Python ndimage.gaussian_gradient_magnitude函数代码示例

def test_multiple_modes():
    # Test that the filters with multiple mode cababilities for different
    # dimensions give the same result as applying a single mode.
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    mode1 = 'reflect'
    mode2 = ['reflect', 'reflect']

    assert_equal(sndi.gaussian_filter(arr, 1, mode=mode1),
                 sndi.gaussian_filter(arr, 1, mode=mode2))
    assert_equal(sndi.prewitt(arr, mode=mode1),
                 sndi.prewitt(arr, mode=mode2))
    assert_equal(sndi.sobel(arr, mode=mode1),
                 sndi.sobel(arr, mode=mode2))
    assert_equal(sndi.laplace(arr, mode=mode1),
                 sndi.laplace(arr, mode=mode2))
    assert_equal(sndi.gaussian_laplace(arr, 1, mode=mode1),
                 sndi.gaussian_laplace(arr, 1, mode=mode2))
    assert_equal(sndi.maximum_filter(arr, size=5, mode=mode1),
                 sndi.maximum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.minimum_filter(arr, size=5, mode=mode1),
                 sndi.minimum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.gaussian_gradient_magnitude(arr, 1, mode=mode1),
                 sndi.gaussian_gradient_magnitude(arr, 1, mode=mode2))
    assert_equal(sndi.uniform_filter(arr, 5, mode=mode1),
                 sndi.uniform_filter(arr, 5, mode=mode2))

    def gaussian_kernel(self,xvalues,yvalues,r200,normalization=100,scale=10,xres=200,yres=220,xmax=6.0,ymax=5000.0,adj=20):
        """
        Uses a 2D gaussian kernel to estimate the density of the phase space.
        As of now, the maximum radius extends to 6Mpc and the maximum velocity allowed is 5000km/s
        The "q" parameter is termed "scale" here which we have set to 10 as default, but can go as high as 50.
        "normalization" is simply H0
        "x/yres" can be any value, but are recommended to be above 150
        "adj" is a custom value and changes the size of uniform filters when used (not normally needed)
        """
        self.x_scale = xvalues/xmax*xres
        self.y_scale = ((yvalues+ymax)/(normalization*scale))/((ymax*2.0)/(normalization*scale))*yres

        img = np.zeros((xres+1,yres+1))
        self.x_range = np.linspace(0,xmax,xres+1)
        self.y_range = np.linspace(-ymax,ymax,yres+1) 

        for j in range(xvalues.size):
            img[self.x_scale[j],self.y_scale[j]] += 1
        
        #Estimate kernel sizes
        #Uniform
        #self.ksize = 3.12/(xvalues.size)**(1/6.0)*((np.var(self.x_scale[xvalues<r200])+np.var(self.y_scale[xvalues<r200]))/2.0)**0.5/adj
        #if self.ksize < 3.5:
        #    self.ksize = 3.5
        #Gaussian
        self.ksize_x = (4.0/(3.0*xvalues.size))**(1/5.0)*np.std(self.x_scale[xvalues<r200])
        self.ksize_y = (4.0/(3.0*yvalues.size))**(1/5.0)*np.std(self.y_scale[xvalues<r200])
        
        #smooth with estimated kernel sizes
        #img = ndi.uniform_filter(img, (self.ksize,self.ksize))#,mode='reflect')
        self.img = ndi.gaussian_filter(img, (self.ksize_y,self.ksize_x),mode='reflect')
        self.img_grad = ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x))
        self.img_inf = ndi.gaussian_gradient_magnitude(ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x)), (self.ksize_y,self.ksize_x))

def gborders(img, alpha=1.0, sigma=1.0):
    """Stopping criterion for image borders."""
    # The norm of the gradient.
    gausgradm = np.zeros(img.shape, np.double)
    gaussian_gradient_magnitude(img, sigma, output=gausgradm, mode='constant')
	
    thr = 2500
    mask = gausgradm > thr
    gausgradm[mask] = thr
    return 1.0/np.sqrt(1.0 + alpha*gausgradm)

    def pruning(self, skeleton_img, sigma):
        skeleton_img_mat = image_conversion.cv2array(skeleton_img)
        # Ausgabe-Array fuer das Ergebnis der Gradientberechnung
        gradient_output = numpy.empty_like(skeleton_img_mat)
        # Gradienten-Berechnung
        ndimage.gaussian_gradient_magnitude(skeleton_img_mat, sigma, gradient_output)
        # Normalisierung
        gradient_output /= gradient_output.max()
        # Array ins Bild umwandeln
        grad_img = image_conversion.array2cv(gradient_output)
        # Schwellwertbasierte Segmentierung des Gradientbildes
        dist_gradient_thresh = cv.CreateImage(cv.GetSize(grad_img), 8, 1)
        cv.InRangeS(grad_img, 0.6, 1, dist_gradient_thresh)

        return dist_gradient_thresh

