def colorize_by_power(self, image):
     """
     Colorize the image mode-by-mode according to the power in each mode.  
     The top third of modes are colored red, the middle third green, and 
     the lower third blue.  For RGB images, a grayscale equivalent is 
     computed and colorized. 
     """
     print "colorizing....."
     if len(image.shape) == 3:
         power = fft2(np.sum(image, axis=2))**2
     elif len(image.shape) == 2:
         power = fft2(image)**2
     else:
         raise Exception("Invalid image shape: {}".foramt(image.shape))
     thirds = (power.max() - power.min())/3.0
     third_cut = power.min() + thirds
     twothird_cut = third_cut + thirds
     lower = power < third_cut
     upper = power > twothird_cut
     middle = ~(lower | upper)
     colorized = np.zeros((power.shape[0], power.shape[1], 3), 
                          dtype=np.uint8)
     for color, region in enumerate([upper, middle, lower]):
         new_channel = ifft2(np.where(region, power, 0.0))
         shifted = (new_channel - new_channel.min())
         scaled = 255.0*shifted/shifted.max()
         colorized[..., color] = ifft2(np.where(region, power, 0.0))
     return colorized

def correlate_layer(pattern_layer, source_layer):
    """
    Normalized Cross-Correlation for a single channel of an RGB image
    (or a greyscale image). Normalization is done as follows:
    normalized = (x - mean(x)) / std(x)

    pattern_layer - Two-dimensional ndarray, single channel of pattern image
    source_layer - Two-dimensional ndarray, single channel of source image

    """
    # http://bit.ly/WsRveH
    if pattern_layer.std() == 0:
        normalized_pattern = pattern_layer
    else:
        normalized_pattern = ((pattern_layer - np.mean(pattern_layer)) /
                              (np.std(pattern_layer) * pattern_layer.size))
    if source_layer.std() == 0:
        normalized_source = source_layer
    else:
        normalized_source = ((source_layer - np.mean(source_layer)) /
                             np.std(source_layer))

    #Take the fft of both Images, padding the pattern out with 0's
    # to be the same shape as the source
    pattern_fft = fftpack.fft2(normalized_pattern, source_layer.shape)
    source_fft = fftpack.fft2(normalized_source)

    # Perform the correlation in the frequency domain, which just the
    # inverse FFT of the pattern matrix's conjugate * the source matrix
    # http://en.wikipedia.org/wiki/Cross-correlation#Properties
    return fftpack.ifft2(pattern_fft.conjugate() * source_fft)

    def shift_inner(arr, nx, ny, window=False, padding='reflect'):
        """
        Shifts an array by nx and ny respectively.

        """

        if ((nx % 1. == 0.) and (ny % 1. ==0)):
            return sp.roll(sp.roll(arr, int(ny), axis=0),
                           int(nx), axis=1)
        else:
            atype = arr.dtype
            if padding:
                x, y = arr.shape
                pwx, pwy = int(pow(2., np.ceil(np.log2(1.5*arr.shape[0])))), int(pow(2., np.ceil(np.log2(1.5*arr.shape[1]))))
                pwx2, pwy2 = (pwx-x)/2, (pwy-y)/2
                if pad=='zero':
                    arr = pad.with_constant(arr, pad_width=((pwx2, pwx2), (pwy2, pwy2)))
                else:
                    arr = pad.with_reflect(arr, pad_width=((pwx2, pwx2), (pwy2, pwy2)))
            phaseFactor = sp.exp(complex(0., -2.*sp.pi)*(ny*spf.fftfreq(arr.shape[0])[:, np.newaxis]+nx*spf.fftfreq(arr.shape[1])[np.newaxis, :]))
            if window:
                window = spf.fftshift(CXData._tukeywin(arr.shape[0], alpha=0.35))
                arr = spf.ifft2(spf.fft2(arr)*phaseFactor*window)
            else:
                arr = spf.ifft2(spf.fft2(arr)*phaseFactor)
            if padding:
                arr = arr[pwx/4:3*pwx/4, pwy/4:3*pwy/4]

        if atype == 'complex':
            return arr
        else:
            return np.real(arr)

