python - 求重叠椭圆的直径和面积(OpenCV，Python)

Question

这是我在 Stackoverflow 上的第一个问题。我有点激动，如果说错了请原谅。我们从颜料中随机绘制有重叠和没有重叠的混合椭圆。我正在分享我正在处理的图像和我的代码。我不是 opencv 模块的专业人士，我编写的代码是受资料来源启发的研究结果。

我的代码的目的是，

使用 cv2.fitEllipse 方法检测随机绘制的有或没有重叠椭圆。接下来，找到检测到的椭圆的长轴、短轴和面积。

我的代码的问题实际上是这样的，

在重叠椭圆中，正常情况下拟合椭圆时，应该拟合2个椭圆，但拟合了大约6-7个椭圆，无法达到我想要计算的值。

我愿意接受您的帮助，提前谢谢您。

示例图片:

import cv2
import numpy as np
import random as rng
import math

img = cv2.imread('overlapping_ellipses.png', 1)
imge= cv2.cvtColor(img,cv2.COLOR_RGB2BGR)
gray = cv2.cvtColor(imge, cv2.COLOR_BGR2GRAY)
blur = cv2.blur(gray, (2,2), 3)
edged = cv2.Canny(blur, 50, 100)
kernel= np.ones((2,2))
edged1 = cv2.dilate(edged, kernel, iterations=2)
edged2 = cv2.erode(edged1, kernel, iterations=2)

