What is the best way to find the closest match to a complex shape, using opencv and c++?

Question

Alright, here is my source code. This code will take an image in a file and compare it against a list of images in another file. In the file of images you must include a .txt file containing the names of all of the images in the file you are trying to compare. The problem i'm having is that these two images are very similar but are not exactly the same. I need a method to refine these matches further. Perhaps even an entire new way to compare these two shapes (in larger chunks, blobs, ect). One way I was considering is actually making an entire keypoint map, and only comparing keypionts if they are at or near a certain point that corespondes to both images. Ie: compare keypoints at point (12,200), +-10 pixels from (x, y) and see if there are similar keypoint on the other image.

All I need is a way to get the best matches possible from: ActualImplant and XrayOfThatSameImplantButASlightlyDifferentSize. Please and thank you!

PS: you will see commented out sections where I was experimenting with Sobel Derivatives and other such things. I ended up just adjusting contrast and brightness on the xray for the best outline. The same has to be done to the image of the implant before it is used to try to match anything.

#include "opencv2\highgui\highgui.hpp"
#include "opencv2\features2d\features2d.hpp"
#include "opencv2\imgproc.hpp"


#include <iostream>
#include <fstream>
#include <ctime>
const string defaultDetector = "ORB";
const string defaultDescriptor = "ORB";
const string defaultMatcher = "BruteForce-Hamming";
const string defaultXrayImagePath = "../../xray.png";
const string defaultImplantImagesTextListPath = "../../implantImage.txt";
const string defaultPathToResultsFolder = "../../results";

static void printIntro(const string& appName)
{
    cout << "/*                                                                                                       *\n"
        << " * Created by: Alex Gatz. 1/11/12. Created for: Xray Implant Identification                               *\n"
        << " * This code was created to scan a file full of images of differnt implants, generate keypoint maps       *\n"
        << " * for each image, and identifywhich image most closely matches a chosen image in another folder          *\n"
        << " */                                                                                                       *\n"
        << endl;

    cout << endl << "Format:\n" << endl;
    cout << "./" << appName << " [detector] [descriptor] [matcher] [xrayImagePath] [implantImagesTextListPath] [pathToSaveResults]" << endl;
    cout << endl;

    cout << "\nExample:" << endl
        << "./" << appName << " " << defaultDetector << " " << defaultDescriptor << " " << defaultMatcher << " "
        << defaultXrayImagePath << " " << defaultImplantImagesTextListPath << " " << defaultPathToResultsFolder << endl;
}

static void maskMatchesByImplantImgIdx(const vector<DMatch>& matches, int trainImgIdx, vector<char>& mask)
{
    mask.resize(matches.size());
    fill(mask.begin(), mask.end(), 0);
    for (size_t i = 0; i < matches.size(); i++)
    {
        if (matches[i].imgIdx == trainImgIdx)
            mask[i] = 1;
    }
}

static void readImplantFilenames(const string& filename, string& dirName, vector<string>& implantFilenames)
{
    implantFilenames.clear();

    ifstream file(filename.c_str());
    if (!file.is_open())
        return;

    size_t pos = filename.rfind('\\');
    char dlmtr = '\\';
    if (pos == String::npos)
    {
        pos = filename.rfind('/');
        dlmtr = '/';
    }
    dirName = pos == string::npos ? "" : filename.substr(0, pos) + dlmtr;

    while (!file.eof())
    {
        string str; getline(file, str);
        if (str.empty()) break;
        implantFilenames.push_back(str);
    }
    file.close();
}

