Code covered by the BSD License  

Highlights from
ECC image alignment algorithm (image registration)

image thumbnail

ECC image alignment algorithm (image registration)

by

 

15 Apr 2010 (Updated )

This is a Matlab implementation of the ECC image alignment (image registration) algorithm.

ecc(image, template, levels, noi, transform, delta_p_init)
function [results, warp, warpedImage] = ecc(image, template, levels, noi, transform, delta_p_init)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ECC image alignment algorithm
%[RESULTS, WARP, WARPEDiMAGE] = ECC(IMAGE, TEMPLATE, LEVELS, NOI, TRANSFORM, DELTA_P_INIT)
%
% This m-file implements the ECC image alignment algorithm as it is
% presented in the paper "G.D.Evangelidis, E.Z.Psarakis, Parametric Image Alignment
% using Enhanced Correlation Coefficient.IEEE Trans. on PAMI, vol.30, no.10, 2008"
%
% ------------------
% Input variables:
% IMAGE:        the profile needs to be warped in order to be similar to TEMPLATE,
% TEMPLATE:     the profile needs to be reached,
% NOI:          the number of iterations per level; the algorithm is executed
%               (NOI-1) times
% LEVELS:       the number of levels in pyramid scheme (set LEVELS=1 for a
%               non pyramid implementation), the level-index 1
%               corresponds to the highest (original) image resolution
% TRANSFORM:    the image transformation {'translation', 'euclidean', 'affine', 'homography'}
% DELTA_P_INIT: the initial transformation matrix for original images (optional); The identity
%               transformation is the default value (see 'transform initialization'
%               subroutine in the code). In case of affine or euclidean transform, 
%               DELTA_P_INIT must be a 2x3 matrix, in homography case it must be a 3x3 matrix,
%               while with translation transform it must be a 2x1 vector.
%
% For example, to initialize the warp with a rotation by x radians, DELTA_P_INIT must
% be [cos(x) sin(x) 0 ; -sin(x) cos(x) 0] for affinity 
% [cos(x) sin(x) 0 ; -sin(x) cos(x) 0 ; 0 0 1] for homography.
%
%
% Output:
%
% RESULTS:   A struct of size LEVELS x NOI with the following fields:
%
% RESULTS(m,n).warp:     the warp needs to be applied to IMAGE at n-th iteration of m-th level,
% RESULTS(m,n).rho:      the enhanced correlation coefficient value at n-th iteration of m-th level,
% WARP :              the final estimated transformation [usually also stored in RESULTS(1,NOI).WARP ].  
% WARPEDiMAGE:        the final warped image (it should be similar to TEMPLATE).
%
% The first stored .warp and .rho values are due to the initialization. In
% case of poor final alignment results check the warp initialization
% and/or the overlap of the images.
% -------------------
% $ Ver: 1.4, 12/2/2013,  released by Georgios D. Evangelidis, INRIA, FRANCE
% For any comment, please contact the author 
% Email: georgios.evangelidis@inria.fr, evagelid@ceid.upatras.gr
%
% This software is provided "as is" without any kind of warranty. Also, it
% is provided for research purposes only. In any case, please cite the above paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;

break_flag=0; 

if nargin<5
    error('-> Not enough input arguments');
end

transform = lower(transform);
if ~(strcmp(transform,'affine')||strcmp(transform,'euclidean')||strcmp(transform,'homography')||strcmp(transform,'translation'))
    error('-> Not a valid transform string')
end

sZi3 = size(image,3);
sZt3 = size(template,3);
initImage = image;
initTemplate = template;
if sZi3>1
    if ((sZi3==2) || (sZi3>3))
        error('Unknown color image format: check the number of channels');
    else
        image=rgb2gray(uint8(image));
    end
end

if sZt3>1
    if ((sZt3==2) || (sZt3>3))
        error('Unknown color image format: check the number of channels');
    else
        template = rgb2gray(uint8(template));
    end
end

template = double(template);
image = double(image);

%% pyramid images
% The following for-loop creates pyramid images in cells IM and TEMP with varying names
% The variables IM{1} and TEMP{1} are the images with the highest resoltuion

% Smoothing of original images
f = fspecial('gaussian',[7 7],.5);
TEMP{1} = imfilter(template,f);
IM{1} = imfilter(image,f);
TEMP{1} = template;
IM{1} = image;

