www.pudn.com > SURF-based-image-stitching.rar > affine_warp.m, change:2010-09-06,size:9721b


function Iout=affine_warp(Iin,M,mode) 
% Affine transformation function (Rotation, Translation, Resize) 
% This function transforms a volume with a 3x3 transformation matrix  
% 
% Iout=affine_warp(Iin,Minv,mode) 
% 
% inputs, 
%   Iin: The input image 
%   Minv: The (inverse) 3x3 transformation matrix 
%   mode: If 0: linear interpolation and outside pixels set to nearest pixel 
%            1: linear interpolation and outside pixels set to zero 
%            (cubic interpolation only support by compiled mex file) 
%            2: cubic interpolation and outsite pixels set to nearest pixel 
%            3: cubic interpolation and outside pixels set to zero 
% 
% output, 
%   Iout: The transformed image 
% 
% example, 
%   % Read image 
%   I=im2double(imread('lenag2.png')) 
%   % Make a transformation matrix 
%   M=make_transformation_matrix([2 3],2,[1.0 1.1]); 
%   % Transform the image 
%   Iout=affine_warp(I,M,0) 
%   % Show the image 
%   figure, imshow(Iout); 
% 
% Function is written by D.Kroon University of Twente (February 2009) 
   
% Make all x,y indices 
[x,y]=ndgrid(0:size(Iin,1)-1,0:size(Iin,2)-1); 
 
% Calculate center of the image 
% mean= size(Iin)/2; 
% Make center of the image coordinates 0,0 
%xd=x-mean(1);  
%yd=y-mean(2); 
xd=x; 
yd=y; 
 
% Calculate the Transformed coordinates 
Tlocalx = mean(1) + M(1,1) * xd + M(1,2) *yd + M(1,3) * 1; 
Tlocaly = mean(2) + M(2,1) * xd + M(2,2) *yd + M(2,3) * 1; 
 
switch(mode) 
	case 0 
		Interpolation='bilinear'; 
		Boundary='replicate'; 
	case 1 
		Interpolation='bilinear'; 
		Boundary='zero'; 
	case 2 
		Interpolation='bicubic'; 
		Boundary='replicate'; 
	otherwise 
		Interpolation='bicubic'; 
		Boundary='zero'; 
end 
Iout=image_interpolation(Iin,Tlocalx,Tlocaly,Interpolation,Boundary); 
 
function Iout = image_interpolation(Iin,Tlocalx,Tlocaly,Interpolation,Boundary,ImageSize) 
% This function is used to transform an 2D image, in a backwards way with an 
% transformation image. 
% 
%   Iout = image_interpolation(Iin,Tlocalx,Tlocaly,Interpolation,Boundary,ImageSize) 
% 
% inputs, 
%	   Iin : 2D greyscale or color input image 
%	   Tlocalx,Tlocaly : (Backwards) Transformation images for all image pixels 
%	   Interpolation: 
%       'nearest'    - nearest-neighbor interpolation 
%       'bilinear'   - bilinear interpolation 
%       'bicubic'    - cubic interpolation; the default method 
%      Boundary: 
%       'zero'       - outside input image are implicilty assumed to be zero 
%       'replicate'  - Input array values outside the bounds of the array 
%                      are assumed to equal the nearest array border value 
%	(optional) 
%	   ImageSize:    - Size of output image 
% outputs, 
%  	   Iout : The transformed image 
% 
% Function is written by D.Kroon University of Twente (September 2010) 
 
if(~isa(Iin,'double')), Iin=double(Iin); end 
if(nargin<6), ImageSize=[size(Iin,1) size(Iin,2)]; end 
if(ndims(Iin)==2), lo=1; else lo=3; end 
 
