Finding circle in objects touching each other

Question

I'm working on an app to measure radius of circle. The cylinder is placed on 2 supporting pieces of plastic. When I threshold the image, the supporting pieces comes on as part of the cylinder. I tried to use contour detecting and use minimum enclosing circles to find the largest and second largest object, but the outer circle is always incorrect due to the supporting objects. I'm using c# and emgucv. Morphological operations wont work because it will damage the edges of the circle and measurement will be incorrect. How do I detect and measure the circle properly? Any solution in Python and OpenCV will also work.

Here is the code:

 VectorOfVectorOfPoint contours = new VectorOfVectorOfPoint();
  CvInvoke.FindContours(binaryimg, contours, null, RetrType.Tree, 
 ChainApproxMethod.ChainApproxNone);
 VectorOfPoint largestContour = null;
 VectorOfPoint secondLargestContour = null;
 double maxArea = 0;
 double secondMaxArea = 0;

 if (contours.Size < 2)
{
Console.WriteLine("Not enough contours found to detect both OD and ID.");

}
      
for (int i = 0; i < contours.Size; i++)
{
   using (VectorOfPoint contour = contours[i])
 

{
    double area = CvInvoke.ContourArea(contour);

    if (area > maxArea)
    {
        secondMaxArea = maxArea;
        secondLargestContour = largestContour;
        maxArea = area;
        largestContour = contour;
    }
    else if (area > secondMaxArea)
    {
        secondMaxArea = area;
        secondLargestContour = contour;
    }
   }
  }


 if (largestContour != null && secondLargestContour !=null)
  {
    CircleF outerCircle = CvInvoke.MinEnclosingCircle(largestContour);
    CircleF innerCircle = CvInvoke.MinEnclosingCircle(secondLargestContour);
  }

This is the output im getting:

Here is the original image:

Answer 1

You can clip the image (red rectangle) or fill the bottom with white color. The contour you will get is not a full circle but you can reconstruct the outline. To do that I see two strategies: Using hough circles or the enclosing circle of the blob object. For Blob see FindContourTree and contour features. Below is the pseudo code:

roi = new Mat(image, region)
blobObjects = CvInvoke.FindContourTree(roi)
blackObject = biggest black object of blobObjects
whiteObjcet = biggest white object inside blackObject 
circle1 = CvInvoke.EnclosingCircle(blackObject.Arc)
circle2 = CvInvoke.EnclosingCircle(whiteObjcet.Arc)

• If the object location is not constant you can first search for the inner circle.
• The enclosing circle will have a pixel accuracy for the Blob and subpixel accuracy for the HoughCircles solution
• FindContours is the fastes solution compared to the execution time

Answer 2

General note: approaches can be combined to use their strengths and compensate their weaknesses. Some are good at finding circles quickly but only approximately. Others can locate a circle precisely, but they take a long time finding a circle, if at all.

You could use Hough to find the circles in a low resolution picture (for speed), then select edge points from the full resolution image that are close to such a circle, and do a model fit on that.

---

The inner circle is perfect and the outer circle is 90% perfect. That means you could try to fit circle models to subsets of the edge points. That is RANSAC.

RANSAC is usually only used when there is a single object, with some outliers, but the object points are the majority (i.e. 50% ore more) of all points. In this case, both of those assumptions are broken. RANSAC can still work, but it will do badly or cost a lot.

RANSAC is a general algorithm:

1. sample a minimal random subset (here: 3 points), hoping they all belong to the object in question
2. fit the model to that sample set
3. determine set of inliers (points matching the model)
4. maybe refit the model on the inlier set
5. score the value of the model (such as how well the circle is covered)
6. repeat N times and take the best model, or repeat until good enough

RANSAC being random means you have to accept one of two downsides:

• if you bound the number of iterations, the program may not find anything good in the allotted time, even if it's there
• if you require a minimum quality result, the algorithm might run for an unpredictable number of iterations, potentially even forever

