找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 4937|回复: 2
打印 上一主题 下一主题
收起左侧

张正友Matlab相机标定完整版原代码,亲测可用

[复制链接]
跳转到指定楼层
楼主
张正友的Matlab相机标定的完整版原代码,亲测可行!!!


Matlab源程序如下:
  1. % function Zhang(M,m)
  2. %
  3. % ***********************************************************************************
  4. % *******          A Flexible New Technique for Camera Calibration            *******
  5. % ***********************************************************************************
  6. %
  7. % Note:    M:2*N  m:2*N
  8. % M        point on the model plane, when using M=[X,Y]' ---> M=[X,Y,1]'
  9. % m        M's image, when using                m=[u,v]' ---> m=[u,v,1]' , so that
  10. %          s*m = H*M , with H=A*[r1,r2,t];                  (2)
  11. % H        homography matrix
  12. %
  13. % REF:           "A Flexible New Technique for Camera Calibration"
  14. %           - Zhengyou Zhang
  15. %           - Microsoft Research
  16. %
  17. function Zhang(M,m)

  18. %  M=[X,Y]' ---> M=[X,Y,1]'  ;   m=[u,v]' ---> m=[u,v,1]'
  19.     [rows,npts]=size(M);
  20.     matrixone=ones(1,npts);
  21.     M=[M;matrixone];
  22.     num=size(m,3);
  23.     for i=1:num
  24.         m(3,:,i)=matrixone;
  25.     end
  26. % Estimate the H
  27.     for i=1:num
  28.         H(:,:,i)=homography2d1(M,m(:,:,i))';
  29.     end
  30. % solve the intrinsic parameters matrix A
  31. % A=[alpha_u skewness u0
  32. %    0       alpha_v  v0
  33. %    0       0        1]
  34. % see Appendix B "Extraction of the Intrisic Parameters from Matrix B", P18
  35.     V=[];
  36.     for flag=1:num
  37.         v12(:,:,flag)=[H(1,1,flag)*H(2,1,flag), H(1,1,flag)*H(2,2,flag)+H(1,2,flag)*H(2,1,flag), H(1,2,flag)*H(2,2,flag), H(1,3,flag)*H(2,1,flag)+H(1,1,flag)*H(2,3,flag), H(1,3,flag)*H(2,2,flag)+H(1,2,flag)*H(2,3,flag), H(1,3,flag)*H(2,3,flag)];
  38.         v11(:,:,flag)=[H(1,1,flag)*H(1,1,flag), H(1,1,flag)*H(1,2,flag)+H(1,2,flag)*H(1,1,flag), H(1,2,flag)*H(1,2,flag), H(1,3,flag)*H(1,1,flag)+H(1,1,flag)*H(1,3,flag), H(1,3,flag)*H(1,2,flag)+H(1,2,flag)*H(1,3,flag), H(1,3,flag)*H(1,3,flag)];
  39.         v22(:,:,flag)=[H(2,1,flag)*H(2,1,flag), H(2,1,flag)*H(2,2,flag)+H(2,2,flag)*H(2,1,flag), H(2,2,flag)*H(2,2,flag), H(2,3,flag)*H(2,1,flag)+H(2,1,flag)*H(2,3,flag), H(2,3,flag)*H(2,2,flag)+H(2,2,flag)*H(2,3,flag), H(2,3,flag)*H(2,3,flag)];
  40.         V=[V;v12(:,:,flag);v11(:,:,flag)-v22(:,:,flag)];
  41.     end
  42.     k=V'*V;      
  43.     [u,v,d]=svd(k);
  44.     [e,d2]=eig(k);
  45.     b=d(:,6);
  46.     v0=(b(2)*b(4)-b(1)*b(5))/(b(1)*b(3)-b(2)^2);
  47.     s=b(6)-(b(4)^2+v0*(b(2)*b(4)-b(1)*b(5)))/b(1);
  48.     alpha_u=sqrt(s/b(1));
  49.     alpha_v=sqrt(s*b(1)/(b(1)*b(3)-b(2)^2));
  50.     skewness=-b(2)*alpha_u*alpha_u*alpha_v/s;
  51.     u0=skewness*v0/alpha_u-b(4)*alpha_u*alpha_u/s;
  52.     A=[alpha_u skewness u0
  53.         0      alpha_v  v0
  54.         0      0        1];
  55. % solve k1 k1 and all the extrisic parameters, P6
  56.     D=[];
  57.     d=[];
  58.     Rm=[];
  59.     for flag=1:num
  60.         s=(1/norm(inv(A)*H(1,:,flag)')+1/norm(inv(A)*H(2,:,flag)'))/2;
  61.         rl1=s*inv(A)*H(1,:,flag)';
  62.         rl2=s*inv(A)*H(2,:,flag)';
  63.         rl3=cross(rl1,rl2);
  64.         RL=[rl1,rl2,rl3];
  65.         %%%%%%%%%%%%%%%%%%%%
  66.         % see Appendix C "Approximating a 3*3 matrix by a Rotation Matrix", P19
  67.         [U,S,V] = svd(RL);
  68.         RL=U*V';
  69.         %%%%%%%%%%%%%%%%%%%%
  70.         TL=s*inv(A)*H(3,:,flag)';
  71.         RT=[rl1,rl2,TL];
  72.         XY=RT*M;
  73.         UV=A*XY;
  74.         UV=[UV(1,:)./UV(3,:); UV(2,:)./UV(3,:); UV(3,:)./UV(3,:)];
  75.         XY=[XY(1,:)./XY(3,:); XY(2,:)./XY(3,:); XY(3,:)./XY(3,:)];
  76.         for j=1:npts
  77.             D=[D; ((UV(1,j)-u0)*( (XY(1,j))^2 + (XY(2,j))^2 )) , ((UV(1,j)-u0)*( (XY(1,j))^2 + (XY(2,j))^2 )^2) ; ((UV(2,j)-v0)*( (XY(1,j))^2 + (XY(2,j))^2 )) , ((UV(2,j)-v0)*( (XY(1,j))^2 + (XY(2,j))^2 )^2) ];
  78.             d=[d; (m(1,j,flag)-UV(1,j)) ; (m(2,j,flag)-UV(2,j))];
  79.         end
  80.         r13=RL(1,3);
  81.         r12=RL(1,2);
  82.         r23=RL(2,3);
  83.         Q1=-asin(r13);
  84.         Q2=asin(r12/cos(Q1));
  85.         Q3=asin(r23/cos(Q1));
  86.         [cos(Q2)*cos(Q1)   sin(Q2)*cos(Q1)   -sin(Q1) ; -sin(Q2)*cos(Q3)+cos(Q2)*sin(Q1)*sin(Q3)    cos(Q2)*cos(Q3)+sin(Q2)*sin(Q1)*sin(Q3)  cos(Q1)*sin(Q3) ; sin(Q2)*sin(Q3)+cos(Q2)*sin(Q1)*cos(Q3)    -cos(Q2)*sin(Q3)+sin(Q2)*sin(Q1)*cos(Q3)  cos(Q1)*cos(Q3)];
  87.         R_new=[Q1,Q2,Q3,TL'];
  88.         Rm=[Rm , R_new];
  89.     end
  90. % using function (13), P8
  91.     k=inv(D'*D)*D'*d;
  92. % Complete Maximun Likelihood Estimation, using function (14), P8
  93.     para=[Rm,k(1),k(2),alpha_u,skewness,u0,alpha_v,v0];
  94.     options = optimset('LargeScale','off','LevenbergMarquardt','on');
  95.     [x,resnorm,residual,exitflag,output]  = lsqnonlin( @simon_HHH, para, [],[],options, m, M);
  96. % display the result
  97.     k1=x(num*6+1)
  98.     k2=x(num*6+2)
  99.     A=[x(num*6+3) x(num*6+4) x(num*6+5); 0 x(num*6+6) x(num*6+7); 0,0,1]
复制代码

所有资料51hei提供下载:
zhang's+calibration.zip (7.25 MB, 下载次数: 30)


分享到:  QQ好友和群QQ好友和群 QQ空间QQ空间 腾讯微博腾讯微博 腾讯朋友腾讯朋友
收藏收藏1 分享淘帖 顶 踩
回复

使用道具 举报

沙发
ID:473532 发表于 2019-4-25 14:19 | 只看该作者
请问我运行为啥报错。。。

出错 homography2d1 (line 78)
    options = optimset('LargeScale','off','LevenbergMarquardt','on');

出错 Zhang (line 31)
        H(:,:,i)=homography2d1(M,m(:,:,i))';

出错 test (line 28)
Zhang(M,m)
回复

使用道具 举报

板凳
ID:1011605 发表于 2022-3-19 20:12 | 只看该作者
Wangqq 发表于 2019-4-25 14:19
请问我运行为啥报错。。。

出错 homography2d1 (line 78)

我也是这个问题,请问解决了吗?
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

手机版|小黑屋|51黑电子论坛 |51黑电子论坛6群 QQ 管理员QQ:125739409;技术交流QQ群281945664

Powered by 单片机教程网

快速回复 返回顶部 返回列表