static bool createDetectorDescriptorMatcher(const string& detectorType, const string& descriptorType, const string& matcherType,
    Ptr<FeatureDetector>& featureDetector,
    Ptr<DescriptorExtractor>& descriptorExtractor,
    Ptr<DescriptorMatcher>& descriptorMatcher)
{
    cout << "< Creating feature detector, descriptor extractor and descriptor matcher ..." << endl;
    featureDetector = ORB::create( //All of these are parameters that can be adjusted to effect match accuracy and process time.
        10000, //int nfeatures = Maxiumum number of features to retain; max vaulue unknown, higher number takes longer to process. Default: 500
        1.4f, //float scaleFactor= Pyramid decimation ratio; between 1.00 - 2.00.                                                 Default: 1.2f
        6,   //int nlevels = Number of pyramid levels used; more levels more time taken to process, but more accurate results.   Default: 8
        40,   //int edgeThreshold = Size of the border where the features are not detected. Should match patchSize roughly.       Default: 31
        0,    //int firstLevel = Should remain 0 for now.                                                                         Default: 0
        4,    //int WTA_K = Should remain 2.                                                                                      Default: 2
        ORB::HARRIS_SCORE,    //int scoreType = ORB::HARRIS_SCORE is the most accurate ranking possible for ORB.                                          Default: HARRIS_SCORE
        33    //int patchSize = size of patch used by the oriented BRIEF descriptor. Should match edgeThreashold.                 Default: 31
        ); 
    //featureDetector = ORB::create(); // <-- Uncomment this and comment the featureDetector above for default detector-
    //OpenCV 3.1 got rid of the dynamic naming of detectors and extractors. 

    //These two are one in the same when using ORB, some detectors and extractors are separate
    // in which case you would set "descriptorExtractor = descriptorType::create();" or its equivilant.
    descriptorExtractor = featureDetector;

    descriptorMatcher = DescriptorMatcher::create(matcherType);

    cout << ">" << endl;

    bool isCreated = !(featureDetector.empty() || descriptorExtractor.empty() || descriptorMatcher.empty());
    if (!isCreated)
        cout << "Can not create feature detector or descriptor extractor or descriptor matcher of given types." << endl << ">" << endl;

    return isCreated;
}

static void manipulateImage(Mat& image) //Manipulates images into only showing an outline!
{
    //Sobel Dirivative edge finder

    //int scale = 1;
    //int delta = 0;
    //int ddepth = CV_16S;
    ////equalizeHist(image, image); //This will equilize the lighting levels in each image.
    //GaussianBlur(image, image, Size(3, 3), 0, 0, BORDER_DEFAULT);

    //Mat grad_x, grad_y;
    //Mat abs_grad_x, abs_grad_y;
    ////For x
    //Sobel(image, grad_x, ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT);
    //convertScaleAbs(grad_x, abs_grad_x);
    ////For y
    //Sobel(image, grad_y, ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT);
    //convertScaleAbs(grad_y, abs_grad_y);

    //addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0, image);

    //Specific Level adjustment (very clean)
    double alpha = 20; //Best Result: 20
    int beta = -300;   //Best Result: -300
    image.convertTo(image, -1, alpha, beta);
}

static bool readImages(const string& xrayImageName, const string& implantFilename,
    Mat& xrayImage, vector <Mat>& implantImages, vector<string>& implantImageNames)
{
    //TODO: Add a funtion call to automatically adjust all images loaded to best settings for matching. 
    cout << "< Reading the images..." << endl;
    xrayImage = imread(xrayImageName, CV_LOAD_IMAGE_GRAYSCALE); //Turns the image gray while loading.
    manipulateImage(xrayImage); //Runs image manipulations

    if (xrayImage.empty())
    {
        cout << "Xray image can not be read." << endl << ">" << endl;
        return false;
    }
    string trainDirName;
    readImplantFilenames(implantFilename, trainDirName, implantImageNames);
    if (implantImageNames.empty())
    {
        cout << "Implant image filenames can not be read." << endl << ">" << endl;
        return false;
    }
    int readImageCount = 0;
    for (size_t i = 0; i < implantImageNames.size(); i++)
    {
        string filename = trainDirName + implantImageNames[i];
        Mat img = imread(filename, CV_LOAD_IMAGE_GRAYSCALE); //Turns imamges gray while loading.
        //manipulateImage(img); //Runs Sobel Dirivitage on implant image.
        if (img.empty())
        {
            cout << "Implant image " << filename << " can not be read." << endl;
        }
        else
        {
            readImageCount++;
        }
        implantImages.push_back(img);
    }
    if (!readImageCount)
    {
        cout << "All implant images can not be read." << endl << ">" << endl;
        return false;
    }
    else
        cout << readImageCount << " implant images were read." << endl;
    cout << ">" << endl;

    return true;
}

