; file: selfaligncube_v2.pro ; init: Sep 22 2020 Rob Rutten Deil ; last: Oct 26 2020 Rob Rutten Deil ; note: selfalignfitscube.pro does the same using assoc ;+ pro selfaligncube_v2,incube,outcube,$ naverref=naverref,smear=smear,flatten=flatten,$ alignshift=alignshift,alignmetcalf=alignmetcalf,$ applypx2asym=applypx2asym,maxmetcalf=maxmetcalf,$ trimbox=trimbox,transformfile=transformfile,plotpath=plotpath,$ show=show ; self-align a cube in memory ; slow version Metcalf-iterating each time step per findalignimages.pro ; ; INPUTS: ; incube: image sequence [nx,ny,nt], any type ; ; OPTIONAL KEYWORD INPUTS: ; naverref: (tricky, much depending on data quality and duration) ; 0 or 1: successive pairwise alignment (default) ; > 1: align each to average last NN already aligned images before ; -1: align all to first image of sequence (bad at evolution) ; -2: align all to nt/2 mid-sequence image (bad at evolution) ; -3: align all to sequence time-average (smear likely useful) ; < -3: align each to temporal |naverref| box-car time-smoothed ; NB: runaway likely for naverref = 0 or 1, less for > 1 and < -3 ; smear: smooth by boxcar [px] smearing ; flatten: boxcar size in px to (x,y) flatten = subtract boxcar-smeared ; alignshift = 1/0: self-align by shifting only ; alignmetcalf = 1/0: self-align per Metcalf iteration (SLOW!) ; applypx2asym: as in findalignimages.pro but not to ref image (may be bad) ; maxmetcalf: max number Metcalf iterations (default 4) ; trimbox: subarea to self-align [xmin,ymin,xmax,ymax] (eg disk center) ; transformfile: filestring writecol results, default '' = none ; plotpath: path/filestart results graphs, default '' = none ; show = 1/0: showex every image pair ; ; OUTPUTS: ; - outcube: same size and type as incube ; - transformfile for application to parallel cubes ; - figures with shifts, (plus scales, rotate when Metcalf) ; ; REMARKS: ; possibly derotate first (for trimbox center), then rotate result ; for too large cubefiles use selfalignfitscube.pro (shift-align only) ; ; HISTORY: ; Oct 22 2020 RR: start ;- ; answer no-parameter query if (n_params(0) lt 2) then begin sp,selfaligncube_v2 return endif ; defaults for keywords if (n_elements(flatten) eq 0) then flatten=0 if (n_elements(smear) eq 0) then smear=0 if (n_elements(naverref) eq 0) then naverref=0 if (n_elements(alignshift) eq 0) then alignshift=0 if (n_elements(alignmetcalf) eq 0) then alignmetcalf=0 if (n_elements(applypx2asym) eq 0) then applypx2asym=0 if (n_elements(maxmetcalf) eq 0) then maxmetcalf=4 ;NB default if (n_elements(trimbox) eq 0) then trimbox=-1 if (n_elements(transformfile) eq 0) then transformfile='' if (n_elements(plotpath) eq 0) then plotpath='' if (n_elements(show) eq 0) then show=0 ; dimensions szcube=size(incube) nx=szcube[1] ny=szcube[2] nt=szcube[3] outcube=incube ; check inputs if (alignshift eq 0 and alignmetcalf eq 0) then begin print,' ##### selfaligncube_v2 abort: no shift nor metcalf align' return endif if (nx lt 1 or ny lt 1 or nt lt 1) then begin print,' ##### selfaligncube_v2 abort: input not dim (2,2,2) or larger' return endif ; first image initialization im_ref=incube[*,*,0] if (naverref eq -2) then im_ref=incube[*,*,nt/2] if (naverref eq -3) then im_ref=avg(incube,2) if (naverref lt -3) then $ smoothcube=smooth(incube,[1,1,abs(naverref)],/edge_truncate) ; define selfalignimages output arrays for writecol shiftx=fltarr(nt) shifty=fltarr(nt) scale1x=fltarr(nt)+1 scale1y=fltarr(nt)+1 rotate=fltarr(nt) ; ----- large loop over it for it=0,nt-1 do begin ; im1 = present image to get aligned to im_ref im1=incube[*,*,it] if (it eq 0 and naverref gt -2) then begin outcube[*,*,0]=im1 goto, SKIPFIRST endif ; naverref < -3 = use time-smoothed cube sample as reference if (naverref lt -3) then im_ref=smoothcube[*,*,it] ; shift-align only (no rotate, scale) if (alignshift eq 1) then begin shiftit=-findimshift_rr(im_ref,im1,/subpix,filter=10) if (show ne 0) then $ print,' ----- selfaligncube_v2 finds shifts it ='+trimd(it)+$ ' shift ='+trimd(shiftit) shiftx[it]=shiftit[0] shifty[it]=shiftit[1] reformimage,im1,im1_al,shift=-[shiftx,shifty] endif ; end shift-only mode ; ---------- Metcalf iteration per findalignimages.pro if (alignmetcalf eq 1) then begin ; make im2 from im_ref smaller than im1 for findalignimages.pro if (trimbox[0] eq -1) then $ reformimage,im_ref,im2,cutcentralx=nx-30,cutcentraly=ny-30 $ else im2=im_ref ; restrict missing to trimbox area if (it eq 0) then missing=avg(im2) ; fresh reset at every it for findalignimages.pro (no stickies) px1=1.0 px2=1.0 angle=0 px2asym=[1.,1.] shiftfor=[0,0] shiftrev=[0,0] inshift12=[0,0] ; run findalign: "forward" meaning im1 on im2 = current im_ref findalignimages,im1,im2,px1,px2,angle,px2asym,$ shiftfor,shiftrev,nxfor,nyfor,nxrev,nyrev,$ nxmuckfor=nxmuckfor,nymuckfor=nymuckfor,$ nxmuckrev=nxmuckrev,nymuckrev=nymuckrev,$ skipmetcalf=skipmetcalf,maxmetcalf=maxmetcalf,$ applypx2asym=applypx2asym,$ trimboxim2=trimbox,inshift12=inshift12,$ smearim1=smear,smearim2=smear,$ histopim1=histopim1,histopim2=histopim2,$ muckdarkim1=muckdarkim1,muckdarkim2=muckdarkim2,$ muckbrightim1=muckbrightim1,muckbrightim2=muckbrightim2,$ flatten=flatten,blink=0,show=0,verbose=0 ; ----- copy @@@@@@ demo stuff from findalignimages.pro ; but no use of nxfor and nyfor output but nx and ny instead and ; keep to im2 = reference pixels for all output images; pxasym to im1 ; copy @@@@@ demo = set scales but only incomplete copy, no scaling to finest pxratio=float(px1)/px2 ; refresh to newest Metcalfed px2, >1 for 1>2 pxscale1=[pxratio,pxratio] if (applypx2asym eq 1) then pxscale1=pxscale1/px2asym ; not to im2 here ; copy @@@@@ demo forward mode (but nxfor,nyfor > nx,ny; no px2asym to im2) reformimage,im1,im1for,$ congridfactor=pxscale1,shift=shiftfor,rotate=angle,$ missingvalue=missing,$ nxlarge=nx,nylarge=ny,cutcentralx=nx,cutcentraly=ny,splinip=1 ; finallly: Metcalf result outcube[*,*,it]=im1for ; store results for writecol shiftx[it]=shiftfor[0] shifty[it]=shiftfor[1] scale1x[it]=pxscale1[0] scale1y[it]=pxscale1[1] rotate[it]=angle ; print results for this timestep print,' ----- it ='+trimd(it)+'/'+trim(nt-1)+$ ': shiftx ='+trimd(shiftx[it],2)+' shifty ='+trimd(shifty[it],2)+$ ' scale1 ='+trimd([scale1x[it],scale1y[it]],4)+$ ' angle ='+trimd(angle,2) ; optionally inspect latest align pair if (show ne 0) then showex,im1for,im_ref,incube[*,*,it],/blink ;; STOP ; insert for checking variables here endif ; end Metcalf mode ; positive naverref options for next im_ref ; pairwise next one to this one if (naverref eq 0 or naverref eq 1) then im_ref=im1for ; next one to average last preceding ones if (naverref ge 2) then begin naver=min([it+1,naverref]) im_ref=total(outcube[*,*,(it-naver+1):it],3)/naver endif SKIPFIRST: endfor ; end large loop over it ; optionally write alignment results if (transformfile ne '') then begin if (alignshift eq 1) then writecol,transformfile,shiftx,shifty,fmt='(2F10.3)' if (alignmetcalf eq 1) then writecol,transformfile,$ shiftx,shifty,scale1x,scale1y,rotate,fmt='(7F10.4)' endif ; ; ----- optionally plot alignment results to show time variation if (plotpath ne '') then begin if (alignmetcalf eq 1) then begin ; Metcalf scales plot (only scale1) psfilename=plotpath+'scales.ps' axrat=1.62 ; golden ratio openpsplot,psfilename,thick=2,fontsize=9,xsize=8.8,ysize=8.8/axrat xtitle='time step' ytitle='scale factors' xrange=[0-0.2*nt,nt-1+0.2*nt] mima=minmax([scale1x,scale1y]) range=mima[1]-mima[0] yrange=[mima[0]-0.1*range,mima[1]+0.1*range] plot,indgen(nt),scale1x,$ position=[0.2,0.2,0.95,0.95],$ ; margins xticklen=0.03,yticklen=0.03/axrat,$ ; same-length ticks psym=1,symsize=1.3,xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 xyouts,-0.1*nt,scale1x[0],/data,charsize=1.3,'X' xyouts,nt-1+0.1*nt,scale1x[nt-1],/data,charsize=1.3,'X' oplot,indgen(nt),scale1y,psym=6,symsize=1.3 xyouts,-0.1*nt,scale1y[0],/data,charsize=1.3,'y' xyouts,nt-1+0.1*nt,scale1y[nt-1],/data,charsize=1.3,'y' closepsplot,psfilename,opengv=0 ; Metcalf rotate plot psfilename=plotpath+'rotate.ps' axrat=1.62 ; golden ratio openpsplot,psfilename,thick=2,fontsize=9,xsize=8.8,ysize=8.8/axrat xtitle='time step' ytitle='rotation angle [deg]' xrange=[0-0.2*nt,nt-1+0.2*nt] mima=minmax(rotate) range=mima[1]-mima[0] yrange=[mima[0]-0.1*range,mima[1]+0.1*range] plot,indgen(nt),rotate,$ position=[0.2,0.2,0.95,0.95],$ ; margins xticklen=0.03,yticklen=0.03/axrat,$ ; same-length ticks psym=1,symsize=1.3,xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 closepsplot,psfilename,opengv=0 endif ; shifts plot psfilename=plotpath+'shifts.ps' axrat=1.62 ; golden ratio openpsplot,psfilename,thick=2,fontsize=9,xsize=8.8,ysize=8.8/axrat xtitle='time step' ytitle='shifts [px]' xrange=[0-0.2*nt,nt-1+0.2*nt] mima=minmax([shiftx,shifty]) range=mima[1]-mima[0] yrange=[mima[0]-0.1*range,mima[1]+0.1*range] plot,indgen(nt),shiftx,$ position=[0.2,0.2,0.95,0.95],$ ; margins xticklen=0.03,yticklen=0.03/axrat,$ ; same-length ticks psym=1,symsize=1.3,xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 xyouts,-0.1*nt,shiftx[0],/data,charsize=1.3,'X' xyouts,nt-1+0.1*nt,shiftx[nt-1],/data,charsize=1.3,'X' oplot,indgen(nt),shifty,psym=6,symsize=1.3 xyouts,-0.1*nt,shifty[0],/data,charsize=1.3,'y' xyouts,nt-1+0.1*nt,shifty[nt-1],/data,charsize=1.3,'y' closepsplot,psfilename,opengv=0 endif ; end plot production end ; end of program ; =============== main for testing per IDLWAVE H-c ====================== ; same heavy-deform test as under findalignimages.pro, test that there ; sharp SST scene from web demo cd,'/home/rutten/rr/web/sdo-demo/2014-06-24/sst' infile='kr_width_6563_cut.fits' ; much Kevinesq detail incube=readfits(infile) ; int [932, 920, 53] nt=50 selectcube=incube[*,*,0:nt-1] ; set trimbox, avoid SST rotation triangles, smaller means faster trimbox=[200,500,400,700] ; nice 200x200 area upper left missing=avg(selectcube[200:400,500:700,*]) ; randomly deform input cube ("put SST on boat in Stockholm") ; (NB: program aborts on crop problem when deformed too much) szcube=size(selectcube) nt=szcube[3] r1=randomn(seed,nt) r2=randomn(seed,nt) ; continues sampling the same series via same seed r3=randomn(seed,nt) r4=randomn(seed,nt) r5=randomn(seed,nt) shiftsall=2.*[[r1],[r2]] angleall=0.5*r3 scaleall=[[1.0+0.005*r4],[1.0+0.005*r4]] ; random, no px2asym scaleall=[[1.0+0.002*r4],[1.0+0.002*r4]] ; random, no px2asym, small smearall=5.*abs((r5+1)>1)-4 ; add bad seeing moments above 1 px smearall=2.*abs((r5+1)>1)-1 ; add bad seeing moments above 1 px, better ;; ; equate all to first for precise showex alignment check ;; ; commented out uses original time evolution of the SST sequence + deform ;; for it=0,nt-1 do selectcube[*,*,it]=selectcube[*,*,0] ; initial reference image no deformation for aligning to it=0 as precise check shiftsall[0,0:1]=0 angleall[0]=0 scaleall[0,0:1]=[1.0,1.0] ; deform the selected cube reformcube,selectcube,deformcube,shift=shiftsall,rotate=angleall,$ scale=scaleall,missingvalue=missing,splinip=1 ; smear no array option in reformcube so smear per image except first for it=1,nt-1 do begin reformimage,deformcube[*,*,it],imout,smear=smearall[it] deformcube[*,*,it]=imout ; can't rewrite in place as I thought endfor ; parameters ;; naverref=-1 ; align all to first image ; BAD ;; naverref=0 ; align successively pairwise naverref=5 ; align each to average last NN already aligned ; OK ;; naverref=-2 ; align each to mid-time image ; BAD ;; naverref=-3 ; align each to temporal sequence mean (add smear) ; soso ;; naverref=-9 ; align each to time average of NN < -3 inputs around it smear=5 ; useful? flatten=50 ; useful? maxmetcalf=4 applypx2asym=0 show=0 alignmetcalf=1 ; output files transformfile='/tmp/selfalign.dat' plotpath='/tmp/selfalign_' ; run program selfaligncube_v2,deformcube,outcube,$ naverref=naverref,smear=smear,flatten=flatten,$ alignshift=alignshift,alignmetcalf=alignmetcalf,$ applypx2asym=applypx2asym,maxmetcalf=maxmetcalf,$ trimbox=trimbox,transformfile=transformfile,plotpath=plotpath,$ show=show ; inspect result (use showex first play button A versus B) ; zoom to center, check the granules under the fast (white) RBEs showex,outcube,deformcube,selectcube end