def convolution_fourier_RGB(img, fil_fft, fftsize):

    channelR = np.zeros((img.shape[0],img.shape[1]), 'double')
    channelG = np.zeros((img.shape[0],img.shape[1]), 'double')
    channelB = np.zeros((img.shape[0],img.shape[1]), 'double')

    for x in range(img.shape[0]):
        for y in range(img.shape[1]):
            channelR[x,y] = img[x,y][0]
            channelG[x,y] = img[x,y][1]
            channelB[x,y] = img[x,y][2]
        
    matrixR_fft = fftpack.fft2(channelR, (fftsize, fftsize))
    matrixG_fft = fftpack.fft2(channelG, (fftsize, fftsize))
    matrixB_fft = fftpack.fft2(channelB, (fftsize, fftsize))
    
    matrixR_fil_fft = matrixR_fft * fil_fft;
    matrixG_fil_fft = matrixG_fft * fil_fft;
    matrixB_fil_fft = matrixB_fft * fil_fft;

    matrixR_fil = np.real(fftpack.ifft2(matrixR_fil_fft))
    matrixG_fil = np.real(fftpack.ifft2(matrixG_fil_fft))
    matrixB_fil = np.real(fftpack.ifft2(matrixB_fil_fft))
    
    img_fil = np.zeros((matrixR_fil.shape[0], matrixR_fil.shape[1], 3), 'double')

    for x in range(matrixR_fil.shape[0]):
        for y in range(matrixR_fil.shape[1]):
            img_fil[x,y,0] = matrixR_fil[x,y]
            img_fil[x,y,1] = matrixG_fil[x,y]
            img_fil[x,y,2] = matrixB_fil[x,y]
            
    return img_fil

    def _storeLensKernel(self):
        '''
        This is an internal function that will store the lensing kernel for the given stamp dimensions
        when the simulator is initialized.
        '''

        N=self.lensN
        scale=self.lensScale

        i=np.array(range(0,N))
        j=np.array(range(0,N))


        xo=(0.5+i-(N)/2.0)  # pixels!
        yo=(0.5+j-(N)/2.0)
        xMesh, yMesh = np.meshgrid(xo,yo)

        r2Mesh = (xMesh**2.+yMesh**2.)

        kX=(1./pi)*xMesh/r2Mesh  # pixels!
        kY=(1./pi)*yMesh/r2Mesh

        
        self.ftKernelX=ff.fft2(kX)
        self.ftKernelY=ff.fft2(kY)
        del(xMesh)
        del(yMesh)
        del(r2Mesh)
        del(kX)
        del(kY)

