pro sonicf3d,first,last,dimx,dimy ;,pore,scam ; ;subsonic filter in 3 dimensions (k_x,k_y,w) ; ; ;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 ; pore : 'name of the object' (STRING) ; scam: camera 'a' or 'b' (STRING) ; 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) ; time_step: mean time distance of the images ; v_ph : maximum phase velocity ; ; output : filtered images ; ; History: ; Based on the J.Hirzsberger's procedure "wave.pro" ; Adaptation made on 6,Nov,1995 by M.Sobotka and J.A.Bonet ; ;PARAMETERS TO BE EDITED------------------------------------------------- ; ; defining paths ; str1='/scratch/msobotka/a_ser23/dp4.' ; input path str2='/scratch/msobotka/a_ser23/dfp4.' ; output path ; str1='im'+scam+'/'+scam+'_d'+pore+'.' ; input path ; str2='im'+scam+'/'+scam+'_df'+pore+'.' ; output path scale=0.062 ;arcsec/pixel im=intarr(dimx,dimy) ; INPUT IMAGE ; ;------------------------------------------------------------------------- ; ; reading input parameters ; print,'' read,'Number of pixels in x-direction of the box: ', xdim read,'Number of pixels in y-direction of the box: ', ydim if xdim mod 2 ne 0 then xdim=xdim-1 if ydim mod 2 ne 0 then ydim=ydim-1 print,'' read,'Lower left corner x : ', x_anf read,'Lower left corner y : ', y_anf print,'' read,'Mean time distance : ', 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 ; ; reading images into imeq ; imeq=fltarr(tdim,xdim,ydim) for i=first,last do begin dcn=strtrim(i,2) print, 'Now reading image : ',str1 + dcn openr,1,str1 + dcn readu,1,im close,1 imeq(i-first,*,*)=float(im(x_anf:x_anf+xdim-1,y_anf:y_anf+ydim-1)) endfor im=0 ; to save memory ; ; 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) ; ; optional apodization with a cosine window ; if ap ne 0 or cut ne 0 then begin 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*perc/100.) ; width of the edge in t-dimension endif if ap ne 0 then begin av=mean(imeq) imeq=imeq-av xmask=fltarr(xdim)+1. ymask=fltarr(ydim)+1. tmask=fltarr(tdim)+1. print,'Now smoothing edges.' ; smoothing with a cosine 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)) for i=0,tdim-1 do imeq(i,*,*)=imeq(i,*,*)*tmask(i) for i=0,xdim-1 do imeq(*,i,*)=imeq(*,i,*)*xmask(i) for i=0,ydim-1 do imeq(*,*,i)=imeq(*,*,i)*ymask(i) imeq=imeq+av endif ; ; apply FFT ; print,'Now do FFT.' imeq=fft(imeq,-1,/overwrite) ; openw,1,'fft.lap' ; saving transformed array ; writeu,1,imeq ; close,1 ; ; prepares the filter ; nx=xdim/2+1 ny=ydim/2+1 nt=tdim/2+1 filter_mask=fltarr(tdim,xdim,ydim) ; ; 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./perc) 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 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. ; ; multiply transformed image by the filter ; imeq=imeq*filter_mask ; openw,1,'filter.lap' ; saving the filter ; writeu,1,filter_mask ; close,1 filter_mask=0 ;to save space ; ; apply the IFFT ; print,'Now do Inverse FFT:' imeq=FFT(imeq,1,/overwrite) imeq=nint(float(imeq)) ;taking real part and converting to integer ; ; saving filtered images ; for k=0,tdim-1 do begin dcn=strtrim(k+first,2) print, 'Now saving image : ', str2 + dcn openw,1,str2+dcn writeu,1,imeq(k,*,*) close,1 endfor ; fin: end