; file: findimshift_rr.pro = find shift between images by cross-correlation ; init: 1992 by Pit Suetterlin as shc.pro ; torr: Aug 23 2016 Rob Rutten = copy of Pit's shc.pro from dotlib ; last: Oct 21 2020 Rob Rutten Deil ;+ FUNCTION findimshift_rr,im1,im2,filter=filter,subpix=subpix,$ trimbox=trimbox,greyborders=greyborders,cropborders=cropborders,$ croptriangles=croptriangles,niter=niter,cutiter=cutiter,$ fitpeak=fitpeak,power2=power2,polfit=polfit,range=range,$ xmaxcc=xmaxcc,ymaxcc=ymaxcc,maxcc=maxcc,showcc=showcc,$ fwhmfit=fwhmfit,itererror=itererror,finalshift=finalshift,$ verbose=verbose ; find the im2 shift needed to put im2 on im1 ; ; INPUTS: ; im1, im2: 2D arrays, same dimensions but may differ in type ; ; OPTIONAL KEYWORD INPUTS: ; filter=filter before Fourier transform ; 0 = no tapering ; 1 = full-length exponential taper (Pit's original) ; integer = cosine window edge taper, value = percentage of length ; I often set filter=10 to ignore edge spill-over from splinip=1 ; subpix: 0/1 interpolate maximum cross correlation function to sub-pixel ; trimbox: use subfield defined by trimmox=[xmin,ymin,xmax,ymax] ; greyborders = 1/0: make flat borders grey, also within iteration ; cropborders = 1/0: autocrop borders im2, also within iteration ; croptriangles = 1/0: autocrop borders and triangles im2 before iteration ; niter = number of iterations (default 10 with break at cutiter px) ; 0 = only one pass, no iteration ; cutiter = break value (default 0.03) ; fitpeak=fitpeak: fit cc peak by 1=Gaussian, 2=Lorentzian, 3=Moffat ; see explanation in mpfit2dpeak.pro ; power2: 0/1 use subarea of 2^n1+2^power-of-2 for faster FFT ; polfit=polfit: if not 0 first subtract polynomial fit with this degree ; range=range: constrain solution to be lower than this (scalar) ; xmaxcc=xmaxcc: output x location of cc max pixel (maybe subpixel) ; ymaxcc=ymaxcc: output y location of cc max pixel (maybe subpixel) ; maxcc=mmaxcc: output max of the cc matrix at pixel resolution ; showcc: 0/1 = show image and plot of cross-correlation ; fwhmfit: 0/1 = badness estimate = 2-D FWHM Gaussian fit to cc peak ; verbose: 0/1 = less/more output ; ; OUTPUTS: ; function result = [shiftx,shifty] for im2 relative to im1 ; xmaxcc = output x location of cc max pixel (maybe subpixel) ; ymaxcc = output y location of cc max pixel (maybe subpixel) ; maxcc = output max of the cc matrix at pixel resolution ; fwhmfit = estimate of cc peak widths = [fwhmx,fwhmy] ; finalshift = last iteration step ; itererror = 0/other: iteration didn't converge to cutiter value; return ; value = maximum of the last shift values ; ; PROCEDURE: ; compute cross-correlation per FFT and locate maximum ; ; MODIFICATION HISTORY: ; 13-Aug-1992 Peter (Pit) Suetterlin (PS) at KIS (later DOT, now SST) ; 29-Aug-1995 PS: Added Edge filter ; 30-Aug-1995 PS: Add Subpix interpolation. Slight rewrite of ; normal maximum finding (Use shift) ; 11-Sep-1995 PS: Forgot the 1-d case. Re-implemented. ; 10-Sep-1996 PS: Auto-Crop datasets to same dimensions. ; Also trim uneven dimensions, it's ugly for FFT ; 08-Jan-2001 PS: Keyword N2 to use good powers of 2 ; 22-Nov-2001 PS: Subtract of polynomial surface. Some formatting. ; Add on_error return. ; 17-Jul-2003 PS: When N2 is used, center the used part in common FOV ; Aug 23 2016 RR: rewrote to grasp; added edge cosine window ; Feb 4 2017 RR: () > []; mpfit2dpeak options ; Mar 19 2017 RR: niter, cutiter, verbose ; Jan 13 2018 RR: options grey or crop borders, crop triangles ; Oct 20 2020 RR: trimbox ;- ; keyword defaults if (n_elements(filter) eq 0) then filter=0 if (n_elements(subpix) eq 0) then subpix=0 if (n_elements(trimbox) eq 0) then trimbox=-1 if (n_elements(greyborders) eq 0) then greyborders=0 if (n_elements(cropborders) eq 0) then cropborders=0 if (n_elements(croptriangles) eq 0) then croptriangles=0 if (n_elements(niter) eq 0) then niter=10 if (n_elements(cutiter) eq 0) then cutiter=0.03 if (n_elements(power2) eq 0) then power2=0 if (n_elements(polfit) eq 0) then polfit=0 if (n_elements(fitpeak) eq 0) then fitpeak=0 if (n_elements(range) eq 0) then range=0 if (n_elements(fwhmfit) eq 0) then fwhmfit=0 if (n_elements(showcc) eq 0) then showcc=0 if (n_elements(verbose) eq 0) then verbose=0 ; initialize itererror=0 ; checks im1size=size(im1) im2size=size(im2) if (im1size[0] ne 2 or im2size[0] ne 2) then begin print,' ##### findimshift_rr abort: image 1 or image 2 not 2D array' return,-1 endif ;; if (im2size[1] ne im1size[1] or im2size[2] ne im1size[1]) then begin ;; print,' ===== findimshift_rr warning: unequal image sizes' ;; endif ; startup im1crop=im1 ; keep im1 and im2 intact although not used anymore here im2crop=im2 thiscrop=-1 ; optional trimbox if (trimbox[0] ne -1) then begin reformimage,im1crop,im1crop,trimbox=trimbox reformimage,im2crop,im2crop,trimbox=trimbox endif ; optionally make flat borders grey in both (better to filter into) if (greyborders ne 0) then begin reformimage,im1crop,im1crop,/greyborders reformimage,im2crop,im2crop,/greyborders endif ; optionally remove flat borders of im2 (implicit in croptriangles) if (cropborders ne 0 and croptriangles eq 0) then $ reformimage,im2crop,im2crop,/cropborders,cropbox=thiscrop ; alternatively optionally remove triangles (includes borders) of im2 if (croptriangles ne 0) then $ reformimage,im2crop,im2crop,/croptriangles,cropbox=thiscrop ; if whole image was found blank then quit returning shift zero if (thiscrop[0] eq -2) then return,[0,0] ; ###### start giant iteration loop defined by niter sumshift=[0,0] if (niter le 0) then nriter=1 else nriter=niter for iter=1,nriter do begin ; check im2crop size (maybe cropped above and also below) sizecrop=size(im2crop) if ((sizecrop[1] lt 0.3*im2size[1] or sizecrop[2] lt 0.3*im2size[2])$ and trimbox[0] eq -1) then $ print,' !!!!! findimshift_rr warning: much border or triangle cropping' ; crop im1 likewise to im2 if (thiscrop[0] ne -1) then $ im1crop=im1crop[thiscrop[0]:thiscrop[2],thiscrop[1]:thiscrop[3]] ;; ; check ;; print,' ----- nx2 ='+trimd(nx2)+' ny2 ='+trimd(ny2) ;; print,' ----- iter ='+trimd(iter) ;; print,' ----- thiscrop ='+trimd(thiscrop) ;; print,' ----- xsize = '+trimd(sizecrop[1])+' ysize = '+trimd(sizecrop[2]) ;; print,' ----- help,im1crop,im2crop' ;; help,im1crop,im2crop ;; flickrr,im1crop,im2crop,duration=4 ;; if (iter eq 4) then STOP ; average subtraction, delete empty dimensions p1=reform(im1crop-avg(im1crop)) p2=reform(im2crop-avg(im2crop)) ; ------ now a (modified) stretch from Pit Suetterlin's src.pro ;PS find common size (PS = Pit Suetterlin, i.e., following src.pro) s1=size(p1) s2=size(p2) ;PS take the smallest, I am surprized this works but it does ;RR but by now they should be equal sx=s1[1] 0:((sx/2+range) < sx)-1, $ (sy/2-range) > 0:((sy/2+range) < sy)-1) ; find max in integer px mx = max(cc, loc) ;RR loc = 1D subscript of max element (integer) maxcc=mx ;RR output value ;PS simple maximum location ccsz = size (cc) xmax = loc mod ccsz[1] ymax = loc/ccsz[1] ;PS if requested: linear sub-pixel interpolation of maximum location ;RR uses only neighbouring pixels if (subpix ne 0) then begin IF (xmax*ymax gt 0) AND (xmax lt (ccsz[1]-1)) $ AND (ymax lt (ccsz[2]-1)) THEN BEGIN ;PS Sep 91 phw try including more points in interpolations denom = mx*2 - cc[xmax-1, ymax] - cc[xmax+1, ymax] xfra = (xmax-.5) + (mx-cc[xmax-1, ymax])/denom denom = mx*2 - cc[xmax, ymax-1] - cc[xmax, ymax+1] yfra = (ymax-.5) + (mx-cc[xmax, ymax-1])/denom xmax = xfra ymax = yfra ENDIF endif ;RR if requested: cc peak fit if (fitpeak ne 0) then begin if (fitpeak eq 1) then yfit=mpfit2dpeak(cc,afit,/gaussian) if (fitpeak eq 2) then yfit=mpfit2dpeak(cc,afit,/lorentzian) if (fitpeak eq 3) then yfit=mpfit2dpeak(cc,afit,/moffat) xmax=afit[4] ymax=afit[5] endif ;RR determine fwhmfit = Gaussian-fit halfwidth if (iter eq 0 and fwhmfit ne 0) then begin yfit=mpfit2dpeak(cc,afit,/gaussian) fwhmfit=[afit[2],afit[3]] endif ;RR cc image + contours and max overlay if (showcc ne 0) then begin cgimage,bytscl(cc),/keep_aspect,charsize=1.5,$ /display,position=[0.1,0.1,0.95,0.95],$ /axes,axkeywords={font:1,ticklen:-0.02} cgcontour,cc,/onimage,nlevels=20 cgplots,[xmax,xmax],[ymax-1,ymax+1] cgplots,[xmax-1,xmax+1],[ymax,ymax] ; cc cross-section profiles window,/free plot,indgen(ccsz[1]),cc[*,ymax]/max(cc),yrange=[0.1,1.1],ystyle=1,psym=1 plots,[xmax,xmax],[0.1,1.1] oplot,indgen(ccsz[2]),cc[xmax,*]/max(cc),psym=2 plots,[ymax,ymax],[0.1,1.1],linestyle=2 endif ; result of his iteration step itershift=[xmax-ccsz[1]/2, ymax-ccsz[2]/2] sumshift=sumshift+itershift xmaxcc=xmax ymaxcc=ymax ; break when below threshold finalshift=itershift if (max(abs(itershift)) lt cutiter) then break ; shift im2crop to new version im2crop=shift_img_rr(im2crop,itershift,splinip=1) ; optionally make new borders resulting from the shifting grey if (max(abs(itershift)) ge 1 and greyborders ne 0) then $ reformimage,im2crop,im2crop,/greyborders ; optionally remove new borders resulting from the shifting if (max(abs(itershift)) ge 1 and cropborders ne 0) then $ reformimage,im2crop,im2crop,/cropborders,cropbox=thiscrop $ else thiscrop=-1 ; print if verbose if (verbose ne 0) then print,' ----- iteration'+trimd(iter)+$ ': shifts ='+trimd(itershift) endfor ; end of large iteration loop for successive improvement ; no convergence if (nriter gt 5 and iter ge nriter and verbose ne 0) then begin print,' ===== findimshift_rr: iteration not converged to cutiter' itererror=1 endif return,sumshift+finalshift END ;RR =============== add test per IDLWAVE hyper-C ======================== ; this test: small 304, 171, 160, 1700, magnetogram, continuum, Dopplergram cubedir='/home/rutten/data/SDO/2014-06-14-small/target/cubes/' infile=cubedir+'hmicont.fits' incube=readfits(infile) inim=incube[*,*,0] ; [102,102] ;; reformimage,inim,inim,rotate=10;; ,/splinip ;; reformimage,inim,inim,/croptriangles ;RR im1 is supposed not to have edges reformimage,inim,muckim,shift=-[3.333,4.6666],/splinip ;; ; Gregal limb data with eclipsed SDO (must fail) ;; cd,'/home/rutten/data/SST/tests/' ;; inim=readfits('gregal-flatsst.fits') ;; muckim=readfits('gregal-flatsdo.fits') ;; ; check ;; showex,inim,muckim ;; STOP trimbox=-1 trimbox=[30,30,90,90] ; worse ; get shift shift=findimshift_rr(inim,muckim,/subpix,filter=10,$ ; fitpeak=1,$ ; /power2,$ ; niter=100,$ cutiter=0.0001,$ trimbox=trimbox,$ greyborders=0,$ cropborders=0,$ croptriangles=0,$ itererror=itererror,$ finalshift=finalshift,$ xmaxcc=xmaxcc,ymaxcc=ymaxcc,maxcc=maxcc,$ verbose=1) ;; print,' ----- finalshift ='+trimd(finalshift) ;; print,'------ xmax, ymax, maxcc ='+trimd([xmaxcc,ymaxcc,maxcc],1) print,' ===== in: shift = -[3.333, 4.6666]' print,' ===== RR: shift = '+trimd(shift) ; check result while !d.window ne -1 do wdelete,!d.window im2new=shift_img_rr(muckim,shift,/splinip) ;RR splinip necessary for subpx ;; showex,inim,im2new,/blink ; compare with original Roberto Molowny and Pit Suetterlin programs shift=align(inim,muckim) print,' ===== RM: shift = '+trimd(shift) shift=shc(inim,muckim,/interpolate,/filter) print,' ===== PS: shift = '+trimd(shift) end