Python ndimage.generate_binary_structure函数代码示例

def adjust_spot_positions(image, label_image, hp, debug=None):
    """Re-evaluate the spot positions based on the segmentation. 
        Parameters: 
        image: The original image (can be masked) that was sent to findspot3d
        label_image: the label image containing two labels 
        hp: the original hotpoints
        debug: set to true to write out an image debugimg.nii.gz with the stuff
        """
        
    struct2 = generate_binary_structure(3, 2)
    struct1 = generate_binary_structure(3, 1)
    peak_points =[] 

    if debug is None:
        temp_path = os.getenv("PYSBR_TEMP")
        if temp_path is not None:
            debug = os.path.join(temp_path, "debug-labels.nii.gz")
    
    if debug is not None:
            debimg = image.copy()

    nlabels = label_image.max()

    if nlabels!=len(hp):
        raise RuntimeError( 'number of labels and hotspots should be the same' )

    tins = []
    for n in range(nlabels):
        label = n+1
        area = binary_closing(label_image == label, struct2)
        thiniter = np.sum(area.reshape(-1)) / 1500 + 1
        csbr.thinning3d(area, thiniter)
        tins.append(area)
   
    for n in range(nlabels):
        label = n+1
        
        #avoid that a single pixel breaks the evaluation by running a closing 
        area = label_image == label
        
        #evaluate the boundary 
        dmask = binary_dilation(area, struct1)
        border = np.bitwise_xor(dmask, area)
        
        p = adjust_spot_position(image, border, image[tuple(hp[n])], tins[n], tins[(n + 1) % 2])
        peak_points.append(p)

        if debug is not None:
            debimg[border>0] = 196
            debimg[p] = 0
            nib.save(nib.Nifti1Image(debimg, global_affine), debug)

    peak_points = np.array( peak_points )
    return peak_points

def getBoundariesOfimage(image):
    """
    find edges by using erosion
    """
    if np.ndim(image) == 2:
        sElement = ndimage.generate_binary_structure(2, 1)
    else:
        sElement = ndimage.generate_binary_structure(3, 1)
    erode_im = scipy.ndimage.morphology.binary_erosion(image, sElement)
    b = image - erode_im
    return b

 def InDecPatch(self,which,amount):
     s = ndimage.generate_binary_structure(2,1) # taxi-cab struct
     if which == 0:
         ras = ndimage.binary_dilation(self.cl_array,s,iterations=amount,border_value=0)
     else:
         ras = ndimage.binary_erosion(self.cl_array,s,iterations=amount,border_value=0)
     return(ras)

def artifact_mask(imdata, airdata, distance, zscore=10.):
    """Computes a mask of artifacts found in the air region"""
    from statsmodels.robust.scale import mad

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    if np.sum((bg_img > 0).astype(np.uint8)) < 100:
        return np.zeros_like(airdata)

    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    bg_location = np.median(bg_img[bg_img > 0])
    bg_spread = mad(bg_img[bg_img > 0])
    bg_img[bg_img > 0] -= bg_location
    bg_img[bg_img > 0] /= bg_spread

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > zscore] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

