methods: data analysis – methods: numerical – techniques: image processing – techniques: imaging spectroscopy – software: data analysis – galaxies: kinematics and dynamics

1 Introduction↩︎

Modern astronomical surveys, particularly those using integral-field spectroscopy (IFS), produce vast, spatially-resolved datasets containing millions of spectra [1]–[3], Sanchez2012?. A common task is to fit these data with complex physical models to extract quantities like stellar kinematics [4], star formation histories, or chemical abundances [5]–[7]. However, the signal-to-noise ratio \(\mathcal{S/N}\) of individual spatial pixels (spaxels) is often too low for reliable model fitting [8]. Fitting a non-linear model to low-\(\mathcal{S/N}\) data typically yields a highly non-Gaussian posterior probability distribution for the model parameters, making it difficult to derive meaningful best-fitting values and uncertainties from Bayesian methods [9]. Crucially, simply averaging the biased results from low-\(\mathcal{S/N}\) fits does not recover the true parameters, because the posterior does not need to be symmetric around the true parameters.

The most effective solution is to combine the data from adjacent spaxels into larger bins to reach a sufficient \(\mathcal{S/N}\) before performing the model fit. This process, known as adaptive binning, is a crucial preprocessing step in the analysis of IFS data and other 2D datasets like X-ray images [10], [11].

1.1 The Voronoi-Binning Method↩︎

To address this need, [8] introduced the Voronoi-binning algorithm, implemented, e.g., in the VorBin package². This method has become a standard tool in astrophysics for partitioning 2D data to meet three key criteria for an optimal binning scheme: (i) the bins must form a complete, non-overlapping tessellation of the data; (ii) the bins should be as compact (round) as possible to preserve spatial resolution; and (iii) a user-defined scalar function of each bin should be as uniform as possible around a target value.

This optimized function is entirely general: while it is often chosen to be the bin’s \(\mathcal{S/N}\), the algorithm places no restriction on its form. It may represent, for example, the fractional error in a physical parameter derived from a spectral fit to the bin’s data, or a composite metric such as a weighted combination of \(\mathcal{S/N}\) values measured in different photometric bands over the same bin area.

The algorithm achieved this through a two-stage process [8]. First, a bin-accretion stage provides an initial tessellation satisfying the target \(\mathcal{S/N}\). This greedy algorithm starts with the unbinned pixel with the highest \(\mathcal{S/N}\), accretes its nearest neighbours until the target \(\mathcal{S/N}\) is reached, and repeats this process until all pixels are binned. Second, an iterative regularization stage improves the bin morphology using a Centroidal Voronoi Tessellation [12]. In each iteration, a Voronoi tessellation is generated from the current bin generators, and these generators are then updated to be the new centroids of their corresponding cells, making the bins more compact.

1.2 Limitations of Existing Methods↩︎

Despite its widespread success, the evolution of astronomical surveys has revealed some limitations of the original method and its subsequent extensions.

• Non-Convexity: An important extension by [13], which is included in VorBin, introduced a multiplicatively-weighted Voronoi tessellation to improve \(\mathcal{S/N}\) uniformity and bin shapes. While effective, the adopted tessellation sacrifices a key morphological property: the guarantee of convex bins of the original CVT method. Although not a common occurrence, this can result in undesirable non-convex or even disconnected bins, which can be problematic in certain situations.

• Computational Speed: The original implementation was not designed for the massive and rapidly growing datasets of modern and future astronomical surveys. Instruments like MUSE [14], employed in large surveys [15]–[18], can already produce mosaics of millions of pixels—for example, the 3 and 9 million spaxel mosaics of M42 and NGC 253 by [19] and [20] respectively—and this trend is set to accelerate. Next-generation concepts like the Wide-field Spectroscopic Telescope (WST) [21] are projected to produce datasets orders of magnitude larger still. Both the bin-accretion and the iterative regularization stages of VorBin scale poorly with the number of input spaxels, creating a significant computational bottleneck. In particular, the multiplicatively-weighted Voronoi diagram has a fundamental time complexity of \(\mathcal{O}(n^2)\), where \(n\) is the number of bins [22], making it impractical for large \(n\), and the same is true for the bin-accretion algorithm.

The goal of this paper is to introduce PowerBin, a fast and robust algorithm that addresses the limitations of previous adaptive-binning methods. I recast adaptive binning as a semi-discrete optimal-transport, or data-quantization, problem whose solutions are Centroidal Power Diagrams [23]–[27]. Building on a simple geometric/physical insight, PowerBin iteratively adjusts the power-diagram weights to enforce per-bin capacity targets while keeping cells convex and compact. The resulting scheme is computationally efficient, stable in the presence of realistic, non-additive capacity measures (for example, when pixel noise is correlated), and scales to the large datasets produced by modern astronomical surveys. A public Python reference implementation accompanies this paper.

This paper is structured as follows. 2 provides practical examples illustrating the necessity of binning. 3 reviews the family of weighted Voronoi diagrams. 4 introduces the optimal transport framework and its connection to Centroidal Power Diagrams. 5 presents the core physical insight behind my new fast regularization algorithm. 6 describes the new fast bin-accretion algorithm. 7 demonstrates the performance of the new method on real and simulated data. 8 presents execution-time benchmarks. Finally, 9 summarizes my findings.

2 Examples Illustrating the Need for Binning↩︎

While 1 described the general motivation for binning, the importance of this preprocessing step is best understood through practical examples. The core issue is that fitting complex, non-linear models to low-\(\mathcal{S/N}\) data can lead to results that are not just uncertain, but systematically biased. This section presents two common scenarios in spectral analysis that demonstrate this effect and motivate the need for an optimal binning strategy. The behaviour shown is generic and applies to many types of data analysis in other fields. An ideal binning scheme should partition the data according to three criteria [8]: (1) a topological criterion, ensuring all data are used without overlap; (2) a morphological criterion, requiring bins to be as compact (round) as possible to preserve spatial resolution; and (3) a uniformity criterion, minimizing the \(\mathcal{S/N}\) scatter around a target value. The following examples highlight why achieving a target \(\mathcal{S/N}\) is paramount.

2.1 Extracting Stellar Kinematics of Multiple Components↩︎

This example examines the problem of separating the kinematics of multiple stellar populations along a single line of sight. The goal is to illustrate a generic characteristic of non-linear problems in general, and astrophysical spectrum fitting in particular: at low \(\mathcal{S/N}\) the model posterior is not informative of the true solution. For this reason, the choice of spectra and kinematics is rather arbitrary and is not meant to represent a specific situation. I construct a synthetic spectrum made from two distinct stellar components and attempt to recover their velocities and dispersions at two representative per-pixel signal-to-noise levels, showing that low \(\mathcal{S/N}\) yields multimodal, biased parameters while higher \(\mathcal{S/N}\) permits reliable parameter recovery.

