; file: reformcube.pro = manipulate (x,y,t) cube in memory ; init: Jun 13 2017 Rob Rutten Deil from reformcubefile.pro ; last: Sep 25 2020 Rob Rutten Deil ;+ pro reformcube,incube,outcube,$ ; ---- startoff keywords for this program verbose=verbose,$ ; ---- most keywords of reformimage.pro for individual images xrange=xrange,yrange=yrange,trimbox=trimbox,$ greyborders=greyborders,greytriangles=greytriangles,$ xreverse=xreverse,yreverse=yreverse,$ rebinfactor=rebinfactor,congridfactor=congridfactor,$ splinip=splinip,$ sqrtint=sqrtint,logint=logint,bytscale=bytscale,flt2int=flt2int,$ nxlarge=nxlarge,nylarge=nylarge,$ cutcentralx=cutcentralx,cutcentraly=cutcentraly,$ shift=shift,rotate=rotate,scale=scale,missingvalue=missingvalue,$ despike=despike,sharpen=sharpen,flatten=flatten,$ histopttop=histopttop,histoptboth=histoptboth,smear=smear,$ muckdark=muckdark,muckbright=muckbright,muckmag=muckmag,muckosc=muckosc,$ ; ---- extra keywords for whole cube trange=trange,treverse=treverse,skip=skip,freezerange=freezerange,$ cubebytscale=cubebytscale,addzero=addzero,missingmean=missingmean,$ cubehistopttop=cubehistopttop,cubehistoptboth=cubehistoptboth,$ subtractrunmean=subtractrunmean,$ alignshift=alignshift,alignmetcalf=alignmetcalf,alignoutfile=alignoutfile ; NAME: ; reformcube.pro ; ; PURPOSE: ; emulates reformcubefile.pro for (smaller! cube in memory ; manipulate cube in memory usually image sequence) ; cut, shift, rescale, rotate, reverse, subsample, rebin, congrid, ; freeze a part ; self-align pairwise (only shifts or full Metcalf) ; take log or sqrt, bytscale, histo_opt, ; despike, sharpen, smear, subtract running mean, ; make grey monochrome borders and (derotation) edges, ; - and more ; ; USAGE EXAMPLES: ; extend duration (with freeze option) and field size (missingmean option) ; pixel-align to another cube (congrid, shift, rotate, resample, cut) ; change intensity scale from linear to log or sqrt ; histo_opt or bytscale whole cube or individual images ; despike or sharpen or smear ; subtract running mean ; grey monochrome borders and derotation triangles for better intscaling ; ; CALL: ; see above ; ; USAGE EXAMPLE: ; see tests below program end at end of this file ; ; INPUT: ; incube = 3D input array ; ; OUTPUT: ; outcube = 3D output array ; ; OPTIONAL KEYWORD INPUTS: ; ---- keywords for this program ; verbose = 2/1/0: print much, less, none ; ---- keywords of reformimage.pro used for individual image reform ; except croptriangles (since border widths may vary between frames) ; xrange, yrange = partial size cutout specification (2-element arrays) ; trimbox = 4-elem px box [xmin,ymin,xmax,ymax] to ignore eg triangle edges ; greyborders = 1/0: make edge borders grey = image mean (default 0) ; greytriangles = 1/0: also make triangular edges grey ; xreverse, yreverse = 1/0: reverse x or y axis ; rebinfactor = multiplier for rebin; integer [x,y] (not (nt,2) array) ; congridfactor = multiplier for congrid; float [x,y] (not (nt,2) array) ; splinip = 0: nearest neighbor (keep coarse pixels) ; = 1: cubic spline interpolation in poly_2d and congrid ; nxlarge,nylarge = increase field size (specify in new px after resize) ; shift, rotate, scale = inputs to reformimage.pro (=> shift_img.pro) ; shift = vector [shiftx,shifty] (time-constant) ; or (nt_in,2) array [shiftx,shifty] at intiming moments ; or (nt_out,2) array [shiftx,shifty] at outtiming moments ; values must be regridded new px sizes ; rotate = counterclockwise angle in degrees to rotate over ; number or array [nt_in] or array [nt_out] ; scale = vector [scalex,scaley] or array [nt_in,2] or [nt_out,2] ; missingvalue: use value for missing corners (see missingmean below) ; despike = value: value = sigma (e.g. IRIS 4.5) for particle hit removal ; sharpen = 1/0 (using fixed ImageMagic parameters) ; flatten = boxcar size in input px to flatten by subracting boxcar-smeared ; smear = smooth spatially with boxcar width = smear = width in old px ; histopttop: histo_opt value per image, uses /top_only ; histoptboth: histo_opt value per image, bottom and top ; muckdark: set intensities to mean below threshold expressed in mean=1 ; muckbright: set intensities to mean above threshold expressed in mean=1 ; muckmag = 1/0: muck HMI magnetogram into white only ; muckosc = 1/0: remove acoustics in AIA 1600 and 1700 images ; NB: muck order is: histopt, muck, smear ; sqrtint = 1/0: convert image intensities into sqrt(int) ; logint = 1/0: convert image intensities into alog10(int) ; cubebytscale = various options for bytscaling the image intensities ; cubebytscale=0: don't bytscale (retain type) - default ; cubebytscale=1: use (min,max) of first image for whole cube ; cubebytscale=2: use (min, max) of the whole image sequence ; cubebytscale=3: bytscale every image individually ; cubebytscale=4: use (min, max) from histo_opt_rr(mean(cube)) ; flt2int = 1/0: float cube > integer with max(abs(min),max) = 30000 ; cutcentralx = x-axis length of central image part to be cut out ; cutcentraly = y-axis length of central image part to be cut out ; ----- additional keywords for temporal business: ; trange = cut time span ; treverse = 1/0: reverse time axis ; skip = skip images (eg skip=1: use every second) ; trange and skip work only if no new time sampling is requested ; freezerange = set images during freezerange=[t1,t2] equal to image[t1] ; addzero = add this value to the zero level (eg +32768 for IRIS) ; missingmean = use cube mean for leftovers or non-image prior, post ; cubehistopttop = value: histo_opt top cutoff (eg 1E-3) for sequence mean ; cubehistoptboth = value: histo_opt cutoff (eg 1E-3) for sequence mean ; subtractrunmean = nr of preceding frames of which to subtract the mean ; alignshift = 1/0: self-align successive pairs by shifting 2nd ; alignmetcalf = 1/0: self-align successive pair per full Metcalf (SLOW!) ; alignoutfile = filestring to write self-alignment params per writecol ; ; METHOD: ; largely a wrapper to reformimage.pro but without keyword inheritance ; Uses the same keywords as reformimage.pro, use that for prior testing. ; ; RESTRICTIONS: ; ; HISTORY: ; Jun 13 2017 RR: start from reformcubefile.pro ; Sep 22 2020 RR: selfalignment ;- ; answer no-parameter query if (n_params() lt 2) then begin sp,reformcube return endif ; ---- top keywords if (n_elements(verbose) eq 0) then verbose=0 ; ---- keyword defaults of reformimage.pro ## keep the same if (n_elements(xrange) eq 0) then xrange=[0,-1] if (n_elements(yrange) eq 0) then yrange=[0,-1] if (n_elements(greyborders) eq 0) then greyborders=0 if (n_elements(greytriangles) eq 0) then greytriangles=0 if (n_elements(trimbox) eq 0) then trimbox=-1 if (n_elements(xreverse) eq 0) then xreverse=0 if (n_elements(yreverse) eq 0) then yreverse=0 if (n_elements(rebinfactor) eq 0) then rebinfactor=[0,0] if (n_elements(congridfactor) eq 0) then congridfactor=[0,0] if (n_elements(splinip) eq 0) then splinip=0 if (n_elements(nxlarge) eq 0) then nxlarge=0 if (n_elements(nylarge) eq 0) then nylarge=0 if (n_elements(shift) eq 0) then shift=[0,0] if (n_elements(rotate) eq 0) then rotate=0 if (n_elements(scale) eq 0) then scale=[1.0,1.0] if (n_elements(missingvalue) eq 0) then missingvalue=0 if (n_elements(despike) eq 0) then despike=0 if (n_elements(sharpen) eq 0) then sharpen=0 if (n_elements(flatten) eq 0) then flatten=0 if (n_elements(smear) eq 0) then smear=0 if (n_elements(histopttop) eq 0) then histopttop=0 if (n_elements(histoptboth) eq 0) then histoptboth=0 if (n_elements(muckdark) eq 0) then muckdark=0 if (n_elements(muckbright) eq 0) then muckbright=0 if (n_elements(muckmag) eq 0) then muckmag=0 if (n_elements(muckosc) eq 0) then muckosc=0 if (n_elements(sqrtint) eq 0) then sqrtint=0 if (n_elements(logint) eq 0) then logint=0 if (n_elements(bytscale) eq 0) then bytscale=0 if (n_elements(flt2int) eq 0) then flt2int=0 if (n_elements(addzero) eq 0) then addzero=0 if (n_elements(cutcentralx) eq 0) then cutcentralx=0 if (n_elements(cutcentraly) eq 0) then cutcentraly=0 ; ---- keyword defaults for temporal cube business if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(treverse) eq 0) then treverse=0 if (n_elements(skip) eq 0) then skip=0 if (n_elements(freezerange) eq 0) then freezerange=[0,-1] if (n_elements(missingmean) eq 0) then missingmean=0 if (n_elements(cubebytscale) eq 0) then cubebytscale=0 if (n_elements(cubehistopttop) eq 0) then cubehistopttop=0 if (n_elements(cubehistoptboth) eq 0) then cubehistoptboth=0 if (n_elements(subtractrunmean) eq 0) then subtractrunmean=0 if (n_elements(alignshift) eq 0) then alignshift=0 if (n_elements(alignmetcalf) eq 0) then alignmetcalf=0 if (n_elements(alignoutfile) eq 0) then alignoutfile='/tmp/selfalignout.dat' ; ========== start ; print pacifier if (verbose gt 0) then print,' ===== reformcube starts' ; set endian bigendian=1 ; get unique temporary file identifier for sharpening rannr=fix(abs(randomn(seed))*10000) strrannr=strmid(string(rannr+1E5,format='(i6)'),2) ; define temporary files for optional sharpening pngincube='/tmp/raw'+strrannr+'.png' pngoutcube='/tmp/sharp'+strrannr+'.png' ; ===== fits file: get cube geometry and file datatype from the fits header size_in=size(incube) if (size_in[0] ne 3) then begin print,' ##### reformcube aborts: incube not 3D cube' return endif nx_in=size_in[1] ny_in=size_in[2] nt_in=size_in[3] type_in=size_in(4) ; 1=byte, 2=int, 3=long, 4=float, 5=double if (type_in eq 1) then bitpix_in=8 if (type_in eq 2) then bitpix_in=16 if (type_in eq 3) then bitpix_in=32 if (type_in eq 4) then bitpix_in=-32 ; xrange and yrange xmin=xrange[0] xmax=xrange[1] ymin=yrange[0] ymax=yrange[1] if (xmax eq -1) then xmax=nx_in-1 if (ymax eq -1) then ymax=ny_in-1 if (trimbox[0] ne -1) then begin if (trimbox[0] gt xmin) then xmin=trimbox[0] if (trimbox[2] lt xmax) then xmax=trimbox[2] if (trimbox[1] gt ymin) then ymin=trimbox[1] if (trimbox[3] lt ymax) then ymax=trimbox[3] endif nx_out=xmax-xmin+1 ny_out=ymax-ymin+1 ; trange itstart=trange[0] itend=trange[1] if (itend eq -1) then itend=nt_in-1 nt_out=itend-itstart+1 ; redefine output image size if regridding requested ; rebinfactor and congridfactor must be fixed values ot dimensions change ; (use first pair if you have rescaling arrays, set scale around mean) resize=[0,0] pxratio=[1.,1.] if (rebinfactor[0] ne 0) then resize=rebinfactor ; 2-element vector if (congridfactor[0] ne 0) then resize=congridfactor ; 2-element vector if (resize[0] ne 0) then begin nxnew=fix((xmax-xmin+1)*resize[0]+0.5) nynew=fix((ymax-ymin+1)*resize[1]+0.5) nx_out=nxnew ny_out=nynew pxratio=1./resize ; 2-element vector endif ; redefine output image size if enlargement requested if (nxlarge ne 0) then nx_out=nxlarge if (nylarge ne 0) then ny_out=nylarge ; redefine output image size if central cutout requested if (cutcentralx ne 0) then nx_out=cutcentralx if (cutcentraly ne 0) then ny_out=cutcentraly ; =========== start output ; set bitpix_out = word type bitpix_out=bitpix_in if (cubebytscale ne 0) then bitpix_out=8 if (sharpen ne 0 and cubebytscale eq 0) then bitpix_out=16 if (flt2int) then bitpix_out=16 ; open output cube outcube=-1 if (bitpix_out eq 8) then outcube=bytarr(nx_out,ny_out,nt_out) if (bitpix_out eq 16) then outcube=intarr(nx_out,ny_out,nt_out) if (bitpix_out eq 32) then outcube=lonarr(nx_out,ny_out,nt_out) if (bitpix_out eq -32) then outcube=fltarr(nx_out,ny_out,nt_out) if (outcube[0] eq -1) then begin print,' ##### reformcube abort: no good bitpix_out' return endif ; find mean value of subcube for missingmean or subtractrunmean rescaling if (missingmean ne 0 or subtractrunmean gt 0) then begin sumim=fltarr(xmax-xmin+1,ymax-ymin+1) for it=itstart,itend do begin im=despike_rr(incube[*,*,it],despike,silent=(verbose ne 2)) sumim=sumim+float(im[xmin:xmax,ymin:ymax]) endfor cubemean=total(sumim)/((float(xmax)-xmin)*(ymax-ymin)*(itend-itstart)) if (missingmean ne 0) then missing=cubemean ; overrides keyword missing endif ; find (minint,maxint) for cubebytscale=1,2 or flt2int if (cubebytscale eq 2 or flt2int ne 0) then begin imstart=despike_rr(incube[*,*,itstart],despike,silent=(verbose ne 2)) imstart=imstart[xmin:xmax,ymin:ymax] minint=min(imstart) maxint=max(imstart) for it=itstart+1,itend do begin im=despike_rr(incube[*,*,it],despike,silent=(verbose ne 2)) im=im[xmin:xmax,ymin:ymax] mmim=minmax(im) if (mmim[0] lt minint) then minint=mmim[0] if (mmim[1] gt maxint) then maxint=mmim[1] endfor endif ; find histo_opted (minint,maxint) for cubebytscale=4 or cubehistop if (cubebytscale eq 4 or cubehistopttop ne 0 or cubehistoptboth ne 0) $ then begin imsum=float(despike_rr(incube[*,*,itstart],despike,silent=(verbose ne 2))) imsumt=imsum[xmin:xmax,ymin:ymax] for it=itstart+1,itend do begin im=float(despike_rr(incube[*,*,it],despike,silent=(verbose ne 2))) im=im[xmin:xmax,ymin:ymax] imsum=imsum+im endfor immean=imsum/float(itend-itstart+1) if (cubehistopttop ne 0) then $ immean=histo_opt_rr(immean,histopttop,/top_only) if (cubehistoptboth ne 0) then $ immean=histo_opt_rr(immean,histoptboth) if (cubehistopttop eq 0 and cubehistoptboth eq 0) then $ immean=histo_opt_rr(immean) mmmean=minmax(immean) minint=mmmean[0] maxint=mmmean[1] endif ;RR what do these do?? if (sqrtint) then begin cubescale=30000./sqrt(maxint) minint=cubescale*sqrt(minint) ; 30000 for integer cube maxint=cubescale*sqrt(maxint) ; 30000 for integer cube endif if (logint) then begin cubescale=30000./alog10(maxint) minint=cubescale*alog10(minint+0.1) ; 30000 for integer cube maxint=cubescale*alog10(maxint) ; 30000 for integer cube endif ;; if (flt2int) then begin ;RR does nothing?? ;; absmax=max([abs(minint),abs(maxint)]) ;; cubescale=30000./absmax ;; endif ; initialize first pair input images im_in_0=despike_rr(incube[*,*,itstart],despike,silent=(verbose ne 2)) im_in_1=despike_rr(incube[*,*,itstart+skip+1],despike,silent=(verbose ne 2)) im_prev=im_in_0 ; define selfalign output arrays alixshift=fltarr(nt_out) aliyshift=fltarr(nt_out) alirot=fltarr(nt_out) alixscale=fltarr(nt_out) aliyscale=fltarr(nt_out) ; start big loop over it_out timesteps = frame by frame processing ; ================================================================ it_in=itstart for it_out=0,nt_out-1 do begin ; no new timing = no interpolation if (it_in lt nt_in) then $ image=despike_rr(incube[*,*,it_in],despike,silent=(verbose ne 2)) $ else $ image=despike_rr(incube[*,*,it_in-1-skip],despike,silent=(verbose ne 2)) ; manipulate this image as requested ; ---------------------------------- ; selfalignment if (alignshift eq 1 or alignmetcalf eq 1) then begin shiftit=[0,0] rotateit=0 scaleit=[1.,1.] if (it_out gt 0) then begin if (alignshift eq 1) then begin shiftit=-findimshift_rr(im_prev,image,/subpix,filter=10) if (verbose ne 0) then $ print,' ----- alignshift it='+trimd(it_out)+$ ' shift ='+trimd(shiftit) endif if (alignmetcalf eq 1) then begin ;RR Oct 5 2020 I think this fails because uses scale instead of congrid below ;RR seems the ancient order problem solved wit iteratioin ;RR findalignimages.pro pin=[[0,0],[1.0,0]] qin=[[0,1.0],[0,0]] ;RR Metcalf "scale" parameter must be larger for worse starting guesses im1m=auto_align_images(im_prev,image,pin,qin,pmet,qmet,/quiet) pq2rss,pmet,qmet,erot,exscl,eyscl,exshft,eyshft,enrss,$ nx_out,ny_out,/center,quiet=1,/rotfirst shiftit=-[exshft,eyshft] rotateit=-erot scaleit=1./[exscl,eyscl] if (verbose ne 0) then $ print,' ----- alignmetcalf it='+trimd(it_out)+$ ' shift ='+trimd(shiftit)+$ ' rotate ='+trimd(rotateit)+$ ' scaleit ='+trimd(scaleit) endif endif ; copy into files for writeout alixshift[it_out]=shiftit[0] aliyshift[it_out]=shiftit[1] alirot[it_out]=rotateit alixscale[it_out]=scaleit[0] aliyscale[it_out]=scaleit[1] endif ; freeze movie if requested if (freezerange[1] ne -1) then begin if (it_out eq freezerange[0]) then frozen=image if (it_out gt freezerange[0] and it_out le freezerange[1]) then $ image=frozen endif ; set shift from shifts parameter shiftit=[-30000,-30000] if (n_elements(shift) eq 2) then shiftit=shift if (n_elements(shift) eq 2*nt_in) then begin xshiftit=shift[it_in,0] yshiftit=shift[it_in,1] shiftit=[xshiftit,yshiftit] endif if (shiftit[0] eq -30000) then begin print,' ##### reformcube abort: wrong shift definition' return endif ; set rotate from rotate parameter rotateit=-30000 if (n_elements(rotate) eq 1) then rotateit=rotate if (n_elements(rotate) eq nt_in) then rotateit=rotate[it_in] if (n_elements(rotate) eq nt_out) then $ rotateit=reform(rotate[it_out]) if (rotateit eq -30000) then begin print,' ##### reformcube abort: wrong rotate definition' return endif ; set scale from scale parameter scaleit=[-30000,-30000] if (n_elements(scale) eq 2) then scaleit=scale if (n_elements(scale) eq 2*nt_in) then begin xscaleit=scale[it_in,0] yscaleit=scale[it_in,1] scaleit=[xscaleit,yscaleit] endif if (n_elements(scale) eq 2*nt_out and nt_out ne nt_in) then $ scaleit=reform(scale[it_out,*]) ;RR ????? error? if (scaleit[0] eq -30000) then begin print,' ##### reformcube abort: wrong scale definition' return endif ; call reformimage with appropriate keywords reformimage,image,outimage,$ xrange=xrange,yrange=yrange,$ trimbox=trimbox,greyborders=greyborders,greytriangles=greytriangles,$ xreverse=xreverse,yreverse=yreverse,$ rebinfactor=rebinfactor,congridfactor=congridfactor,$ splinip=splinip,$ sqrtint=sqrtint,logint=logint,bytscale=0,flt2int=flt2int,$ nxlarge=nxlarge,nylarge=nylarge,$ cutcentralx=cutcentralx,cutcentraly=cutcentraly,$ shift=shiftit,rotate=rotateit,scale=scaleit,missingvalue=missingvalue,$ despike=despike,sharpen=sharpen,flatten=flatten,$ histopttop=histopttop,histoptboth=histoptboth,smear=smear,$ muckdark=muckdark,muckbright=muckbright,muckmag=muckmag,muckosc=muckosc ;; ; check: sv,image & sv,outimage flickrr,image,outimage ;; print, ' ====== reformcube: shift =',shiftit ;; STOP image=outimage im_prev=outimage ; previous for for selfalign ; cut at cubehistop values if requested if (cubehistopttop ne 0) then image=image(minint)