def manual_split(probs, seg, body, seeds, connectivity=1, boundary_seeds=None):
    """Manually split a body from a segmentation using seeded watershed.

    Input:
        - probs: the probability of boundary in the volume given.
        - seg: the current segmentation.
        - body: the label to be split.
        - seeds: the seeds for the splitting (should be just two labels).
        [-connectivity: the connectivity to use for watershed.]
        [-boundary_seeds: if not None, these locations become inf in probs.]
    Value:
        - the segmentation with the selected body split.
    """
    struct = generate_binary_structure(seg.ndim, connectivity)
    body_pixels = seg == body
    bbox = find_objects(body_pixels)[0]
    body_pixels = body_pixels[bbox]
    body_boundary = binary_dilation(body_pixels, struct) - body_pixels
    non_body_pixels = True - body_pixels - body_boundary
    probs = probs.copy()[bbox]
    probs[non_body_pixels] = probs.min()-1
    if boundary_seeds is not None:
        probs[boundary_seeds[bbox]] = probs.max()+1
    probs[body_boundary] = probs.max()+1
    seeds = label(seeds.astype(bool)[bbox], struct)[0]
    outer_seed = seeds.max()+1 # should be 3
    seeds[non_body_pixels] = outer_seed
    seg_new = watershed(probs, seeds, 
        dams=(seg==0).any(), connectivity=connectivity, show_progress=True)
    seg = seg.copy()
    new_seeds = unique(seeds)[:-1]
    for new_seed, new_label in zip(new_seeds, [0, body, seg.max()+1]):
        seg[bbox][seg_new == new_seed] = new_label
    return seg

    def __init__(self, label_image=None, connectivity=1, data=None, **attr):

        super(RAG, self).__init__(data, **attr)
        if self.number_of_nodes() == 0:
            self.max_id = 0
        else:
            self.max_id = max(self.nodes_iter())

        if label_image is not None:
            fp = ndi.generate_binary_structure(label_image.ndim, connectivity)
            # In the next ``ndi.generic_filter`` function, the kwarg
            # ``output`` is used to provide a strided array with a single
            # 64-bit floating point number, to which the function repeatedly
            # writes. This is done because even if we don't care about the
            # output, without this, a float array of the same shape as the
            # input image will be created and that could be expensive in
            # memory consumption.
            ndi.generic_filter(
                label_image,
                function=_add_edge_filter,
                footprint=fp,
                mode='nearest',
                output=as_strided(np.empty((1,), dtype=np.float_),
                                  shape=label_image.shape,
                                  strides=((0,) * label_image.ndim)),
                extra_arguments=(self,))

 def f_returnInternalEdge(self,cl_array):
     # Internal edge: Count of neighboring non-zero cell       
     kernel = ndimage.generate_binary_structure(2, 1) # Make a kernel
     kernel[1, 1] = 0
     b = ndimage.convolve(cl_array, kernel, mode="constant")
     n_interior = b[cl_array != 0].sum() # Number of interiror edges
     return n_interior

def edge_matrix(labels, connectivity=1):
    """Generate a COO matrix containing the coordinates of edge pixels.

    Parameters
    ----------
    labels : array of int
        An array of labeled pixels (or voxels).
    connectivity : int in {1, ..., labels.ndim}
        The square connectivity for considering neighborhood.

    Returns
    -------
    edges : sparse.coo_matrix
        A COO matrix where (i, j) indicate neighboring labels and the
        corresponding data element is the linear index of the edge pixel
        in the labels array.
    """
    conn = ndi.generate_binary_structure(labels.ndim, connectivity)
    eroded = ndi.grey_erosion(labels, footprint=conn).ravel()
    dilated = ndi.grey_dilation(labels, footprint=conn).ravel()
    labels = labels.ravel()
    boundaries0 = np.flatnonzero(eroded != labels)
    boundaries1 = np.flatnonzero(dilated != labels)
    labels_small = np.concatenate((eroded[boundaries0], labels[boundaries1]))
    labels_large = np.concatenate((labels[boundaries0], dilated[boundaries1]))
    n = np.max(labels_large) + 1
    data = np.concatenate((boundaries0, boundaries1))
    sparse_graph = sparse.coo_matrix((data, (labels_small, labels_large)),
                                     dtype=np.int_, shape=(n, n))
    return sparse_graph

