; file: rotatefitscube.pro ; init: Jun 18 2014 Rob Rutten Deil start from cube_rotate.pro ; last: Jul 27 2014 Rob Rutten Deil pro rotatefitscube,infile,outfile,angle=angle,$ fullfield=fullfield,missing=missing ;+ ; NAME: ; rotatefitscube.pro ; ; PURPOSE: ; rotate a fitscube, eg image cube = movie ; ; DESCRIPTION: ; rotates a fitscube file using assoc to permit large cube size ; ; CALL: ; rotatefitscube,infile,outfile,angle=angle,$ ; fullfield=fullfield,missing=missing ; ; MANDATORY INPUTS: ; infile: string with path/name input fits file = 3D data cube of any type ; eg movie [nx,ny,nt] ; outfile: string with path/name output fits file ; angle = angle, value to rotate over in degrees clockwise ; either a constant or a [nt] vector specifying angle per frame ; ; OPTIONAL INPUTS: ; fullfield = 1/0: extend array size to accommodate rotated corners ; missing = data value for rotated-in blank pixels ; default: missing=avg(incube) ; NB: missing=0 gives the default; use eg missing=0.001 to get black ; ; OUTPUTS: ; fits file with rotated cube, same type as incube, same size when ; /fullfield is not set, but symmetrically extended in x and y when ; /fullfield is set ; ; METHOD: ; uses assoc to permit cube sizes exceeding memory ; uses keyword for parameter angle to permit call by doallfiles.pro ; in order to process multiple similar files ; ; RESTRICTIONS: ; does not update rotation (WCS matrix) keywords ; ; HISTORY: ; Jun 19 2014 RR: start from cube_rotate.pro ; Jul 27 2014 RR: make suited for doallfiles.pro ;- ; answer a no-parameter query if (n_params() lt 2) then begin print,' -- rotatefitscube,infile,outfile,angle=angle,$' print,' fullfield=fullfield,missing=missing' print,' -- angle: degrees clockwise, either constant or [nt] vector' print,' -- optional fullfield: expand image size to accommodate corners' print,' -- optional missing = data value rotated-in blank px' print,' default = mean(incube)' print,' use missing = 0.001 or such to get black, not missing=0' return endif ; set endian bigendian=1 ; get cube geometry and file datatype from the fits header inheader=headfits_rr(infile) inheadersize=(1+fix(n_elements(inheader)/36.))*2880 bitpix=fxpar(inheader,'bitpix') nx=fxpar(inheader,'naxis1') ny=fxpar(inheader,'naxis2') nt=fxpar(inheader,'naxis3') ; check angle specification if (n_elements(angle) ne 1 and n_elements(angle) ne nt) then begin print,' ===== cube_rotate ERROR: angle not a constant nor [nt] vector' return endif ; open input file for assoc get_lun, unit_in if (bigendian) then openr,unit_in,infile,/swap_if_little_endian $ else openr,unit_in,infile if (bitpix eq -32) then inassoc=assoc(unit_in,fltarr(nx,ny),inheadersize) if (bitpix eq 16) then inassoc=assoc(unit_in,intarr(nx,ny),inheadersize) if (bitpix eq 8) then inassoc=assoc(unit_in,bytarr(nx,ny),inheadersize) ; default for missing (= grey value of pixels outside frame) if (n_elements(missing) eq 0) then missing=avg(inassoc[0]) ; extend image plane for fullfield option (no cutoff for rotated corners) if (n_elements(fullfield) eq 0) then begin nxout=nx nyout=ny dxmax=0 dymax=0 endif else begin dx=intarr(nt) dy=intarr(nt) semidiag=sqrt(float(nx-1)^2+float(ny-1)^2)/2. angdiag=acos((nx-1)/2./semidiag)*!radeg for it=0,nt-1 do begin dx[it]=fix(cos((angle[it]-angdiag)*!dtor)*semidiag-(nx-1)/2.+1) dy[it]=fix(sin((angle[it]+angdiag)*!dtor)*semidiag-(ny-1)/2.+1) if (n_elements(angle) eq 1) then break endfor dxmax=max(dx) dymax=max(dy) nxout=nx+2*dxmax nyout=ny+2*dymax endelse ; define outheader outheader=inheader outheadersize=inheadersize sxaddpar,outheader,'naxis1',nxout sxaddpar,outheader,'naxis2',nyout ; open output file for assoc, write header get_lun, unit_out if (bigendian) then openw,unit_out,outfile,/swap_if_little_endian $ else openw,unit_out,outfile if (bitpix eq -32) then $ outassoc=assoc(unit_out,fltarr(nxout,nyout),outheadersize) if (bitpix eq 16) then $ outassoc=assoc(unit_out,intarr(nxout,nyout),outheadersize) if (bitpix eq 8) then $ outassoc=assoc(unit_out,bytarr(nxout,nyout),outheadersize) if (outheadersize ne 0) then begin rec=assoc(unit_out, bytarr(outheadersize)) rec[0]=byte(outheader) endif ; rotate frame by frame ang=angle[0] for it=0,nt-1 do begin if (n_elements(angle) ne 1) then ang=angle[it] inimage=inassoc[it] if (bitpix eq -32) then newimage=fltarr(nxout,nyout)+float(missing) if (bitpix eq 16) then newimage=intarr(nxout,nyout)+fix(missing) if (bitpix eq 8) then newimage=bytarr(nxout,nyout)+byt(missing) newimage[dxmax:nx-1+dxmax,dymax:ny-1+dymax]=inimage[*,*] outassoc[it]=rot(newimage,ang,cubic=-0.5,missing=missing) endfor ; free the input and output files free_lun,unit_in,unit_out print,' ===== rotatefitscube wrote ',outfile end ; ================== main for IDLWAVE test ====================== infile='/home/rutten/data/SDO/2014-06-14-small/target/cubes/aia304.fits' outfile='/tmp/aia304_rot.fits' angle=30. rotatefitscube,infile,outfile,angle=angle ;; ,/fullfield showex,infile,outfile end