ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Ransac拟合椭圆

2021-05-09 12:58:27  阅读:230  来源: 互联网

标签:Ransac 椭圆 self np inliers 拟合 ay ax


一、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 cv2
	import math
	import random
	import numpy as np
	from numpy.linalg import inv, svd, det

	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 = 999

	        self.count = 0  # 内点计数器
	        self.best_model = ((0, 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 Geometric2Conic(self, ellipse):
	    # 这个部分参考了GitHub中的一位大佬的,但是时间太久,忘记哪个人的了
	        """计算椭圆方程系数"""
	        # Ax ^ 2 + Bxy + Cy ^ 2 + Dx + Ey + F
	        (x0, y0), (bb, aa), phi_b_deg = ellipse
	
	        a, b = aa / 2, bb / 2  # Semimajor and semiminor axes
	        phi_b_rad = phi_b_deg * np.pi / 180.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
	        a2 = a * a
	        b2 = b * b
	
	        # Conic parameters
	        if a2 > 0 and b2 > 0:
	            A = ax * ax / a2 + ay * ay / b2
	            B = 2 * ax * ay / a2 - 2 * ax * ay / b2
	            C = ay * ay / a2 + ax * ax / b2
	            D = (-2 * ax * ay * y0 - 2 * ax * ax * x0) / a2 + (2 * ax * ay * y0 - 2 * ay * ay * x0) / b2
	            E = (-2 * ax * ay * x0 - 2 * ay * ay * y0) / a2 + (2 * ax * ay * x0 - 2 * ax * ax * y0) / b2
	            F = (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 - 1
	        else:
	            # Tiny dummy circle - response to a2 or b2 == 0 overflow warnings
	            A, B, C, D, E, F = (1, 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.Geometric2Conic(ellipse)
	        E = 4 * a * c - b * b
	        if E <= 0:
	            # print('this is not an ellipse')
	            return 0, 0
	
	        #  which long axis ?
	        (x, y), (LAxis, SAxis), Angle = ellipse
	        LAxis, SAxis = LAxis / 2, SAxis / 2
	        if SAxis > LAxis:
	            temp = SAxis
	            SAxis = LAxis
	            LAxis = temp
	
	        # calculate focus
	        Axis = 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 points
	        f1, 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 set
	        inliers = np.nonzero(Z < 0.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):
	            # 1.select N points at random
	            select_points = self.random_sampling(self.N)
	          
	            # 2.fitting N ellipse points
	            ellipse = cv2.fitEllipse(select_points)
	
	            # 3.assess model and calculate inliers points
	            inliers_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 points
	                self.count = inliers_count  # Update inliers set
	                self.best_model = ellipse_  # Update best ellipse
	
	                # 5.number of inliers points reach the expected value
	                if self.count > self.max_inliers:
	                    print('the number of inliers: ', self.count)
	                    break
	
	                # Update items
	                self.items = math.log(1 - self.P) / math.log(1 - pow(inliers_count/self.length, self.N))
	
	        return self.best_model
	
	if __name__ == '__main__':
	# 这个是根据我的工程实际问题写的寻找椭圆轮廓点,你们可以根据自己实际来该
	    # 1.find ellipse edge line 
	    contours, hierarchy = cv2.findContours(img, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)

	    # 2.Ransac fit ellipse param
	    points_data = np.reshape(contours, (-1, 2))  # ellipse edge points set
	    Ransac = RANSAC(data=points_data, threshold=1., P=.99, S=.9, N=5)
	    (X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()

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

四、部分对ransac的改进论文

链接:https://pan.baidu.com/s/1_P3rQJhRHMTE8sFmqRsxRg
提取码:0lia

标签:Ransac,椭圆,self,np,inliers,拟合,ay,ax
来源: https://blog.csdn.net/qq_41994220/article/details/116562339

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有