我们从Python开源项目中，提取了以下27个代码示例，用于说明如何使用scipy.ndimage.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 __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 _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
