《基于双经度模型的鱼眼图像畸变校正方法》code

网友说现在用得多的是双经度法来校正鱼眼图像，理论部分来自题目说的这个期刊论文。

主要分5步：

1，和期刊论文《一种鱼眼图象到透视投影图象的变换模型》中一样求光学中心和真实半径。。。循环三次求平均值得到鱼眼图像的中心和半径

2，目标图像坐标转为双经度坐标

3，双经度坐标转为球面坐标

4，球面坐标转鱼眼图像坐标

5，双线性插值

其中第4步，有两种投影方式：正交投影和等距投影

第1步：

%《一种鱼眼图象到透视投影图象的变换模型》

% 对于齐次线性方程 A*X =0;当A的秩大于列数时，就需要求解最小二乘解 使用SVD分解矩阵A，[U S V] = svd(A); U 由 A*A'的特征向量组成，V 由

%A'*A的特征向量组成，因此，奇异值矩阵S中最小的奇异值对应的V中的奇异向量即为最小二乘解。

function result=doublejingdu(A)

A=imread('F:\fisheye\14.JPG');

[A,R]=kuaisusaomiao(A,40);

[width,height,kk]=size(A);

%xo=R;

%yo=R;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

imshow(A)

hold on;

[x1 y1]=ginput(1);

plot(x1,y1,'g.');

drawnow;

[x2 y2]=ginput(1);

plot(x2,y2,'g.');

drawnow;

[x3 y3]=ginput(1);

plot(x3,y3,'g.');

drawnow;

[x4 y4]=ginput(1);

plot(x4,y4,'g.');

drawnow;

[x5 y5]=ginput(1);

plot(x5,y5,'g.');

drawnow;

[x6 y6]=ginput(1);

plot(x6,y6,'g.');

drawnow;

[x7 y7]=ginput(1);

plot(x7,y7,'g.');

drawnow;

hold off;

%[x8 y8]=ginput(1);

%plot(x8,y8,'g.');

%drawnow;

%[x9 y9]=ginput(1);

%plot(x9,y9,'g.');

%drawnow;

%[x10 y10]=ginput(1);

%plot(x10,y10,'g.');

%drawnow;

%hold off;

M0=[x1,y1;x2,y2;x3,y3;x4,y4;x5,y5;x6,y6;x7,y7];

M=zeros(7,6);

j=M0(:,2);

i=M0(:,1);

M(:,1)=j.^2;

M(:,2)=2*i.*j;

M(:,3)=i.^2;

M(:,4)=j;

M(:,5)=i;

M(:,6)=1;

[U,S,V]=svd(M);

%选S里最小对应的V即为所求的六个参数A B C D E F

S =

1.0e+007 *

2.9628 0 0 0 0 0

0 0.6137 0 0 0 0

0 0 0.0733 0 0 0

0 0 0 0.0000 0 0

0 0 0 0 0.0000 0

0 0 0 0 0 0.0000

0 0 0 0 0 0

V =

-0.5589 -0.7725 0.3016 0.0002 -0.0004 0.0000

-0.7728 0.3533 -0.5272 -0.0002 -0.0000 0.0000

-0.3007 0.5277 0.7944 -0.0000 -0.0002 0.0000

-0.0002 -0.0003 0.0003 -0.6256 0.7802 -0.0008

-0.0002 0.0001 0.0000 0.7802 0.6256 -0.0001

-0.0000 -0.0000 0.0000 -0.0004 0.0007 1.0000

A=0.0002;

B=-A;

C=0;

D=-0.6256;

E=0.7802;

F=-0.0004;

u1=(C*D-B*E)/2/(B^2-A*C);

v1=(A*E-B*D)/2/(B^2-A*C);

a1=sqrt(((C*D^2-2*B*D*E+A*E^2)/4/(B^2-A*C)+F)*((A+C+sqrt((A-C)^2+4*B^2))/2/(B^2-A*C)));

u1 =

1.9505e+003

v1 =

386.5000

a1 =

0 +1.9277e+003i

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u2=(C*D-B*E)/2/(B^2-A*C);

v2=(A*E-B*D)/2/(B^2-A*C);

a2=sqrt(((C*D^2-2*B*D*E+A*E^2)/4/(B^2-A*C)+F)*((A+C+sqrt((A-C)^2+4*B^2))/2/(B^2-A*C)));

u2 =

2.4340e+003

v2 =

646.8333

a2 =

1.4902e+003

u3=(C*D-B*E)/2/(B^2-A*C);

v3=(A*E-B*D)/2/(B^2-A*C);

a3=sqrt(((C*D^2-2*B*D*E+A*E^2)/4/(B^2-A*C)+F)*((A+C+sqrt((A-C)^2+4*B^2))/2/(B^2-A*C)));

u3 =

337

v3 =

649.8333

a3 =

634.9031

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




