PRO Conefilt_file, file, sx, sy, n, x, t, vel, $ RESFILE=resfile, PARTSIZE=partsize, PLOTOUT=plotout, _extra=extra ;+ ; NAME: ; CONEFILT_FILE ; PURPOSE: ; k-w filtering of a 2D time series (on disk) ; CALLING SEQUENCE: ; Conefilt_file, datafile, sx, sy, num, dx, dt, Vel ; INPUTS: ; Datafile: (string) name of a file with the datacube (integer!) ; SX, SY: (integer) (spatial) size of the images ; num: (integer) number of images in time ; dx, dt: (float) spatial and temporal binsize ; vel: (num) cutoff (in same units as dx/dt) ; OPTIONAL PARAMETERS: ; ; KEYWORDS: ; RESFILE: (string) filename for the filtered data, default sum.dat ; PARTSIZE: (2-el. integer) size of the sub-parts, default 384x256 ; PLOTOUT: switch to have images on screen or not (default yes) ; _EXTRA: _extra keywords of conefilt,pro ; OUTPUTS: ; ; COMMON BLOCKS: ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; PROCEDURE: ; ; MODIFICATION HISTORY: ; 20-Oct-1999 P.Suetterlin, SIU ; 26-Oct-1999 First and last image are somehow "smeared". Try ; to put an artificial interpolated image before ; and after the series. ; 06-Mar-2001 automatic split- and rejoin, use conefilt for parts. ; 25-Feb-2002 It could happen that parts do not overlapp. Fixed. ; Some error checking. ; 16-Apr-2004 Add the (pass-through) keywords for normal conefilt ;- IF n_params() LT 7 THEN BEGIN message, "Use: Conefilt_file, datafile, sx, sy, num, dx, dt, Vel", /info return ENDIF IF NOT keyword_set(resfile) THEN resfile = 'sum.dat' IF NOT keyword_set(partsize) THEN BEGIN dx = 384 dy = 256 ENDIF ELSE BEGIN IF n_elements(partsize) EQ 1 THEN BEGIN dx = partsize dy = partsize ENDIF ELSE BEGIN dx = partsize(0) dy = partsize(1) ENDELSE ENDELSE IF NOT keyword_set(plotout) THEN plotout = 0 IF (sx LT dx) OR (sy LT dy) THEN BEGIN message, "Part dimensions must be smaller than image dimensions!" return ENDIF ;;; ;;; Compute the number of subfields ;;; special case: exactly one subimage fits in, then stay with n=1, ;;; else add one to assure an overlapp nx = float(sx)/dx IF nx EQ 1. THEN nx = 1 ELSE nx = fix(nx+1) ny = float(sy)/dy IF ny EQ 1. THEN ny = 1 ELSE ny = fix(ny+1) ;;; The startpoints of the subfields in X and Y xsta = intarr(nx) ysta = intarr(ny) xsta(0) = 0 FOR i=1, nx-2 DO xsta(i) = i*(sx-dx)/(nx-1) xsta(nx-1) = sx-dx ysta(0) = 0 FOR i=1, ny-2 DO ysta(i) = i*(sy-dy)/(ny-1) ysta(ny-1) = sy-dy mov = intarr(dx, dy, n) print, 'Splitting into '+strtrim(nx*ny, 2)+' parts' pa = 1 openr, unit00, file, /get p = assoc(unit00, intarr(sx, sy)) get_lun, unit01 FOR i=0, nx-1 DO BEGIN FOR j=0, ny-1 DO BEGIN writeu, -1, 'Part '+nnumber(pa, 2)+$ string(xsta(i), xsta(i)+dx-1, ysta(j), ysta(j)+dy-1, $ form="('[',i4,',',i4,',',i4,',',i4,'] ....')") FOR k=0, n-1 DO $ mov(*, *, k)=p(xsta(i):xsta(i)+dx-1, ysta(j):ysta(j)+dy-1, k) conefilt, mov, x, t, vel, _extra=extra openw, unit01, 'part.'+nnumber(pa, 2) writeu, unit01, mov close, unit01 print, ' done' pa = pa+1 ENDFOR ENDFOR close, unit00 openr, unit01, 'part.01' FOR i=2, nx*ny DO $ OK=execute("openr, unit"+nnumber(i, 2)+", 'part."+nnumber(i, 2)+"',/get") FOR i=1, nx*ny DO $ OK=execute("p"+nnumber(i, 2)+"=assoc(unit"+nnumber(i, 2)+",intarr(dx,dy))") openw, unit00, resfile FOR img=0, n-1 DO BEGIN pa = 1 wp1 = fltarr(sx, sy) pp = fltarr(sx, sy) FOR i=0, nx-1 DO BEGIN FOR j=0, ny-1 DO BEGIN cmd = 'paste,pp,p'+nnumber(pa, 2)+'(img),wp1,xsta(i),ysta(j)' ok = execute(cmd) pa = pa+1 ENDFOR ENDFOR pp = fix(pp/wp1) writeu, unit00, pp if (plotout) then tvscl, rescale(pp, sc=.5) ENDFOR close, unit00 FOR i=0, nx*ny DO ok = execute("free_lun, unit"+nnumber(i, 2)) END