Python numpy.fft.fftn() Examples

The following are code examples for showing how to use numpy.fft.fftn(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: powerfit   Author: haddocking   File: test_gpyfft.py    Apache License 2.0 6 votes vote down vote up
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())) 
Example 2
Project: powerfit   Author: haddocking   File: test_gpyfft.py    Apache License 2.0 6 votes vote down vote up
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())) 
Example 3
Project: powerfit   Author: haddocking   File: test_gpyfft.py    Apache License 2.0 6 votes vote down vote up
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))) 
Example 4
Project: AcousticNLOS   Author: computational-imaging   File: AcousticNLOSReconstruction.py    MIT License 6 votes vote down vote up
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 
Example 5
Project: AcousticNLOS   Author: computational-imaging   File: lct.py    MIT License 6 votes vote down vote up
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 
Example 6
Project: AcousticNLOS   Author: computational-imaging   File: lct.py    MIT License 6 votes vote down vote up
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 
Example 7
Project: tomominer   Author: alberlab   File: funcs.py    GNU General Public License v3.0 6 votes vote down vote up
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 
Example 8
Project: indigo   Author: mbdriscoll   File: test_operators.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 9
Project: indigo   Author: mbdriscoll   File: test_operators.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 10
Project: indigo   Author: mbdriscoll   File: test_operators.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 11
Project: indigo   Author: mbdriscoll   File: test_operators.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 12
Project: fftw-cffi   Author: ghisvail   File: test_fftw.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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)
        ) 
Example 13
Project: mridata-recon   Author: MRSRL   File: fftc.py    GNU General Public License v3.0 5 votes vote down vote up
def fftnc(x, axes):
    tmp = fft.fftshift(x, axes=axes)
    tmp = fft.fftn(tmp, axes=axes)
    return fft.ifftshift(tmp, axes=axes) 
Example 14
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: bench_basic.py    GNU General Public License v3.0 5 votes vote down vote up
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() 
Example 15
Project: ocelot   Author: ocelot-collab   File: sc.py    GNU General Public License v3.0 5 votes vote down vote up
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] 
Example 16
Project: Computable   Author: ktraunmueller   File: bench_basic.py    MIT License 5 votes vote down vote up
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() 
Example 17
Project: TurboGenPY   Author: saadgroup   File: filters.py    MIT License 5 votes vote down vote up
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 
Example 18
Project: ProxImaL   Author: comp-imaging   File: utils.py    MIT License 5 votes vote down vote up
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 
Example 19
Project: dl-cs   Author: MRSRL   File: fftc.py    MIT License 5 votes vote down vote up
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) 
Example 20
Project: aitom   Author: xulabs   File: band_pass.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 21
Project: aitom   Author: xulabs   File: faml.py    GNU General Public License v3.0 5 votes vote down vote up
def fourier_transform(v):
    return fftshift(fftn(v)) 
Example 22
Project: aitom   Author: xulabs   File: vol.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 23
Project: aitom   Author: xulabs   File: ctf.py    GNU General Public License v3.0 5 votes vote down vote up
def apply_ctf(v, ctf):
    vc = N.real(     ifftn(  ifftshift(  ctf * fftshift(fftn(v)) ) )     )           # convolute v with ctf
    return vc 
Example 24
Project: aitom   Author: xulabs   File: gaussian.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 25
Project: aitom   Author: xulabs   File: util.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 26
Project: senior-design   Author: james-tate   File: bench_basic.py    GNU General Public License v2.0 5 votes vote down vote up
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() 
Example 27
Project: AcousticNLOS   Author: computational-imaging   File: ADMMReconstruction.py    MIT License 5 votes vote down vote up
def fconv(x, otf):
    return np.real(ifftn(fftn(x) * otf)) 
Example 28
Project: tomominer   Author: alberlab   File: worker_funcs.py    GNU General Public License v3.0 5 votes vote down vote up
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) 
Example 29
Project: tomominer   Author: alberlab   File: utils.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 30
Project: tomominer   Author: alberlab   File: filters.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 31
Project: tomominer   Author: alberlab   File: gaussian.py    GNU General Public License v3.0 5 votes vote down vote up
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 
Example 32
Project: tomominer   Author: alberlab   File: worker_funcs.py    GNU General Public License v3.0 4 votes vote down vote up
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