For this test, I took two model spectra from the MILES library [28], both with solar metallicity but with different ages (12.6 Gyr and 1.0 Gyr). The spectra were logarithmically rebinned to a velocity scale of 70 \(\text{km s}^{-1}\), typical of large surveys like SDSS [29], and normalized to contribute equally to the total flux in the fitted region (354–741 nm). I assigned distinct kinematics to each component: \((V_1, \sigma_1) = (-150, 100)\) \(\text{km s}^{-1}\) for the old population and \((V_2, \sigma_2) = (150, 200)\) \(\text{km s}^{-1}\) for the young one.

I then used the Penalized PiXel Fitting (ppxf) software³ [30]–[32] to recover the four kinematic parameters \((V_1, \sigma_1, V_2, \sigma_2)\). To explore the parameter space, I coupled ppxfwith the Adaptive Metropolis MCMC sampler by [33] as implemented in the AdaMet package⁴ [34]. Assuming a uniform prior, the posterior probability is \(P\propto\exp(-\chi^2/2)\). I ran a chain of \(10^5\) steps for two cases: a low-\(\mathcal{S/N}\) case (\(\mathcal{S/N}=2\) per pixel) and a high-\(\mathcal{S/N}\) case (\(\mathcal{S/N}=20\)).

The results are shown in 1. At \(\mathcal{S/N}=2\) (left panel), the posterior distribution is highly complex, non-Gaussian, and biased away from the true input values (magenta lines). The posterior is hard to interpret and one cannot meaningfully quote a single best-fitting value or its error. Averaging such biased results from many low-\(\mathcal{S/N}\) spaxels would not recover the true average kinematics. In contrast, at \(\mathcal{S/N}=20\) (right panel), the posterior is unimodal, symmetric, and correctly centred on the true values. The right panel unambiguously reveals the two kinematic components and the parameters are well-constrained. This demonstrates that reaching a sufficient \(\mathcal{S/N}\) by binning is essential before attempting such a measurement.

2.2 Star Formation History from Full-Spectrum Fitting↩︎

My second example concerns the recovery of a galaxy’s star formation history (SFH), a problem that involves fitting a spectrum with a large linear combination of template spectra representing stellar populations of different ages and metallicities.

Here, I constructed a synthetic galaxy spectrum with a simple SFH, dominated by a single burst of star formation 0.3 Gyr ago. The mock spectrum was built from a linear combination of 25 solar-metallicity templates from the MILES models [28], with ages spaced logarithmically between 0.063 and 15.8 Gyr. I then attempted to recover the SFH using ppxfat three different \(\mathcal{S/N}\) levels: 5, 50, and 500. For each \(\mathcal{S/N}\) level, I ran 1000 Monte Carlo realisations, adding appropriate Gaussian noise to the spectrum in each run and fitting for the template weights. The SFH recovery was performed with a low level of regularization (‘regul=10’; see [32, Fig. 5] for visual reference). This example is not intended to illustrate the expected performance of ppxfin a realistic setting; for that, the reader is referred to [35]. Rather, its purpose is to demonstrate a generic feature of non-linear models, particularly those used in astrophysical problems where binning is common: at low \(\mathcal{S/N}\), the model posteriors can be biased and unrepresentative of the true solution.

2 summarizes the Monte Carlo SFH recoveries. At \(\mathcal{S/N}=5\) the solutions are dominated by noise and a clear systematic bias: the fitter spuriously assigns weight to old populations and fails to recover the true single-burst history. At \(\mathcal{S/N}=50\) the correct peak is recovered but with substantial scatter, while only at \(\mathcal{S/N}=500\) does the ensemble reliably reproduce the input SFH with high fidelity. These tests demonstrate that reaching a problem-dependent minimum \(\mathcal{S/N}\) is essential for trustworthy spectral inference; adaptive binning is therefore a necessary preprocessing step, not merely a cosmetic choice.

3 Generalizations of Voronoi Diagrams↩︎

The original Voronoi-binning algorithm [8] used the ordinary Voronoi tessellation, the simplest member of a broader family of weighted Voronoi diagrams. Later, [13] extended the regularization phase by adopting the multiplicatively weighted variant. In order to motivate the adoption of a different tessellation for adaptive binning, it is useful to place these and other variants side-by-side. This section therefore reviews the principal types of weighted Voronoi diagrams, following the classic treatments of [36] and [37], and uses 3 to give a geometric interpretation of each. This comparison will make clear that one particular member of the family occupies a special position for our purposes.

3.1 Ordinary Voronoi Diagram↩︎

Given a set of \(n\) distinct points \(\mathbf{G}=\{\mathbf{g}_1,\ldots,\mathbf{g}_n\}\) in a \(d\)-dimensional Euclidean space \(\mathbb{R}^d\), called generators, the ordinary Voronoi diagram is a partition of the space into regions based on the nearest-neighbour rule. The Voronoi cell \(\mathcal{V}(\mathbf{g}_j)\) associated with a generator \(\mathbf{g}_j\) contains all points in \(\mathbb{R}^d\) that are closer to \(\mathbf{g}_j\) than to any other generator \(\mathbf{g}_k\). Using the Euclidean distance \[d(\mathbf{x},\mathbf{g}_j)=\lVert\mathbf{x} - \mathbf{g}_j\rVert\] the cell is defined as: \[\mathcal{V}(\mathbf{g}_j)=\{\mathbf{x} \mid \lVert\mathbf{x} - \mathbf{g}_j\rVert\leq \lVert\mathbf{x} - \mathbf{g}_k\rVert, \forall k \neq j\}.\] The boundaries of the Voronoi cells are formed by segments of perpendicular bisectors between pairs of generators. Consequently, the cells are always convex polytopes. Efficient algorithms exist for computing ordinary Voronoi diagrams, with a worst-case time complexity of \(\mathcal{O}(n\log n)\) in 2D [37].

Geometrically, the ordinary Voronoi diagram can be visualized as the projection onto the 2D plane of the intersections of a set of identical 3D cones, whose apices are located at the generator positions (3a). The diagram can be seen by looking at the cones from above.

3.2 Multiplicatively Weighted Voronoi Diagram↩︎

