Ransac拟合椭圆

一、Ransac算法介绍

RANSAC(RAndom SAmple Consensus,随机采样一致)最早是由Fischler和Bolles在SRI上提出用来解决LDP(Location Determination Proble)问题的，该算法是从一组含有“外点”(outliers)的数据中正确估计数学模型参数的迭代算法。“外点”一般指的的数据中的噪声点，比如说匹配中的误匹配和估计曲线中的离群点。除去“外点”后剩下的数据点则为“内点”，即符合拟合模型的数据点。RANSAC算法的主要处理思想就是通过不断的迭代数据点，通过决断条件剔除“外点”，筛选出“内点”，最后将筛选出的全部内点作为好的点集去重新拟合去估计匹配模型。但RANSAC算法具有不稳定性，它只能在一种概率下产生结果，并且这个概率会随着迭代次数的增加而加大。

二、Ransac算法实现示例

这一部分基本上网上代码及讲解都非常多，就不过多讲解了。
附上我在知乎上看到的一篇写的比较完整的一个拟合直线模型示例：

https://zhuanlan.zhihu.com/p/62238520

三、基于pyhton的Ransac拟合椭圆实现

1.Ransac拟合椭圆流程：
①. 找到可能是椭圆的轮廓线，可以使用cv2.findContours去找，返回的轮廓线数据点是[[x1, y1], [x2, y2], … , [xn, yn]]格式的数据
②. 随机取5个数据点进行椭圆方程拟合，并求出椭圆方程参数。椭圆一般方程：Ax ^ 2 + Bxy + Cy ^ 2 + Dx + Ey + F = 0
③. 利用求得的椭圆模型进行数据点评估，判断条件一般为：①代数距离二范数；②几何距离二范数，判断出符合该椭圆模型的内点集
④. 判断3求出的内点集数量是否大于上一次内点集数量，如果是，则更新最优内点集及最优椭圆模型
⑤. 判断最优内点集是否达到预期数量，如果是，则退出循环，输出最优椭圆模型
⑥. 循环2-5步骤
⑦. 完成椭圆拟合

