| FGG_3d(f,knots,N,accuracy,GridListx, GridListy, GridListz)
|
function F = FGG_3d(f,knots,N,accuracy,GridListx, GridListy, GridListz)
%Description:
%This code implements the 3D version of the "accelerated"
%Gaussian-gridding-based NUFFT described in Greengard and Lee [1]. The
%gridding approach is very similar to previous work by Nguyen and Liu [2];
%the only difference is the use of a different convolution kernel ([1]
%claims that the Gaussian kernel has computational advantages). Both
%algorithms allow the user to specify the numerical precision of the
%routine, but [1] provides a nice summary that would allow one to tabulate
%the appropriate variable values for each desired numerical precision. Code
%for [2] is also available upon request.
%
%This code performs NUFFTs that form rectangular cubes of size N=[Nx,Ny,Nz]
%(not counting frequency padding for image interpolation). The approximate
%DFT attains errors on the order of 1e-6 (for more or less accuracy, vary
%the optional "accuracy" parameter).
%
%Inputs:
% f: frequency-domain data (a complex Mx1 vector) unwrapped from a
% matrix knots: k-space locations at which the data were measured
% (an Mx3 vector). This data must be in double format.
% knots: the frequency locations of the data points. These locations
% should be normalized to correspond to the grid boundaries
% [-N/2, N/2 -1/N], N even. If the knots are not scaled
% properly, they will be shifted and scaled into this normalized
% form.
% N = [Nx,Ny,Nz]: the 1x3 vector denoting the size of the spatial
% grid in the image domain. This parameter will determine the
% spatial extent of the image.
% accuracy: (optional input parameter) a positive integer indicating
% the desired number of digits of accuracy
%Optional Inputs (highly recommended):
% GridListx:(column vector) The x-locations of the
% frequency grid (not yet scaled by c or 2*pi) onto which the
% data should be interpolated
% GridListy:(column vector) The y-locations of the
% frequency grid (not yet scaled by c or 2*pi) onto which the
% data should be interpolated
% GridListz:(column vector) The z-locations of the
% frequency grid (not yet scaled by c or 2*pi) onto which the
% data should be interpolated
%Outputs:
% F: the 3D NUFFT (approximate DFT) of f, with dimension [Nx,Ny,Nz].
%
%
%Usage Notes:
%In order for this function to work, the C
%file "FGG_Convolution3D.c" must be compiled into a Matlab executable
%(cmex) with the following command in the command prompt:
%
%mex FGG_Convolution3D.c
%
%A note on the effect of M_sp on the algorithm's accuracy:
%[R,M_sp]=[2,3] ==> 1e-3 accuracy
%[R,M_sp]=[2,6] ==> single precision
%[R,M_sp]=[2,9] ==> 1e-9 accuracy
%[R,M_sp]=[2,12] ==> double precision
%
%References:
%[1] L. Greengard and J. Lee, "Accelerating the Nonuniform Fast Fourier
%Transform," SIAM Review, Vol. 46, No. 3, pp. 443-454.
%[2] N. Nguyen and Q. H. Liu, "Nonuniform fast fourier transforms," SIAM J.
%Sci. Comput., 1999.
%
%Please note:
%This code is free to use, but we ask that you please reference the source,
%as this will encourage future funding for more free AFRL products. This
%code was developed through the AFOSR Lab Task "Moving-Target Radar Feature
%Extraction."
%Project Manager: Arje Nachman
%Principal Investigator: Matthew Ferrara
%Date: November 2008
%
%Code by (send correspondence to):
%Matthew Ferrara, Research Mathematician
%AFRL Sensors Directorate Innovative Algorithms Branch (AFRL/RYAT)
%Matthew.Ferrara@wpafb.af.mil
%Explanation of variables used
%bw The bandwidth of the input data (fmax-fmin)
%D A [Nx,Ny,Nz] deconvolution matrix formed from the 3D
% kronecker product of E4_x, E4_y, and E4_z
%dgx,dgy,dgz The separation (in the frequency domain) of the
% user-defined k-space grid onto which the data should be
% interpolated
%E1, E2, E3 factors of the Gaussian filter in the frequency domain
% (the Gaussian is factored to eliminate redundant
% exponential calculations)
%E4 The Gaussian deconvolution filter
%f The Mx1 vector of nonuniformly-spaced frequency data given
% as input to the NUFFT routine
%f_tau The [R*Nx,R*Ny] matrix of uniformly-spaced fequency-domain
% data values after "gridding"
%F_tau The FFT of f_tau
%F The approximate DFT of f (the deconvolved FFT of f_tau)
%fmean The average knot values of the input data (1x3 vector)
%j Index variable that indexes the convolution loop through
% the data points (1 <= j <= M)
%kmin The minimum knot values of the input data (1x3 vector)
%kmax The maximum knot values of the input data (1x3 vector)
%knots The Mx3 vector of frequency locations given as input to
% the NUFFT.
%M The number of (k-space) data points (the length of the
% input data vector f)
%M_sp Width of the frequency-domain box used in the approximate
% interpolation of each data point onto the frequency grid
% (M_sp=6 for single precision and M_sp=12 for double
% precision)
%N The image size of the NUFFT output (N=[Nx, Ny, Nz])
%Nx The length of the image in the x dimension
%Ny The length of the image in the y dimension
%Nz The length of the image in the z dimension
%R Oversampling ratio for gridding in the frequency (data)
% domain
%S The [Nz,Nx*Ny]-sized kronecker product of E4z, E4_x, and
% E4_y
%scale 1x3 vector used to scale the input data locations into the
% normalized form
%shift 1x3 vector used to shift the input data locations into the
% normalized form
%tau The 3x1 Gaussian kernel spreading factor
%End Explanation of variables
%% Step 1: Initialize constant variables:
M=length(f);%number of frequency-domain data points
if nargin<4, accuracy=6; end
if nargin<3, N=M; accuracy=6; end
if length(N)<2, N=N(1)*[1,1,1];accuracy=6; end
%The size parameters [Nx,Ny] determine the spatial extent of the image
Nx=N(1); Ny=N(2); Nz=N(3);
%R is the oversampling ratio(>1) for gridding in the frequency (data)
%domain. There are diminishing returns in accuracy after R=2 (M_sp has a
%more direct effect on accuracy).
R=2;
%M_sp is the length of the convolution kernel
M_sp=accuracy;%This gives roughly 6 digits of accuracy
%The variance, tau, of the Gaussian filter may be different in each
%dimension
%tau = M_sp./(N.^2);%I was initially using this value
tau = (pi*M_sp./(N.*N*R*(R-.5)));%Suggested value of tau by Greengard [1]
%The length of the oversampled grid
M_r = R*N;
%Scale knots (data locations) to [-N/2,N/2-1] (I don't initially use
%Greengard's [0,2*pi] convention, but instead assume users will most likely
%be thinking in terms of fs<=f<=fe)
kmin=min(knots);
kmax=max(knots);
if nargin <=4,%Simply choose the most convenient frequency-domain grid
bw=kmax-kmin;
scale=(N-1)./bw;
shift=-N/2-kmin.*scale;
knots=repmat(scale,[M,1]).*knots + repmat(shift,[M,1]);
else%we specify knot locations in terms of the user-defined grid (this
%consequently specifies the image pixel locations)
fmean=(kmin+kmax)/2;
dgx=GridListx(2)-GridListx(1);
dgy=GridListy(2)-GridListy(1);
dgz=GridListz(2)-GridListz(1);
%We choose to ensure that the data are perfectly centered on the grid:
kminx=fmean(1)-(Nx/2)*dgx;
kmaxx=kminx+(Nx-1)*dgx;
kminy=fmean(2)-(Ny/2)*dgy;
kmaxy=kminy+(Ny-1)*dgy;
kminz=fmean(3)-(Nz/2)*dgz;
kmaxz=kminz+(Nz-1)*dgz;
kmin=[kminx, kminy, kminz];
kmax=[kmaxx, kmaxy, kmaxz];
%The BW that covers the whole k-space region when confined to the
%specified grid spacing
bw=[kmaxx-kminx, kmaxy-kminy, kmaxz-kminz];
scale=(N-1)./bw;
shift=-.5*[Nx,Ny,Nz]-[kminx,kminy,kminz].*scale;
knots=repmat(scale,[M,1]).*knots + repmat(shift,[M,1]);
end
%Switch knot locations to [0,2*pi] convention (used by Greengard):
knots=mod(2*pi*knots./repmat(N,[M,1]),2*pi);%Makes NUFFT implementation
%more straightforward when notation is the same!
%Precompute E_3, the constant component of the (truncated) Gaussian:
E_3x(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(1)).^2)/tau(1));
%don't waste (slow) exponential calculations
E_3x=[fliplr(E_3x(1:(M_sp-1))),1,E_3x];
E_3y(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(2)).^2)/tau(2));
%don't waste (slow) exponential calculations
E_3y=[fliplr(E_3y(1:(M_sp-1))),1,E_3y];
E_3z(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(3)).^2)/tau(3));
%don't waste (slow) exponential calculations
E_3z=[fliplr(E_3z(1:(M_sp-1))),1,E_3z];
%Precompute E_4 (for devonvolution after the FFT)
kx_vec = (-Nx/2):(Nx/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4x(1:Nx,1)=sqrt(pi/tau(1))*exp(tau(1)*(kx_vec.^2));
ky_vec = (-Ny/2):(Ny/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4y(1,1:Ny)=sqrt(pi/tau(2))*exp(tau(2)*(ky_vec.^2));%
kz_vec = (-Nz/2):(Nz/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4z(1:Nz,1)=sqrt(pi/tau(3))*exp(tau(3)*(kz_vec.^2));%
%End initialization of constant variables
%% Step 2: Approximate convolution for each datum location, (x_j,y_j). This
% step is implemented in C and compiled into a Matlab-executable (cmex)
% file.
%Initialize convolved data matrix
f_tau=zeros(M_r(1),M_r(2),M_r(3));
f_taui=zeros(M_r(1)*M_r(2)*M_r(3),1);%Imaginary components of f_tau
f_taur=zeros(M_r(1)*M_r(2)*M_r(3),1);%Real components of f_tau
%Perform convolution onto finely-spaced grid
[f_taur,f_taui]=...
FGG_Convolution3D(double(real(f(:))),double(imag(f(:))),...
double(knots(:)),double(E_3x), double(E_3y), double(E_3z), ...
[double(M_sp), double(tau(1)), double(tau(2)), double(tau(3)),...
double(M_r(1)), double(M_r(2)), double(M_r(3))]);
f_tau = reshape(f_taur+sqrt(-1)*f_taui,[M_r(1), M_r(2), M_r(3)]);
%% Step 3: Perform FFT and deconvolve the result
%Perform FFT
F_tau=fftshift(fftn(ifftshift(f_tau)));
%Chop off excess pixels(the fine spacing in the frequency domain expanded
%the image in the transform domain)
F_tau(1:round(.5*(R-1)*Nx),:,:)=[];
F_tau(Nx+1:end,:,:)=[];
F_tau(:,1:round(.5*(R-1)*Ny),:)=[];
F_tau(:,Ny+1:end,:)=[];
F_tau(:,:,1:round(.5*(R-1)*Nz))=[];
F_tau(:,:,Nz+1:end)=[];
%Deconvolve the FFT by Hadamard multiplication with the Gaussian
%To form the deconvolution matrix S, we could do this with a loop:
%for i=1:Nx
% for j=1:Ny
% for k=1:Nz
% S(i,j,k)=E_4x(i)*E_4y(j)*E_4z(k);
% end
% end
%end
%HOWEVER, Matlab is extremely slow with loops so we instead use the kron()
%function and reshape the result (alternatively, we could have written
%another mex function):
S=kron(E_4z,E_4x*E_4y);%Scalars for 3D kronecker product (matrix must be
%reshaped)
D=permute(reshape(S,[Nx,Nz,Ny]),[1,3,2]);
F=F_tau.*D/(M*R*R*R); %Note that the scaling factor is 1/M, which is
%apparently not 1/M_r as shown in eqn (9) in [1]. This makes sense because
%the scaling factor for the DFT we are attempting to approximate (in eqn
%(1) of [1]) is 1/M.
%% If desired, compare to DFT
% F_true=zeros(Nx,Ny,Nz);
% for m=(-Nx/2):(Nx/2 -1)
% for n=(-Ny/2):(Ny/2 -1)
% for o=(-Nz/2):(Nz/2 -1)
% for k=1:M
% F_true(m+Nx/2+1,n+Ny/2+1,o+Nz/2+1)=...
% F_true(m+Nx/2+1,n+Ny/2+1,o+Nz/2+1)+...
% f(k)*exp(sqrt(-1)*(m*knots(k,1) + n*knots(k,2) + ...
% o*knots(k,3)));
% end
% end
% end
% F_true=F_true/M;
%
% figure
% isosurface(abs(F))
% title('F via NUFFT')
%
% figure
% isosurface(abs(F_true))
% title('F via DFT')
%
%error('done!')
|
|
