问题:
我们在拥有一系列散列的点(x1,y1),(x2,y2)... (xm,ym),这些点在一条直线附近,通过点拟合直线。
我在工程中是要拟合一系列线段,其实一条线段就对应着两个要拟合的点,算法上稍有区别,原理完全一致。
思路一:
用OpenCV自带的最小二乘法拟合
函数:cvFitLine()
void cvFitLine( const CvArr* points, int dist_type, double param, double reps, double aeps, float* line );
points
2D 或 3D 点集,32-比特整数或浮点数坐标 ,存放输入点。
Char* cvSeqPush(CvSeq* seq,void* element=NULL)
功能:存放元素到序列尾部。
dist_type
拟合的距离类型.
dist_type=CV_DIST_L2 (L2): ρ(r)=r2/2 (最小二乘法,我们要用的)
dist_type=CV_DIST_L1 (L1): ρ(r)=r
dist_type=CV_DIST_L12 (L1-L2): ρ(r)=2?[sqrt(1+r2/2) - 1]
dist_type=CV_DIST_FAIR (Fair): ρ(r)=C2?[r/C - log(1 + r/C)], C=1.3998
dist_type=CV_DIST_WELSCH (Welsch): ρ(r)=C2/2?[1 - exp(-(r/C)2)], C=2.9846
dist_type=CV_DIST_HUBER (Huber): ρ(r)= r2/2, if r < C; C?(r-C/2), otherwise; C=1.345
param
对某些距离的数字参数,如果是 0, 则选择某些最优值
reps, aeps
半径 (坐标原点到直线的距离) 和角度的精度,一般设为0.01。
line
输出的直线参数。2D 拟合情况下,它是包含 4 个浮点数的数组 (vx, vy, x0, y0),其中 (vx, vy) 是线的单位向量而 (x0, y0) 是线上的某个点.
对 3D 拟合,它是包含 6 个浮点数的数组 (vx, vy, vz, x0, y0, z0), 其中 (vx, vy, vz) 是线的单位向量,而 (x0, y0, z0) 是线上某点。
float *line = new float[4];
CvMemStorage* storage = cvCreateMemStorage(0);
//往点序列中存放需要参与拟合直线的点
CvSeq* point_seq = cvCreateSeq( CV_32FC2, sizeof(CvSeq), sizeof(CvPoint2D32f), storage );
for(int j=0;j<temp.size();j++)
{
cvSeqPush(point_seq, &cvPoint2D32f(temp[j].stPot().x,temp[j].stPot().y));
cvSeqPush(point_seq, &cvPoint2D32f(temp[j].enPot().x,temp[j].enPot().y));
}
cvFitLine(point_seq,CV_DIST_L2,0,0.01,0.01,line); //CV_DIST_L2表示最小二乘法
//line[0],line[1]为x,y的单位方向向量
//line[2],line[3]为直线经过某点的X,Y值
//正好组成参数方程
//X=line[0]t+line[2]
//Y=line[1]t+line[3]
//根据line[4]计算自己需要的直线或线段
Point A(0,line[3]-(line[1]*line[2]/line[0]));
Point B(src->width,(src->width-line[2])*line[1]/line[0]+line[3]);
cvLine(src,A,B,CV_RGB(255,0,0),3,8,0);
cvClearSeq(point_seq);
cvReleaseMemStorage(&storage);
思路二:
自己写算法,但是需要注意,openCV中原点在屏幕的左上角,与我们数学中的直角坐标系不太一样。
//最小二乘法拟合直线
void line_fit(vector<MyLine>& h,IplImage* src)
{
int n = 2*h.size();
int k; //目标直线斜率
int b; //目标直线截距
vector<MyLine>::iterator it = h.begin();
int sumx=0,sumy=0,sumxy=0,sumxsq=0;
while(it != h.end())
{
sumx += it->stPot().x;
sumx += it->enPot().x;
sumy += it->stPot().y;
sumy += it->enPot().y;
sumxy += it->stPot().x*it->stPot().y;
sumxy += it->enPot().x*it->enPot().y;
sumxsq += it->stPot().x*it->stPot().x;
sumxsq += it->enPot().x*it->enPot().x;
it++;
}
if(sumxsq == (sumx*sumx/n))
k = atan2(0,1.0);
else
k = (sumxy-((sumx*sumy)/n))/(sumxsq-(sumx*sumx/n));
b = (sumy-k*sumx)/n;
cvLine(src,Point(0,b),Point(src->width,k*src->width+b),CV_RGB(255,0,0),1,8,0);
}