To find several circles, first find one, remove its points from the set of points, repeat, until the algorithm takes "long enough" to find nothing more. This too can fail, in that it might not find all the circles.

There are simple equations to determine the number of iterations, if you can estimate the probability of a number of points belonging to a circle, and you can say with what probability you want to get a result.

If we have 11000 points in total, and the big circle consists of ~3000 points, and we need three from the same circle, that's (3000/11000)**3 = 0.02, and we want a result with 99% probability, then log(1-0.99) / log(1-0.02) = 228 iterations would be needed. If you wanted to find a circle of just 1000 points (the small one) from the remaining set of 11-3 = 8000 points with 99% probability, you'd have to run over 2300 iterations.

---

I can imagine improvements in a few places. I haven't investigated those deeply.

Instead of picking points randomly, one could pick points from the same contour, and perhaps pick consecutive points with some distance between them. That is a lot more likely to grab points from the same circle.

Instead of randomly, that could be done "exhaustively", just walking along every contour with a sliding window, and trying the fit-evaluate approach of RANSAC.

If the circles are known to be concentric, one could "bias" the random hypotheses towards circle centers near known ones.

If this is supposed to run on a video and the circles don't move (much), you could see if the known circles are "still" there.

If the points had an orientation each (from the gradient), that would allow discarding more points, thus rejecting circle models that merely intersect a bunch of points that aren't actually aligned with the circle.

---

I gained edge points from a Scharr-filtered image.

gx = cv.Scharr(im, cv.CV_32F, 1, 0)
gy = cv.Scharr(im, cv.CV_32F, 0, 1)
gm = cv.magnitude(gx, gy)
gm /= gm.max()
edge_mask = (gm >= 0.5)
edge_indices = np.argwhere(edge_mask) # (i,j) pairs
edge_points = edge_indices[:,::-1] # (x,y) pairs

---

Some general functions for the task:

def fit_circle(X, Y):
    A = np.c_[X, Y, np.ones(len(X))] # design matrix A with columns for X, Y, and a constant term
    f = X**2 + Y**2 # get the function as x²+y²
    C, _, _, _ = np.linalg.lstsq(A, f, rcond=None) # optimize
    cx = C[0] / 2 # get the centre x coordinate
    cy = C[1] / 2 # get the centre y coordinate
    radius = np.sqrt(C[2] + cx**2 + cy**2) # calculate the radius
    return cx, cy, radius # return everythin as int

(fit_circle: code from here and math detailed there)

def circle_point_distance(model, points):
    # returns (absolute) distance from circle for each point
    (cx, cy, radius) = model
    return np.abs(np.linalg.norm(points - (cx, cy), axis=1) - radius)

def score_distance_to_circle(model, points):
    (cx, cy, radius) = model
    # score is the sum of squared distances from the circle
    # the lower the score, the better the circle fits the points
    distances = circle_point_distance(model, points)
    return np.sum(distances**2)

# this scores a set of points by the "gaps" it leaves in the coverage of the circle
def score_circle_coverage(model, points, threshold=2):
    (cx, cy, radius) = model
    angles = np.arctan2(points[:,1] - cy, points[:,0] - cx)
    sorted_indices = np.argsort(angles)
    angles = angles[sorted_indices]
    distances = np.roll(angles, -1) - angles
    distances %= 2*np.pi
    distances *= radius
    gap_mask = (distances > threshold)
    cost = distances[gap_mask].sum() / (2*np.pi*radius)
    return 1 - cost

---

The RANSAC function:

def ransac_circle(points, threshold, max_iters=1000, early_stop=None, n_refinements=1):
    n_points = 3 # more just make the set unlikelier to be from the same circle
    best_model = None
    best_score = 0
    best_coverage = 0
    best_numinliers = 0

    for i in range(max_iters):
        sample_indices = np.random.choice(len(points), n_points, replace=False)
        sample = points[sample_indices]

        for _ in range(1 + n_refinements):
            model = fit_circle(sample[:,0], sample[:,1])
            distances = circle_point_distance(model, points)
            sample_mask = (distances <= threshold)
            num_inliers = np.count_nonzero(sample_mask)
            sample = points[sample_mask]

        coverage = score_circle_coverage(model, sample)
        score = num_inliers * coverage

        if score > best_score:
            best_model = model
            best_score = score
            best_coverage = coverage
            best_numinliers = num_inliers

        if early_stop is not None:
            if early_stop(coverage, num_inliers):
                break
    else:
        i = max_iters

    return i, best_model, best_coverage, best_numinliers

