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)*npoints

Mat 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)*npoints
Mat W = temp3*src_y;//(order+1)*npoints mul npoints*(order+1)
W.copyTo(dst);
}

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