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 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 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 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 invert_vort(uc, dx=dx, dy=dy, nx=nx, ny=ny, geom=tad.geom):

    ucv = geom.validview(uc)
    vort = ucv[0]
    u = ucv[1]
    v = ucv[2]

    f = fft2(vort)

    nx, ny = vort.shape

    scal_y = 2*pi/dy/ny
    scal_x = 2*pi/dx/nx

    k = fftfreq(nx, 1/nx)[:,None] * 1j * scal_x
    l = fftfreq(ny, 1/ny)[None,:] * 1j * scal_y


    lapl = k**2 + l**2
    lapl[0,0] = 1.0

    psi = f/lapl
    u[:] = -real(ifft2(psi * l))
    v[:] = real(ifft2(psi * k))

    return uc

 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 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 spec2grid(sfield):
    """
    Transform one frame of SQG model
    output to a grided (physical) representation.  Assumes 'sfield'
    to be up-half plane, and specifies lower half plane by conjugate
    sym (since physical field is assumed real-valued).  Input field
    should have dimensions  (...,kmax+1,2*kmax+1,nz), where
    kmax=2^n-1, hence physical resolution will be 2^(n+1) x 2^(n+1).
    NOTE: top row of the input field corresponds to ky = 0, the
    kx<0 part is NOT assumed a priori to be conjugate- symmetric
    with the kx>0 part.  NOTE: grid2spec(spec2grid(fk)) = fk.
    OPTIONAL: da = true pads input with 0s before transfoming to
    gridspace, for dealiased products.  Default is da = false.
    
    Args:
        sfield: complex spectrum field with shape (t(optional), ky, kx, z(optional))
    """        
    if not _is_single_layer(sfield):
        hres = sfield.shape[-2] + 1
        fk = fullspec(sfield)
        fk = fftpack.ifftshift(fk, axes=(-2,-3))
        return hres*hres*np.real(fftpack.ifft2(fk, axes=(-2,-3)))
    else:
        hres = sfield.shape[-1] + 1
        fk = fullspec(sfield, True)
        fk = fftpack.ifftshift(fk, axes=(-1,-2))
        return hres*hres*np.real(fftpack.ifft2(fk, axes=(-1,-2)))

