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)
python类distance_transform_edt()的实例源码
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)