def filters_bank(M, N, J, L=8):
    filters = {}
    filters['psi'] = []

    offset_unpad = 0
    for j in range(J):
        for theta in range(L):
            psi = {}
            psi['j'] = j
            psi['theta'] = theta
            psi_signal = morlet_2d(M, N, 0.8 * 2**j, (int(L - L / 2 - 1) - theta) * np.pi / L, 3.0 / 4.0 * np.pi / 2**j,offset=offset_unpad)  # The 5 is here just to match the LUA implementation :)
            psi_signal_fourier = fft.fft2(psi_signal)
            for res in range(j + 1):
                psi_signal_fourier_res = crop_freq(psi_signal_fourier, res)
                psi[res] = tf.constant(np.stack((np.real(psi_signal_fourier_res), np.imag(psi_signal_fourier_res)), axis=2))
                psi[res] = tf.div(psi[res], (M * N // 2**(2 * j)), name="psi_theta%s_j%s" % (theta, j))
            filters['psi'].append(psi)

    filters['phi'] = {}
    phi_signal = gabor_2d(M, N, 0.8 * 2**(J - 1), 0, 0, offset=offset_unpad)
    phi_signal_fourier = fft.fft2(phi_signal)
    filters['phi']['j'] = J
    for res in range(J):
        phi_signal_fourier_res = crop_freq(phi_signal_fourier, res)
        filters['phi'][res] = tf.constant(np.stack((np.real(phi_signal_fourier_res), np.imag(phi_signal_fourier_res)), axis=2))
        filters['phi'][res] = tf.div(filters['phi'][res], (M * N // 2 ** (2 * J)), name="phi_res%s" % res)

    return filters

def filters_bank(M, N, J, L=8):
    filters = {}
    filters['psi'] = []

    offset_unpad = 0
    for j in range(J):
        for theta in range(L):
            psi = {}
            psi['j'] = j
            psi['theta'] = theta
            psi_signal = morlet_2d(M, N, 0.8 * 2**j, (int(L - L / 2 - 1) - theta) * np.pi / L, 3.0 / 4.0 * np.pi / 2**j,offset=offset_unpad)  # The 5 is here just to match the LUA implementation :)
            psi_signal_fourier = fft.fft2(psi_signal)
            for res in range(j + 1):
                psi_signal_fourier_res = crop_freq(psi_signal_fourier, res)
                psi[res] = torch.FloatTensor(np.stack((np.real(psi_signal_fourier_res), np.imag(psi_signal_fourier_res)), axis=2))
                # Normalization to avoid doing it with the FFT!
                psi[res].div_(M * N // 2**(2 * j))
            filters['psi'].append(psi)

    filters['phi'] = {}
    phi_signal = gabor_2d(M, N, 0.8 * 2**(J - 1), 0, 0, offset=offset_unpad)
    phi_signal_fourier = fft.fft2(phi_signal)
    filters['phi']['j'] = J
    for res in range(J):
        phi_signal_fourier_res = crop_freq(phi_signal_fourier, res)
        filters['phi'][res] = torch.FloatTensor(np.stack((np.real(phi_signal_fourier_res), np.imag(phi_signal_fourier_res)), axis=2))
        filters['phi'][res].div_(M * N // 2 ** (2 * J))

    return filters

def drawPic():
    m = 200
    n = 200
    a = 20
    b = 20
    r = 1
    vp0 = 1500
    variance = 0.08
    dx = 3
    C = np.zeros((m, n))

    au = autocorrelation_func(a, b)

    for i in xrange(m):
        for j in xrange(n):
            C[i, j] = au.gaussian((i-(m-1)/2.)*dx, (j-(n-1)/2.)*dx)
    S = fft2(C)  # power spectral density
    z = np.random.randn(m, n)
    Z = fft2(z)

    G = np.sqrt(dx * S)

    GZ = G * Z
    gz = ifft2(GZ)

    A = np.real(gz)

    K = np.real(A)
    f = plt.figure()
    plt.imshow(K)
    return f

def frankotchellappa(dzdx, dzdy):
    rows, cols = dzdx.shape

    # The following sets up matrices specifying frequencies in the x and y
    # directions corresponding to the Fourier transforms of the gradient
    # data.  They range from -0.5 cycles/pixel to + 0.5 cycles/pixel.
    # The scaling of this is irrelevant as long as it represents a full
    # circle domain. This is functionally equivalent to any constant * pi
    pi_over_2 = np.pi / 2.0
    row_grid = np.linspace(-pi_over_2, pi_over_2, rows)
    col_grid = np.linspace(-pi_over_2, pi_over_2, cols)
    wy, wx = np.meshgrid(row_grid, col_grid, indexing='ij')

    # Quadrant shift to put zero frequency at the appropriate edge
    wx = ifftshift(wx)
    wy = ifftshift(wy)

    # Fourier transforms of gradients
    DZDX = fft2(dzdx)
    DZDY = fft2(dzdy)

    # Integrate in the frequency domain by phase shifting by pi/2 and
    # weighting the Fourier coefficients by their frequencies in x and y and
    # then dividing by the squared frequency
    denom = (wx ** 2 + wy ** 2)
    Z = (-1j * wx * DZDX - 1j * wy * DZDY) / denom
    Z = np.nan_to_num(Z)

    return np.real(ifft2(Z))

 def getGeneralStatistics(self, hara=False, zern=False, tamura=False, only1D=None):
     generalStatistics = []
     if self.rows == 1 and self.columns == 1:
         for index in range(3):
             generalStatistics.append(self.image[0, 0, index])
         return generalStatistics
     if not only1D is None:
         im = only1D
         generalStatistics.extend(self._calculateStatistics(im, haralick=hara, zernike=zern))
         fourierTransform = np.abs(fftpack.fft2(im))  # fourierTransform
         generalStatistics.extend(self._calculateStatistics(fourierTransform))
         waveletTransform = pywt.dwt2(im, "sym5")[0]
         generalStatistics.extend(self._calculateStatistics(waveletTransform))
         waveletFourierTransform = pywt.dwt2(fourierTransform, "sym5")[0]
         generalStatistics.extend(self._calculateStatistics(waveletFourierTransform))
         if tamura:
             generalStatistics.extend(self.get3Dstatistics(tamura=True))
         return generalStatistics
     for index in range(3):
         im = self.image[:, :, index]
         generalStatistics.extend(self._calculateStatistics(im, haralick=hara, zernike=zern))
         fourierTransform = np.abs(fftpack.fft2(im))  # fourierTransform
         generalStatistics.extend(self._calculateStatistics(fourierTransform))
         waveletTransform = pywt.dwt2(im, "sym5")[0]
         generalStatistics.extend(self._calculateStatistics(waveletTransform))
         waveletFourierTransform = pywt.dwt2(fourierTransform, "sym5")[0]
         generalStatistics.extend(self._calculateStatistics(waveletFourierTransform))
     if tamura:
         generalStatistics.extend(self.get3Dstatistics(tamura=True))
     return generalStatistics

def phase_cor( A, B ):
  """
  find the 2D correlation between images A and B. This is much faster than any
  of the other scipy correlation methods which insist on working in the spatial
  domain.
  """
  return ( ifft2( fft2(A)*numpy.conj(fft2(B)) ) ).real

def inverseFilter(img,fftsize):
    im = np.mean(img,2)/255.
    
    fftsize = 1024
    im_fft = fftpack.fft2(im, (fftsize, fftsize))
    
    #Complementary of a Gaussian filter
    SZ = 1024
    sigma = 0.25
    [xx,yy]=np.meshgrid(np.linspace(-4,4,SZ),np.linspace(-4,4,SZ))
    gaussian = np.exp(-0.5*(xx*xx+yy*yy)/(sigma*sigma))
    fil =1.-fftpack.fftshift(gaussian/np.max(gaussian))
    
    fil_fft =  fil
    
    im_fil_fft = im_fft * fil_fft
    
    im_fil = np.real(fftpack.ifft2(im_fil_fft))
    
    hs=np.floor(SZ/2.)
    #Careful with the crop. Because we work directly in the Fourier domain there is no padding.
    im_crop = im_fil[0:im.shape[0], 0:im.shape[1]]     
    F=fftpack.fft2(im_crop,(1024,1024))
    H=fil_fft
    tol= 1e-2
    I = F/H
    print np.min(I)
    I=np.where(np.abs(H)<tol,0,I)
    i_reconstructed = np.real(fftpack.ifft2(I))
    plt.imshow(i_reconstructed[:im.shape[0],:im.shape[1]],cmap="gray")

    def solve(self, u, v, dx, dy):
        import numexpr as ne
        nx, ny = u.shape
        assert u.shape == tuple(self.shape)

        fu = fft2(u)
        fv = fft2(v)

        mpx = self.mpx
        mmx = self.mmx
        dpx = self.dpx
        dmx = self.dmx
        mpy = self.mpy
        mmy = self.mmy
        dpy = self.dpy
        dmy = self.dmy

        d = ne.evaluate("fu*mmy * dmx + fv * mmx * dmy")
        lapl = ne.evaluate("mpy  * mmy * dpx * dmx + mpx*mmx *dpy *dmy")
        lapl[0, 0] = 1.0

        p = d / lapl
        px = np.real(ifft2(mpy * dpx * p))
        py = np.real(ifft2(mpx * dpy * p))

        # self.p = np.real(ifft2(p))

        u -= px
        v -= py

        return px, py

    def _getCrossCorrelation(self, ref, mask, fft_ref = False, fft_mask = False):
        """ Computes the cross correlation between reference and mask images.
        For parameter description, refer to <self._getDriftValue()> """
        # Images should be square and of same dimensions at this point.
        assert(ref.shape==mask.shape)

        if not fft_ref:
            four_ref = fft2(ref)
        else:
            four_ref = ref

        if not fft_mask:
            # Crop the mask and replace the edges with 0.0 values. Helps the crosscorrelation.
            if self.cropping:
                size = min(mask.shape)
                crop = self.cropping
                mask_cropped = np.copy(mask[crop:(size-crop), crop:(size-crop)])
                mask_padded = np.pad(mask_cropped, crop, mode='constant')
                four_mask = fft2(mask_padded)
            else:
                four_mask = fft2(mask)           
        else:
            four_mask = mask

        # Conjugate the mask.
        four_mask_conj = np.conjugate(four_mask)
        # Compute pointwise product of reference and mask.
        product = np.multiply(four_mask_conj, four_ref)
        # Compute ifft of this product
        xcorr = ifft2(product)
        # Take the absolute value
        xcorr_abs = np.absolute(xcorr)
        return xcorr_abs

 def step4(self):
   '''
   Perform a 4th order timestep
   '''
   
   def order2(c):
     Vc = np.exp( -1j * c * self.dt / 2. * 
                   ( self.V - self.gravity() + 
                     self.g * abs( self.psi ) ** 2
                   )
                 )
     Tc = self.expksquare ** c
     return Vc, Tc
   
   p = 1/(4.-4.**(1/3.))
   q = 1 - 4 * p
   
   Vp,Tp = order2( p )
   Vq,Tq = order2( q )
   
   return Vp * ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp ** 2 * 
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp * Vq *
               ff.fftshift( ff.ifft2( Tq * ff.fft2( ff.fftshift( Vq * Vp * 
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp ** 2 *  
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp * self.psi 
               ) ) ) )
               ) ) ) )
               ) ) ) )
               ) ) ) )
               ) ) ) )

