function [f,iprep] = nufft1_3D(iprep,alpha,c2,c1,c3,X,Y,Z,m,q);
%
%
% NUFFT1 algorithm for forward 3D FFT of nonuniform time samples
%
% Input: iprep --- Preprocessing flag:
% iprep=0 if preprocessing has not been done;
% iprep=1 if preprocessing has already been done.
% alpha --- Input complex function to be transformed
% (c1,c2,c3) --- Nonuniform locations of tk in [-M/2,M/2],[-N/2,N/2]
% X,Y,Z --- Dimension of the uniform domain (after transform)
% m --- Over sampling rate (usually 2)
% q --- Number of interpolation points (even)
% Output: f --- Array of NUFFT1 output
% iprep --- Same meaning as input, but has changed value after
% initialization.
%
%==============================================================================
% Author: Jiayu Song
% Date: Feb 25, 2004
% Modified: Sep 30, 2006
%==============================================================================
%
global P1 P2 P3 mukx muky mukz sij
q = ceil(q);
mX = round(m*X/2)*2;
mY = round(m*Y/2)*2;
mZ = round(m*Z/2)*2;
n = length(c1);
if iprep == 0;
%ct1 = cputime;
iprep = pre_nufft1_3D(c1,c2,c3,X,Y,Z,m,q);
%ct2 = cputime;
%fprintf('CPU TIME for pre_nufft is %f s\n',ct2-ct1)
end;
clear c1 c2 c3;
% calculate Fourier coefficients
b = (ceil(-q/2):ceil(q/2))+mX/2;
c = (ceil(-q/2):ceil(q/2))+mY/2;
d = (ceil(-q/2):ceil(q/2))+mZ/2;
g1 = mod(repmat(mukx,1,q+1)+repmat(b,n,1),mX)+1;
clear mukx;
g2 = mod(repmat(muky,1,q+1)+repmat(c,n,1),mY)+1;
clear muky;
g3 = mod(repmat(mukz,1,q+1)+repmat(d,n,1),mZ)+1;
clear mukx;
fc = zeros(mX,mY,mZ);
t1=cputime;
G1 = zeros(q+1,q+1);
G = zeros(q+1,q+1,q+1);
for k=1:n
G1 = P1(:,k)*P2(k,:);
for i = 1: q+1
G(:,:,i) = G1*P3(k,i);
end
fc(g1(k,:),g2(k,:),g3(k,:)) = fc(g1(k,:),g2(k,:),g3(k,:))+alpha(k)*G;
end
clear g1 g2 g3 G1 G alpha P1 P2 P3
% calculate uniform FFT
% fc = fftshift(fc,1);
% fc = fftshift(fc,2);
% fc = fftshift(fc,3);
fc = fftshift(fc);
FFTs = ifftn(fc);
clear fc
FFTs = fftshift(FFTs);
% FFTs = fftshift(FFTs,1);
% FFTs = fftshift(FFTs,2);
% FFTs = fftshift(FFTs,3);
%scale FFT to approximate NUFFT
FFTs = FFTs((mX/2-X/2+1):(mX+X)/2,(mY/2-Y/2+1):(mY+Y)/2,(mZ/2-Z/2+1):(mZ+Z)/2);
f = FFTs.*sij*m*m*m;
%f = FFTs;
t2=cputime;
fprintf(['CPU TIME without pre-processing is', ...
' %f s\n'],t2-t1)