Python fftpack.ifftn函数代码示例

    def psf_calc(self, psf, kz, data_size):
        '''Pre calculate OTFs etc ...'''
        g = psf;

        self.height = data_size[0]
        self.width  = data_size[1]
        self.depth  = data_size[2]

        (x,y,z) = mgrid[-floor(self.height/2.0):(ceil(self.height/2.0)), -floor(self.width/2.0):(ceil(self.width/2.0)), -floor(self.depth/2.0):(ceil(self.depth/2.0))]
        
       
        gs = shape(g);
        g = g[int(floor((gs[0] - self.height)/2)):int(self.height + floor((gs[0] - self.height)/2)), int(floor((gs[1] - self.width)/2)):int(self.width + floor((gs[1] - self.width)/2)), int(floor((gs[2] - self.depth)/2)):int(self.depth + floor((gs[2] - self.depth)/2))]
	
        g = abs(ifftshift(ifftn(abs(fftn(g)))));
        g = (g/sum(sum(sum(g))));
	
        self.g = g;
        
        self.H = cast['f'](fftn(g));
        self.Ht = cast['f'](ifftn(g));

        tk = 2*kz*z
        
        t = g*exp(1j*tk)
        self.He = cast['F'](fftn(t));    
        self.Het = cast['F'](ifftn(t));

        tk = 2*tk
        
        t = g*exp(1j*tk)
        self.He2 = cast['F'](fftn(t));    
        self.He2t = cast['F'](ifftn(t));

def fft_correlation(in1, in2, normalize=False):
    """Correlation of two N-dimensional arrays using FFT.

    Adapted from scipy's fftconvolve.

    Parameters
    ----------
    in1, in2 : array
    normalize: bool
        If True performs phase correlation

    """
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    size = s1 + s2 - 1
    # Use 2**n-sized FFT
    fsize = 2 ** np.ceil(np.log2(size))
    IN1 = fftn(in1, fsize)
    IN1 *= fftn(in2, fsize).conjugate()
    if normalize is True:
        ret = ifftn(np.nan_to_num(IN1 / np.absolute(IN1))).real.copy()
    else:
        ret = ifftn(IN1).real.copy()
    del IN1
    return ret

 def invert(self, estimate):
     """Invert the estimate to produce slopes.
     
     Parameters
     ----------
     estimate : array_like
         Phase estimate to invert.
     
     Returns
     -------
     xs : array_like
         Estimate of the x slopes.
     ys : array_like
         Estimate of the y slopes.
     
     
     """
     if self.manage_tt:
         estimate, ttx, tty = remove_tiptilt(self.ap, estimate)
     
     est_ft = fftpack.fftn(estimate) / 2.0
     
     xs_ft = self.gx * est_ft
     ys_ft = self.gy * est_ft
     
     xs = np.real(fftpack.ifftn(xs_ft))
     ys = np.real(fftpack.ifftn(ys_ft))
     
     if self.manage_tt and not self.suppress_tt:
         xs += ttx
         ys += tty
     
     return (xs, ys)

 def test_definition(self):
     x = [[1,2,3],[4,5,6],[7,8,9]]
     y = ifftn(x)
     assert_array_almost_equal(y,direct_idftn(x))
     x = random((20,26))
     assert_array_almost_equal(ifftn(x),direct_idftn(x))
     x = random((5,4,3,20))
     assert_array_almost_equal(ifftn(x),direct_idftn(x))

 def test_definition(self):
     x = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=self.dtype)
     y = ifftn(x)
     assert_(y.dtype == self.cdtype)
     assert_array_almost_equal_nulp(y,direct_idftn(x),self.maxnlp)
     x = random((20,26))
     assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp)
     x = random((5,4,3,20))
     assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp)

    def test_invalid_sizes(self):
        with assert_raises(ValueError,
                           match="invalid number of data points"
                           r" \(\[1 0\]\) specified"):
            ifftn([[]])

        with assert_raises(ValueError,
                           match="invalid number of data points"
                           r" \(\[ 4 -3\]\) specified"):
            ifftn([[1, 1], [2, 2]], (4, -3))

