Code covered by the BSD License  

Highlights from
Active Shape Model (ASM) and Active Appearance Model (AAM)

image thumbnail

Active Shape Model (ASM) and Active Appearance Model (AAM)

by

Dirk-Jan Kroon (view profile)

 

16 Feb 2010 (Updated )

Cootes 2D/3D Active Shape & Appearance Model for automatic image object segmentation and recognition

image_interpolation_backward(Iin,Tlocal,Interpolation,Boundary,ImageSize)
function Iout = image_interpolation_backward(Iin,Tlocal,Interpolation,Boundary,ImageSize)
% This function is used to transform an 2D image, in a backwards way with an
% transformation image.
%
%  Iout = image_interpolation_backward(Iin,T,Interpolation,Boundary,ImageSize)
%
%
% inputs,
%	   Iin : 2D/3D greyscale or color input image
%	   T : (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<5), ImageSize=[size(Iin,1) size(Iin,2)]; end

if((ndims(Iin)==3)&&(size(Tlocal,ndims(Tlocal))==2)), lo=3; else lo=1; end

if(size(Tlocal,ndims(Tlocal))==3), Is3D=true; else Is3D=false; end
 

switch(lower(Interpolation))
    case 'nearest'
        xBas0=round(Tlocal(:,1));
        yBas0=round(Tlocal(:,2));
        if(Is3D)
            zBas0=round(Tlocal(:,3));
        end
    case 'bilinear'
        xBas0=floor(Tlocal(:,1));
        yBas0=floor(Tlocal(:,2));
        if(Is3D)
            zBas0=floor(Tlocal(:,3));
        end
        xBas1=xBas0+1;
        yBas1=yBas0+1;
        if(Is3D)
            zBas1=zBas0+1;
        end
        % Linear interpolation constants (percentages)
        tx=Tlocal(:,1)-xBas0;
        ty=Tlocal(:,2)-yBas0;
        if(Is3D)
            tz=Tlocal(:,3)-zBas0;
        end
        
        if(~Is3D)
            perc0=(1-tx).*(1-ty);
            perc1=(1-tx).*ty;
            perc2=tx.*(1-ty);
            perc3=tx.*ty;
        else
            perc0=(1-tx).*(1-ty).*(1-tz);
            perc1=(1-tx).*(1-ty).*tz;
            perc2=(1-tx).*ty.*(1-tz);
            perc3=(1-tx).*ty.*tz;
            perc4=tx.*(1-ty).*(1-tz);
            perc5=tx.*(1-ty).*tz;
            perc6=tx.*ty.*(1-tz);
            perc7=tx.*ty.*tz;
        end
        
    case 'bicubic'
        if(Is3D)
            error('image_interpolation_backward:inputs','tricubic not supported');
        end
        xBas0=floor(Tlocal(:,1));
        yBas0=floor(Tlocal(:,2));
        tx=Tlocal(:,1)-xBas0;
        ty=Tlocal(:,2)-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));
        if(Is3D)
            check_zBas0=(zBas0<0)|(zBas0>(size(Iin,3)-1));
        end
        xBas0=min(max(xBas0,0),size(Iin,1)-1);
        yBas0=min(max(yBas0,0),size(Iin,2)-1);
        if(Is3D)
            zBas0=min(max(zBas0,0),size(Iin,3)-1);
        end
         
    case 'bilinear'
        check_xBas0=(xBas0<0)|(xBas0>(size(Iin,1)-1));
        check_yBas0=(yBas0<0)|(yBas0>(size(Iin,2)-1));
        if(Is3D)
            check_zBas0=(zBas0<0)|(zBas0>(size(Iin,3)-1));
        end
        check_xBas1=(xBas1<0)|(xBas1>(size(Iin,1)-1));
        check_yBas1=(yBas1<0)|(yBas1>(size(Iin,2)-1));
        if(Is3D)
            check_zBas1=(zBas1<0)|(zBas1>(size(Iin,3)-1));
        end
        xBas0=min(max(xBas0,0),size(Iin,1)-1);
        yBas0=min(max(yBas0,0),size(Iin,2)-1);
        if(Is3D)
            zBas0=min(max(zBas0,0),size(Iin,3)-1);
        end
        xBas1=min(max(xBas1,0),size(Iin,1)-1);
        yBas1=min(max(yBas1,0),size(Iin,2)-1);
        if(Is3D)
            zBas1=min(max(zBas1,0),size(Iin,3)-1);
        end
    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

