How do I sample a line across a blob at a perpendicular angle? (in Python/OpenCV unless you suggest switching to something else)

Question

I am working on a program that will trace dots down the center of binary blobs resembling curved confetti pieces. Later I will fit these points with a cubic spline tracing the curve.

As part of the program, I need to:

-create a 2D vector sampling an angled line across a binary image,

-calculate the angle to use at each position along the confetti blob.

Here are some examples of the images, and sketches of what the points traced might look like:

Finding the center of a vertical section of a black confetti is straightforward. Provided you are in a black pixel, find the left and right white edge, the middle is half the distance between these. Doing this is easy because the 2d vector used to make the above calculation is just a row of the image.

But the confetti pieces do not always line up straight and vertically! Sometimes they are curved, or aligned horizontally. What's needed here is a 2d vector that cuts a section through the confetti at an angle. What is the most efficient way to sample this angled vector from the image? Practically, in an image processing library such as Python PIL or OpenCV, are there operations that can get vectors of lines at angles through an image? If I make one, how can I take care to make sure it is efficient?

What is the most efficient way to calculate the angle of the vector needed? One way to get the appropriate angle is to find the angle that results in the minimum width black segment in the returned 2d vector. I don't need to do this exhaustively, only cycle through 360 degrees at 30 degree increments Another way to get the appropriate angle might be to find the tangent of the curve of the confetti piece, and use the line perpendicular to that - but that might be more complicated.

Any thoughts on how to better tackle the problem would be very helpful. Any specific suggestions regarding how to fetch a 2d line across an image, and an efficient way to get the perpendicular angle would also be great.

Answer 1

It would appear that you are interested in "medial axis fitting" - along with an estimate of orientation (if you have the axis then usually the tangent of the axis at any point suffices).

Technically, with OpenCV, you might consider using a distance transform (cv.PointPolygonTest), using Voronoi cells (cv.cvCalcSubdivVoronoi2D), or - as suggested by @remi morphological thinning...

But, if you didn't want to use the scikits-image package, and simply had to use OpenCV - here's a starter attempt using some skeletonization code (based on a quick and easy technique).

You could then follow this by some spline fitting along the discovered points to work out your samples and tangents (but this would require a bit more work to discover the ends and to remove any stray points/gaps...)

import cv

# get images 
orig = cv.LoadImage('o1Mlp.png')

# create storage images
grey = cv.CreateImage(cv.GetSize(orig), 8, 1) 
skel = cv.CreateImage(cv.GetSize(orig),8, 1) 
temp = cv.CreateImage(cv.GetSize(orig),8,1)

# convert image to pure binary B&W based on threshold
cv.CvtColor(orig, grey, cv.CV_RGB2GRAY)
cv.Threshold(grey,grey,200,255,cv.CV_THRESH_BINARY_INV)

# Create cross operator - good for skeletonization 
elem = cv.CreateStructuringElementEx(3,3,1,1,cv.CV_SHAPE_CROSS)

# Loop until we have consumed all the values in the original image
while True:
  cv.MorphologyEx(grey,temp,None,elem,cv.CV_MOP_OPEN) # Shrink..
  cv.Not(temp,temp) # ...invert...
  cv.And(grey,temp,temp) # ...intersect with original...
  cv.Or(skel,temp,skel) # ... add to current skeleton...
  cv.Erode(grey,grey,elem) # and reduce original ready for next.

  (minVal,maxVal,minLoc,maxLoc)= cv.MinMaxLoc(grey)
  if (maxVal==0): # Any pixels left?
    break

# show result 
cv.ShowImage("orig", orig)
cv.ShowImage("skel", skel)
cv.WaitKey(-1)

Answer 2

Concerning the last part of the problem, finding the normals: I use my own algorithm. It seems to work. If you find a standard solution or improve on mine, please let us know!

    /// <summary>
    /// Returns absolute angle between points at offset length, or double.MinValue when not found.
    /// </summary>
    /// <param name="sequence">Ordered array of points (e.g., from OpenCV Contour)</param>
    /// <param name="length">Number of points used to calculate angle</param>
    /// /// <param name="increment">number of points between each angle calculation (e.g., 1 to attempt to determine angles for all points)</param>
    /// <returns></returns>
    public static double[] FindAbsoluteAngles(Point[] sequence, int length, int increment)
    {
        double[] angles = new double[sequence.Length];
        for (int i = 0; i < sequence.Length; i++)
            angles[i] = double.MinValue;

        double last = double.MinValue;
        for (int i = length; i < sequence.Length; i += increment)
        {
            int i1 = i - length;
            int i2 = i - ((int)length / 2);
            int i3 = i;

            Point p1 = sequence[i1];
            Point p2 = sequence[i2];
            Point p3 = sequence[i3];

            if (p1.X != p3.X & p1.Y != p3.Y)//Is a diagonal
            {
                angles[i2] = 180 - Math.Atan(1.0 * (p1.X - p3.X) / (p1.Y - p3.Y)) * 180 / Math.PI;
            }
            else if (last != double.MinValue)
            {
                //USe previous angle to determine non-diagonals (which can be: 0 or 180; 90 or 270) 
                double error;
                if (p1.X == p3.X)//Is a vertical
                {
                    error = Math.Abs(last - 180);
                    if (Math.Min(error, 180 - error) == error)
                        angles[i2] = 180;
                    else
                        angles[i2] = 0;
                }
                else if (p1.Y == p3.Y)//Is a horizontal
                {
                    error = Math.Abs(last - 270);
                    if (Math.Min(error, 180 - error) == error)
                        angles[i2] = 270;
                    else
                        angles[i2] = 90;
                }
            }

            last = angles[i2];
        }

        return angles;
    }