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

曲线拟合的最小二乘法基于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;
}
拟合后的结果如下:
其中左边是离散的一些点，右边绿色的线是拟合后的直线。