def _cvt_equal_mass(x, y, signal, noise, xnode, ynode, pixelsize, quiet, sn_func, wvt):
"""
Implements the modified Lloyd algorithm
in section 4.1 of Cappellari & Copin (2003).
NB: When the keyword WVT is set this routine includes
the modification proposed by Diehl & Statler (2006).
"""
dens2 = (signal/noise)**4 # See beginning of section 4.1 of CC03
scale = np.ones_like(xnode) # Start with the same scale length for all bins
for it in range(1, xnode.size): # Do at most xnode.size iterations
xnode_old, ynode_old = xnode.copy(), ynode.copy()
classe = voronoi_tessellation(x, y, xnode, ynode, scale)
# Computes centroids of the bins, weighted by dens**2.
# Exponent 2 on the density produces equal-mass Voronoi bins.
# The geometric centroids are computed if WVT keyword is set.
#
good = np.unique(classe)
if wvt:
for k in good:
index = np.flatnonzero(classe == k) # Find subscripts of pixels in bin k.
xnode[k], ynode[k] = np.mean(x[index]), np.mean(y[index])
sn = sn_func(index, signal, noise)
scale[k] = np.sqrt(index.size/sn) # Eq. (4) of Diehl & Statler (2006)
else:
mass = ndimage.sum(dens2, labels=classe, index=good)
xnode = ndimage.sum(x*dens2, labels=classe, index=good)/mass
ynode = ndimage.sum(y*dens2, labels=classe, index=good)/mass
diff2 = np.sum((xnode - xnode_old)**2 + (ynode - ynode_old)**2)
diff = np.sqrt(diff2)/pixelsize
if not quiet:
print('Iter: %4i Diff: %.4g' % (it, diff))
if diff < 0.1:
break
# If coordinates have changed, re-compute (Weighted) Voronoi Tessellation of the pixels grid
#
if diff > 0:
classe = voronoi_tessellation(x, y, xnode, ynode, scale)
good = np.unique(classe) # Check for zero-size Voronoi bins
# Only return the generators and scales of the nonzero Voronoi bins
return xnode[good], ynode[good], scale[good], it
#-----------------------------------------------------------------------
评论列表
文章目录