static void detectKeypoints(const Mat& xrayImage, vector<KeyPoint>& xrayKeypoints,
    const vector<Mat>& implantImages, vector<vector<KeyPoint> >& implantKeypoints,
    Ptr<FeatureDetector>& featureDetector)
{
    cout << endl << "< Extracting keypoints from images..." << endl;
    featureDetector->detect(xrayImage, xrayKeypoints);
    featureDetector->detect(implantImages, implantKeypoints);
    cout << ">" << endl;
}

static void computeDescriptors(const Mat& xrayImage, vector<KeyPoint>& implantKeypoints, Mat& implantDescriptors,
    const vector<Mat>& implantImages, vector<vector<KeyPoint> >& implantImageKeypoints, vector<Mat>& implantImageDescriptors,
    Ptr<DescriptorExtractor>& descriptorExtractor)
{
    cout << "< Computing descriptors for keypoints..." << endl;
    descriptorExtractor->compute(xrayImage, implantKeypoints, implantDescriptors);
    descriptorExtractor->compute(implantImages, implantImageKeypoints, implantImageDescriptors);

    int totalTrainDesc = 0;
    for (vector<Mat>::const_iterator tdIter = implantImageDescriptors.begin(); tdIter != implantImageDescriptors.end(); tdIter++)
        totalTrainDesc += tdIter->rows;

    cout << "Query descriptors count: " << implantDescriptors.rows << "; Total train descriptors count: " << totalTrainDesc << endl;
    cout << ">" << endl;
}

static void matchDescriptors(const Mat& xrayDescriptors, const vector<Mat>& implantDescriptors,
    vector<DMatch>& matches, Ptr<DescriptorMatcher>& descriptorMatcher)
{
    cout << "< Set implant image descriptors collection in the matcher and match xray descriptors to them..." << endl;
    //time_t timerBegin, timerEnd;

    //time(&timerBegin);
    descriptorMatcher->add(implantDescriptors);
    descriptorMatcher->train();
    //time(&timerEnd);
    //double buildTime = difftime(timerEnd, timerBegin);

    //time(&timerBegin);
    descriptorMatcher->match(xrayDescriptors, matches);
    //time(&timerEnd);
    //double matchTime = difftime(timerEnd, timerBegin);

    CV_Assert(xrayDescriptors.rows == (int)matches.size() || matches.empty());

    cout << "Number of imageMatches: " << matches.size() << endl;
    //cout << "Build time: " << buildTime << " ms; Match time: " << matchTime << " ms" << endl;
    cout << ">" << endl;
}

static void saveResultImages(const Mat& xrayImage, const vector<KeyPoint>& xrayKeypoints,
    const vector<Mat>& implantImage, const vector<vector<KeyPoint> >& implantImageKeypoints,
    const vector<DMatch>& matches, const vector<string>& implantImagesName, const string& resultDir)
{
    cout << "< Save results..." << endl;
    Mat drawImg;
    vector<char> mask;
    for (size_t i = 0; i < implantImage.size(); i++)
    {
        if (!implantImage[i].empty())
        {
            maskMatchesByImplantImgIdx(matches, (int)i, mask);
            drawMatches(xrayImage, xrayKeypoints, implantImage[i], implantImageKeypoints[i],
                matches, drawImg, Scalar::all(-1), Scalar(0, 0, 255), mask, 4);
            string filename = resultDir + "/result_" + implantImagesName[i];
            if (!imwrite(filename, drawImg))
                cout << "Image " << filename << " can not be saved (may be because directory " << resultDir << " does not exist)." << endl;
        }
    }
    cout << ">" << endl;

    //After all results have been saved, another function will scan and place the final result in a separate folder. 
    //For now this save process is required to manually access each result and determine if the current settings are working well. 
}

