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)

    def icwt2d(self, da=0.25):
        '''
        Inverse bi-dimensional continuous wavelet transform as in Wang and
        Lu (2010), equation [5].

        Parameters
        ----------
        da : float, optional
            Spacing in the frequency axis.
        '''
        if self.Wf is None:
            raise TypeError("Run cwt2D before icwt2D")
        m0, l0, k0 = self.Wf.shape

        if m0 != self.scales.size:
            raise Warning('Scale parameter array shape does not match\
                           wavelet transform array shape.')
        # Calculates the zonal and meridional wave numters.
        L, K = 2 ** int(np.ceil(np.log2(l0))), 2 ** int(np.ceil(np.log2(k0)))
        # Calculates the zonal and meridional wave numbers.
        l, k = fftfreq(L, self.dy), fftfreq(K, self.dx)
        # Creates empty inverse wavelet transform array and fills it for every
        # discrete scale using the convolution theorem.
        self.iWf = np.zeros((m0, L, K), 'complex')
        for i, an in enumerate(self.scales):
            psi_ft_bar = an * self.wavelet.psi_ft(an * k, an * l)
            W_ft = fftn(self.Wf[i, :, :], s=(L, K))
            self.iWf[i, :, :] = ifftn(W_ft * psi_ft_bar, s=(L, K)) *\
                da / an ** 2.

        self.iWf = self.iWf[:, :l0, :k0].real.sum(axis=0) / self.wavelet.cpsi

        return self

    def _make_noise_input(self, init):
        """
            Creates an initial input (generated) image.
        """

        # specify dimensions and create grid in Fourier domain
        dims = tuple(self.net.blobs["data"].data.shape[2:]) + \
               (self.net.blobs["data"].data.shape[1], )
        grid = np.mgrid[0:dims[0], 0:dims[1]]

        # create frequency representation for pink noise
        Sf = (grid[0] - (dims[0]-1)/2.0) ** 2 + \
             (grid[1] - (dims[1]-1)/2.0) ** 2
        Sf[np.where(Sf == 0)] = 1
        Sf = np.sqrt(Sf)
        Sf = np.dstack((Sf**int(init),)*dims[2])

        # apply ifft to create pink noise and normalize
        ifft_kernel = np.cos(2*np.pi*np.random.randn(*dims)) + \
                      1j*np.sin(2*np.pi*np.random.randn(*dims))
        img_noise = np.abs(ifftn(Sf * ifft_kernel))
        img_noise -= img_noise.min()
        img_noise /= img_noise.max()

        # preprocess the pink noise image
        x0 = self.transformer.preprocess("data", img_noise)

        return x0

def fftconvolve(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.
    
    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    if (s1.dtype.char in ['D','F']) or (s2.dtype.char in ['D', 'F']):
        cmplx=1
    else: cmplx=0
    size = s1+s2-1
    IN1 = fftn(in1,size)
    IN1 *= fftn(in2,size)
    ret = ifftn(IN1)
    del IN1
    if not cmplx:
        ret = real(ret)
    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 apply(self, image, freq_image=False, invert=True):
        """
        apply(image, freq_image=False, invert=True)

        Apply the current filter.

        Parameters
        ----------
        image : array_like
            Image to apply filter to.

        freq_image : bool optional
            Flag to indicate if the input image is already the DFT of the
            input image.

        invert : bool optional
            Flag to indicate if the output image should be inverted from the
            frequency domain.
        """
        try:
            H = self.H
        except AttributeError:
            raise AttributeError, 'No filter currently set.'

        if not freq_image:
            image = fftpack.fftn(image)

        E = image * H

        if invert:
            E = np.real(fftpack.ifftn(E))

        return E

def rolloff(grid_shape, nhood_shape=[6,6], ncrop=None, osfactor=2, threshold=0.01, axes=[-1, -2], combi=False, slice_prof_coef=1):
    if combi:
        nz = spatial_dims[-1]
        spatial_dims = grid_shape[:-1]
        ndims = len(spatial_dims)
        ro_data = ones(nz, dtype='complex64')
        ro_K = zeros((nz, ndims), dtype='float32')
        ro_K[:,2] = linspace(-nz/2, nz/2, nz, endpoint=False)
    else:
        spatial_dims = grid_shape
        ndims = len(spatial_dims)
        ro_data = array([1.0], dtype='complex64')
        ro_K = array([[0]*ndims], dtype='float32')
    G = sparse_matrix_operator(ro_K, grid_shape=spatial_dims, nhood_shape=nhood_shape, osfactor=osfactor, combi=combi, slice_prof_coef=slice_prof_coef)
    ro = G.T * ro_data
    n = sqrt(G.shape[1])
    ro = reshape(ro, (n,n))
    #if osfactor > 1 and ncrop == None:
    #    ncrop = int((spatial_dims[0]/osfactor) * (osfactor - 1) / 2)
    ro = fftshift(abs(ifftn(fftshift(ro, axes=axes), axes=axes)), axes=axes)  # transform to image
    if ncrop > 0:
        ro = ro[ncrop:-ncrop, ncrop:-ncrop]
    #print 'rolloff shape:', ro.shape
    ro = ro / ro.max()  # normalize
    ro[ro < threshold] = 1.0
    #ro_max = ro.max()
    #ro[ro < threshold*ro_max] = 1.0
    ro = 1.0 / ro
    #ro = ro**2
    #print 'TOOK OUT SQUARED RO'
    #ro = ro / ro.max()
    return ro