switch(lower(Interpolation)) 
    case 'nearest' 
        xBas0=round(Tlocalx); 
        yBas0=round(Tlocaly); 
    case 'bilinear' 
        xBas0=floor(Tlocalx); 
        yBas0=floor(Tlocaly); 
        xBas1=xBas0+1; 
        yBas1=yBas0+1; 
         
        % Linear interpolation constants (percentages) 
        tx=Tlocalx-xBas0; 
        ty=Tlocaly-yBas0; 
        perc0=(1-tx).*(1-ty); 
        perc1=(1-tx).*ty; 
        perc2=tx.*(1-ty); 
        perc3=tx.*ty; 
    case 'bicubic' 
        xBas0=floor(Tlocalx); 
        yBas0=floor(Tlocaly); 
        tx=Tlocalx-xBas0; 
        ty=Tlocaly-yBas0; 
         
        % Determine the t vectors 
        vec_tx0= 0.5; vec_tx1= 0.5*tx; vec_tx2= 0.5*tx.^2; vec_tx3= 0.5*tx.^3; 
        vec_ty0= 0.5; vec_ty1= 0.5*ty; vec_ty2= 0.5*ty.^2;vec_ty3= 0.5*ty.^3; 
         
        % t vector multiplied with 4x4 bicubic kernel gives the to q vectors 
        vec_qx0= -1.0*vec_tx1 + 2.0*vec_tx2 - 1.0*vec_tx3; 
        vec_qx1=  2.0*vec_tx0 - 5.0*vec_tx2 + 3.0*vec_tx3; 
        vec_qx2=  1.0*vec_tx1 + 4.0*vec_tx2 - 3.0*vec_tx3; 
        vec_qx3= -1.0*vec_tx2 + 1.0*vec_tx3; 
         
        vec_qy0= -1.0*vec_ty1 + 2.0*vec_ty2 - 1.0*vec_ty3; 
        vec_qy1=  2.0*vec_ty0 - 5.0*vec_ty2 + 3.0*vec_ty3; 
        vec_qy2=  1.0*vec_ty1 + 4.0*vec_ty2 - 3.0*vec_ty3; 
        vec_qy3= -1.0*vec_ty2 + 1.0*vec_ty3; 
         
        % Determine 1D neighbour coordinates 
        xn0=xBas0-1; xn1=xBas0; xn2=xBas0+1; xn3=xBas0+2; 
        yn0=yBas0-1; yn1=yBas0; yn2=yBas0+1; yn3=yBas0+2; 
    otherwise 
        error('image_interpolation:inputs','unknown interpolation method'); 
end 
 
% limit indexes to boundaries 
switch(lower(Interpolation)) 
    case 'nearest' 
        check_xBas0=(xBas0<0)|(xBas0>(size(Iin,1)-1)); 
        check_yBas0=(yBas0<0)|(yBas0>(size(Iin,2)-1)); 
        xBas0=min(max(xBas0,0),size(Iin,1)-1); 
        yBas0=min(max(yBas0,0),size(Iin,2)-1); 
    case 'bilinear' 
        check_xBas0=(xBas0<0)|(xBas0>(size(Iin,1)-1)); 
        check_yBas0=(yBas0<0)|(yBas0>(size(Iin,2)-1)); 
        check_xBas1=(xBas1<0)|(xBas1>(size(Iin,1)-1)); 
        check_yBas1=(yBas1<0)|(yBas1>(size(Iin,2)-1)); 
        xBas0=min(max(xBas0,0),size(Iin,1)-1); 
        yBas0=min(max(yBas0,0),size(Iin,2)-1); 
        xBas1=min(max(xBas1,0),size(Iin,1)-1); 
        yBas1=min(max(yBas1,0),size(Iin,2)-1); 
    case 'bicubic' 
        check_xn0=(xn0<0)|(xn0>(size(Iin,1)-1)); 
        check_xn1=(xn1<0)|(xn1>(size(Iin,1)-1)); 
        check_xn2=(xn2<0)|(xn2>(size(Iin,1)-1)); 
        check_xn3=(xn3<0)|(xn3>(size(Iin,1)-1)); 
        check_yn0=(yn0<0)|(yn0>(size(Iin,2)-1)); 
        check_yn1=(yn1<0)|(yn1>(size(Iin,2)-1)); 
        check_yn2=(yn2<0)|(yn2>(size(Iin,2)-1)); 
        check_yn3=(yn3<0)|(yn3>(size(Iin,2)-1)); 
        xn0=min(max(xn0,0),size(Iin,1)-1); 
        xn1=min(max(xn1,0),size(Iin,1)-1); 
        xn2=min(max(xn2,0),size(Iin,1)-1); 
        xn3=min(max(xn3,0),size(Iin,1)-1); 
        yn0=min(max(yn0,0),size(Iin,2)-1); 
        yn1=min(max(yn1,0),size(Iin,2)-1); 
        yn2=min(max(yn2,0),size(Iin,2)-1); 
        yn3=min(max(yn3,0),size(Iin,2)-1); 
end 
 