def convolve_turb(image, fwhm, get_psf=False):
    """
    Convolve the input image with a turbulent psf

    parameters
    ----------
    image:
        A numpy array
    fwhm:
        The FWHM of the turbulent psf.
    get_psf:
        If True, return a tuple (im,psf)

    The images dimensions should be square, and even so the psf is
    centered.
    """

    dims = array(image.shape)
    if dims[0] != dims[1]:
        raise ValueError("only square images for now")

    # add padding for PSF in real space
    # sigma is approximate
    kdims=dims.copy()
    kdims += 2*4*fwhm/TURB_SIGMA_FAC

    # Always use 2**n-sized FFT
    kdims = 2**ceil(log2(kdims))
    kcen = kdims/2.

    imfft = fftn(image,kdims)

    k0 = 2.92/fwhm
    # in fft units
    k0 *= kdims[0]/(2*pi)

    otf = pixmodel.ogrid_turb_kimage(kdims, kcen, k0)
    otf = fftshift(otf) 

    ckim = otf*imfft
    cim = ifftn(ckim)[0:dims[0], 0:dims[1]]
    cim = cim.real

    if get_psf:
        psf = ifftn(otf)
        psf = fftshift(psf)
        psf = sqrt(psf.real**2 + psf.imag**2)
        psf = pixmodel._centered(psf, dims)
        return cim, psf
    else:
        return cim

def preWhitenCube(**kwargs):
	'''
	Pre-whitenening using noise estimates from a cube taken from the difference map.
	Returns a the pre-whitened volume and various spectra. (Alp Kucukelbir, 2013)

	'''
	print '\n= Pre-whitening the Cubes'
	tStart = time()

	n             = kwargs.get('n', 0)
	vxSize        = kwargs.get('vxSize', 0)
	elbowAngstrom = kwargs.get('elbowAngstrom', 0)
	rampWeight    = kwargs.get('rampWeight',1.0)
	dataF         = kwargs.get('dataF', 0)
	dataBGF       = kwargs.get('dataBGF', 0)
	dataBGSpect   = kwargs.get('dataBGSpect', 0)

	epsilon = 1e-10

	pWfilter = createPreWhiteningFilter(n             = n,
										spectrum      = dataBGSpect,
										elbowAngstrom = elbowAngstrom,
										rampWeight    = rampWeight,
										vxSize        = vxSize)

	# Apply the pre-whitening filter to the inside cube
	dataF       = np.multiply(pWfilter['pWfilter'],dataF)

	dataPWFabs  = np.abs(dataF)
	dataPWFabs  = dataPWFabs-np.min(dataPWFabs)
	dataPWFabs  = dataPWFabs/np.max(dataPWFabs)
	dataPWSpect = sphericalAverage(dataPWFabs**2) + epsilon

	dataPW = np.real(fftpack.ifftn(fftpack.ifftshift(dataF)))
	del dataF

	# Apply the pre-whitening filter to the outside cube
	dataBGF       = np.multiply(pWfilter['pWfilter'],dataBGF)

	dataPWBGFabs  = np.abs(dataBGF)
	dataPWBGFabs  = dataPWBGFabs-np.min(dataPWBGFabs)
	dataPWBGFabs  = dataPWBGFabs/np.max(dataPWBGFabs)
	dataPWBGSpect = sphericalAverage(dataPWBGFabs**2) + epsilon

	dataBGPW = np.real(fftpack.ifftn(fftpack.ifftshift(dataBGF)))
	del dataBGF

	m, s = divmod(time() - tStart, 60)
	print "  :: Time elapsed: %d minutes and %.2f seconds" % (m, s)

	return {'dataPW':dataPW, 'dataBGPW':dataBGPW, 'dataPWSpect': dataPWSpect, 'dataPWBGSpect': dataPWBGSpect, 'peval': pWfilter['peval'], 'pcoef': pWfilter['pcoef'] }

    def Ahfunc(self, f):
        fs = reshape(f, (self.height, self.width, self.depth))

        F = fftn(fs)

        d_1 = ifftshift(ifftn(F*self.Ht));

        d_e = ifftshift(ifftn(F*self.Het));

        d_e2 = ifftshift(ifftn(F*self.He2t));

        d = (1.5*d_1 + 2*real(d_e*exp(1j*self.alpha)) + 0.5*real(d_e2*exp(2*1j*self.alpha)));
        
        d = real(d);
        return ravel(d)

    def Afunc(self, f):
        fs = reshape(f, (self.height, self.width, self.depth))

        F = fftn(fs)

        d_1 = ifftshift(ifftn(F*self.H));
        
        d_e = ifftshift(ifftn(F*self.He));
        
        d_e2 = ifftshift(ifftn(F*self.He2));
        
        d = (1.5*real(d_1) + 2*real(d_e*self.e1) + 0.5*real(d_e2*self.e2))
        
        d = real(d);
        return ravel(d)