2.python拟合椭圆模型代码：

	import cv2import mathimport randomimport numpy as npfrom numpy.linalg import inv, svd, detclass 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 = 999self.count = 0  # 内点计数器self.best_model = ((0, 0), (1e-6, 1e-6), 0)  # 椭圆模型存储器def random_sampling(self, n):# 这个部分有修改的空间，这样循环次数太多了，可以看看别人改进的ransac拟合椭圆的论文"""随机取n个数据点"""all_point = self.point_dataselect_point = np.asarray(random.sample(list(all_point), n))return select_pointdef Geometric2Conic(self, ellipse):# 这个部分参考了GitHub中的一位大佬的，但是时间太久，忘记哪个人的了"""计算椭圆方程系数"""# Ax ^ 2 + Bxy + Cy ^ 2 + Dx + Ey + F(x0, y0), (bb, aa), phi_b_deg = ellipsea, b = aa / 2, bb / 2  # Semimajor and semiminor axesphi_b_rad = phi_b_deg * np.pi / 180.0  # Convert phi_b from deg to radax, ay = -np.sin(phi_b_rad), np.cos(phi_b_rad)  # Major axis unit vector# Useful intermediatesa2 = a * ab2 = b * b# Conic parametersif a2 > 0 and b2 > 0:A = ax * ax / a2 + ay * ay / b2B = 2 * ax * ay / a2 - 2 * ax * ay / b2C = ay * ay / a2 + ax * ax / b2D = (-2 * ax * ay * y0 - 2 * ax * ax * x0) / a2 + (2 * ax * ay * y0 - 2 * ay * ay * x0) / b2E = (-2 * ax * ay * x0 - 2 * ay * ay * y0) / a2 + (2 * ax * ay * x0 - 2 * ax * ax * y0) / b2F = (2 * ax * ay * x0 * y0 + ax * ax * x0 * x0 + ay * ay * y0 * y0) / a2 + \(-2 * ax * ay * x0 * y0 + ay * ay * x0 * x0 + ax * ax * y0 * y0) / b2 - 1else:# Tiny dummy circle - response to a2 or b2 == 0 overflow warningsA, B, C, D, E, F = (1, 0, 1, 0, 0, -1e-6)# Compose conic parameter arrayconic = np.array((A, B, C, D, E, F))return conicdef eval_model(self, ellipse):# 这个地方也有很大修改空间，判断是否内点的条件在很多改进的ransac论文中有说明，可以多看点论文"""评估椭圆模型，统计内点个数"""# this an ellipse ?a, b, c, d, e, f = self.Geometric2Conic(ellipse)E = 4 * a * c - b * bif E <= 0:# print('this is not an ellipse')return 0, 0#  which long axis ?(x, y), (LAxis, SAxis), Angle = ellipseLAxis, SAxis = LAxis / 2, SAxis / 2if SAxis > LAxis:temp = SAxisSAxis = LAxisLAxis = temp# calculate focusAxis = math.sqrt(LAxis * LAxis - SAxis * SAxis)f1_x = x - Axis * math.cos(Angle * math.pi / 180)f1_y = y - Axis * math.sin(Angle * math.pi / 180)f2_x = x + Axis * math.cos(Angle * math.pi / 180)f2_y = y + Axis * math.sin(Angle * math.pi / 180)# identify inliers pointsf1, f2 = np.array([f1_x, f1_y]), np.array([f2_x, f2_y])f1_distance = np.square(self.point_data - f1)f2_distance = np.square(self.point_data - f2)all_distance = np.sqrt(f1_distance[:, 0] + f1_distance[:, 1]) + np.sqrt(f2_distance[:, 0] + f2_distance[:, 1])Z = np.abs(2 * LAxis - all_distance)delta = math.sqrt(np.sum((Z - np.mean(Z)) ** 2) / len(Z))# Update inliers setinliers = np.nonzero(Z < 0.8 * delta)[0]inlier_pnts = self.point_data[inliers]return len(inlier_pnts), inlier_pntsdef execute_ransac(self):Time_start = time.time()while math.ceil(self.items):# 1.select N points at randomselect_points = self.random_sampling(self.N)# 2.fitting N ellipse pointsellipse = cv2.fitEllipse(select_points)# 3.assess model and calculate inliers pointsinliers_count, inliers_set = self.eval_model(ellipse)# 4.number of new inliers points more than number of old inliers points ?if inliers_count > self.count:ellipse_ = cv2.fitEllipse(inliers_set)  # fitting ellipse for inliers pointsself.count = inliers_count  # Update inliers setself.best_model = ellipse_  # Update best ellipse# 5.number of inliers points reach the expected valueif self.count > self.max_inliers:print('the number of inliers: ', self.count)break# Update itemsself.items = math.log(1 - self.P) / math.log(1 - pow(inliers_count/self.length, self.N))return self.best_modelif __name__ == '__main__':# 这个是根据我的工程实际问题写的寻找椭圆轮廓点，你们可以根据自己实际来该# 1.find ellipse edge line contours, hierarchy = cv2.findContours(img, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)# 2.Ransac fit ellipse parampoints_data = np.reshape(contours, (-1, 2))  # ellipse edge points setRansac = RANSAC(data=points_data, threshold=1., P=.99, S=.9, N=5)(X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()

整篇Ransac拟合椭圆的代码我感觉其实还是有不少小问题的，毕竟我菜鸡，写不出啥好东西，大家可以参考参考编程思路我感觉就可以了，我感觉我这个编程思路应该问题不大

或者你觉得你想顺便练习下矩阵的运算，也可以自己写椭圆方程的参数拟合过程，网上有很多关于椭圆的参数计算的，可以去参考参考：

