def dilated(self, structure=None, connectivity=2, iterations=100):
"""
This function ...
:param connectivity:
:param iterations:
:return:
"""
# Define the structure for the expansion
if structure is None: structure = ndimage.generate_binary_structure(2, connectivity=connectivity)
# Make the new mask, made from 100 iterations with the structure array
data = ndimage.binary_dilation(self, structure, iterations)
# Return the dilated mask
#data, name=None, description=None
return Mask(data, name=self.name, description=self.description)
# -----------------------------------------------------------------
python类generate_binary_structure()的实例源码
def dilated(self, structure=None, connectivity=2, iterations=100):
"""
This function ...
:param connectivity:
:param iterations:
:return:
"""
# Define the structure for the expansion
if structure is None: structure = ndimage.generate_binary_structure(2, connectivity=connectivity)
# Make the new mask, made from 100 iterations with the structure array
data = ndimage.binary_dilation(self, structure, iterations)
# Return the dilated mask
#data, name=None, description=None
return Mask(data, name=self.name, description=self.description)
# -----------------------------------------------------------------
def _run_interface(self, runtime):
in_file = nb.load(self.inputs.in_file)
wm_mask = nb.load(self.inputs.wm_mask).get_data()
wm_mask[wm_mask < 0.9] = 0
wm_mask[wm_mask > 0] = 1
wm_mask = wm_mask.astype(np.uint8)
if self.inputs.erodemsk:
# Create a structural element to be used in an opening operation.
struc = nd.generate_binary_structure(3, 2)
# Perform an opening operation on the background data.
wm_mask = nd.binary_erosion(wm_mask, structure=struc).astype(np.uint8)
data = in_file.get_data()
data *= 1000.0 / np.median(data[wm_mask > 0])
out_file = fname_presuffix(self.inputs.in_file,
suffix='_harmonized', newpath='.')
in_file.__class__(data, in_file.affine, in_file.header).to_filename(
out_file)
self._results['out_file'] = out_file
return runtime
def _prepare_mask(mask, label, erode=True):
fgmask = mask.copy()
if np.issubdtype(fgmask.dtype, np.integer):
if isinstance(label, string_types):
label = FSL_FAST_LABELS[label]
fgmask[fgmask != label] = 0
fgmask[fgmask == label] = 1
else:
fgmask[fgmask > .95] = 1.
fgmask[fgmask < 1.] = 0
if erode:
# Create a structural element to be used in an opening operation.
struc = nd.generate_binary_structure(3, 2)
# Perform an opening operation on the background data.
fgmask = nd.binary_opening(fgmask, structure=struc).astype(np.uint8)
return fgmask
def run(self, ips, imgs, para = None):
k, unit = ips.unit
strc = generate_binary_structure(3, 1 if para['con']=='4-connect' else 2)
lab, n = label(imgs==0 if para['inv'] else imgs, strc, output=np.uint16)
idx = (np.ones(n+1)*(0 if para['inv'] else para['front'])).astype(np.uint8)
ls = regionprops(lab)
for i in ls:
if para['vol'] == 0: break
if para['vol']>0:
if i.area*k**3 < para['vol']: idx[i.label] = para['back']
if para['vol']<0:
if i.area*k**3 >= -para['vol']: idx[i.label] = para['back']
for i in ls:
if para['dia'] == 0: break
d = norm(np.array(i.bbox[:3]) - np.array(i.bbox[3:]))
if para['dia']>0:
if d*k < para['dia']: idx[i.label] = para['back']
if para['dia']<0:
if d*k >= -para['dia']: idx[i.label] = para['back']
idx[0] = para['front'] if para['inv'] else 0
imgs[:] = idx[lab]
def generate_markers_3d(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal_labels = np.zeros(image.shape).astype(np.int16)
for i in range(marker_internal.shape[0]):
marker_internal[i] = segmentation.clear_border(marker_internal[i])
marker_internal_labels[i] = measure.label(marker_internal[i])
#areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas = [r.area for i in range(marker_internal.shape[0]) for r in measure.regionprops(marker_internal_labels[i])]
for i in range(marker_internal.shape[0]):
areas = [r.area for r in measure.regionprops(marker_internal_labels[i])]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels[i]):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[i, coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
# 3x3 structuring element with connectivity 1, used by default
struct1 = ndimage.generate_binary_structure(2, 1)
struct1 = struct1[np.newaxis,:,:] # expand by z axis .
external_a = ndimage.binary_dilation(marker_internal, structure=struct1, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, structure=struct1, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
#marker_watershed = np.zeros((512, 512), dtype=np.int) # origi
marker_watershed = np.zeros((marker_external.shape), dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def generate_markers_3d(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal_labels = np.zeros(image.shape).astype(np.int16)
for i in range(marker_internal.shape[0]):
marker_internal[i] = segmentation.clear_border(marker_internal[i])
marker_internal_labels[i] = measure.label(marker_internal[i])
#areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas = [r.area for i in range(marker_internal.shape[0]) for r in measure.regionprops(marker_internal_labels[i])]
for i in range(marker_internal.shape[0]):
areas = [r.area for r in measure.regionprops(marker_internal_labels[i])]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels[i]):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[i, coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
# 3x3 structuring element with connectivity 1, used by default
struct1 = ndimage.generate_binary_structure(2, 1)
struct1 = struct1[np.newaxis,:,:] # expand by z axis .
external_a = ndimage.binary_dilation(marker_internal, structure=struct1, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, structure=struct1, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
#marker_watershed = np.zeros((512, 512), dtype=np.int) # origi
marker_watershed = np.zeros((marker_external.shape), dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def __voxel_first_layer(self, keep_background=True):
"""
Extract the first layer of voxels at the surface of the biological object.
"""
print "Extracting the first layer of voxels..."
mask_img_1 = (self.image == self.background())
struct = nd.generate_binary_structure(3, 1)
dil_1 = nd.binary_dilation(mask_img_1, structure=struct)
layer = dil_1 - mask_img_1
if keep_background:
return self.image * layer + mask_img_1
else:
return self.image * layer
def eroded(self, structure=None, connectivity=2, iterations=100):
"""
This function ...
:param structure:
:param connectivity:
:param iterations:
:return:
"""
# Define the structure for the expansion
if structure is None: structure = ndimage.generate_binary_structure(2, connectivity=connectivity)
try:
# Make the new mask, made from 100 iterations with the structure array
data = ndimage.binary_erosion(self, structure, iterations)
except:
#print(self)
#print(structure)
data = np.zeros((self.ysize,self.xsize), dtype=bool)
# Reassign this object
#data, name=None, description=None
return Mask(data, name=self.name, description=self.description)
# -----------------------------------------------------------------
def eroded(self, structure=None, connectivity=2, iterations=100):
"""
This function ...
:param structure:
:param connectivity:
:param iterations:
:return:
"""
# Define the structure for the expansion
if structure is None: structure = ndimage.generate_binary_structure(2, connectivity=connectivity)
try:
# Make the new mask, made from 100 iterations with the structure array
data = ndimage.binary_erosion(self, structure, iterations)
except:
#print(self)
#print(structure)
data = np.zeros((self.ysize,self.xsize), dtype=bool)
# Reassign this object
#data, name=None, description=None
return Mask(data, name=self.name, description=self.description)
# -----------------------------------------------------------------
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 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 get_morphological_patch(dimension, shape):
"""
:param dimension: dimension of the image (NOT the shape).
:param shape: circle or square.
:return: morphological patch as ndimage
"""
if shape == 'circle':
morpho_patch = ndimage.generate_binary_structure(dimension, 1)
elif shape == 'square':
morpho_patch = ndimage.generate_binary_structure(dimension, 3)
else:
return IOError
return morpho_patch
def weights_map(ys):
"""Compute corresponding weight map when use cross entropy loss.
Argument:
ys: [depth, height, width]
Return:
weights_map: [depth, height, width]
"""
weights = ys.astype(np.float64)
# Balance class frequencies.
cls_ratio = np.sum(1 - ys) / np.sum(ys)
weights *= cls_ratio
# Generate boundaries using morphological operation.
se = generate_binary_structure(3, 1)
bigger = binary_dilation(ys, structure=se).astype(np.float64)
small = binary_erosion(ys, structure=se).astype(np.float64)
edge = bigger - small
# Balance edge frequencies.
edge_ratio = np.sum(bigger) / np.sum(edge)
edge *= np.exp(edge_ratio) * 10
# `weights` should > 0
# `targets * -log(sigmoid(logits)) * pos_weight + (1 - targets) * -log(1 - sigmoid(logits))`
return weights + edge + 1
def floodfill(img, x, y, thr, con):
color = img[int(y), int(x)]
buf = np.subtract(img, color, dtype=np.int16)
msk = np.abs(buf)<=thr
if buf.ndim==3:
msk = np.min(msk, axis=2)
buf = buf[:,:,0]
strc = generate_binary_structure(2, con+1)
label(msk, strc, output = buf)
msk = buf == buf[int(y), int(x)]
#msk[[0,-1],:], msk[:,[0,-1]] = 0, 0
#imsave('test.png', msk.astype(np.uint8))
return msk
def neighbors(shape, conn=1):
dim = len(shape)
block = generate_binary_structure(dim, conn)
block[tuple([1]*dim)] = 0
idx = np.where(block>0)
idx = np.array(idx, dtype=np.uint8).T
idx = np.array(idx-[1]*dim)
acc = np.cumprod((1,)+shape[::-1][:-1])
return np.dot(idx, acc[::-1])
def neighbors(shape):
dim = len(shape)
block = generate_binary_structure(dim, 1)
block[tuple([1]*dim)] = 0
idx = np.where(block>0)
idx = np.array(idx, dtype=np.uint8).T
idx = np.array(idx-[1]*dim)
acc = np.cumprod((1,)+shape[::-1][:-1])
return np.dot(idx, acc[::-1])
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, 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
def dilate(y):
mask = ndimage.generate_binary_structure(2, 2)
y_f = ndimage.binary_dilation(y, structure=mask).astype(y.dtype)
return y_f
##########################
### Array routine Operation
##########################
def wall_voxels_between_two_cells(self, label_1, label_2, bbox=None, verbose=False):
"""
Return the voxels coordinates defining the contact wall between two labels.
Args:
image: (ndarray of ints) - Array containing objects defined by labels
label_1: (int) - object id #1
label_2: (int) - object id #2
bbox: (dict, optional) - If given, contain a dict of slices
Returns:
- xyz 3xN array.
"""
if bbox is not None:
if isinstance(bbox, dict):
label_1, label_2 = sort_boundingbox(bbox, label_1, label_2)
boundingbox = bbox[label_1]
elif isinstance(bbox, tuple) and len(bbox)==3:
boundingbox = bbox
else:
try:
boundingbox = find_smallest_boundingbox(self.image, label_1, label_2)
except:
print "Could neither use the provided value of `bbox`, nor gess it!"
boundingbox = tuple([(0,s-1,None) for s in self.image.shape])
dilated_bbox = dilation( boundingbox )
dilated_bbox_img = self.image[dilated_bbox]
else:
try:
boundingbox = find_smallest_boundingbox(self.image, label_1, label_2)
except:
dilated_bbox_img = self.image
mask_img_1 = (dilated_bbox_img == label_1)
mask_img_2 = (dilated_bbox_img == label_2)
struct = nd.generate_binary_structure(3, 2)
dil_1 = nd.binary_dilation(mask_img_1, structure=struct)
dil_2 = nd.binary_dilation(mask_img_2, structure=struct)
x,y,z = np.where( ( (dil_1 & mask_img_2) | (dil_2 & mask_img_1) ) == 1 )
if bbox is not None:
return np.array( (x+dilated_bbox[0].start, y+dilated_bbox[1].start, z+dilated_bbox[2].start) )
else:
return np.array( (x, y, z) )
def cells_voxel_layer(self, labels, region_boundingbox = False, single_frame = False):
"""
This function extract the first layer of voxel surrounding a cell defined by `label`
Args:
label: (int|list) - cell-label for which we want to extract the first layer of voxel;
region_boundingbox: (bool) - if True, consider a boundingbox surrounding all labels, instead of each label alone.
single_frame: (bool) - if True, return only one array with all voxels position defining cell walls.
:Output:
returns a binary image: 1 where the cell-label of interest is, 0 elsewhere
"""
if isinstance(labels,int):
labels = [labels]
if single_frame:
region_boundingbox=True
if not isinstance(region_boundingbox,bool):
if sum([isinstance(s,slice) for s in region_boundingbox])==3:
bbox = region_boundingbox
else:
print "TypeError: Wong type for 'region_boundingbox', should either be bool or la tuple of slices"
return None
elif isinstance(region_boundingbox,bool) and region_boundingbox:
bbox = self.region_boundingbox(labels)
else:
bboxes = self.boundingbox(labels, real=False)
# Generate the smaller eroding structure possible:
struct = nd.generate_binary_structure(3, 2)
if single_frame:
vox_layer = np.zeros_like(self.image[bbox], dtype=int)
else:
vox_layer = {}
for clabel in labels:
if region_boundingbox:
bbox_im = self.image[bbox]
else:
bbox_im = self.image[bboxes[clabel]]
# Creating a mask (1 where the cell-label of interest is, 0 elsewhere):
mask_bbox_im = (bbox_im == clabel)
# Erode the cell using the structure:
eroded_mask_bbox_im = nd.binary_erosion(mask_bbox_im, structure=struct)
if single_frame:
vox_layer += np.array(mask_bbox_im - eroded_mask_bbox_im, dtype=int)
else:
vox_layer[clabel] = np.array(mask_bbox_im - eroded_mask_bbox_im, dtype=int)
if len(labels)==1:
return vox_layer[clabel]
else:
return vox_layer
def median_otsu(unweighted_volume, median_radius=4, numpass=4, dilate=1):
""" Simple brain extraction tool for dMRI data.
This function is inspired from the ``median_otsu`` function from ``dipy``
and is copied here to remove a dependency.
It uses a median filter smoothing of the ``unweighted_volume``
automatic histogram Otsu thresholding technique, hence the name
*median_otsu*.
This function is inspired from Mrtrix's bet which has default values
``median_radius=3``, ``numpass=2``. However, from tests on multiple 1.5T
and 3T data. From GE, Philips, Siemens, the most robust choice is
``median_radius=4``, ``numpass=4``.
Args:
unweighted_volume (ndarray): ndarray of the unweighted volumes brain volumes
median_radius (int): Radius (in voxels) of the applied median filter (default 4)
numpass (int): Number of pass of the median filter (default 4)
dilate (None or int): optional number of iterations for binary dilation
Returns:
ndarray: a 3D ndarray with the binary brain mask
"""
b0vol = unweighted_volume
logger = logging.getLogger(__name__)
logger.info('We will use a single precision float type for the calculations.'.format())
for env in mot.configuration.get_load_balancer().get_used_cl_environments(mot.configuration.get_cl_environments()):
logger.info('Using device \'{}\'.'.format(str(env)))
m = MedianFilter(median_radius)
b0vol = m.filter(b0vol, nmr_of_times=numpass)
thresh = _otsu(b0vol)
mask = b0vol > thresh
if dilate is not None:
cross = generate_binary_structure(3, 1)
mask = binary_dilation(mask, cross, iterations=dilate)
return mask
def run(self, ips, imgs, para = None):
inten = WindowsManager.get(para['inten']).ips
if not para['slice']:
imgs = [inten.img]
msks = [ips.img]
else:
msks = ips.imgs
if len(msks)==1:
msks *= len(imgs)
buf = imgs[0].astype(np.uint16)
strc = ndimage.generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
idct = ['Max','Min','Mean','Variance','Standard','Sum']
key = {'Max':'max','Min':'min','Mean':'mean',
'Variance':'var','Standard':'std','Sum':'sum'}
idct = [i for i in idct if para[key[i]]]
titles = ['Slice', 'ID'][0 if para['slice'] else 1:]
if para['center']: titles.extend(['Center-X','Center-Y'])
if para['extent']: titles.extend(['Min-Y','Min-X','Max-Y','Max-X'])
titles.extend(idct)
k = ips.unit[0]
data, mark = [], []
for i in range(len(imgs)):
n = ndimage.label(msks[i], strc, output=buf)
index = range(1, n+1)
dt = []
if para['slice']:dt.append([i]*n)
dt.append(range(n))
xy = ndimage.center_of_mass(imgs[i], buf, index)
xy = np.array(xy).round(2).T
if para['center']:dt.extend([xy[1]*k, xy[0]*k])
boxs = [None] * n
if para['extent']:
boxs = ndimage.find_objects(buf)
boxs = [(i[0].start, i[1].start, i[0].stop, i[1].stop) for i in boxs]
for j in (0,1,2,3):
dt.append([i[j]*k for i in boxs])
if para['max']:dt.append(ndimage.maximum(imgs[i], buf, index).round(2))
if para['min']:dt.append(ndimage.minimum(imgs[i], buf, index).round(2))
if para['mean']:dt.append(ndimage.mean(imgs[i], buf, index).round(2))
if para['var']:dt.append(ndimage.variance(imgs[i], buf, index).round(2))
if para['std']:dt.append(ndimage.standard_deviation(imgs[i], buf, index).round(2))
if para['sum']:dt.append(ndimage.sum(imgs[i], buf, index).round(2))
mark.append([(center, cov) for center,cov in zip(xy.T, boxs)])
data.extend(list(zip(*dt)))
IPy.table(inten.title+'-region statistic', data, titles)
inten.mark = Mark(mark)
inten.update = True
def run(self, ips, imgs, para = None):
if not para['slice']:imgs = [ips.img]
k = ips.unit[0]
titles = ['Slice', 'ID'][0 if para['slice'] else 1:]
if para['center']:titles.extend(['Center-X','Center-Y'])
if para['area']:titles.append('Area')
if para['l']:titles.append('Perimeter')
if para['extent']:titles.extend(['Min-Y','Min-X','Max-Y','Max-X'])
if para['ed']:titles.extend(['Diameter'])
if para['ca']:titles.extend(['ConvexArea'])
if para['holes']:titles.extend(['Holes'])
if para['fa']:titles.extend(['FilledArea'])
if para['solid']:titles.extend(['Solidity'])
if para['cov']:titles.extend(['Major','Minor','Ori'])
buf = imgs[0].astype(np.uint16)
data, mark = [], []
strc = generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
for i in range(len(imgs)):
label(imgs[i], strc, output=buf)
ls = regionprops(buf)
dt = [[i]*len(ls), list(range(len(ls)))]
if not para['slice']:dt = dt[1:]
if not para['cov']: cvs = [None] * len(ls)
else: cvs = [(i.major_axis_length, i.minor_axis_length, i.orientation) for i in ls]
centroids = [i.centroid for i in ls]
mark.append([(center, cov) for center,cov in zip(centroids, cvs)])
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])
if para['area']:
dt.append([i.area*k**2 for i in ls])
if para['l']:
dt.append([round(i.perimeter*k,1) for i in ls])
if para['extent']:
for j in (0,1,2,3):
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['ca']:
dt.append([i.convex_area*k**2 for i in ls])
if para['holes']:
dt.append([1-i.euler_number for i in ls])
if para['fa']:
dt.append([i.filled_area*k**2 for i in ls])
if para['solid']:
dt.append([round(i.solidity, 2) for i in ls])
if para['cov']:
dt.append([round(i.major_axis_length*k, 1) for i in ls])
dt.append([round(i.minor_axis_length*k, 1) for i in ls])
dt.append([round(i.orientation*k, 1) for i in ls])
data.extend(list(zip(*dt)))
ips.mark = Mark(mark)
IPy.table(ips.title+'-region', data, titles)
# center, area, l, extent, cov
def run(self, ips, snap, img, para = None):
k, unit = ips.unit
strc = generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
lab, n = label(snap==0 if para['inv'] else snap, strc, output=np.uint16)
idx = (np.ones(n+1)*(0 if para['inv'] else para['front'])).astype(np.uint8)
ls = regionprops(lab)
for i in ls:
if para['area'] == 0: break
if para['area']>0:
if i.area*k**2 < para['area']: idx[i.label] = para['back']
if para['area']<0:
if i.area*k**2 >= -para['area']: idx[i.label] = para['back']
for i in ls:
if para['l'] == 0: break
if para['l']>0:
if i.perimeter*k < para['l']: idx[i.label] = para['back']
if para['l']<0:
if i.perimeter*k >= -para['l']: idx[i.label] = para['back']
for i in ls:
if para['holes'] == 0: break
if para['holes']>0:
if 1-i.euler_number < para['holes']: idx[i.label] = para['back']
if para['holes']<0:
if 1-i.euler_number >= -para['holes']: idx[i.label] = para['back']
for i in ls:
if para['solid'] == 0: break
if para['solid']>0:
if i.solidity < para['solid']: idx[i.label] = para['back']
if para['solid']<0:
if i.solidity >= -para['solid']: idx[i.label] = para['back']
for i in ls:
if para['e'] == 0: break
if para['e']>0:
if i.minor_axis_length>0 and i.major_axis_length/i.minor_axis_length < para['e']:
idx[i.label] = para['back']
if para['e']<0:
if i.minor_axis_length>0 and i.major_axis_length/i.minor_axis_length >= -para['e']:
idx[i.label] = para['back']
idx[0] = para['front'] if para['inv'] else 0
img[:] = idx[lab]
def labelHealpix(pixels, values, nside, threshold=0, xsize=1000):
"""
Label contiguous regions of a (sparse) HEALPix map. Works by mapping
HEALPix array to a Mollweide projection and applying scipy.ndimage.label
Assumes non-nested HEALPix map.
Parameters:
pixels : Pixel values associated to (sparse) HEALPix array
values : (Sparse) HEALPix array of data values
nside : HEALPix dimensionality
threshold : Threshold value for object detection
xsize : Size of Mollweide projection
Returns:
labels, nlabels
"""
proj = healpy.projector.MollweideProj(xsize=xsize)
vec = healpy.pix2vec(nside,pixels)
xy = proj.vec2xy(vec)
ij = proj.xy2ij(xy)
xx,yy = proj.ij2xy()
# Convert to Mollweide
searchims = []
if values.ndim < 2: iterate = [values]
else: iterate = values.T
for i,value in enumerate(iterate):
logger.debug("Labeling slice %i...")
searchim = numpy.zeros(xx.shape,dtype=bool)
select = (value > threshold)
yidx = ij[0][select]; xidx = ij[1][select]
searchim[yidx,xidx] |= True
searchims.append( searchim )
searchims = numpy.array(searchims)
# Full binary structure
s = ndimage.generate_binary_structure(searchims.ndim,searchims.ndim)
### # Dilate in the z-direction
logger.info(" Dilating image...")
searchims = ndimage.binary_dilation(searchims,s,1)
# Do the labeling
logger.info(" Labeling image...")
labels,nlabels = ndimage.label(searchims,structure=s)
# Convert back to healpix
pix_labels = labels[:,ij[0],ij[1]].T
pix_labels = pix_labels.reshape(values.shape)
pix_labels *= (values > threshold) # re-trim
return pix_labels, nlabels