for nol=2:levels
    IM{nol} = imresize(IM{nol-1},.5);
    TEMP{nol} = imresize(TEMP{nol-1},.5);
end


%% transform initialization

% In case of translation transform the initialiation matrix is of size 2x1:
%  delta_p_init = [p1;
%                  p2]
% In case of affine transform the initialiation matrix is of size 2x3:
%
%  delta_p_init = [p1, p3, p5;
%                  p2, p4, p6]
%
% In case of euclidean transform the initialiation matrix is of size 2x3:
%
%  delta_p_init = [p1, p3, p5;
%                  p2, p4, p6]
%
% where p1=cos(theta), p2 = sin(theta), p3 = -p2, p4 =p1
%
% In case of homography transform the initialiation matrix is of size 3x3:
%  delta_p_init = [p1, p4, p7;
%                 p2, p5, p8;
%                 p3, p6,  1]
if strcmp(transform,'translation')
    nop=2; %number of parameteres
    if nargin==5;
        warp=zeros(2,1);
    else
        if (size(delta_p_init,1)~=2)|(size(delta_p_init,2)~=1)
            error('-> In translation case the size of initialization matrix must be 2x1 ([deltaX;deltaY])');
        else
            warp=delta_p_init;
        end
    end
end
if strcmp(transform,'euclidean')
    nop=3; %number of parameteres
    if nargin==5;
        warp=[1 0 0; 0 1 0; 0 0 0];
    else
        if (size(delta_p_init,1)~=2)||(size(delta_p_init,2)~=3)
            error('-> In euclidean case the size of initialization matrix must be 2x3');
        else
            warp=[delta_p_init;zeros(1,3)];
        end
    end
end
if strcmp(transform,'affine')
    nop=6; %number of parameters
    if nargin==5;
        warp=[1 0 0; 0 1 0; 0 0 0];
    else
        if (size(delta_p_init,1)~=2)|(size(delta_p_init,2)~=3)
            error('-> In affine case the size of initialization matrix must be 2x3');
        else
            warp=[delta_p_init;zeros(1,3)];
        end
    end
end

if strcmp(transform,'homography')
    nop=8; %number of parameteres
    if nargin==5;
        warp=eye(3);
    else
        if (size(delta_p_init,1)~=3)|(size(delta_p_init,2)~=3)
            error('-> In homography case the size of initialization matrix must be 3x3');
        else
            warp=delta_p_init;
            if warp(3,3)~=1
                error('The ninth element of homography must be equal to 1');
            end
        end
    end
end

% in case of pyramid implementation, the initial transformation must be
% appropriately modified
for ii=1:levels-1
    warp=next_level(warp, transform, 0);
end

