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,
%   % 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

```