; file: sdo_shiftaligncubefiles.pro = find shifts(t) and align two fits cubes ; init: Nov 11 2013 Rob Rutten Sac Peak (then crossalignfitscubes.pro) ; last: Oct 31 2020 Rob Rutten Deil ;+ pro sdo_shiftaligncubefiles,infile1,infile2,$ driftshifts,wobbleshifts,$ outfile1=outfile1,outfile2=outfile2,$ alignmean=alignmean,trange=trange,trimbox=trimbox,$ niter=niter,cutiter=cutiter,$ nxtile=nxtile,tilecut=tilecut,tileprecision=tileprecision,$ driftspline=driftspline,wobblespline=wobblespline,nspliter=nspliter,$ cutminmax=cutminmax,scale=scale,rotate=rotate,$ smearref=smearref,plotdir=plotdir,$ saveshifts=saveshifts,missing=missing,verbose=verbose ; find cross-align for two image-sequence fitscube files ; their timing must be identical ; INPUTS: ; infile1 = string with path/name first input file = reference file ; infile2 = string with path/name second input file = file to align ; OPTIONAL KEYWORD INPUTS: ; alignmean =1/0: align with mean shifts after spline etc (default 0) ; trange: 2-element vector with it-start and it-end (default [0,-1]) ; trimbox = 4-elem integer vector to define borders to ignore ; niter = nr iterations findimshift_rr.pro (with break at cutiter px) ; 0 = only one pass, no iteration ; cutiter = break value in call findimshift_rr.pro (default 0.01) ; nxtile = 0 no tiling ; > 0 = tile x-length (e.g., 50 for AIA = 30 arcsec) ; tilecut = cut tile outliers at tilecut*rms (default 2.) ; tileprecision = stop tile loop at this increment (default 0.05) ; driftspline: smooth value for cubic spline approximation to drifts ; = 0: none (default) ; = 1E8: advised ; = 1E-5: strict through all x-shift and y-shift points ; = 1E+10: straight line fit ; wobblespline: smooth value for wobble (default 0, e.g. driftspline*1.E-6) ; nspliter = iterations to remove outliers from spline fit (default 0) ; cutminmax = set min,max cutoff (fraction range 1st image) and reverse ; e.g. [0.0,0.5] puts weight on the darkest areas ; scale = scale factor for 2nd cube (before shift determination) ; rotate = rotate value for 2nd cube (before shift determination) ; smearref = smear refcube images temporarily (use e.g. 3.*px1/px2) ; missing = intensity for blank side excess (default -1 = mean 1st image) ; verbose =1/0: extra printout ; OUTPUTS: ; driftshifts = drift shift values (spline solutions for driftspline ne 0) ; wobbleshifts = wobble shift values (idem, modulation from driftshifts)) ; OPTIONAL KEYWORD OUTPUTS: ; outfile1 = string with path/name first result file (may be cut to common) ; outfile2 = string with path/name second result file (aligned to first) ; saveshifts = optional filename to write shifts (ASCII, default in /tmp) ; plotdir = dir for shifts plots (default none) ; RESTRICTIONS: ; the cubes must have identical timing sequences ; smearing useful only if refcube is finer-scale ; REMARKS: ; finds only shifts, no rotation or scale changes (use findalignimages.pro) ; uses assoc so can handle larger files than memory size ; the shifts are printed when verbose and plotted when plotdir ne '' ; writing result files can take long ; HISTORY: ; Nov 11 2013 RR: start (as crossalignfitscubes.pro) ; Nov 12 2013 RR: added iteration ; Nov 24 2013 RR: added alignmean ; Nov 25 2013 RR: added cutminmax,scale,rotate ; Jan 11 2014 RR: ps plotshifts instead of cgplot window ; Aug 18 2014 RR: output shifts ; Feb 11 2015 RR: option smearref ; Sep 28 2015 RR: smear the refcube, not 2nd cube ; Mar 22 2017 RR: revamp including tile, spline, outlier removal ; Dec 28 2019 RR: add wobbleshifts ; Jan 22 2020 RR: cross-aligned outfiles optional ; Oct 31 2020 RR: name shiftalign > sdo_shiftalign ;- ; answer a no-parameter query if (n_params() lt 4) then begin sp,sdo_shiftaligncubefiles return endif ; defaults for keywords if (n_elements(outfile1) eq 0) then outfile1='' if (n_elements(outfile2) eq 0) then outfile2='' if (n_elements(alignmean) eq 0) then alignmean=0 if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(niter) eq 0) then niter=5 if (n_elements(cutiter) eq 0) then cutiter=0.01 if (n_elements(nxtile) eq 0) then nxtile=0 if (n_elements(tilecut) eq 0) then tilecut=2. if (n_elements(tileprecision) eq 0) then tileprecision=0.05 if (n_elements(driftspline) eq 0) then driftspline=0 if (n_elements(wobblespline) eq 0) then wobblespline=0 if (n_elements(nspliter) eq 0) then nspliter=0 if (n_elements(trimbox) eq 0) then trimbox=[-1,-1,-1,-1] if (n_elements(cutminmax) eq 0) then cutminmax=0 if (n_elements(scale) eq 0) then scale=0 if (n_elements(rotate) eq 0) then rotate=0 if (n_elements(smearref) eq 0) then smearref=0 if (n_elements(plotdir) eq 0) then plotdir='' if (n_elements(saveshifts) eq 0) then saveshifts='/tmp/saveshifts' if (n_elements(missing) eq 0) then missing=-1 if (n_elements(verbose) eq 0) then verbose=0 ; print pacifier print,' ----- sdo_shiftaligncubefiles starts on '$ +file_basename(infile1)+' and '+file_basename(infile2) ; check parameters if (driftspline eq 0) then wobblespline=0 if (n_elements(trimbox) ne 4) then begin print,' ##### shiftalignfiles abort: wrong trimbox' return endif ; set endian bigendian=1 ; get cube1 dimensions and file datatype from fits header inheader=headfits_rr(infile1) inheadersize=(1+fix(n_elements(inheader)/36.))*2880 nx1=fxpar(inheader,'naxis1') ny1=fxpar(inheader,'naxis2') nt1=fxpar(inheader,'naxis3') bitpix=fxpar(inheader,'bitpix') ; open input file 1 for assoc get_lun, unit_in1 if (bigendian) then openr,unit_in1,infile1,/swap_if_little_endian $ else openr,unit_in1,infile1 if (bitpix eq -32) then inassoc1=assoc(unit_in1,fltarr(nx1,ny1),inheadersize) if (bitpix eq 16) then inassoc1=assoc(unit_in1,intarr(nx1,ny1),inheadersize) if (bitpix eq 8) then inassoc1=assoc(unit_in1,bytarr(nx1,ny1),inheadersize) ; get cube2 dimensions and file datatype from the fits header inheader=headfits_rr(infile2) inheadersize=(1+fix(n_elements(inheader)/36.))*2880 nx2=fxpar(inheader,'naxis1') ny2=fxpar(inheader,'naxis2') nt2=fxpar(inheader,'naxis3') bitpix2=fxpar(inheader,'bitpix') ; check they have the same bitpix if (bitpix2 ne bitpix) then begin print,' #### sdo_shiftaligncubefiles abort: cubes do not have same bitpix' return endif ; open input file 2 for assoc get_lun, unit_in2 if (bigendian) then openr,unit_in2,infile2,/swap_if_little_endian $ else openr,unit_in2,infile2 if (bitpix eq -32) then inassoc2=assoc(unit_in2,fltarr(nx2,ny2),inheadersize) if (bitpix eq 16) then inassoc2=assoc(unit_in2,intarr(nx2,ny2),inheadersize) if (bitpix eq 8) then inassoc2=assoc(unit_in2,bytarr(nx2,ny2),inheadersize) ; find common dimensions nx=min([nx1,nx2]) ny=min([ny1,ny2]) nt=min([nt1,nt2]) ; set itstart and itend itstart=trange[0] if (trange[1] eq -1) then itend=nt-1 else itend=trange[1] ;RR???? Aug 13 2022 nt1-1 => nt-1 ; define optional output header outheader=inheader sxaddpar,outheader,'NAXIS3',itend-itstart+1 sizeoutheader=size(outheader) ; fits outheader = Nx36 "card images" = Nx2880 bytes outheadersize=(1+fix(sizeoutheader[1]/36.))*2880 ; open optional output file 1 for assoc, write header if (outfile1 ne '') then begin get_lun, unit_out1 if (bigendian) then openw,unit_out1,outfile1,/swap_if_little_endian $ else openw,unit_out1,outfile1 if (bitpix eq -32) then outassoc1=assoc(unit_out1,fltarr(nx,ny),outheadersize) if (bitpix eq 16) then outassoc1=assoc(unit_out1,intarr(nx,ny),outheadersize) if (bitpix eq 8) then outassoc1=assoc(unit_out1,bytarr(nx,ny),outheadersize) if (outheadersize ne 0) then begin rec=assoc(unit_out1, bytarr(outheadersize)) rec[0]=byte(outheader) endif endif ; open optional output file 2 for assoc, write header if (outfile2 ne '') then begin get_lun, unit_out2 if (bigendian) then openw,unit_out2,outfile2,/swap_if_little_endian $ else openw,unit_out2,outfile2 if (bitpix eq -32) then outassoc2=assoc(unit_out2,fltarr(nx,ny),outheadersize) if (bitpix eq 16) then outassoc2=assoc(unit_out2,intarr(nx,ny),outheadersize) if (bitpix eq 8) then outassoc2=assoc(unit_out2,bytarr(nx,ny),outheadersize) if (outheadersize ne 0) then begin rec=assoc(unit_out2, bytarr(outheadersize)) rec[0]=byte(outheader) endif endif ; specify cutout size specified by trimbox if (trimbox[0] ne -1) then begin xmincut=trimbox[0] xmaxcut=trimbox[2] ymincut=trimbox[1] ymaxcut=trimbox[3] endif else begin xmincut=0 xmaxcut=nx-1 ymincut=0 ymaxcut=ny-1 endelse ; define intensity ranges if cutints requested if (cutminmax ne 0) then begin im1=inassoc1[0] ;; mm1=minmax(histo_opt_rr(im1)) mm1=minmax(im1) im2=inassoc2[0] ;; mm2=minmax(histo_opt_rr(im2)) mm2=minmax(im2) endif ; big loop = shift determination between successive image pairs ; ------------------------------------------------------------- shiftx=fltarr(itend-itstart+1) shifty=fltarr(itend-itstart+1) for it=itstart,itend do begin ; get the images (MB: this routine requires identical cube timings) im1=inassoc1[it] im2=inassoc2[it] ; set greyscale value missing for filling blank edges if (it eq itstart) then begin if (missing eq -1) then miss=avg(im2) else miss=0 endif ; scale and/or rotate the second image if (scale ne 0 or rotate ne 0) then $ im2=shift_img(im2,xscale=scale,yscale=scale,rot=rotate,missing=miss) ; get common area im1com=im1[0:nx-1,0:ny-1] im2com=im2[0:nx-1,0:ny-1] im1cut=im1com[xmincut:xmaxcut,ymincut:ymaxcut] im2cut=im2com[xmincut:xmaxcut,ymincut:ymaxcut] ; apply intensity cuts when specified and bytscale if (cutminmax ne 0) then begin im1byte=bytscl(im1cut,$ min=mm1[0]+cutminmax[0]*(mm1[1]-mm1[0]),$ max=mm1[0]+cutminmax[1]*(mm1[1]-mm1[0])) im2byte=bytscl(im2cut,$ min=mm2[0]+cutminmax[0]*(mm2[1]-mm2[0]),$ max=mm2[0]+cutminmax[1]*(mm2[1]-mm2[0])) endif else begin im1byte=bytscl(im1cut) im2byte=bytscl(im2cut) endelse ; smear im1byte temporarily if requested if (smearref ne 0) then reformimage,im1byte,im1smear,smear=fix(smearref) $ else im1smear=im1byte ; get shift between these images (for SST croptriangles may be important) if (nxtile eq 0) then begin ;; shift=align(im1smear,im2byte) ; Molowny Horas ;; shift=shc(im1smear,im2byte,/interpolate,/filter,/n2) ; Suetterlin ;; im2bytenew=shift_sub(im2byte,shift[0],shift[1]) ; Chae shift=findimshift_rr(im1smear,im2byte,/subpix,niter=niter,filter=10,$ cropborders=1,croptriangles=1,cutiter=cutiter) endif else begin shift=findtiledimshift(im1smear,im2byte,nxtile,tilecut,tileprecision) endelse shiftx[it-itstart]=shift[0] shifty[it-itstart]=shift[1] if (verbose) then print,' ==== time step '+trim(it)+$ ' shift (x,y) = '+trim(shift[0])+' '+trim(shift[1]) ; optional output first = only common area if (outfile1 ne '') then outassoc1[it-itstart]=im1com endfor ; end big loop over timesteps ; define initial shift parameters nshifts=n_elements(shiftx) itarr=indgen(nshifts) nshifts0=nshifts shiftx0=shiftx ; originally measured offsets every time step shifty0=shifty ; (outliers get cropped out in interation) itarr0=itarr ; ==================== spline determinations ; if requested get cubic slice approximations with optional outlier removal if (driftspline ne 0) then begin ; iteratively fit to handle outliers (tricky dicky) sdevx=1. sdevy=1. cutsdev=3. splinx0=fltarr(nt) spliny0=fltarr(nt) for spliter=0,nspliter-1 do begin ; first fine spline to cut out single outliers, then smoother for drift if (spliter eq 0) then begin drxsmooth=wobblespline drysmooth=wobblespline woxsmooth=wobblespline woysmooth=wobblespline endif if (spliter gt 0) then begin drxsmooth=driftspline drysmooth=driftspline endif ; adapt array size after earlier outlier deletions, with check > 1 left nshifts=n_elements(shiftx) if (nshifts lt 3) then begin print,' ##### shift_aligfiles OOPS: '+$ ' too few values left after outlier removals; all ignored' nshifts=nt shiftx=shiftx0 shifty=shifty0 itarr=itarr0 endif splinx=fltarr(nshifts) spliny=fltarr(nshifts) ; define driftspline approximation (modified Vitas program) splinapp_set,itarr,shiftx,ax,bx,cx,dx,smooth=drxsmooth splinapp_set,itarr,shifty,ay,by,cy,dy,smooth=drysmooth ; get solutions at the same it's for it=0,nshifts-1 do splinx[it]=$ splinapp_get(itarr[it],itarr,ax,bx,cx,dx) for it=0,nshifts-1 do spliny[it]=$ splinapp_get(itarr[it],itarr,ay,by,cy,dy) ; define wobble spline approximation as modulation from driftspline if (spliter gt 0) then begin splinapp_set,itarr,shiftx-splinx,woax,wobx,wocx,wodx,smooth=woxsmooth splinapp_set,itarr,shifty-spliny,woay,woby,wocy,wody,smooth=woysmooth endif ; get standard deviations diffx=abs(shiftx-splinx) momx=moment(diffx,sdev=sdevx) diffy=abs(shifty-spliny) momy=moment(diffy,sdev=sdevy) ; find and treat outliers for next pass if (spliter lt nspliter-1) then begin ; find outliers in x inddelete=-1 indexcessx=where(diffx gt cutsdev*sdevx) if (indexcessx[0] ne -1) then begin inddelete=indexcessx if (verbose) then print, $ ' X outliers at it values ',trim(indexcessx) endif ; find outliers in y indexcessy=where(diffy gt cutsdev*sdevy) if (indexcessy[0] ne -1) then begin inddelete=[indexcessx,indexcessy] if (verbose) then print, $ ' Y outliers at it values ',trim(indexcessy) endif ; remove new outliers if (inddelete[0] ne -1 and n_elements(inddelete) lt nt) $ then remove,inddelete,itarr,shiftx,shifty endif ; check if (verbose) then $ print,' ----- spliter,sdevx,sdevy,drxsmooth,drysmooth = '+$ trim(spliter)+' '+trim(sdevx)+' '+trim(sdevy)+$ ' '+trim(drxsmooth)+' '+trim(drysmooth) endfor ; end of spline iteration endif ; ====================== end driftspline determination ; get drift spline solution at all original it's splinx0=fltarr(nshifts0) spliny0=fltarr(nshifts0) if (driftspline ne 0) then begin for it=0,nshifts0-1 do splinx0[it]=$ splinapp_get(it,itarr,ax,bx,cx,dx) for it=0,nshifts0-1 do spliny0[it]=$ splinapp_get(it,itarr,ay,by,cy,dy) endif ; get wobble spline solution at all original it's wosplinx0=fltarr(nshifts0) wospliny0=fltarr(nshifts0) if (wobblespline ne 0) then begin for it=0,nshifts0-1 do wosplinx0[it]=$ splinapp_get(it,itarr,woax,wobx,wocx,wodx) for it=0,nshifts0-1 do wospliny0[it]=$ splinapp_get(it,itarr,woay,woby,wocy,wody) endif ; get mean shifts if (alignmean ne 0) then begin if (driftspline eq 0) then begin meanx=avg(shiftx0) meany=avg(shifty0) endif else begin meanx=avg(splinx0) meany=avg(spliny0) endelse if (verbose) then print,' ----- mean shift = '+trim(meanx)+' '+trim(meany) endif ; now write cropped and shifted images on optional outfile2 if (outfile2 ne '') then begin for it=0,nshifts0-1 do begin im2=inassoc2[it] im2com=im2[0:nx-1,0:ny-1] if (driftspline eq 0 and alignmean eq 0) then begin im2new=shift_img(im2com,[shiftx[it],shifty[it]],missing=miss,$ xscale=scale,yscale=scale,rot=rotate) endif if (driftspline ne 0 and alignmean eq 0) then begin im2new=shift_img(im2com,[splinx0[it]+wosplinx0[it],$ spliny0[it]+wospliny0[it]],$ missing=miss,xscale=scale,yscale=scale,rot=rotate) endif if (alignmean eq 1) then begin im2new=shift_img(im2com,[meanx,meany],missing=miss,$ xscale=scale,yscale=scale,rot=rotate) endif outassoc2[it]=im2new endfor endif ; free the input and optional output files free_lun,unit_in1,unit_in2 if (outfile1 ne '') then free_lun,unit_out1 if (outfile2 ne '') then free_lun,unit_out2 ; make shifts plot if (plotdir ne '') then begin filestr1=file_basename(infile1,'.fits') filestr2=file_basename(infile2,'.fits') psfilename=filestr1+'--'+filestr2+'.ps' psfile=plotdir+'/'+psfilename openpsplot,psfile,thick=2,fontsize=9,xsize=8.8,ysize=8.8 ; yrange set neglecting extreme samples or by driftspline extremes if splined if (driftspline eq 0) then begin meanx=avg(shiftx0) meany=avg(shifty0) dummy=moment(shiftx0,sdev=sdevx) dummy=moment(shifty0,sdev=sdevy) yrange=[min([meanx-3*sdevx,meany-3*sdevy]),$ max([meanx+3*sdevx,meany+3*sdevy])] endif else begin yrange=[min([min(splinx0+wosplinx0),min(spliny0+wospliny0)]),$ max([max(splinx0+wosplinx0),max(spliny0+wospliny0)])] yrange=[yrange[0]-0.2*(yrange[1]-yrange[0]),$ yrange[1]+0.2*(yrange[1]-yrange[0])] endelse plot,itarr0,shiftx0,$ position=[0.2,0.1,0.95,0.92],$ psym=2,symsize=4./sqrt(nshifts),$ xrange=[0,itend],yrange=yrange,xstyle=1,ystyle=1,$ xtitle='time step',ytitle='shift in px' oplot,itarr0,shifty0,$ psym=6,symsize=4./sqrt(nshifts) ; overplot driftspline solution if there is one if (driftspline ne 0) then begin oplot,itarr0,splinx0 oplot,itarr0,spliny0 endif ; overplot drift+wobble solution if there is one if (wobblespline ne 0) then begin oplot,itarr0,splinx0+wosplinx0 oplot,itarr0,spliny0+wospliny0 endif ; add x and y labels labeloffset=(yrange[1]-yrange[0])*0.05 itlabel=fix(0.1*nshifts) xyouts,itlabel,shiftx[itlabel]+1.2*labeloffset,'x',charsize=1.5 xyouts,itlabel,shifty[itlabel]-labeloffset,'y',charsize=1.5 itlabel=fix(0.9*nshifts) xyouts,itlabel,shiftx[itlabel]+1.2*labeloffset,'x',charsize=1.5 xyouts,itlabel,shifty[itlabel]-labeloffset,'y',charsize=1.5 ; figtitle = main work dir + filename ; get current dir path; \ flag spaces in path if present (bloody Mac?) ; https://groups.google.com/forum/#!topic/comp.lang.idl-pvwave/Lo10H5XtN80 cd,current=thisdir thisdir=strjoin(strsplit(thisdir,' ',/extract),'\ ') workdir=file_basename(thisdir) xyouts,0.05,0.95,/norm,workdir+' '+psfilename closepsplot,psfile,opengv=1 endif if (verbose and outfile1 ne '') then $ print,' ----- sdo_shiftaligncubefiles wrote '+outfile1+' and '+outfile2 ; write the spline coefficients and shifts to file if requested if (n_elements(saveshifts) ne 0) then begin get_lun, unit_shifts openw,unit_shifts,saveshifts for it=itstart,itend do printf,unit_shifts,$ ntostr(shiftx0[it]),' ',ntostr(shifty0[it]),$ ' ',ntostr(splinx0[it]),' ',ntostr(spliny0[it]),$ ' ',ntostr(wosplinx0[it]),' ',ntostr(wospliny0[it]) free_lun,unit_shifts if (verbose) then print,' ----- selfalignfile wrote ',saveshifts endif ; output the shifts as variable for further use driftshifts=fltarr(nshifts0,2) wobbleshifts=fltarr(nshifts0,2) for it=0,nshifts0-1 do begin if (driftspline eq 0) then driftshifts[it,0:1]=[shiftx0[it],shifty0[it]] if (driftspline ne 0) then begin driftshifts[it,0:1]=[splinx0[it],spliny0[it]] if (wobblespline ne 0) then $ wobbleshifts[it,0:1]=[wosplinx0[it],wospliny0[it]] endif endfor end ; end of program ; =============== main for testing per IDLWAVE H-c ======================= ; test dir has brief 304, 171, 160, 1700, magnetogram, continuum, Dopplergram cubedir='/home/rutten/data/SDO/2014-06-14-small/target/cubes/' infile1=cubedir+'hmicont.fits' header_in=headfits_rr(infile1) nt_in=fxpar(header_in,'naxis3') ; construct shifts shifts=0.1*[[indgen(nt_in)],[indgen(nt_in)]] ; increasing per time step shifts=shifts+[[sin(indgen(nt_in)/2.)],[sin(indgen(nt_in)/2.)]] ; wobble shifts[12,0:1]=10.*shifts[12,0:1] ; add outlier shifts[16,0:1]=0.1*shifts[12,0:1] ; add outlier ; check ; print,shifts ; muck the input file muckfile1=cubedir+'hmimuck1.fits' muckfile=cubedir+'hmimuck.fits' reformcubefile,infile1,muckfile1,shift=[-10,-15],/splinip,missingvalue=1 reformcubefile,muckfile1,muckfile,shift=shifts,/splinip ; get alignment infile2=cubedir+'hmimuck.fits' outfile1='/tmp/hmicont_xal.fits' outfile2='/tmp/hmimuck_xal.fits' driftspline=1E6 wobblespline=1. driftspline=0 wobblespline=0 sdo_shiftaligncubefiles,infile1,infile2,$ driftshifts,wobbleshifts,$ outfile1=outfile1,outfile2=outfile2,$ ;; trange=[2,22],$ ;; alignmean=1,$ niter=5,cutiter=0.01,$ ; nxtile=10,tileprecision=0.001,$ driftspline=driftspline,wobblespline=wobblespline,nspliter=2,$ plotdir='/tmp/' ; check result if (outfile1 ne '') then showex,outfile2,outfile1 ; pull zoom, A jumps while B is stable end