DBSCAN+椭圆拟合
首先,这篇文章使用的方法,已经有论文出版,如果需要使用,请注意版权的保护,论文名为Binocular Vision System Integrated with Transcranial Magnetic Stimulation Automatic Therapeutic Apparatus。
优点:
1.计算精度高
2.算法原理简单,
缺点:
1.计算复杂度较高,实时性不好
代替方案:可以先寻找边缘,然后findconter,之后再拟合,或者计算点集的矩,获得中心点。但是这个缺点就是精度稍微差一点。
原理概述
至于DBSCAN在中已经有介绍,大家可以到文章中提到的链接直接观察到DBSCAN算法的过程。第二个就是椭圆拟合,最小二乘椭圆拟合,网上有很多算法,在这里就不介绍了,本文采用的是Opencv中自带的椭圆拟合方法。cv.elipese()。
具体过程可以概述为 : 图像前处理----->得到图像边缘----》对边缘进行分类-------》得到分类结果------》对每个分类结果进行椭圆拟合------》评价那个椭圆中心是我们需要的
评价方法:因为特征是近圆的,所以第一个条件就是长轴和短轴近似相等。
其次由于特征是由两个同心圆组成的,故可以判断两个中心近似相等的即为我们所需要的。
python代码:
#__aditor:dwy
#__date: 2019/3/30
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sklearn as sk
from sklearn.cluster import DBSCAN
import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
import math
#stp1: read img
def getData(adress1):
#adress1 = "D:\\My code\\C_C++\\test_opencv\\data\\20190411\\left1\\16.tiff"
img = cv.imread(adress1)
img = cv.GaussianBlur(img,(5,5),0)
cv.namedWindow("img",cv.WINDOW_NORMAL)
cv.imshow("img",img)
cv.waitKey(10)
img_gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
thresh = cv.threshold(img_gray, 100, 255, type=cv.THRESH_BINARY)
thresh = thresh[1]
cv.namedWindow("thresh1", cv.WINDOW_NORMAL)
cv.imshow("thresh1", thresh)
cv.waitKey(10)
kernal1 = cv.getStructuringElement(cv.MORPH_RECT, (50, 50))
thresh = cv.erode(thresh, kernal1)
thresh = cv.dilate(thresh, kernal1)
cv.namedWindow("thresh", cv.WINDOW_NORMAL)
cv.imshow("thresh",thresh)
cv.waitKey(10)
kernal2 = cv.getStructuringElement(cv.MORPH_RECT, (20, 20))
img_gray2 = cv.dilate(thresh, kernal2)
img_gray2 = cv.erode(img_gray2, kernal2)
cv.namedWindow("open", cv.WINDOW_NORMAL)
cv.imshow("open",img_gray2)
kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])
img_gray2 = cv.filter2D(img_gray2, -1, kernel)
kernal2 = cv.getStructuringElement(cv.MORPH_RECT, (10, 10))
img_gray2 = cv.dilate(img_gray2, kernal2)
cv.namedWindow("laplace", cv.WINDOW_NORMAL)
cv.imshow("laplace",~img_gray2)
cv.waitKey(10)
return img_gray2
#stp2:data handle
def dataHandle(data):
weight,height = data.shape
output=[]
for i in range(height):
for j in range(weight):
if data[j][i]>0:
point = [i,j]
output.append(point)
output = np.array(output)
return output
#椭圆拟合
def circleMatch(points):
ellips = cv.fitEllipse(points)
return ellips
#stp3:using DBSCAN to classification
def testDBSCAN(input,address1):
img = cv.imread(address1)
#print(input.shape)
db = DBSCAN(eps=8,min_samples=10).fit(input)
y_pred = DBSCAN(eps=8,min_samples=10).fit_predict(input)
#绘制分类后的结果
plt.xlabel("u direction pixel coordinates")
plt.ylabel("v direction pixel coordinates")
plt.title("Result of DBSCAN",fontsize=10.5)
plt.scatter(input[:, 0], input[:, 1], c=y_pred, marker='.')
plt.show()
labels = db.labels_# 和X同一个维度,labels对应索引序号的值 为她所在簇的序号。若簇编号为-1,表示为噪声
# print('每个样本的簇标号:')
# print(labels)
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0) # 获取分簇的数目
elipses=[]
for i in range(n_clusters_):
one_cluster = input[labels == i]
plt.xlabel("u direction pixel coordinates")
plt.ylabel("v direction pixel coordinates")
string = "Result of "+str(i)+"th category"
plt.title(string,fontsize=10.5)
if i>=0:
ellips=circleMatch(one_cluster)
elipses.append(ellips)
#cv.ellipse(img,ellips,(0,0,255),10)
plt.plot(one_cluster[:, 0], one_cluster[:, 1],'.')
plt.show()
output_elipses=[]
points =[]
ranges = len(elipses)
j=0
#提取目标椭圆部分
#暴力查找,如果两个中心近似相等我们认为这是我们需要的,此处可以用dict优化,先排序,后搜索
while(j<ranges):
i=0
while(i<ranges):
if(i==j):
i+=1
else:
if abs((elipses[j][0][0]-elipses[i][0][0]))<5 and abs((elipses[j][0][1]-elipses[i][0][1]))<5:
output_elipses.append(elipses[j])
break
else:
i+=1
j+=1
range2=len(output_elipses)
points=[]
m=0
#将目标椭圆绘制到图上,并且计算出中心位置像素坐标
while(m<range2):
print("第%d个椭圆的参数为:"%(m+1),output_elipses[m])
cv.ellipse(img, output_elipses[m], (0, 0, 255), 10)
n=m+1
while(n<range2):
if abs((output_elipses[m][0][0]-output_elipses[n][0][0]))<5 and abs((output_elipses[m][0][1]-output_elipses[n][0][1]))<5:
pointx=(output_elipses[m][0][0]+output_elipses[n][0][0])/2
pointy=(output_elipses[m][0][1]+output_elipses[n][0][1])/2
point=[pointx,pointy]
points.append(point)
break
else:
n+=1
m+=1
po=0
title=['A','B','C']
for p in points:
print("第%d个提取中心点位:" % po , p)
t=tuple(map(int,p))
cv.circle(img ,t,10,(0,0,255),-1,4)
cv.putText(img,title[po],t,cv.FONT_HERSHEY_COMPLEX,5,(0,0,255),5)
po=po+1
cv.namedWindow("imshow",cv.WINDOW_NORMAL)
cv.imshow("imshow",img)
cv.waitKey(10)
return points
#stp4:距离计算,因为此处用到了立体标定,计算两点之间的距离,不用可以忽略。
def caculateDistance(pointL,pointR,Q):
PLx = np.array([pointL[0][0],pointL[1][0],pointL[2][0]])
PRx = np.array([pointR[0][0],pointR[1][0],pointR[2][0]])
PLy = np.array([pointL[0][1],pointL[1][1],pointL[2][1]])
output = []
for i in range(3):
d = PLx[i] - PRx[i]
input = np.array([PLx[i], PLy[i], d, 1])
out = np.dot(Q, input)
out2 = [out[0] / out[3], out[1] / out[3], out[2] / out[3]]
output.append(out2)
print("相机坐标系下第%d个坐标值为:(%f,%f,%f)" % (i + 1, out[0] / out[3], out[1] / out[3], out[2] / out[3]))
d1 = math.sqrt(np.square(output[0][0] - output[1][0]) + np.square(output[0][1] - output[1][1]) + np.square(
output[0][2] - output[1][2]))
d2 = math.sqrt(np.square(output[0][0] - output[2][0]) + np.square(output[0][1] - output[2][1]) + np.square(
output[0][2] - output[2][2]))
d3 = math.sqrt(np.square(output[1][0] - output[2][0]) + np.square(output[1][1] - output[2][1]) + np.square(
output[1][2] - output[2][2]))
print("点1和点2两点之间的距离单位mm:%f" % (d1 * 1000))
print("点1和点3两点之间的距离单位mm:%f" % (d2 * 1000))
print("点2和点3两点之间的距离单位mm:%f" % (d3 * 1000))
if __name__=="__main__":
address1 = "G:\\cameraData\\0422-first\\left\\rectifyImageL.jpg"
address2 = "G:\\cameraData\\0422-first\\right\\rectifyImageR.jpg"
data = getData(address1)
db = dataHandle(data)
pointL=testDBSCAN(db,address1)
#print(pointL)
print("---------------------------------")
data = getData(address2)
db = dataHandle(data)
pointR=testDBSCAN(db,address2)
#print(pointR)
print("-------------计算结果--------------------")
Q=np.array([[1,0,0,2138.81555],[0,1,0,-1868.5245],[0,0,0,20079.8178],[0,0,4.55778,39719.39836]])
caculateDistance(pointL,pointR,Q)
cv.waitKey(0)
cv.destroyAllWindows()
plt.close("all")
结果:
输入图像
分类后结果
拟合后结果
小结
本文主要是介绍一种可行的特征检测方案,供大家参考学习。简单记录自己的学习生涯。如有不当请多多指教。如果需要更详细的解释,更好的算法,可以联系作者本人。定将全力帮助。