Python ndimage.binary_opening函数代码示例

def analyse(im, scales=[4,10,14,40], verbose=True):
    try:
        mask = im > im.mean()
    except:
        im = im.matrix
        mask = im > im.mean()

    granulo = granulometry(mask, sizes=np.arange(2, 19, 4))
    print 'granulo:', granulo

    plt.figure(figsize=(6, 2.2))
    plt.subplot(121)
    plt.imshow(mask, cmap=plt.cm.gray, origin='lower')
    opened_less = ndimage.binary_opening(mask, structure=disk_structure(scales[0]))
    opened = ndimage.binary_opening(mask, structure=disk_structure(scales[1]))
    opened_more = ndimage.binary_opening(mask, structure=disk_structure(scales[2]))
    opened_even_more = ndimage.binary_opening(mask, structure=disk_structure(scales[3]))
    plt.contour(opened_less, [0.5], colors='g', linewidths=2)
    plt.contour(opened, [0.5], colors='b', linewidths=2)
    plt.contour(opened_more, [0.5], colors='r', linewidths=2)
    plt.contour(opened_even_more, [0.5], colors='k', linewidths=2)

    plt.axis('off')
    plt.subplot(122)
    plt.plot(np.arange(2, 19, 4), granulo, 'ok', ms=8)
    plt.subplots_adjust(wspace=0.02, hspace=0.15, top=0.95, bottom=0.15, left=0, right=0.95)
    if verbose:
        plt.show()
    
    return opened_less, opened, opened_more, opened_even_more

def open_isl(isls, crms):
    import pylab as pl
    import scipy.ndimage as nd
    import numpy as N

    thr = crms
    ft1 = N.array(((1,0,1), (0,1,0), (1,0,1)), int)
    ft2 = N.array(((0,1,0), (1,1,1), (0,1,0)), int)
    ft3 = N.ones((3,3), int)
    ft5 = N.ones((5,5), int)
    for isl in isls:
      ma = ~isl.mask_active
      open1 = nd.binary_opening(ma, ft1)
      open2 = nd.binary_opening(ma, ft2)
      open3 = nd.binary_opening(ma, ft3)
      open5 = nd.binary_opening(ma, ft5)

      pl.figure()
      pl.suptitle('Island '+str(isl.island_id))
      pl.subplot(2,2,1); pl.imshow(N.transpose(isl.image), origin='lower', interpolation='nearest'); pl.title('Image')
      pl.subplot(2,2,2); pl.imshow(N.transpose(ma), origin='lower', interpolation='nearest'); pl.title('mask')
      pl.subplot(2,2,3); pl.imshow(N.transpose(open3), origin='lower', interpolation='nearest'); pl.title('open 3x3')
      pl.subplot(2,2,4); pl.imshow(N.transpose(open5), origin='lower', interpolation='nearest'); pl.title('open 5x5')
      #pl.subplot(2,2,3); pl.imshow(N.transpose(open1), origin='lower', interpolation='nearest'); pl.title('open diag')
      #pl.subplot(2,2,4); pl.imshow(N.transpose(open2), origin='lower', interpolation='nearest'); pl.title('open str')
      pl.savefig('cyga_p_w12_bigisl_'+str(isl.island_id)+'_open.png')

