Non-rigid image alignment in twenty lines of Matlab

We have previously discussed rigid and non-rigid image alignment algorithms for calcium imaging data. If you have the image processing toolbox, a particularly compact solution for non-rigid image registration can be written in ~20 lines of code or so.  The code below returns a displacement field for each frame in the image sequence (output variable ‘disp’) that needs top be used to align the image.  The mean image after alignment is returned as well (output variable ‘m’).  The strategy follows the same recursive approach we used previously for rigid alignment.  The main Matlab function doing the calculation is imregdemons.


function [m,disp] = sbxalign_nonrigid(fname,idx)

if(length(idx)==1)
   A = squeeze(sbxread(fname,idx(1),1)); % just one frame... easy!
   m = A;
   disp = {zeros([size(A) 2])};
elseif (length(idx)==2) % align two frames
   A = squeeze(sbxread(fname,idx(1),1)); % read the frames
   B = squeeze(sbxread(fname,idx(2),1));
   [D,Ar] = imregdemons(A,B,[32 16 8 4],'AccumulatedFieldSmoothing',2.5,'PyramidLevels',4,'DisplayWaitBar',false);
   m = (Ar/2+B/2);
   disp = {D zeros([size(A) 2])};
else
   idx0 = idx(1:floor(end/2)); % split dataset in two
   idx1 = idx(floor(end/2)+1 : end); % recursive alignment
   [A,D0] = sbxalign_nonrigid(fname,idx0);
   [B,D1] = sbxalign_nonrigid(fname,idx1);
   [D,Ar] = imregdemons(A,B,[32 16 8 4],'AccumulatedFieldSmoothing',2.5,'PyramidLevels',4,'DisplayWaitBar',false);
   m = (Ar/2+B/2);
   D0 = cellfun(@(x) (x+D),D0,'UniformOutput',false); % concatenate distortions
   disp = [D0 D1];
end

Once the displacement field is computed one can obtain the aligned images by applying the distortion to each frame individually using imwarp.

The example below shows the original (top) and stabilized sequences in one case where imaging at high magnification clearly leads to non-rigid deformations that need correction.