The multiplicatively weighted Voronoi (MWV) diagram assigns a weight \(w_j > 0\) to each generator \(\mathbf{g}_j\), with the weight acting as a scaling factor on the distance. This allows one to change the influence of each generator on the partitioning of space. The weighted distance is \[d_{\rm MW}(\mathbf{x}, \mathbf{g}_j) = \lVert\mathbf{x} - \mathbf{g}_j\rVert / w_j.\] The cell is defined as: \[\mathcal{V}_{\rm MW}(\mathbf{g}_j)=\{\mathbf{x} \mid \lVert\mathbf{x} - \mathbf{g}_j\rVert/w_j \leq \lVert\mathbf{x} - \mathbf{g}_k\rVert/w_k, \forall k \neq j\}.\] The boundary between two cells, \(\lVert\mathbf{x} - \mathbf{g}_j\rVert/\lVert\mathbf{x} - \mathbf{g}_k\rVert = w_j/w_k\), is a hypersphere (known as a Circle of Apollonius in 2D). A key characteristic of MWV diagrams is that their cells are not necessarily convex and in general can be disconnected. 3 (b) shows an example of a MWV diagram with a cell contained inside another. Crucially, the optimal computation time for its generation scales as \(\mathcal{O}(n^2)\) [22], [37].

The geometric interpretation of the boundaries of a MWV diagram is the projection of the intersections of cones with different inclinations (slopes), determined by their weights \(w_j\) (3b).

3.3 Additively Weighted Voronoi Diagram↩︎

The additively weighted Voronoi (AWV) diagram also assigns a real-valued weight \(w_j\) to each generator \(\mathbf{g}_j\). The distance metric is modified by subtracting this weight, defining a weighted distance \[d_{\rm AW}(\mathbf{x}, \mathbf{g}_j) = \lVert\mathbf{x} - \mathbf{g}_j\rVert - w_j.\] The corresponding cell is: \[\mathcal{V}_{\rm AW}(\mathbf{g}_j)=\{\mathbf{x} \mid \lVert\mathbf{x} - \mathbf{g}_j\rVert - w_j \leq \lVert\mathbf{x} - \mathbf{g}_k\rVert - w_k, \forall k \neq j\}.\] The boundary between two cells, defined by \(\lVert\mathbf{x} - \mathbf{g}_j\rVert - \lVert\mathbf{x} - \mathbf{g}_k\rVert = w_j - w_k\), is a sheet of a hyperboloid of revolution. In 2D, the boundaries are hyperbolic arcs. Like for MWV diagrams, the cells of AWV diagrams are not guaranteed to be convex. The computational complexity is like for the ordinary Voronoi diagrams, scaling as \(\mathcal{O}(n \log n)\) [37], [38].

This diagram can be visualized as the projection of the intersections of cones with identical slopes but different heights, where the apex of each cone is shifted vertically by its weight \(w_j\) (3c).

3.4 Power Diagram: The Ideal Candidate↩︎

Among the family of weighted Voronoi diagrams, the power diagram (also known as the Laguerre-Voronoi diagram) is not just another variation; it possesses a unique combination of geometric and computational properties that make it the ideal mathematical structure for the adaptive binning problem. It is defined using the ‘power distance’, where each generator \(\mathbf{g}_j\) is associated with a real-valued weight \(w_j\). The power of a point \(\mathbf{x}\) with respect to a generator \(\mathbf{g}_j\) is given by: \[\mathrm{pow}(\mathbf{x}, \mathbf{g}_j) = \lVert\mathbf{x} - \mathbf{g}_j\rVert^2 - w_j.\] The power cell \(\mathcal{V}_{\rm pow}(\mathbf{g}_j)\) consists of all points whose power with respect to \(\mathbf{g}_j\) is less than or equal to their power with respect to any other generator: \[\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\{\mathbf{x} \mid \lVert\mathbf{x} - \mathbf{g}_j\rVert^2 - w_j \leq \lVert\mathbf{x} - \mathbf{g}_k\rVert^2 - w_k, \forall k \neq j\}.\] The boundary condition, \(\lVert\mathbf{x} - \mathbf{g}_j\rVert^2 - \lVert\mathbf{x} - \mathbf{g}_k\rVert^2 = w_j - w_k\), simplifies to a linear equation. This means the boundaries are hyperplanes (straight lines in 2D), which in turn guarantees that the cells are always convex polytopes—a critical property that both AWV and MWV diagrams lack.

Furthermore, this linear boundary property allows power diagrams to be computed with optimal efficiency. By transforming the problem into a convex hull construction in one higher dimension, the tessellation can be found in \(\mathcal{O}(n\log n)\) time in 2D [23], [39], matching the speed of the simplest ordinary Voronoi diagram.

Geometrically, the boundaries of a power diagram are the vertical projection of the intersections of a set of upward-opening paraboloids \(z = \lVert\mathbf{x} - \mathbf{g}_j\rVert^2 - w_j\). Equivalently, the boundaries of a 2D power diagram can be viewed as the vertical projection of the intersection of a set of 3D spheres with different radii, whose centres lie on the \(z=0\) 2D plane (3d). To see this, consider a sphere \(j\) centred at \((\mathbf{g}_j, 0)\) with radius \(r_j\). Its equation is \(\lVert\mathbf{x} - \mathbf{g}_j\rVert^2 + z^2 = r_j^2\). The intersection of two spheres \(j\) and \(k\) lies on a plane (the radical plane) defined by equating their equations, which yields \(\lVert\mathbf{x} - \mathbf{g}_j\rVert^2 - r_j^2 = \lVert\mathbf{x} - \mathbf{g}_k\rVert^2 - r_k^2\). This is precisely the boundary condition for a power diagram with weights \(w_j = r_j^2\). This sphere-based interpretation provides a direct link to the physical analogy of packed soap bubbles (4), which inspires the fast heuristic algorithm presented in this paper. As I will show in 4, this specific geometric form makes the power diagram the natural solution to the problem of optimal transport, providing a rigorous mathematical foundation for capacity-constrained binning. This unique combination of guaranteed convexity, computational efficiency, and a direct link to optimal transport theory sets the power diagram apart as the ideal choice for our application.

4 Centroidal Power Diagrams for Optimal-Transport Binning↩︎

In this section, I show that among the various generalizations of Voronoi diagrams, power diagrams are uniquely suited for the adaptive binning of empirical data. This is because they provide a natural solution to a class of problems known as optimal transport, which offers a rigorous mathematical foundation for the binning criteria I outlined in 1.

4.1 The Optimal Transport Problem↩︎

The theory of optimal transport, first formulated by [40], provides a mathematical framework for finding the most efficient way to remap one distribution of ‘mass’ (or any density) to another, given a specified transport cost. For comprehensive reviews of the theory and its computational methods, see [41] and [42].