def objextract(Fg):

    s = nd.generate_binary_structure(2,2)
    labeled_array, num_features = nd.measurements.label(Fg, structure=s)
    coor = []
    cnt = []

    if num_features == 0:
       idx = []
    else:    
        lth = 200   # label pixel number less than lth will be removed
        Lth = 6500
        for i in range(1,num_features+1):
            coor.append(np.where(labeled_array==i))
            cnt.append(len(np.where(labeled_array==i)[1]))

        cnt = array(cnt)
        idx = arange(num_features)
        idx = idx[(cnt<Lth)&(cnt>lth)]
      
        if len(idx)==0:
            idx = []
        elif len(idx)>1:              
            #idx = [idx[cnt[idx].argmax()]]
            idx = sorted(range(len(cnt)),key=lambda x:cnt[x])[::-1][0:2]

    return idx,labeled_array,coor,cnt

def find_local_max(img, d_rad, threshold=1e-15):
    """
    This is effectively a replacement for pkfnd in the matlab/IDL code.

    The output of this function is meant to be feed into :py:func:`~subpixel_centroid`

    The magic of numpy means this should work for any dimension data.

    :param img: an ndarray representing the data to find the local maxes
    :param d_rad: the radius of the dilation, the smallest possible spacing between local maximum
    :param threshold: optional, voxels < threshold are ignored.

    :rtype: (d,N) array of the local maximums.
    """
    d_rad = int(d_rad)
    img = np.array(np.squeeze(img))       # knock out singleton dimensions
    img[img < threshold] = -np.inf        # mask out pixels below threshold
    dim = img.ndim                        # get the dimension of data

    # make structuring element
    s = ndimage.generate_binary_structure(dim, 1)
    # scale it up to the desired size
    d_struct = ndimage.iterate_structure(s, int(d_rad))
    dilated_img = ndimage.grey_dilation(img,
                                        footprint=d_struct,
                                        cval=0,
                                        mode='constant')   # do the dilation

    # find the locations that are the local maximum
    # TODO clean this up
    local_max = np.where(np.exp(img - dilated_img) > (1 - 1e-15))
    # the extra [::-1] is because matplotlib and ndimage disagree an xy vs yx
    return np.vstack(local_max[::-1])

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 remove_small_objects(ar, min_size=64, connectivity=1, in_place=False):
    # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7
    if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)):
        raise TypeError("Only bool or integer image types are supported. "
                        "Got %s." % ar.dtype)

    if in_place:
        out = ar
    else:
        out = ar.copy()

    if min_size == 0:  # shortcut for efficiency
        return out

    if out.dtype == bool:
        selem = nd.generate_binary_structure(ar.ndim, connectivity)
        ccs = np.zeros_like(ar, dtype=np.int32)
        nd.label(ar, selem, output=ccs)
    else:
        ccs = out

    try:
        component_sizes = np.bincount(ccs.ravel())
    except ValueError:
        raise ValueError("Negative value labels are not supported. Try "
                         "relabeling the input with `scipy.ndimage.label` or "
                         "`skimage.morphology.label`.")

    too_small = component_sizes < min_size
    too_small_mask = too_small[ccs]
    out[too_small_mask] = 0

    return out

def artifact_mask(imdata, airdata, distance):
    """Computes a mask of artifacts found in the air region"""
    import nibabel as nb

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    # CHANGED - to the 75 percentile
    bg_threshold = np.percentile(bg_img[airdata > 0], 75)

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > bg_threshold] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

def compute_sparsity(im):
    l_x = len(im)
    X, Y = np.ogrid[:l_x, :l_x]
    mask = ((X - l_x/2)**2 + (Y - l_x/2)**2 <= (l_x/2)**2)
    grad1 = ndimage.morphological_gradient(im, footprint=np.ones((3, 3)))
    grad2 = ndimage.morphological_gradient(im, footprint=ndimage.generate_binary_structure(2, 1))
    return (grad1[mask] > 0).mean(), (grad2[mask] > 0).mean() 

def ClusterizeImage(image,thresh=None,connectivity=3):
    if thresh is None:
        thresh = 0
    image[np.where(image<=thresh)] = 0
    s = generate_binary_structure(3,connectivity)
    larray, nf = label(image,s)
    return larray