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.