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

最新推荐文章于 2022-07-07 14:49:48 发布

万能的七爷

最新推荐文章于 2022-07-07 14:49:48 发布

阅读量2.6w

点赞数 29

文章标签： opencv 圆检测椭圆检测几何特征

本文链接：https://blog.csdn.net/txlqy9041/article/details/83240897

版权

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

参考论文：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 算法流程图

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

三、算法分析

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

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

将与处理过程剥离开，增加了参数bool _isSmooth。_isSmooth为ture时，PrePeocessing使用高斯滤波，_isSmooth为false时,不使用高斯滤波；根据实际应用场景做不同预处理，预处理后的图像直接检测椭圆。
边界检测自适应canny变换时，参数percent_of_pixels_not_edges 和 threshold_ratio是默认值。将整两个参数放出来可以调整,_dNotEdges, _dThMinRatioMax。
变量、类名称的更改。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) 为凹函数。

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

那么边界上的点如何判定凹凸性呢？图像分别对方向和方向求一阶偏导得到 $\frac{\partial _{I}}{\partial _{x}}$ 和 $\frac{\partial _{I}}{\partial _{y}}$ 。当 $\frac{\partial _{I}}{\partial _{x}} > 0$ 时， $I$ 沿着方向递增，当 $\frac{\partial _{I}}{\partial _{y}} < 0$ 时， $I$ 沿着方向递增（图像原点为左上点，因此轴反向）。 $I$ 沿着方向和方向同时递增（递减）时，当前点判定为凸（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.$

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

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

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段弧为例，寻找椭圆中心点。

弧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}}$ 。
取出弧b的前半段 $\overset{\frown}{L_{b}M_{b}}$ 和弧a，连接 $L_{a}$ 和 $M_{b}$ 作为基准弦，利用夹逼方法求出一组平行弦。
利用这组平行弦的中点计算出中点质心，中点连线斜率和斜率中值。
同理利用弧a的后半段 $\overset{\frown}{M_{a}R_{a}}$ 和弧b，连接 $M_{a}$ 和 $R_{b}$ 作为基准弦，利用夹逼方法求出一组平行弦。
利用这组平行弦的中点计算出中点质心，中点连线斜率和斜率中值。
两组平行弦中点连线的交点为椭圆中心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
		);

利用平行弦斜率和平行弦中点连线斜率，计算长轴短轴比例值 $N=\frac{B}{A}$ ，和椭圆倾角 $\varphi$ 。
证明过程有些繁琐，大体思路是椭圆的标准方程和参考方程之间的互换。详细过程请参考链接证明过程和作者论文中的3.2.3章节。
利用椭圆弧上点坐标计算A和B。这里是利用椭圆一般方程和标准方程之间的转换。
拟和椭圆，给椭圆评分。评分准则：(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。在拟和椭圆参数是用到了椭圆标准方程、一般方程和参数方程的相互转换。