【opencv】鱼眼图像畸变校正——双经度法

5月23日更新代码
参考：基于双经度模型的鱼眼图像畸变矫正方法_魏利胜
http://download.youkuaiyun.com/detail/qq_15947787/9583010
注意：原文中公式14，15，17，18，22，23有误，均未考虑符号问题。双经度法适用与视场接近180°以及小于180°情况

#define Q_R 311
#define X0 392
#define Y0 282

IplImage* Double( IplImage* img)
{
	float R = Q_R;
	IplImage* imgCalib = cvCreateImage( cvSize( 2*R, 2*R ),IPL_DEPTH_8U, 3);		//创建校正后图像大小
	float sita = 0;
	float fi = 0;
	float tempu = 0;
	float tempv = 0;
	float x0 = X0 ;//扩展边界后圆心发生偏移
	float y0 = Y0 ;
	float flag1 =1;
	float flag2 =1;
	CvPoint2D32f xy;
	CvPoint2D32f uv;
	
	for (int i = 0; i < imgCalib->width; i++)
	{
		for (int j = 0; j < imgCalib->height; j++)
		{
			sita = PI - (PI/2/R)*i;
			fi = PI - (PI/2/R)*j;
			tempu = tan(fi)*tan(fi) + 1 + (tan(fi)*tan(fi))/(tan(sita)*tan(sita));
			tempv = tan(sita)*tan(sita) + 1 + (tan(sita)*tan(sita))/(tan(fi)*tan(fi));
			if (fi >PI/2)    //文中公式22，23 ，17，18有符号错误
				flag1 = 1;
			else
				flag1 = -1;

			if (sita <PI/2)
				flag2 = 1;
			else
				flag2 = -1;
			uv.x = flag1*R/sqrt(tempu) + x0;//u-fi 下正上负
			uv.y = flag2*R/sqrt(tempv) + y0;//v-sita 右正左负

			int iu = (int)uv.x;
			int iv = (int)uv.y;
			//!此处取的是左上点，为了双线性方便
			CvPoint3D32f f1;
			CvPoint3D32f f2;
			if (iu>=0&&iv>=0&&iu<img->width&&iv<img->height)
			{
				//双线性插值法
				f1.x = ((uchar*)(img->imageData + img->widthStep*iv))[iu*3] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*iv))[(iu+1)*3] * (1-abs(uv.x-(iu+1)));

				f1.y = ((uchar*)(img->imageData + img->widthStep*iv))[iu*3+1] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*iv))[(iu+1)*3+1] * (1-abs(uv.x-(iu+1)));

				f1.z = ((uchar*)(img->imageData + img->widthStep*iv))[iu*3+2] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*iv))[(iu+1)*3+2] * (1-abs(uv.x-(iu+1)));

				f2.x = ((uchar*)(img->imageData + img->widthStep*(iv+1)))[iu*3] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*(iv+1)))[(iu+1)*3] * (1-abs(uv.x-(iu+1)));

				f2.y = ((uchar*)(img->imageData + img->widthStep*(iv+1)))[iu*3+1] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*(iv+1)))[(iu+1)*3+1] * (1-abs(uv.x-(iu+1)));

				f2.z = ((uchar*)(img->imageData + img->widthStep*(iv+1)))[iu*3+2] * (1-abs(uv.x-iu))+
					((uchar*)(img->imageData + img->widthStep*(iv+1)))[(iu+1)*3+2] * (1-abs(uv.x-(iu+1)));

				((uchar*)(imgCalib->imageData + imgCalib->widthStep*j))[i*3] = f1.x*(1-abs(uv.y-iv))+f2.x*(abs(uv.y-iv));
				((uchar*)(imgCalib->imageData + imgCalib->widthStep*j))[i*3+1] = f1.y*(1-abs(uv.y-iv))+f2.y*(abs(uv.y-iv));
				((uchar*)(imgCalib->imageData + imgCalib->widthStep*j))[i*3+2] = f1.z*(1-abs(uv.y-iv))+f2.z*(abs(uv.y-iv));
			}			
		}
	}
	return imgCalib;
}

可以看到在四个极点图像拉伸现象非常严重。