def inverse_gaussian_gradient(image, alpha=100.0, sigma=5.0):
    """Inverse of gradient magnitude.

    Compute the magnitude of the gradients in the image and then inverts the
    result in the range [0, 1]. Flat areas are assigned values close to 1,
    while areas close to borders are assigned values close to 0.

    This function or a similar one defined by the user should be applied over
    the image as a preprocessing step before calling
    `morphological_geodesic_active_contour`.

    Parameters
    ----------
    image : (M, N) or (L, M, N) array
        Grayscale image or volume.
    alpha : float, optional
        Controls the steepness of the inversion. A larger value will make the
        transition between the flat areas and border areas steeper in the
        resulting array.
    sigma : float, optional
        Standard deviation of the Gaussian filter applied over the image.

    Returns
    -------
    gimage : (M, N) or (L, M, N) array
        Preprocessed image (or volume) suitable for
        `morphological_geodesic_active_contour`.
    """
    gradnorm = ndi.gaussian_gradient_magnitude(image, sigma, mode='nearest')
    return 1.0 / np.sqrt(1.0 + alpha * gradnorm)

def VesicleEdge_phc(img, x0, y0, r0, N=100, phi1=0, phi2=2 * np.pi, sigma=1):
    Xedge = np.empty(N)
    Yedge = np.empty(N)
    for i, phi in enumerate(np.linspace(phi1, phi2, N)):
        x = x0 + r0 * np.cos(phi)
        y = y0 + r0 * np.sin(phi)
        if x < 0:
            x = 0
            y = y0 + (x - x0) * np.tan(phi)
        elif x > img.shape[1] - 1:
            x = img.shape[1] - 1
            y = y0 + (x - x0) * np.tan(phi)
        if y < 0:
            y = 0
            x = x0 + (y - y0) / np.tan(phi)
        elif y > img.shape[0] - 1:
            y = img.shape[1] - 1
            x = x0 + (y - y0) / np.tan(phi)

        point1 = np.asarray(((y0, x0), (PIX_ERR, PIX_ERR)))
        point2 = np.asarray(((y, x), (PIX_ERR, PIX_ERR)))
        metric, metric_err, line = section_profile(img, point1, point2)
        grad = gaussian_gradient_magnitude(line, sigma)
        pos = np.argmax(grad)
        Xedge[i] = x0 + pos * np.cos(phi) * metric
        Yedge[i] = y0 + pos * np.sin(phi) * metric
    return Xedge, Yedge

 def remove_dust(self):
     self.seq = []
     for i in xrange(0,self.n_frames):
         # greyscale conversion
         img = np.average(self.reader.get_data(i),axis=2)
         self.seq.append(img)
     self.seq = np.array(self.seq)
     #var = np.var(self.seq, axis=0)
     #min = np.min(self.seq, axis=0)
     #max = np.max(self.seq, axis=0)
     #delta = max - min
     #var = stats.variation(self.seq, axis=0)
     #gmean = stats.gmean(self.seq, axis=0)
     a = np.average(self.seq, axis=0)
     #grad = ndimage.gaussian_gradient_magnitude(a , 0.25)
     #map = ndimage.prewitt(a)
     map = ndimage.gaussian_laplace(a,2.5) * ndimage.gaussian_gradient_magnitude(a , 0.25)
     cutoff = np.percentile(map,99.9)
     map[map<cutoff]=0
     map[map>0]=1
     #map = grad
     #map[map>300]=300
     fig = plt.figure(figsize=(20,8), frameon=False)
     fig.subplots_adjust(hspace=0)
     fig.subplots_adjust(wspace=0)
     ax1 = fig.add_subplot(1, 2, 1)
     ax1.imshow(map,interpolation='nearest')
     ax1.set_title('variance')
     ax2 = fig.add_subplot(1, 2, 2)
     ax2.imshow(self.seq[0], cmap='Greys_r',interpolation='nearest')
     ax2.set_title('img')
     fig.set_tight_layout(True)
     plt.show()

 def __init__(self, *args,**kwargs):
     im = kwargs.pop('image')
     super(Variance_no_D_Cmap,self).__init__(*args, **kwargs)
     
     fim = ndimage.gaussian_gradient_magnitude(im,2)
     fim = fim/np.std(fim)
     alpha = 1e0
     self.cmap = np.exp(-fim.flat[self.anchor['imask']]*alpha) + 1e-10

	def gaussian_kernel(self,xvalues,yvalues,r200,normalization=100,scale=10,res=200,adj=20,see=False):
		yres = 220
		#x_scale = (xvalues-np.min(xvalues))/np.max(xvalues-np.min(xvalues))*res
		#y_scale = ((yvalues-np.min(yvalues))/(normalization*scale))/np.max(xvalues-np.min(xvalues))*res
		self.x_scale = xvalues/6.0*res
		self.y_scale = ((yvalues+5000)/(normalization*scale))/(10000.0/(normalization*scale))*yres

		#img = np.zeros((int(np.max(x_scale))+1,int(np.max(y_scale))+1))
		img = np.zeros((res+1,yres+1))
		#x_range = np.linspace(np.min(xvalues),np.max(xvalues),int(np.max(x_scale))+1)
		#y_range = np.linspace(np.min(yvalues),np.max(yvalues),int(np.max(y_scale))+1)
		x_range = np.linspace(0,6,res+1)
		y_range = np.linspace(-5000,5000,yres+1) 

		for j in range(xvalues.size):
			img[self.x_scale[j],self.y_scale[j]] += 1

		#pcolormesh(img.T)
		#find ksize
		#xval = xvalues[np.where((xvalues<3) & (yvalues<2000) & (yvalues > -2000))]
		#yval = yvalues[np.where((xvalues<3) & (yvalues<2000) & (yvalues > -2000))]
		#x_scale2 = (xval-np.min(xval))/np.max(xval-np.min(xval))*res
		#y_scale2 = ((yval-np.min(yval))/(normalization*scale))/np.max(xval-np.min(xval))*res
		#xksize = 3.12/(xvalues.size)**(1.0/6.0)*((np.var(x_scale))/2.0)**0.5/adj
		#yksize = 3.12/(xvalues.size)**(1.0/6.0)*((np.var(y_scale))/2.0)**0.5/adj
	
		self.ksize = 3.12/(xvalues.size)**(1/6.0)*((np.var(self.x_scale[xvalues<r200])+np.var(self.y_scale[xvalues<r200]))/2.0)**0.5/adj
                self.ksize_x = (4.0/(3.0*xvalues.size))**(1/5.0)*np.std(self.x_scale[xvalues<r200])
                self.ksize_y = (4.0/(3.0*yvalues.size))**(1/5.0)*np.std(self.y_scale[xvalues<r200])
                if self.ksize < 3.5:
                    self.ksize = 3.5
		#ksize = 6.77588630223
		#print 'kernel size',ksize
		#img = ndi.uniform_filter(img, (self.ksize,self.ksize))#,mode='reflect')
                img = ndi.gaussian_filter(img, (self.ksize_y,self.ksize_x))#,mode='reflect')
                img_grad = ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x))
                img_inf = ndi.gaussian_gradient_magnitude(ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x)), (self.ksize_y,self.ksize_x))
		
