曲线拟合的最小二乘法（基于OpenCV实现）的，拟合图像中离散点的拟合直线

今天使用拟合的最小二乘法，求出了给定的一组坐标系上的点对最接近的直线的。

　　其具体理论如下：

在科学实验数据处理中，往往要根据一组给定的实验数据

，求出自变量x与因变量y的函数关系

，这是

为待定参数，由于观测数据总有误差，且待定参数ai的数量比给定数据点的数量少(即n＜m)，因此它不同于插值问题.这类问题不要求

通过点

，而只要求在给定点

上的误差

的平方和

最小.当

时，即

　　　　

　　　　　(5.8.1)

这里

是线性无关的函数族，假定在

上给出一组数据

，

以及对应的一组权

，这里

为权系数，要求

使

最小，其中

　　　　

　　　　　　　(5.8.2)

这就是最小二乘逼近，得到的拟合曲线为y=s(x)，这种方法称为曲线拟合的最小二乘法.

　　(5.8.2)中

实际上是关于

的多元函数，求I的最小值就是求多元函数I的极值，由极值必要条件，可得

　　

　　　　(5.8.3)

根据内积定义(见第三章)引入相应带权内积记号

　　　　　　

　　　　　　　(5.8.4)

则(5.8.3)可改写为

　　　　　


这是关于参数

的线性方程组，用矩阵表示为

　　　　　　

　　　　　　(5.8.5)

(5.8.5)称为法方程.当

线性无关，且在点集

上至多只有n个不同零点，则称

在X上满足Haar条件，此时(5.8.5)的解存在唯一(证明见［3］).记(5.8.5)的解为

　　　　　　　　


从而得到最小二乘拟合曲线

　　　　　

　　　　　(5.8.6)

可以证明对

，有

　　　　　


故(5.8.6)得到的

即为所求的最小二乘解.它的平方误差为

　　　　　　　　　　

　　　　　　　　(5.8.7)

均方误差为

　　　　　　　　　　


　　在最小二乘逼近中，若取

，则

，表示为

　　　　　　　

　　　　　　(5.8.8)

此时关于系数

的法方程(5.8.5)是病态方程，通常当n≥3时都不直接取

作为基。

源代码如下：

基于OpenCV

int calAngle(IplImage **lpSrc, IplImage **lpRes, double* angle)

{

long double a,b,q;

long double u11,u12,u21,u22,c1,c2,p;

int i,j;

MYPOINT *points = 0;

IplImage *src = *lpSrc;

IplImage *res = *lpRes;

int step = src->widthStep/sizeof(char); //排列的图像行大小，以字节为单位

uchar *srcData = (uchar*)src->imageData;

int elementCount = 0;

MYPOINT tempPoint;

MYPOINT *interatorPoint;

CvPoint startPoint, endPoint;

res = cvCreateImage( cvGetSize(src), 8, 3 );

memset(res->imageData, 0x00, (src->width)*(src->height));

u11 = 0.0;

u12 = 0.0;

u21 = 0.0;

u22 = 0.0;

c1 = 0.0;

c2 = 0.0;

p = 0.0;

points = (MYPOINT*)malloc((src->height) * sizeof(MYPOINT));

interatorPoint = points;

for (j = 0; j < src->height; ++j)

{

for (i = 0; i < src->width; ++i)

{

if (srcData[j*step + i] == 255)

{

tempPoint.x = i;

tempPoint.y = j;

*interatorPoint = tempPoint;

++interatorPoint;

++elementCount;

break;

}

}

}

interatorPoint = points;

for (i = 0; i < elementCount; i++)//列法方程

{

if (0 == i)

{

startPoint.x = (*interatorPoint).x;

startPoint.y = (*interatorPoint).y;

}

u11+=1.0;

u12+=(*interatorPoint).x;

u22+=(*interatorPoint).x*(*interatorPoint).x;

c1+=(*interatorPoint).y;

c2+=(*interatorPoint).y*(*interatorPoint).x;

++interatorPoint;

//printf("##:%d:/t %lf, %lf, %lf, %lf, %lf%\n", i, u11, u12, u22, c1, c2);

}

u21=u12;

a=(c1*u22-c2*u12)*1.0/(u11*u22-u12*u21);//求a,b的解

b=(c1*u21-c2*u11)*1.0/(u21*u12-u22*u11);

interatorPoint = points;

for(i = 0; i < elementCount; i++)

{

(*interatorPoint).z=a+b*(*interatorPoint).x;

p+=pow(((*interatorPoint).z-(*interatorPoint).y),2);

++interatorPoint;

}

q=sqrt(p);

//输出线性方程

printf("answer:s(x)=%lf+%lf(x)\n",a,b);

endPoint.y = src->height - 1;

endPoint.x = ( endPoint.y -a ) / b;

printf("endPoint:%d, %d\n", endPoint.x, endPoint.y);

//输出均方误差

printf("均方误差:q=%lf\n",q);

*angle = atan(b);

printf("Angle is: %f", *angle*180.0/PI);

cvLine(res, startPoint, endPoint, CV_RGB(0,255,0), 1, 8 );

cvReleaseImage(lpRes);

*lpRes = res;

free(points);

return 0;

}

拟合后的结果如下：

其中左边是离散的一些点，右边绿色的线是拟合后的直线。
原文地址：http://blog.chinaunix.net/u1/58901/showart_2085724.html