pro large3d_fft,pm,first,last,dimx,dimy ; ;PURPOSE: ; Computes by parts the fft of a large 3D-array I(x,y,t) which cannot ; be handled in one go by the cpu. ; ;METHOD: ; Let us assume that we have a series of n=first-last+1 images of ; dimensions dimx*dimy, and we are interested in the fft of a sub-array ; of dimensions xdim*ydim*n. The first step consist in computing the ; fft of separated sub-images(xdim*ydim), thus creating in the hard ; disk a series of n complex Fourier-transformed-sub-images. ; From this series, 3D-complex-sub-arrays of dimensions bxdim*bydim*n, ; able to enter in the cpu in one go, are created, and then for each ; spatial position is extracted a temporal array and computed its 1D-fft. ; The procedure is based on the splitting of the 3D Fourier integral ; as follows: ; S{exp(-2*pi*j*nt*t)[SS{f(x,y,t)exp(-2*pi*j*(nx*x+ny*y))dxdx}]dt} ; ; WARNING! bxdim & bydim should be as larger as possible to minimize ; the intermediate read-write processes in the execution of the program. ; ;HISTORY: ; Edited by J.A.Bonet on Jul,7,1995. Modified on Nov,17,1995. ; ; ;INPUT (some of the inputs are interacting): ; ; pm : direction of the FFT: direct (=-1), invers (=1) ; first: number of the first image of the series ; last : number of the last image of the series ; dimx : x-dimension of the original images of the series ; dimy : y-dimension of the original images of the series ; xdim : x-dimension of the box extracted in the origi.images (even number) ; ydim : y-dimension of the box extracted in the origi.images (even number) ; x_inf: lower left corner of the box (x-dimension) ; y_inf: lower left corner of the box (y-dimension) ; bxdim: nu_x-dimension of the sub-box in the Fourier domain partition ; bydim: nu_y-dimension of the sub-box in the Fourier domain partition ; ;OUTPUTS: ; ; The resulting Fourier Transform consist on a series of last-first+1 ; separated complex images, directly stored on the hard disk, with ; names fft.n where n=0,...,first-last. ; Subsequent filtering processes requires sequencial manipulation (not ; in one go) of these complex images. ; The invers Fourier Transform will also produce a complex formatted ; array ;____________________________________________________________________________ ; ; defining paths TO BE CHANGED BY EDITOR ; str1=string('/home/stefans/images/') ; input path ; ; it is assumed that the filenames are 0, 1, 2, 3 etc. ;____________________________________________________________________________ ; selection: case pm of -1: begin strR=str1 strW=str1+'fft.' end 1: begin strR=str1+'fft.' strW=str1+'F' print,'' read,'x-dimension of the extracted box : ', xdim read,'y-dimension of the extracted box : ', ydim if xdim mod 2 ne 0 then xdim=xdim-1 if ydim mod 2 ne 0 then ydim=ydim-1 x_inf=0 & y_inf=0 goto, inv end else: begin read, 'valid pm values: -1 (forward), 1 (inverse) : ',pm print, pm goto, selection end endcase ; ; ; ; reading input parameters ; print,'' read,'x-dimension of the extracted box : ', xdim read,'y-dimension of the extracted 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 of the box: ', x_inf read,'Lower left corner y of the box: ', y_inf inv: print,'' read,'x-dimens.of the sub-box in the partition of the F. domain:',bxdim read,'y-dimens.of the sub-box in the partition of the F. domain:',bydim ; ;.............................................................................. ; ; reads separated images and computes their the 2D-Fourier Transform. ; for i=first,last do begin dcn=strtrim(i,1) if pm eq -1 then im=fltarr(dimx,dimy) else im=complexarr(xdim,ydim) print, 'Now reading image : '+ strR + dcn openr,1,strR + dcn readu,1,im close,1 ;window,i+6,xs=256,ys=256 ;tvscl,im(x_inf:x_inf+xdim-1,y_inf:y_inf+ydim-1) fim=fft(im(x_inf:x_inf+xdim-1,y_inf:y_inf+ydim-1),pm,/overwrite) print, '.... and writing its (x,y)-Fourier Transform' openw,1,strW + dcn writeu,1,fim close,1 ;window,i,xs=256,ys=256 ;img=float(fft(fim,-pm,/overwrite)) ;tvscl,img ;img=0 endfor fim=0 im=0 ;goto,fin ; ;creates 3d_subarrays (nu_x,nu_y,t) with dimensions in spatial freq. small ;enough to enter the 3d_subarray in the cpu, and computes the fft of ;1D temporal arrays correponding to separated x,y pixels ; a=xdim mod bxdim b=ydim mod bydim n=fix(xdim/bxdim) m=fix(ydim/bydim) if a ne 0 then begin dix=intarr(n+1)+bxdim dix(n)=a endif else begin dix=intarr(n)+bxdim endelse if b ne 0 then begin diy=intarr(m+1)+bydim diy(m)=b endif else begin diy=intarr(m)+bydim endelse nelx=fix(n_elements(dix)) nely=fix(n_elements(diy)) print, 'Number of sub-arrays = ',strtrim(nelx*nely,1) fim=complexarr(xdim, ydim) for jbox=0,nely-1 do begin for ibox=0,nelx-1 do begin box3d=complexarr(dix(ibox),diy(jbox),last-first+1) print, 'Now reading in the Fourier domain the sub-array: ',ibox,jbox for i=first,last do begin dcn=strtrim(i,1) openr,1,strW + dcn readu,1,fim close,1 box3d(*,*,i-first)=fim(ibox*bxdim:ibox*bxdim+dix(ibox)-1,$ jbox*bydim:jbox*bydim+diy(jbox)-1) endfor for y=0,diy(jbox)-1 do begin for x=0,dix(ibox)-1 do begin tt=box3d(x,y,*) box3d(x,y,*)=fft(tt,pm,/overwrite) endfor endfor print,'....and writing its t-Fourier transform' for i=first,last do begin dcn=strtrim(i,1) openr,1,strW + dcn readu,1,fim close,1 fim(ibox*bxdim:ibox*bxdim+dix(ibox)-1,$ jbox*bydim:jbox*bydim+diy(jbox)-1)=box3d(*,*,i-first) openw,1,strW + dcn writeu,1,fim close,1 endfor endfor endfor fin: print, '' print, '*** END OF PROGRAM ***' end