#		if see == True:
#			s = figure()
#			ax = s.add_subplot(111)
#			ax.pcolormesh(x_range,y_range,img.T)
#			show()
		return (x_range,y_range,img,np.abs(img_grad),np.abs(img_inf))

	def onEdgeFinding(self,evt):	
		if not self.panel.selectiontool.isTargeting('Auto Create Contours'):
			return 
		point = self.panel.view2image((evt.m_x,evt.m_y))

		ndarray = mrc.read(os.path.join(self.appionloop.params['rundir'], self.appionloop.imgtree[self.index]['filename']+'.dwn.mrc'))
		mrc.write(ndarray, os.path.join(self.appionloop.params['rundir'], 'beforefilter'+'.dwn.mrc'))

		negative = False	
		if self.filters:
			ndarray = ndimage.gaussian_filter(ndarray,1)
			ndarray = ndimage.gaussian_gradient_magnitude(ndarray,2)
			markers = []
			for i in range(3):
				for j in range(3):
					if i!=0 or j!=0:
						markers.append((point[0]-1+i,point[1]-1+j))
			markers = (1,2,3,4,5,6,7,8)
			#ndarray = ndimage.watershed_ift(ndarray,markers)
			ndarray = ndimage.laplace(ndarray)
			ndarray = ndimage.gaussian_filter(ndarray,1)
			#ndarray = apImage.preProcessImage(ndarray,params=self.appionloop.params)
			negative = True 
		mrc.write(ndarray, os.path.join(self.appionloop.params['rundir'], 'afterfilter'+'.dwn.mrc'))

		delta = .1
		targets = []
		radius = 20
		size = 50
		rangeSize = 50
		maker = PixelCurveMaker()
		maker._init_(size,rangeSize);
		for theta in range(size):
			theta +=0
			theta*=math.pi*2/rangeSize
			for rad in range(size):
				try:
					if negative:
						maker.addData(theta,rad,127-ndarray[int(point[1]+rad*math.sin(theta))][int(point[0]+rad*math.cos(theta))])
					else:
						maker.addData(theta,rad,ndarray[int(point[1]+rad*math.sin(theta))][int(point[0]+rad*math.cos(theta))])
				except IndexError:
					maker.addData(theta,rad,0)
		maker.makeCalculations()
		s = self.filterSelectorChoices[self.filterSelector.GetSelection()]	
		dilate = 2
		if s == 'Latex Bead':
			dilate = 0
		for theta in range(size):
			theta += 0
			theta*=math.pi*2/rangeSize
			targets.append((point[0]+(dilate+maker.getData(theta))*math.cos(theta),point[1]+(dilate+maker.getData(theta))*math.sin(theta)))
		self.addPolyParticle(targets)
		#this section draws all of the contours that the algorithm considers - useful for debugging
		'''

def getGradientVideo(I, IDims, sigma = 1):
    GV = np.zeros(I.shape)
    for i in range(I.shape[0]):
        X = np.reshape(I[i, :], IDims)
        G = rgb2gray(X, False)
        GM = gaussian_gradient_magnitude(G, sigma)
        F = np.zeros(IDims)
        for k in range(F.shape[2]):
            F[:, :, k] = GM
        GV[i, :] = F.flatten()
    return GV

def hipass(image,gaussrad,dilrat):
    # 
    # A derivative of gaussian (suggested radii of 0.5) to get edges
    # followed by a morphological dilate to widen the effect and pick up extra pixels on edges
    # This can be used as a mask for 

    timg = np.copy(image).reshape(96,96)
    
    gradim = ndimage.gaussian_gradient_magnitude(timg,gaussrad)
    
    digradim= ndimage.morphology.grey_dilation(gradim,size=(dilrat,dilrat))
    
    return digradim.flatten()

def test_multiple_modes_gaussian_gradient_magnitude():
    # Test gaussian_gradient_magnitude filter for multiple
    # extrapolation modes
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    expected = np.array([[0.04928965, 0.09745625, 0.06405368],
                         [0.23056905, 0.14025305, 0.04550846],
                         [0.19894369, 0.14950060, 0.06796850]])

    modes = ['reflect', 'wrap']

    calculated = sndi.gaussian_gradient_magnitude(arr, 1, mode=modes)

    assert_almost_equal(expected, calculated)

def task2_4():
    img = normalize_intensity(imread(CAMERAMAN))
    img = img[30:95, [i for i in range(80, 160)]]  # select subsection of image
    vel_x, vel_y = gradient(f=img)

    magn_img = gaussian_gradient_magnitude(img, 3)
    output_path = os.path.join(OUTPUT_DIR, "2_4_gradien_magnitude_" + os.path.split(CAMERAMAN)[-1])
    imsave(output_path, magn_img)

    dim_x, dim_y = len(img[0]), len(img)
    x, y = range(dim_x), range(dim_y)
    x, y = meshgrid(x, y)
    plt.figure()
    imgplot = plt.imshow(img)
    imgplot.set_cmap('gray')
    plt.ylim(dim_y, 0)
    plt.quiver(x, y, vel_x, vel_y, pivot='middle')
    plt.show()

def wall_points_pix_3(img, refsx, axisy, sigma):
    N = 2
    refs = np.array([])
    axisy = np.array([113, 115])
    for i, refx in enumerate(refsx):
        prof = img[:, refx]
        mid = axisy[i]
        #        mid = len(prof)/2
        filtered = ndimage.gaussian_gradient_magnitude(ndimage.sobel(prof), sigma)

        refy = np.asarray((np.argmax(filtered[:mid]), np.argmax(filtered[mid:]) + mid))
        dref = np.asarray([PIX_ERR, 0])
        rx = np.tile(refx, N)
        xy = np.column_stack((refy, rx))  # .flatten()
        drefe = np.repeat(np.expand_dims(dref, 0), N, 0)
        ref = np.concatenate((xy, drefe), 1)
        refs = np.append(refs, ref)
    return refs.reshape(-1, 2, 2)