int main(int argc, char** argv)
{
    //Intialize variables to global defaults.
    string detector = defaultDetector;
    string descriptor = defaultDescriptor;
    string matcher = defaultMatcher;
    string xrayImagePath = defaultXrayImagePath;
    string implantImagesTextListPath = defaultImplantImagesTextListPath;
    string pathToSaveResults = defaultPathToResultsFolder;

    //As long as you have 7 arguments, you can procede
    if (argc != 7 && argc != 1)
    {
        //This will be called if the incorrect amount of commands are used to start the program.
        printIntro(argv[1]);
        system("PAUSE");
        return -1;
    }

    //As long as you still have 7 arguments, I will set the variables for this
    // to the arguments you decided on. 
    //If testing using XrayID --> Properties --> Debugging --> Command Arguments, remember to start with [detector] as the first command
    // C++ includes the [appName] command as the first argument automantically. 

    if (argc != 1) //I suggest placing a break here and stepping through this to ensure the proper commands were sent in. With a
                   // GUI this would nto matter because the GUI would structure the input and use a default if no input was used. 
    {
        detector = argv[1]; 
        descriptor = argv[2]; 
        matcher = argv[3];
        xrayImagePath = argv[4]; 
        implantImagesTextListPath = argv[5];
        pathToSaveResults = argv[6];
    }

    //Set up cv::Ptr's for tools. 
    Ptr<FeatureDetector> featureDetector;
    Ptr<DescriptorExtractor> descriptorExtractor;
    Ptr<DescriptorMatcher> descriptorMatcher;

    //Check to see if tools are created, if not true print intro and close program.
    if (!createDetectorDescriptorMatcher(detector, descriptor, matcher, featureDetector, descriptorExtractor, descriptorMatcher))
    {
        printIntro(argv[0]);
        system("PAUSE");
        return -1;
    }

    Mat testImage;
    vector<Mat> implantImages;
    vector<string> implantImagesNames;

    //Check to see if readImages completes properly, if not true print intro and close program. 
    if (!readImages(xrayImagePath, implantImagesTextListPath, testImage, implantImages, implantImagesNames))
    {
        printIntro(argv[0]);
        system("PAUSE");
        return -1;
    }

    vector<KeyPoint> xrayKeypoints;
    vector<vector<KeyPoint> > implantKeypoints;
    detectKeypoints(testImage, xrayKeypoints, implantImages, implantKeypoints, featureDetector);

    Mat xrayDescriptors;
    vector<Mat> implantTestImageDescriptors;
    computeDescriptors(testImage, xrayKeypoints, xrayDescriptors, implantImages, implantKeypoints, implantTestImageDescriptors,
        descriptorExtractor);

    vector<DMatch> imageMatches;
    matchDescriptors(xrayDescriptors, implantTestImageDescriptors, imageMatches, descriptorMatcher);
    saveResultImages(testImage, xrayKeypoints, implantImages, implantKeypoints, imageMatches, implantImagesNames, pathToSaveResults);

    system("PAUSE");
    return 0;
}

Answer 1