选择的是这三条弧线 第一条是最右边这条 第二条是金字塔那条 第三条是柱子那条 我在想是不是不能选第一条这样横着的 或者说我这个图不适合双经度校正？ 第一个算出来的u1 v1和a1那么奇怪 只有u3 v3 a3算出来是看起来比较正常的 还有SVD分解那里 最小S对应的V就是参数A B C D E F 可是比如第一次算出的S有三个都是0那我选哪个呢？？ 奇怪 然后我主要想看后面的效果
我就暂且把第三次算的变化一两个像素作为平均值好了

%final center and R

u0=338;

v0=648;

R=635;

后面的几步：

function resultimg=calibration(fisheyeimg,Rtrue,u0,v0,flag)

%《基于双经度模型的鱼眼图像畸变校正方法》

% Detailed explanation goes here

R=Rtrue+10;

for i=1:2*R

for j=1:2*R

arfa=(pi/2/R)*i;

beta=(pi/2/R)*j;

sita=pi-arfa;

fi=pi-beta;

fenmu1=sqrt((tan(fi))^2+1+((tan(fi))^2/(tan(sita))^2));

x=R/fenmu1;

fenmu2=sqrt((tan(sita))^2+1+((tan(sita))^2/(tan(fi))^2));

y=R/fenmu2;

fenmu3=sqrt(1/((tan(fi))^2)+1+(1/(tan(sita))^2));

z=R/fenmu3;

if(flag==1)

u=x+u0;

v=y+v0;

else

r=R*arccos(z/R);

u=r*x/sqrt(x^2+y^2)+u0;

v=r*y/sqrt(x^2+y^2)+v0;

end

[ER,EG,EB]=doublelinechazhi(fisheyeimg,v,u);

resultimg(i,j,1)=ER;

resultimg(i,j,2)=EG;

resultimg(i,j,3)=EB;

end

end

end

双线性插值自己写的感觉很慢啊：

function [ER,EG,EB]=doublelinechazhi(fisheyeimg,u,v)

%《基于双经度模型的鱼眼图像畸变校正方法》双线性插值

% Detailed explanation goes here

fisheyeimgR=fisheyeimg(:,:,1);

fisheyeimgG=fisheyeimg(:,:,2);

fisheyeimgB=fisheyeimg(:,:,3);

%R channel

FR=(u-floor(u))*(fisheyeimgR(floor(u),ceil(v))-fisheyeimgR(floor(u),floor(v)))+fisheyeimgR(floor(u),floor(v));

GR=(u-floor(u))*(fisheyeimgR(ceil(u),ceil(v))-fisheyeimgR(ceil(u),floor(v)))+fisheyeimgR(ceil(u),floor(v));

ER=(v-floor(v))*(GR-FR)+FR;

%G channel

FG=(u-floor(u))*(fisheyeimgG(floor(u),ceil(v))-fisheyeimgG(floor(u),floor(v)))+fisheyeimgG(floor(u),floor(v));

GG=(u-floor(u))*(fisheyeimgG(ceil(u),ceil(v))-fisheyeimgG(ceil(u),floor(v)))+fisheyeimgG(ceil(u),floor(v));

EG=(v-floor(v))*(GG-FG)+FG;

%B channel

FB=(u-floor(u))*(fisheyeimgB(floor(u),ceil(v))-fisheyeimgB(floor(u),floor(v)))+fisheyeimgB(floor(u),floor(v));

GB=(u-floor(u))*(fisheyeimgB(ceil(u),ceil(v))-fisheyeimgB(ceil(u),floor(v)))+fisheyeimgB(ceil(u),floor(v));

EB=(v-floor(v))*(GB-FB)+FB;

end

太慢了 就不用双线性插值 用最近邻的像素代替不是整数的点 看下效果 结果：


我只能呵呵了 看来我错了

重新看论文 看到一个地方错了 空间直线投影成球面大圆 球面大圆投影成鱼眼图像的椭圆 算出的长轴a只是直径 应该再除以2 才是半径Rtrue 可是光学中心和鱼眼图像中心是一个意思吗 即使改了 出来也不对啊

我重新选一幅图 试试这次出来的S：


这次的S出来分别是3.9489、1.1605、0.2224、0.0003、0.0001、0那么最小的就是最后一个0 它对应的V向量 A=0 、B=0、C=0.0001、D=-0.0093、E=-0.0174、F=0.9998可是这样一来 那个光学中心u0、v0分母都是0了啊 长轴的分母也是0啊？？？为什么 难道选弧线还有什么特别要求？

然后看到http://blog.csdn.net/qq_15947787/article/details/51441128这位兄台也写了双经度法 果然论文里还是有笔误的 只是这个兄台的中心和半径是直接给出的。 所以我也不知道他会不会遇到我这种问题 可是我按照他的改了

先不看左右反了 单纯看校正效果 很差啊 虽然半径和光学中心都有偏差 但也不该这样房子上的柱子还是那么弯 不过可以看出是双极点
四个角拉伸果然严重。哎。。。还是不行。。。话说这个论文是在《电子测量与仪器学报》今年的CSCD核心上发表的 还可以有笔误？