def get_subvolume(self, bounds):
subvol_shape = bounds.stop - bounds.start
label_shape = subvol_shape - 2 * bounds.label_margin
parent_bounds = SubvolumeBounds(self.world_coord_to_parent(bounds.start),
self.world_coord_to_parent(bounds.stop),
label_margin=self.world_coord_to_parent(bounds.label_margin))
subvol = self.parent.get_subvolume(parent_bounds)
subvol.image = subvol.image.reshape(
[subvol_shape[0], self.scale[0],
subvol_shape[1], self.scale[1],
subvol_shape[2], self.scale[2]]).mean(5).mean(3).mean(1)
if subvol.label_mask is not None:
# Downsample body mask by considering blocks where the majority
# of voxels are in the body to be in the body. Alternatives are:
# - Conjunction (tends to introduce false splits)
# - Disjunction (tends to overdilate and merge)
# - Mode label (computationally expensive)
if CONFIG.volume.label_downsampling == 'conjunction':
subvol.label_mask = subvol.label_mask.reshape(
[label_shape[0], self.scale[0],
label_shape[1], self.scale[1],
label_shape[2], self.scale[2]]).all(5).all(3).all(1)
else:
subvol.label_mask = subvol.label_mask.reshape(
[label_shape[0], self.scale[0],
label_shape[1], self.scale[1],
label_shape[2], self.scale[2]]).mean(5).mean(3).mean(1) > 0.5
# Note that this is not a coordinate xform to parent in the typical
# sense, just a rescaling of the coordinate in the subvolume-local
# coordinates. Hence no similar call in VolumeView.get_subvolume.
subvol.seed = self.parent_coord_to_world(subvol.seed)
return subvol
python类label()的实例源码
def get_largest_component(self, closing_shape=None):
mask, bounds = self._get_bounded_mask(closing_shape)
label_im, num_labels = ndimage.label(mask)
label_sizes = ndimage.sum(mask, label_im, range(num_labels + 1))
label_im[(label_sizes < label_sizes.max())[label_im]] = 0
label_im = np.minimum(label_im, 1)
if label_im[tuple(self.seed - bounds[0])] == 0:
logging.warning('Seed voxel ({}) is not in connected component.'.format(np.array_str(self.seed)))
return label_im, bounds
def get_seeded_component(self, closing_shape=None):
mask, bounds = self._get_bounded_mask(closing_shape)
label_im, _ = ndimage.label(mask)
seed_label = label_im[tuple(self.seed - bounds[0])]
if seed_label == 0:
raise ValueError('Seed voxel (%s) is not in body.', np.array_str(self.seed))
label_im[label_im != seed_label] = 0
label_im[label_im == seed_label] = 1
return label_im, bounds
def compute_centroids(image):
"""Find the centroids of all nonzero connected blobs in `image`.
Parameters
----------
image : ndarray
The input image.
Returns
-------
label_image : ndarray of int
The input image, with each connected region containing a different
integer label.
Examples
--------
>>> image = np.array([[1, 0, 1, 0, 0, 1, 1],
... [1, 0, 0, 1, 0, 0, 0]])
>>> labels, centroids = compute_centroids(image)
>>> print(labels)
[[1 0 2 0 0 3 3]
[1 0 0 2 0 0 0]]
>>> centroids
array([[ 0.5, 0. ],
[ 0.5, 2.5],
[ 0. , 5.5]])
"""
connectivity = np.ones((3,) * image.ndim)
labeled_image = ndi.label(image, connectivity)[0]
nz = np.nonzero(labeled_image)
nzpix = labeled_image[nz]
sizes = np.bincount(nzpix)
coords = np.transpose(nz)
grouping = np.argsort(nzpix)
sums = np.add.reduceat(coords[grouping], np.cumsum(sizes)[:-1])
means = sums / sizes[1:, np.newaxis]
return labeled_image, means
def mesh_sizes(skeleton):
"""Compute the area in pixels of the spaces between skeleton branches.
This only makes sense for 2D images.
Parameters
----------
skeleton : array, shape (M, N)
An image containing a single-pixel-wide closed skeleton.
Returns
-------
sizes : array of int, shape (P,)
The sizes of all spaces delineated by the skeleton *not* touching
the borders.
Examples
--------
>>> image = np.array([[0, 0, 1, 0, 0],
... [0, 0, 1, 1, 1],
... [0, 0, 1, 0, 0],
... [0, 1, 0, 1, 0]])
>>> print(mesh_sizes(image))
[]
>>> from skan.nputil import pad
>>> image2 = pad(image, 1) # make sure mesh not touching border
>>> print(mesh_sizes(image2)) # sizes in row order of first pixel in space
[7 2 3 1]
"""
spaces = ~skeleton.astype(bool)
labeled = ndi.label(spaces)[0]
touching_border = np.unique(np.concatenate((labeled[0], labeled[-1],
labeled[:, 0],
labeled[:, -1])))
sizes = np.bincount(labeled.flat)
sizes[touching_border] = 0
sizes = sizes[sizes != 0]
return sizes
def image_summary(skeleton, *, spacing=1):
"""Compute some summary statistics for an image.
Parameters
----------
skeleton : array, shape (M, N)
The input image.
Other Parameters
----------------
spacing : float or array of float, shape (`skeleton.ndim`,)
The resolution along each axis of `skeleton`.
Returns
-------
stats : pandas.DataFrame
Selected statistics about the image.
"""
stats = pd.DataFrame()
stats['scale'] = [spacing]
g, coords, degimg = csr.skeleton_to_csgraph(skeleton, spacing=spacing)
degrees = np.diff(g.indptr)
num_junctions = np.sum(degrees > 2)
stats['number of junctions'] = num_junctions
pixel_area = (spacing ** skeleton.ndim if np.isscalar(spacing) else
np.prod(spacing))
stats['area'] = np.prod(skeleton.shape) * pixel_area
stats['junctions per unit area'] = (stats['number of junctions'] /
stats['area'])
sizes = mesh_sizes(skeleton)
stats['average mesh area'] = np.mean(sizes)
stats['median mesh area'] = np.median(sizes)
stats['mesh area standard deviation'] = np.std(sizes)
structure = np.ones((3,) * skeleton.ndim)
stats['number of disjoint skeletons'] = ndi.label(skeleton, structure)[1]
return stats
def find_peaks_minmax(z, separation=5., threshold=10.):
"""
Method to locate the positive peaks in an image by comparing maximum
and minimum filtered images.
Parameters
----------
z : ndarray
Matrix of image intensities.
separation : float
Expected distance between peaks.
threshold : float
Minimum difference between maximum and minimum filtered images.
Returns
-------
peaks: array with dimensions (npeaks, 2) that contains the x, y coordinates
for each peak found in the image.
"""
data_max = ndi.filters.maximum_filter(z, separation)
maxima = (z == data_max)
data_min = ndi.filters.minimum_filter(z, separation)
diff = ((data_max - data_min) > threshold)
maxima[diff == 0] = 0
labeled, num_objects = ndi.label(maxima)
peaks = np.array(
ndi.center_of_mass(z, labeled, range(1, num_objects + 1)))
return clean_peaks(peaks)
def candidates_to_image(cands,radius):
image_names = []
for subset in xrange(0,10):
subset_names = glob.glob("../data/subset{0}/*.mhd".format(subset))
names = [os.path.split(x)[1].replace('.mhd','') for x in subset_names]
image_names.append(names)
previous_candidate = ""
images = []
image = []
origin = []
spacing = []
number = 0
for candidate in tqdm(cands.values):
if candidate[0] != previous_candidate:
number = 0
previous_candidate = candidate[0]
for image_subset in xrange(0,10):
if candidate[0] in image_names[image_subset]:
image,origin,spacing = load_itk_image("../data/subset{0}/{1}.mhd".format(image_subset,candidate[0]))
break
coords = world_2_voxel([candidate[3],candidate[2],candidate[1]],origin,spacing)
im = image_part_from_candidate(image,coords,radius)
#images.append(im)
if candidate[4]:
label = "true"
else:
label = "false"
scipy.misc.imsave('../data/samples/{0}_{1}_{2}.jpg'.format(candidate[0],number,label), im)
number += 1
return images
def average_boxes(hot_windows, old_boxes,
image_shape):
hot_boxes = initialize_center_box(hot_windows)
new_boxes = add_center_box(hot_boxes, old_boxes)
# print('new_boxes', new_boxes)
filtered_boxes = []
for new_box in new_boxes:
# Draw boxes only if the confidence level is above 1
if new_box[-1] > 2:
filtered_boxes.append(new_box)
new_windows = []
# convert center-width-height to lefttop-rightbottom format
for filtered_box in filtered_boxes:
new_center, new_width, new_height,new_move, new_prob = filtered_box
new_windows.append(((int(new_center[0]-new_width), int(new_center[1]-new_height)),
(int(new_center[0]+new_width), int(new_center[1]+new_height))))
# Create a heatmap
heatmap = create_heatmap(new_windows, image_shape)
# Check if there is any overlap of windows
if np.unique(heatmap)[-1] >= 2:
labels = ndi.label(heatmap)[0]
heatmap_2 = np.zeros_like(heatmap)
heatmap_2[heatmap>=2] = 1
labels_2 = ndi.label(heatmap_2)
array_2 = np.argwhere(labels_2[0])
for car_number in range(1, labels_2[1]+1):
# Find pixels with each car_number label value
nonzero = (labels_2[0] == car_number).nonzero()
# Identify x and y values of those pixels
num = labels[nonzero[0][0], nonzero[1][0]]
labels[labels == num] = 0
heatmap = labels + heatmap_2
new_windows = find_windows_from_heatmap(heatmap)
# return the boxes with high confidence and new set of probability array
return new_windows, new_boxes
def _run_interface(self, runtime):
in_file = nb.load(self.inputs.in_file)
data = in_file.get_data()
mask = np.zeros_like(data, dtype=np.uint8)
mask[data <= 0] = 1
# Pad one pixel to control behavior on borders of binary_opening
mask = np.pad(mask, pad_width=(1,), mode='constant', constant_values=1)
# Remove noise
struc = nd.generate_binary_structure(3, 2)
mask = nd.binary_opening(mask, structure=struc).astype(
np.uint8)
# Remove small objects
label_im, nb_labels = nd.label(mask)
if nb_labels > 2:
sizes = nd.sum(mask, label_im, list(range(nb_labels + 1)))
ordered = list(reversed(sorted(zip(sizes, list(range(nb_labels + 1))))))
for _, label in ordered[2:]:
mask[label_im == label] = 0
# Un-pad
mask = mask[1:-1, 1:-1, 1:-1]
# If mask is small, clean-up
if mask.sum() < 500:
mask = np.zeros_like(mask, dtype=np.uint8)
out_img = in_file.__class__(mask, in_file.affine, in_file.header)
out_img.header.set_data_dtype(np.uint8)
out_file = fname_presuffix(self.inputs.in_file,
suffix='_rotmask', newpath='.')
out_img.to_filename(out_file)
self._results['out_file'] = out_file
return runtime
def makeMask(self):
#Narrow down to pixels which measured something every frame
counts = np.sum(self.Zs > 0, 2)
self.Mask = (counts == self.Zs.shape[2])
#Find biggest connected component out of remaining pixels
ILabel, NLabels = ndimage.label(self.Mask)
idx = np.argmax(ndimage.sum(self.Mask, ILabel, range(NLabels+1)))
self.Mask = (ILabel == idx)
plt.imshow(self.Mask)
plt.show()
#Actually sets up all of the vertex, face, and adjacency structures in the
#PolyMeshObject
def binary_volume_opening(vol, minvol):
lb_vol, num_objs = label(vol)
lbs = np.arange(1, num_objs + 1)
v = labeled_comprehension(lb_vol > 0, lb_vol, lbs, np.sum, np.int, 0)
ix = np.isin(lb_vol, lbs[v >= minvol])
newvol = np.zeros(vol.shape)
newvol[ix] = vol[ix]
return newvol
def autofind_pois(self, neighborhood_size=1, min_threshold=10000, max_threshold=1e6):
"""Automatically search the xy scan image for POIs.
@param neighborhood_size: size in microns. Only the brightest POI per neighborhood will be found.
@param min_threshold: POIs must have c/s above this threshold.
@param max_threshold: POIs must have c/s below this threshold.
"""
# Calculate the neighborhood size in pixels from the image range and resolution
x_range_microns = np.max(self.roi_map_data[:, :, 0]) - np.min(self.roi_map_data[:, :, 0])
y_range_microns = np.max(self.roi_map_data[:, :, 1]) - np.min(self.roi_map_data[:, :, 1])
y_pixels = len(self.roi_map_data)
x_pixels = len(self.roi_map_data[1, :])
pixels_per_micron = np.max([x_pixels, y_pixels]) / np.max([x_range_microns, y_range_microns])
# The neighborhood in pixels is nbhd_size * pixels_per_um, but it must be 1 or greater
neighborhood_pix = int(np.max([math.ceil(pixels_per_micron * neighborhood_size), 1]))
data = self.roi_map_data[:, :, 3]
data_max = filters.maximum_filter(data, neighborhood_pix)
maxima = (data == data_max)
data_min = filters.minimum_filter(data, 3 * neighborhood_pix)
diff = ((data_max - data_min) > min_threshold)
maxima[diff is False] = 0
labeled, num_objects = ndimage.label(maxima)
xy = np.array(ndimage.center_of_mass(data, labeled, range(1, num_objects + 1)))
for count, pix_pos in enumerate(xy):
poi_pos = self.roi_map_data[pix_pos[0], pix_pos[1], :][0:3]
this_poi_key = self.add_poi(position=poi_pos, emit_change=False)
self.rename_poi(poikey=this_poi_key, name='spot' + str(count), emit_change=False)
# Now that all the POIs are created, emit the signal for other things (ie gui) to update
self.signal_poi_updated.emit()
def candidates_to_image(cands,radius):
image_names = []
for subset in xrange(0,10):
subset_names = glob.glob("../data/subset{0}/*.mhd".format(subset))
names = [os.path.split(x)[1].replace('.mhd','') for x in subset_names]
image_names.append(names)
previous_candidate = ""
images = []
image = []
origin = []
spacing = []
number = 0
for candidate in tqdm(cands.values):
if candidate[0] != previous_candidate:
number = 0
previous_candidate = candidate[0]
for image_subset in xrange(0,10):
if candidate[0] in image_names[image_subset]:
image,origin,spacing = load_itk_image("../data/subset{0}/{1}.mhd".format(image_subset,candidate[0]))
break
coords = world_2_voxel([candidate[3],candidate[2],candidate[1]],origin,spacing)
im = image_part_from_candidate(image,coords,radius)
#images.append(im)
if candidate[4]:
label = "true"
else:
label = "false"
scipy.misc.imsave('../data/samples/{0}_{1}_{2}.jpg'.format(candidate[0],number,label), im)
number += 1
return images
def count(n):
a = [(n>>i) & 1 for i in range(8)]
if sum(a)<=1:return False
if a[1] & a[3] & a[5] & a[7]:return False
a = np.array([[a[0],a[1],a[2]],
[a[7], 0 ,a[3]],
[a[6],a[5],a[4]]])
n = label(a, strc)[1]
return n<2
def run(self, ips, snap, img, para = None):
if para == None: para = self.para
ips.lut = self.lut
msk = img>para['thr']
con = 1 if para['con']=='4-Connect' else 2
strc = generate_binary_structure(2, con)
msk = label(msk, strc, output = np.int16)
IPy.show_img([msk[0]], ips.title+'-label')
def run(self, ips, snap, img, para = None):
intenimg = WindowsManager.get(para['inten']).ips.img
strc = ndimage.generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
buf, n = ndimage.label(snap, strc, output=np.uint16)
index = range(1, n+1)
idx = (np.ones(n+1)*para['front']).astype(np.uint8)
msk = np.ones(n, dtype=np.bool)
if para['mean']>0: msk *= ndimage.mean(intenimg, buf, index)>=para['mean']
if para['mean']<0: msk *= ndimage.mean(intenimg, buf, index)<-para['mean']
if para['max']>0: msk *= ndimage.maximum(intenimg, buf, index)>=para['max']
if para['max']<0: msk *= ndimage.maximum(intenimg, buf, index)<-para['max']
if para['min']>0: msk *= ndimage.minimum(intenimg, buf, index)>=para['min']
if para['min']<0: msk *= ndimage.minimum(intenimg, buf, index)<-para['min']
if para['sum']>0: msk *= ndimage.sum(intenimg, buf, index)>=para['sum']
if para['sum']<0: msk *= ndimage.sum(intenimg, buf, index)<-para['sum']
if para['std']>0: msk *= ndimage.standard_deviation(intenimg, buf, index)>=para['std']
if para['std']<0: msk *= ndimage.standard_deviation(intenimg, buf, index)<-para['std']
xy = ndimage.center_of_mass(intenimg, buf, index)
xy = np.array(xy).round(2).T
idx[1:][~msk] = para['back']
idx[0] = 0
img[:] = idx[buf]
WindowsManager.get(para['inten']).ips.mark = RGMark((xy.T, msk))
WindowsManager.get(para['inten']).ips.update = True
def run(self, ips, snap, img, para = None):
ips.lut = self.buflut
lab, n = ndimg.label(snap>para['thr2'], np.ones((3,3)), output=np.uint16)
sta = ndimg.sum(snap>para['thr1'], lab, index=range(n+1)) > 0
img[:] = (sta*255)[lab]
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, imgs, para = None):
k = ips.unit[0]
titles = ['ID']
if para['center']:titles.extend(['Center-X','Center-Y','Center-Z'])
if para['vol']:titles.append('Volume')
if para['extent']:titles.extend(['Min-Z','Min-Y','Min-X','Max-Z','Max-Y','Max-X'])
if para['ed']:titles.extend(['Diameter'])
if para['fa']:titles.extend(['FilledArea'])
buf = imgs.astype(np.uint16)
strc = generate_binary_structure(3, 1 if para['con']=='4-connect' else 2)
label(imgs, strc, output=buf)
ls = regionprops(buf)
dt = [range(len(ls))]
centroids = [i.centroid for i in ls]
if para['center']:
dt.append([round(i.centroid[1]*k,1) for i in ls])
dt.append([round(i.centroid[0]*k,1) for i in ls])
dt.append([round(i.centroid[2]*k,1) for i in ls])
if para['vol']:
dt.append([i.area*k**3 for i in ls])
if para['extent']:
for j in (0,1,2,3,4,5):
dt.append([i.bbox[j]*k for i in ls])
if para['ed']:
dt.append([round(i.equivalent_diameter*k, 1) for i in ls])
if para['fa']:
dt.append([i.filled_area*k**3 for i in ls])
IPy.table(ips.title+'-region', list(zip(*dt)), titles)
# center, area, l, extent, cov