def thresh_callback(val):
 threshold = val

 canny_output = cv2.Canny(edged2, threshold, threshold * 4)
 contours, _ = cv2.findContours(canny_output, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
 minRect = [None]*len(contours)
 minEllipse = [None]*len(contours)
 for i, c in enumerate(contours):
    minRect[i] = cv2.minAreaRect(c)
    if c.shape[0] > 5:
        minEllipse[i] = cv2.fitEllipse(c)
        (x1,y1),(d1,d2),angle = minEllipse[i]
        print('\nX1: ', round(x1,4), '\nY1: ', round(y1,4), '\nD1:',round(d1,4), '\nD2',round(d2,4), '\nAngle:', round(angle,4))
        long= x1-d2
        small= y1-d1
        major= long/2
        minor= small/2
        pixel= 37.795275591
        major1= major/pixel
        minor1= minor/pixel
        print('--------------------------------')
        print('Major axis is: ', abs(round(major1,4)), 'cm')
        print('Minor axis is: ', abs(round(minor1,4)), 'cm')
        print('--------------------------------')
drawing = np.zeros((canny_output.shape[1], canny_output.shape[1], 3), dtype=np.uint8)

for i, c in enumerate(contours):
    color = (rng.randint(0,256), rng.randint(0,256), rng.randint(0,256))
    cv2.drawContours(drawing, contours, i, color)
    if c.shape[0] > 5:
        cv2.ellipse(drawing, minEllipse[i], color, 1)
        
cv2.imshow('Fitting Ellips', drawing)

source_window = 'Source'
cv2.namedWindow(source_window)
cv2.imshow(source_window, img)
max_thresh = 255
thresh = 100
cv2.createTrackbar('Canny Thresh:', source_window, thresh, max_thresh, thresh_callback)
thresh_callback(thresh)
cv2.waitKey()

Answer 1

第 1 步:识别并分离输入图像中的 Blob 。

由于我们在这里不关心颜色信息，因此我们可以直接将图像加载为灰度。

image = cv2.imread('input.png', cv2.IMREAD_GRAYSCALE)

输入图像包含白色背景上的黑色椭圆。 我们只需要 Blob 的外部轮廓，并且 cv2.findContours期望黑色背景上出现白色 Blob 。 因此我们需要反转图像。同时我们需要一个二值图像。我们可以使用cv2.threshold完成这两项任务。

一旦我们检测到 Blob 轮廓，我们就可以将每个 Blob 的一些有用信息收集到一个简单的基于 map 的数据结构中。

def detect_blobs(image):
    _,img_binary = cv2.threshold(image, 127, 255, cv2.THRESH_BINARY_INV)
    contours, _ = cv2.findContours(img_binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    blobs = []
    for i, contour in enumerate(contours):
        orig_x, orig_y, width, height = cv2.boundingRect(contour)
        roi_image = image[orig_y:orig_y+height,orig_x:orig_x+width]
        blobs.append({
            "i" : i
            , "contour" : contour
            , "origin" : (orig_x, orig_y)
            , "size" : (width, height)
            , "roi_image" : roi_image
        })
    return blobs

---

第 2 步:处理每个 blob

首先我们需要确定 Blob 是否是单个椭圆，或者是否是相交椭圆上的一对。 实现此目的的一种方法是寻找凸性缺陷。

由于轮廓的坐标是用整数表示的，因此即使是单椭圆场景也会表现出一些凸性缺陷。 然而，它们的大小(轮廓上最远点与封闭凸包段之间的距离)将非常小，通常低于 1 像素。 另一方面，一对相交椭圆的轮廓将具有较大的凸性缺陷，曲线相交的四个点中的每一个都有一个。

这种区别可以在以下两张图像上看到(轮廓为蓝色，凸包为红色，识别出的大凸面缺陷的交点/位置为橙色圆圈):

<表类=“s-表”> <标题> <正文> 单个椭圆 两个相交的椭圆

因此，我们过滤掉任何小的凸性缺陷，并记下大的凸性缺陷的位置。现在我们只剩下 3 种可能的情况。

---

场景 A:未检测到交叉点

仅发现了小的凸度缺陷，这意味着这很可能是单个椭圆。我们只需将椭圆拟合到轮廓上即可继续。

场景 B:恰好检测到 4 个交叉点

在本例中，我们有 2 个相交的椭圆。我们使用交点将轮廓分成 4 段，每个段对应 Blob 的每个“叶”。每个线段都应包含界定它的两个交点。

在下图中，线段以绿色、黄色、青色和洋红色显示，交点为橙色圆圈:

现在，我们可以组合彼此相对的线段对(即绿色+青色和黄色+洋红色)来获得两个点列表，每个椭圆对应一个点列表。同样，我们只是简单地为每个点列表拟合一个椭圆。

场景 C:检测到其他数量的交叉点

这被认为是无效的情况。

def process_blob(blob):
    MAJOR_DEFECT_THRESHOLD = 2.0
    
    contour = blob["contour"]
    blob["hull"] = cv2.convexHull(contour)
    
    hull_idx = cv2.convexHull(contour, returnPoints=False)
    defects = cv2.convexityDefects(contour, hull_idx)
    
    intersections = []
    for i,defect in enumerate(np.squeeze(defects, 1)):
        _, _, far_idx, far_dist = defect
        real_far_dist = far_dist / 256.0
        if real_far_dist >= MAJOR_DEFECT_THRESHOLD:
            intersections.append(far_idx)
    
    if len(intersections) == 0:
        print("One ellipse")
        blob["ellipses"] = [cv2.fitEllipse(contour)]
    elif len(intersections) == 4:
        print("Two ellipses")
        blob["segments"] = [
            contour[intersections[0]:intersections[1]+1]
            , contour[intersections[1]:intersections[2]+1]
            , contour[intersections[2]:intersections[3]+1]
            , np.vstack([contour[intersections[3]:],contour[:intersections[0]+1]])
        ]
        split_contours = [
            np.vstack([blob["segments"][0], blob["segments"][2]])
            , np.vstack([blob["segments"][1], blob["segments"][3]])
        ]
        blob["ellipses"] = [cv2.fitEllipse(c) for c in split_contours]
    else:
        print("Invalid scenario")
        blob["ellipses"] = []
        
    return blob["ellipses"]

---

此时，计算所需的参数很简单 - 我将把它作为练习留给读者。

作为奖励，这里有一些用于调试目的的简单可视化:

def visualize_blob(blob):
    PADDING = 20
    
    orig_x, orig_y = blob["origin"]
    offset = (orig_x - PADDING, orig_y - PADDING)
    
    input_img = cv2.copyMakeBorder(blob["roi_image"]
        , PADDING, PADDING, PADDING, PADDING
        , cv2.BORDER_CONSTANT, None, 255)

    adjusted_img = cv2.add(input_img, 127) - 63
    output_img_ch = cv2.cvtColor(adjusted_img, cv2.COLOR_GRAY2BGR)
    output_img_seg = output_img_ch.copy()
    output_img_el = output_img_ch.copy()
    
    cv2.drawContours(output_img_ch, [blob["hull"] - offset], 0, (127,127,255), 4)
    cv2.drawContours(output_img_ch, [blob["contour"] - offset], 0, (255,127,127), 2)
    
    SEGMENT_COLORS = [(0,255,0),(0,255,255),(255,255,0),(255,0,255)]
    if "segments" in blob:
        for i in range(4):
            cv2.polylines(output_img_seg, [blob["segments"][i] - offset], False, SEGMENT_COLORS[i], 4)
        for i in range(4):
            center = (blob["segments"][i] - offset)[0][0]
            cv2.circle(output_img_ch, center, 4, (0,191,255), -1)
            cv2.circle(output_img_seg, center, 4, (0,191,255), -1)
        
    
    for ellipse in blob["ellipses"]:
        offset_ellipse = ((ellipse[0][0] - offset[0], ellipse[0][1] - offset[1]), ellipse[1], ellipse[2])
        cv2.ellipse(output_img_el, offset_ellipse, (0,0,255), 2)
    
    cv2.imshow('', np.hstack([output_img_ch,output_img_seg, output_img_el]))
    cv2.imwrite('output_%d_ch.png' % blob["i"], output_img_ch)
    cv2.imwrite('output_%d_seg.png' % blob["i"], output_img_seg)
    cv2.imwrite('output_%d_el.png' % blob["i"], output_img_el)
    cv2.waitKey()

---

把它们放在一起:

import cv2
import numpy as np

## INSERT THE FUNCTIONS LISTED ABOVE IN THE QUESTION ##

image = cv2.imread('input.png', cv2.IMREAD_GRAYSCALE)

blobs = detect_blobs(image)
print("Found %d blob(s)." % len(blobs))

for blob in blobs:
    process_blob(blob)
    visualize_blob(blob)