def fftconvolve(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1 + s2 - 1

    # Always use 2**n-sized FFT
    fsize = 2 ** np.ceil(np.log2(size))
    IN1 = fftn(in1, fsize)
    IN1 *= fftn(in2, fsize)
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1, axis=0) > product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret, osize)
    elif mode == "valid":
        return _centered(ret, abs(s2 - s1) + 1)

def fast_multinomial(pii, nsum, thresh):
    """generate multinomial distribution for given probability tuple pii.

    *nsum* is the overal number of atoms of a fixed element, pii is a tuple holding the
    distribution of the isotopes.

    this generator yields all combinations and their probabilities which are above *thresh*.

    Remark: the count of the first isotope of all combinations is not computed and "yielded", it is
    automatically *nsum* minus the sum of the elemens in the combinations. We could compute this
    value in this generator, but it is faster to do this later (so only if needed).

    Example:: given three isotopes ([n1]E, [n2]E, [n3]E) of an element E which have
              probabilities 0.2, 0.3 and 0.5.

    To generate all molecules consisting of 5 atoms of this element where the overall probability
    is abore 0.1 we can run:

        for index, pi in gen_n((0.2, 0.3, 0.5), 5, 0.1):
            print(index, pi)

    which prints:

        (1, 3) 0.15
        (2, 2) 0.135
        (2, 3) 0.1125

    the first combination refers to (1, 1, 3) (sum is 5), the second to (1, 2, 2) and the last to
    (0, 2, 3).

    So the probability of an molecule with the overall formula [n1]E1 [n2]E1 [n3]E3 is 0.15, for
    [n1]E1 [n2]E2 [n3]E2 is 0.135, and for [n2]E2 [n3]E3 is 0.1125.

    Implementation:: multinomial distribution can be described as the n times folding (convolution)
    of an underlying simpler distribution. convolution can be fast computed with fft + inverse
    fft as we do below.

    This is often 100 times faster than the old implementatation computing the full distribution
    using its common definition.
    """
    n = len(pii)

    if n == 1:
        yield (0,), pii[0]
        return

    dim = n - 1
    a = np.zeros((nsum + 1,) * dim)
    a[(0,) * dim] = pii[0]
    for i, pi in enumerate(pii[1:]):
        idx = [0] * dim
        idx[i] = 1
        a[tuple(idx)] = pi

    probs = ifftn(fftn(a) ** nsum).real
    mask = probs >= thresh
    pi = probs[mask]
    ii = zip(*np.where(mask))
    for iii, pii in zip(ii, pi):
        yield iii, pii

  def FourierToSpaceTimeNorway(P12,fx1,fy1,f1,**kwargs):

    opt = dotdict({'resolution' : np.r_[6e-5, 6e-5, 6e-5],
                   'c'          : 1540.0})
    opt.update(**kwargs)

    resolution = opt.resolution
    c          = opt.c

    testPlot = False

    dfx=(fx1[1]-fx1[0]);
    dfy=(fy1[1]-fy1[0]);
    df=(f1[1]-f1[0]);

    Nfx=int(2*np.round(c/dfx/resolution[0]/2)+1)
    Nfy=int(2*np.round(c/dfy/resolution[1]/2)+1)
    Nfs=int(2*np.round(c/2/df/resolution[2]/2)+1)

    xax=np.linspace(-1/dfx*c/2,1/dfx*c/2,Nfx);
    yax=np.linspace(-1/dfy*c/2,1/dfy*c/2,Nfy);
    zax=np.linspace(-1/df*c/4,1/df*c/4,Nfs);

    [Nx,Ny,Nf]=P12.shape
    P12zp=np.zeros((Nfx,Nfy,Nfs),dtype=np.complex128)
    ix1=int(np.round(1+(Nfx-1)/2-(Nx-1)/2))
    iy1=int(np.round(1+(Nfy-1)/2-(Ny-1)/2))
    if1=int(np.round(1+(Nfs-1)/2+1+f1[0]/df) - 1)
    P12zp[ix1:ix1+Nx,iy1:iy1+Ny,if1:if1+Nf]=P12;
    P12zp=np.fft.fftshift(P12zp);
    p12=ifftn(P12zp);
    p12=np.fft.fftshift(p12);

    return (p12,xax,yax,zax)

