dftregistration

 dftregistration : Efficient subpixel image registration by crosscorrelation. 
 This code gives the same precision as the FFT upsampled cross correlation in a
 small fraction of the computation time and with reduced memory 
 requirements. It obtains an initial estimate of the crosscorrelation peak
 by an FFT and then refines the shift estimation by upsampling the DFT
 only in a small neighborhood of that estimate by means of a 
 matrix-multiply DFT. With this procedure all the image points are used to
 compute the upsampled crosscorrelation.
 Manuel Guizar - Dec 13, 2007

0001 function [output, Greg] = dftregistration(buf1ft,buf2ft,usfac)
0002 % dftregistration : Efficient subpixel image registration by crosscorrelation.
0003 % This code gives the same precision as the FFT upsampled cross correlation in a
0004 % small fraction of the computation time and with reduced memory
0005 % requirements. It obtains an initial estimate of the crosscorrelation peak
0006 % by an FFT and then refines the shift estimation by upsampling the DFT
0007 % only in a small neighborhood of that estimate by means of a
0008 % matrix-multiply DFT. With this procedure all the image points are used to
0009 % compute the upsampled crosscorrelation.
0010 % Manuel Guizar - Dec 13, 2007
0011 
0012 % Portions of this code were taken from code written by Ann M. Kowalczyk
0013 % and James R. Fienup.
0014 % J.R. Fienup and A.M. Kowalczyk, "Phase retrieval for a complex-valued
0015 % object by using a low-resolution image," J. Opt. Soc. Am. A 7, 450-458
0016 % (1990).
0017 
0018 % Citation for this algorithm:
0019 % Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
0020 % "Efficient subpixel image registration algorithms," Opt. Lett. 33,
0021 % 156-158 (2008).
0022 
0023 % Inputs
0024 % buf1ft    Fourier transform of reference image,
0025 %           DC in (1,1)   [DO NOT FFTSHIFT]
0026 % buf2ft    Fourier transform of image to register,
0027 %           DC in (1,1) [DO NOT FFTSHIFT]
0028 % usfac     Upsampling factor (integer). Images will be registered to
0029 %           within 1/usfac of a pixel. For example usfac = 20 means the
0030 %           images will be registered within 1/20 of a pixel. (default = 1)
0031 
0032 % Outputs
0033 % output =  [error,diffphase,net_row_shift,net_col_shift]
0034 % error     Translation invariant normalized RMS error between f and g
0035 % diffphase     Global phase difference between the two images (should be
0036 %               zero if images are non-negative).
0037 % net_row_shift net_col_shift   Pixel shifts between images
0038 % Greg      (Optional) Fourier transform of registered version of buf2ft,
0039 %           the global phase difference is compensated for.
0040 
0041 % Default usfac to 1
0042 if exist('usfac')~=1, usfac=1; end
0043 
0044 % Compute error for no pixel shift
0045 if usfac == 0,
0046     CCmax = sum(sum(buf1ft.*conj(buf2ft))); 
0047     rfzero = sum(abs(buf1ft(:)).^2);
0048     rgzero = sum(abs(buf2ft(:)).^2); 
0049     error = 1.0 - CCmax.*conj(CCmax)/(rgzero*rfzero); 
0050     error = sqrt(abs(error));
0051     diffphase=atan2(imag(CCmax),real(CCmax)); 
0052     output=[error,diffphase];
0053         
0054 % Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the
0055 % peak
0056 elseif usfac == 1,
0057     [m,n]=size(buf1ft);
0058     CC = ifft2(buf1ft.*conj(buf2ft));
0059     [max1,loc1] = max(CC);
0060     [max2,loc2] = max(max1);
0061     rloc=loc1(loc2);
0062     cloc=loc2;
0063     CCmax=CC(rloc,cloc); 
0064     rfzero = sum(abs(buf1ft(:)).^2)/(m*n);
0065     rgzero = sum(abs(buf2ft(:)).^2)/(m*n); 
0066     error = 1.0 - CCmax.*conj(CCmax)/(rgzero(1,1)*rfzero(1,1));
0067     error = sqrt(abs(error));
0068     diffphase=atan2(imag(CCmax),real(CCmax)); 
0069     md2 = fix(m/2); 
0070     nd2 = fix(n/2);
0071     if rloc > md2
0072         row_shift = rloc - m - 1;
0073     else
0074         row_shift = rloc - 1;
0075     end
0076 
0077     if cloc > nd2
0078         col_shift = cloc - n - 1;
0079     else
0080         col_shift = cloc - 1;
0081     end
0082     output=[error,diffphase,row_shift,col_shift];
0083     
0084 % Partial-pixel shift
0085 else
0086     
0087     % First upsample by a factor of 2 to obtain initial estimate
0088     % Embed Fourier data in a 2x larger array
0089     [m,n]=size(buf1ft);
0090     mlarge=m*2;
0091     nlarge=n*2;
0092     CC=zeros(mlarge,nlarge);
0093     CC(m+1-fix(m/2):m+1+fix((m-1)/2),n+1-fix(n/2):n+1+fix((n-1)/2)) = ...
0094         fftshift(buf1ft).*conj(fftshift(buf2ft));
0095   
0096     % Compute crosscorrelation and locate the peak
0097     CC = ifft2(ifftshift(CC)); % Calculate cross-correlation
0098     [max1,loc1] = max(CC);
0099     [max2,loc2] = max(max1);
0100     rloc=loc1(loc2);cloc=loc2;
0101     CCmax=CC(rloc,cloc);
0102     
0103     % Obtain shift in original pixel grid from the position of the
0104     % crosscorrelation peak
0105     [m,n] = size(CC); md2 = fix(m/2); nd2 = fix(n/2);
0106     if rloc > md2 
0107         row_shift = rloc - m - 1;
0108     else
0109         row_shift = rloc - 1;
0110     end
0111     if cloc > nd2
0112         col_shift = cloc - n - 1;
0113     else
0114         col_shift = cloc - 1;
0115     end
0116     row_shift=row_shift/2;
0117     col_shift=col_shift/2;
0118 
0119     % If upsampling > 2, then refine estimate with matrix multiply DFT
0120     if usfac > 2,
0121         %%% DFT computation %%%
0122         % Initial shift estimate in upsampled grid
0123         row_shift = round(row_shift*usfac)/usfac; 
0124         col_shift = round(col_shift*usfac)/usfac;     
0125         dftshift = fix(ceil(usfac*1.5)/2); %% Center of output array at dftshift+1
0126         % Matrix multiply DFT around the current shift estimate
0127         CC = conj(dftups(buf2ft.*conj(buf1ft),ceil(usfac*1.5),ceil(usfac*1.5),usfac,...
0128             dftshift-row_shift*usfac,dftshift-col_shift*usfac))/(md2*nd2*usfac^2);
0129         % Locate maximum and map back to original pixel grid
0130         [max1,loc1] = max(CC);   
0131         [max2,loc2] = max(max1); 
0132         rloc = loc1(loc2); cloc = loc2;
0133         CCmax = CC(rloc,cloc);
0134         rg00 = dftups(buf1ft.*conj(buf1ft),1,1,usfac)/(md2*nd2*usfac^2);
0135         rf00 = dftups(buf2ft.*conj(buf2ft),1,1,usfac)/(md2*nd2*usfac^2);  
0136         rloc = rloc - dftshift - 1;
0137         cloc = cloc - dftshift - 1;
0138         row_shift = row_shift + rloc/usfac;
0139         col_shift = col_shift + cloc/usfac;    
0140 
0141     % If upsampling = 2, no additional pixel shift refinement
0142     else    
0143         rg00 = sum(sum( buf1ft.*conj(buf1ft) ))/m/n;
0144         rf00 = sum(sum( buf2ft.*conj(buf2ft) ))/m/n;
0145     end
0146     error = 1.0 - CCmax.*conj(CCmax)/(rg00*rf00);
0147     error = sqrt(abs(error));
0148     diffphase=atan2(imag(CCmax),real(CCmax));
0149     % If its only one row or column the shift along that dimension has no
0150     % effect. We set to zero.
0151     if md2 == 1,
0152         row_shift = 0;
0153     end
0154     if nd2 == 1,
0155         col_shift = 0;
0156     end
0157     output=[error,diffphase,row_shift,col_shift];
0158 end  
0159 
0160 % Compute registered version of buf2ft
0161 if (nargout > 1)&&(usfac > 0),
0162     [nr,nc]=size(buf2ft);
0163     Nr = ifftshift([-fix(nr/2):ceil(nr/2)-1]);
0164     Nc = ifftshift([-fix(nc/2):ceil(nc/2)-1]);
0165     [Nc,Nr] = meshgrid(Nc,Nr);
0166     Greg = buf2ft.*exp(i*2*pi*(-row_shift*Nr/nr-col_shift*Nc/nc));
0167     Greg = Greg*exp(i*diffphase);
0168 elseif (nargout > 1)&&(usfac == 0)
0169     Greg = buf2ft*exp(i*diffphase);
0170 end
0171 return
0172 
0173 function out=dftups(in,nor,noc,usfac,roff,coff)
0174 % function out=dftups(in,nor,noc,usfac,roff,coff);
0175 % Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just
0176 % a small region.
0177 % usfac         Upsampling factor (default usfac = 1)
0178 % [nor,noc]     Number of pixels in the output upsampled DFT, in
0179 %               units of upsampled pixels (default = size(in))
0180 % roff, coff    Row and column offsets, allow to shift the output array to
0181 %               a region of interest on the DFT (default = 0)
0182 % Recieves DC in upper left corner, image center must be in (1,1)
0183 % Manuel Guizar - Dec 13, 2007
0184 % Modified from dftus, by J.R. Fienup 7/31/06
0185 
0186 % This code is intended to provide the same result as if the following
0187 % operations were performed
0188 %   - Embed the array "in" in an array that is usfac times larger in each
0189 %     dimension. ifftshift to bring the center of the image to (1,1).
0190 %   - Take the FFT of the larger array
0191 %   - Extract an [nor, noc] region of the result. Starting with the
0192 %     [roff+1 coff+1] element.
0193 
0194 % It achieves this result by computing the DFT in the output array without
0195 % the need to zeropad. Much faster and memory efficient than the
0196 % zero-padded FFT approach if [nor noc] are much smaller than [nr*usfac nc*usfac]
0197 
0198 [nr,nc]=size(in);
0199 % Set defaults
0200 if exist('roff')~=1, roff=0; end
0201 if exist('coff')~=1, coff=0; end
0202 if exist('usfac')~=1, usfac=1; end
0203 if exist('noc')~=1, noc=nc; end
0204 if exist('nor')~=1, nor=nr; end
0205 % Compute kernels and obtain DFT by matrix products
0206 kernc=exp((-i*2*pi/(nc*usfac))*( ifftshift([0:nc-1]).' - floor(nc/2) )*( [0:noc-1] - coff ));
0207 kernr=exp((-i*2*pi/(nr*usfac))*( [0:nor-1].' - roff )*( ifftshift([0:nr-1]) - floor(nr/2)  ));
0208 out=kernr*in*kernc;
0209 return