def fraction_positive(bin_im, positive_im, erode=2, overlap_thresh=0.9,
                      bin_name='nuclei', positive_name='tf'):
    """Compute fraction of objects in bin_im overlapping positive_im.

    The purpose of this function is to compute the fraction of nuclei
    that express a particular transcription factor. By providing the
    thresholded DAPI channel as `bin_im` and the thresholded TF channel
    as `positive_im`, this fraction can be computed.

    Parameters
    ----------
    bin_im : 2D array of bool
        The image of objects being tested.
    positive_im : 2D array of bool
        The image of positive objects.
    erode : int, optional
        Radius of structuring element used to smooth input images.
    overlap_thresh : float, optional
        The minimum amount of overlap between an object in `bin_im` and
        the `positive_im` to consider that object "positive".
    bin_name : string, optional
        The name of the objects being tested.
    positive_name : string, optional
        The name of the property being measured.

    Returns
    -------
    f : 1D array of float, shape (1,)
        The feature vector.
    name : list of string, length 1
        The name of the feature.

    Examples
    --------
    >>> bin_im = np.array([[1, 1, 0],
    ...                    [0, 0, 0],
    ...                    [1, 1, 1]], dtype=bool)
    >>> pos_im = np.array([[1, 0, 0],
    ...                    [0, 1, 1],
    ...                    [0, 1, 1]], dtype=bool)
    >>> f = fraction_positive(bin_im, pos_im, erode=0, overlap_thresh=0.6)
    >>> f[0]
    array([0.5])
    >>> f[1][0]
    'frac-nuclei-pos-tf-erode-0-thresh-0.60'
    """
    selem = skmorph.disk(erode)
    if erode > 0:
        bin_im = nd.binary_opening(bin_im, selem)
        positive_im = nd.binary_opening(positive_im, selem)
    lab_im = nd.label(bin_im)[0].ravel()
    pos_im = positive_im.ravel().astype(int)
    counts = sparse.coo_matrix((np.ones(lab_im.size),
                                (lab_im, pos_im))).todense()
    means = counts[:, 1] / np.sum(counts, axis=1)
    f = np.array([np.mean(means[1:] > overlap_thresh)])
    name = ['frac-%s-pos-%s-erode-%i-thresh-%.2f' %
            (bin_name, positive_name, erode, overlap_thresh)]
    return f, name

  def findCom(self,data): 
    data = data.astype(int)
    if self.update_counter >= 5:
      self.update_counter = 0

      ##########################################################################
      ## Update the background image, adding a new image and removing the oldest.
      ##########################################################################
      self.background_list.insert(0,data)
      self.background_list.pop()
      background = np.zeros((480, 640, 3), dtype=int)
      for b in self.background_list:
        background += b
      self.background = background/len(self.background_list)
    
    ############################################################################
    ## Detect foreground by looking at difference from mean.
    ############################################################################
    foreground = np.sum(np.abs(np.subtract(self.background,data)),axis=2)
    falseImage = foreground
    ## clean foreground image
    falseImage[falseImage > 100] = 255
    falseImage[falseImage < 101] = 0
    falseImage = ndimage.binary_opening(falseImage)
    falseImage = ndimage.binary_closing(falseImage)
    com = ndimage.measurements.center_of_mass(falseImage)
    self.update_counter += 1
    return com

def get_difference_spots(pix):
    bpix = pix > 20
    bpix = ndimage.binary_opening(bpix)
    bpix = ndimage.binary_closing(bpix)
    labels, n = ndimage.measurements.label(bpix)
    clicks = ndimage.measurements.center_of_mass(pix, labels, range(1, n+1))
    return clicks