def early_stop(coverage, numinliers):
    return (coverage >= 0.8) and (numinliers >= 500)

---

To remove the points fitting one model from the overall set:

distances = circle_point_distance(model, remaining_points)
remaining_points = remaining_points[distances > distance_threshold]

---

Finding all (two) circles would give such a result:

took 0.169695 s
found circle with 90.0% coverage, 3059 inliers for (635.6666319482904, 591.8268363474424, 225.5372141978907)
took 0.232468 s
found circle with 100.0% coverage, 1342 inliers for (630.5985030266783, 601.595411392499, 91.78840587862769)
took 1.161267 s
found nothing

Answer 3

This is a rather simple approach in Python/OpenCV that requires a little estimation for doing Hough Circles. I took the liberty to crop the image first.

The idea is rather simple. Read the input, convert to grayscale, threshold and do HoughCircles for the inner circle. Make the minDist parameter very big as we only want one circle. Estimate the range on the radius. For the second outer circle, first invert the image, then repeat, except after inverting, apply a large morphology open kernel to separate the circle from the rest of the white regions. The repeat as before with new approximate estimates for the radius range.

Input (cropped):

import cv2
import numpy as np

# Read image
img = cv2.imread('2circles.png')
hh, ww = img.shape[:2]

# Convert to grayscale
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# threshold
thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)[1]

# get Hough circles
circles = cv2.HoughCircles(thresh, cv2.HOUGH_GRADIENT, 1, minDist=5000, param1=150, param2=10, minRadius=50, maxRadius=150)
#print(circles)

# draw circles
result = img.copy()
for circle in circles[0]:
    # draw the circle in the output image, then draw a rectangle
    # corresponding to the center of the circle
    (x,y,r) = circle
    print('center1:', (x,y), 'radius1:', r)
    x = int(x)
    y = int(y)
    r = int(r)
    cv2.circle(result, (x, y), r, (0, 0, 255), 1)

# invert thresh
thresh2 = 255 - thresh

# apply morphology to separate outer circle
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (100,100))
morph = cv2.morphologyEx(thresh2, cv2.MORPH_OPEN, kernel)
#morph = cv2.morphologyEx(morph, cv2.MORPH_CLOSE, kernel)

# get Hough circles
circles2 = cv2.HoughCircles(morph, cv2.HOUGH_GRADIENT, 1, minDist=5000, param1=150, param2=10, minRadius=r+20, maxRadius=300)
#print(circles2)

# draw circles
for circle in circles2[0]:
    # draw the circle in the output image, then draw a rectangle
    # corresponding to the center of the circle
    (x,y,r) = circle
    print('center2:', (x,y), 'radius2:', r)
    x = int(x)
    y = int(y)
    r = int(r)
    cv2.circle(result, (x, y), r, (0, 0, 255), 1)

# save results
cv2.imwrite('2circles_thresh.png', thresh)
cv2.imwrite('2circles_thresh2.png', thresh2)
cv2.imwrite('2circles_morph.png', morph)
cv2.imwrite('2circles_circles.png', result)

# show images
cv2.imshow('thresh', thresh)
cv2.imshow('thresh2', thresh2)
cv2.imshow('morph', morph)
cv2.imshow('result', result)
cv2.waitKey(0)
cv2.destroyAllWindows()

Threshold for inner circle:

Threshold for outer circle:

Morphology image for outer circle:

Resulting 2 circles:

Data For 2 Circles:

center1: (589.5, 559.5) radius1: 91.4
center2: (594.5, 549.5) radius2: 225.7