椭圆拟合算法实现

椭圆的目标函数：

F(A,B,C,D,E)=XiGeMa(xi^2+A*xiyi+B*yi^2+C*xi+D*yi+E)^2

分别对A,B,C,D,E 求一阶偏导并令其等于0

得到线性方程组：

|A1B1C1D1E1|

|A2B2C2D2E2|

|A3B3C3D3E3|

|A4B4C4D4E4|

|A5B5C5D5E5||A|=|resul1||B|=|resul2||C|=|resul3||D|=|resul4||E|=|resul5|

求得A,B,C,D,E.

椭圆的五个参数：

center.x =(2*B*C-A*D)/(A*A-4*B);

center.y =(2*D-A*D)/(A*A-4*B);

fenzi =2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);

fenmu =(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B)) +1);

femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B)) +1);

long =sqrt(fabs(fenzi/fenmu));

short =sqrt(fabs(fenzi/femmu2));

theta =atan(sqrt((center.x*center.x-center.y*center.y*B)/(center.x*B-center.y*center.y))+0.0001)*180/cv_pi;;*center.x

vectorgetEllipsepar(vectorvec_point)

{

vectorvec_result;

double x3y1=0,x1y3=0,x2y2=0,yyy4=0, xxx3=0,xxx2=0,x2y1=0,yyy3=0,x1y2=0,yyy2=0,x1y1=0,xxx1=0,yyy1=0;

int N =vec_point.size();

for (intm_i=0;m_i

{

double xi =vec_point[m_i].x;

double yi =vec_point[m_i].y;

x3y1+=xi*xi*xi*yi;

x1y3+=xi*yi*yi*yi;

x2y2+=xi*xi*yi*yi;;

yyy4+=yi*yi*yi*yi;

xxx3+=xi*xi*xi;

xxx2+=xi*xi;

x2y1+=xi*xi*yi;

x1y2+=xi*yi*yi;

yyy2+=yi*yi;

x1y1+=xi*yi;

xxx1+=xi;

yyy1+=yi;

yyy3+=yi*yi*yi;

}

long double resul1=-(x3y1);

long double resul2=-(x2y2);

long double resul3=-(xxx3);

long double resul4=-(x2y1);

long double resul5=-(xxx2);

long double B1=x1y3, C1=x2y1, D1=x1y2,

long double B2=yyy4, C2=x1y2, D2=yyy3,

long double B3=x1y2, C3=xxx2, D3=x1y1,

long double B4=yyy3, C4=x1y1, D4=yyy2,

long double B5=yyy2, C5=xxx1, D5=yyy1,

//

CvMat*Ma =cvCreateMat(5,5,CV_64FC1);

CvMat*Md =cvCreateMat(5,1,CV_64FC1);

CvMat*Mb =cvCreateMat(5,1,CV_64FC1);

//

cvmSet(Mb,0,0,resul1);

cvmSet(Mb,1,0,resul2);

cvmSet(Mb,2,0,resul3);

cvmSet(Mb,3,0,resul4);

cvmSet(Mb,4,0,resul5);

cvmSet(Ma,0,0,A1);

cvmSet(Ma,0,1,B1);

cvmSet(Ma,0,2,C1);

cvmSet(Ma,0,3,D1);

cvmSet(Ma,0,4,E1);

cvmSet(Ma,1,0,A2);

cvmSet(Ma,1,1,B2);

cvmSet(Ma,1,2,C2);

cvmSet(Ma,1,3,D2);

cvmSet(Ma,1,4,E2);

cvmSet(Ma,2,0,A3);

cvmSet(Ma,2,1,B3);

cvmSet(Ma,2,2,C3);

cvmSet(Ma,2,3,D3);

cvmSet(Ma,2,4,E3);

cvmSet(Ma,3,0,A4);

cvmSet(Ma,3,1,B4);

cvmSet(Ma,3,2,C4);

cvmSet(Ma,3,3,D4);E1=x1y1, E2=yyy2, E3=xxx1, E4=yyy1, E5=N, A1=x2y2; A2=x1y3; A3=x2y1; A4=x1y2; A5=x1y1;

cvmSet(Ma,3,4,E4);

cvmSet(Ma,4,0,A5);

cvmSet(Ma,4,1,B5);

cvmSet(Ma,4,2,C5);

cvmSet(Ma,4,3,D5);

cvmSet(Ma,4,4,E5);

cvSolve(Ma,Mb, Md);

long double A =cvmGet(Md,0,0);

long double B =cvmGet(Md,1,0);

long double C =cvmGet(Md,2,0);

long double D =cvmGet(Md,3,0);

long double E =cvmGet(Md,4,0);

double XC =(2*B*C-A*D)/(A*A-4*B);

double YC =(2*D-A*D)/(A*A-4*B);

long double fenzi =2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);

long double fenmu =(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B)) +1);

long double femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B)) +1);

double XA =sqrt(fabs(fenzi/fenmu));

double XB =sqrt(fabs(fenzi/femmu2));

double Xtheta =atan(sqrt((XA*XA-XB*XB*B)/(XA*XA*B-XB*XB))+0.0001)*180/3.1415926;vec_result.push_back(XC);

vec_result.push_back(YC);

vec_result.push_back(XA);

vec_result.push_back(XB);

vec_result.push_back(Xtheta);

return vec_result;

}