Fitting ellipses using RANSAC


Fitting ellipses using RANSAC algorithm

Basic knowledge about ellipse equation http://mathworld.wolfram.com/Ellipse.html

In this tutorial, I am going to use the general quadratic curve which came from the above link. Please note that the below equations are copied from WolframMathWorld. I just put it below just in case you could not access to the above link. But the equation is not pretty (hard to read).

ax^2+2bxy+cy^2+2dx+2fy+g=0

From the above equation, we could find the center of ellipse as:

x_0 = (cd-bf)/(b^2-ac)

y_0 = (af-bd)/(b^2-ac)

the semi-axes lengths are:

a^' = sqrt( ( 2(af^2+cd^2+gb^2-2bdf-acg ) )/( (b^2-ac)[sqrt( (a-c)^2+4b^2)-(a+c)]) )

b^' = sqrt( ( 2(af^2+cd^2+gb^2-2bdf-acg ) )/( (b^2-ac)[-sqrt( (a-c)^2+4b^2)-(a+c)]) )

and the counterclockwise angle of rotation from the x-axis to the major axis of the ellipse is

Φ={0 for b=0 and a<c; 1/2pi for b=0 and a>c; 1/2cot^(-1)( (a-c)/(2b) ) for b!=0 and a<c; pi/2+1/2cot^(-1)( (a-c)/(2b) ) for b!=0 and a>c.

Ellipse implementation

First, create a structure to store information of ellipse as following code:

typedef struct{
  double x,y,majorAxis,minorAxis,angle;
  double f1_x,f1_y,f2_x,f2_y;
} Ellipse;

This is a function which pass the coefficient a,b,c,d,f,g as and input parameter and get ellipse struct back as output. This function will calculate the center of the ellipse, semi-axes, angle and two foci.

void getEllipseParam(double a,double b,double c,double d,double f,double g,Ellipse& ellipse){
  ellipse.x = (c * d - b * f)/(b * b - a * c);
  ellipse.y = (a * f - b * d)/(b * b - a * c);
 
  ellipse.majorAxis = sqrt( (2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g))/((b*b-a*c)*(sqrt((a-c)*(a-c)+4*b*b)-(a+c))));
  ellipse.minorAxis = sqrt( (2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g))/((b*b-a*c)*(sqrt((a-c)*(a-c)+4*b*b)+(a+c))));
 
  ellipse.angle=0;
  if(b == 0 && a < c){
    ellipse.angle = 0;
  }
  else if(b == 0 && a > c){
    ellipse.angle = 90;
  }
  else if(b != 0 && a < c){
    ellipse.angle = 0.5 * acot_d( (a-c)/(2*b) );
  }
  else if(b != 0 && a > c){
    ellipse.angle = 90 + 0.5 * acot_d( (a-c)/(2*b) );
  }
  if(ellipse.minorAxis > ellipse.majorAxis){
    double temp = ellipse.majorAxis;
    ellipse.majorAxis = ellipse.minorAxis;
    ellipse.minorAxis = temp;
    ellipse.angle += 90;
  }
 
  double temp_c;
  if(ellipse.majorAxis > ellipse.minorAxis)
    temp_c = sqrt(ellipse.majorAxis * ellipse.majorAxis - ellipse.minorAxis * ellipse.minorAxis);
  else
    temp_c = sqrt(ellipse.minorAxis * ellipse.minorAxis - ellipse.majorAxis * ellipse.majorAxis);
  ellipse.f1_x = ellipse.x - temp_c * cos(ellipse.angle*PI/180);
  ellipse.f1_y = ellipse.y - temp_c * sin(ellipse.angle*PI/180);
  ellipse.f2_x = ellipse.x + temp_c * cos(ellipse.angle*PI/180);
  ellipse.f2_y = ellipse.y + temp_c * sin(ellipse.angle*PI/180);
}

Function arc-cotan:

double acot_d(double val){
  double acot = atan(1/val);
 
  return acot*180/PI;
}

Function to check whether a given point is inside ellipse or not. This function will find the total distance from the input point to two focus points of ellipse. If the total distance is less than or equal to 2*major axis, this point is inside the ellipse. Otherwise, it is outside ellipse.