def detect_vortices(cloud, radius=70, showplots=False):
    """
    Detects whether there are vortex-like features within a given radius
    of the peak density in the TOF image of an expanded BEC
    """
    OD = cloud.get_OD()    
    peak_coord = cloud.results['peak coordinates']
    center_region = ROI(center=peak_coord,
                        size=(1.5 * radius, 1.5 * radius)).slices
    smooth_cloud = ndi.median_filter(OD[center_region], size=4)
    minOD = smooth_cloud.min()
    maxOD = smooth_cloud.max()
    cloud_median = ndi.median_filter(smooth_cloud, size=10)
    belowthresh = where(smooth_cloud < cloud_median * 0.75, 1, 0)
    opened = ndi.binary_opening(belowthresh, iterations=1)
    closed = ndi.binary_closing(opened, iterations=1)
    vort_found = ndi.label(closed)[1]
    cloud.results['vort_found'] = vort_found
    if showplots == True:
        fig = plt.figure(1999)
        fig.add_subplot(221, xticks=[], yticks=[])   
        plt.imshow(smooth_cloud, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(222, xticks=[], yticks=[]) 
        plt.imshow(cloud_median, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(223, xticks=[], yticks=[]) 
        plt.imshow(closed, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
        fig.add_subplot(224, xticks=[], yticks=[]) 
        plt.imshow(belowthresh, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
    return vort_found

	def cell_to_image(self):
	
		# Find x, y coordinate bounds
		x_res = max(self.pix_list, key=itemgetter(0))[0]
		y_res = max(self.pix_list, key=itemgetter(1))[1]

		# Creating labeled_img
		self.cell_img = NP.zeros([x_res+2, y_res+2], dtype=NP.int_)
		
		for (x_pix, y_pix) in self.pix_list:
			self.cell_img[x_pix-1, y_pix-1] = 1

		# Find the pixels that make up the perimeter
		eroded_image = NDI.binary_erosion(self.cell_img)

		eroded_image_open = NDI.binary_opening(eroded_image, structure=NP.ones((3,3)))
		eroded_image_open2 = NDI.binary_erosion(eroded_image_open)

		# self.perim_img = self.cell_img - eroded_image
		self.eroded_img = eroded_image_open - eroded_image_open2
		self.perim_img = self.cell_img - eroded_image

		# Create a list of the coordinates of the pixels (use the center of the pixels)
		perim_image_ind = NP.where(self.perim_img == 1)
		perim_image_coord = NP.array([perim_image_ind[0], perim_image_ind[1]])
		self.perim_coord = NP.transpose(perim_image_coord)

		return

def segmentAlg(I,frame):
    n = 10
    np.random.seed(1)

    Nx, Ny = I.size

    im = np.asarray(I)

    mask = (im > im.mean()).astype(np.float)
    mask += 1 * im


    # Thresholding algorithm
    #binary_img = white.white(I,Nx,Ny,20,1.15)#mask < 80
    binary_img = im > 100
    #scipy.misc.imsave('images/outfile'+str(frame)+'_m.jpg', binary_img)

    binary_img = ndimage.binary_opening(binary_img, structure=np.ones((2,2))).astype(np.int)
    #scipy.misc.imsave('images/outfile'+str(frame)+'_m2.jpg', binary_img)

    label_im, nb_labels = ndimage.label(binary_img)

    sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))
    mask_size = sizes > 150000
    remove_pixel = mask_size[label_im]
    label_im[remove_pixel] = 0

    #scipy.misc.imsave('images/outfile'+str(frame)+'_m3.jpg', label_im)

    return binary_img, label_im, nb_labels, sizes

def binaryOpening(binarydata, structure=None, iterations=3):
    '''
    Opening [R52] is a mathematical morphology operation [R53] that consists
    in the succession of an erosion and a dilation of the input with the same
    structuring element. Opening therefore removes objects smaller than the
    structuring element.

    Together with closing (binary_closing), opening can be used for noise removal.
    '''
    result = np.empty_like(binarydata)
    
    if structure is None:
        ndimage.binary_opening(binarydata, iterations=iterations, output=result)
    else:
        ndimage.binary_opening(binarydata, structure, iterations=iterations, output=result)
    return result

def og_cb(og_msg, param_list):
    global occupancy_difference_threshold, connected_comonents_size_threshold
    rospy.loginfo('og_cb called')
    diff_og = param_list[0]
    curr_og = rog.og_msg_to_og3d(og_msg, to_binary = False)

    if diff_og == None:
        param_list[0] = curr_og
        return

    pog.subtract(diff_og, curr_og)
    param_list[0] = curr_og

    diff_og.to_binary(occupancy_difference_threshold)
    # filter the noise
    connect_structure = np.zeros((3,3,3), dtype=int)
    connect_structure[1,1,:] = 1
#    connect_structure[1,1,0] = 0
    diff_og.grid = ni.binary_opening(diff_og.grid, connect_structure,
                                     iterations = 1)
    #    diff_og.grid, n_labels = diff_og.connected_comonents(connected_comonents_size_threshold)

    print 'np.all(diff_og == 0)', np.all(diff_og.grid == 0)
    diff_og_msg = rog.og3d_to_og_msg(diff_og)
    diff_og_msg.header.frame_id = og_msg.header.frame_id
    diff_og_msg.header.stamp = og_msg.header.stamp
    param_list[1].publish(diff_og_msg)

def smoothImage(image,disk_size):
    '''
    Smooths the outline of a binary object by using morphological opening with
    a user-defined kernel size.
    '''
    smoothed = ndimage.binary_opening(image,structure=np.ones((disk_size,disk_size))).astype('uint8') * 255
    return smoothed

def find_feature_mask_simple(s_msk, sigma=1, ax=None, x_values=None):
    """
    sigma is estimated globally.
    """
    # find emission features from observed spec.

    filtered_spec = s_msk - ni.median_filter(s_msk, 15)
    #filtered_spec = ni.gaussian_filter1d(filtered_spec, 0.5)
    #smoothed_std = get_smoothed_std(filtered_spec,
    #                                rad=3, smooth_length=3)
    std = np.nanstd(filtered_spec)
    for i in [0, 1]:
        std = filtered_spec[np.abs(filtered_spec)<3*std].std()

    emission_feature_msk_ = filtered_spec > sigma*std
    #emission_feature_msk_ = ni.binary_closing(emission_feature_msk_)
    emission_feature_msk = ni.binary_opening(emission_feature_msk_,
                                             iterations=1)

    if ax is not None:
        if x_values is None:
            x_values = np.arange(len(s_msk))

        #ax.plot(x_values, s_msk)
        ax.plot(x_values, filtered_spec)
        #ax.plot(x_values, smoothed_std)
        ax.axhline(sigma*std)
        ax.plot(x_values[emission_feature_msk],
                emission_feature_msk[emission_feature_msk],
                "ys", mec="none")

    return emission_feature_msk

def granulometry(data, sizes=None):
    s = max(data.shape)
    if sizes == None:
        sizes = range(1, s/2, 2)
    granulo = [ndimage.binary_opening(data, \
            structure=disk_structure(n)).sum() for n in sizes]
    return granulo

def thresholding(img, thresh, size=9):
    """
    Segment using a thresholding algorithm
    
    Input:
     - img  ndarray : Image array (ndim=2)
     - thresh float : Threshold value for pixels selectino
     - size     int : Minimum size a group of pixels must have
    
    Output:
     - regions : Binary array for each segmented region
    
    ---
    """

    logging.debug("Threshold: %.2f", thresh)
    logging.debug("Objects min size: %d", size)

    # Take the binary image thresholded
    img_bin = img > thresh

    # And use (MO) binary opening (erosion + dilation) for cleaning spurious Trues
    strct = ndi.generate_binary_structure(2, 2)
    img_bin = ndi.binary_opening(img_bin, strct)

    # Label each group/region (value==True) of pixels
    regions, nlbl = ndi.label(img_bin)
    for i in xrange(1, nlbl + 1):
        inds = np.where(regions == i)
        if inds[0].size < size:
            regions[inds] = 0

    logging.debug("Threshold labels: %s", np.unique(regions))

    return regions.astype(np.bool)

def find_feature_mask(s_msk, sigma=1, ax=None, x_values=None):
    # find emission features from observed spec.

    filtered_spec = s_msk - ni.median_filter(s_msk, 15)
    filtered_spec = ni.gaussian_filter1d(filtered_spec, 0.5)
    smoothed_std = get_smoothed_std(filtered_spec,
                                    rad=3, smooth_length=3)
    emission_feature_msk_ = filtered_spec > sigma*smoothed_std
    #emission_feature_msk_ = ni.binary_closing(emission_feature_msk_)
    emission_feature_msk = ni.binary_opening(emission_feature_msk_,
                                             iterations=1)

    if ax is not None:
        if x_values is None:
            x_values = np.arange(len(s_msk))

        #ax.plot(x_values, s_msk)
        ax.plot(x_values, filtered_spec)
        ax.plot(x_values, smoothed_std)

        ax.plot(x_values[emission_feature_msk],
                emission_feature_msk[emission_feature_msk],
                "ys", mec="none")

    return emission_feature_msk