def convolve2D_fft(theta,F):
    func = KS_kernel_real
    #func = gaussian_real

    assert theta.shape == F.shape

    N1,N2 = theta.shape
    
    dtheta1 = theta[1,1].real - theta[0,0].real
    theta1 = dtheta1*(numpy.arange(2*N1)-N1)
    theta1 = fftpack.ifftshift(theta1)
    
    dtheta2 = theta[1,1].imag - theta[0,0].imag
    theta2 = dtheta2*(numpy.arange(2*N2)-N2)
    theta2 = fftpack.ifftshift(theta2)
    
    theta_kernel = numpy.zeros((2*N1,2*N2),dtype=complex)
    theta_kernel += theta1.reshape((2*N1,1))
    theta_kernel += 1j*theta2

    kernel = func(theta_kernel)
    
    dA = dtheta1 * dtheta2

    F_fft = fftpack.fft2(F, (2*N1,2*N2) ) * dA
    F_fft *= fftpack.fft2(kernel,(2*N1,2*N2) ) 
    
    out = fftpack.ifft2(F_fft)
    
    return out[:N1,:N2]

def mvd_wiener(initImg, imgList, psfList, iterNum, mu, positiveOnly=True):
    if positiveOnly:
        initImg[initImg < 0.0] = 0.0
    viewNum = len(imgList)
    fftFactor = np.sqrt(initImg.shape[0]*initImg.shape[1])
    mu = mu * fftFactor
    I = np.sum(np.abs(initImg))
    e = fft2(initImg)
    e_img_old = initImg
    e_img = initImg
    if iterNum == 0:
        return e_img
    # pre-compute spectra
    ijList = [fft2(img) for img in imgList]
    pjList = [fft2(pad_and_center_psf(psf, initImg.shape)) for psf in psfList]
    for i in xrange(iterNum):
        c_all = np.zeros(e.shape, dtype=float)
        for j in xrange(viewNum):
            ij = ijList[j]
            pj = pjList[j]
            sj = e * pj
            cj = (np.conj(pj) * (ij - sj))/(np.square(np.abs(pj)) + mu**2)
            c_all = c_all + cj / float(viewNum)
        e = e + c_all
        e_img = np.real(ifft2(e))
        if positiveOnly:
            e_img[e_img < 0.0] = 0.0
        e_img = e_img / np.sum(np.abs(e_img)) * I
        e = fft2(e_img)
        print 'iter #%d, total change: %f.' %\
            (i+1, np.sum(np.abs(e_img_old-e_img))/I)
        e_img_old = e_img
    return e_img