def preparePSF(md, PSSize):
    global PSFFileName, cachedPSF, cachedOTF2, cachedOTFH, autocorr
    
    PSFFilename = md.PSFFile
                
    if (not (PSFFileName == PSFFilename)) or (not (cachedPSF.shape == PSSize)):
        try:
            ps, vox = md.taskQueue.getQueueData(md.dataSourceID, 'PSF')
        except:
            fid = open(getFullExistingFilename(PSFFilename), 'rb')
            ps, vox = pickle.load(fid)
            fid.close()
            
        ps = ps.max(2)
        ps = ps - ps.min()
        #ps = ps*(ps > 0)
        ps = ps*scipy.signal.hanning(ps.shape[0])[:,None]*scipy.signal.hanning(ps.shape[1])[None,:]
        ps = ps/ps.sum()
        PSFFileName = PSFFilename
        pw = (numpy.array(PSSize) - ps.shape)/2.
        pw1 = numpy.floor(pw)
        pw2 = numpy.ceil(pw)
        cachedPSF = pad.with_constant(ps, ((pw2[0], pw1[0]), (pw2[1], pw1[1])), (0,))
        cachedOTFH = ifftn(cachedPSF)*cachedPSF.size
        cachedOTF2 = cachedOTFH*fftn(cachedPSF)

    def recon_dm_trans(self):

        for i, (x_start, x_end, y_start, y_end) in enumerate(self.point_info):
            prb_obj =  self.prb[:,:] * self.obj[x_start:x_end,y_start:y_end]
            tmp = 2. * prb_obj - self.product[i]

            if self.sf_flag:
                tmp_fft = sf.fftn(tmp) / npy.sqrt(npy.size(tmp))
            else:
                tmp_fft = npy.fft.fftn(tmp) / npy.sqrt(npy.size(tmp))
    
            amp_tmp = npy.abs(tmp_fft)
            ph_tmp = tmp_fft / (amp_tmp+self.sigma1)
            (index_x,index_y) = npy.where(self.diff_array[i] >= 0.)
            dev = amp_tmp - self.diff_array[i]
            power = npy.sum(npy.sum((dev[index_x,index_y])**2))/(self.nx_prb*self.ny_prb)
    
            if power > self.sigma2: 
                amp_tmp[index_x,index_y] = self.diff_array[i][index_x,index_y] + dev[index_x,index_y] * npy.sqrt(self.sigma2/power)

            if self.sf_flag:
                tmp2 =  sf.ifftn(amp_tmp*ph_tmp) *  npy.sqrt(npy.size(tmp))
            else:
                tmp2 = npy.fft.ifftn(amp_tmp*ph_tmp) * npy.sqrt(npy.size(tmp))
                    
            self.product[i] += self.beta*(tmp2 - prb_obj)

        del(prb_obj)
        del(tmp)
        del(amp_tmp)
        del(ph_tmp)
        del(tmp2)