;+ ;NAME: ivm_desee ;PURPOSE: wrap-around program for ivm_destretchdeblur, to reduce ; effects of seeing variations in ivm data ;CATEGORY: image processing ;CALLING SEQUENCE: ivm_desee,fgin,fdin,hstr ;INPUTS: ; hstr: variable holding the IVM header structure, i.e. oblk ; gin: array(nn,nn, nimages) or (nn,nn,nf,nw*ns) of input geometry ; din: array(nn,nn, nimages) or (nn,nn,nf,nw*ns) of input data ; fgin: string, geometry-file name for input ; fdin: string, data-file name for input ; filein: file input flag ; fileout: file output flag ; vb: verbose flag -- 0 is errors-only, 1 is start/stop times, etc. ;OPTIONAL INPUT PARAMETERS: ;OPTIONAL KEYWORD PARAMETERS: ; smoothsize: integer, number of pixels for smoothing. Default is ; 2 * nn/256 + 1. Much more doesn't help, much smaller ; retains the pixellation noise this is trying to defeat. ; seeingout: if set, only computes/writes seeing array and does ; *not* continue with the deblur. Useful (nay, required) ; if using image-substitution but not deblur. ; nodeblur: passed to ivm_destretchdeblur, to do only destretch ; standard: output "standard" array ; err: is set when there is a problem with input files ;OUTPUTS: writes seeing-mitigated fgin and fdin files, if fileout ; is set. Outputs gout and dout arrays as 4D (see gin/din). ; also writes out seeing-array file and destretch/deblur ; coefficients array ;COMMON BLOCKS: none. ;SIDE EFFECTS: writes files same size as input files ;RESTRICTIONS: ;PROCEDURE: ;MODIFICATION HISTORY: written 9/01 KDL as upgrade from Rivm_v3 approaches ; 20091210 ELW - Array I/O, convert between 3D/4D data ; 20100106 ELW - Logical operators, not bitwise (AND=&&, NOT=~) ; 05/2010 KDL small commentary and testing for nterms being used for derivatives ; 20101220 ELW - Error flag keyword added ;- pro ivm_desee, hstr, filein, fileout, vb, gin=gin, din=din, dout=dout, gout=gout, $ rmssee=rmssee, fgin=fgin, fgout=fgout, fdin=fdin, fdout=fdout, $ seef=seef, smoothsize=smoothsize,seeingout=seeingout, nodeblur=nodeblur, $ standard=standard, err=err if (vb ge 1) then begin if keyword_set(seeingout) THEN $ PRINT,'ivm_desee start [seeingout set]: ',SYSTIME() $ else $ PRINT,'ivm_desee start: ',SYSTIME() endif error=0 ; Set-up nn = hstr.ppix nw = hstr.N_IMAGE_SETS ns = hstr.NTIMES nf = hstr.NFRAMES nimages=nf*nw*ns nwr=nw*ns nsmooth = 2 * nn/256 + 1 ;basic smoothing to counter pixellation noise. if (keyword_set(smoothsize)) then nsmooth = smoothsize nodblr=0 if (keyword_set(nodeblur)) then nodblr=1 ;first, construct reference image from sqrt(nimage) best ;according to seeing data, which is obtained here: if (vb ge 2) then PRINT,' ** computing seeing array **' rmssee = FLTARR(nimages) ; Use file input for geometry data if 'filein' set if (filein) then begin if ~keyword_set(fgin) then begin print, 'ivm_desee: ERROR! fgin filename not set. Stopping.' stop endif OPENR, lun1, fgin, error=err, /GET if (err ne 0) then begin print, 'ivm_desee: Exiting. Unable to open ', fgin return endif gin = FLTARR(nn,nn, nimages) READU, lun1, gin CLOSE, lun1 & FREE_LUN, lun1 endif else if ~KEYWORD_SET(gin) then begin print, 'ivm_desee: ERROR! gin array not set. Stopping.' stop endif else begin ; Make sure input data arrays are 3D for computation r = data_4Dto3D(gin, nn, nn, nf, nwr) if (r ne 1) then begin err=1 print, 'ivm_desee: unable to translate to 3D data.' & return endif endelse FOR i=0,nimages-1 DO BEGIN refimg = gin(*,*,i) ;rmssee(i) = SEERMS(refimg, NSMOOTH=nsmooth, /HIGH, /CENTER) rmssee(i) = SEERMS(refimg, /HIGH, /CENTER) ENDFOR IF KEYWORD_SET(seef) THEN BEGIN OPENW, lunx, seef, /GET_LUN WRITEU, lunx, rmssee CLOSE, lunx & FREE_LUN, lunx if (vb ge 2) then PRINT,' ** wrote '+seef+', = fltarr(',strtrim(nimages,2),')' endif ; If seeingout is set, we're done. if keyword_set(seeingout) THEN begin if (vb ge 2) then PRINT,' ** end [seeing array only]: ',SYSTIME() return endif if (vb ge 2) then PRINT,' ** Continuing with deblur **' ; Use file input for data if 'filein' set, otherwise use 'fdin' array if (filein) then begin if ~keyword_set(fdin) then begin print, 'ivm_desee: ERROR! fdin filename not set. Stopping.' stop endif din = FLTARR(nn,nn, nimages) OPENR, lun2, fdin,error=err, /GET if (err ne 0) then begin print, 'ivm_desee: Exiting. Unable to open ', fdin return endif READU, lun2, din CLOSE,lun2 & FREE_LUN,lun2 endif else if ~KEYWORD_SET(din) then begin print, 'ivm_desee: ERROR! din array not set. Stopping.' stop endif else begin ; Make sure input data arrays are 3D for computation r = data_4Dto3D(din, nn, nn, nf, nwr) if (r ne 1) then begin print, 'ivm_desee: unable to translate to 3D data.' err=1 return endif endelse ntiles=16. ; ELW-20100308 re-floated this ;ntiles=16 ; number of tiles in x,y direction, for building reference ; image. 32 is tilesize to start for 512^2 images, 16 is ; tilesize for 256^2, same patch of sky size for both. ; thought - perhaps a minimum variation in these sqrt(nimages)'s seeing ; should be instituted, so reference image doesn't get screwed up?? if (vb ge 2) then PRINT,' ** Constructing reference image **' ibest = REVERSE(SORT(rmssee) ) ibest = ibest(0:SQRT(nimages)-1) ibest = ibest(SORT(ibest)) nbest = N_ELEMENTS(ibest) dmid = ABS(ibest - nimages/2) ixmid = WHERE(dmid EQ MIN(dmid)) ixmid = ixmid(0) initref = gin(*,*,ibest(ixmid)) bigrefim = FLTARR(nn,nn,nbest) bigrefshifts = FLTARR(nn,nn,2,nbest) ; only correct for shifts for refim. FOR i=0,nbest-1 DO BEGIN ncoefs = IVM_DESTRETCHDEBLUR(ntiles,gin(*,*,ibest(i)),initref,newim,$ SMOOTH=nsmooth,nodeblur=nodblr) bigrefshifts(*,*,*,i) = ncoefs(*,*,1:2) ENDFOR ; in fact, need to subtract mean of bigrefshifts from bigrefshifts ; before interpolation. avdx = REFORM(REBIN(bigrefshifts(*,*,0,*),nn,nn,1,1)) avdy = REFORM(REBIN(bigrefshifts(*,*,1,*),nn,nn,1,1)) xx = REBIN( FINDGEN(nn), nn, nn ) yy = REBIN( TRANSPOSE( FINDGEN(nn) ), nn, nn ) FOR i=0,nbest-1 DO BEGIN bigrefshifts(*,*,0,i) = bigrefshifts(*,*,0,i)-avdx bigrefshifts(*,*,1,i) = bigrefshifts(*,*,1,i)-avdy bigrefim(*,*,i) = INTERPOLATE(gin(*,*,ibest(i)), xx-bigrefshifts(*,*,0,i), $ yy-bigrefshifts(*,*,1,i), CUBIC=-0.5 ) ENDFOR ; Average bigrefim into one. standard = REBIN(bigrefim,nn,nn) ; Start to do full data & geometry data. ; number of tiles in x,y direction, finer here than in computing reference above. ntiles = 2*ntiles ; use standard as ref for rest of images. if (vb ge 2) then PRINT,' ** Correcting seeing... **' outcoefs = FLTARR(ntiles,ntiles,4,nimages) ; WARNING 4 IS HARDWIRED FOR NTERMS=3 as in IVM_DESTRETCHDEBLUR KDL 05/2010 gout = fltarr(nn, nn, nimages) dout = fltarr(nn, nn, nimages) FOR k=0,nimages-1 DO BEGIN coefs = IVM_DESTRETCHDEBLUR(ntiles,gin(*,*,k),standard,newgeom, $ DATIM=din(*,*,k), NEWDATIM=newdat,SMOOTH=nsmooth,/SMALLOUT,nodeblur=nodblr) ;stop outcoefs(*,*,*,k) = coefs gout(*,*,k) = newgeom dout(*,*,k) = newdat ENDFOR if fileout then begin if keyword_set(fgout) then begin OPENW, lun3, fgout,/GET_LUN ; Write geom WRITEU,lun3, gout CLOSE,lun3 & FREE_LUN,lun3 endif else begin print, 'ivm_desee: ERROR! fgout filename not set. Stopping.' stop endelse if keyword_set(fdout) then begin OPENW, lun4, fdout,/GET_LUN ; Write data WRITEU,lun4, dout CLOSE,lun4 & FREE_LUN,lun4 endif else begin print, 'ivm_desee: ERROR! fdout filename not set. Stopping.' stop endelse datbase = STRMID(fdout,5,13) OPENW,lun5,'blur_coef_'+datbase,/GET_LUN WRITEU,lun5,outcoefs if (vb ge 1) then PRINT,' wrote blur_coef_'+datbase+', fltarr(',ntiles,ntiles,4,nimages,')' CLOSE,lun5 & FREE_LUN,lun5 endif ;======================== ; Set all arrays to 4D ;======================== if keyword_set(din) then r=data_3Dto4D(din, nn, nn, nf, nwr) if keyword_set(gin) then r=data_3Dto4D(gin, nn, nn, nf, nwr) if keyword_set(dout) then r=data_3Dto4D(dout, nn, nn, nf, nwr) if keyword_set(gout) then r=data_3Dto4D(gout, nn, nn, nf, nwr) if (vb ge 2) then print,' ** end: ',systime() END