基于opencv的一种快速有效椭圆检测方法

 

 

      本篇介绍的椭圆检测方法来自以下论文,论文作者提供了测试代码。本文主要是对这个方法做出详解。

       参考论文:A fast and effective ellipse detector for embedded vision applications

       代码链接:fast_ellipse_detector

 一、调用及参数说明


Size	szPreProcessingGaussKernelSize = Size(5,5); //高斯滤波窗体大小
double	dPreProcessingGaussSigma = 1.0;//高斯滤波 Sigma
float 	fThPosition = 1.0;//依据相关凸性(弧的相对位置)移除某些点
float	fMaxCenterDistance = 100*0.05;// 两个中心点最大距离
int		iMinEdgeLength = 16;//单段弧最小边界尺寸
float	fMinOrientedRectSide = 3.0;//最小包含圆弧方向包围盒尺寸
float	fDistanceToEllipseContour = 1.0;//椭圆上的点到椭圆边界的最大容忍距离
float	fMinScore = 0.4;//椭圆有效的最低分数
float	fMinReliability= 0.4;//椭圆有效的最低辅助分数
int     iNs = 16;//此参数决定 遍历还是抽样弧上的点


int main()
{
    Mat1b input;//输入图像
    vector<Ellipse> ellipseDst;//输出椭圆

    CEllipseDetectorYaed yaed;//初始化椭圆检测算子

    //初始化参数
    yaed.SetParameters( szPreProcessingGaussKernelSize,
	                    dPreProcessingGaussSigma,
	                    fThPosition,
	                    fMaxCenterDistance,
	                    iMinEdgeLength,
	                    fMinOrientedRectSide,
	                    fDistanceToEllipseContour,
	                    fMinScore,
	                    fMinReliability,
	                    iNs
	                    );
    //调用椭圆检测
    yaed.Detect(input, ellipseDst);

    return 0;
}


 

 二、算法流程

                                                                       

                                                                     图1  算法流程图

  1. 作者是利用高斯滤波做预处理。我将这一部分摘出来,根据实际需要做预处理,预处理后的图像作为输入。
  2. 边界检测部分用到了自适应Canny检测。
  3. 将边界分为凹弧和凸弧,根据输入参数筛选弧段。
  4. 利用弧段来估计椭圆参数,交叉验算得出椭圆中心点,计算出椭圆得分。
  5. 利用圆心、长短轴和旋转角度聚类。

 

三、算法分析

下面结合论文和代码对每个步骤分析。

本文在作者提供的源码基础上做了些许调整,使用起来更加方便。调整如下:

  1. 将与处理过程剥离开,增加了参数bool _isSmooth。_isSmooth为ture时,PrePeocessing使用高斯滤波,_isSmooth为false时,不使用高斯滤波;根据实际应用场景做不同预处理,预处理后的图像直接检测椭圆。
  2. 边界检测自适应canny变换时,参数percent_of_pixels_not_edges 和 threshold_ratio是默认值。将整两个参数放出来可以调整,_dNotEdges, _dThMinRatioMax。
  3. 变量、类名称的更改。CEllipseDetectorYaed 改为CEllipseDetector;Ellipse改为Ellipse_Common。

 

Detect(Mat1b& I, vector<Ellipse_Common>& ellipses)
{
	// Set the image size
	//获取图像size
	_szImg = I.size();

	// Initialize temporary data structures
	//初始化 凸弧和凹弧
	Mat1b DP = Mat1b::zeros(_szImg);		// arcs along positive diagonal
	Mat1b DN = Mat1b::zeros(_szImg);		// arcs along negative diagonal

	// Initialize accumulator dimensions
	//初始化累加器的size
	ACC_N_SIZE = 101;
	ACC_R_SIZE = 180;
	ACC_A_SIZE = max(_szImg.height, _szImg.width);

	// Allocate accumulators
	// 配置累加器
	accN = new int[ACC_N_SIZE];
	accR = new int[ACC_R_SIZE];
	accA = new int[ACC_A_SIZE];

	// Other temporary
	VVP points_1, points_2, points_3, points_4;		//属于1,2,3,4象限的点
	//vector of points, one for each convexity class
	map<uint, EllipseData> centers;		//已经计算好的椭圆数据重构的哈希图
	//hash map for reusing already computed EllipseData

	// Preprocessing
	// From input image I, find edge point with coarse convexity along positive (DP) or negative (DN) diagonal
	//预处理部分 计算出图像中属于凸弧、凹弧部分
	PrePeocessing(I, DP, DN, _dNotEdges, _dThMinRatioMax);

	if(is_Ellipsedebug)
	{
		imwrite("DP.bmp", DP);
		imwrite("DN.bmp", DN);
	}

	// Detect edges and find convexities
	// 凸弧中属于1,3象限 凹弧中属于2,4象限的点分类
	DetectEdges13(DP, points_1, points_3);
	DetectEdges24(DN, points_2, points_4);

	if(is_Ellipsedebug)
	{
		//把点plot出来保存数据
		Mat3b out(I.rows, I.cols, Vec3b(0,0,0));
		for(unsigned i=0; i<points_1.size(); ++i)
		{
			Vec3b color(255,0,0);
			for(unsigned j=0; j<points_1[i].size(); ++j)
				out(points_1[i][j]) = color;
		}

		for(unsigned i=0; i<points_2.size(); ++i)
		{
			Vec3b color(0,255,0);
			for(unsigned j=0; j<points_2[i].size(); ++j)
				out(points_2[i][j]) = color;
		}

		for(unsigned i=0; i<points_3.size(); ++i)
		{
			Vec3b color(0,0,255);
			for(unsigned j=0; j<points_3[i].size(); ++j)
				out(points_3[i][j]) = color;
		}

		for(unsigned i=0; i<points_4.size(); ++i)
		{
			Vec3b color(255,0,255);
			for(unsigned j=0; j<points_4[i].size(); ++j)
				out(points_4[i][j]) = color;
		}

		imwrite("out.bmp", out);//chenglm  2017.3.6
	}

	//find triplets
	//通过3个象限的弧估算出椭圆
	Triplets124(points_1, points_2, points_4, centers, ellipses);
	Triplets231(points_2, points_3, points_1, centers, ellipses);
	Triplets342(points_3, points_4, points_2, centers, ellipses);
	Triplets413(points_4, points_1, points_3, centers, ellipses);


	// Sort detected ellipses with respect to score
	//椭圆排序
	sort(ellipses.begin(), ellipses.end());

	// Free accumulator memory
	//释放累加器资源
	delete[] accN;
	delete[] accR;
	delete[] accA;


	// Cluster detections
	//椭圆聚类
	ClusterEllipses(ellipses);
}

3.1 预处理

预处理分为2个部分:滤波和边界检测。

