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
评论列表
文章目录