Python scipy.ndimage 模块，distance_transform_edt() 实例源码

def d_H(im1, im2, lab, return_mm3):
    """
    Asymmetric component of the Hausdorff distance.
    :param im1: first image
    :param im2: second image
    :param lab: label in the image
    :param return_mm3: final unit of measures of the result.
    :return: max(d(x, contourY)), x: point belonging to the first contour,
                                  contourY: contour of the second segmentation.
    """
    arr1 = im1.get_data() == lab
    arr2 = im2.get_data() == lab
    if np.count_nonzero(arr1) == 0 or np.count_nonzero(arr2) == 0:
        return np.nan
    if return_mm3:
        dt2 = nd.distance_transform_edt(1 - arr2, sampling=list(np.diag(im1.affine[:3, :3])))
    else:
        dt2 = nd.distance_transform_edt(1 - arr2, sampling=None)
    return np.max(dt2 * arr1)

def calculate_distance(centers_of_mass, image):
    """
    takes the centers of each blob, and an image to be segmented. Divides the image according to the center of masses
    by a random walk

    :param image: a binarised image to be segmented
    :param centers_of_mass: the centres that will define the maxima of the watershed segmentation
    :return segmentation_labels: a labelled image/segmentation, where each index belongs do a different center of mass

    """
    # random walk segmentation of 2D image-mask array
    distance = ndimage.distance_transform_edt(np.abs(image-1))
    local_maxi = np.zeros_like(image)

    for c in centers_of_mass:
        local_maxi[int(c[0]), int(c[1])] = 1

    markers = ndimage.label(local_maxi)[0]
    segmentation_labels = segmentation.random_walker(distance, markers, beta=60)

    return segmentation_labels

def blob_labels(centers_of_mass, blob_image):
    """
    label nuclei with segmentation - so labels are in the same order as the outer layer

    :param list centers_of_mass: centers of target blobs/cells
    :param np.array blob_image: image of the target cells, or nuclei of cells
    :return segmented_blobs: a labelled image where each index is a different cell
    :return distance: image where each pixel's value is related to how far away from a blob center it is
    """

    image = np.abs(blob_image-1)
    distance = ndimage.distance_transform_edt(np.abs(image-1))
    local_maxi = np.zeros_like(image)
    for c in centers_of_mass:
        local_maxi[int(c[0]), int(c[1])] = 1
    markers = ndimage.label(local_maxi)[0]
    segmented_blobs = segmentation.random_walker(distance, markers, beta=20)

    return segmented_blobs, distance

def __init__(self, map_msg):
        self.memoize_disks = True
        self.map_msg = map_msg
        self.map_info = map_msg.info
        #  # 0: permissible, -1: unmapped, 100: blocked
        self.raw_matrix = np.array(map_msg.data).reshape((map_msg.info.height, map_msg.info.width))
        self.exploration_coeff = float(rospy.get_param("~exploration_coeff", 0.75))

        # reversed from expectation since this is what distance_transform_edt requires
        self.occupancy_grid = np.ones_like(self.raw_matrix, dtype=bool)
        self.occupancy_grid[self.raw_matrix>50] = 0

        # 0: not permissible, 1: permissible
        self.permissible_region = np.zeros_like(self.raw_matrix, dtype=bool)
        self.permissible_region[self.raw_matrix==0] = 1

        self.distmap = ndimage.distance_transform_edt(self.occupancy_grid)
        self.exploration_buffer = np.zeros_like(self.raw_matrix, dtype=bool)

        if self.memoize_disks:
            self.memo_table = {}
            self.memoized = 0
            self.unmemoized = 0

def distancetransf(image):
    if image.dtype=='bool':
        return ndimage.distance_transform_edt(1-image.astype(np.int8))
        else:
        return ndimage.distance_transform_edt(1-image)