def tsvd_fft(B, PSF, center=None, tol=None,i=None):
    """TSVD_FFT Truncated SVD image deblurring using the FFT algorithm.

    X, tol = tsvd_fft(B, PSF, center, tol)

  Compute restoration using an FFT-based truncated spectral factorization.

  Input:
        B  Array containing blurred image.
      PSF  Array containing the point spread function.
  
  Optional Inputs:
   center  [row, col] = indices of center of PSF.
      tol  Regularization parameter (truncation tolerance).
             Default parameter chosen by generalized cross validation.

  Output:
        X  Array containing computed restoration.
      tol  Regularization parameter used to construct restoration."""
    #
    # compute the center of the PSF if it is not provided
    #
    if center is None:
        center = array(PSF.shape) / 2

    #
    # if PSF is smaller than B, pad it to the same size as B
    #
    if PSF.size < B.size:
        PSF = padarray(PSF, array(B.shape) - array(PSF.shape),
                       direction='post')
    
    #
    # Use the FFT to compute the eigenvalues of the BCCB blurring matrix.
    #
    S = fft2( circshift(PSF, 0-center) )
    #
    # If a regularization parameter is not given, use GCV to find one.
    #
    bhat = fft2(B)
    if i !=None:
        ev=abs(S.flatten())
        ev.sort()
        ev=ev[::-1]
        tol=ev[i] 
    elif tol is None:
        tol = gcv_tsvd(S.flatten('f'), bhat.flatten('f'))
	

    #
    # Compute the TSVD regularized solution.
    #
    
    idx = where(abs(S) >= tol)
    Sfilt = zeros(shape(bhat),'d')
    Sfilt[idx] = 1 / S[idx]
    X = real(ifft2(bhat * Sfilt))

    return X, sorted(S.flatten())[::-1]