//预处理函数 得到凸弧和凹弧
void PrePeocessing(Mat1b& I, Mat1b& DP, Mat1b& DN,double dNotEdges,
		           double dThRatio, bool isSmooth=false);

3.1.1滤波

滤波部分没有太多的可以分析,原文就是高斯滤波。也可以用其他方法。

  3.1.2边界检测

边界检测使用的是adaptCanny函数,与cv::Canny大体相同,有差异的地方是自动计算low_thresh和high_thresh。

//自适应Canny变换
void adaptCanny(cv::Mat image, cv::Mat& _edges,cv::Mat& _sobel_x, cv::Mat& _sobel_y,
	            int apertureSize, bool L2gradient,
	            double dNot_edges,
	            double dThreshold_ratio
	            )

//Canny变换
void cv::Canny( InputArray _src, OutputArray _dst,
                double low_thresh, double high_thresh,
                int aperture_size, bool L2gradient )

首先是计算双向的soble变换。

cvSobel( src, dx, 1, 0, aperture_size );
cvSobel( src, dy, 0, 1, aperture_size );

然后是计算梯度的峰值。

//% Calculate Magnitude of Gradient
//magGrad = hypot(dx, dy);

//计算梯度峰值

Mat1f magGrad(size.height, size.width, 0.f);
float maxGrad(0);
float val(0);
for(i=0; i<size.height; ++i)
{
	float* _pmag = magGrad.ptr<float>(i);
	const short* _dx = (short*)(dx->data.ptr + dx->step*i);
	const short* _dy = (short*)(dy->data.ptr + dy->step*i);
	for(j=0; j<size.width; ++j)
	{
		val = float(abs(_dx[j]) + abs(_dy[j]));
		_pmag[j] = val;
		maxGrad = (val > maxGrad) ? val : maxGrad;
	}
}

接着计算梯度直方图,将梯度压缩到NUM_BINS阶,作者源码中NUM_BINS=64。

        //% Normalize for threshold selection
	//normalize(magGrad, magGrad, 0.0, 1.0, NORM_MINMAX);

	//% Determine Hysteresis Thresholds

	//set magic numbers
	//在现有的灰度上分为64阶
	const int NUM_BINS = 64;
	const double percent_of_pixels_not_edges = dNot_edges;
	const double threshold_ratio = dThreshold_ratio;

	//compute histogram
	//计算梯度直方图
	int bin_size = cvFloor(maxGrad / float(NUM_BINS) + 0.5f) + 1;
	if (bin_size < 1) bin_size = 1;
	int bins[NUM_BINS] = { 0 };
	for (i=0; i<size.height; ++i)
	{
		float *_pmag = magGrad.ptr<float>(i);
		for(j=0; j<size.width; ++j)
		{
			int hgf = int(_pmag[j]);
			bins[int(_pmag[j]) / bin_size]++;
		}
	}

根据直方图计算出low_thresh和high_thresh。

        //% Select the thresholds
	//选择阈值:直方图内前N个灰度阶梯点数累计和到达阈值
	float total(0.f);
	float target = float(size.height * size.width * percent_of_pixels_not_edges);
	int low_thresh, high_thresh(0);

	while(total < target)
	{
		total+= bins[high_thresh];
		high_thresh++;
	}
	high_thresh *= bin_size;
	low_thresh = cvFloor(threshold_ratio * float(high_thresh));

接下来的部分和canny变换类似,就不再分析了。如果有想更进一步了解的,可以参考matlab的edge.m文件。其中有以下部分代码,与作者提供的源码是一样的。参考matlab参数,设置_dNotEdges=0.7, _dThMinRatioMax=0.4;

if strcmp(method,'canny')
    % Magic numbers
    PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
    ThresholdRatio = .4;          % Low thresh is this fraction of the high.
    
    % Calculate gradients using a derivative of Gaussian filter
    [dx, dy] = smoothGradient(a, sigma);
    
    % Calculate Magnitude of Gradient
    magGrad = hypot(dx, dy);
    
    % Normalize for threshold selection
    magmax = max(magGrad(:));
    if magmax > 0
        magGrad = magGrad / magmax;
    end
    
    % Determine Hysteresis Thresholds
    [lowThresh, highThresh] = selectThresholds(thresh, magGrad, PercentOfPixelsNotEdges, ThresholdRatio, mfilename);
    
    % Perform Non-Maximum Suppression Thining and Hysteresis Thresholding of Edge
    % Strength
    e = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh);
    thresh = [lowThresh highThresh];

小结:

adaptCanny变换中根据图像梯度直方图,自适应计算了low_thresh和high_thresh。当使用cv:Canny函数时,low_thresh和high_thresh取值不恰当,往往得出糟糕的结果。adaptCanny显得更加的灵活。

canny变换中low_thresh和high_thresh的意义为:

(1)当前点梯度值小于low_thresh,被标记为1,非边界点;

(2)当前点梯度值大于high_thresh,标记为2,边界点;

(3)当前点梯度值大于low_thresh且小于high_thresh,标记为0,可能的边界点。

当然被认为是边界还需要非极大值抑制抑制。

图像边界理解为像素值跳变处,sobel算子是求取图像梯度取较大值处为边界。adaptCanny函数中统计各梯度对应像素点数,由低到高各梯度值对应像素点数累加和大于非边界占比像素点数时。此时的梯度值为cv::Canny中高门限阈值high_thresh。低门限阈值low_thresh=high_thresh*_dThMinRatioMax,应当是由大量统计得出的经验值。

 

3.2够成椭圆的弧判定方法

3.2.1判定弧的凹凸性

在判定边界属于凹弧还是凸弧前,先考虑几个问题。边界定义为凹或者凸的依据是什么?判定凹凸性有什么作用呢?

首先参考函数的凹凸性。设函数f(x)在区间I上定义,若对I中任意2点x_{1}x_{2}都有:   

                                         f(\lambda x_{1}+(1-\lambda)x_{2}) <=\lambda f(x_{1}) + (1-\lambda) f(x_{2}) ,\lambda\in (0,1)                                                    

那么f(x)为凸函数;

设函数f(x)在区间I上定义,若对I中任意2点x_{1}x_{2}都有:   

                                         f(\lambda x_{1}+(1-\lambda)x_{2}) >=\lambda f(x_{1}) + (1-\lambda) f(x_{2}) ,\lambda\in (0,1)                                                    

那么f(x)为凹函数。

函数的凹凸性是描述函数图像弯曲方向的一个重要性质。函数的凹凸性还可以描述为:当且仅当一元可微函数在某个区间上的导数单调不减,此函数是凸的;反之则是凹函数。

