How to watershed two connected circles

Question

There are the two images I am working with:

(Here the images are not the same size, but in my programme they are the same)

After taking the skimage.metrics.structural_similarity() of the two images above, I have the following thresh:

As you can see, it consists of 2 shapes which are almost circles but not quite (the extra part at the bottom right is the circle's shadow)

I want to watershed this thresh so that I obtain two circles, but my current code gives me this:

Instead, I want something that looks like this in blue:

# import the necessary packages
from skimage.feature import peak_local_max
from skimage.segmentation import watershed
from scipy import ndimage
import numpy as np
import cv2
from skimage.metrics import structural_similarity

imageA = cv2.imread("frames/thing150.png") #the left image
imageB = cv2.imread("frames/thing180.png") #the right image
grayA = cv2.cvtColor(imageA, cv2.COLOR_BGR2GRAY)
grayB = cv2.cvtColor(imageB, cv2.COLOR_BGR2GRAY)

(score, diff) = structural_similarity(grayA, grayB, full=True)
diff = (diff * 255).astype("uint8")

thresh = cv2.threshold(diff, 0, 255, cv2.THRESH_BINARY_INV | cv2.THRESH_OTSU)[1]
# cv2.namedWindow("Thresh", cv2.WINDOW_NORMAL)
# cv2.imshow("Thresh", thresh)

# compute the exact Euclidean distance from every binary
# pixel to the nearest zero pixel, then find peaks in this
# distance map
D = ndimage.distance_transform_edt(thresh)
localMax = peak_local_max(D, indices=False, min_distance=100, labels=thresh)
# perform a connected component analysis on the local peaks,
# using 8-connectivity, then appy the Watershed algorithm
markers = ndimage.label(localMax, structure=np.ones((3, 3)))[0]
labels = watershed(-D, markers, mask=thresh)
print(f"[INFO] {len(np.unique(labels)) - 1} unique segments found")

# loop over the unique labels returned by the Watershed
# algorithm
for label in np.unique(labels):
    # if the label is zero, we are examining the 'background'
    # so simply ignore it
    if label == 0:
        continue
    # otherwise, allocate memory for the label region and draw
    # it on the mask
    mask = np.zeros(grayB.shape, dtype="uint8")
    mask[labels == label] = 255
    # detect contours in the mask and grab the largest one
    (cnts, _) = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    c = max(cnts, key=cv2.contourArea)
    # draw a circle enclosing the object
    ((x, y), r) = cv2.minEnclosingCircle(c)
    cv2.circle(imageB, (int(x), int(y)), int(r), (0, 255, 0), 2)
    cv2.putText(imageB, "#{}".format(label), (int(x) - 10, int(y)),
        cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 0, 255), 2)
    print(len(cnts))
print("----")
print(np.unique(labels))
# show the output imageB
cv2.namedWindow("Output", cv2.WINDOW_NORMAL)
cv2.imshow("Output", imageB)
cv2.waitKey(0)

I am not familiar with watershed, so I copied the code from this website. I tried varying the parameter for min_distance in localMax = peak_local_max(D, indices=False, min_distance=100, labels=thresh), but it didn't solve my problem.

I also tried using OpenCV's watershed algorithm, but for some reason it didn't work. If that's better than skimage's, then I will try that.

Any suggestions would be appreciated.

P.S. Is it thresh = cv2.threshold(diff, 0, 255, cv2.THRESH_BINARY_INV | cv2.THRESH_OTSU)[1]

or thresh = cv2.threshold(diff, 0, 255, cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)?

I have seen various sources using both and both of them work for me, are there any differences?

Answer 1

This solution is gonna start from the thresh image you uploaded:

1. Fill contours: Since the thresh image is not really connected (there are a lot of black spaces in between the white spaces), we have to find a way to fill those "holes". We could use some opening kernels, but the size will tend to vary inbetween images. I figured it would be better to find the external contours and create some sort of mask by filling them.

contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
# get largest contour (snowman figurine)
idx = np.argmax([cv2.contourArea(cnt) for cnt in contours])
cnt = contours[idx]
# create mask
thresh = cv2.fillPoly(thresh, np.int32([cnt]), 0, (255, 255, 255)) #Bug with fillPoly, needs explict cast to 32bit

2. Watershed algorithm: Use the watershed algorithm implementation in order to separate the two connected objects. I did some minor modifications since the tutorial is not as updated as one would hope to.

kernel = np.ones((3, 3), np.uint8)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)
background = cv2.dilate(opening, kernel, iterations=3)

# extract foreground
dst = cv2.distanceTransform(opening, cv2.DIST_L2, 5, dstType=cv2.CV_32F)
_, foreground = cv2.threshold(dst, 0.7 * dst.max(), 255, cv2.THRESH_BINARY)
foreground = np.uint8(foreground)
unknown = cv2.subtract(background, foreground)

# get markers
_, markers = cv2.connectedComponents(foreground)
markers += 1
markers[unknown == 255] = 0

# Since the original image is in fact two different images, use
# instead the thresh image, which is the combination of both.
thresh = cv2.cvtColor(thresh, cv2.COLOR_GRAY2BGR)
markers = cv2.watershed(thresh, markers)
# normalize markers so they can be visualized
markers = cv2.normalize(
  markers, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U
)

3. Threshold image (again) + contours: Now that we finally separated the contours, we can now find the enclosing circles. In order to do this, I would recommend to threshold the markers so we don't get more circles that we intend to.

NOTE: I used the inverse binary so as to get that nice middle line in between the two contours, otherwise we would be back to square one.

_, thresh = cv2.threshold(markers, 150, 255, cv2.THRESH_BINARY_INV)
contours, _ = cv2.findContours(thresh, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)

4. Approximate circles: Now we can find the minimum enclosing circles.

NOTE: I used a filter for the radius since the number of fetched contours might be more than two; you might need to use this filter as well just to make sure you always get the circles you are looking for and not bigger (or smaller) ones.

circles = []
for cnt in contours:
  (cx, cy), radius = cv2.minEnclosingCircle(cnt)
  if radius < 100:
    circles.append((int(cx), int(cy), int(radius)))

This procedure can probably be optimized/reduced since I improvised a little bit and didn't put much thought into it. But anyway, I hope this can be of help :)

PS: It doesn't matter if you use either "+" or "|" in opencv, they mean the same thing. Although I'd recommend you use "+" since it is more readable, but it is completely up to you.