Python基于随机采样一至性实现拟合椭圆

Python
256
0
0
2023-06-16
标签   Python实践

检测这些圆,先找轮廓后通过轮廓点拟合椭圆

import cv
import numpy as np
import matplotlib.pyplot as plt
import math
from Ransac_Process import RANSAC
 
def lj_img(img):
    wlj, hlj = img.shape[], img.shape[0]
    lj_dis =      # 连接白色区域的判定距离
    for ilj in range(wlj):
        for jlj in range(hlj):
            if img[jlj, ilj] ==:    # 判断上下左右是否存在白色区域并连通
                for im in range(, lj_dis):
                    for jm in range(, lj_dis):
                        if ilj - im >= and jlj - jm >= 0  and img[jlj - jm, ilj - im] == 255:
                            cv.line(img, (jlj, ilj), (jlj - jm, ilj - im), (255, 255, 255), thickness=1)
                        if ilj + im < wlj and jlj + jm < hlj and img[jlj + jm, ilj + im] ==:
                            cv.line(img, (jlj, ilj), (jlj + jm, ilj + im), (255, 255, 255), thickness=1)
    return img
 
def cul_area(x_mask, y_mask, r_circle, mask):
    mask_label = mask.copy()
    num_area =
    for xm in range(x_mask+r_circle-, x_mask+r_circle+10):
        for ym in range(y_mask+r_circle-, y_mask+r_circle+10):
            # print(mask[ym, xm])
            if (pow((xm-x_mask),) + pow((ym-y_mask), 2) - pow(r_circle,  2)) == 0 and mask[ym, xm][0] == 255:
                num_area +=
                mask_label[ym, xm] = (, 0, 255)
    cv.imwrite('./test2/mask_label.png', mask_label)
    print(num_area)
    return num_area
 