def find_foci(evt, type,key,type2,key2,minPhase=-500000, maxPhase=500000, steps=101, field_of_view_rad=100, wavelength=1.053, CCD_S_DIST=0.375, PX_SIZE=75e-6):
    img = evt[type][key].data
    centroids = evt[type2][key2].data

    Nfoci = centroids.shape[0]
    Xrange, Yrange = img.shape
    Npixel = field_of_view_rad
    
    p = numpy.linspace(-Xrange/2, Xrange/2-1, Xrange)
    q = numpy.linspace(-Yrange/2, Yrange/2-1, Yrange)
    pp, qq = numpy.meshgrid(p, q)
   
    phase_matrix = (2*numpy.pi/wavelength)*numpy.sqrt(1-((PX_SIZE/CCD_S_DIST)**2)*(qq**2 + pp**2))
    prop_length = numpy.linspace(minPhase, maxPhase, steps)
    
    variance = numpy.zeros([steps, Nfoci])
    # shift stuff for performance reasons
    img_shifted = fftshift(img)
    phase_matrix_shifted = fftshift(phase_matrix)
    
    for idx, phase in enumerate(prop_length):
        
        img_propagated = img_shifted * numpy.exp(1.j*phase*phase_matrix_shifted)
        recon = fftshift(ifft2(img_propagated))
        
        for CC in numpy.arange(Nfoci):
            centerx, centery = centroids[CC, :]
            ###print  centerx, centery
            reconcut = numpy.abs(recon[numpy.max([0, centerx-Npixel-1]).astype(int): numpy.min([Xrange-1, centerx+Npixel]).astype(int), numpy.max([0, centery-Npixel-1]).astype(int): numpy.min([Yrange-1, centery+Npixel]).astype(int)])
            variance[idx, CC] = reconcut.var()
    
    focus_distance = numpy.zeros(Nfoci)
    CC_size = numpy.zeros(Nfoci)
    focused_CC = numpy.zeros(4*Npixel**2 * Nfoci).reshape(Nfoci, 2*Npixel, 2*Npixel)
    
    for CC in numpy.arange(Nfoci):
        ind_max = numpy.argmax(variance[:, CC])
        tmp = variance[:, CC]
        # get max which is not at border
        loc_max_bool = numpy.r_[True, tmp[1:] > tmp[:-1]] & numpy.r_[tmp[:-1] > tmp[1:], True]
        loc_max_bool[0] = False
        loc_max_bool[-1] = False
        ind_max = numpy.argmax(tmp*loc_max_bool)
        
        focus_distance[CC] = prop_length[ind_max]
        img_propagated = img_shifted * numpy.exp(1.j * focus_distance[CC]  * phase_matrix_shifted)
        recon = fftshift(ifft2(img_propagated))
        
        centerx, centery = centroids[CC, :]
        
        reconcut = numpy.real(recon[numpy.max([0, centerx-Npixel]).astype(int): numpy.min([Xrange-1, centerx+Npixel]).astype(int), numpy.max([0, centery-Npixel]).astype(int): numpy.min([Yrange-1, centery+Npixel]).astype(int)])
        focused_CC[CC, 0:reconcut.shape[0], 0:reconcut.shape[1]] = reconcut
        CC_size[CC] = numpy.sum(get_CC_size(reconcut))
    
    if len(focused_CC):
        add_record(evt["analysis"], "analysis", "focused_CC", focused_CC[0])
        add_record(evt["analysis"], "analysis", "focus distance", focus_distance)
        add_record(evt["analysis"], "analysis", "CC_size", CC_size)
        add_record(evt["analysis"], "analysis", "propagation length", prop_length)

def calc_gradients(inputvar, kxgrid, kygrid):
    var_k = fft2(np.complex128(inputvar))
    gradx_k = 1j * kxgrid * var_k
    grady_k = 1j * kygrid * var_k
    dvar_dx = np.real(ifft2(gradx_k))
    dvar_dy = np.real(ifft2(grady_k))
    grad2 = dvar_dx**2 + dvar_dy**2
    return dvar_dx, dvar_dy, grad2 

def ifft(a, overwrite=False, shift=True):
    if shift:
        res = fftpack.ifft2(fftpack.fftshift(a, axes=[0,1]), axes=[0, 1],
                            overwrite_x=overwrite)
    else:
        res = fftpack.ifft2(a, overwrite_x=overwrite)

    return res

 def calcAcovf2d(self):
     """Calculate the 2d auto covariance function. """
     # See Wiener-Kinchine theorem
     if self.shift:
         # Note, the ACovF needs the unshifted 2d PSD for inverse FFT, so unshift.
         #  Then shift back again. 
         self.acovf = fftpack.fftshift(fftpack.ifft2(fftpack.ifftshift(self.psd2d)))
     else:
         self.acovf = fftpack.ifft2(self.psd2d)
     return

