Perfect resampling of 3D Data at non-uniform points

Let f(x,y,z) is a set of discrete data in 3D, namely a volume image. I want to map this image onto a deformed configuration, i.e.:

g(xx,yy,zz)=f(x-u,y-v,z-w) where u,v and w are functions of space coordinates (x,y,z).

I don't want to use interpolation for this purpose as it leads to phase errors attributed to interpolation methods. Therefore, I need to perform perfect resampling in Fourier space, for instance.

In 2D, I do this as follows:

 DIM=size(f);

 F=fftshift(fft2(f));

 [x,y]=meshgrid(.5:DIM(2)-.5,.5:DIM(1)-.5);

 g=0;

 dim0=ceil((DIM+1)/2);

 for i=(1:DIM(1))-dim0(1)

    for j=(1:DIM(2))-dim0(2)

        m=i+dim0(1);

        n=j+dim0(2);

        g=g+F(m,n)*exp(2*pi*1i*(j*(x-u)/DIM(2)+i*(y-v)/DIM(1)));

    end
 end

 g=real(f);

However, when it comes to 3D, this will take extremely large computational time, which is almost impossible to run even on strong computers.

I am no expert in non-uniform discrete fourier transform (NDFT), anyone knows if this problem can be solved using NDFT more efficiently, or any other better idea which leads to perfect resampling in 3D?

Your ideas are very much appreciated

[Farhad M]

Sorry, in the last line:

 g=real(g);