For the adaptive binning problem, I consider the ‘semi-discrete’ case: transporting a continuous density distribution \(\rho(\mathbf{x})\), which is approximated by our \(N\) pixels of data, to a discrete set of \(n\) target locations, the bin generators \(\{\mathbf{g}_j\}_{j=1}^n\). I define the transport cost as the squared Euclidean distance, \(\lVert\mathbf{x} - \mathbf{g}_j\rVert^2\). The goal is to find a partition of the data into a set of \(n\) bins \(\{\mathcal{V}_j\}\) that minimizes the total transport cost, \[\mathcal{E}(\{\mathbf{g}_j\}, \{\mathcal{V}_j\}) = \sum_{j=1}^n \int_{\mathcal{V}_j} \lVert\mathbf{x}-\mathbf{g}_j\rVert^2\,\rho(\mathbf{x})\,\mathrm{d}\mathbf{x}, \label{eq:primal95energy}\tag{1}\] while ensuring each bin \(j\) contains a prescribed amount of mass, or ‘capacity’, \(\nu_j\). That is, the partition must satisfy the constraint \[m_j = \int_{\mathcal{V}_j} \rho(\mathbf{x})\,\mathrm{d}\mathbf{x} = \nu_j, \quad \forall j. \label{eq:capacity95constraint}\tag{2}\] It can be shown that the optimal partition \(\{\mathcal{V}_j\}\) for this problem is a power diagram [24].

4.2 Energy Functional and Centroidal Power Diagrams↩︎

While the primal energy \(\mathcal{E}\) in 1 is intuitive, finding the partition \(\{\mathcal{V}_j\}_{j=1}^n\) that minimizes it under capacity constraints is difficult. The problem becomes tractable by considering the dual problem, which can be expressed using a Lagrangian functional \(\mathcal{F}\) that depends on the generator positions \(\{\mathbf{g}_j\}\) and a set of real-valued weights \(\{w_j\}\) acting as Lagrange multipliers [24]–[27]: \[\mathcal{F}(\{\mathbf{g}_j\}, \{w_j\}) = \mathcal{E}(\{\mathbf{g}_j\}, \{\mathcal{V}_j\}) - \sum_{j=1}^n w_j (m_j - \nu_j). \label{eq:dual95energy}\tag{3}\] Here, \(\mathcal{E}\) is the primal energy from 1 , where the partition \(\{\mathcal{V}_j\}\) is now explicitly shown to be the power diagram defined by the weights \(\{w_j\}\). The term \(m_j\) is the current capacity of cell \(\mathcal{V}_j\) from 2 , and \(\nu_j\) is the target capacity, which for the Voronoi-binning problem I generally assume to be constant \(\nu\), although this is not a requirement for the method.

The key insight is that finding the optimal binning is equivalent to finding a saddle point of \(\mathcal{F}\). The gradients of \(\mathcal{F}\) with respect to the weights and generator positions reveal its utility:

1. Gradient w.r.t. weights: The gradient with respect to a weight \(w_j\) is simply the difference between the cell’s target capacity and its current capacity [24], [26]: \[\nabla_{w_j} \mathcal{F} = \nu_j - m_j.\] For a fixed set of generators, finding the weights \(\{w_j\}\) that maximize the dual functional \(\mathcal{F}\) is a convex optimization problem [24]. This is a crucial property, as it guarantees the existence of a unique global maximum. Standard and efficient algorithms, such as Newton’s method, can be used to find this solution by driving the gradient to zero, thus ensuring the capacity constraints \(m_j = \nu_j\) are satisfied.

2. Gradient w.r.t. generators: The gradient with respect to a generator position \(\mathbf{g}_j\) is [26]: \[\nabla_{\mathbf{g}_j} \mathcal{F} = 2 m_j (\mathbf{g}_j - \mathbf{b}_j),\] where \(\mathbf{b}_j\) is the barycentre (density-weighted centroid) of the cell \(\mathcal{V}_j\). Setting this gradient to zero implies that the generator must coincide with its cell’s barycentre: \(\mathbf{g}_j = \mathbf{b}_j\).

A configuration that is a stationary point of \(\mathcal{F}\)—simultaneously satisfying the capacity constraints and the barycentric condition—is called a Centroidal Power Diagram (CPD). A CPD corresponds to a (local) minimum of the original transport energy \(\mathcal{E}\), thus providing a complete and principled solution to the adaptive binning problem.

4.3 Challenges in Applying to Astronomical Data↩︎

While the CPD framework is theoretically ideal, its direct application to the binning of astronomical data faces two main practical challenges:

1. Discrete Data vs. Continuous Theory: The optimal transport theory is formulated for continuous density functions. When applied to discrete data, integrals are replaced by sums over pixels. This approximation is valid when bins are large, but it breaks down when bins contain only a few pixels, leading to numerical instabilities. While one could interpolate the discrete data to create a continuous density, this is only rigorously applicable when the capacity function is additive.

2. Non-Additive Capacity: The most significant challenge is that the bin ‘capacity’ is often not a simple additive quantity. For example, when binning to a target signal-to-noise ratio (\(\mathcal{S/N}\)), the bin’s total \((\mathcal{S/N})^2\) is only the sum of the pixel \((\mathcal{S/N})^2\) if the noise is uncorrelated [8]. In practice, instrumental effects and data reduction steps (e.g., dithering, resampling) introduce significant covariance between pixels [4]. This makes the capacity \(m_j\) a non-linear, non-additive function of its constituent pixels. As a consequence, the dual functional \(\mathcal{F}\) loses its convenient convexity at fixed generator positions, and the analytic gradients become invalid. As a result, standard gradient-based optimization methods become unstable and fail to converge.

I confirmed this limitation through numerical experiments using the formalisms of [24] and [26]. While these variational approaches perform well for additive capacities in the continuum limit (i.e., large bins with many pixels), they fail catastrophically for the non-additive capacities typical of real data with correlated noise.

Because of these issues, a direct implementation of a mathematically exact CPD solver is not robust for real-world data binning. In the following section, I introduce a new algorithm, PowerBin, which is inspired by the optimal transport framework but uses a fast and robust heuristic to handle these complexities.

5 Fast centroidal power-diagram solver↩︎

The previous section established that while Centroidal Power Diagrams (CPDs) provide a theoretically ideal framework for adaptive binning, formal solvers based on gradient-based optimization of the dual energy functional are impractical for real astronomical data. The non-additive nature of capacity measures like \(\mathcal{S/N}\) with correlated noise violates the assumptions required for these methods to converge reliably. This section introduces the core of the PowerBin algorithm: a fast, robust, and physically-motivated heuristic that bypasses these problems. Instead of relying on complex and fragile numerical optimization, I develop a simple iterative scheme inspired by the geometry of packed cells, which proves highly effective at enforcing capacity constraints while maintaining computational efficiency and bin convexity.