def f(vort, t, k=k, l=l, lapl=lapl):
    fv = fft2(vort)
    psi = fv/lapl

    u = ifft2(-psi * l)
    v = ifft2(psi * k)

    adv = -(u* ifft2(fv*k) + v*ifft2(fv*l)) #+ ifft2(lapl * fv/R)

    return adv

	def FilterElectrons(self,sign, Psi):
 		'''
		Routine that uses the Fourier transform to filter positrons/electrons
		Options:
			sign=1   Leaves electrons
			sign=-1	 Leaves positrons
		'''
		print '  '
		print '  	Filter Electron routine '
		print '  '

		px = self.c*self.Px
		py = self.c*self.Py
		
		m = self.mass
		c= self.c

		energy = np.sqrt(  (m*c**2)**2 + px**2 + py**2 )

		EP_11 = 1. + sign*m*c**2/energy
		EP_12 = 0.
		EP_13 = 0.
		EP_14 = sign*(px - 1j*py)/energy
		
		EP_21 = 0.
		EP_22 = 1. + sign*m*c**2/energy
		EP_23 = sign*(px + 1j*py)/energy
		EP_24  = 0.

		EP_31 = 0.
		EP_32 = sign*(px - 1j*py)/energy
		EP_33 = 1. - sign*m*c**2/energy
		EP_34 = 0.

		EP_41 = sign*(px + 1j*py)/energy
		EP_42 = 0.
		EP_43 = 0.
		EP_44 = 1. - sign*m*c**2/energy	
		
		#Psi1, Psi2, Psi3, Psi4 = Psi

		psi1_fft = fftpack.fft2( Psi[0]  ) 
		psi2_fft = fftpack.fft2( Psi[1]  ) 
		psi3_fft = fftpack.fft2( Psi[2]  ) 
		psi4_fft = fftpack.fft2( Psi[3]  ) 		
		
		psi1_fft_electron = EP_11*psi1_fft + EP_12*psi2_fft + EP_13*psi3_fft + EP_14*psi4_fft
		psi2_fft_electron = EP_21*psi1_fft + EP_22*psi2_fft + EP_23*psi3_fft + EP_24*psi4_fft		
		psi3_fft_electron = EP_31*psi1_fft + EP_32*psi2_fft + EP_33*psi3_fft + EP_34*psi4_fft
		psi4_fft_electron = EP_41*psi1_fft + EP_42*psi2_fft + EP_43*psi3_fft + EP_44*psi4_fft
						
		return np.array([ fftpack.ifft2( psi1_fft_electron   ),
				  fftpack.ifft2( psi2_fft_electron   ),
	   			  fftpack.ifft2( psi3_fft_electron   ),
				  fftpack.ifft2( psi4_fft_electron   )  ])

	def _FilterElectrons(self,sign):
 		'''
		Routine that uses the Fourier transform to filter positrons/electrons
		Options:
			sign=1   Leaves electrons
			sign=-1	 Leaves positrons
		'''
		print '  '
		print '  	Filter Electron routine '
		print '  '
                min_Px = np.pi*self.X_gridDIM/(2*self.min_X)
		dPx = 2*np.abs(min_Px)/self.X_gridDIM
		px_Vector  = fftpack.fftshift ( np.linspace(min_Px, np.abs(min_Px) - dPx, self.X_gridDIM ))

		min_Py = np.pi*self.Y_gridDIM/(2*self.min_Y)
		dPy = 2*np.abs(min_Py)/self.Y_gridDIM
		py_Vector  = fftpack.fftshift ( np.linspace(min_Py, np.abs(min_Py) - dPy, self.Y_gridDIM ))


		px = px_Vector[np.newaxis,:]
		py = py_Vector[:,np.newaxis]


		sqrtp = sign*2*np.sqrt( self.mass*self.mass*self.c**4 + self.c*self.c*px*px + self.c*self.c*py*py )
		aa = sign*self.mass*self.c*self.c/sqrtp
		bb = sign*(px/sqrtp - 1j*py/sqrtp)
		cc = sign*(px/sqrtp + 1j*py/sqrtp)
	        
		ElectronProjector = np.matrix([ [0.5+aa , 0.  , 0.  , bb  ],
						[0. , 0.5+aa  , cc  , 0.  ],
						[0. , bb  , 0.5-aa  , 0.  ],
						[cc , 0.  , 0.  , 0.5-aa] ])

		psi1_fft = fftpack.fft2( self.Psi1_init   ) 
		psi2_fft = fftpack.fft2( self.Psi2_init   ) 
		psi3_fft = fftpack.fft2( self.Psi3_init   ) 
		psi4_fft = fftpack.fft2( self.Psi4_init   ) 		
		
		psi1_fft_electron = ElectronProjector[0,0]*psi1_fft + ElectronProjector[0,1]*psi2_fft +\
		ElectronProjector[0,2]*psi3_fft + ElectronProjector[0,3]*psi4_fft	

		psi2_fft_electron = ElectronProjector[1,0]*psi1_fft + ElectronProjector[1,1]*psi2_fft +\
		ElectronProjector[1,2]*psi3_fft + ElectronProjector[1,3]*psi4_fft

                psi3_fft_electron = ElectronProjector[2,0]*psi1_fft + ElectronProjector[2,1]*psi2_fft +\
		ElectronProjector[2,2]*psi3_fft + ElectronProjector[2,3]*psi4_fft	

                psi4_fft_electron = ElectronProjector[3,0]*psi1_fft + ElectronProjector[3,1]*psi2_fft +\
		ElectronProjector[3,2]*psi3_fft + ElectronProjector[3,3]*psi4_fft

                self.Psi1_init  = fftpack.ifft2( psi1_fft_electron   ) 
		self.Psi2_init  = fftpack.ifft2( psi2_fft_electron   ) 
		self.Psi3_init  = fftpack.ifft2( psi3_fft_electron   ) 
		self.Psi4_init  = fftpack.ifft2( psi4_fft_electron   ) 					

