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

来源：互联网发布：什么是java并发编辑：程序博客网时间：2024/04/29 18:31

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

主要分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核心上发表的还可以有笔误？

0 0