Python numpy.fft.fftn() Examples

def test_simple_plan_1d(self):
        """Test 1D plan with multiple batches"""
        shape = (10,)
        plan = self.G.create_plan(self.ctx, shape)
        plan.bake(self.queue)
        plan.inplace = False
        plan.layouts = (self.CLFFT_COMPLEX_INTERLEAVED, self.CLFFT_COMPLEX_INTERLEAVED)
        plan.batch_size = 2

        np_in = np.zeros(20, dtype=np.complex64)
        np_out = np.zeros(20, dtype=np.complex64)
        np_in.real = np.random.rand(20)

        cl_in = cl_array.to_device(self.queue, np_in)
        cl_out = cl_array.to_device(self.queue, np_out)

        plan.enqueue_transform(self.queue, cl_in.data, cl_out.data)
        self.queue.finish()
        out = fft.fftn(np_in[:10])
        self.assertTrue(np.allclose(out[:10], cl_out[:10].get())) 

def test_simple_plan_2d(self):
        shape = (8, 8)
        plan = self.G.create_plan(self.ctx, shape)
        plan.inplace = False
        plan.layouts = (self.CLFFT_COMPLEX_INTERLEAVED, self.CLFFT_COMPLEX_INTERLEAVED)
        plan.batch_size = 2
        plan.distances = (64, 64)
        plan.bake(self.queue)

        np_in = np.zeros(128, dtype=np.complex64)
        np_out = np.zeros(128, dtype=np.complex64)
        np_in[:64].real = np.random.rand(64)
        np_in[:64].imag = np.random.rand(64)
        np_in[64:].real = np.random.rand(64)
        np_in[64:].imag = np.random.rand(64)

        cl_in = cl_array.to_device(self.queue, np_in)
        cl_out = cl_array.to_device(self.queue, np_out)

        plan.enqueue_transform(self.queue, cl_in.data, cl_out.data)
        out = fft.fftn(np_in[64:].reshape(8, 8)).ravel()

        self.assertTrue(np.allclose(out, cl_out[64:].get())) 

def test_rfft_2d(self):
        shape = (8, 8)
        plan = self.G.create_plan(self.ctx, shape)
        plan.inplace = False
        plan.layouts = (self.CLFFT_REAL, self.CLFFT_HERMITIAN_INTERLEAVED)
        plan.batch_size = 2
        plan.distances = (64, 5 * 8)
        plan.strides_in = (8, 1)
        plan.strides_out = (8, 1)
        plan.bake(self.queue)

        np_in = np.random.rand(64 * 2).astype(np.float32)
        np_out = np.zeros(40 * 2, dtype=np.complex64)

        answer1 = fft.fftn(np_in[:64].reshape(8, 8))[:5]
        answer2 = fft.fftn(np_in[64:].reshape(8, 8))[:5]

        cl_in = cl_array.to_device(self.queue, np_in)
        cl_out = cl_array.to_device(self.queue, np_out)
        plan.enqueue_transform(self.queue, cl_in.data, cl_out.data)
        self.assertTrue(np.allclose(answer1, cl_out[:40].get().reshape(5, 8)))
        self.assertTrue(np.allclose(answer2, cl_out[40:].get().reshape(5, 8))) 