def tik_fft(B, PSF, center=None, alpha=None,sigma=None):
    """TIK_FFT Tikhonov image deblurring using the FFT algorithm.

    X, alpha = tik_fft(B, PSF, center, alpha)

    Compute restoration using an FFT-based Tikhonov filter, 
    with the identity matrix as the regularization operator.

    Input:
    B  Array containing blurred image.
    PSF  Array containing the point spread function.
  
    Optional Inputs:
    center  [row, col] = indices of center of PSF.
    alpha  Regularization parameter.
    Default parameter chosen by generalized cross validation.
    
    Output:
    X  Array containing computed restoration.
    alpha  Regularization parameter used to construct restoration.
    """
    #
    # compute the center of the PSF if it is not provided
    #
    if center is None:
        center = array(PSF.shape) / 2
    #
    # if PSF is smaller than B, pad it to the same size as B
    #
    if PSF.size < B.size:
        PSF = padarray(PSF, array(B.shape) - array(PSF.shape),
                       direction='post')
    #
    # Use the FFT to compute the eigenvalues of the BCCB blurring matrix.
    #
    S = fft2( circshift(PSF, 0-center) )
    s = S.flatten('f')
    bhat = fft2(B)
    bhat = bhat.flatten('f')
    #
    # If a regularization parameter is not given, use GCV to find one.
    #
##    if alpha is None:
    def tik_fft_obj(alpha,B,PSF,sigma):
        return abs(norm(tik_fft(B, PSF,alpha=abs(alpha)),2)-sigma) 
    
    if alpha is None:
        alpha = gcv_tik(s, bhat)
        print "alpha:", alpha,log10(alpha)

    #
    # Compute the Tikhonov regularized solution.
    #
    D = conj(s)*s + abs(alpha**2)
    bhat = conj(s) * bhat
    xhat = bhat / D
    xhat = reshape(asmatrix(xhat), shape(B), 'f')
    X = real(ifft2(xhat))
    return X

def gauss_smooth(arr, w):
    """
    Smooths an input array by w pixels
    """
    alpha = w/2.3548
    pix = arr.shape[0]
    ker = sp.exp(-(alpha**-2.)*sp.hypot(*sp.ogrid[-pix/2:pix/2,-pix/2:pix/2])**2.)
    return sp.real(spf.fftshift(spf.ifft2( spf.fft2(arr)*spf.fft2(ker) )))