imtransfpor.m

图像对齐（inage registration）工具箱

function [imret1,imret2,pr2,pc2,slu2,slv2]=imtransform(M1,M2,par,intyp,culpc,culpr);
% transform the image (translation + rotation) using the intyp transform
%image rotate about the left top corner of the image, if the angle is for example pi/4, 
%it rotates counterclockwise
%par are [cos(angle) sin(angle) shift_left_right shift_up_down ]
% shift: if the number is positive, shift is up or left, negative down or right
% pr and pc correspond to the position of [1,1] pixel according to the coordinates

vel2=size(M2);
vel1=size(M1);

slx=[1+culpc-1 1+culpc-1       vel2(2)+culpc-1 vel2(2)+culpc-1];
rdy=[1+culpr-1 vel2(1)+culpr-1 1+culpr-1       vel2(1)+culpr-1];

stx=(culpc+ vel2(2)+culpc-1)/2;
sty=(culpr+ vel2(1)+culpr-1)/2;


s=sqrt(par(1)^2+par(2)^2);
u1=par(1)/s;
u2=par(2)/s;
u3=1/s;

%*****************************

p1=u3*(u1*par(3)-u2*par(4));
p2=u3*(u2*par(3)+u1*par(4));



par(3)=p1;
par(4)=p2;
%*****************************


slu2= par(1)*(slx-par(3))+par(2)*(rdy-par(4));
slv2=-par(2)*(slx-par(3))+par(1)*(rdy-par(4));

sttx= par(1)*(stx-par(3))+par(2)*(sty-par(4));
stty=-par(2)*(stx-par(3))+par(1)*(sty-par(4));


d1u=floor(min(slu2));
d1v=floor(min(slv2));

d2u=ceil(max(slu2));
d2v=ceil(max(slv2));

[yu,xv]=meshgrid(d1u:1:d2u,d1v:1:d2v);

xxu=u3*(u1*yu-u2*xv)+par(3);
yyv=u3*(u2*yu+u1*xv)+par(4);

maska=find(xxu>(vel2(2)+culpc-1)|xxu<(1+culpc-1)|yyv>(vel2(1)+culpr-1)|yyv<(1+culpr-1));

trans=interp2((culpc):(vel2(2)+culpc-1),(culpr):(vel2(1)+culpr-1),M2,xxu,yyv,intyp);
trans(maska)=0;
vel3=size(trans);

if d1v<(culpr) 
   idy=1:vel3(1);
   idu=(1:vel1(1))+culpr-d1v;
   pr2=d1v;
else
   idy=(d1v:d2v)-culpr+1;
   idu=1:vel1(1);
   %   pr=d1v-culpr+1;
   pr2=culpr;
end;
sx= max([vel1(1)+culpr-d1v ,d2v-culpr+1,vel3(1),vel1(1)]);

if d1u < (culpc)
   idx=1:vel3(2);
   idv=(1:vel1(2))+culpc-d1u;
   pc2=d1u;
else
   idx=(d1u:d2u)-culpc+1;
   idv=1:vel1(2);
   pc2=culpc;
%   pc=d1u-culpc+1;
end;

sy=max([-d1u+culpc+vel1(2),d2u-culpc+1,vel3(2),vel1(2)]);

imret1(sx,sy)=0;
imret2(sx,sy)=0;

imret1(idy,idx)=trans;
imret2(idu,idv)=M1;