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类binary_dilation()的实例源码
def disk_dilation(self, radius=5, iterations=1):
"""
This function ...
:param radius:
:return:
"""
structure = morphology.disk(radius, dtype=bool)
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 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 disk_dilation(self, radius=5, iterations=1):
"""
This function ...
:param radius:
:return:
"""
structure = morphology.disk(radius, dtype=bool)
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 generate_markers(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal = segmentation.clear_border(marker_internal)
marker_internal_labels = measure.label(marker_internal)
areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
external_a = ndimage.binary_dilation(marker_internal, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
marker_watershed = np.zeros(image.shape, dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def generate_markers(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal = segmentation.clear_border(marker_internal)
marker_internal_labels = measure.label(marker_internal)
areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
external_a = ndimage.binary_dilation(marker_internal, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
marker_watershed = np.zeros(image.shape, dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def generate_markers(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal = segmentation.clear_border(marker_internal)
marker_internal_labels = measure.label(marker_internal)
areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
external_a = ndimage.binary_dilation(marker_internal, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
marker_watershed = np.zeros(image.shape, dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def generate_markers(image):
#Creation of the internal Marker
marker_internal = image < -400
marker_internal = segmentation.clear_border(marker_internal)
marker_internal_labels = measure.label(marker_internal)
areas = [r.area for r in measure.regionprops(marker_internal_labels)]
areas.sort()
if len(areas) > 2:
for region in measure.regionprops(marker_internal_labels):
if region.area < areas[-2]:
for coordinates in region.coords:
marker_internal_labels[coordinates[0], coordinates[1]] = 0
marker_internal = marker_internal_labels > 0
#Creation of the external Marker
external_a = ndimage.binary_dilation(marker_internal, iterations=10)
external_b = ndimage.binary_dilation(marker_internal, iterations=55)
marker_external = external_b ^ external_a
#Creation of the Watershed Marker matrix
marker_watershed = np.zeros(image.shape, dtype=np.int)
marker_watershed += marker_internal * 255
marker_watershed += marker_external * 128
return marker_internal, marker_external, marker_watershed
def separate_data_mask(self):
def ret_data(image):
if isinstance(image, np.ma.MaskedArray):
image_data = image.data
else:
image_data = image
return image_data
badpixmask = None
if _has_mask(self.refimage):
badpixmask = ndimage.binary_dilation(
self.refimage.mask.astype('uint8'),
structure=np.ones(self.k_shape)).astype('bool')
if _has_mask(self.image):
badpixmask += self.image.mask
elif _has_mask(self.image):
badpixmask = self.image.mask
return ret_data(self.image), ret_data(self.refimage), badpixmask
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 wall(mask_img, label_id):
"""
Function detecting wall position for a given cell-id (`label_id`) within a segmented image (`mask_img`).
"""
img = (mask_img == label_id)
dil = nd.binary_dilation(img)
contact = dil - img
return mask_img[contact]
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 color_blend(rgb1, depth1, rgb2, depth2):
mask = np.all(rgb1 == 0, axis=2)
mask = ndimage.binary_dilation(mask).astype(mask.dtype)
depth1[mask] = 0
rgb1[mask, :] = 0
mask = mask.astype(np.uint8)
new_depth = depth2 * mask + depth1
new_color = rgb2 * mask[:, :, np.newaxis] + rgb1
return new_color.astype(np.uint8), new_depth
def depth_blend(rgb1, depth1, rgb2, depth2):
new_depth2 = depth2.copy()
new_depth1 = depth1.copy()
rgb1_mask = np.all(rgb1 == 0, axis=2)
rgb2_mask = np.all(rgb2 == 0, axis=2)
rgb1_mask = ndimage.binary_dilation(rgb1_mask)
new_depth2[rgb2_mask] = -100000
new_depth1[rgb1_mask] = -100000
mask = (new_depth1 < new_depth2)
pos_mask = mask.astype(np.uint8)
neg_mask = (mask == False).astype(np.uint8)
masked_rgb_occluder = rgb1 * pos_mask[:, :, np.newaxis]
masked_rgb_object = rgb2 * neg_mask[:, :, np.newaxis]
masked_depth_occluder = depth1 * pos_mask
masked_depth_object = depth2 * neg_mask
blend_rgb = masked_rgb_occluder + masked_rgb_object
blend_depth = masked_depth_occluder + masked_depth_object
return blend_rgb, blend_depth
def resize_save(path, size):
path = os.path.join(os.getcwd(), path)
images_files = os.listdir(path)
el = (np.array([[0,1],[1,1]])).astype(np.float32)
for i, image in enumerate(images_files):
image_file = os.path.join(path,image)
try:
im = (ndimage.imread(image_file)).astype(np.float32)
np.putmask(im, im < 100, 0)
im = ndimage.binary_dilation(im, structure=el)
im = (misc.imresize(im, (size, size))).astype(np.uint8)
new_im = Image.fromarray(im)
new_im.save(image_file)
except Exception as e:
print('Could not read:', image_file, ':', e, '- it\'s ok, skipping.')
def __init__(self, tps, obs, pivots, dtps=None, **kwargs):
super(freq_est_clipped, self).__init__()
tmp_obs = np.array(sorted(zip(tps, obs), key=lambda x:x[0]))
self.tps = tmp_obs[:,0]
self.obs = np.array(tmp_obs[:,1], dtype=bool)
self.pivots = pivots
pivot_dt = self.pivots[1]-self.pivots[0]
if dtps==None:
self.dtps = 6.0*pivot_dt
else:
self.dtps = np.max(dtps, pivot_dt)
cum_obs = np.diff(self.obs).cumsum()
first_obs = max(pivots[0], self.tps[cum_obs.searchsorted(cum_obs[0]+1)])
last_obs = min(pivots[-1], max(first_obs,self.tps[min(len(self.tps)-1, 20+cum_obs.searchsorted(cum_obs[-1]))]))
tps_lower_cutoff = first_obs - self.dtps
tps_upper_cutoff = last_obs + self.dtps
self.good_tps = (self.tps>=tps_lower_cutoff)&(self.tps<tps_upper_cutoff)&(self.tps<self.pivots[-1])&(self.tps>=self.pivots[0])
self.valid = True
if self.good_tps.sum()<3:
print("too few valid time points:", self.good_tps.sum())
self.valid=False
return None
if self.good_tps.sum()<7:
from scipy.ndimage import binary_dilation
self.good_tps = binary_dilation(self.good_tps, iterations=5)
reduced_obs = self.obs[self.good_tps]
reduced_tps = self.tps[self.good_tps]
self.pivot_lower_cutoff = min(reduced_tps[0], tps_lower_cutoff)-pivot_dt
self.pivot_upper_cutoff = max(reduced_tps[-1], tps_upper_cutoff)+pivot_dt
self.good_pivots = (self.pivots>=self.pivot_lower_cutoff)\
&(self.pivots<self.pivot_upper_cutoff)
if self.good_pivots.sum()<2:
from scipy.ndimage import binary_dilation
self.good_pivots = binary_dilation(self.good_pivots, iterations=2)
self.fe = frequency_estimator(reduced_tps, reduced_obs,
self.pivots[self.good_pivots], **kwargs)
def get_morphological_mask(point, omega, radius=5, shape='circle', morpho_patch=None):
if morpho_patch is None:
d = omega.ndim
morpho_patch = get_morphological_patch(d, shape=shape)
mask = np.zeros(omega, dtype=np.bool)
mask.itemset(tuple(point), 1)
for _ in range(radius):
mask = ndimage.binary_dilation(mask, structure=morpho_patch).astype(mask.dtype)
return mask
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 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 run(self, ips, snap, img, para = None):
strc = np.ones((para['h'], para['w']), dtype=np.uint8)
ndimg.binary_dilation(snap, strc, output=img)
img *= 255
def run(self, ips, snap, img, para = None):
ndimg.binary_dilation(snap, output=img)
img *= 255
img -= snap
def run(self, ips, imgs, para = None):
strc = np.ones((para['r'], para['r'],para['r']), dtype=np.uint8)
imgs[:] = ndimg.binary_dilation(imgs, strc)
imgs *= 255
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 calc_fwhm(distribution, is_log=True):
"""
Assess the width of the probability distribution. This returns
full-width-half-max
"""
if isinstance(distribution, interp1d):
if is_log:
ymin = distribution.y.min()
prob = np.exp(-(distribution.y-ymin))
else:
prob = distribution.y
xvals = distribution.x
elif isinstance(distribution, Distribution):
# Distribution always stores log-prob
xvals = distribution._func.x
prob = distribution.prob_relative(xvals)
else:
raise TypeError("Error in computing the FWHM for the distribution. "
" The input should be either Distribution or interpolation object");
x_idxs = binary_dilation(prob >= 0.4*(prob.max() - prob.min())+prob.min(), iterations=1)
xs = xvals[x_idxs]
if xs.shape[0] < 2:
print ("Not enough points to compute FWHM: returning zero")
return min(TINY_NUMBER, distribution.xmax - distribution.xmin)
else:
return xs.max() - xs.min()
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 cell_wall_area(self, label_id, neighbors, real = True):
"""
Return the area of contact between a label and its neighbors.
A list or a unique id can be given as neighbors.
:Examples:
>>> import numpy as np
>>> a = np.array([[1, 2, 7, 7, 1, 1],
[1, 6, 5, 7, 3, 3],
[2, 2, 1, 7, 3, 3],
[1, 1, 1, 4, 1, 1]])
>>> from vplants.tissue_analysis.spatial_image_analysis import SpatialImageAnalysis
>>> analysis = SpatialImageAnalysis(a)
>>> analysis.cell_wall_area(7,2)
1.0
>>> analysis.cell_wall_area(7,[2,5])
{(2, 7): 1.0, (5, 7): 2.0}
"""
resolution = self.get_voxel_face_surface()
try:
dilated_bbox = dilation(self.boundingbox(label_id))
dilated_bbox_img = self.image[dilated_bbox]
except:
#~ dilated_bbox = tuple( [slice(0,self.image.shape[i]-1) for i in xrange(len(self.image.shape))] ) #if no slice can be found we use the whole image
dilated_bbox_img = self.image
mask_img = (dilated_bbox_img == label_id)
xyz_kernels = self.neighbor_kernels()
unique_neighbor = not isinstance(neighbors,list)
if unique_neighbor:
neighbors = [neighbors]
wall = {}
for a in xrange(len(xyz_kernels)):
dil = nd.binary_dilation(mask_img, structure=xyz_kernels[a])
frontier = dilated_bbox_img[dil-mask_img]
for n in neighbors:
nb_pix = len(frontier[frontier==n])
if real: area = float(nb_pix*resolution[a//2])
else : area = nb_pix
i,j = min(label_id,n), max(label_id,n)
wall[(i,j)] = wall.get((i,j),0.0) + area
if unique_neighbor: return wall.itervalues().next()
else : return wall
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 gradient_threshold(in_file, in_segm, thresh=1.0, out_file=None):
""" Compute a threshold from the histogram of the magnitude gradient image """
import os.path as op
import numpy as np
import nibabel as nb
from scipy import ndimage as sim
struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 2)
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_gradmask{}'.format(fname, ext))
imnii = nb.load(in_file)
hdr = imnii.get_header().copy()
hdr.set_data_dtype(np.uint8) # pylint: disable=no-member
data = imnii.get_data().astype(np.float32)
mask = np.zeros_like(data, dtype=np.uint8) # pylint: disable=no-member
mask[data > 15.] = 1
segdata = nb.load(in_segm).get_data().astype(np.uint8)
segdata[segdata > 0] = 1
segdata = sim.binary_dilation(segdata, struc, iterations=2, border_value=1).astype(np.uint8) # pylint: disable=no-member
mask[segdata > 0] = 1
mask = sim.binary_closing(mask, struc, iterations=2).astype(np.uint8) # pylint: disable=no-member
# Remove small objects
label_im, nb_labels = sim.label(mask)
artmsk = np.zeros_like(mask)
if nb_labels > 2:
sizes = sim.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
artmsk[label_im == label] = 1
mask = sim.binary_fill_holes(mask, struc).astype(np.uint8) # pylint: disable=no-member
nb.Nifti1Image(mask, imnii.get_affine(), hdr).to_filename(out_file)
return out_file
def spikes_mask(in_file, in_mask=None, out_file=None):
"""
Utility function to calculate a mask in which check
for :abbr:`EM (electromagnetic)` spikes.
"""
import os.path as op
import nibabel as nb
import numpy as np
from nilearn.image import mean_img
from nilearn.plotting import plot_roi
from scipy import ndimage as nd
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_spmask{}'.format(fname, ext))
out_plot = op.abspath('{}_spmask.pdf'.format(fname))
in_4d_nii = nb.load(in_file)
orientation = nb.aff2axcodes(in_4d_nii.affine)
if in_mask:
mask_data = nb.load(in_mask).get_data()
a = np.where(mask_data != 0)
bbox = np.max(a[0]) - np.min(a[0]), np.max(a[1]) - \
np.min(a[1]), np.max(a[2]) - np.min(a[2])
longest_axis = np.argmax(bbox)
# Input here is a binarized and intersected mask data from previous section
dil_mask = nd.binary_dilation(
mask_data, iterations=int(mask_data.shape[longest_axis]/9))
rep = list(mask_data.shape)
rep[longest_axis] = -1
new_mask_2d = dil_mask.max(axis=longest_axis).reshape(rep)
rep = [1, 1, 1]
rep[longest_axis] = mask_data.shape[longest_axis]
new_mask_3d = np.logical_not(np.tile(new_mask_2d, rep))
else:
new_mask_3d = np.zeros(in_4d_nii.shape[:3]) == 1
if orientation[0] in ['L', 'R']:
new_mask_3d[0:2, :, :] = True
new_mask_3d[-3:-1, :, :] = True
else:
new_mask_3d[:, 0:2, :] = True
new_mask_3d[:, -3:-1, :] = True
mask_nii = nb.Nifti1Image(new_mask_3d.astype(np.uint8), in_4d_nii.get_affine(),
in_4d_nii.get_header())
mask_nii.to_filename(out_file)
plot_roi(mask_nii, mean_img(in_4d_nii), output_file=out_plot)
return out_file, out_plot