• Try below code.Hope this will help you.

  #include <opencv2/nonfree/nonfree.hpp>
  #include <iostream>
  #include <dirent.h>
  #include <ctime>  
  #include <stdio.h>
  using namespace cv;
  using namespace std;
  
  int main(int argc, const char *argv[])
  {
     double ratio = 0.9;
  
     Mat image1 = imread("Image1_path);
     Mat image2 = cv::imread("Image2_path");
  
     Ptr<FeatureDetector> detector;
     Ptr<DescriptorExtractor> extractor;
  
    // TODO default is 500 keypoints..but we can change
    detector = FeatureDetector::create("ORB");
    extractor = DescriptorExtractor::create("ORB");
  
    vector<KeyPoint> keypoints1, keypoints2;
    detector->detect(image1, keypoints1);
    detector->detect(image2, keypoints2);
  
    cout << "# keypoints of image1 :" << keypoints1.size() << endl;
    cout << "# keypoints of image2 :" << keypoints2.size() << endl;
  
    Mat descriptors1,descriptors2;
    extractor->compute(image1,keypoints1,descriptors1);
    extractor->compute(image2,keypoints2,descriptors2);
  
    cout << "Descriptors size :" << descriptors1.cols << ":"<< descriptors1.rows << endl;
  
    vector< vector<DMatch> > matches12, matches21;
    Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");
    matcher->knnMatch( descriptors1, descriptors2, matches12, 2);
    matcher->knnMatch( descriptors2, descriptors1, matches21, 2);
  
    //BFMatcher bfmatcher(NORM_L2, true);
    //vector<DMatch> matches;
    //bfmatcher.match(descriptors1, descriptors2, matches);
    double max_dist = 0; double min_dist = 100;
    for( int i = 0; i < descriptors1.rows; i++)
    {
        double dist = matches12[i].data()->distance;
        if(dist < min_dist)
           min_dist = dist;
        if(dist > max_dist)
           max_dist = dist;
    }
    printf("-- Max dist : %f \n", max_dist);
    printf("-- Min dist : %f \n", min_dist);
    cout << "Matches1-2:" << matches12.size() << endl;
    cout << "Matches2-1:" << matches21.size() << endl;
  
    std::vector<DMatch> good_matches1, good_matches2;
    for(int i=0; i < matches12.size(); i++)
    {
        if(matches12[i][0].distance < ratio * matches12[i][1].distance)
           good_matches1.push_back(matches12[i][0]);
    }
  
    for(int i=0; i < matches21.size(); i++)
    {
        if(matches21[i][0].distance < ratio * matches21[i][1].distance)
           good_matches2.push_back(matches21[i][0]);
    }
  
    cout << "Good matches1:" << good_matches1.size() << endl;
    cout << "Good matches2:" << good_matches2.size() << endl;
  
   // Symmetric Test
   std::vector<DMatch> better_matches;
   for(int i=0; i<good_matches1.size(); i++)
   {
       for(int j=0; j<good_matches2.size(); j++)
       {
           if(good_matches1[i].queryIdx == good_matches2[j].trainIdx && good_matches2[j].queryIdx == good_matches1[i].trainIdx)
           {
               better_matches.push_back(DMatch(good_matches1[i].queryIdx, good_matches1[i].trainIdx, good_matches1[i].distance));
          break;
           }
       }
   }
  
   cout << "Better matches:" << better_matches.size() << endl;
   double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
  
   // show it on an image
   Mat output;
   drawMatches(image1, keypoints1, image2, keypoints2, better_matches, output);
   imshow("Matches result",output);
   waitKey(0);
  
   return 0;
  }
  

Answer 2

That image looks rather like an artificial hip. If you're dealing with medical images, you should definitely check out The Insight Toolkit (ITK) which has many special features designed for the particular needs of this domain. You could do a simple Model-Image Registration between your real-world image and your template data to find the best result. I think you would get much better results with this approach than with the point-based testing described above.

This sort of registration performs an iterative optimisation of a set of parameters (in this case, an affine transform) which seeks to find the best mapping of the model to the image data.

• ITK Affine Registration example

The example above takes a fixed image and attempts to find a transform that maps the moving image onto it. The transform is a 2D affine transform (rotation and translation in this case) and its parameters are the result of running the optimiser which maximises the matching metric. The metric measures how well the fixed image and the transformed moving image match. The interpolator is what takes the moving image and applies the transform to map it onto the fixed image.

In your sample images, fixed image could be the original X-ray and the moving image the actual implant. You will probably need to add scaling to make a full affine transform since the size of the two differs.

The metric is a measure of how well the transformed moving image matches the fixed image, so you would need to determine a tolerance or minimum metric for a match to be valid. If the images are very different, the metric would be very low and can be rejected.

The output is a set of transformation parameters and the output image is the final optimal transform applied to the moving image (not a combination of the images). The result is basically telling you where the implant is found in the X-ray.