Iout=zeros([ImageSize(1:2) lo]); 
for i=1:lo; % Loop incase of RGB 
    Iin_one=Iin(:,:,i); 
    switch(lower(Interpolation)) 
        case 'nearest' 
            % Get the intensities 
            intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1)); 
             
            % Set pixels outside the image 
            switch(lower(Boundary)) 
                case 'zero' 
                    intensity_xyz0(check_xBas0|check_yBas0)=0; 
                otherwise 
            end 
             
            % Combine the weighted neighbour pixel intensities 
            Iout_one=intensity_xyz0; 
        case 'bilinear' 
            % Get the intensities 
            intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1)); 
            intensity_xyz1=Iin_one(1+xBas0+yBas1*size(Iin,1)); 
            intensity_xyz2=Iin_one(1+xBas1+yBas0*size(Iin,1)); 
            intensity_xyz3=Iin_one(1+xBas1+yBas1*size(Iin,1)); 
             
            % Set pixels outside the image 
            switch(lower(Boundary)) 
                case 'zero' 
                    intensity_xyz0(check_xBas0|check_yBas0)=0; 
                    intensity_xyz1(check_xBas0|check_yBas1)=0; 
                    intensity_xyz2(check_xBas1|check_yBas0)=0; 
                    intensity_xyz3(check_xBas1|check_yBas1)=0; 
                otherwise 
            end 
             
            % Combine the weighted neighbour pixel intensities 
            Iout_one=intensity_xyz0.*perc0+intensity_xyz1.*perc1+intensity_xyz2.*perc2+intensity_xyz3.*perc3; 
        case 'bicubic' 
            % Get the intensities 
            Iy0x0=Iin_one(1+xn0+yn0*size(Iin,1));Iy0x1=Iin_one(1+xn1+yn0*size(Iin,1)); 
            Iy0x2=Iin_one(1+xn2+yn0*size(Iin,1));Iy0x3=Iin_one(1+xn3+yn0*size(Iin,1)); 
            Iy1x0=Iin_one(1+xn0+yn1*size(Iin,1));Iy1x1=Iin_one(1+xn1+yn1*size(Iin,1)); 
            Iy1x2=Iin_one(1+xn2+yn1*size(Iin,1));Iy1x3=Iin_one(1+xn3+yn1*size(Iin,1)); 
            Iy2x0=Iin_one(1+xn0+yn2*size(Iin,1));Iy2x1=Iin_one(1+xn1+yn2*size(Iin,1)); 
            Iy2x2=Iin_one(1+xn2+yn2*size(Iin,1));Iy2x3=Iin_one(1+xn3+yn2*size(Iin,1)); 
            Iy3x0=Iin_one(1+xn0+yn3*size(Iin,1));Iy3x1=Iin_one(1+xn1+yn3*size(Iin,1)); 
            Iy3x2=Iin_one(1+xn2+yn3*size(Iin,1));Iy3x3=Iin_one(1+xn3+yn3*size(Iin,1)); 
             
            % Set pixels outside the image 
            switch(lower(Boundary)) 
                case 'zero' 
                    Iy0x0(check_yn0|check_xn0)=0;Iy0x1(check_yn0|check_xn1)=0; 
                    Iy0x2(check_yn0|check_xn2)=0;Iy0x3(check_yn0|check_xn3)=0; 
                    Iy1x0(check_yn1|check_xn0)=0;Iy1x1(check_yn1|check_xn1)=0; 
                    Iy1x2(check_yn1|check_xn2)=0;Iy1x3(check_yn1|check_xn3)=0; 
                    Iy2x0(check_yn2|check_xn0)=0;Iy2x1(check_yn2|check_xn1)=0; 
                    Iy2x2(check_yn2|check_xn2)=0;Iy2x3(check_yn2|check_xn3)=0; 
                    Iy3x0(check_yn3|check_xn0)=0;Iy3x1(check_yn3|check_xn1)=0; 
                    Iy3x2(check_yn3|check_xn2)=0;Iy3x3(check_yn3|check_xn3)=0; 
                otherwise 
            end 
             
            % Combine the weighted neighbour pixel intensities 
            Iout_one=vec_qy0.*(vec_qx0.*Iy0x0+vec_qx1.*Iy0x1+vec_qx2.*Iy0x2+vec_qx3.*Iy0x3)+... 
                vec_qy1.*(vec_qx0.*Iy1x0+vec_qx1.*Iy1x1+vec_qx2.*Iy1x2+vec_qx3.*Iy1x3)+... 
                vec_qy2.*(vec_qx0.*Iy2x0+vec_qx1.*Iy2x1+vec_qx2.*Iy2x2+vec_qx3.*Iy2x3)+... 
                vec_qy3.*(vec_qx0.*Iy3x0+vec_qx1.*Iy3x1+vec_qx2.*Iy3x2+vec_qx3.*Iy3x3); 
    end 
    Iout(:,:,i)=reshape(Iout_one, ImageSize); 
end