; file: shift_img_rr.pro = slight modification of Metcalf shift_img.pro ; init: Apr 25 2016 Rob Rutten Deil ; last: Mar 20 2022 Rob Rutten Deil ;+ function shift_img_rr,data,shifts,missing=missing, $ xscalein=xscalein,yscalein=yscalein, $ rot12in=rot12in,anchorin=anchorin,splinip=splinip ; Shift images in a data cube with optional scale and rotation. ; ;RR This RR modification: option to keep original pixelation ;RR RR modifications marked by ;RR comments ; ; INPUTS: ; data = image data cube (or a single image) ; shifts = array of x and y pixel shifts: fltarr(2,nimages). A ; positive value for the shifts moves the image up and ; right. ; ; KEYWORD INPUT PARAMETERS: ; missing = value for array for areas that shift out of the FOV ; (def = 0) ; xscale, yscale = scale change for images: fltarr(nimages) (def = 1.0) ; rot12 = rotate images CCW degrees: fltarr(nimages) (def = 0.0) ; anchor = center of rotation [x,y]: fltarr(2,nimages). ; (def = center of image) ;RR splinip = 0: nearest neighbour ;RR = 1: cubic spline interpolation with parameter -0.5 (Metcalf) ;RR = 2: bilinear interpolation ; ; OUTPUT: ; outdata = array of shifted-scaled-rotated images, same size as data ; ; PROCEDURE: ; Uses rss2pq to define the transformation and poly_2d to do the ; shifting. ; ; EXAMPLE: ; Here is an example of making a dejittered TRACE movie: ; itref = 0 ; itrace = trace_cube_pointing(tindex,tdata,itref) ; correct pointing ; xscale = tindex[itref].cdelt1/tindex.cdelt1 ; yscale = tindex[itref].cdelt2/tindex.cdelt2 ; xshift0 = ((tindex.naxis1-1)-(tindex[itref].naxis1-1)*xscale)/2. ; yshift0 = ((tindex.naxis2-1)-(tindex[itref].naxis2-1)*yscale)/2. ; xshift = -xshift0 - (tindex[itref].xcen-tindex.xcen)/tindex.cdelt1 ; yshift = -yshift0 - (tindex[itref].ycen-tindex.ycen)/tindex.cdelt2 ; shifts = transpose([[xshift],[yshift]]) ; dnew = shift_img(tdata,shifts) ; dejittered data cube ; This example does not correct for solar rotation and the ; pointing in the index structure will not apply to the shifted data. ; ; MODIFICATION HISTORY: ; T. Metcalf 2003-07-29 ; Mar 20 2022 RR: splinip > choice 0,1,2 commented ;RR ; splinip=1 gives overshoot negatives black outside AIA304 ; splinip=0 seems to introduce shift => added splinip=2 ;- ;RR default interpolation = nearest neighbor if (n_elements(splinip) eq 0) then splinip=0 dsize = size(data) if dsize[0] LT 2 OR dsize[0] GT 3 then $ message,'ERROR: data must have 2 or 3 dimensions' nx = dsize[1] ny = dsize[2] if dsize[0] EQ 3 then nimage=dsize[3] else nimage = 1L if n_elements(missing) NE 1 then missing = 0. xscale = fltarr(nimage) yscale = fltarr(nimage) rot12 = fltarr(nimage) anchor = fltarr(2,nimage) if keyword_set(xscalein) then begin if n_elements(xscalein) EQ nimage then xscale[*] = xscalein[*] $ else if n_elements(xscalein) EQ 1 then xscale[*] = xscalein $ else message,'ERROR: xscale must have either 1 element or '+string(nimage)+' elements' endif else xscale[*] = 1.0 if keyword_set(yscalein) then begin if n_elements(yscalein) EQ nimage then yscale[*] = yscalein[*] $ else if n_elements(yscalein) EQ 1 then yscale[*] = yscalein $ else message,'ERROR: yscale must have either 1 element or '+string(nimage)+' elements' endif else yscale[*] = 1.0 if keyword_set(rot12in) then begin if n_elements(rot12in) EQ nimage then rot12[*] = rot12in[*] $ else if n_elements(rot12in) EQ 1 then rot12[*] = rot12in $ else message,'ERROR: rot12 must have either 1 element or '+string(nimage)+' elements' endif else rot12[*] = 0.0 if keyword_set(anchorin) then begin if n_elements(anchorin) EQ 2*nimage then begin anchor[0,*] = anchorin[0,*] anchor[1,*] = anchorin[1,*] endif else if n_elements(anchorin) EQ 2 then begin anchor[0,*] = anchorin[0] anchor[1,*] = anchorin[1] endif else message,'ERROR: anchor must have either 2 elements or '+string(2*nimage)+' elements' endif else begin anchor[0,*] = (nx-1)/2. ;RR Dec 16 2017 funny top row at rotate=180 anchor[1,*] = (ny-1)/2. endelse outdata = data for i=0L,nimage-1L do begin xshift = shifts[0,i] yshift = shifts[1,i] t = rss2pq(nx,ny,xscale=xscale[i],yscale=yscale[i], $ xshift=xshift,yshift=yshift, $ rot12=rot12[i],p=p,q=q,anchor=anchor[*,i]) ;RR add splinip choice if (splinip eq 0) then $ outdata[*,*,i]=poly_2d(data[*,*,i],p,q,0,nx,ny,missing=missing) if (splinip eq 1) then $ outdata[*,*,i]=poly_2d(data[*,*,i],p,q,2,nx,ny,cubic=-0.5,missing=missing) if (splinip eq 2) then $ outdata[*,*,i]=poly_2d(data[*,*,i],p,q,1,nx,ny,missing=missing) endfor return,outdata end ; ==================== main test per IDLWAVE Hyper C ====================== ; SST file available https://robrutten.nl/rrweb/sdo-demo/instruction.html cd,'/media/rutten/RRHOME/alldata/SST/2014-06-24-solnet/sst' sstfile='wb.6563.corrected.aligned_cut.icube' im1full=lp_get(sstfile,1) szim1=size(im1full) nx=szim1[1] ny=szim1[2] im1=im1full[nx/4:nx*3./4,ny/4:ny*3./4] ; central cutout for smaller im2 xscalein=2 yscalein=2 shifts=[10.3,20.6] rot12in=0. missing=avg(im1) ;; ;RR classic all together: no good at all in findalignimages.pro ;; im2=shift_img_rr(im1,shifts,missing=missing, $ ;; xscalein=xscalein,yscalein=yscalein, $ ;; rot12in=rot12in,anchorin=anchorin,splinip=splinip) ;RR so unroll tests but beware: findalignimages has its own order ;; ; unroll: shift+rotate, scale = no good at all ;; im2a=shift_img_rr(im1,shifts,$ ;; rot12in=rot12in,$ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2=shift_img_rr(im2a,[0,0],$ ;; xscalein=xscalein,yscalein=yscalein, $ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; ; unroll: rotate, scale, shift = nearly good ; try in reformimage.pro? ;; im2a=shift_img_rr(im1,[0,0],$ ;; rot12in=rot12in,$ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2b=shift_img_rr(im2a,[0,0],$ ;; xscalein=xscalein,yscalein=yscalein, $ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2=shift_img_rr(im2b,shifts,$ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; ; unroll: rotate, shift, scale =; no good at all ;; im2a=shift_img_rr(im1,[0,0],$ ;; rot12in=rot12in,$ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2b=shift_img_rr(im2a,shifts,$ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2=shift_img_rr(im2b,[0,0],$ ;; xscalein=xscalein,yscalein=yscalein, $ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; ; unroll: rotate,shift+scale: no good at all ;; im2a=shift_img_rr(im1,[0,0],$ ;; rot12in=rot12in,$ ;; xscalein=1.,yscalein=1., $ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; im2=shift_img_rr(im2a,shifts,$ ;; xscalein=xscalein,yscalein=yscalein, $ ;; missing=missing,anchorin=anchorin,splinip=splinip) ;; showex,im1,im2,/blink px1=1 px2=1.05/avg([xscalein,yscalein]) angle=rot12in-2. findalignimages,im1full,im2,px1,px2,angle,$ px2asym,shiftfor,shiftrev,nxfor,nyfor,nxrev,nyrev,$ checkmetcalf=1,blink=0,show=3,verbose=1 ;; cd,'/home/rutten/data/SST/tests' ;; im=readfits('sst-demo-best.fits') ;; im[100,100]=min(im) ;; im[101,101]=max(im) ;; ;; im2=shift_img_rr(im,[3.3,3,3],xscale=3,yscale=3) ;;,/splinip) ;; im2=shift_img_rr(im,[1.3,1,7],/splinip) ;; ,xscale=3,yscale=3;;,/splinip) ;; showex,im,im2 end