曲線擬合的最小二乘法(基於OpenCV實現)
阿新 • • 發佈:2019-01-02
在科學實驗資料處理中,往往要根據一組給定的實驗資料,求出自變數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;
}
擬合後的結果如下:
其中左邊是離散的一些點,右邊綠色的線是擬合後的直線。
(5.8.1)
這裡是線性無關的函式族,假定在上給出一組資料,以及對應的一組權,這裡為權係數,要求使最小,其中
(5.8.2)
這就是最小二乘逼近,得到的擬合曲線為y=s(x),這種方法稱為曲線擬合的最小二乘法.
(5.8.2)中實際上是關於的多元函式,求I的最小值就是求多元函式I的極值,由極值必要條件,可得
(5.8.3)
根據內積定義(見第三章)引入相應帶權內積記號
則(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;
}
擬合後的結果如下:
其中左邊是離散的一些點,右邊綠色的線是擬合後的直線。