椭圆拟合之最小二乘法

原理：可能以后的篇幅再详细写了，

代码：

QImage inputImage ：输入图像；

float ellipse ： 输出椭圆的各个参数

bool Vertical ： 是否垂直

bool total ： 椭圆是否有上下部分，还是只是椭圆的一段曲线

void Ellipse(QImage *inputImage,float* ellipse,bool Vertical,bool total)

{

vector> XX_mats;

// double x[2];

vector x_mats;

vector y_mats;

int Min_x = 100000,Max_x = 0;

float b = 100000;

float Max_a = 0;

float a;

int y_Min = 0, y_Max = 0;

float X0,Y0;

//初始值

float sum_x = 0;

float sum_y = 0;

int width = inputImage->width();

int height = inputImage->height();

for(int i = 25; i

{

for(int j = 25; j

{

if(QColor(inputImage->pixel(j,i)).red() == 255)

{

if(j>Max_x)

{

Max_x = j;

y_Max = i;

}

if(j

{

Min_x = j;

y_Min = i;

}

sum_y += i;

sum_x += j;

x_mats.push_back(j);

y_mats.push_back(i);

}

}

}

if(total)

{

X0 = sum_x /x_mats.size();

Y0 = sum_y / y_mats.size();

}

else {

X0 = ((Max_x + Min_x) / 2);

Y0 = ((y_Max + y_Min) / 2);

}

cout<

a = sqrtf(powf((Max_x-Min_x),2)+powf((y_Max - y_Min),2)) / 2.0 ;

for(int k=0; k

{

float juli = sqrtf(powf((x_mats[k]-X0),2) + powf((y_mats[k] - Y0),2));

if(juli

{

b = juli;

}

if(juli>Max_a)

{

Max_a = juli;

}

}

//a=Max_a;

float tem;

if(Vertical)

{

tem = a;

a = b;

b = tem;

}

float A = a*a;

float B = b*b;

float sum_0 = 0;

for(int j = 0; j

{

float y = powf((y_mats[j] - Y0),2) / B;

float x = powf((x_mats[j] - X0),2) / A;

sum_0 += powf((1-x-y),2);

}

cout<

float theta =0;

float s1 = 1.0f / x_mats.size();

float learn = 1.0f;

float* hmatrix = new float[x_mats.size()];

memset(hmatrix,0,x_mats.size() * sizeof(float));

//梯度下降

float sum_min = 1000.0;

int countes = 0;

for(int i = 0;i<13140;i++)

{

bool kaiguang = true;

float sum_1 = 0,sum_2 = 0,sum_3 = 0,sum_4 = 0,sum_5 = 0,sum = 0;

for(int j = 0; j

{

float cx = x_mats[j] - X0;

float cy = y_mats[j] - Y0;

float x = cx * cos(theta) - cy * sin(theta);

float y = cx * sin(theta) + cy * cos(theta);

x = (x * x) / powf(a,2);

y = (y * y) / powf(b,2);

//float y = powf((y_mats[j] - Y0),2) / B;

//float x = powf((x_mats[j] - X0),2) / A;

hmatrix[j] = 1 - x -y;

}

for(int j = 0; j

{

float cx = x_mats[j] - X0;

float cy = y_mats[j] - Y0;

float x = cx * cos(theta) - cy * sin(theta);

float y = cx * sin(theta) + cy * cos(theta);

sum_1 += hmatrix[j] * (2) * (x*cos(theta) / powf(a,2) + y*sin(theta) / powf(b,2));

sum_2 += hmatrix[j] * (2) * (x*sin(theta) / powf(a,2) + y*cos(theta) / powf(b,2));

sum_3 += hmatrix[j] * (2) * (powf(x,2) / powf(a,3));

sum_4 += hmatrix[j] * (2) * (powf(y,2) / powf(b,3));

sum_5 += hmatrix[j] * (2) * ((x*(-y)) / powf(a,2) + (y*x) / powf(b,2));

sum += powf(hmatrix[j],2);

}

if(sum

{

sum_min = sum;

kaiguang = false;

}

if(kaiguang)

{

if(fabs(sum - sum_min)<= 0.00001)

{

countes +=1;

if(countes >= 11)

{

break;

}

}

}

sum_1 = learn*s1*sum_1;

sum_2 = learn*s1*sum_2;

sum_3 = learn*s1*sum_3;

sum_4 = learn*s1*sum_4;

sum_5 = learn*s1*sum_5;

// 更新 参数theta

X0 = X0 - sum_1;

Y0 = Y0 - sum_2;

a = a - sum_3;

b = b - sum_4;

theta = theta + sum_5;

cout<

}

/**/

for(int i = 0; iheight()/2;i++)

{

for(int j = 25; jwidth()-10; j++)

{

/*

float y = sqrtf(b*b - ((b*b)/(a*a))*(j-X0)*(j-X0)) + Y0;

y1 = Y0 - sqrtf(b*b - ((b*b)/(a*a))*(j-X0)*(j-X0));

sum_0 += sqrtf(powf((y_mats[j] - y),2.0));

if(y>0&&y<890)

{

inputImage->setPixel(j,int(y+0.3),qRgb(0,255,0));

}

if(y1>0&&y1<900)

inputImage->setPixel(j,int(y1),qRgb(0,255,0));

*/

float cx = j - X0;

float cy = i - Y0;

float x = cx * cos(theta) - cy * sin(theta);

float y = cx * sin(theta) + cy * cos(theta);

x = (x * x) / powf(a,2);

y = (y * y) / powf(b,2);

float result = x + y;

if(abs(result -1)<0.01)

{

inputImage->setPixel(j,i,qRgb(0,255,0));

}

}

}

Show_Image(*inputImage);

/*

*/

ellipse[0] = X0;

ellipse[1] = Y0;

ellipse[2] = a;

ellipse[3] = b;

ellipse[4] = theta;

}

效果图：

误差值：