5.1 A physical heuristic for the weight–capacity relation↩︎

The central challenge in constructing a capacity-constrained power diagram is to find the set of weights \(\{w_j\}\) that yields cells with the desired capacities \(\{\nu_j\}\). As noted by experts in the field, ‘the relation between the weights and the measures of the power cells is non-trivial’ [27]. This complexity arises from two main factors. First, the weights are defined only up to a common additive constant; adding the same value to all weights leaves the tessellation unchanged. Second, for an arbitrary arrangement of generators, there is no simple, direct relationship between a weight \(w_j\) and the area \(A_j\) of its corresponding cell. The cells can be highly elongated, and some generators may even have empty cells.

However, the problem simplifies dramatically if one considers the specific geometry of a centroidal tessellation, where each generator is close to the centre of its cell. My approach is based on two key insights:

1. Power weights as squared radii. As shown in 3, a power diagram can be defined by associating a circle of radius \(r_j\) with each generator \(\mathbf{g}_j\), such that the weight is \(w_j = r_j^2\). This provides a direct geometric interpretation: the cell boundaries are the radical axes of these circles, and their sizes control the tessellation.

2. The ‘packed bubbles’ approximation. In a centroidal configuration, the generators are close to the cell centroids, which naturally produces compact, nearly-round cells. In this limit, the tessellation resembles a foam of packed bubbles (4). For such a packing, the area \(A_j\) of a cell is well-approximated by the area of its defining circle: \[A_j \approx \pi r_j^2. \label{eq:area95approx}\tag{4}\] This simple approximation provides the crucial link between a cell’s area (\(A_j\)) and its generator’s weight (\(w_j = r_j^2\)).

With this physical model, I can derive a simple update rule. The goal is to adjust the radii \(\{r_j\}\) until the measured capacity \(m_j\) of each cell matches a target value \(\nu\). I start by assuming that a cell’s capacity is, to first order, proportional to its area: \(m_j \approx \rho_j A_j\), where \(\rho_j\) is an effective local capacity density.

To achieve the target capacity \(\nu\), cell \(j\) would need a target area \(A_j^{\star} \approx \nu / \rho_j\). Using our approximation from 4 , the corresponding target squared radius would be \(r_j^{2\,\star} \approx A_j^{\star}/\pi = \nu/(\pi\,\rho_j)\). The unknown density \(\rho_j\) can be eliminated by substituting its value from the current iteration, \(\rho_j \approx m_j/A_j\). This yields a simple update rule for the target squared radius based entirely on measurable quantities from the current tessellation: \[r_j^{\text{new}} \leftarrow \sqrt{\frac{\nu}{m_j} \frac{A_j}{\pi}} = \sqrt{\frac{f_j A_j}{\pi}}, \quad \text{where} \quad f_j \equiv \frac{\nu}{m_j}. \label{eq:multiplicative95rule}\tag{5}\]

In terms of the power weights themselves, the update is \(w_j^{\text{new}} \leftarrow f_j A_j / \pi\).

This leads to the simple yet powerful iterative algorithm summarized in 5. The update rule in 5 is formally similar to the WVT rule of [13], but it operates on the radii of a power diagram, not the weights of an MWV diagram. This is a crucial distinction: for MWV diagrams, multiplying all weights by a common factor has no effect, whereas for a power diagram, it changes the tessellation. This implies that for a MWV the constant factors are irrelevant, while for a Power Diagram they are essential. By combining the natural area–radius relation (\(A\approx \pi r^2\)) with centroidal recentring, this heuristic achieves stable, few-iteration convergence for capacity equalization while preserving the convexity and compactness of the bins.

5.2 Implementation Details↩︎

The success and speed of the regularization stage of the PowerBin algorithm rely on a few crucial implementation choices, which are detailed in 5, 6 and 7.

Efficient Power Tessellation:

A key advantage of power diagrams over other weighted Voronoi diagrams is their computational efficiency. A power diagram can be computed with \(\mathcal{O}(N\log N)\) complexity [23], whereas an MWV diagram requires \(\mathcal{O}(N^2)\) operations [22]. I achieve this efficiency for my discrete dataset by implementing the geometric lifting technique described in [39], as detailed in 7. Each 2D generator \(\mathbf{g}_j\) with radius \(r_j\) is ‘lifted’ to a 3D point \((\mathbf{g}_j, z_j)\) where \(z_j^2 = r_{\max}^2 - r_j^2\). The power-diagram assignment for a 2D spaxel \(\mathbf{x}_i\) is then equivalent to finding the nearest 3D neighbour to the point \((\mathbf{x}_i, 0)\) among the lifted generators. This 3D nearest-neighbour search is performed efficiently using a standard KD-Tree (scipy.spatial.KDTree) in the SciPy library [43], which implements the algorithm of [44].

Geometric Centroids vs. Barycentres:

A key choice in my algorithm is the update of the generators. Instead of moving them to the capacity-weighted barycentre of each cell, I update them to the cell’s unweighted, geometric centroid (the mean position of its constituent pixels). This means I do not strictly minimize the optimal transport energy functional from 1 . However, this choice is deliberate and crucial for robustness. The main reason is that it allows the method to handle data with negative values (e.g., background-subtracted X-ray data), where the definition of a barycentre breaks down. This choice was adopted in the standard VorBin algorithm for the same reason. Moreover, this choice makes the algorithm a more direct implementation of the physical ‘soap bubble’ analogy (4) that inspired my area-weight relation, as the cells are centred geometrically rather than by mass. Furthermore, for real data with non-additive capacities, the energy functional itself is ill-defined, making its formal minimization moot. For positive data, I find that my algorithm is not particularly sensitive to this choice, unlike the standard unweighted CVT which relies on the barycentric condition.

Failure of Formal Optimization:

The use of a heuristic update rule and geometric centroids is a direct consequence of the practical failure of formal optimization methods. During development, I implemented a version of the algorithm that directly minimized the energy functional from [24] and [26] using their analytic gradient and a quasi-Newton optimizer. This approach is not just slower; it fails completely. Even in the continuum limit with many spaxels per bin, the method fails to converge to a sensible result (i.e., the capacity is not equalized) when the capacity function is non-additive, as is generally the case with real data. The fast, physically-motivated heuristic of PowerBin proved to be far superior in all practical tests, converging quickly and reliably to a high-quality solution where formal methods could not.

Robust Convergence:

Iterative methods can sometimes stall or enter a cycle. To prevent this, I use under-relaxation when I detect cycling, and employ a robust early-stopping heuristic. This method monitors the sequence of generator shifts and terminates the loop if the improvement stagnates or if it detects that the values are oscillating without any significant downward trend. This ensures reliable termination even in difficult cases.

