图像处理中二次曲线拟合

来源：互联网发布：国内网络婚纱品牌排行编辑：程序博客网时间：2024/05/18 00:51

2016/7/16

在一次提取发光管的中心线程序中，由于我们只拍到了断续而弯曲的发光管，所以无法使用光带中心线提取的方法进行提取。

在此背景下，我想到了拟合。之前有学过直线拟合的方法，名为最小二乘法。其基本步骤如下：

(1) 设需要拟合的直线为y=a*x+b。

(2) 首先选取进行拟合的点集，选取方法可以为阈值分割，模板匹配等，设最后选出的点集为。

(3) 求该点集到直线的距离平方和，

(4) 对Sum分别求关于x,y的偏导函数。

(5) 根据偏导求出该距离平方和最小时的a，b值即为拟合的曲线的参数。

为了与之后的二次曲线拟合做对比，我写了下直线拟合的函数，先从源图像中根据阈值分割选取拟合点集（这里我认为灰度值超过45即被选取），再根据上述步骤计算a，b的值。源代码如下：

/*输入为三通道图像*/

/*对图像中的亮点进行直线拟合*/

void cvLineFit1D(IplImage* src_getin)

{

IplImage*src =cvCloneImage(src_getin);

IplImage*image_threshold = cvCreateImage(cvGetSize(src),8,1);

cvCvtColor(src,image_threshold,CV_BGR2GRAY);

cvThreshold(image_threshold,image_threshold,45,255,CV_THRESH_BINARY);

cvShowImage("cvLineFit1D[Threshold]",image_threshold);

//设拟合的二次曲线方程为y=ax+b;

//先求出各点到拟合直线上的距离的平方和;

//求出使得该平方和最小的a,b的值 ;

long long int k1=0;

long long int k2=0;

long long int k3=0;

long long int k4=0;

long long int k5=0;

long long int k6=0;

for (inti=0;i<image_threshold->height;i++)

{

for (intj=0;j<image_threshold->width;j++)

{

if(cvGetReal2D(image_threshold,i,j)==255)

{

k1+=2*i*i;

k2+=2*i;

k3+=2*i*j;

k4+=2*i;

k5+=2;

k6+=2*j;

}

double a,b;

a = (double)(k3*k5-k6*k2)/(k1*k5-k4*k2);

b = (double)(k3-k1*a)/k2;

int bottom = image_threshold->height;

int top = 0;

for(int i =0;i<image_threshold->height;i++)

for(int j =0;j<image_threshold->width;j++)

{

if(cvGetReal2D(image_threshold,i,j) == 255)

{

if(i>top)

top= i;

if(i<bottom)

bottom= i;

}

for (inti=bottom;i<top;i++)

{

int temp_y = a*i+b;

if(temp_y<src->width&&temp_y>0)

{

cvCircle(src,cvPoint(temp_y,i),1,cvScalar(255,255,255));

//cvSet2D(src,i,temp_y,cvScalar(255,255,255));

}

cvShowImage("直线拟合",src);

}

cvReleaseImage(&src);

}

拟合效果如下：

而曲线拟合，则是在直线拟合的基础上，进行相似的偏导计算，只是设所拟合的曲线方程为y=a*x*x+b*x+c

代码如下：

/*对图像中的亮点进行二次曲线拟合;*/

void cvCurveFit2D(IplImage* src)

{

IplImage*image_threshold = cvCreateImage(cvGetSize(src),8,1);

cvCvtColor(src,image_threshold,CV_BGR2GRAY);

cvThreshold(image_threshold,image_threshold,50,255,CV_THRESH_BINARY);

//cvShowImage("cvCurveFit2D[Threshold]",image_threshold);

//设拟合的二次曲线方程为y=ax2+bx+c;

//先求出各点到拟合曲线上的距离的平方和;

//求出使得该平方和最小的a,b,c的值 ;

double a,b,c;

long long int k1=0;

long long int k2=0;

long long int k3=0;

long long int k4=0;

long long int k5=0;

long long int k6=0;

long long int k7=0;

long long int k8=0;

long long int k9=0;

long long int k10=0;

long long int k11=0;

long long int k12=0;

for (inti=0;i<image_threshold->height;i+=5)

{

for (intj=0;j<image_threshold->width;j+=5)

{

if(cvGetReal2D(image_threshold,i,j)==255)

{

//cvCircle(src,cvPoint(j,i),1,cvScalar(0,0,0));

//如果该点为亮点，则将其加入拟合点;

//cout<<"k1="<<k1<<",i="<<i<<",j="<<j<<endl;

k1= k1+(long longint)i*i*i*i; //k1==2x^4

//cout<<"k1="<<k1<<endl;

//getchar();

k2= k2 + (long longint)i*i*i; //k2==2x^3

k3= k3 + (long longint)i*i; //k3==2x^2

k4= k4 + (long longint)i*i*j; //k4==-2y*x^2

k5= k5 + (long longint)i*i*i;

k6= k6 + (long longint)i*i;

k7= k7 + (long longint)i;

k8= k8 + (long longint)i*j;

k9= k9 + (long longint)i*i;

k10= k10 + (long longint)i;

k11= k11 + 1;

k12= k12 + (long longint)j;

}

//根据偏导求出a,b,c的值;

long long int t1 = k1/k3-k5/k7;

long long int t2 = k2/k3-k6/k7;

long long int t3 = k4/k3-k8/k7;

long long int t4 = k1/k3-k9/k11;

long long int t5 = k2/k3-k10/k11;

long long int t6 = k4/k3-k12/k11;

a = (double)(t3*t5-t6*t2)/(t1*t5-t4*t2);

b = (double)(t3-t1*a)/t2;

double b1=(double)(t6-t4*a)/t5;

c = (double)(k4-k1*a-k2*b)/k3;