if(~Is3D)
    Iout=zeros([ImageSize(1:2) lo]);
else
    Iout=zeros(ImageSize);
end

for i=1:lo; % Loop incase of RGB
    if(~Is3D)
        Iin_one=Iin(:,:,i);
    else
        Iin_one=Iin;
    end
    switch(lower(Interpolation))
        case 'nearest'
            % Get the intensities
            if(~Is3D)
                intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1));
            else
                intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1)+zBas0*size(Iin,1)*size(Iin,2));
            end
            
            % Set pixels outside the image
            switch(lower(Boundary))
                case 'zero'
                    if(~Is3D)
                        intensity_xyz0(check_xBas0|check_yBas0)=0;
                    else
                        intensity_xyz0(check_xBas0|check_yBas0|check_zBas0)=0;
                    end
                otherwise
            end
            
            % Combine the weighted neighbour pixel intensities
            Iout_one=intensity_xyz0;
        case 'bilinear'
            % Get the intensities
            if(~Is3D)
                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));
            else
                intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1)+zBas0*size(Iin,1)*size(Iin,2));
                intensity_xyz1=Iin_one(1+xBas0+yBas0*size(Iin,1)+zBas1*size(Iin,1)*size(Iin,2));
                intensity_xyz2=Iin_one(1+xBas0+yBas1*size(Iin,1)+zBas0*size(Iin,1)*size(Iin,2));
                intensity_xyz3=Iin_one(1+xBas0+yBas1*size(Iin,1)+zBas1*size(Iin,1)*size(Iin,2));
                intensity_xyz4=Iin_one(1+xBas1+yBas0*size(Iin,1)+zBas0*size(Iin,1)*size(Iin,2));
                intensity_xyz5=Iin_one(1+xBas1+yBas0*size(Iin,1)+zBas1*size(Iin,1)*size(Iin,2));
                intensity_xyz6=Iin_one(1+xBas1+yBas1*size(Iin,1)+zBas0*size(Iin,1)*size(Iin,2));
                intensity_xyz7=Iin_one(1+xBas1+yBas1*size(Iin,1)+zBas1*size(Iin,1)*size(Iin,2));
            end
            
            % Set pixels outside the image
            switch(lower(Boundary))
                case 'zero'
                    if(~Is3D)
                        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;
                    else
                        intensity_xyz0(check_xBas0|check_yBas0|check_zBas0)=0;
                        intensity_xyz1(check_xBas0|check_yBas0|check_zBas1)=0;
                        intensity_xyz2(check_xBas0|check_yBas1|check_zBas0)=0;
                        intensity_xyz3(check_xBas0|check_yBas1|check_zBas1)=0;
                        intensity_xyz4(check_xBas1|check_yBas0|check_zBas0)=0;
                        intensity_xyz5(check_xBas1|check_yBas0|check_zBas1)=0;
                        intensity_xyz6(check_xBas1|check_yBas1|check_zBas0)=0;
                        intensity_xyz7(check_xBas1|check_yBas1|check_zBas1)=0;
                    end
                otherwise
            end
            
            % Combine the weighted neighbour pixel intensities
            if(~Is3D)
                Iout_one=intensity_xyz0.*perc0+intensity_xyz1.*perc1+intensity_xyz2.*perc2+intensity_xyz3.*perc3;
            else
                Iout_one=intensity_xyz0.*perc0+intensity_xyz1.*perc1+intensity_xyz2.*perc2+intensity_xyz3.*perc3+...
                    intensity_xyz4.*perc4+intensity_xyz5.*perc5+intensity_xyz6.*perc6+intensity_xyz7.*perc7;
            end
        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





Contact us