Python fftpack.ifft2函数代码示例

 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   )