Generality of Capacity Function:

The combination of my robust heuristic and the use of geometric centroids makes the PowerBin algorithm extremely versatile. The capacity function is not limited to the specific forms discussed in this paper and can be adapted to various applications. This flexibility allows for the incorporation of complex, non-additive constraints, making the algorithm applicable to a broad range of problems.

6 A linear-time bin-accretion algorithm↩︎

A crucial, and perhaps under-appreciated, aspect of all successful adaptive-binning schemes is the quality of the initial tessellation. The iterative refinement stages, whether based on a Centroidal Voronoi Tessellation (CVT), a Weighted Voronoi Tessellation (WVT), or the Centroidal Power Diagram (CPD) presented here, are all local optimizers. They are variants of [45] algorithm, which is known to be sensitive to the initial placement of the generators. Unlike the optimization of weights at fixed generator’s location, the optimization of the energy functional of 3 with respect to the generator positions is not convex and presents a large number of secondary minima [27]. This means that the final result can vary significantly depending on the starting configuration.

One cannot, for instance, initialize the generators with points drawn randomly from the underlying signal or \(\mathcal{S/N}\) distribution and expect the iterations to converge to a satisfactory result. The discrete nature of the data and the non-convexity of the optimization landscape mean that such an approach will invariably become trapped in a poor local minimum, yielding a tessellation with suboptimal bin shapes and poor capacity uniformity. Consequently, the bin-accretion algorithm, first introduced in [8], has always been the indispensable foundation of the entire procedure. It provides an excellent initial guess that already satisfies the capacity constraint, allowing the subsequent refinement to focus solely on improving the bin morphology.

With the development of the fast CPD solver, which has a time complexity of \(\mathcal{O}(N\log N)\) for \(N\) pixels, the original bin-accretion algorithm became the computational bottleneck. To fully realize the performance gains of the new method, it was essential to devise an accretion algorithm with a comparable, near-linear time complexity. I achieve this through four key improvements, while the complete revised method is described in 8:

1. Delaunay Adjacency: I begin by pre-computing a single Delaunay triangulation of all input pixel coordinates. This is an \(\mathcal{O}(N\log N)\) operation that provides a static adjacency graph for the entire dataset. For any given pixel, its neighbours are instantly known without requiring any further geometric searches. The computation is done using scipy.spatial.Delaunay, which is based on the QHull library [46].

2. Frontier-Based Growth: To efficiently find the next pixel to add, the algorithm restricts its search to the bin’s ‘frontier’. This frontier is defined as the set of unbinned pixels that are Delaunay neighbours to any pixel already in the current bin. By only considering this small subset of pixels at each step, the search for the one closest to the bin’s centroid becomes extremely fast.

3. Incremental Updates: All quantities required to assess the validity of adding a new pixel to a bin—namely its centroid and mean squared radius—are updated incrementally using [47] algorithm. The total capacity is also updated incrementally if it is additive; for an additive capacity, this simply means adding the capacity of the new pixel to the running total for the bin. A recomputation for the entire bin is still required in the case of a black-box generic capacity function. Adding a pixel involves a simple update to a running sum, an \(\mathcal{O}(1)\) operation, rather than a full re-computation over all pixels in the growing bin.

4. Heap-Managed Seeding: To efficiently select the starting point for each new bin, I use a max-heap data structure containing all unbinned pixels, prioritized by their initial capacity density (or ‘brightness’, or \(\mathcal{S/N}\)). After a bin is completed, the next seed is retrieved from the top of the heap in an \(\mathcal{O}(\log N)\) operation. This approach is more efficient than the original method of [8], which started a new bin from the unbinned pixel closest to the centroid of all remaining unbinned pixels. The brightness of a pixel is a static property, allowing for a single heap construction at the start, whereas the centroid of unbinned pixels changes after every bin is created, requiring a costly recalculation.

In many common astronomical applications, such as binning the surface brightness of a galaxy, this change in seeding strategy has a minimal impact on the final result. In such cases, the brightest unbinned pixels are naturally located near the galaxy’s centre, which is also where the centroid of the remaining light is. The old method would accrete bins in nearly spherical shells around the galaxy’s core, while the new ‘brightest-pixel’ approach will still start bins near the centre but will tend to grow them along isophotes. In general, both methods provide an excellent starting point for the subsequent regularization. However, the new approach appears superior in the case of multiple distinct objects, as illustrated in 7.2.

Apart from these significant algorithmic optimizations, the new implementation aims to reproduce the logic of the original bin-accretion algorithm from [8]. There are, however, two minor differences in the acceptance criteria. (a) I employ a different definition of roundness, based on the mean squared radius of the pixel coordinates, which is faster to update incrementally. (b) The precise conditions for accepting a new pixel into a bin have been slightly adjusted. The new algorithm therefore represents a much faster, but functionally very similar, version of the well-tested accretion method. It is worth noting that the roundness criterion is not critical for the success of the overall algorithm; the final tessellation after the regularization stage is quite similar, though not identical, even if this test is omitted entirely from the accretion phase.

7 Power Binning Examples↩︎

In this section, I demonstrate the compactness and uniformity performance, and the versatility of the PowerBin algorithm. I apply it to a range of test cases, including simulated galaxies with different morphologies and noise properties, real integral-field spectroscopic data, and a large, complex image to showcase its scalability and applicability beyond standard astronomical use cases.

7.1 Application to Mock Single Galaxies↩︎

I first test the algorithm on mock galaxy data to assess its quality performance under controlled conditions. I created two types of galaxy images: one with an exponential disk profile (Sérsic \(n_{\rm Ser}=1\)) and another with a highly concentrated, elliptical-like profile [48] described by a [49] profile (\(n_{\rm Ser}=8\)). These two cases test the algorithm’s ability to handle both shallow and steep signal gradients. For each galaxy, I consider two scenarios for the noise: uncorrelated or correlated.