def integration_and_analysis(zeta0, phi0, f_zeta_t, f_phi_t, M,
                             kxgrid, kygrid, ktgrid, dt, time, damping,
                             storage, wmax, wmin):
    '''
    Perform the 4th order Runge-Kutta integration scheme for
    surface and potential.

    zeta0: surface at time step n.
    phi0: potential at time step n.
    f_zeta_t, f_phi_t: functions with Euler eqs of M-order to be solved
    found with derive_euler_equation_functions(M).
    time: time in which the integration takes place
    damping: non-linear damping factor in Euler eqs.

    Return: surface and potential at time step n+1. 
    Other variables like the orders of phi can easily be returned!
    '''

    zeta = zeta0
    phi = phi0

    rk1_zeta, rk1_phi, phi_m = tderiv_surface_potential(zeta, phi, f_zeta_t, f_phi_t, M,
                                                 kxgrid, kygrid, ktgrid, time, damping,
                                                 return_phi_m = 1)
    zeta = zeta0 + rk1_zeta*dt/2
    phi = phi0 + rk1_phi*dt/2

    rk2_zeta, rk2_phi = tderiv_surface_potential(zeta, phi, f_zeta_t, f_phi_t, M,
                                                 kxgrid, kygrid, ktgrid, time + dt/2, damping)
    zeta = zeta0 + rk2_zeta*dt/2
    phi = phi0 + rk2_phi*dt/2

    rk3_zeta, rk3_phi = tderiv_surface_potential(zeta, phi, f_zeta_t, f_phi_t, M,
                                                 kxgrid, kygrid, ktgrid, time + dt/2, damping)
    zeta = zeta0 + rk3_zeta*dt
    phi = phi0 + rk3_phi*dt

    rk4_zeta, rk4_phi = tderiv_surface_potential(zeta, phi, f_zeta_t, f_phi_t, M,
                                                 kxgrid, kygrid, ktgrid, time + dt, damping)

    dzeta_dt = 1/6 * (rk1_zeta + rk4_zeta + 2*(rk2_zeta + rk3_zeta))
    dphi_dt = 1/6 * (rk1_phi + rk4_phi + 2*(rk2_phi + rk3_phi))

    
    kernel = monitor_conserved_quantities(phi0, zeta0, dzeta_dt, kxgrid, kygrid)
    print('Total Energy: ' + str(kernel['kin'] + kernel['poten']) + ' Total Mass:' +  str(kernel['mass']))

    zeta_next = np.real(ifft2(dealias(fft2(zeta0 + dt*dzeta_dt), M))) #We transform, dealias, then transform back as real. 
    phi_next = np.real(ifft2(dealias(fft2(phi0 + dt*dphi_dt), M)))

    storage, wmax, wmin = detect_rogue_waves(zeta0, zeta_next, wmax, wmin, sig_h, storage, time)

    return zeta_next, phi_next, phi_m, storage, wmax, wmin

