Python scipy.ndimage.morphology.binary_erosion() Examples

def decouple_volumes(v1, v2, mode, se=None, iterations=1):
    """
    
    mode : {inner-from-outer, outer-from-inner, neighbors}
        inner-from-outer: this changes v1 by removing voxels
        outer-from-inner: this changes v2 by adding voxels
        neighbors: this changes v2 by removing voxels
    
    """
    assert mode in ["inner-from-outer","outer-from-inner","neighbors"]
    
    if isinstance(v1, str) and os.path.isfile(v1):
        v1 = nib.load(v1)
    assert isinstance(v1, nib.Nifti1Image) or isinstance(v1, nib.Nifti2Image)
    d1 = v1.get_data()
    if isinstance(v2, str) and os.path.isfile(v2):
        v2 = nib.load(v2)
    assert isinstance(v2, nib.Nifti1Image) or isinstance(v2, nib.Nifti2Image)
    d2 = v2.get_data()
    
    assert d1.ndim is d2.ndim
    
    
    if se is None:
        se = mrph.generate_binary_structure(d1.ndim,1)
    
    if mode == "inner-from-outer":
        # make v2/d2 the inner volume
        d1, d2 = d2, d1
        v1, v2 = v2, v1        
        d2 = d2 & mrph.binary_erosion(d1, se, iterations)
        
    if mode == "outer-from-inner":
        d2 = d2 | mrph.binary_dilation(d1, se, iterations)
        
    if mode == "neighbors":
        d2 = d2 & ~mrph.binary_dilation(d1, se, iterations)
    
    d2 = nib.Nifti1Image(d2, v2.affine, header=v2.header)
    d2.set_filename(v2.get_filename())
    return d2 

def surface_distances(y_true, y_pred, voxel_shape=None):
    check_bool(y_pred, y_true)
    check_shapes(y_pred, y_true)

    pred_border = np.logical_xor(y_pred, binary_erosion(y_pred))
    true_border = np.logical_xor(y_true, binary_erosion(y_true))

    dt = distance_transform_edt(~true_border, sampling=voxel_shape)
    return dt[pred_border] 

def local_maxima(arr):
    # http://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array/3689710#3689710
    """
    Takes an array and detects the troughs using the local maximum filter.
    Returns a boolean mask of the troughs (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """
    # define an connected neighborhood
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#generate_binary_structure
    neighborhood = morphology.generate_binary_structure(len(arr.shape),2)
    # apply the local maximum filter; all locations of maximal value 
    # in their neighborhood are set to 1
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.filters.html#maximum_filter
    local_max = (filters.maximum_filter(arr, footprint=neighborhood)==arr)
    # local_max is a mask that contains the peaks we are 
    # looking for, but also the background.
    # In order to isolate the peaks we must remove the background from the mask.
    # 
    # we create the mask of the background
    background = (arr==arr.min())           # mxu: in the original version, was         background = (arr==0)
    # 
    # a little technicality: we must erode the background in order to 
    # successfully subtract it from local_max, otherwise a line will 
    # appear along the background border (artifact of the local maximum filter)
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#binary_erosion
    eroded_background = morphology.binary_erosion(
        background, structure=neighborhood, border_value=1)
    # 
    # we obtain the final mask, containing only peaks, 
    # by removing the background from the local_max mask
    #detected_maxima = local_max - eroded_backround             # mxu: this is the old version, but the boolean minus operator is deprecated
    detected_maxima = np.bitwise_and(local_max, np.bitwise_not(eroded_background))          # Material nonimplication, see http://en.wikipedia.org/wiki/Material_nonimplication
    return np.where(detected_maxima) 

def calcImaMaskBrd(self):
        """Calculate borders of image mask slice."""
        return self.imaSlcMsk - binary_erosion(self.imaSlcMsk) 

def find_edges(self, member='total'):
        self.find_footprint(member=member)

        edges = morphology.binary_erosion(self.footprint).astype(np.int16)
        self.edges = self.footprint - edges 

def detect(self, image):
        # define an 8-connected neighborhood
        neighborhood = generate_binary_structure(2, 2)

        # apply the local maximum filter; all pixel of maximal value
        # in their neighborhood are set to 1
        local_max = maximum_filter(image, footprint=neighborhood) == image
        # local_max is a mask that contains the peaks we are
        # looking for, but also the background.
        # In order to isolate the peaks we must remove the background from the mask.

        # we create the mask of the background
        background = (image < self.min_th)

        # a little technicality: we must erode the background in order to
        # successfully subtract it form local_max, otherwise a line will
        # appear along the background border (artifact of the local maximum filter)
        eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)

        # we obtain the final mask, containing only peaks,
        # by removing the background from the local_max mask (xor operation)
        detected_peaks = local_max ^ eroded_background

        detected_peaks[image < self.min_th] = False
        peaks = np.array(np.nonzero(detected_peaks)).T

        if len(peaks) == 0:
            return peaks, np.array([])

        # nms
        if len(peaks) == 1:
            clusters = [0]
        else:
            clusters = fclusterdata(peaks, self.min_dist, criterion="distance")
        peak_groups = {}
        for ind_junc, ind_group in enumerate(clusters):
            if ind_group not in peak_groups.keys():
                peak_groups[ind_group] = []
                peak_groups[ind_group].append(peaks[ind_junc])
        peaks_nms = []
        peaks_score = []
        for peak_group in peak_groups.values():
            values = [image[y, x] for y, x in peak_group]
            ind_max = np.argmax(values)
            peaks_nms.append(peak_group[int(ind_max)])
            peaks_score.append(values[int(ind_max)])

        return np.float32(np.array(peaks_nms)), np.float32(np.array(peaks_score)) 

def distance_to_surface(structure_mask, target_mask, voxel_pts):
    """For every voxel in structure, find minimum distance to target surface.
    The target surface is calculated by removing the binary erosion of the target mask.
    
    Args:
        structure_mask: a 3D numpy array of boolean values
        target_mask: a 3D numpy array of boolean values
        voxel_pts: a 4D numpy array of voxel location points (x,y,z) for each voxel in masks

    Returns:
        a 3D numpy array of minimum distances to target_mask surface for each True voxel in structure_mask, 
        numpy.NaN otherwise.
    """
    
    # TODO: make multitreaded

    # compute surface voxels
    target_shell = np.logical_and(target_mask,np.logical_not(binary_erosion(target_mask,iterations=1)))        
    surface_pts = voxel_pts[np.where(target_shell)]

    # index voxels inside structure but outside target
    idx_outside = np.where(np.logical_and(structure_mask,np.logical_not(target_mask)))
    
    # index voxels inside structure and inside target
    idx_inside = np.where(np.logical_and(structure_mask,target_mask))
    
    structure_pts_outside = voxel_pts[idx_outside]
    structure_pts_inside = voxel_pts[idx_inside]
    
    dists_to_target = np.empty(structure_mask.shape) 
    dists_to_target[:] = np.NaN

    # only compute if there are structure voxels ouside the target
    if len(structure_pts_outside.shape) != 0:
        dists_to_target[idx_outside] = list(map(lambda pt: _min_dist(pt,surface_pts),structure_pts_outside))    
    
    # only compute if there are structure voxels inside the target
    if len(structure_pts_inside) != 0:
        # pts inside have negative dist
        dists_to_target[idx_inside] = list(map(lambda pt: - _min_dist(pt,surface_pts),structure_pts_inside))
    
    return dists_to_target