import cv2
import math
import random
import numpy as np
from numpy.linalg import inv, svd, detclass RANSAC:def __init__(self, data, items, threshold, P, S):self.point_data = data  # 数据点self.items = items  # 迭代次数self.threshold = threshold  # 数据和模型之间可接受的差值self.P = P  # 希望的得到正确模型的概率self.size = len(self.point_data) * S  # 判断是否当前模型已经符合超过预设点self.label = Falseself.count = 0self.best_model = Nonedef random_sampling(self, n):"""随机取n个数据点"""all_point = self.point_datanp.random.shuffle(all_point)select_point = all_point[:n]return select_pointdef make_model(self, select_point):# create ellipse modelpoint1, point2, point3, point4, point5 = select_pointMatA = np.matrix([[point1[0] * point1[0], point1[0] * point1[1], point1[1] * point1[1], point1[0], point1[1], 1],[point2[0] * point2[0], point2[0] * point2[1], point2[1] * point2[1], point2[0], point2[1], 1],[point3[0] * point3[0], point3[0] * point3[1], point3[1] * point3[1], point3[0], point3[1], 1],[point4[0] * point4[0], point4[0] * point4[1], point4[1] * point4[1], point4[0], point4[1], 1],[point5[0] * point5[0], point5[0] * point5[1], point5[1] * point5[1], point5[0], point5[1], 1]])U, D, V = svd(MatA)VT = np.transpose(V)V05 = VT[0, 5]if math.fabs(V05) < 0.000001:return 0A = 1B = VT[1, 5] / V05C = VT[2, 5] / V05D = VT[3, 5] / V05E = VT[4, 5] / V05F = VT[5, 5] / V05thres1 = 4 * A * C - B * BellParamMat = np.matrix([[2 * A, B, D], [B, 2 * C, E], [D, E, 2 * F]])thres2 = det(ellParamMat) * (A + C)if (thres1 <= 0) or (thres2 >= 0):return 0# get ellipse paramx = (B * E - 2 * C * D) / (4 * A * C - B * B)y = (B * D - 2 * A * E) / (4 * A * C - B * B)if abs(B) <= 0.0001 and A < C:Angle = 0elif abs(B) <= 0.0001 and A > C:Angle = 90elif A < C:Angle = 0.5 * math.atan(B / (A - C)) * 180 / math.pielse:Angle = 90 + 0.5 * math.atan(B / (A - C)) * 180 / math.piepTemp1 = A * x * x + C * y * y + B * x * y - FepTemp2 = A + CepTemp3 = math.sqrt((A - C) * (A - C) + B * B)LAxis = math.sqrt(2 * epTemp1 / (epTemp2 - epTemp3))SAxis = math.sqrt(2 * epTemp1 / (epTemp2 + epTemp3))return x, y, LAxis, SAxis, Angle, [A, B, C, D, E, F]def eval_model(self, model):Count = 0x, y, LAxis, SAxis, Angle, EllipseFunc_Param = model  # model paramA, B, C, D, E, F = EllipseFunc_Param  # elliptic equation coefficient# calculate the ellipse focusAxis = math.sqrt(LAxis * LAxis - SAxis * SAxis)f1_x = x - Axis * math.cos(Angle * math.pi / 180)f1_y = y - Axis * math.sin(Angle * math.pi / 180)f2_x = x + Axis * math.cos(Angle * math.pi / 180)f2_y = y + Axis * math.sin(Angle * math.pi / 180)F1F2 = math.sqrt((f1_x - f2_x) ** 2 + (f1_y - f2_y) ** 2)if (2 * LAxis) <= F1F2:  # should: 2a > |F1F2|print('The model does not satisfy the definition of elliptic equation')return 0for point in self.point_data:dist1 = math.sqrt((point[0] - f1_x) ** 2 + (point[1] - f1_y) ** 2)dist2 = math.sqrt((point[0] - f2_x) ** 2 + (point[1] - f2_y) ** 2)dist = dist1 + dist2# |PF1| + |PF2| = 2aif math.fabs(dist - 2 * LAxis) < self.threshold:Count += 1return Countdef execute_ransac(self):while self.items > 0:# random select 5 pointsselect_point = self.random_sampling(n=5)# create model and get ellipse parammodel = self.make_model(select_point)if model == 0:continue# eval model and calculate number of inter pointsCount = self.eval_model(model)# number of new inter points more than number of old inter pointsif Count > self.count:self.count = Count  # Update internal pointsself.best_model = model  # Save the best model# Update the number of iterations dynamicallyself.items = int(math.log(1 - self.P) / math.log(1 - pow(Count / len(self.point_data), 5)))print(self.items)# inter points reach the expected valueif self.count > self.size:breakself.items -= 1return self.best_modeldef data_generator(contour):x_data = []y_data = []for point in contour:x_data.append(point[0])y_data.append(point[1])return x_data, y_dataif __name__ == '__main__':# 这个是根据我的工程实际问题写的寻找椭圆轮廓点，你们可以根据自己实际来该# 1.find ellipse edge line contours, hierarchy = cv2.findContours(img, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)points_data = np.reshape(contours, (-1, 2))  # ellipse edge points set# 2.Ransac fit ellipse param Ransac = RANSAC(data=points_data, items=999, threshold=1, P=.96, S=.56)(X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()

