Python基于随机采样一至性实现拟合椭圆(优化版)
作者:天人合一peng
这篇文章主要对上一版的Python基于随机采样一至性实现拟合椭圆的优化,文中的示例代码讲解详细,具有一定的借鉴价值,感兴趣的可以了解一下
上次版本如果在没有找到轮廓或轮廓的点集数很小无法拟合椭圆或在RANSAC中寻找最优解时会死循环中,优化后的代码
import cv2 import os import numpy as np import matplotlib.pyplot as plt import math from Ransac_Process import RANSAC def cul_area(x_mask, y_mask, r_circle, mask): mask_label = mask.copy() num_area = 0 for xm in range(x_mask+r_circle-10, x_mask+r_circle+10): for ym in range(y_mask+r_circle-10, y_mask+r_circle+10): # print(mask[ym, xm]) if (pow((xm-x_mask), 2) + pow((ym-y_mask), 2) - pow(r_circle, 2)) == 0 and mask[ym, xm][0] == 255: num_area += 1 mask_label[ym, xm] = (0, 0, 255) cv2.imwrite('./test2/mask_label.png', mask_label) print(num_area) return num_area def mainFigure(img, point0): point_center = [] # cv2.imwrite('./test2/img_source.png', img) img_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV) # cv2.imwrite('./test2/img_hsv.png', img_hsv) w, h = img.shape[1], img.shape[0] w_hsv, h_hsv = img_hsv.shape[1], 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][0] < 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][0] < 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, 255 else: img_hsv[j_hsv, i_hsv] = 0, 0, 0 # cv2.imwrite('./test2/img_hsvhb.png', img_hsv) # cv2.imshow("hsv", img_hsv) # cv2.waitKey() # 灰度化处理图像 grayImage = cv2.cvtColor(img_hsv, cv2.COLOR_BGR2GRAY) # mask = np.zeros((grayImage.shape[0], grayImage.shape[1]), np.uint8) # mask = cv2.cvtColor(mask, cv2.COLOR_GRAY2BGR) # cv2.imwrite('./mask.png', mask) # 尝试寻找轮廓 contours, hierarchy = cv2.findContours(grayImage, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE) # 合并轮廓 if len(contours) > 1: contours_merge = np.vstack([contours[0], contours[1]]) for i in range(2, len(contours)): contours_merge = np.vstack([contours_merge, contours[i]]) # cv2.drawContours(img, contours_merge, -1, (0, 255, 255), 1) # # cv2.imwrite('./test2/img_res.png', img) # cv2.imshow("contours_merge", img) # cv2.waitKey() elif len(contours) == 1: contours_merge = contours[0] else: print("No contours!") return 0,0 # RANSAC拟合 points_data = np.reshape(contours_merge, (-1, 2)) # ellipse edge points set # print("points_data", len(points_data)) # 2.Ransac fit ellipse param # Ransac = RANSAC(data=points_data, threshold=0.1, P=.99, S=.5, N=10) Ransac = RANSAC(data=points_data, threshold=0.5, P=.98, S=.6, N=10) ellipse_values = Ransac.execute_ransac() # 检测到轮廓里数量太少(<5)则无法拟合椭圆 if ellipse_values is None: return 0,0 (X, Y), (LAxis, SAxis), Angle = ellipse_values # print( (X, Y), (LAxis, SAxis)) # 拟合圆 cv2.ellipse(img, ((X, Y), (LAxis, SAxis), Angle), (0, 0, 255), 1, cv2.LINE_AA) # 画圆 cv2.circle(img, (int(X), int(Y)), 3, (0, 0, 255), -1) # 画圆心 point_center.append(int(X)) point_center.append(int(Y)) # 直接拟合 # rrt = cv2.fitEllipse(contours_merge) # x, y)代表椭圆中心点的位置(a, b)代表长短轴长度,应注意a、b为长短轴的直径,而非半径,angle 代表了中心旋转的角度 # # print("rrt", rrt) # cv2.ellipse(img, rrt, (255, 0, 0), 1, cv2.LINE_AA) # 画圆 # x, y = rrt[0] # cv2.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)) # # # 两种方法坐标的距离 # dis_two_method = math.sqrt(math.pow(X - x, 2) + math.pow(Y - y, 2)) # print("两种方法坐标的距离", dis_two_method) cv2.imshow("fit circle", img) cv2.waitKey(3) # cv2.imwrite("./test2/fitcircle.png", img) return point_center[0], point_center[1] if __name__ == "__main__": # 测试所有图片 mainFolder = "./Images/save_img" myFolders = os.listdir(mainFolder) print("myFolders", myFolders) myImageList = [] path = '' for folder in myFolders: path = mainFolder + '/' + folder myImageList = os.listdir(path) # print(myImageList) # print(f'Tatal images deteted is {len(myImageList)}') i = 0 for imagN in myImageList: curImg = cv2.imread(f'{path}/{imagN}') # images.append(curImg) print(f'{path}/{imagN}') point0 = [0, 0] cir_x, cir_y = mainFigure(curImg, point0) print("This is ", i, "圆心为",(cir_x, cir_y)) i += 1 # # 测试2 # for i in range(1,6): # imageName = "s" # imageName += str(i) # path = './Images/danHoles/' + imageName + '.png' # print(path) # img = cv2.imread(path) # point0 = [0, 0] # cir_x, cir_y = mainFigure(img, point0) # # 测试1 # img = cv2.imread('./Images/danHoles/s6.png') # point0 = [0, 0] # cir_x, cir_y = mainFigure(img, point0)
Ransac_Process.py
import cv2 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 = 8 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 if len(all_point) >= n: select_point = np.asarray(random.sample(list(all_point), n)) return select_point else: print("轮廓点数太少,数量为", len(all_point)) return None 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): # print(self.max_inliers) # 1.select N points at random select_points = self.random_sampling(self.N) # 当从轮廓中采集的点不够拟合椭圆时跳出循环 if select_points is None or len(select_points) < 5: # print(select_points) return None else: # 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: if len(inliers_set) > 4: ellipse_ = cv2.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) # 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 # print(math.log(1 - pow(inliers_count / self.length, self.N))) if math.log(1 - pow(inliers_count / self.length, self.N)) != 0: self.items = math.log(1 - self.P) / math.log(1 - pow(inliers_count / self.length, self.N)) Time_end = time.time() # print(Time_end - Time_start ) if Time_end - Time_start >= 4: # print("time is too long") break return self.best_model if __name__ == '__main__': # 1.find ellipse edge line contours, hierarchy = cv2.findContours(grayImage, 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=0.5, P=.99, S=.618, N=10) (X, Y), (LAxis, SAxis), Angle = Ransac.execute_ransac()
以上就是Python基于随机采样一至性实现拟合椭圆(优化版)的详细内容,更多关于Python拟合椭圆的资料请关注脚本之家其它相关文章!