%% Run ECC algorithm for each level of pyramid
for nol=levels:-1:1
    
    im = IM{nol};
    [vx,vy]=gradient(im);
    
    temp = TEMP{nol};
    
    [A,B]=size(temp);
    % Warning for tiny images
    if prod([A,B])<400
        disp(' -> ECC Warning: The size of images in high pyramid levels is quite small and it may cause errors.');
        disp(' -> To avoid such errors you could try fewer levels or larger images.');
        disp(' -> Press any key to continue.')
        pause
    end
        
    % Define the rectangular Region of Interest (ROI) by nx and ny (you can modify
    % the ROI). Here we just ignore some image margins. Margin is equal to 
    % 5 percent of the mean of [height,width].
    m0=mean([A,B]);
    margin=floor(m0*.05/(2^(nol-1)));
    margin =0;%no-margin - modify these two lines if you want to exclude a margin
    nx=margin+1:B-margin;
    ny=margin+1:A-margin;
    temp=double(temp(ny,nx,:));
    
    %%%%   temp=temp-mean(temp(:)); % zero-mean image; is useful for brightness change compensation, otherwise you can comment this line
    % MODIFIED 12/1/2013 (we subtract the zero-mean inside the loop by taking
    % into account only the overlapping area
    
    
    %% ECC, Forwards Additive Algorithm -------------------------------
    for i=1:noi
        
        disp(['Level: ' num2str(nol) ', Iteration: ' num2str(i)])
        %Image interpolation method
        str='linear'; % bilinear interpolation (you may also choose cubic)
        
        
        wim = spatial_interp(im, warp, str, transform, nx, ny); %inverse (backward) warping
        
        %ADDED 7/8/2012; MODIFIED 12/2/2013
        % define a mask to deal with warping outside the image borders 
        % (they may have negative values due to the subtraction of the mean value)
        ones_map = spatial_interp(ones(size(im)), warp, 'nearest', transform, nx, ny); %inverse (backward) warping
        numOfElem = sum(sum(ones_map~=0));
               
        meanOfWim = sum(sum(wim.*(ones_map~=0)))/numOfElem;
        meanOfTemp = sum(sum(temp.*(ones_map~=0)))/numOfElem;
        
        
        wim = wim-meanOfWim;% zero-mean image; is useful for brightness change compensation, otherwise you can comment this line
        tempzm = temp-meanOfTemp; % zero-mean template
        
        wim(ones_map==0) = 0; % for pixels outside the overlapping area
        tempzm(ones_map==0)=0;
        
        
        %Save current transform
        if (strcmp(transform,'affine')||strcmp(transform,'euclidean'))
            results(nol,i).warp = warp(1:2,:);
        else
            results(nol,i).warp = warp;
        end
        
        results(nol,i).rho = dot(temp(:),wim(:)) / norm(tempzm(:)) / norm(wim(:));
          
        if (i == noi) % the algorithm is executed (noi-1) times
            break;
        end
        
        % Gradient Image interpolation (warped gradients)
        wvx = spatial_interp(vx, warp, str, transform, nx, ny);
        wvy = spatial_interp(vy, warp, str, transform, nx, ny);
        
        % Compute the jacobian of warp transform
        J = warp_jacobian(nx, ny, warp, transform);
        
        % Compute the jacobian of warped image wrt parameters (matrix G in the paper)
        G = image_jacobian(wvx, wvy, J, nop);
        
        % Compute Hessian and its inverse
        C= G' * G;% C: Hessian matrix
        con=cond(C);
        if con>1.0e+15
            disp('->ECC Warning: Badly conditioned Hessian matrix. Check the initialization or the overlap of images.')
        end
        i_C = inv(C);
       
        % Compute projections of images into G
        Gt = G' * tempzm(:);
        Gw = G' * wim(:);
        
        
        %% ECC closed form solution
        
        % Compute lambda parameter
        num = (norm(wim(:))^2 - Gw' * i_C * Gw);
        den = (dot(tempzm(:),wim(:)) - Gt' * i_C * Gw);
        lambda = num / den;
        
        % Compute error vector
        imerror = lambda * tempzm - wim;
        
        % Compute the projection of error vector into Jacobian G
        Ge = G' * imerror(:);
        
        % Compute the optimum parameter correction vector
        delta_p = i_C * Ge;
        
        if (sum(isnan(delta_p)))>0 %Hessian is close to singular
            disp([' -> Algorithms stopped at ' num2str(i) '-th iteration of ' num2str(nol) '-th level due to bad condition of Hessian matrix.']);
            disp([' -> Current results have been saved at results(' num2str(nol) ',' num2str(i) ').warp and results(' num2str(nol) ',' num2str(i) ').rho.']);
            disp([' -> If you enabled a multilevel running, the output variables (warp, warpedImage) have been computed after mapping the current warp into the high-resolution level']);
            
            break_flag=1;
            break;
        end
        
        % Update parmaters
        warp = param_update(warp, delta_p, transform);
        
        
    end
    
    if break_flag==1
        break;
    end
    
    % modify the parameteres appropriately for next pyramid level
    if (nol>1)&(break_flag==0)
        warp = next_level(warp, transform,1);
    end
    
end

toc

if break_flag==1 % this conditional part is only executed when algorithm stops due to Hessian singularity
    for jj=1:nol-1
        warp = next_level(warp, transform,1);
        %m0=2*m0;
    end
    %margin=floor(m0*.05);
    %nx=margin+1:size(template,2)-margin;
    %ny=margin+1:size(template,1)-margin;
end
   
if break_flag == 1
    final_warp = warp;
else
    final_warp = results(1,end).warp;
end

% return the final warped image using the whole support area (include
% margins)

nx2 = 1:B;
ny2 = 1:A;

for ii = 1:sZi3
    warpedImage(:,:,ii) = spatial_interp(double(initImage(:,:,ii)), final_warp, str, transform, nx2, ny2);
end

warpedImage = uint8(warpedImage);
warp = final_warp;
%%%

Contact us