def main():
    fftsize=1024
    #img=io.imread("torre.jpg")
    #print "Shape: ",img.shape
    #im = np.mean(img,axis=2)/255.    
    #im_fft=fftpack.fft2(im,(fftsize, fftsize))
    #F = np.log(1+np.abs(im_fft))
    #recovered = np.real(fftpack.ifft2(im_fft))
    a=np.zeros((3,3),dtype=float)
    b=np.ones((3,3),dtype=float)
    a[0][0]=2    
    a[0][1]=3
    a[0][2]=1
    #plt.show()
    #plt.imshow(im, cmap='gray')
    #plt.title('Imagen en gris')
    #blurImg(im,fftsize)

    imA = io.imread("torre.jpg")


    im = np.mean(imA,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_sqg(self):
        import scipy.fftpack as fft
        import numpy as np
        from math import pi
        
        self.precondition()
        
        dx = self.dx
        dy = self.dy
        rhos = self.ssd
        
        bhat = fft.fft2( - 9.81 * rhos / self.rho0)  # calculate bouyance 
        ny, nx = rhos.shape
        nz = self.nz
        k = 2 * pi * fft.fftfreq(nx)
        l = 2 * pi * fft.fftfreq(ny)

        ipsihat = np.zeros((nz+3, ny, nx))*complex(0, 0)
        
        Q = np.zeros((nz + 1, 1), dtype='float64'); Q[[0,-1]] = 0.0 # for interior PV, no used in this version
        
        # cutoff value
        ck, cl = 2 * pi / self.filterL, 2 * pi / self.filterL
        # loop through wavenumbers
        bhats = np.zeros_like(bhat)
        for ik in np.arange(k.size):
            for il in np.arange(l.size):
                wv2 = ((k[ik] / dx[il, 0]) ** 2 + (l[il] / dy[0, ik]) ** 2)
                if wv2 > (ck * ck + cl * cl):
                    bhats[il,ik] = bhat[il,ik]
                    right = - bhat[il, ik] / self.f0 * self.Rp
                    left = self.M - wv2 * np.eye(self.nz+1)
                    ipsihat[1:-1, il, ik] = np.linalg.solve(left, right).flatten()
                else:
                    print 'skip k(ik,il)', ik, il, "wv2 = ", wv2
        
        for k in range(1,nz+2):
            ipsihat[k, :, :] = (fft.ifft2(ipsihat[k, :, :]))
        
        if self.bottomboundary == 'psi=0':
            self.psis = np.r_[(np.real(ipsihat)), np.zeros((1,ny,nx))]
        else:
            self.psis = np.real(ipsihat)
            self.psis[0,:,:]= self.psis[1,:,:]
            self.psis[-1,:,:]=self.psis[-2,:,:]-self.dzc[-1]*np.real(fft.ifft2(-bhats))/self.f0
            
        self.rhos = self.psi2rho(self.psis)
        self.us, self.vs = psi2uv(self.lon, self.lat, self.psis)
        
        return

def make_initial_surface(spectrum, ktgrid, dkx_dky):
    '''Make an initial surface and potential from a spectrum
    by giving each mode a random phase'''

    random_phase = (2*pi)*np.random.random(np.shape(spectrum))
    amplitude_k = np.sqrt(2*spectrum*dkx_dky)
    wk = np.sqrt(ktgrid*grav)
    wk[0,0] = 1
    surface_k = amplitude_k * np.exp(1j * random_phase)
    potential_k = 1j * grav / wk * surface_k
    potential_k[0,0] = 0.0
    surface = np.real(ifft2(surface_k))
    potential = np.real(ifft2(potential_k))

    return surface, potential

    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