《基于双经度模型的鱼眼图像畸变校正方法》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核心上发表的  还可以有笔误？