PRO bigsonic,first,last,dimx,dimy,bxdim,bydim ; ;Subsonic filtering in 3 dimensions (k_x,k_y,w), ;for large data volumes. ; ;The data volume is limited by the memory requirement ;of the 3-D array of filter (float), which is computed ;as a whole. Then, it is saved in a form of 2-D frames ;'filter.n', where n=first,...,last. ; ;Images (in SAVE2 format) are read, apodized in space and time, ;converted to float, and stored as unformatted files 1,2,3,... . ; ;For the direct FFT is used BIG3DFFT.PRO ;!! NOTICE !! Edit path in BIG3DFFT.PRO !! ;The FT transformed frames are stored as 'fft.n'. ; ;The 'fft' frames are multiplied by the filter frames, ;one by one. ; ;Then, the inverse FFT is done by BIG3DFFT.PRO, ;storing the frames inverse FT frames as 'F.n'. ; ;Finally, the 'F' frames are converted to integer and saved ;by SAVE2. ; ;Routines called: SAVE2, RESTORE2 (can be replaced by other ; read/write routines), NINT. ; ; ;INPUT PARAMETERS: ; ; first: number of the first image of the series ; last : number of the last image of the series ; dimx : x-dimension of the input images ; dimy : y-dimension of the input images ; bxdim: x-dimension of the sub-box in the Fourier domain partition ; bydim: y-dimension of the sub-box in the Fourier domain partition ; ; TO BE CHANGED IN EDITOR: ; PATHS ; xdim : x-dimension of the box to be extracted ; ydim : y-dimension of the box to be extracted ; x_anf: lower left corner of the box (x-dimension) ; y_anf: lower left corner of the box (y-dimension) ; scale: image scale arcsec/pixel ; ; PROMPTED: ; time_step: mean time distance of the images (seconds) ; v_ph : maximum phase velocity (km/s) ; ; OUTPUT: filtered images ; ; History: ; Some parts on the J.Hirzberger's procedure "wave.pro" ; and its adaptation "sonicf3d.pro" by M.Sobotka and J.A.Bonet ; are used. Written by M.S., Apr.18.2000 ; ; ;--------------------------------------------------------------------- ;----------------- PARAMETERS TO BE EDITED---------------------------- ; defining paths str0='/scratch/work/' ; general str1='/scratch/inpt_sav2.' ; input str2='/scratch/sonfilt.' ; output strf='/scratch/work/filter.' ; filter ; other parameters xdim=dimx ydim=dimy x_anf=0 y_anf=0 scale=0.12 ;arcsec/pixel ;------------------------------------------------------------------ ;------------------------------------------------------------------ if xdim mod 2 ne 0 then xdim=xdim-1 if ydim mod 2 ne 0 then ydim=ydim-1 ; ; reading input parameters ; print,'' read,'Mean time distance [s]: ', time_step read,'Max. Phase Velocity [km/s] : ', v_ph read,'Shall we apodize?: YES(=1), NO(=0): ',ap print,'Filter cutoff:' read,'abrupt (=0), smoothed WEDGE-MODE(=1), CONSTANT-MODE(=2) ',cut if (cut gt 2) or (cut lt 0) then begin print,'UNEXPECTED VALUE. Program terminated.' goto,fin endif if ap ne 0 or cut ne 0 then begin read,'length of the damped edge (% of total length)=',perc endif ;-------------------------------------------------------------------------- ; ; number of images should be even ; if (last-first+1) mod 2 ne 0 then last=last-1 tdim=last-first+1 print,last-first+1,' images to be filtered ;SUBSONIC FILTER CONSTRUCTION --------------------------------------------- ; ; defining unit in the Fourier domain ; print,'Spatial Resolution =',scale,' arcsec/pixel' print,'' kx_step=1./(scale*725.*xdim) ky_step=1./(scale*725.*ydim) w_step=1./(time_step*tdim) ; ; edge-smoothing parameters used for the filter and apodization ; if ap ne 0 or cut ne 0 then begin perct=perc*float(xdim)/float(tdim) ; % of damping in time smooth_x=fix(xdim*perc/100.) ; width of the edge in x-dimension smooth_y=fix(ydim*perc/100.) ; width of the edge in y-dimension smooth_t=fix(tdim*perct/100.) ; width of the edge in t-dimension endif ; ; prepares the filter (t-dimension is first here!) ; nx=xdim/2+1 ny=ydim/2+1 nt=tdim/2+1 filter_mask=fltarr(tdim,xdim,ydim) ;max.memory required here ; ; calculate subsonic filter with abrupt cutoff: ; if cut eq 0 or cut eq 1 then begin print,'Now calculating filter.' for j=0,ny-1 do begin for i=0,nx-1 do begin k_by_v=sqrt((i*kx_step)^2+(j*ky_step)^2)*v_ph for k=1,nt-1 do if k*w_step le k_by_v then filter_mask(k,i,j)=1. endfor endfor endif ; ;optional smoothing of the cutoff in the filter: WEDGE-MODE transition ; if cut eq 1 then begin ;optional smoothing of the cutoff. Wedge-mode five=nint(50./perct)+1 smoo=fltarr(smooth_t,smooth_t) for i=0,smooth_t-1 do begin for j=0,smooth_t-i-1 do begin smoo(i,j)=(1.+cos(!pi*float(j)/float(smooth_t-i)))/2. endfor endfor for i=0,nx-1 do begin for j=0,ny-1 do begin k=0 repeat begin if filter_mask(k+1,i,j) ne 1. then begin if k+1 ge five then trans=(k+1)/five else trans=1 if k ge 1 then $ filter_mask(k-trans+1:k-trans+smooth_t,i,j)=smoo(smooth_t-trans,*) k=nt-1 endif else begin k=k+1 endelse endrep until k eq nt-1 endfor endfor smoo=0 endif ; ; subsonic filter with optional smoothing of the cutoff in the filter: ; CONSTANT-WIDTH transition ; if cut eq 2 then begin print,'Now calculating filter.' ntrans=smooth_t-1 for j=0,ny-1 do begin for i=0,nx-1 do begin k_by_v=sqrt((i*kx_step)^2+(j*ky_step)^2)*v_ph for k=1,nt-1 do begin dif=k*w_step-k_by_v if dif le 0 then begin filter_mask(k,i,j)=1. endif else begin dif=dif/w_step in=0 for kk=k,k+ntrans do begin if kk le nt-1 then begin filter_mask(kk,i,j)=0.5+0.5*cos(!pi*(dif+in)/smooth_t) in=in+1 endif else begin goto,label endelse endfor goto,label endelse endfor label: endfor endfor endif ; ; calculate the other 7 octants of the filter ; filter_mask(*,*,ny:ydim-1)=reverse(filter_mask(*,*,1:ny-2),3) filter_mask(*,nx:xdim-1,*)=reverse(filter_mask(*,1:nx-2,*),2) filter_mask(nt:tdim-1,*,*)=reverse(filter_mask(1:nt-2,*,*),1) filter_mask(0,*,*)=1. ; ; writing filter frames ; for i=first,last do begin dcn=strtrim(i,2) openw,1,strf + dcn writeu,1,reform(filter_mask(i-first,*,*)) close,1 endfor print,'the filter was written to ' + strf filter_mask=0 ; to release memory space ; ;READING AND APODIZING IMAGES ----------------------------------------------- ; ; apodization functions (cosine window) ; if ap ne 0 then begin xmask=fltarr(xdim)+1. ymask=fltarr(ydim)+1. tmask=fltarr(tdim)+1. for i=0,smooth_x do xmask(i)=(1.-cos(!pi*float(i)/float(smooth_x)))/2. for i=0,smooth_y do ymask(i)=(1.-cos(!pi*float(i)/float(smooth_x)))/2. for i=0,smooth_t do tmask(i)=(1.-cos(!pi*float(i)/float(smooth_t)))/2. xmask(xdim-smooth_x:xdim-1)=reverse(xmask(1:smooth_x)) ymask(ydim-smooth_y:ydim-1)=reverse(ymask(1:smooth_y)) tmask(tdim-smooth_t:tdim-1)=reverse(tmask(1:smooth_t)) endif ; ; loop of reading, optional apodization, and writing of images ; for n=first,last do begin dcn=strtrim(n,2) print, 'apodizing image : ',str1 + dcn RESTORE2,str1+dcn,ima ima=float(ima(x_anf:x_anf+xdim-1,y_anf:y_anf+ydim-1)) ; ; apodization ; if ap ne 0 then begin av=mean(ima) ima=ima-av for i=0,xdim-1 do ima(i,*)=ima(i,*)*xmask(i) ; in x for i=0,ydim-1 do ima(*,i)=ima(*,i)*ymask(i) ; in y ima=ima*tmask(n-first) ; in t ima=ima+av endif ; ; writing ; openw,1,str0+dcn writeu,1,ima close,1 endfor print,'apodized images were written to ' + str0 xmask=0 ymask=0 tmask=0 ima=0 ; ;DIRECT FFT --------------------------------------------------------------- ; print,'' print,'Calling BIG3DFFT' BIG3DFFT,-1,first,last,dimx,dimy,bxdim,bydim,xdim,ydim,x_anf,y_anf ; ;MULTIPLY TRANSFORMED IMAGES BY FILTER IMAGES ----------------------------- ; ima=complexarr(xdim,ydim) filter=fltarr(xdim,ydim) for n=first,last do begin dcn=strtrim(n,2) openr,1,str0+'fft.'+dcn readu,1,ima close,1 openr,1,strf+dcn readu,1,filter close,1 print,'multiplying by the filter ',str0+'fft.'+dcn ima=ima*filter openw,1,str0+'fft.'+dcn writeu,1,ima close,1 endfor ima=0 filter=0 ; ;INVERSE FFT --------------------------------------------------------------- ; print,'' print,'Calling BIG3DFFT' BIG3DFFT,1,first,last,dimx,dimy,bxdim,bydim,xdim,ydim,x_anf,y_anf ; ;SAVING RESULTS ------------------------------------------------------------ ; ima=complexarr(xdim,ydim) for n=first,last do begin dcn=strtrim(n,2) openr,1,str0+'F.'+dcn readu,1,ima close,1 im=NINT(float(ima)) ;taking real part and converting to integer print,'saving image : ', str2 + dcn SAVE2,str2+dcn,im endfor fin: END