def mainFigure(img, point):
    # params = cv.SimpleBlobDetector_Params()  # 黑色斑点面积大小:1524--1581--1400--周围干扰面积: 1325--1695--1688--
    # # Filter by Area.   设置斑点检测的参数
    # params.filterByArea = True  # 根据大小进行筛选
    # params.minArea =e2
    # params.maxArea =e4
    # params.minDistBetweenBlobs =  # 设置两个斑点间的最小距离 10*7.5
    # # params.filterByColor = True             # 跟据颜色进行检测
    # params.filterByConvexity = False  # 根据凸性进行检测
    # params.minThreshold =  # 二值化的起末阈值,只有灰度值大于当前阈值的值才会被当成特征值
    # params.maxThreshold = * 2.5  # 75
    # params.filterByColor = True  # 检测颜色限制,黑色,255白色
    # params.blobColor =
    # params.filterByCircularity = True
    # params.minCircularity =.3
 
    point_center = []
    # cv.imwrite('./test2/img_source.png', img)
    img_hsv = cv.cvtColor(img, cv2.COLOR_BGR2HSV)
    # cv.imwrite('./test2/img_hsv.png', img_hsv)
    w, h = img.shape[], img.shape[0]
    w_hsv, h_hsv = img_hsv.shape[], img_hsv.shape[0]
    for i_hsv in range(w_hsv):
        for j_hsv in range(h_hsv):
            if img_hsv[j_hsv, i_hsv][] < 200 and img_hsv[j_hsv, i_hsv][1] < 130 and img_hsv[j_hsv, i_hsv][2] > 120:
                # if hsv[j_hsv, i_hsv][] < 100 and hsv[j_hsv, i_hsv][1] < 200 and hsv[j_hsv, i_hsv][2] > 80:
                img_hsv[j_hsv, i_hsv] =, 255, 255
            else:
                img_hsv[j_hsv, i_hsv] =, 0, 0
    # cv.imwrite('./test2/img_hsvhb.png', img_hsv)
    # cv.imshow("hsv", img_hsv)
    # cv.waitKey()
 
    # 灰度化处理图像
    grayImage = cv.cvtColor(img_hsv, cv2.COLOR_BGR2GRAY)
    # mask = np.zeros((grayImage.shape[], grayImage.shape[1]), np.uint8)
    # mask = cv.cvtColor(mask, cv2.COLOR_GRAY2BGR)
    # cv.imwrite('./mask.png', mask)
 
    # 尝试寻找轮廓
    contours, hierarchy = cv.findContours(grayImage, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
 
    # 合并轮廓
    if len(contours) >:
        # print(contours)
        # 去掉离图中心最远的圆
        max_idex, dis_max =, 0
        for c_i in range(len(contours)):
            c = contours[c_i]
            cx, cy, cw, ch = cv.boundingRect(c)
            dis = math.sqrt(pow((cx + cw / - w / 2), 2) + pow((cy + ch / 2 - h / 2), 2))
            if dis > dis_max:
                dis_max = dis
                max_idex = c_i
        contours.pop(max_idex)
        # print(contours)
 
        if len(contours) >:
            contours_merge = np.vstack([contours[], contours[1]])
            for i in range(, len(contours)):
                contours_merge = np.vstack([contours_merge, contours[i]])
            cv.drawContours(img, contours_merge, -1, (0, 255, 255), 1)
            cv.imwrite('./test2/img_res.png', img)
            # cv.imshow("contours_merge", img)
            # cv.waitKey()
        else:
            contours_merge = contours[]
    else:
        contours_merge = contours[]
 
    # RANSAC拟合
    points_data = np.reshape(contours_merge, (-, 2))  # ellipse edge points set
 
    print("points_data", len(points_data))
    #.Ransac fit ellipse param
    Ransac = RANSAC(data=points_data, threshold=.5, P=.99, S=.5, N=20)
    # Ransac = RANSAC(data=points_data, threshold=.05, P=.99, S=.618, N=25)
 
    (X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()
    # print( (X, Y), (LAxis, SAxis))
    # 拟合圆
    cv.ellipse(img, ((X, Y), (LAxis, SAxis), Angle), (0, 0, 255), 1, cv2.LINE_AA)  # 画圆
    cv.circle(img, (int(X), int(Y)), 3, (0, 0, 255), -1)  # 画圆心
    point_center.append(int(X))
    point_center.append(int(Y))
 
    rrt = cv.fitEllipse(contours_merge)  # x, y)代表椭圆中心点的位置(a, b)代表长短轴长度,应注意a、b为长短轴的直径,而非半径,angle 代表了中心旋转的角度
    # print("rrt", rrt)
    cv.ellipse(img, rrt, (255, 0, 0), 1, cv2.LINE_AA)  # 画圆
    x, y = rrt[]
    cv.circle(img, (int(x), int(y)), 3, (255, 0, 0), -1)  # 画圆心
    point_center.append(int(x))
    point_center.append(int(y))
    # print("no",(x,y))
 
    cv.imshow("fit circle", img)
    cv.waitKey()
    # cv.imwrite("./test2/fitcircle.png", img)
 
    # # 尝试寻找轮廓
    # contours, hierarchy = cv.findContours(grayImage, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    # # print('初次检测数量: ', len(contours))
    # if len(contours) ==:
    #     cv.drawContours(mask, contours[0], -1, (255, 255, 255), 1)
    #     cv.imwrite('./mask.png', mask)
    #     x, y, w, h = cv.boundingRect(contours[0])
    #     cv.circle(img, (int(x+w/2), int(y+h/2)), 1, (0, 0, 255), -1)
    #     cv.rectangle(img, (x, y), (x + w + 1, y + h + 1), (0, 255, 255), 1)
    #     point_center.append(x + w / + point0[0])
    #     point_center.append(y + h / + point0[1])
    #     cv.imwrite('./center1.png', img)
    # else:
    #     # 去除小面积杂点, 连接轮廓,求最小包围框
    #     kernel = np.ones((3, 3), dtype=np.uint8)
    #     kernel = np.ones((2, 2), dtype=np.uint8)
    #     grayImage = cv.dilate(grayImage, kernel1, 1)  # 1:迭代次数,也就是执行几次膨胀操作
    #     grayImage = cv.erode(grayImage, kernel2, 1)
    #     cv.imwrite('./img_dilate_erode.png', grayImage)
    #     contours, hierarchy = cv.findContours(grayImage, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    #     if len(contours) ==:
    #         cv.drawContours(mask, contours[0], -1, (255, 255, 255), 1)
    #         cv.imwrite('./mask.png', mask)
    #         x, y, w, h = cv.boundingRect(contours[0])
    #         cv.circle(img, (int(x + w / 2), int(y + h / 2)), 1, (0, 0, 255), -1)
    #         cv.rectangle(img, (x, y), (x + w + 1, y + h + 1), (0, 255, 255), 1)
    #         point_center.append(x + w / + point0[0])
    #         point_center.append(y + h / + point0[1])
    #         cv.imwrite('./center1.png', img)
    #     else:
    #         gray_circles = cv.HoughCircles(grayImage, cv2.HOUGH_GRADIENT, 4, 10000, param1=100, param2=81, minRadius=10, maxRadius=19)
    #         # cv.imwrite('./img_gray_circles.jpg', gray_circles)
    #         if len(gray_circles[]) > 0:
    #             print('霍夫圆个数:', len(gray_circles[]))
    #             for (x, y, r) in gray_circles[]:
    #                 x = int(x)
    #                 y = int(y)
    #                 cv.circle(grayImage, (x, y), int(r), (255, 255, 255), -1)
    #             cv.imwrite('./img_hf.jpg', grayImage)
    #
    #             detector = cv.SimpleBlobDetector_create(params)
    #             keypoints = list(detector.detect(grayImage))
    #             for poi in keypoints:  # 回归到原大图坐标系
    #                 x_poi, y_poi = poi.pt[], poi.pt[1]
    #                 cv.circle(img, (int(x_poi), int(y_poi)), 20, (255, 255, 255), -1)
    #                 point_center.append(poi.pt[] + point0[0])
    #                 point_center.append(poi.pt[] + point0[1])
    #                 cv.imwrite('./img_blob.png', img)
    #         else:
    #             for num_cont in range(len(contours)):
    #                 cont = cv.contourArea(contours[num_cont])
    #                 # if cont >:
    #                 #     contours.append(contours[num_cont])
    #                 if cont <=:
    #                     x, y, w, h = cv.boundingRect(contours[num_cont])
    #                     cv.rectangle(grayImage, (x, y), (x + w, y + h), (0, 0, 0), -1)
    #             cv.imwrite('./img_weilj.png', grayImage)
    #             grayImage = lj_img(grayImage)
    #             cv.imwrite('./img_lj.png', grayImage)
    #             contours, hierarchy = cv.findContours(grayImage, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    #             # print('再次检测数量: ', len(contours))
    #
    #             cv.drawContours(mask, contours[0], -1, (255, 255, 255), 1)
    #             cv.imwrite('./mask.png', mask)
    #             x, y, w, h = cv.boundingRect(contours[0])
    #             cv.circle(img, (int(x + w / 2), int(y + h / 2)), 1, (0, 0, 255), -1)
    #             cv.rectangle(img, (x, y), (x + w + 1, y + h + 1), (0, 255, 255), 1)
    #             point_center.append(x + w / + point0[0])
    #             point_center.append(y + h / + point0[1])
    #             cv.imwrite('./center1.png', img)
 
    return point_center[], point_center[1]
 
if __name__ == "__main__":
 
    for i in range(,6):
        imageName = "s"
        imageName += str(i)
        path = './Images/danHoles/' + imageName + '.png'
        print(path)
        img = cv.imread(path)
        point = [0, 0]
        cir_x, cir_y = mainFigure(img, point)
 
    # img = cv.imread('./Images/danHoles/s2.png')
    # point = [0, 0]
    # cir_x, cir_y = mainFigure(img, point)

Ransac_Process.py

import cv
import math
import random
import numpy as np
from numpy.linalg import inv, svd, det
import time
 
class RANSAC:
    def __init__(self, data, threshold, P, S, N):
        self.point_data = data  # 椭圆轮廓点集
        self.length = len(self.point_data)  # 椭圆轮廓点集长度
        self.error_threshold = threshold  # 模型评估误差容忍阀值
 
        self.N = N  # 随机采样数
        self.S = S  # 设定的内点比例
        self.P = P  # 采得N点去计算的正确模型概率
        self.max_inliers = self.length * self.S  # 设定最大内点阀值
        self.items =
 
        self.count =  # 内点计数器
        self.best_model = ((, 0), (1e-6, 1e-6), 0)  # 椭圆模型存储器
 
    def random_sampling(self, n):
        # 这个部分有修改的空间,这样循环次数太多了,可以看看别人改进的ransac拟合椭圆的论文
        """随机取n个数据点"""
        all_point = self.point_data
        select_point = np.asarray(random.sample(list(all_point), n))
        return select_point
 
    def GeometricConic(self, ellipse):
        # 这个部分参考了GitHub中的一位大佬的,但是时间太久,忘记哪个人的了
        """计算椭圆方程系数"""
        # Ax ^ + Bxy + Cy ^ 2 + Dx + Ey + F
        (x, y0), (bb, aa), phi_b_deg = ellipse
 
        a, b = aa /, bb / 2  # Semimajor and semiminor axes
        phi_b_rad = phi_b_deg * np.pi /.0  # Convert phi_b from deg to rad
        ax, ay = -np.sin(phi_b_rad), np.cos(phi_b_rad)  # Major axis unit vector
 
        # Useful intermediates
        a = a * a
        b = b * b
 
        # Conic parameters
        if a > 0 and b2 > 0:
            A = ax * ax / a + ay * ay / b2
            B = * ax * ay / a2 - 2 * ax * ay / b2
            C = ay * ay / a + ax * ax / b2
            D = (- * ax * ay * y0 - 2 * ax * ax * x0) / a2 + (2 * ax * ay * y0 - 2 * ay * ay * x0) / b2
            E = (- * ax * ay * x0 - 2 * ay * ay * y0) / a2 + (2 * ax * ay * x0 - 2 * ax * ax * y0) / b2
            F = ( * ax * ay * x0 * y0 + ax * ax * x0 * x0 + ay * ay * y0 * y0) / a2 + \
                (- * ax * ay * x0 * y0 + ay * ay * x0 * x0 + ax * ax * y0 * y0) / b2 - 1
        else:
            # Tiny dummy circle - response to a or b2 == 0 overflow warnings
            A, B, C, D, E, F = (, 0, 1, 0, 0, -1e-6)
 
        # Compose conic parameter array
        conic = np.array((A, B, C, D, E, F))
        return conic
 
    def eval_model(self, ellipse):
        # 这个地方也有很大修改空间,判断是否内点的条件在很多改进的ransac论文中有说明,可以多看点论文
        """评估椭圆模型,统计内点个数"""
        # this an ellipse ?
        a, b, c, d, e, f = self.GeometricConic(ellipse)
        E = * a * c - b * b
        if E <=:
            # print('this is not an ellipse')
            return, 0
 
        #  which long axis ?
        (x, y), (LAxis, SAxis), Angle = ellipse
        LAxis, SAxis = LAxis /, SAxis / 2
        if SAxis > LAxis:
            temp = SAxis
            SAxis = LAxis
            LAxis = temp
 
        # calculate focus
        Axis = math.sqrt(LAxis * LAxis - SAxis * SAxis)
        f_x = x - Axis * math.cos(Angle * math.pi / 180)
        f_y = y - Axis * math.sin(Angle * math.pi / 180)
        f_x = x + Axis * math.cos(Angle * math.pi / 180)
        f_y = y + Axis * math.sin(Angle * math.pi / 180)
 
        # identify inliers points
        f, f2 = np.array([f1_x, f1_y]), np.array([f2_x, f2_y])
        f_distance = np.square(self.point_data - f1)
        f_distance = np.square(self.point_data - f2)
        all_distance = np.sqrt(f_distance[:, 0] + f1_distance[:, 1]) + np.sqrt(f2_distance[:, 0] + f2_distance[:, 1])
 
        Z = np.abs( * LAxis - all_distance)
        delta = math.sqrt(np.sum((Z - np.mean(Z)) **) / len(Z))
 
        # Update inliers set
        inliers = np.nonzero(Z <.8 * delta)[0]
        inlier_pnts = self.point_data[inliers]
 
        return len(inlier_pnts), inlier_pnts
 
    def execute_ransac(self):
        Time_start = time.time()
        while math.ceil(self.items):
            # print(self.max_inliers)
 
            #.select N points at random
            select_points = self.random_sampling(self.N)
 
            #.fitting N ellipse points
            ellipse = cv.fitEllipse(select_points)
 
            #.assess model and calculate inliers points
            inliers_count, inliers_set = self.eval_model(ellipse)
 
            #.number of new inliers points more than number of old inliers points ?
            if inliers_count > self.count:
                ellipse_ = cv.fitEllipse(inliers_set)  # fitting ellipse for inliers points
                self.count = inliers_count  # Update inliers set
                self.best_model = ellipse_  # Update best ellipse
                # print("self.count", self.count)
 
                #.number of inliers points reach the expected value
                if self.count > self.max_inliers:
                    print('the number of inliers: ', self.count)
                    break
 
                # Update items
                # print(math.log( - pow(inliers_count / self.length, self.N)))
                self.items = math.log( - self.P) / math.log(1 - pow(inliers_count / self.length, self.N))
 
        return self.best_model
 
 
if __name__ == '__main__':
 
 
    #.find ellipse edge line
    contours, hierarchy = cv.findContours(grayImage, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)
 
    #.Ransac fit ellipse param
    points_data = np.reshape(contours, (-, 2))  # ellipse edge points set
    Ransac = RANSAC(data=points_data, threshold=.5, P=.99, S=.618, N=10)
    (X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()

检测对象

检测结果

蓝色是直接椭圆拟合的结果

红色是Ransc之后的结果