The results are shown in 9. To aid visualization, the tessellation is coloured using the networkx.coloring.greedy_color algorithm from the NetworkX package [50] to approximately four-colour the Delaunay graph of the bin generators. The top two panels show the ideal case of uncorrelated, Poissonian noise. In this scenario, the bin signal-to-noise is calculated using the standard capacity function for uncorrelated noise \[\label{eq:capacity95uncorrelated} \mathcal{C}(\{i\})=(\mathcal{S/N}_\mathrm{bin})^2 = \frac{(\sum_i \mathcal{S}_i)^2}{\sum_i \mathcal{N}_i^2},\tag{6}\] where \(\mathcal{S}_i\) and \(\mathcal{N}_i\) are the signal and noise of the individual pixels in the bin. For Poissonian noise, \(\mathcal{N}_i^2 = \mathcal{S}_i\), and the bin capacity becomes \((\mathcal{S/N}_\mathrm{bin})^2 = \sum_i \mathcal{S}_i\), which is additive. The left sub-panels show the resulting power diagram tessellation. The bins are compact and convex, and their sizes adapt smoothly to the underlying signal gradient, becoming larger in the low-\(\mathcal{S/N}\) outskirts. The right sub-panels confirm that the final bin \(\mathcal{S/N}\) (blue points) clusters tightly around the target value (dashed line), demonstrating excellent uniformity.

The bottom two panels illustrate a more realistic scenario where noise covariance is present. I simulate this by making the capacity function non-additive, using an empirical formula derived from real IFS data to penalize the \(\mathcal{S/N}\) of bins with many pixels. Specifically, the bin \(\mathcal{S/N}\) is modified as \[\mathcal{S/N}_\mathrm{bin}\rightarrow \frac{\mathcal{S/N}_\mathrm{bin}}{1 + 1.07 \lg N_{\rm pix}},\] where \(N_{\rm pix}\) is the number of spaxels in the bin. This formula is not general but depends on the data under study. It was derived for CALIFA data [51] and is used here for illustrative purposes. This test highlights a key strength of my heuristic approach: despite the non-linear and non-additive nature of the capacity, the algorithm converges robustly and still produces bins with excellent \(\mathcal{S/N}\) uniformity. A formal gradient-based optimizer would fail to converge in this regime, but the physically-motivated update rule of PowerBin, combined with the bin-accretion initialization, handles it with ease.

7.2 Application to Mock Sky-Subtracted Data↩︎

A crucial consideration for modern IFS data is that they are often background-limited rather than purely Poissonian. In this regime, the noise is roughly constant across the field, meaning that accreting spaxels with very low signal can decrease a bin’s \(\mathcal{S/N}\) instead of increasing it. This raises the question of how the algorithm performs on large mosaics that may contain multiple sources of interest within a vast field of mostly ‘empty sky’.

To address this scenario, I simulated a mock galaxy group image, shown in 10. The simulation generates an image by adding four exponential profile galaxies onto a constant sky background of 100 counts. The galaxy models have peak surface brightness ranging from 500 to 1000 counts. The combined image is then subjected to Poisson noise, and a slightly imperfect sky subtraction is performed to mimic a realistic data reduction process. This creates a background-limited image where the noise is nearly constant, dominated by the sky contribution. The results demonstrate that PowerBin handles this challenging case robustly. It correctly identifies and tessellates the regions with galaxy signal, with the resulting bins clustering tightly around the target \(\mathcal{S/N}\). The empty, noise-dominated regions are efficiently grouped into very large bins, which a user would typically discard in a subsequent analysis step.

The algorithm’s success in this regime is due to two key design choices. First, the bin-accretion stage seeds each new bin from the unbinned pixel with the highest \(\mathcal{S/N}\). This allows the algorithm to effectively grow bins ‘in parallel’ around multiple, disconnected sources. This contrasts with the classic VorBin approach, which started new bins from the barycentre of all remaining pixels and could risk a large noise-bin growing to engulf a faint, distant source. Second, unlike some previous implementations of VorBin, the new accretion algorithm does not stop if adding a pixel temporarily decreases the bin’s \(\mathcal{S/N}\). Instead, it continues to accrete neighbours until the target capacity is met or the bin’s roundness criterion is violated. This makes the process more robust in the background-dominated regime.

This example suggests that one could apply PowerBin to an entire data cube and select the scientifically useful bins afterward. However, for fields with extremely sparse signals, a conservative approach of pre-selecting regions of interest before binning may still be a sensible choice.

7.3 Application on Real IFS Data↩︎

To provide a direct comparison with previous work, I apply PowerBin to the SAURON integral-field data of the galaxy NGC 2273, shown in 11. This dataset served as the primary test case for the original Voronoi-binning method and its subsequent multiplicatively-weighted extension, making it a benchmark for two decades. The application of those previous methods to this same dataset is shown in [8] and [13] respectively, and can be directly compared with my 11. The figure shows that the new algorithm performs flawlessly on this classic dataset. It produces a clean, convex tessellation that adapts to the galaxy’s morphology, and the resulting bin \(\mathcal{S/N}\) is highly uniform around the target value. The example included in the public VorBin package, using the WVT regularization, gives a fractional rms scatter of 7.3%, which is slightly larger than the 6.5% scatter produced by PowerBin on the same input data. This confirms that PowerBin successfully matches the results of the original method on real astronomical data.

This paper focuses on the algorithm and does not present new scientific applications, such as the generation of kinematic maps. This is because, while PowerBin resolves the critical issues of non-convex bins and slow computation, the resulting tessellations are visually very similar to those from the classic VorBin method. The latter has been successfully applied to hundreds of datasets over many years.

Instead of reproducing similar science results, I refer the reader to existing work for examples of the high-quality scientific products that can be derived from this binning approach. For beautiful maps of stellar kinematics and populations, particularly from high-quality MUSE data, see, for instance, [16], [52]–[55]. For the largest applications of VorBin to date on the ever-increasing samples from major IFS surveys, see the results from the ATLAS^3D, CALIFA, SAMI, and MaNGA surveys in [1], [56], [57], and [4], respectively.

7.4 General-Purpose Tessellation↩︎

Finally, to demonstrate the scalability and versatility of PowerBin on a non-astronomical problem, I apply it to a task common in computer graphics: creating a stipple drawing from an image. I took a \(512 \times 512\) pixel greyscale self-portrait (12) and tasked the algorithm with partitioning it into \(10^4\) bins of equal integrated flux. To achieve the desired artistic effect, the image was inverted so that dark regions correspond to high signal.

The result, shown in 12, plots the final positions of the bin generators. The algorithm successfully handles this large and complex input, which features sharp, irregular discontinuities. As expected from the connection to optimal transport, the generators form a ‘blue noise’ point distribution, where their density traces the underlying signal structure. This example showcases the computational efficiency and robustness of PowerBin on large datasets. The entire process took just 20 seconds on a standard laptop, with 12 seconds for the bin-accretion phase and 8 seconds for the CPD regularization. This highlights its potential as a general-purpose tool for capacity-constrained tessellation in a wide variety of scientific and technical fields.

8 Execution Time Benchmarks↩︎

To quantify the performance improvement of the new algorithm, I conducted a series of benchmark tests comparing the execution time of PowerBin against the classic VorBin package. All tests were performed on a standard laptop with an Intel i7-1355U processor. The process was running on a single core at a sustained frequency of about 3GHz. The results are presented in 13.