那么边界上的点如何判定凹凸性呢?图像分别对x方向和y方向求一阶偏导得到\frac{\partial _{I}}{\partial _{x}}\frac{\partial _{I}}{\partial _{y}}。当\frac{\partial _{I}}{\partial _{x}} > 0时,I沿着x方向递增,当\frac{\partial _{I}}{\partial _{y}} < 0时,I沿着y方向递增(图像原点为左上点,因此y轴反向)。I沿着x方向和y方向同时递增(递减)时,当前点判定为凸(convex );反之判定为凹(concave )。

当前点的凹凸性有以下关系式确定。                                                      

                                                        \left\{\begin{matrix} \frac{\partial _{I}}{\partial _{x}} > 0 ,\frac{\partial _{I}}{\partial _{y}} > 0 , concave \\ \frac{\partial _{I}}{\partial _{x}} > 0 ,\frac{\partial _{I}}{\partial _{y}} < 0 ,convex \\ \frac{\partial _{I}}{\partial _{x}} < 0 ,\frac{\partial _{I}}{\partial _{y}} > 0 ,convex \\ \frac{\partial _{I}}{\partial _{x}} < 0 ,\frac{\partial _{I}}{\partial _{y}} < 0 ,concave \end{matrix}\right.

进一步的可以表达为:

                                                         \left\{\begin{matrix} \frac{\partial _{I}}{\partial _{x}} *\frac{\partial _{I}}{\partial _{y}} > 0 , concave \\ \frac{\partial _{I}}{\partial _{x}} *\frac{\partial _{I}}{\partial _{y}} < 0 ,convex \end{matrix}\right.

当判定了边界的凹凸性后,如何利用凹凸性呢?观察一个椭圆,以椭圆中心为原点构造一个直角坐标系。会发现第一象限和第三象限的弧是凸的,第二象限和第四象限的弧是凹的。y轴反向,所以与笛卡尔坐标系下函数图像看起来是相反的。

那么凹弧位于第一象限和第三象限,凸弧位于第二象限和第四象限可以构成椭圆吗?显然是无法构成椭圆的。

 

3.2.2判定弧属于哪个象限

接着如何判定能够构成椭圆的弧属于哪个象限呢?以弧段的起点和终点为对角线,构建一个包含矩形。包含矩形被弧划分为2个部分。凸弧上方面积小于下方面积时,此弧处于第一象限;凸弧上方面积大于下方面积时,此弧处于第三象限;凹弧上方面积小于下方面积时,此弧处于第二象限;凹弧上方面积大于下方面积时,此弧处于第四象限。总结如下

                                             \left\{\begin{matrix}if convex \cap U < O , frist \: quadrant & &\\ if concave \cap U < O ,second \: quadrant & & \\ if convex \cap U > O , third \: quadrant & & \\ if concave \cap U > O , fourth\: quadrant \end{matrix}\right.

 

此时边界分类完毕,全部的边界被分为4类:

(1)属于第一象限的凸弧;

(2)属于第二象限的凹弧;

(3)属于第三象限的凸弧;

(4)属于第四象限的凹弧。

下面结合代码进一步理解。作者在PrePeocessing函数中边界检测步骤后对弧段分类:凹弧和凸弧。

          // For each edge points, compute the edge direction
	//针对每个边界点计算边界方向
	for (int i = 0; i<_szImg.height; ++i)
	{
		short* _dx = DX.ptr<short>(i);
		short* _dy = DY.ptr<short>(i);
		uchar* _e = E.ptr<uchar>(i);
		uchar* _dp = DP.ptr<uchar>(i);
		uchar* _dn = DN.ptr<uchar>(i);

		for (int j = 0; j<_szImg.width; ++j)
		{
			if (!((_e[j] <= 0) || (_dx[j] == 0) || (_dy[j] == 0)))
			{
				// Angle of the tangent
				float phi = -(float(_dx[j]) / float(_dy[j]));//边界点的正切值  注意图像坐标的y轴于笛卡尔直角坐标的y轴相反

				// Along positive or negative diagonal
				if (phi > 0)	//属于凸弧
					_dp[j] = (uchar)255;
				else if (phi < 0)	//属于凹弧
					_dn[j] = (uchar)255;
			}
		}
	}

调用以下2各函数将凸弧分到第一象限和第三象限,将凹弧分配到第二象限和第四象限。

        //凸弧中区分第一象限和第三象限的部分
	void DetectEdges13(Mat1b& DP, VVP& points_1, VVP& points_3);
	//凹弧中区分第二象限和第四象限部分
	void DetectEdges24(Mat1b& DN, VVP& points_2, VVP& points_4);

以将凸弧为例

DetectEdges13(Mat1b& DP,      //输入凸弧
               VVP& points_1, //输出第一象限点集
               VVP& points_3) //输出第三象限点集
{
	// Vector of connected edge points
	VVP contours;

	// Labeling 8-connected edge points, discarding edge too small
	//抛弃过短的弧单个边界弧长度小于_iMinEdgeLength被抛弃
	Labeling(DP, contours, _iMinEdgeLength);
	int iContoursSize = int(contours.size());

	// For each edge
	for (int i = 0; i < iContoursSize; ++i)
	{
		VP& edgeSegment = contours[i];

#ifndef DISCARD_CONSTRAINT_OBOX   //利用最小包含矩形的尺寸约束

		// Selection strategy - Step 1 - See Sect [3.1.2] of the paper
		// Constraint on axes aspect ratio
		RotatedRect oriented = minAreaRect(edgeSegment);
		float o_min = min(oriented.size.width, oriented.size.height);

		if (o_min < _fMinOrientedRectSide)
		{
			continue;
		}
#endif

		// Order edge points of the same arc
		sort(edgeSegment.begin(), edgeSegment.end(), SortTopLeft2BottomRight);//将边界点按坐标从左上到右下排序(凹弧是左下到右上)
		int iEdgeSegmentSize = unsigned(edgeSegment.size());

		// Get extrema of the arc
		Point& left = edgeSegment[0];
		Point& right = edgeSegment[iEdgeSegmentSize - 1];

		// Find convexity - See Sect [3.1.3] of the paper
		int iCountTop = 0;
		int xx = left.x;
		for (int k = 1; k < iEdgeSegmentSize; ++k)//可以理解为梯形图积分求面积
		{
			if (edgeSegment[k].x == xx) continue;

			iCountTop += (edgeSegment[k].y - left.y);
			xx = edgeSegment[k].x;
		}

		int width = abs(right.x - left.x) + 1;
		int height = abs(right.y - left.y) + 1;
		int iCountBottom = (width * height) - iEdgeSegmentSize - iCountTop;

		if (iCountBottom > iCountTop)
		{	//弧下方面积大于上方面积,属于第一象限
			points_1.push_back(edgeSegment);
		}
		else if (iCountBottom < iCountTop)
		{	//弧下方面积小于上方面积,属于第三象限
			points_3.push_back(edgeSegment);
		}
	}
}

 

小结:

在获取边界后如何利用这些边界?通过一系列的筛选,去除干扰,将边界按照凹凸性和所在象限分为4个部分。只有这4个部分中的弧相互组合,才有可能构成椭圆。接下来就是利用这些弧拟和椭圆。

 

3.3 拟和椭圆参数

3.3.1估算椭圆中心

根据椭圆的性质:一组平行弦的中点连线必过椭圆中心。

 

证明如下:

设椭圆的方程为:\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1,(a>b>1),当弦垂直于x轴或y轴时,其中点连线必然过椭圆中心。当弦与x轴和y轴均不垂直时,设其方程为:y=kx+m

联立方程\left\{\begin{matrix} \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1,(a>b>1) \\ y=kx+m \end{matrix}\right.可得:

b^{2}x^{2}+a^{2}(kx+m)^{2}=a^{2}b^{2}

(b^{2}+a^{2}k^{2})x^{2}+2kma^{2}x+(a^{2}m^{2}-a^{2}b^{2})=0

b^{2}(y-m)^{2}+a^{2}k^{2}y^{2}=a^{2}b^{2}k^{2}

(b^{2}+a^{2}k^{2})y^{2}-2mb^{2}y+(b^{2}m^{2}-a^{2}b^{2}k^{2})=0

 

设弦与椭圆相交的两点中点是(x_{0},y_{0}),利用求根公式可得:

x_{0}=-\frac{kma^{2}}{b^{2}+a^{2}k^{2}}

y_{0}=\frac{mb^{2}}{b^{2}+a^{2}k^{2}}

其中b^{2}+a^{2}k^{2}\neq 0

\frac{y_{0}}{x_{0}}=-\frac{mb^{2}}{kma^{2}}=-\frac{b^{2}}{ka^{2}},与m无关。因此(x_{0},y_{0})轨迹是过椭圆中心的直线y=-\frac{b^{2}}{ka^{2}}x在椭圆内部的部分。

 

如何判断段弧属于同一椭圆呢?作者利用3段不同象限的弧,两两组合求出椭圆的中心点,若中心点距离小于阈值则认为这三段弧可能属于同一椭圆。

这里使用2段弧为例,寻找椭圆中心点。

  1. 弧a和弧b按照坐标排定方向,弧a分为\overset{\frown}{L_{a}M_{a}}\overset{\frown}{M_{a}R_{a}},弧b分为\overset{\frown}{L_{b}M_{b}}\overset{\frown}{M_{b}R_{b}}

  2. 取出弧b的前半段\overset{\frown}{L_{b}M_{b}}和弧a,连接L_{a}M_{b}作为基准弦,利用夹逼方法求出一组平行弦。

  3. 利用这组平行弦的中点计算出中点质心,中点连线斜率和斜率中值。

  4. 同理利用弧a的后半段\overset{\frown}{M_{a}R_{a}}和弧b,连接M_{a}R_{b}作为基准弦,利用夹逼方法求出一组平行弦。

  5. 利用这组平行弦的中点计算出中点质心,中点连线斜率和斜率中值。

  6. 两组平行弦中点连线的交点为椭圆中心C。

 

一下结合代码分析如何利用两段弧计算出椭圆中心。

        //快速计算中心点
	void GetFastCenter(vector<Point>& e1, vector<Point>& e2, EllipseData& data);
GetFastCenter(vector<Point>& e1, vector<Point>& e2, EllipseData& data)
{
	data.isValid = true;//数据有效

	unsigned size_1 = unsigned(e1.size());
	unsigned size_2 = unsigned(e2.size());

	unsigned hsize_1 = size_1 >> 1;//求取size的1/2 将弧分为2个部分
	unsigned hsize_2 = size_2 >> 1;

	Point& med1 = e1[hsize_1];//chenglm 2017.3.6    弧的中点
	Point& med2 = e2[hsize_2];

	Point2f M12, M34;//平行弦中点的质心
	float q2, q4;//平行弦中点构造的直线的斜率

	{
		// First to second  第一段弧e1 到  第二段弧e2

		// Reference slope //参考斜率    e1的起点 和 e2中点连线

		float dx_ref = float(e1[0].x - med2.x);
		float dy_ref = float(e1[0].y - med2.y);

		if (dy_ref == 0) dy_ref = 0.00001f;//防止0/0

		float m_ref = dy_ref / dx_ref;//平行弦的参考斜率
		data.ra = m_ref;

		// Find points with same slope as reference//找到一组与参考弦平行的弦
		vector<Point2f> med;// 平行弦的中点
		med.reserve(hsize_2);

		unsigned minPoints = (_uNs < hsize_2) ? _uNs : hsize_2;//  最小点数

                //针对e2一半弧是抽样还是遍历,由参数_uNs 决定,目的是在确保有足够精度的同时还要确保计算速度
		vector<uint> indexes(minPoints);//e2上参与计算的点
		if (_uNs < hsize_2)//抽样   
		{
			unsigned iSzBin = hsize_2 / unsigned(_uNs);//抽样间隔
			unsigned iIdx = hsize_2 + (iSzBin / 2);

			for (unsigned i = 0; i<_uNs; ++i)
			{
				indexes[i] = iIdx;
				iIdx += iSzBin;
			}
		}
		else//遍历
		{
			for (int i = 0; i < hsize_2; i++)
				indexes[i] = i + hsize_2;

			//iota(indexes.begin(), indexes.end(), hsize_2);

		}


                //e2前半段上的点和e1上的点连线,求取斜率。是否和参考弦平行
		for (uint ii = 0; ii<minPoints; ++ii)
		{
			uint i = indexes[ii];

			float x1 = float(e2[i].x);
			float y1 = float(e2[i].y);

			uint begin = 0;
			uint end = size_1 - 1;

			float xb = float(e1[begin].x);
			float yb = float(e1[begin].y);
			float res_begin = ((xb - x1) * dy_ref) - ((yb - y1) * dx_ref);//chenglm 2017.3,6(x1-x2)/x   =    (y1-y2)/y
			                                                                // (x1-x2)*y=(y1-y2)*x
			int sign_begin = sgn(res_begin);
			if (sign_begin == 0)//利用夹逼方法求解    已知:f(x)在区间D上连续,区间D上有x1<x2,f(x1)>0,f(x2)<0;必然存在x3在(x1,x2)之间,使f(x3)=0;
			{
				//found
				med.push_back(Point2f((xb + x1)* 0.5f, (yb + y1)* 0.5f));
				continue;
			}

			float xe = float(e1[end].x);
			float ye = float(e1[end].y);
			float res_end = ((xe - x1) * dy_ref) - ((ye - y1) * dx_ref);
			int sign_end = sgn(res_end);
			if (sign_end == 0)
			{
				//found
				med.push_back(Point2f((xe + x1)* 0.5f, (ye + y1)* 0.5f));
				continue;
			}

			if ((sign_begin + sign_end) != 0)//如果不是整数解,解必然存在与这两个点之间
			{
				continue;
			}

			uint j = (begin + end) >> 1;//二分法

			while (end - begin > 2)
			{
				float x2 = float(e1[j].x);
				float y2 = float(e1[j].y);
				float res = ((x2 - x1) * dy_ref) - ((y2 - y1) * dx_ref);
				int sign_res = sgn(res);

				if (sign_res == 0)
				{
					//found
					med.push_back(Point2f((x2 + x1)* 0.5f, (y2 + y1)* 0.5f));
					break;
				}

				if (sign_res + sign_begin == 0)
				{
					sign_end = sign_res;
					end = j;
				}
				else
				{
					sign_begin = sign_res;
					begin = j;
				}
				j = (begin + end) >> 1;
			}

			med.push_back(Point2f((e1[j].x + x1)* 0.5f, (e1[j].y + y1)* 0.5f));
		}

		if (med.size() < 2)//两段弧线没有更多的平行线,无法构成圆或者椭圆
		{
			data.isValid = false;//数据无效
			return;
		}

		q2 = GetMedianSlope(med, M12, data.Sa);//求取平均斜率     平行弦中点连线    必过中心点
	}


	{   //同理
		// Second to first

		// Reference slope
		float dx_ref = float(med1.x - e2[0].x);
		float dy_ref = float(med1.y - e2[0].y);

		if (dy_ref == 0) dy_ref = 0.00001f;

		float m_ref = dy_ref / dx_ref;
		data.rb = m_ref;

		// Find points with same slope as reference
		vector<Point2f> med;
		med.reserve(hsize_1);

		uint minPoints = (_uNs < hsize_1) ? _uNs : hsize_1;

		vector<uint> indexes(minPoints);
		if (_uNs < hsize_1)
		{
			unsigned iSzBin = hsize_1 / unsigned(_uNs);
			unsigned iIdx = hsize_1 + (iSzBin / 2);

			for (unsigned i = 0; i<_uNs; ++i)
			{
				indexes[i] = iIdx;
				iIdx += iSzBin;
			}
		}
		else
		{
			for (int i = 0; i < hsize_1; i++)
				indexes[i] = i + hsize_1;

			//iota(indexes.begin(), indexes.end(), hsize_1);

		}


		for (uint ii = 0; ii<minPoints; ++ii)
		{
			uint i = indexes[ii];

			float x1 = float(e1[i].x);
			float y1 = float(e1[i].y);

			uint begin = 0;
			uint end = size_2 - 1;

			float xb = float(e2[begin].x);
			float yb = float(e2[begin].y);
			float res_begin = ((xb - x1) * dy_ref) - ((yb - y1) * dx_ref);
			int sign_begin = sgn(res_begin);
			if (sign_begin == 0)
			{
				//found
				med.push_back(Point2f((xb + x1)* 0.5f, (yb + y1)* 0.5f));
				continue;
			}

			float xe = float(e2[end].x);
			float ye = float(e2[end].y);
			float res_end = ((xe - x1) * dy_ref) - ((ye - y1) * dx_ref);
			int sign_end = sgn(res_end);
			if (sign_end == 0)
			{
				//found
				med.push_back(Point2f((xe + x1)* 0.5f, (ye + y1)* 0.5f));
				continue;
			}

			if ((sign_begin + sign_end) != 0)
			{
				continue;
			}

			uint j = (begin + end) >> 1;

			while (end - begin > 2)
			{
				float x2 = float(e2[j].x);
				float y2 = float(e2[j].y);
				float res = ((x2 - x1) * dy_ref) - ((y2 - y1) * dx_ref);
				int sign_res = sgn(res);

				if (sign_res == 0)
				{
					//found
					med.push_back(Point2f((x2 + x1)* 0.5f, (y2 + y1)* 0.5f));
					break;
				}

				if (sign_res + sign_begin == 0)
				{
					sign_end = sign_res;
					end = j;
				}
				else
				{
					sign_begin = sign_res;
					begin = j;
				}
				j = (begin + end) >> 1;
			}

			med.push_back(Point2f((e2[j].x + x1)* 0.5f, (e2[j].y + y1)* 0.5f));
		}

		if (med.size() < 2)
		{
			data.isValid = false;
			return;
		}
		q4 = GetMedianSlope(med, M34, data.Sb);//求取平均斜率     平行弦中点连线    必过中心点
	}


	if (q2 == q4)//斜率相同,无法相交   (此处可设置阈值,如果斜率相近,可忽略????)
	{
		data.isValid = false;
		return;
	}

	float invDen = 1 / (q2 - q4);
	data.Cab.x = (M34.y - q4*M34.x - M12.y + q2*M12.x) * invDen;
	data.Cab.y = (q2*M34.y - q4*M12.y + q2*q4*(M12.x - M34.x)) * invDen;//中心点
	data.ta = q2;//一中心连线斜率
	data.tb = q4;
	data.Ma = M12;//一中心连线过的一点
	data.Mb = M34;
}

利用3段弧估算出椭圆中心、平行弦斜率等参数后,进一步交叉计算求出更准确的椭圆中心。设椭圆中心点(x_{0},y_{0})(x_{0},y_{0}),一组平行弦中点的方程为y-y_{0}=k(x-x_{0}),将上一步获取的4组参数代入,两两联立方程求解可以得到4组椭圆中心点(x_{0},y_{0})(x_{1},y_{1})(x_{2},y_{2})(x_{3},y_{3}),再加上已经计算得出的2组椭圆中心点(x_{4},y_{4})(x_{5},y_{5}),最后加上2组椭圆中心点的均值(x_{6},y_{6})共计7组椭圆中心点。分别求x和y的中值为最终椭圆中心。

作者用GetCenterCoordinates函数用来进一步估算椭圆中心点。

        //获取椭圆中心点的坐标
	Point2f GetCenterCoordinates(EllipseData& data_ij, EllipseData& data_ik);

3.3.2估算椭圆倾角和长短轴

计算长轴2A、短轴2B、椭圆倾角\varphi并拟和椭圆,作者使用了FindEllipses函数。

        //估算椭圆
	void FindEllipses( Point2f& center,VP& edge_i,VP& edge_j,VP& edge_k,
		EllipseData& data_ij,EllipseData& data_ik,
		vector<Ellipse_Common>& ellipses
		);
  1. 利用平行弦斜率和平行弦中点连线斜率,计算长轴短轴比例值N=\frac{B}{A},和椭圆倾角\varphi

    证明过程有些繁琐,大体思路是椭圆的标准方程和参考方程之间的互换。详细过程请参考链接证明过程和作者论文中的3.2.3章节。

  2. 利用椭圆弧上点坐标计算A和B。这里是利用椭圆一般方程和标准方程之间的转换。
  3. 拟和椭圆,给椭圆评分。评分准则:(1)边界点在拟和椭圆上的比例。(2)弧的端点距离与参数(A+B)的比值,意味着弧尽可能长,椭圆离心率尽可能大。
FindEllipses(Point2f& center,VP& edge_i,VP& edge_j,VP& edge_k,
	EllipseData& data_ij,EllipseData& data_ik,
	vector<Ellipse_Common> &ellipses
	)
{
	// Find ellipse parameters

	// 0-initialize accumulators
	memset(accN, 0, sizeof(int)*ACC_N_SIZE);
	memset(accR, 0, sizeof(int)*ACC_R_SIZE);
	memset(accA, 0, sizeof(int)*ACC_A_SIZE);


	// Get size of the 4 vectors of slopes (2 pairs of arcs)
	int sz_ij1 = int(data_ij.Sa.size());
	int sz_ij2 = int(data_ij.Sb.size());
	int sz_ik1 = int(data_ik.Sa.size());
	int sz_ik2 = int(data_ik.Sb.size());

	// Get the size of the 3 arcs
	size_t sz_ei = edge_i.size();
	size_t sz_ej = edge_j.size();
	size_t sz_ek = edge_k.size();

	// Center of the estimated ellipse
	float a0 = center.x;
	float b0 = center.y;


	// Estimation of remaining parameters
	// Uses 4 combinations of parameters. See Table 1 and Sect [3.2.3] of the paper.

	//chenglm 2017.3.7
	/*
	直接引用论文 ON USING DIRECTIONAL INFORMATION FOR PARAMETER SPACE DECOMPOSITION IN ELLIPSE DETECTION
	应该是利用椭圆标准方程和椭圆参考方程相互转换

	*/



	{
		float q1 = data_ij.ra;
		float q3 = data_ik.ra;
		float q5 = data_ik.rb;

		for (int ij1 = 0; ij1 < sz_ij1; ++ij1)
		{
			float q2 = data_ij.Sa[ij1];

			float q1xq2 = q1*q2;

			for (int ik1 = 0; ik1 < sz_ik1; ++ik1)
			{
				float q4 = data_ik.Sa[ik1];

				float q3xq4 = q3*q4;

				// See Eq. [13-18] in the paper

				float a = (q1xq2 - q3xq4);
				float b = (q3xq4 + 1)*(q1 + q2) - (q1xq2 + 1)*(q3 + q4);
				float Kp = (-b + sqrt(b*b + 4 * a*a)) / (2 * a);
				float zplus = ((q1 - Kp)*(q2 - Kp)) / ((1 + q1*Kp)*(1 + q2*Kp));

				if (zplus >= 0.0f)
				{
					continue;
				}

				float Np = sqrt(-zplus);
				float rho = atan(Kp);
				int rhoDeg;
				if (Np > 1.f)
				{
					Np = 1.f / Np;
					rhoDeg = cvRound((rho * 180 / CV_PI) + 180) % 180; // [0,180)
				}
				else
				{
					rhoDeg = cvRound((rho * 180 / CV_PI) + 90) % 180; // [0,180)
				}

				int iNp = cvRound(Np * 100); // [0, 100]

				if (0 <= iNp	&& iNp < ACC_N_SIZE &&
					0 <= rhoDeg	&& rhoDeg < ACC_R_SIZE
					)
				{
					++accN[iNp];	// Increment N accumulator
					++accR[rhoDeg];	// Increment R accumulator
				}
			}


			for (int ik2 = 0; ik2 < sz_ik2; ++ik2)
			{
				float q4 = data_ik.Sb[ik2];

				float q5xq4 = q5*q4;

				// See Eq. [13-18] in the paper

				float a = (q1xq2 - q5xq4);
				float b = (q5xq4 + 1)*(q1 + q2) - (q1xq2 + 1)*(q5 + q4);
				float Kp = (-b + sqrt(b*b + 4 * a*a)) / (2 * a);
				float zplus = ((q1 - Kp)*(q2 - Kp)) / ((1 + q1*Kp)*(1 + q2*Kp));

				if (zplus >= 0.0f)
				{
					continue;
				}

				float Np = sqrt(-zplus);
				float rho = atan(Kp);
				int rhoDeg;
				if (Np > 1.f)//旋转超过90度
				{
					Np = 1.f / Np;
					rhoDeg = cvRound((rho * 180 / CV_PI) + 180) % 180; // [0,180)
				}
				else//旋转不超过90度
				{
					rhoDeg = cvRound((rho * 180 / CV_PI) + 90) % 180; // [0,180)
				}

				int iNp = cvRound(Np * 100); // [0, 100]

				if (0 <= iNp	&& iNp < ACC_N_SIZE &&
					0 <= rhoDeg	&& rhoDeg < ACC_R_SIZE
					)
				{
					++accN[iNp];		// Increment N accumulator
					++accR[rhoDeg];		// Increment R accumulator
				}
			}

		}
	}


	{
		float q1 = data_ij.rb;
		float q3 = data_ik.rb;
		float q5 = data_ik.ra;

		for (int ij2 = 0; ij2 < sz_ij2; ++ij2)
		{
			float q2 = data_ij.Sb[ij2];

			float q1xq2 = q1*q2;

			for (int ik2 = 0; ik2 < sz_ik2; ++ik2)
			{
				float q4 = data_ik.Sb[ik2];

				float q3xq4 = q3*q4;

				// See Eq. [13-18] in the paper

				float a = (q1xq2 - q3xq4);
				float b = (q3xq4 + 1)*(q1 + q2) - (q1xq2 + 1)*(q3 + q4);
				float Kp = (-b + sqrt(b*b + 4 * a*a)) / (2 * a);
				float zplus = ((q1 - Kp)*(q2 - Kp)) / ((1 + q1*Kp)*(1 + q2*Kp));

				if (zplus >= 0.0f)
				{
					continue;
				}

				float Np = sqrt(-zplus);
				float rho = atan(Kp);
				int rhoDeg;
				if (Np > 1.f)
				{
					Np = 1.f / Np;
					rhoDeg = cvRound((rho * 180 / CV_PI) + 180) % 180; // [0,180)
				}
				else
				{
					rhoDeg = cvRound((rho * 180 / CV_PI) + 90) % 180; // [0,180)
				}

				int iNp = cvRound(Np * 100); // [0, 100]

				if (0 <= iNp	&& iNp < ACC_N_SIZE &&
					0 <= rhoDeg	&& rhoDeg < ACC_R_SIZE
					)
				{
					++accN[iNp];		// Increment N accumulator
					++accR[rhoDeg];		// Increment R accumulator
				}
			}


			for (int ik1 = 0; ik1 < sz_ik1; ++ik1)
			{
				float q4 = data_ik.Sa[ik1];

				float q5xq4 = q5*q4;

				// See Eq. [13-18] in the paper

				float a = (q1xq2 - q5xq4);
				float b = (q5xq4 + 1)*(q1 + q2) - (q1xq2 + 1)*(q5 + q4);
				float Kp = (-b + sqrt(b*b + 4 * a*a)) / (2 * a);
				float zplus = ((q1 - Kp)*(q2 - Kp)) / ((1 + q1*Kp)*(1 + q2*Kp));

				if (zplus >= 0.0f)
				{
					continue;
				}

				float Np = sqrt(-zplus);
				float rho = atan(Kp);
				int rhoDeg;
				if (Np > 1.f)
				{
					Np = 1.f / Np;
					rhoDeg = cvRound((rho * 180 / CV_PI) + 180) % 180; // [0,180)
				}
				else
				{
					rhoDeg = cvRound((rho * 180 / CV_PI) + 90) % 180; // [0,180)
				}

				int iNp = cvRound(Np * 100); // [0, 100]

				if (0 <= iNp	&& iNp < ACC_N_SIZE &&
					0 <= rhoDeg	&& rhoDeg < ACC_R_SIZE
					)
				{
					++accN[iNp];		// Increment N accumulator
					++accR[rhoDeg];		// Increment R accumulator
				}
			}

		}
	}


	// Find peak in N and K accumulator
	int iN = FindMaxN(accN);
	int iK = FindMaxK(accR);

	// Recover real values
	float fK = float(iK);
	float Np = float(iN) * 0.01f;
	float rho = fK * float(CV_PI) / 180.f;	//deg 2 rad
	float Kp = tan(rho);

	// Estimate A. See Eq. [19 - 22] in Sect [3.2.3] of the paper




	//chenglm    2017.3.7   此处是标准椭圆方程和一般椭圆方程转换
	/*
	椭圆一般方程:
	A*x^2+B*x*y+C*y^2+D*x+E*y+F=0;
	变形为:
	A*(x-x0)^2+B*(x-x0)*(y-y0)+C*(y-y0)^2+f=0
	令:x'=x-x0    y'=y-y0
	x0和b0是椭圆中心

	则:A*x'^2+B*x'*y'+C*y'^2+f=0

	标准椭圆方程:(x^2)/(a^2)+(y^2)/(b^2)=1


	若椭圆旋转角度为    -θ

	x=x'*cosθ-y'*sinθ
	y=x'*sinθ+y'*cosθ

	上式变形为:
	((x'*cosθ-y'sinθ)^2)/(a^2)+((x'*sinθ+y'cosθ)^2)/(b^2)=1

	化简后得到:
	(a^2*sinθ^2+b^2*cosθ^2)*x'+(a^2*cosθ^2+b^2*sinθ^2)*y'+2(a^2-b^2)*sinθ*cosθ*x'*y'-a^2*b^2=0


	A=a^2*sinθ^2+b^2*cosθ^2
	B=2(a^2-b^2)*sinθ*cosθ
	C=a^2*cosθ^2+b^2*sinθ^2
	f=-a^2*b^2

	*/



	//chenglm 2017.3.7   此处的变形
	/*

	step 1
	Np=A*B
	Kp = tan(rho);
	sk = 1.f / sqrt(Kp*Kp + 1.f)=cos(rho);
	sin(rho)=Kp*sk
	cos(rho)=1.f / sqrt(Kp^2 + 1.f)
	step2

	x0=x'*cosθ-y'*sinθ
	=(x-a0)*sk+(((y-b0)*Kp)*sk)

	y0=x'*sinθ+y'*cosθ
	=(((x-a0)*Kp)*sk)+(y-b0)*sk

	step3
	椭圆标准方程

	(x0^2)/(a^2)+(y0^2)/(b^2)=1
	(x0^2*b^2)/(a^2)+y0^2=(a^2*b^2)/(a^2)

	用Np=b/a 得到

	x0^2*Np^2+y0^2=a^2*Np^2

	用a=Ax/cos(rho)得到

	x0^2*Np^2+y0^2=Ax^2*Np^2/(cos(rho)^2)
	x0^2*Np^2+y0^2=Ax^2*Np^2*(1+Kp^2)
	Ax=sqrt((x0^2*Np^2+y0^2)/(Np^2*(1+Kp^2)))

	*/




	for (ushort l = 0; l < sz_ei; ++l)
	{
		Point& pp = edge_i[l];
		float sk = 1.f / sqrt(Kp*Kp + 1.f);
		float x0 = ((pp.x - a0) * sk) + (((pp.y - b0)*Kp) * sk);
		float y0 = -(((pp.x - a0) * Kp) * sk) + ((pp.y - b0) * sk);
		float Ax = sqrt((x0*x0*Np*Np + y0*y0) / ((Np*Np)*(1.f + Kp*Kp)));
		int A = cvRound(abs(Ax / cos(rho)));
		if ((0 <= A) && (A < ACC_A_SIZE))
		{
			++accA[A];
		}
	}

	for (ushort l = 0; l < sz_ej; ++l)
	{
		Point& pp = edge_j[l];
		float sk = 1.f / sqrt(Kp*Kp + 1.f);
		float x0 = ((pp.x - a0) * sk) + (((pp.y - b0)*Kp) * sk);
		float y0 = -(((pp.x - a0) * Kp) * sk) + ((pp.y - b0) * sk);
		float Ax = sqrt((x0*x0*Np*Np + y0*y0) / ((Np*Np)*(1.f + Kp*Kp)));
		int A = cvRound(abs(Ax / cos(rho)));
		if ((0 <= A) && (A < ACC_A_SIZE))
		{
			++accA[A];
		}
	}




	for (ushort l = 0; l < sz_ek; ++l)
	{
		Point& pp = edge_k[l];
		float sk = 1.f / sqrt(Kp*Kp + 1.f);
		float x0 = ((pp.x - a0) * sk) + (((pp.y - b0)*Kp) * sk);
		float y0 = -(((pp.x - a0) * Kp) * sk) + ((pp.y - b0) * sk);
		float Ax = sqrt((x0*x0*Np*Np + y0*y0) / ((Np*Np)*(1.f + Kp*Kp)));
		int A = cvRound(abs(Ax / cos(rho)));
		if ((0 <= A) && (A < ACC_A_SIZE))
		{
			++accA[A];
		}
	}

	// Find peak in A accumulator
	int A = FindMaxA(accA);
	float fA = float(A);

	// Find B value. See Eq [23] in the paper
	float fB = abs(fA * Np);

	// Got all ellipse parameters!
	Ellipse_Common ell(a0, b0, fA, fB, fmod(rho + float(CV_PI)*2.f, float(CV_PI)));//确保发布的角度在  0~CV_PI  之间,而不是  -CV_PI/2~CV_PI/2  之间

	//	Toc(3); //estimation
	//	Tac(4); //validation

	// Get the score. See Sect [3.3.1] in the paper

	// Find the number of edge pixel lying on the ellipse

    //统计边界上的点落在拟和椭圆上的比例  
	float _cos = cos(-ell._rad);
	float _sin = sin(-ell._rad);

	float invA2 = 1.f / (ell._a * ell._a);
	float invB2 = 1.f / (ell._b * ell._b);

	float invNofPoints = 1.f / float(sz_ei + sz_ej + sz_ek);
	int counter_on_perimeter = 0;

	for (ushort l = 0; l < sz_ei; ++l)
	{
		float tx = float(edge_i[l].x) - ell._xc;
		float ty = float(edge_i[l].y) - ell._yc;
		float rx = (tx*_cos - ty*_sin);
		float ry = (tx*_sin + ty*_cos);

		float h = (rx*rx)*invA2 + (ry*ry)*invB2;
		if (abs(h - 1.f) < _fDistanceToEllipseContour)
		{
			++counter_on_perimeter;
		}
	}

	for (ushort l = 0; l < sz_ej; ++l)
	{
		float tx = float(edge_j[l].x) - ell._xc;
		float ty = float(edge_j[l].y) - ell._yc;
		float rx = (tx*_cos - ty*_sin);
		float ry = (tx*_sin + ty*_cos);

		float h = (rx*rx)*invA2 + (ry*ry)*invB2;
		if (abs(h - 1.f) < _fDistanceToEllipseContour)
		{
			++counter_on_perimeter;
		}
	}

	for (ushort l = 0; l < sz_ek; ++l)
	{
		float tx = float(edge_k[l].x) - ell._xc;
		float ty = float(edge_k[l].y) - ell._yc;
		float rx = (tx*_cos - ty*_sin);
		float ry = (tx*_sin + ty*_cos);

		float h = (rx*rx)*invA2 + (ry*ry)*invB2;
		if (abs(h - 1.f) < _fDistanceToEllipseContour)
		{
			++counter_on_perimeter;
		}
	}

	//no points found on the ellipse
	if (counter_on_perimeter <= 0)
	{
		//		Toc(4); //validation
		return;
	}


	// Compute score   计算椭圆得分
	float score = float(counter_on_perimeter) * invNofPoints;
	if (score < _fMinScore)
	{
		//	Toc(4); //validation
		return;
	}

	// Compute reliability
	// this metric is not described in the paper, mostly due to space limitations.
	// The main idea is that for a given ellipse (TD) even if the score is high, the arcs
	// can cover only a small amount of the contour of the estimated ellipse.
	// A low reliability indicate that the arcs form an elliptic shape by chance, but do not underlie
	// an actual ellipse. The value is normalized between 0 and 1.
	// The default value is 0.4.

	// It is somehow similar to the "Angular Circumreference Ratio" saliency criteria
	// as in the paper:
	// D. K. Prasad, M. K. Leung, S.-Y. Cho, Edge curvature and convexity
	// based ellipse detection method, Pattern Recognition 45 (2012) 3204-3221.


    //计算每段弧端点距离与长短轴比例关系
    //即弧尽可能长,椭圆离心率尽可能大

	float di, dj, dk;
	{
		Point2f p1(float(edge_i[0].x), float(edge_i[0].y));
		Point2f p2(float(edge_i[sz_ei - 1].x), float(edge_i[sz_ei - 1].y));
		p1.x -= ell._xc;
		p1.y -= ell._yc;
		p2.x -= ell._xc;
		p2.y -= ell._yc;
		Point2f r1((p1.x*_cos - p1.y*_sin), (p1.x*_sin + p1.y*_cos));
		Point2f r2((p2.x*_cos - p2.y*_sin), (p2.x*_sin + p2.y*_cos));
		di = abs(r2.x - r1.x) + abs(r2.y - r1.y);
	}
	{
		Point2f p1(float(edge_j[0].x), float(edge_j[0].y));
		Point2f p2(float(edge_j[sz_ej - 1].x), float(edge_j[sz_ej - 1].y));
		p1.x -= ell._xc;
		p1.y -= ell._yc;
		p2.x -= ell._xc;
		p2.y -= ell._yc;
		Point2f r1((p1.x*_cos - p1.y*_sin), (p1.x*_sin + p1.y*_cos));
		Point2f r2((p2.x*_cos - p2.y*_sin), (p2.x*_sin + p2.y*_cos));
		dj = abs(r2.x - r1.x) + abs(r2.y - r1.y);
	}
	{
		Point2f p1(float(edge_k[0].x), float(edge_k[0].y));
		Point2f p2(float(edge_k[sz_ek - 1].x), float(edge_k[sz_ek - 1].y));
		p1.x -= ell._xc;
		p1.y -= ell._yc;
		p2.x -= ell._xc;
		p2.y -= ell._yc;
		Point2f r1((p1.x*_cos - p1.y*_sin), (p1.x*_sin + p1.y*_cos));
		Point2f r2((p2.x*_cos - p2.y*_sin), (p2.x*_sin + p2.y*_cos));
		dk = abs(r2.x - r1.x) + abs(r2.y - r1.y);
	}

	// This allows to get rid of thick edges
	float rel = min(1.f, ((di + dj + dk) / (3 * (ell._a + ell._b))));

	if (rel < _fMinReliability)
	{
		//		Toc(4); //validation
		return;
	}

	// Assign the new score!
	ell._score = (score + rel) * 0.5f;

	// The tentative detection has been confirmed. Save it!
	ellipses.push_back(ell);

	//	Toc(4); // Validation
}

小结:估算椭圆中心利用到了椭圆的性质:一组平行弦中点连线必过椭圆中心。交叉验算椭圆中心时用到了robust Theil-Sen estimator。在拟和椭圆参数是用到了椭圆标准方程、一般方程和参数方程的相互转换。

3.4椭圆聚类

按照以下步骤给椭圆聚类:

  1. 椭圆按照得分由高到低排序。
  2. 椭圆的中心点C、长轴A、短轴B和倾角\varphi接近的椭圆归为一类;用得分最高的椭圆代表这一类。要注意的是A\approx B时,此时可能是圆被识别为椭圆了,倾角\varphi不再作为分类依据。

 

至此椭圆检测完成。

 

 

4总结

 

 

 

 

 

 

 

 

评论 12
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值