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")

结果:

输入图像

opencv 拟合圆 opencv椭圆拟合获取圆心_ci

分类后结果

opencv 拟合圆 opencv椭圆拟合获取圆心_拟合_02

拟合后结果

opencv 拟合圆 opencv椭圆拟合获取圆心_拟合_03

小结

本文主要是介绍一种可行的特征检测方案,供大家参考学习。简单记录自己的学习生涯。如有不当请多多指教。如果需要更详细的解释,更好的算法,可以联系作者本人。定将全力帮助。