For these tests, I generated a sequence of mock galaxy images of increasing size. The input signal for all tests was a simulated galaxy with an exponential surface brightness profile (Sérsic \(n_{\rm Ser}=1\)), a fixed axial ratio of 4/3, and Poisson noise, corresponding to the uncorrelated noise case shown in the top panel of 9. I created five images, starting at \(320\times240\) pixels and progressively doubling the total number of pixels up to \(1280\times960\). For each image, the target signal-to-noise was chosen to make the number of bins, \(n\), scale proportionally with the total number of pixels, \(N\). This approach, which keeps the average number of pixels per bin constant, simulates a common scientific goal: exploiting a larger number of pixels to increase the spatial sampling of the binned map (i.e., more bins). This setup provides a realistic benchmark of how the algorithm’s performance scales as both the number of input pixels and output bins increase. I adopted a number of bins logarithmically spaced from \(n=1600\) to \(n=25600\).

13 plots the execution time versus the number of input pixels on a log-log scale for the two main computational stages: the accretion and regularization steps for both the old VorBin and the new PowerBin algorithms. The performance difference is dramatic and confirms the expected theoretical scaling laws.

The classic VorBin method, shown by the upper two curves, exhibits a computational time that scales significantly more steeply than \(\mathcal{O}(N\log N)\). At large \(N\), both its accretion (blue triangles) and regularization (orange diamonds) stages approach a scaling consistent with \(\mathcal{O}(N^2)\) (dot-dashed brown line). This is the expected behaviour, as the regularization is based on a MWV diagram, which has a quadratic time complexity. Similarly, the classic bin-accretion algorithm also performs operations on average linear in \(N\) for every pixel, leading to a \(\mathcal{O}(N^2)\) time complexity.

In stark contrast, the new PowerBin algorithm, shown by the lower two curves, demonstrates vastly superior performance. Both the new fast bin-accretion stage (green circles) and the Centroidal Power Diagram regularization (red squares) follow a trend that is nearly perfectly described by the theoretical \(\mathcal{O}(N\log N)\) scaling (dashed purple line). This is the optimal complexity for this class of geometric problem and is a direct result of the algorithmic improvements described in 5 and 6. The bottom line is that for a dataset with one million pixels, the new algorithm is approximately two orders of magnitude faster than the previous standard, turning a computation that would take 6 hours into one that takes 3 minutes. This efficiency gain is critical for the practical analysis of large-scale astronomical surveys.

It is important to emphasize that this benchmark compares the relative performance and algorithmic scaling of the two methods. Both the classic VorBin and the new PowerBin are implemented entirely in Python, relying on standard scientific libraries like NumPy [58] and SciPy [43]. The absolute execution times could be substantially reduced by porting the computationally-intensive parts to a compiled language or specialized hardware like GPUs. However, such optimizations would not change the fundamental time complexity of the algorithms. The key result of this comparison is the difference in scaling—\(\mathcal{O}(N\log N)\) versus \(\mathcal{O}(N^2)\)—which demonstrates the inherent efficiency of the new approach, independent of the specific implementation.

9 Conclusions↩︎

In this paper, I have introduced PowerBin, a new algorithm for the adaptive binning of two-dimensional data. This work was motivated by the increasing scale of modern astronomical surveys and the limitations of existing methods, which are either too slow or lack guarantees of bin convexity. The main contributions of this work can be summarized as follows:

1. A New Theoretical Framework: I have framed the adaptive binning problem within the mathematical theory of optimal transport. The natural solution in this framework is a Centroidal Power Diagram (CPD), a generalization of a Centroidal Voronoi Tessellation that rigorously accommodates capacity constraints while guaranteeing convex bins.

2. A Fast and Robust Heuristic Solver: Formal CPD solvers, based on gradient-based optimization of a dual energy functional, are ill-suited to real astronomical data, which is discrete and often has non-additive noise properties. I have introduced a novel heuristic algorithm that circumvents these issues. It is based on a simple physical insight into the geometry of packed cells, which provides a direct, non-linear update rule for the power diagram weights to enforce the target capacity. This approach is fast, robust, and converges reliably even when the capacity function is non-additive, a regime where formal methods fail.

3. An Optimized Bin-Accretion Algorithm: The bin-accretion stage, which provides the crucial starting point for the iterative refinement, would have become the computational bottleneck of the current methods. Therefore, I developed a new implementation with near-linear time complexity, \(\mathcal{O}(N\log N)\). This is achieved by using a pre-computed Delaunay triangulation for adjacency information, a frontier-based growth strategy, incremental updates for bin properties, and a heap-managed frontier to efficiently select pixels.

4. Superior Performance and Scalability: The combination of the fast CPD solver and the optimized bin-accretion algorithm results in a dramatic performance improvement. Benchmark tests show that the entire PowerBin algorithm scales as \(\mathcal{O}(N\log N)\), in stark contrast to the \(\mathcal{O}(N^2)\) scaling of previous methods. For a dataset of one million pixels, the new algorithm is approximately two orders of magnitude faster than the widely-used VorBin package.

I have demonstrated through a series of tests on simulated and real data that PowerBin produces high-quality, convex tessellations with excellent capacity uniformity. It successfully handles a wide range of signal distributions and is robust to the challenges of correlated noise. Its robustness was further demonstrated on challenging background-subtracted mock observations containing multiple objects, where the noise is highly non-Poissonian. Its performance on the classic SAURON data of NGC 2273 confirms that it reproduces and improves upon the results of the original Voronoi-binning method.

By addressing the key limitations of speed and convexity, PowerBin provides a powerful and scalable tool for the analysis of the massive datasets generated by current and future astronomical surveys. Its versatility, demonstrated on a non-astronomical imaging problem, also suggests its potential for broad application in other scientific and technical fields requiring capacity-constrained tessellation. Notably, both the bin-accretion and regularization stages of PowerBin extend naturally to higher dimensions with minimal code changes and no change in the underlying concept. The Python implementation of the algorithm is publicly available to the community in the PowerBin package⁵.

Acknowledgements↩︎

I thank the referee for a constructive and insightful report that helped improve the clarity and completeness of this paper. I am grateful for the serene setting of the island of Pellestrina, whose tranquil location between the Venetian Lagoon and the Adriatic Sea proved to be the perfect environment for the writing of this paper. A special thanks to my family, who demonstrated saintly patience while enduring countless updates about the PowerBin algorithm.

Data Availability↩︎

No new data were generated in support of this research. The PowerBin package is available at https://pypi.org/project/powerbin/.

References↩︎