bool pointInEllipse(CvPoint point,Ellipse ellipse){
  double dist1 = sqrt((point.x - ellipse.f1_x) * (point.x - ellipse.f1_x) + 
		      (point.y - ellipse.f1_y) * (point.y - ellipse.f1_y));
  double dist2 = sqrt((point.x - ellipse.f2_x) * (point.x - ellipse.f2_x) + 
		      (point.y - ellipse.f2_y) * (point.y - ellipse.f2_y));
  double max;
  if(ellipse.majorAxis > ellipse.minorAxis)
    max = ellipse.majorAxis;
  else
    max = ellipse.minorAxis;
  if(dist1+dist2 <= 2*max)
    return true;
  else
    return false;
}

RANSAC implementation

This part is an RANSAC algorithm. First randomly select 5 points. Then construct a matrix A. Then find SVD of A. The solution is the last column vector of V. The vector V is the coefficient a,b,c,d,f,g. Then use the coefficient to get ellipse information. After having ellipse information, count number of points that are inside of the ellipse. If there is number of points inside ellipse, stop RANSAC algorithm. Otherwise, re-run the RANSAC algorithm.

This is the Function of RANSAC algorithm:

Ellipse fitEllipseRANSAC(vector<CvPoint> points,int &count){
  Ellipse ellipse;
  count=0;
  int index[5];
  bool match=false;
  for(int i=0;i<5;i++){
    do {
      match = false;
      index[i]=rand()%points.size();
      for(int j=0;j<i;j++){
	if(index[i] == index[j]){
	  match=true;
	}
      }
    }
    while(match);
  }
  double aData[] = {
    points[index[0]].x * points[index[0]].x, 2 * points[index[0]].x * points[index[0]].y, points[index[0]].
    y * points[index[0]].y, 2 * points[index[0]].x, 2 * points[index[0]].y,
 
    points[index[1]].x * points[index[1]].x, 2 * points[index[1]].x * points[index[1]].y, points[index[1]].
    y * points[index[1]].y, 2 * points[index[1]].x, 2 * points[index[1]].y,
 
    points[index[2]].x * points[index[2]].x, 2 * points[index[2]].x * points[index[2]].y, points[index[2]].
    y * points[index[2]].y, 2 * points[index[2]].x, 2 * points[index[2]].y,
 
    points[index[3]].x * points[index[3]].x, 2 * points[index[3]].x * points[index[3]].y, points[index[3]].
    y * points[index[3]].y, 2 * points[index[3]].x, 2 * points[index[3]].y,
 
    points[index[4]].x * points[index[4]].x, 2 * points[index[4]].x * points[index[4]].y, points[index[4]].
    y * points[index[4]].y, 2 * points[index[4]].x, 2 * points[index[4]].y };
  CvMat matA=cvMat(5,5,CV_64F,aData);
  CvMat *D,*U,*V;
  D=cvCreateMat(5,5,CV_64F);
  U=cvCreateMat(5,5,CV_64F);
  V=cvCreateMat(5,5,CV_64F);
 
  cvSVD(&matA,D,U,V,CV_SVD_MODIFY_A);
 
  double a,b,c,d,f,g;
  a=cvmGet(V,0,4);
  b=cvmGet(V,1,4);
  c=cvmGet(V,2,4);
  d=cvmGet(V,3,4);
  f=cvmGet(V,4,4);
  g=1;
 
  getEllipseParam(a,b,c,d,f,g,ellipse);
 
  vector<CvPoint>::iterator point_iter;
  if(ellipse.majorAxis > 0 && ellipse.minorAxis > 0){
    for(point_iter=points.begin();point_iter!=points.end();point_iter++){
      CvPoint point = *point_iter;
      if(pointInEllipse(point,ellipse)){
	count++;
      }
    }
  }
 
  return ellipse;
}

Result

This is the result that is quite good:

However, the result might not be as good as expected:



Improvements:

1) As we see, the fitted ellipses tend to be overly large compared to their underneath irregular shapes.  It is caused by the way of evaluating consensus points.  The metric should be changed into computing the error of a new point fitting the model estimated by sample subset , since points to be fitted are contour points  of object rather than their interior points.  

2) I think the matrix should be 5*6 rather than 5*5, since there are actually 6 parameters in a generic ellipse equation. Even though the last coefficient could be set to a constant, say, 1, it should be included and involved in the computation of SVD.


  • 1
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值