补充 - 2021.06.13（差点忘记补充了，写完一直忘记）

class RANSAC:def __init__(self, img, data, n, max_refines, max_inliers, threshold, show):self.img, self.orig_img = img  # 分割后图片、label图片self.point_data = data  # 椭圆轮廓点集self.length = len(self.point_data)  # 椭圆轮廓点集长度self.show = showself.n = n  # 随机选取n个点self.max_items = 999  # 迭代次数self.inliers_num = 0  # 内点计数器self.max_refines = max_refines  # 最大重迭代次数self.max_perc_inliers = max_inliers  # 最大内点比例self.error_threshold = threshold  # 模型误差容忍阀值self.best_model = ((0, 0), (1e-6, 1e-6), 0)  # 椭圆模型存储器def random_sampling(self, n):"""随机取n个数据点"""select_point = []all_point = self.point_datak, m = divmod(len(all_point), n)list_split = [all_point[i * k + min(i, m):(i + 1) * k + min(i + 1, m)] for i in list(range(n))]select_patch = random.sample(list_split, n // 2)for list_ in select_patch:temp = random.sample(list(list_), n // 2 // 2)select_point.append(temp[0])select_point.append(temp[1])return np.array(select_point)# return np.asarray(random.sample(list(self.point_data), n))def Geometric2Conic(self, ellipse):""" Calculate the elliptic equation coefficients """# Ax ^ 2 + Bxy + Cy ^ 2 + Dx + Ey + F(x0, y0), (bb, aa), phi_b_deg = ellipse# Semimajor and semiminor axesa, b = aa / 2, bb / 2# Convert phi_b from deg to radphi_b_rad = phi_b_deg * np.pi / 180.0# Major axis unit vectorax, ay = -np.sin(phi_b_rad), np.cos(phi_b_rad)# Useful intermediatesa2 = a * ab2 = b * b# Conic parametersif a2 > 0 and b2 > 0:A = ax * ax / a2 + ay * ay / b2B = 2 * ax * ay / a2 - 2 * ax * ay / b2C = ay * ay / a2 + ax * ax / b2D = (-2 * ax * ay * y0 - 2 * ax * ax * x0) / a2 + (2 * ax * ay * y0 - 2 * ay * ay * x0) / b2E = (-2 * ax * ay * x0 - 2 * ay * ay * y0) / a2 + (2 * ax * ay * x0 - 2 * ax * ax * y0) / b2F = (2 * ax * ay * x0 * y0 + ax * ax * x0 * x0 + ay * ay * y0 * y0) / a2 + \(-2 * ax * ay * x0 * y0 + ay * ay * x0 * x0 + ax * ax * y0 * y0) / b2 - 1else:# Tiny dummy circle - response to a2 or b2 == 0 overflow warningsA, B, C, D, E, F = (1, 0, 1, 0, 0, -1e-6)# Compose conic parameter arrayconic = np.array((A, B, C, D, E, F))return conicdef ConicFunctions(self, points, ellipse):"""Calculate various conic quadratic curve support functionsParameters----points : n x 2 array of floatsellipse : tuple of tuplesReturns----distance : array of floatsgrad : array of floatsabsgrad : array of floatsnormgrad : array of floats"""# Convert from geometric to conic ellipse parametersconic = self.Geometric2Conic(ellipse)# Row vector of conic parameters (Axx, Axy, Ayy, Ax, Ay, A1) (1 x 6)C = np.array(conic)# Extract vectors of x and y valuesx, y = points[:, 0], points[:, 1]# Construct polynomial array (6 x n)X = np.array((x * x, x * y, y * y, x, y, np.ones_like(x)))# Calculate Q/distance for all points (1 x n)distance = C.dot(X)# Construct conic gradient coefficients vector (2 x 3)Cg = np.array(((2 * C[0], C[1], C[3]), (C[1], 2 * C[2], C[4])))# Construct polynomial array (3 x n)Xg = np.array((x, y, np.ones_like(x)))# Gradient array (2 x n)grad = Cg.dot(Xg)# Normalize gradient -> unit gradient vectorabsgrad = np.sqrt(np.sqrt(grad[0, :] ** 2 + grad[1, :] ** 2))normgrad = grad / absgradreturn distance, grad, absgrad, normgraddef EllipseError(self, points, ellipse):# Calculate algebraic distances and gradients of all points from fitted ellipsedistance, grad, absgrad, normgrad = self.ConicFunctions(points, ellipse)# eps = 2.220446049250313e-16eps = np.finfo(distance.dtype).epsabsgrad = np.maximum(absgrad, eps)# Calculate error from distance and gradienterr = distance / absgradreturn errdef EllipseNormError(self, points, ellipse):(x0, y0), (bb, aa), phi_b_deg = ellipse# Semiminor axisb = bb / 2# Convert phi_b from deg to radphi_b_rad = phi_b_deg * np.pi / 180.0# Minor axis vectorbx, by = np.cos(phi_b_rad), np.sin(phi_b_rad)# Point one pixel out from ellipse on minor axisp1 = np.array((x0 + (b + 1) * bx, y0 + (b + 1) * by)).reshape(1, 2)# Error at this pointerr_p1 = self.EllipseError(p1, ellipse)# Errors at provided pointserr_pnts = self.EllipseError(points, ellipse)return err_pnts / err_p1def OverlayRANSACFit(self, overlay, outlier_points, inlier_points, best_ellipse):""" Ransac IMG """# outlier_points in greenfor col, row in outlier_points:overlay[row, col] = [0, 255, 0]# inlier_points in redfor col, row in inlier_points:overlay[row, col] = [0, 0, 255]""" Seg_Orign """# all points in greenfor col, row in self.point_data:self.img[row, col] = [0, 255, 0]# inlier fitted ellipse in redcv2.ellipse(self.img, best_ellipse, (0, 0, 255), 1)""" Label_Orign """# all points in greenfor col, row in self.point_data:self.orig_img[row, col] = [0, 255, 0]# # inlier fitted ellipse in redcv2.ellipse(self.orig_img, best_ellipse, (0, 0, 255), 1)def execute_ransac(self):count = 0while count < int(self.max_items) + 1:count += 1# 1.select n points at randomselect_points = self.random_sampling(self.n)# 2.fitting selected ellipse pointsellipse = cv2.fitEllipse(select_points)# 3.Refine inliers iterativelyfor refine in range(self.max_refines):# Calculate normalized errors for all pointsnorm_err = self.EllipseNormError(self.point_data, ellipse)# Identify inliers and outlierinliers = np.nonzero(norm_err ** 2 < self.error_threshold)[0]outlier = np.nonzero(norm_err ** 2 >= self.error_threshold)[0]# Update inliers and outlier setinlier_points = self.point_data[inliers]outlier_points = self.point_data[outlier]if inliers.size < 5:break# Fit ellipse to refined inliers setellipse = cv2.fitEllipse(np.asarray(inlier_points))# 4.Count inliers ratioinliers_lenght = inliers.sizeperc_inliers = inliers_lenght / self.lengthif inliers_lenght > self.inliers_num:# 5.Update paramself.best_model = ellipseself.inliers_num = inliers_lenght# 6.Judge whether the current inliers is greater than the maximum inliersif perc_inliers * 100.0 > self.max_perc_inliers:breakif perc_inliers > 0.08:self.max_items = count + math.log(1 - 0.99) / math.log(1 - pow(perc_inliers, self.n))# Update fitted image and displayif self.show:overlay = np.zeros((512, 512, 3), np.uint8)self.OverlayRANSACFit(overlay, outlier_points, inlier_points, self.best_model)cv2.imshow('Ransac', overlay)cv2.imshow('Seg_Orign', self.img)cv2.imshow('Label_Orign', self.orig_img)cv2.waitKey(0)cv2.destroyAllWindows()return self.best_model