Python基于随机采样一至性实现拟合椭圆(优化版)
天人合一peng 人气:0上次版本如果在没有找到轮廓或轮廓的点集数很小无法拟合椭圆或在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()加载全部内容