Susan Point Matching

[fdncred]

What would you like to submit? (put an 'x' inside the bracket that applies)

• [x] question
• [ ] bug report
• [ ] feature request

Issue description @cesarsouza, I have found that SusanCornersDetector() seems to work well with my matching endeavor. However, correlation matching is extremely slow. I tried parallelizing the code but it just messed up things due to my ignorace with parallelism. Are there any other alternatives to matching these points that are faster? I tried to use KNearestNeighborMatching() but it wants a feature with descriptors versus just points.

var susan = new SusanCornersDetector();
var detector = new CornerFeaturesDetector(susan);
var susanPoints1 = detector.Transform(needle).Select(x => new IntPoint((int)x.X, (int)x.Y));
var susanPoints2 = detector.Transform(hayCrop).Select(x => new IntPoint((int)x.X, (int)x.Y));
var matcher = new CorrelationMatching(9, needle, hayCrop);
var matches = matcher.Match(susanPoints1, susanPoints2);

Thanks, Darren

[cesarsouza]

Hi Darren,

Sorry for the delay in answering! For the moment no, if you only have a points detector like Susan the current only approach is to use correlation matching. It could be possible to use KNearestNeighborMatching, but I suspect it would be just as slow because first we would need to extract pixel windows around the corner points anyways. However, it might be possible to parallelize it to some extent. I will add this to the backlog as an enhancement/feature request, but I can't promise when it would be done.

(if you would like to share your work on trying to parallelize it so far by the way, I could try to give some hints on why it is not working so well)

Regards, Cesar

[fdncred]

Hey @cesarsouza, Good to hear from you. Unfortunately I've had to move on from Accord.Net due to speed issues. I'm using OpenCvSharp now for TemplateMatching and KazeFeatures. I haven't tried Susan Corners & CorrelationMatching yet in OpenCvSharp but I assume it'll be faster as well.

However, for the off chance that this thread may help others, I'm posting the code below that is taking the most time.

We start with the method that looks for a "needle" in a "haystack". The hayCrop variable is the haystack but cropped to a smaller ROI to help searching and matching go faster.

var modelSubList = modelFull.Take(10).ToList();
for (int i = 0; i < modelSubList.Count; i++)
{
    var curTemplate = modelSubList[i];
    using (needle = modelSubList[i].TemplateROIBitmap.Resize(scale))
    using (hayCrop = haystackData.SeedBitmap.Crop(modelSubList[i].SearchROI).Resize(scale))
    {
        // Step 1: Detect feature points using Harris Corners Detector
        var susanPoints1 = detector.Transform(needle).Select(x => new IntPoint((int)x.X, (int)x.Y));
        var susanPoints2 = detector.Transform(hayCrop).Select(x => new IntPoint((int)x.X, (int)x.Y));
        if (susanPoints1.Count() < 4 || susanPoints2.Count() < 4) continue;
        // Step 2: Match feature points using a correlation measure
        var matcher = new CorrelationMatching(7, /*400,*/ needle, hayCrop);
        var matches = matcher.Match(susanPoints1, susanPoints2);
        if (matches[0].Length < 4 || matches[1].Length < 4) continue;
        // Step 3: Create the homography matrix using a robust estimator
        var ransac = new RansacHomographyEstimator(0.001, 0.99);
        var homography = ransac.Estimate(matches[0], matches[1]);
        var foundAtX = (int)Math.Round(homography.OffsetX / scale + modelSubList[i].SearchROI.X);
        var foundAtY = (int)Math.Round(homography.OffsetY / scale + modelSubList[i].SearchROI.Y);
        var offsetPt = new Accord.Point(foundAtX - modelSubList[i].TemplateROI.X, foundAtY - modelSubList[i].TemplateROI.Y);
        offsetPoints.Add(offsetPt);
    }
}

The "new CorrelationMatching()" statement eventually calls computeCorrelationMatrix and this is where all the time is spent.

// Generate correlation matrix
double[,] correlationMatrix =
    computeCorrelationMatrix(grayImage1, points1, grayImage2, points2, windowSize, dmax);

Below is my modified computeCorrelationMatrix() method. You can see where I have commented out stop watches to instrument things as well as the Parllel.ForEach'es. You may also note CalcWindow1 and CalcWindow2 methods commented out. They just pulled that loop below them into a separate function where I Parallel.For'd through them. Surprisingly, the CalcWindow methods ended up making computeCorrelationMatrix many times slower. Most of the time in this method is what I labeled as loop3, contructing the matrix. I just couldn't get it to work faster and when I did get things working (read: not throwing exceptions) the matrix was all wrong. I chocked that up to threads writing things in the wrong order. I considered putting all this mess in a threadsafe ordered dictionary but at that point I was just tired of messing with it. So I moved on.

