; file: cross_cor.pro ; init: Oct 14 2013 RR: from George H. Fischer, found at ; http://solarmuri.ssl.berkeley.edu/~fisher/public/software/FLCT/ pro cross_cor,f,g,ccor,shiftx,shifty,quiet=quiet,hires=hires ; ; - compute cross-correlation between 2 functions f & g and find the shift ; - in x and y which maximizes the cross-correlation function ; ; - INPUT: f, g are 2-d floating point arrays with x and y dimensions that ; - are the same ; ; OUTPUT: ccor,shiftx,shifty - ccor is the cross-correlation ; function, and shiftx, shifty are the shifts in f that are necessary ; to maximize the cross-correlation function. Note that ccor is ; shifted to be centered at nx/2, ny/2 to avoid aliasing problems when ; doing interpolation near the maximum ; nparms=N_Params() if (nparms ne 5) then begin Message, 'Usage: cross_cor,f,g,ccor,shiftx,shifty,/quiet,/hires', /info Message, 'Purpose: Compute cross-correlation function between', /info Message, ' images f,g, and find the position of the maximum', /info Message, ' to 0.1 (or .02 w/hires kwd) pixel resolution.', /info Message, 'Input: images f,g of equal size. ', /info Message, 'Output: ccor - (complex) 2-d cross-correlation function', /info Message, ' ccor has been shifted so that max occurs', /info Message, ' near the center of the array', /info Message, ' shiftx - the shift in x pixels of f needed to', /info Message, ' best correlate with g', /info Message, ' shifty - the shift in y pixels of f needed to', /info Message, ' best correlate with g', /info Message, '03/30/05 search for max confined to +/- 1gp instead of 5', /info Message, '01/15/04 hires kwd added for shifts to .02 pix res (GHF)', /info Message, '01/15/04 quiet kwd added to suppress warning msg', /info Message, 'Written: George H. Fisher UCB/SSL 11/7/2002 (version 0.1)', /info Message, 'Require: idl version 5.4 or more recent', /info return endif arrinfo=size(f) nx=arrinfo(1) ny=arrinfo(2) ; - - SHOULD check to make sure g and f have same size (but not done now) finv=fft(f,-1) ginv=fft(g,-1) spec=ginv*conj(finv) ; ccor=fft(spec,1) ; ; - - shift the origin to the middle of the box to avoid aliasing problems ; - - for small shifts. ; ccor=shift(ccor,nx/2,ny/2) absccor=abs(ccor) ; ; - - get coordinates of maximum cross-correlation to 1 pixel accuracy: ; cmax=max(absccor,indmax) ix = indmax mod nx iy = indmax/nx ; ; - - Now define a 10x10 pixel square region near the maximum ; around which to do higher order interpolation, to locate the ; - - maximum to 0.1 pixel resolution: ; shiftx0=float(ix) shifty0=float(iy) if keyword_set(hires) then begin nxinterp=101 & nyinterp = 101 & nfgppergp = 50 endif else begin nxinterp=21 & nyinterp = 21 & nfgppergp = 10 endelse range=float(nxinterp-1)/float(nfgppergp) xinterp=findgen(nxinterp)/float(nfgppergp)-0.5*range+shiftx0 yinterp=findgen(nyinterp)/float(nfgppergp)-0.5*range+shifty0 ; - - warning message if < 10 gp in some direction (but it keeps running) if keyword_set(quiet) then begin endif else begin if (min(xinterp) lt 0.) or (min(yinterp) lt 0.) or $ (max(xinterp) gt float(nx-1)) or (max(yinterp) gt float(ny-1)) $ then print,'WARNING: < 10 gp in at least 1 dimension in cross_cor or huge shift' endelse ; ; - - the array peakarea is interpolated abs. value of cross-correlation ; - - function ; peakarea=interpolate(absccor,xinterp,yinterp,cubic=-0.5,/grid,missing=0.) cmaxmax=max(peakarea,indmaxmax) ixx=indmaxmax mod nyinterp iyy=indmaxmax / nyinterp ; ; return the location of the maximum to subpixel resolution ; shiftxx = xinterp(ixx) shiftyy = yinterp(iyy) ; ; - - now shift result back to original coordinate system ; shiftx=shiftxx-float(nx/2) shifty=shiftyy-float(ny/2) return end