def run_lct(self, meas, X, Y, Z, S, x_max, y_max):
        self.max_dist = int(self.T * self.v/2 / self.channels[0])
        slope_x = self.dx * x_max / (self.fend/self.B * self.max_dist) * (1 + ((self.fstart/self.B)/(self.fend/self.B))**2)
        slope_y = self.dy * y_max / (self.fend/self.B * self.max_dist) * (1 + ((self.fstart/self.B)/(self.fend/self.B))**2)
        
        # params
        psf, fpsf = lct.getPSF(X, Y, Z, S, slope_x, slope_y)
        mtx, mtxi = lct.interpMtx(Z, S, self.fstart/self.B * self.max_dist, self.fend/self.B * self.max_dist)

        def pad_array(x, S, Z, X, Y):
            return np.pad(x, ((S*Z//2, S*Z//2), (Y//2, Y//2), (X//2, X//2)), 'constant')

        def trim_array(x, S, Z, X, Y):
            return x[S*int(np.floor(Z/2))+1:-S*int(np.ceil(Z/2))+1, Y//2+1:-Y//2+1, X//2+1:-X//2+1]

        invpsf = np.conj(fpsf) / (abs(fpsf)**2 + 1 / self.snr)
        tdata = np.matmul(mtx, meas.reshape((Z, -1))).reshape((-1, Y, X))

        fdata = fftn(pad_array(tdata, S, Z, X, Y))
        tvol = abs(trim_array(ifftn(fdata * invpsf), S, Z, X, Y))
        out = np.matmul(mtxi, tvol.reshape((S*Z, -1))).reshape((-1, Y, X))

        return out 

def getPSF(X, Y, Z, S, slope_x, slope_y):
    x = np.linspace(-1, 1, 2*X)
    y = np.linspace(-1, 1, 2*Y)
    z = np.linspace(0,2,2*S*Z)
    grid_z, grid_y, grid_x = np.meshgrid(z, y, x, indexing='ij')
    psf = np.abs(slope_x**2 * grid_x**2 + slope_y**2 * grid_y**2 - grid_z)

    psf = psf == np.tile(np.min(psf, axis=0, keepdims=True), (2*S*Z, 1, 1)) 
    psf = psf.astype(np.float32)
    psf = psf / np.sum(psf)

    psf = np.roll(psf, X, axis=2)
    psf = np.roll(psf, Y, axis=1)
    fpsf = fftn(psf)

    return psf, fpsf 

def lct(x1, y1, t1, v, vol, snr):
    X = len(x1)
    Y = len(y1)
    Z = len(t1)
    S = 2
    slope = np.max(x1) / (np.max(t1) * v/2)
    slope = np.max(y1) / (np.max(t1) * v/2)
    psf, fpsf = getPSF(X, Y, Z, S, slope)
    mtx, mtxi = interpMtx(Z, S, 0, np.max(t1)*v)

    def pad_array(x, S, Z, X):
        return np.pad(x, ((S*Z//2, S*Z//2), (X//2, X//2), (Y//2, Y//2)), 'constant')

    def trim_array(x, S, Z, X):
        return x[S*(Z//2)+1:-S*(Z//2)+1, X//2+1:-X//2+1, Y//2+1:-Y//2+1]

    invpsf = np.conj(fpsf) / (abs(fpsf)**2 + 1 / snr)
    tdata = np.matmul(mtx, vol)
    fdata = fftn(pad_array(tdata, S, Z, X))
    tvol = abs(trim_array(ifftn(fdata * invpsf), S, Z, X))
    vol = np.matmul(mtxi, tvol)
    return vol 

def dimension_reduction_randomized_pca_stack_diff(avg_vol_key, data):
  """
  :TODO: documentation
  :TODO: rename
  :TODO: move somewhere else
  """
  vol_avg = get_mrc(avg_vol_key)
  vol_avg_fft = fftshift(fftn(vol_avg))

  mat = np.zeros((len(data), vol_avg.size))

  for i, (vk, mk, ang, loc) in enumerate(data):
    # load vol/mask.
    vol  = get_mrc(vk)
    mask = get_mrc(mk)

    # rotate vol and mask according to angle/loc.
    vol_r  = core.rotate_vol_pad_mean(vol, ang, loc)
    mask_r = core.rotate_mask(mask, ang)

    v_dif = np.real(ifftn(ifftshift(fftshift(fftn(vol_r) - vol_avg_fft) * mask_r)))

    mat[i, :] = v_dif.flatten()

  return mat 

def test_UnscaledFFT_3d(backend, M, N, K, B ):
    b = backend()

    # forward
    x = b.rand_array( (M*N*K, B) )
    y = b.rand_array( (M*N*K, B) )
    x_h = x.to_host().reshape( (M,N,K,B), order='F' )

    A = b.UnscaledFFT( (M,N,K), dtype=x.dtype )

    A.eval(y, x)

    y_exp = np.fft.fftn( x_h, axes=(0,1,2) )
    y_act = y.to_host().reshape( (M,N,K,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2)

    # adjoint
    x = b.rand_array( (M*N*K, B) )
    y = b.rand_array( (M*N*K, B) )
    x_h = x.to_host().reshape( (M,N,K,B), order='F' )

    A.H.eval(y, x)

    y_exp = np.fft.ifftn( x_h, axes=(0,1,2) ) * (M*N*K)
    y_act = y.to_host().reshape( (M,N,K,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2) 

def test_UnscaledFFT_2d(backend, M, N, B ):
    b = backend()

    # forward
    x = b.rand_array( (M*N, B) )
    y = b.rand_array( (M*N, B) )
    x_h = x.to_host().reshape( (M,N,B), order='F' )

    A = b.UnscaledFFT( (M,N), dtype=x.dtype )

    A.eval(y, x)

    y_exp = np.fft.fftn( x_h, axes=(0,1) )
    y_act = y.to_host().reshape( (M,N,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2)

    # adjoint
    x = b.rand_array( (M*N, B) )
    y = b.rand_array( (M*N, B) )
    x_h = x.to_host().reshape( (M,N,B), order='F' )

    A.H.eval(y, x)

    y_exp = np.fft.ifftn( x_h, axes=(0,1) ) * (M*N)
    y_act = y.to_host().reshape( (M,N,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2) 

def test_UnscaledFFT_1d(backend, M, B ):
    b = backend()

    # forward
    x = b.rand_array( (M, B) )
    y = b.rand_array( (M, B) )
    x_h = x.to_host().reshape( (M,B), order='F' )

    A = b.UnscaledFFT( (M,), dtype=x.dtype )

    A.eval(y, x)

    y_exp = np.fft.fftn( x_h, axes=(0,) )
    y_act = y.to_host().reshape( (M,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2)

    # adjoint
    x = b.rand_array( (M, B) )
    y = b.rand_array( (M, B) )
    x_h = x.to_host().reshape( (M,B), order='F' )

    A.H.eval(y, x)

    y_exp = np.fft.ifftn( x_h, axes=(0,) ) * M
    y_act = y.to_host().reshape( (M,B), order='F' )
    npt.assert_allclose(y_act, y_exp, rtol=1e-2) 

def test_CenteredFFT(backend, M, N, K, B ):
    from numpy.fft import fftshift, ifftshift, fftn, ifftn

    b = backend()
    A = b.FFTc( (M,N,K), dtype=np.dtype('complex64') )

    # forward
    ax = (0,1,2)
    x = b.rand_array( (M*N*K,B) )
    y = b.rand_array( (M*N*K,B) )
    x_h = x.to_host().reshape( (M,N,K,B), order='F' )

    A.eval(y, x)

    y_act = y.to_host().reshape( (M,N,K,B), order='F' )
    y_exp = fftshift( fftn( ifftshift(x_h, axes=ax), axes=ax, norm='ortho'), axes=ax)
    npt.assert_allclose(y_act, y_exp, rtol=1e-2)

    # adjoint
    x = b.rand_array( (M*N*K,B) )
    y = b.rand_array( (M*N*K,B) )
    x_h = x.to_host().reshape( (M,N,K,B), order='F' )

    A.H.eval(y, x)

    y_act = y.to_host().reshape( (M,N,K,B), order='F' )
    y_exp = fftshift( ifftn( ifftshift(x_h, axes=ax), axes=ax, norm='ortho'), axes=ax)
    npt.assert_allclose(y_act, y_exp, rtol=1e-2) 

def test_plan_call():
    for shape in tested_shapes:
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.forward,
        )
        testing.assert_allclose(
            plan(),
            fft.fftn(plan.input_array)
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.fftn(plan.input_array) / plan.input_array.size
        )
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.backward
        )
        testing.assert_allclose(
            plan(),
            fft.ifftn(plan.input_array) * plan.input_array.size
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.ifftn(plan.input_array)
        ) 

def fftnc(x, axes):
    tmp = fft.fftshift(x, axes=axes)
    tmp = fft.fftn(tmp, axes=axes)
    return fft.ifftshift(tmp, axes=axes) 

def bench_random(self):
        from numpy.fft import fftn as numpy_fftn
        print()
        print('    Multi-dimensional Fast Fourier Transform')
        print('===================================================')
        print('          |    real input     |   complex input    ')
        print('---------------------------------------------------')
        print('   size   |  scipy  |  numpy  |  scipy  |  numpy ')
        print('---------------------------------------------------')
        for size,repeat in [((100,100),100),((1000,100),7),
                            ((256,256),10),
                            ((512,512),3),
                            ]:
            print('%9s' % ('%sx%s' % size), end=' ')
            sys.stdout.flush()

            for x in [random(size).astype(double),
                      random(size).astype(cdouble)+random(size).astype(cdouble)*1j
                      ]:
                y = fftn(x)
                #if size > 500: y = fftn(x)
                #else: y = direct_dft(x)
                assert_array_almost_equal(fftn(x),y)
                print('|%8.2f' % measure('fftn(x)',repeat), end=' ')
                sys.stdout.flush()

                assert_array_almost_equal(numpy_fftn(x),y)
                print('|%8.2f' % measure('numpy_fftn(x)',repeat), end=' ')
                sys.stdout.flush()

            print(' (secs for %s calls)' % (repeat))

        sys.stdout.flush() 

def potential(self, q, steps):
        hx = steps[0]
        hy = steps[1]
        hz = steps[2]
        Nx = q.shape[0]
        Ny = q.shape[1]
        Nz = q.shape[2]
        out = np.zeros((2*Nx-1, 2*Ny-1, 2*Nz-1))
        out[:Nx, :Ny, :Nz] = q
        K1 = self.sym_kernel(q.shape, steps)
        K2 = np.zeros((2*Nx-1, 2*Ny-1, 2*Nz-1))
        K2[0:Nx, 0:Ny, 0:Nz] = K1
        K2[0:Nx, 0:Ny, Nz:2*Nz-1] = K2[0:Nx, 0:Ny, Nz-1:0:-1] #z-mirror
        K2[0:Nx, Ny:2*Ny-1,:] = K2[0:Nx, Ny-1:0:-1, :]        #y-mirror
        K2[Nx:2*Nx-1, :, :] = K2[Nx-1:0:-1, :, :]             #x-mirror
        t0 = time.time()
        if pyfftw_flag:
            nthread = multiprocessing.cpu_count()
            K2_fft = pyfftw.builders.fftn(K2, axes=None, overwrite_input=False, planner_effort='FFTW_ESTIMATE',
                                       threads=nthread, auto_align_input=False, auto_contiguous=False, avoid_copy=True)
            out_fft = pyfftw.builders.fftn(out, axes=None, overwrite_input=False, planner_effort='FFTW_ESTIMATE',
                                          threads=nthread, auto_align_input=False, auto_contiguous=False, avoid_copy=True)
            out_ifft = pyfftw.builders.ifftn(out_fft()*K2_fft(), axes=None, overwrite_input=False, planner_effort='FFTW_ESTIMATE',
                                          threads=nthread, auto_align_input=False, auto_contiguous=False, avoid_copy=True)
            out = np.real(out_ifft())
        else:
            out = np.real(ifftn(fftn(out)*fftn(K2)))
        t1 = time.time()
        logger.debug('fft time:' + str(t1-t0) + ' sec')
        out[:Nx, :Ny, :Nz] = out[:Nx,:Ny,:Nz]/(4*pi*epsilon_0*hx*hy*hz)
        return out[:Nx, :Ny, :Nz] 

def bench_random(self):
        from numpy.fft import fftn as numpy_fftn
        print()
        print('    Multi-dimensional Fast Fourier Transform')
        print('===================================================')
        print('          |    real input     |   complex input    ')
        print('---------------------------------------------------')
        print('   size   |  scipy  |  numpy  |  scipy  |  numpy ')
        print('---------------------------------------------------')
        for size,repeat in [((100,100),100),((1000,100),7),
                            ((256,256),10),
                            ((512,512),3),
                            ]:
            print('%9s' % ('%sx%s' % size), end=' ')
            sys.stdout.flush()

            for x in [random(size).astype(double),
                      random(size).astype(cdouble)+random(size).astype(cdouble)*1j
                      ]:
                y = fftn(x)
                #if size > 500: y = fftn(x)
                #else: y = direct_dft(x)
                assert_array_almost_equal(fftn(x),y)
                print('|%8.2f' % measure('fftn(x)',repeat), end=' ')
                sys.stdout.flush()

                assert_array_almost_equal(numpy_fftn(x),y)
                print('|%8.2f' % measure('numpy_fftn(x)',repeat), end=' ')
                sys.stdout.flush()

            print(' (secs for %s calls)' % (repeat))

        sys.stdout.flush() 

def spectralcutoff(u, kappa, cx, cy, cz):
  nx = len(u[:,0,0])
  ny = len(u[0,:,0])
  nz = len(u[0,0,:])  
  nt= nx*ny*nz  
  uh = fftn(u)/nt  
  for i in range(0,nx):
    for j in range (0,ny):
      for k in range (0,nz):      
        rk = sqrt(cx*i*i + cy*j*j + cz*k*k)
        #rk = int(np.round(rk))
        if(rk >= kappa):
          uh[i,j,k] = 0.0
  ureal = ifftn(uh)*nt
  return ureal.real

#def spectralcutoff(u, kappa):
#  nx = len(u[:,0,0])
#  ny = len(u[0,:,0])
#  nz = len(u[0,0,:])  
#  nt= nx*ny*nz  
#  uh = fftn(u)/nt  
#  for i in range(0,nx):
#    for j in range (0,ny):
#      for k in range (0,nz):      
#        rk = sqrt(i*i + j*j + k*k)
#        #rk = int(np.round(rk))
#        if(rk >= kappa):
#          uh[i,j,k] = 0.0
#  ureal = ifftn(uh)*nt
#  return ureal.real 

def fftd(I, dims=None):

    # Compute fft
    if dims is None:
        X = fftn(I)
    elif dims == 2:
        X = fft2(I, axes=(0, 1))
    else:
        X = fftn(I, axes=tuple(range(dims)))

    return X 

def fftnc(x, axes, ortho=True):
    tmp = fft.fftshift(x, axes=axes)
    tmp = fft.fftn(tmp, axes=axes, norm="ortho" if ortho else None)
    return fft.ifftshift(tmp, axes=axes) 

def filter_given_curve(v, curve):
    grid = GV.grid_displacement_to_center(v.shape, GV.fft_mid_co(v.shape))
    rad = GV.grid_distance_to_center(grid)
    rad = N.round(rad).astype(N.int)
    b = N.zeros(rad.shape)
    for (i, a) in enumerate(curve):
        b[(rad == i)] = a
    vf = ifftn(ifftshift((fftshift(fftn(v)) * b)))
    vf = N.real(vf)
    return vf 

def fourier_transform(v):
    return fftshift(fftn(v)) 

def fsc(v1, v2, band_width_radius=1.0):
    
    siz = v1.shape
    assert(siz == v2.shape)

    origin_co = GV.fft_mid_co(siz)
    
    x = N.mgrid[0:siz[0], 0:siz[1], 0:siz[2]]
    x = x.astype(N.float)

    for dim_i in range(3):      x[dim_i] -= origin_co[dim_i]

    rad = N.sqrt( N.square(x).sum(axis=0) )

    vol_rad = int( N.floor( N.min(siz) / 2.0 ) + 1)

    v1f = NF.fftshift( NF.fftn(v1) )
    v2f = NF.fftshift( NF.fftn(v2) )

    fsc_cors = N.zeros(vol_rad)

    # the interpolation can also be performed using scipy.ndimage.interpolation.map_coordinates()
    for r in range(vol_rad):

        ind = ( abs(rad - r) <= band_width_radius )

        c1 = v1f[ind]
        c2 = v2f[ind]

        fsc_cor_t = N.sum( c1 * N.conj(c2) ) / N.sqrt( N.sum( N.abs(c1)**2 ) * N.sum( N.abs(c2)**2) )
        fsc_cors[r] = N.real( fsc_cor_t )

    return fsc_cors 

def apply_ctf(v, ctf):
    vc = N.real(     ifftn(  ifftshift(  ctf * fftshift(fftn(v)) ) )     )           # convolute v with ctf
    return vc 

def dog_smooth__large_map(v, s1, s2=None):

    if s2 is None:      s2 = s1 * 1.1       # the 1.1 is according to a DoG particle picking paper
    assert      s1 < s2

    size = v.shape


    pad_width = int(N.round(s2*2))
    vp = N.pad(array=v, pad_width=pad_width, mode='reflect')

    v_fft = fftn(vp).astype(N.complex64)
    del v;      GC.collect()


    g_small = difference_of_gauss_function(size=N.array([int(N.round(s2 * 4))]*3), sigma1=s1, sigma2=s2)
    assert      N.all(N.array(g_small.shape) <= N.array(vp.shape))       # make sure we can use CV.paste_to_whole_map()

    g = N.zeros(vp.shape)
    paste_to_whole_map(whole_map=g, vol=g_small, c=None)

    g_fft_conj = N.conj(   fftn(ifftshift(g)).astype(N.complex64)   )    # use ifftshift(g) to move center of gaussian to origin
    del g;      GC.collect()

    prod_t = (v_fft * g_fft_conj).astype(N.complex64)
    del v_fft;      GC.collect()
    del g_fft_conj;      GC.collect()

    prod_t_ifft = ifftn( prod_t ).astype(N.complex64)
    del prod_t;      GC.collect()

    v_conv = N.real( prod_t_ifft )
    del prod_t_ifft;      GC.collect()
    v_conv = v_conv.astype(N.float32)

    v_conv = v_conv[(pad_width+1):(pad_width+size[0]+1), (pad_width+1):(pad_width+size[1]+1), (pad_width+1):(pad_width+size[2]+1)]
    assert      size == v_conv.shape

    return v_conv 

def translation_align_given_rotation_angles(v1, m1, v2, m2, angs):
    v1f = fftn(v1)
    v1f[0,0,0] = 0.0
    v1f = fftshift(v1f)

    a = [None] * len(angs)
    for i, ang in enumerate(angs):
        v2r = GR.rotate_pad_mean(v2, angle=ang)
        v2rf = fftn(v2r)
        v2rf[0,0,0] = 0.0
        v2rf = fftshift(v2rf)

        m2r = GR.rotate_pad_zero(m2, angle=ang)
        m1_m2r = m1 * m2

        # masked images
        v1fm = v1f * m1_m2r
        v2rfm = v2rf * m1_m2r

        # normalize values
        v1fmn = v1fm / N.sqrt(N.square(N.abs(v1fm)).sum())
        v2rfmn = v2rfm / N.sqrt(N.square(N.abs(v2rfm)).sum())

        lc = translation_align__given_unshifted_fft(ifftshift(v1fmn), ifftshift(v2rfmn))

        a[i] = {'ang':ang, 'loc':lc['loc'], 'score':lc['cor']}

    return a 

def bench_random(self):
        from numpy.fft import fftn as numpy_fftn
        print
        print '    Multi-dimensional Fast Fourier Transform'
        print '==================================================='
        print '          |    real input     |   complex input    '
        print '---------------------------------------------------'
        print '   size   |  scipy  |  numpy  |  scipy  |  numpy '
        print '---------------------------------------------------'
        for size,repeat in [((100,100),100),((1000,100),7),
                            ((256,256),10),
                            ((512,512),3),
                            ]:
            print '%9s' % ('%sx%s'%size),
            sys.stdout.flush()

            for x in [random(size).astype(double),
                      random(size).astype(cdouble)+random(size).astype(cdouble)*1j
                      ]:
                y = fftn(x)
                #if size > 500: y = fftn(x)
                #else: y = direct_dft(x)
                assert_array_almost_equal(fftn(x),y)
                print '|%8.2f' % measure('fftn(x)',repeat),
                sys.stdout.flush()

                assert_array_almost_equal(numpy_fftn(x),y)
                print '|%8.2f' % measure('numpy_fftn(x)',repeat),
                sys.stdout.flush()

            print ' (secs for %s calls)' % (repeat)

        sys.stdout.flush() 

def fconv(x, otf):
    return np.real(ifftn(fftn(x) * otf)) 

def masked_difference_given_vol_avg_fft(v, m, ang, loc, vol_avg_fft, vol_mask_avg, smoothing_gauss_sigma=0):
  """
  :TODO: documentation

  # mxu: IMPORTANT: correction of dimension reduction procedure according to the paper  Heumann12

  :param v:
  :param m:
  :param ang:
  :param loc:
  :param vol_avg_fft:
  :param vol_mask_avg:
  :param smoothing_gauss_sigma:
  """
  v_r = core.rotate_vol_pad_mean(v,ang,loc)
  m_r = core.rotate_mask(m, ang)

  v_r_fft = fftshift(fftn(v_r))
  v_r_msk_dif = np.real(ifftn(ifftshift( (v_r_fft - vol_avg_fft) * m_r * vol_mask_avg )))

  if smoothing_gauss_sigma > 0:
    # mxu:
    # when the subtomograms are very noisy, PCA alone may not work well.
    # Because PCA only process vectors but does consider spetial
    # structures. In such case we can make use spatial information as an
    # additional layer of constrain (or approsimation model) by apply
    # certain amount of gaussian smoothing
    v_r_msk_dif = filtering.gaussian_smoothing(v=v_r_msk_dif, sigma=smoothing_gauss_sigma)

  return (v_r_msk_dif, m_r) 

def fourier_shell_correlation(v1, v2, bandwidth_radius = 1.0):
  """
  from: Min

  :todo: verify.
  """


  size = v1.shape

  assert( v1.shape == v2.shape )

  x = np.mgrid[0:size[0], 0:size[1], 0:size[2]]
  for i in range(3):
    x[i] -= size[i]//2

  R = np.sqrt(np.square(x).sum(axis=0))

  v1f = fftshift(fftn(v1))
  v2f = fftshift(fftn(v2))

  fsc = np.zeros(np.min(size)//2)

  for i in range(0,np.min(size)//2):
    mask = (abs(R-(i+1)) < 0.5)
    c1 = v1f[mask]
    c2 = v2f[mask]

    fsc[i] = np.real(sum(c1 * np.conj(c2))) / np.sqrt( sum(abs(c1)**2) * sum(abs(c2)**2))

  return fsc 

def gaussian_smoothing(v, sigma):
  """Smooth volume v, with spherical gaussian filter.

  Gaussian filter has single parameter sigma, which is symmetric standard
  deviation.

  Returns the volume v, after Gaussian smoothing is applied.
  """

  g = gauss_function(size=v.shape, sigma=sigma)
  # use ifftshift(g) to move center of gaussian to origin
  g_fft = fftn(ifftshift(g))
  v_conv = np.real( ifftn( fftn(v) * np.conj( g_fft ) ) )
  return v_conv 

def gaussian_smoothing(v, sigma):
 
  g = gauss_function(size=v.shape, sigma=sigma)


  g_fft = fftn(ifftshift(g));   # use ifftshift(g) to move center of gaussian to origin

  v_conv = np.real( ifftn( fftn(v) * np.conj( g_fft ) ) )

  return v_conv 

def vol_avg_fft_map(data, vol_shape, pass_dir):
  """
  Local work for computing average of all volumes.  This calculates a sum
  over the given subset of volumes.  Several of these are done and
  combined with vol_avg_global.  This is the map() stage of the global
  averaging step.

  :param vol_shape: The dimensions of the subtomograms we will be processing
  :param vmal_in: The data for incoming data.  A list of (volume, mask, angle, disp) tuples.
  :param v_out_key: The file location to write our local averaged volume to.
  :param m_out_key: The file location to write our local averaged mask to.
  """

  # temporary collection of local volume, and mask.
  vol_sum  = np.zeros(vol_shape, dtype=np.complex128, order='F')
  mask_sum = np.zeros(vol_shape, dtype=np.float64,  order='F')

  # iterate over all volumes/masks and incorporate data into averages.
  # volume_key, mask_key, angle offset, location offset.
  for vk, mk, ang, loc in data:

    # load vol/mask.
    vol  = get_mrc(vk)
    mask = get_mrc(mk)


    # rotate vol and mask according to angle/loc.
    vol  = core.rotate_vol_pad_mean(vol, ang, loc)
    mask = core.rotate_mask(mask, ang)

    vol_sum  += (fftshift(fftn(vol)) * mask)
    mask_sum += mask

  # volume/mask temporary accumulation locations.
  (v_fh, v_name) = tempfile.mkstemp(prefix='tm_tmp_vafmv_', suffix='.npy', dir=pass_dir)
  (m_fh, m_name) = tempfile.mkstemp(prefix='tm_tmp_vafmm_', suffix='.npy', dir=pass_dir)
  os.close(v_fh)
  os.close(m_fh)

  np.save(v_name, vol_sum)
  np.save(m_name, mask_sum)

  return v_name, m_name