private static double[,] computeCorrelationMatrix(
    Bitmap image1, IntPoint[] points1,
    Bitmap image2, IntPoint[] points2,
    int windowSize, double maxDistance)
{
    var LockObj = new Object();
    // Create the initial correlation matrix
    double[,] matrix = Matrix.Create(points1.Length, points2.Length, Double.NegativeInfinity);

    // Gather some information
    int width1 = image1.Width;
    int width2 = image2.Width;
    int height1 = image1.Height;
    int height2 = image2.Height;

    int r = (windowSize - 1) / 2; //  'radius' of correlation window
    double m = maxDistance * maxDistance; // maximum considered distance
    double[,] w1 = new double[windowSize, windowSize]; // first window
    double[,] w2 = new double[windowSize, windowSize]; // second window

    // Lock the images
    BitmapData bitmapData1 = image1.LockBits(new Rectangle(0, 0, width1, height1),
        ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);
    BitmapData bitmapData2 = image2.LockBits(new Rectangle(0, 0, width2, height2),
        ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);

    int stride1 = bitmapData1.Stride;
    int stride2 = bitmapData2.Stride;

    // We will ignore points at the edge
    int[] idx1 = Matrix.Find(points1, p => p.X >= r && p.X < width1 - r &&
                                           p.Y >= r && p.Y < height1 - r);

    int[] idx2 = Matrix.Find(points2, p => p.X >= r && p.X < width2 - r &&
                                           p.Y >= r && p.Y < height2 - r);

    //var loop1 = new List<long>();
    //var loop2 = new List<long>();
    //var loop3 = new List<long>();
    //var sw = new Stopwatch();
    //var overall = new Stopwatch();
    //overall.Start();
    // For each index in the first set of points
    foreach (int n1 in idx1)
    //Parallel.ForEach(idx1, (n1) =>
    {
        //sw.Reset();
        //sw.Start();
        // Get the current point
        IntPoint p1 = points1[n1];
        double sum = 0;

        unsafe // Create the first window for the current point
        {
            byte* src = (byte*)bitmapData1.Scan0 + (p1.X - r) + (p1.Y - r) * stride1;
            //var res1 = CalcWindow1(windowSize, src, stride1);
            for (int j = 0; j < windowSize; j++)
            {
                for (int i = 0; i < windowSize; i++)
                {
                    double w = (byte)(*(src + i));
                    w1[i, j] = w;
                    sum += w * w;
                }
                src += stride1;
            }
            //w1 = res1.W1;
            //sum = res1.Sum;
        }

        // Normalize the window
        w1.Divide(Math.Sqrt(sum), result: w1);
        //sw.Stop();
        ////Console.WriteLine($"Loop1:[{sw.ElapsedMilliseconds}]");
        //loop1.Add(sw.ElapsedMilliseconds);
        //sw.Reset();
        //sw.Start();
        // Identify the indices of points in p2 that we need to consider.
        int[] candidates;
        if (maxDistance == 0)
        {
            // We should consider all points
            candidates = idx2;
        }
        else
        {
            // We should consider points that are within
            //  distance maxDistance apart

            // Compute distances from the current point
            //  to all points in the second image.
            double[] distances = new double[idx2.Length];
            for (int i = 0; i < idx2.Length; i++)
            {
                double dx = p1.X - points2[idx2[i]].X;
                double dy = p1.Y - points2[idx2[i]].Y;
                distances[i] = dx * dx + dy * dy;
            }

            candidates = idx2.Get(Matrix.Find(distances, d => d < m));
        }
        //sw.Stop();
        ////Console.WriteLine($"Loop2: Distance:[{sw.ElapsedMilliseconds}]");
        //loop2.Add(sw.ElapsedMilliseconds);
        //sw.Reset();
        //sw.Start();

        // Calculate normalized correlation measure
        foreach (int n2 in candidates)
        //Parallel.ForEach(candidates, n2 =>
        {
            IntPoint p2 = points2[n2];

            unsafe // Generate window in 2nd image
            {
                byte* src = (byte*)bitmapData2.Scan0 + (p2.X - r) + (p2.Y - r) * stride2;
                //w2 = CalcWindow2(windowSize, src, stride2);
                for (int j = 0; j < windowSize; j++)
                {
                    for (int i = 0; i < windowSize; i++)
                        w2[i, j] = (byte)(*(src + i));
                    src += stride2;
                }
            }

            double sum1 = 0, sum2 = 0;
            for (int i = 0; i < windowSize; i++)
            {
                for (int j = 0; j < windowSize; j++)
                {
                    sum1 += w1[i, j] * w2[i, j];
                    sum2 += w2[i, j] * w2[i, j];
                }
            }

            matrix[n1, n2] = sum1 / Math.Sqrt(sum2);
        }//);
        //sw.Stop();
        //loop3.Add(sw.ElapsedMilliseconds);
    }//);
    //overall.Stop();
    //Console.WriteLine($"Loop1: Window:[{loop1.Average()}]");
    //Console.WriteLine($"Loop2: Distance:[{loop2.Average()}]");
    //Console.WriteLine($"Loop3: Candidates:[{loop3.Average()}]");
    //Console.WriteLine($"Overall:[{overall.ElapsedMilliseconds}]");
    // Release the images
    image1.UnlockBits(bitmapData1);
    image2.UnlockBits(bitmapData2);

    return matrix;
}