def _neighbor_distances(ndim):
    center = np.ones((3,) * ndim, dtype=bool)
    center.ravel()[center.size//2] = False
    out = ndi.distance_transform_edt(center).ravel()
    return out

def one_country(maparray, country, background):
    cmask = mask_by_color(maparray, colorfn(country))
    edtmask = cmask
    for c in background:
        edtmask = edtmask | mask_by_color(maparray, colorfn(c))
    _, (indx,indy) = ndimage.distance_transform_edt(edtmask, return_indices=True)
    edtvals = maparray[indx,indy]
    for i in range(edtvals.shape[-1]): edtvals[...,i] *= cmask
    return Image.fromarray(edtvals)

def bwdist(a):
    """
    Intermediary function. 'a' has only True/False vals,
    so we convert them into 0/1 values - in reverse.
    True is 0, False is 1, distance_transform_edt wants it that way.
    """
    return nd.distance_transform_edt(a == 0)


# Displays the image with curve superimposed

def mask_to_dist(mask):
  """Returns distance transform on a mask

  Arguments:
    mask (array): the mask for which to calculate distance transform

  Returns:
    array: distance transform of the mask
  """
  return nd.distance_transform_edt(mask == 0)

def main(args):
    if args.threshold is None:
        print("Please provide a binarization threshold")
        return 1
    data, hdr = read(args.input, inc_header=True)
    mask = data >= args.threshold
    if args.minvol is not None:
        mask = binary_volume_opening(mask, args.minvol)
    if args.fill:
        mask = binary_fill_holes(mask)
    if args.extend is not None and args.extend > 0:
        if args.relion:
            se = binary_sphere(args.extend, False)
            mask = binary_dilation(mask, structure=se, iterations=1)
        else:
            dt = distance_transform_edt(~mask)
            mask = mask | (dt <= args.edge_width)
    if args.close:
        se = binary_sphere(args.extend, False)
        mask = binary_closing(mask, structure=se, iterations=1)
    final = mask.astype(np.single)
    if args.edge_width is not None:
        dt = distance_transform_edt(~mask)  # Compute *outward* distance transform of mask.
        idx = (dt <= args.edge_width) & (dt > 0)  # Identify edge points by distance from mask.
        x = np.arange(1, args.edge_width + 1)  # Domain of the edge profile.
        if "sin" in args.edge_profile:
            y = np.sin(np.linspace(np.pi/2, 0, args.edge_width + 1))  # Range of the edge profile.
        f = interp1d(x, y[1:])
        final[idx] = f(dt[idx])  # Insert edge heights interpolated at distance transform values.
    write(args.output, final, psz=hdr["xlen"] / hdr["nx"])
    return 0

def create_roi_masks(centers_of_mass, putative_nuclei_image, putative_somata_image=None, radius=3):
        """
        create roi masks for the outer segment of the cell (i.e. soma)

        :param radius: limits the size of the mask
        :return roi_mask_list: a list of pixels for each cell, for further analysis
        """
        roi_mask_list = []

        if putative_somata_image is None:
            putative_somata_image = np.zeros_like(putative_nuclei_image)
        putative_nuclei_image = remove_small_blobs(centers_of_mass, putative_nuclei_image)[0]
        watershed_image = np.logical_or(putative_nuclei_image, putative_somata_image)
        labelled_watershed = calculate_distance(centers_of_mass, watershed_image)
        labelled_putative_somata = putative_somata_image*labelled_watershed
        labelled_putative_nuclei = calculate_distance(centers_of_mass, putative_nuclei_image)*putative_nuclei_image  # nuclei need their own watershed

        for i in range(np.max(labelled_putative_somata)):  # for each nucleus

            # calculate the distance away from the nucleus boundary

            distance_from_blob_centre = ndimage.distance_transform_edt(labelled_putative_nuclei != i+1)
            bool_mask = np.ones_like(labelled_putative_nuclei)
            bool_mask[distance_from_blob_centre > radius] = 0
            bool_mask[distance_from_blob_centre == 0] = 0

            # take all indices within the radius number of pixels of the nucleus boundary

            idx = np.where(labelled_putative_somata*bool_mask == i+1)
            x = np.vstack((idx[0], idx[1]))
            m = np.reshape(x.T, (len(idx[0]), 2))
            roi_mask_list.append(m)

        return roi_mask_list

def mid_axis(img):
    dis = ndimg.distance_transform_edt(img)
    idx = np.argsort(dis.flat).astype(np.int32)
    medial_axis(dis, idx, lut)
    return dis

def run(self, ips, snap, img, para = None):
        return ndimg.distance_transform_edt(snap)

def run(self, ips, snap, img, para = None):
        img[:] = snap
        dist = -ndimg.distance_transform_edt(snap)
        pts = find_maximum(dist, para['tor'], False)
        buf = np.zeros(ips.size, dtype=np.uint16)
        buf[pts[:,0], pts[:,1]] = 1
        markers, n = ndimg.label(buf, np.ones((3,3)))
        line = watershed(dist, markers, line=True, conn=para['con']+1)
        img[line==0] = 0

def run(self, ips, snap, img, para = None):
        dist = ndimg.distance_transform_edt(snap)
        markers, n = ndimg.label(snap==0, np.ones((3,3)))

        line = watershed(dist, markers)
        if para['type']=='segment with ori':
            img[:] = np.where(line==0, 0, snap)
        if para['type']=='segment only':
            img[:] = (line>0) * 255
        if para['type']=='white line':
            img[:] = (line==0) * 255
        if para['type']=='gray line':
            img[:] = np.where(line==0, dist, 0)

def extend_to_global(var, src_grid, global_grid, arctic_filler=None):
    """
    Use nearest neighbour to extend obs/source over the whole globe, including land.
    """

    # Create masked array of the correct shape.
    tmp = np.zeros((global_grid.num_levels, global_grid.num_lat_points,
                    global_grid.num_lon_points))
    new_data = np.ma.array(tmp, mask=global_grid.mask, copy=True)
    # Mask everything by default. The code below fills masked values with
    # nearest neighbour.
    new_data.mask[:] = True

    # Drop obs data into new grid at correct location
    lat_min_idx = find_nearest_index(global_grid.y_t[:, 0], np.min(src_grid.y_t[:]))
    if np.max(global_grid.y_t[:]) <= np.max(src_grid.y_t[:]):
        new_data[:, lat_min_idx:, :] = var[:]
    else:
        lat_max_idx = find_nearest_index(global_grid.y_t[:, 0], np.max(src_grid.y_t[:]))
        new_data[:, lat_min_idx:lat_max_idx+1, :] = var[:]

    # Fill in missing values on each level with nearest neighbour
    for l in range(var.shape[0]):
        ind = nd.distance_transform_edt(new_data[l, :, :].mask,
                                        return_distances=False,
                                        return_indices=True)
        tmp = new_data[l, :, :]
        tmp = tmp[tuple(ind)]
        new_data[l, :, :] = tmp[:, :]

    return new_data

def symmetric_contour_distance_l(im1, im2, lab, return_mm3, formula='normalised'):
    """
    Generalised normalised symmetric contour distance.
    On the set {d(x, contourY)) | x in contourX}, several statistics can be computed.
     Mean, median and standard deviation can be useful, as well as a more robust normalisation
     Formula can be
    :param im1:
    :param im2:
    :param lab:
    :param return_mm3:
    :param formula: 'normalised', 'averaged', 'median', 'std'
    'normalised' = (\sum_{x in contourX} d(x, contourY))  + \sum_{y in contourY} d(y, contourX))) / (|contourX| + |contourY|)
    'averaged'   = 0.5 (mean({d(x, contourY)) | x in contourX}) + mean({d(y, contourX)) | y in contourY}))
    'median'     = 0.5 (median({d(x, contourY)) | x in contourX}) + median({d(y, contourX)) | y in contourY}))
    'std'        = 0.5 \sqrt(std({d(x, contourY)) | x in contourX})^2 + std({d(y, contourX)) | y in contourY})^2)
    :return:
    """
    arr1 = im1.get_data() == lab
    arr2 = im2.get_data() == lab

    if np.count_nonzero(arr1) == 0 or np.count_nonzero(arr2) == 0:
        return np.nan

    arr1_contour = contour_from_array_at_label_l(arr1, 1)
    arr2_contour = contour_from_array_at_label_l(arr2, 1)

    if return_mm3:
        dtb1 = nd.distance_transform_edt(1 - arr1_contour, sampling=list(np.diag(im1.affine[:3, :3])))
        dtb2 = nd.distance_transform_edt(1 - arr2_contour, sampling=list(np.diag(im1.affine[:3, :3])))
    else:
        dtb1 = nd.distance_transform_edt(1 - arr1_contour)
        dtb2 = nd.distance_transform_edt(1 - arr2_contour)

    dist_border1_array2 = arr2_contour * dtb1
    dist_border2_array1 = arr1_contour * dtb2

    if formula == 'normalised':
        return (np.sum(dist_border1_array2) + np.sum(dist_border2_array1)) / (np.count_nonzero(arr1_contour) + np.count_nonzero(arr2_contour))
    elif formula == 'averaged':
        return .5 * (np.mean(dist_border1_array2) + np.mean(dist_border2_array1))
    elif formula == 'median':
        return .5 * (np.median(dist_border1_array2) + np.median(dist_border2_array1))
    elif formula == 'std':
        return np.sqrt( .5 * (np.std(dist_border1_array2)**2 + np.std(dist_border2_array1)**2))
    elif formula == 'average_std':
        return .5 * (np.mean(dist_border1_array2) + np.mean(dist_border2_array1)), \
               np.sqrt(.5 * (np.std(dist_border1_array2) ** 2 + np.std(dist_border2_array1) ** 2))
    else:
        raise IOError


# --- distances - (segm, segm) |-> pandas.Series (indexed by labels)