OpenCV(2.4.4及之前版本)polyfit函数存在潜在的bug

1.Problem description

前几天做tracking点的polyfit时，发现OpenCV在contrib模块下有社区成员贡献的polyfit函数，很高兴地觉得不用自己写代码了，结果那一整个下午和晚上都花在寻找bug上了（因为polyfit的曲线和MATLAB拟合出来的相差非常大，而且非常明显原来自带的polyfit函数有bug，见图1）。

图 1.  默认的polyfit函数结果

2.Introduction to Least Squares Fitting

假设拟合的目标曲线为(kth order) 

那么residual记作

拟合的目标是使residual最小，求上式对各个系数求偏导，并将偏导置为0

...

将上式写成矩阵形式

（3）

可以将上面的矩阵拆分为如下形式

3.Possible reason and Solution

问题其实很简单，当输入x变量的范围过大，且进行高阶拟合时（比如2阶。。。），由于默认的精度(CV_32F)不足导致结果偏差比较大，将数据换成CV_64F后解决掉了。但是不确定的是，在高阶拟合时，如果X变量的输入足够大，CV_64F恐怕精度也不会够吧。

已将问题提交给OpenCV社区 ，问题的详细描述及解决办法见http://code.opencv.org/issues/2887

临时补救的代码如下（其实还是有很多地方可以再优化的，本人比较懒。。。）

原代码在/modules/contrib/src/polyfit.cpp内

CODE

/**@ Function : Polyfit for OpenCV/**@ Auhtor : chouclee/**@ Date : Mar.15/2013**/#include <opencv2/opencv.hpp>typedef double polyfit_type;void mypolyfit2(const cv::Mat& _src_x, const cv::Mat& _src_y, cv::Mat& dst, int order){CV_Assert(_src_x.channels() == 1);dst.create(order + 1, 1, CV_MAKETYPE(DataType<polyfit_type>::depth, 1));//目前简单起见，只做单通道的数据拟合int npoints = _src_x.rows;Mat X = Mat::zeros(order + 1, npoints, CV_MAKETYPE(DataType<polyfit_type>::depth, 1));//(order+1)*npointsMat src_x, src_y;_src_y.convertTo(src_y, CV_MAKETYPE(DataType<polyfit_type>::depth, 1));_src_x.convertTo(src_x,CV_MAKETYPE(DataType<polyfit_type>::depth, 1));polyfit_type* pSrcX = (polyfit_type*)src_x.data;//重写了矩阵X的计算部分polyfit_type* pXData = (polyfit_type*)X.data;int stepX = (int)(X.step/X.elemSize1());for (int y = 0; y < order + 1; y++){for (int x = 0; x < npoints; x++){if (y == 0)pXData[x] = 1;else if (y == 1)pXData[x + stepX] = pSrcX[x];else pXData[x + y*stepX] = pSrcX[x]* pXData[x + (y-1)*stepX];}}Mat X_t, X_inv;transpose(X,X_t);//X_t --> npoints*(order+1)Mat temp = X*X_t;//(order+1)*npoints mul npoints*(order+1) --> (order+1)*(order+1)Mat temp2;invert (temp,temp2);Mat temp3 = temp2*X;//(order+1)*(order+1) mul (order+1)*npointsMat W = temp3*src_y;//(order+1)*npoints mul npoints*(order+1)W.copyTo(dst);}

图2. 重写函数之后的拟合结果

从2张图的结果也可以看出，当X输入比较小的，CV_32和CV_64相比拟合的结果差不多，但是当输入达到较高的level时，结果糟透了～～