Centroidal Power Diagram (CPD) Overview

• Centroidal Power Diagrams (CPDs) are geometric partitions where each generator is the mass centroid of its power cell defined by weighted points.
• They extend classical Voronoi tessellations by incorporating weight adjustments via a mass-cost relation, enabling optimal transport and adaptive binning.
• CPDs are computed efficiently using generalized Lloyd algorithms and Hamilton–Jacobi formulations, with applications in material science, data compression, and astronomy.

A centroidal power diagram (CPD) is a power diagram in which each generator coincides with the centroid of its power cell. In the Euclidean formulation, the diagram is defined by weighted points $(\mathbf{x}_i,w_i)$, with cells

$$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$

and the centroidal condition requires

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$

Within the framework developed for optimal location and quantization, critical points of the relevant energy are centroidal power diagrams, and centroidal Voronoi tessellations arise as the special case of constant weights (Bourne et al., 2014). The concept extends beyond Euclidean geometry: geodesic centroidal power diagrams replace Euclidean distance by a geodesic distance $d_C$, admit capacity constraints, and can be characterized by a system of Hamilton–Jacobi equations (Camilli et al., 2021). More recent work uses CPDs as the geometric solution class for capacity-constrained adaptive binning, emphasizing convex cells and computational scalability (Cappellari, 8 Sep 2025).

1. Definition and geometric structure

Power diagrams generalize Voronoi diagrams by assigning a real weight to each generator. In the Euclidean setting, a power cell is determined by the comparison of the quantities $|\mathbf{x}-\mathbf{x}_i|^2-w_i$, so the tessellation depends on both locations and weights (Bourne et al., 2014). In the formulation used for adaptive binning, the same construction is written as

$$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$

and its boundaries are hyperplanes; thus, cells are convex polytopes (Cappellari, 8 Sep 2025).

A centroidal power diagram is obtained when the generator of each power cell is also its centroid. In the Euclidean mass-weighted formulation, this is expressed by

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$

with the associated mass $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$ (Bourne et al., 2014). The same paper states that, at a critical point of the energy, the weights satisfy

$w_i=-f'(m_i),$

up to adding the same constant to all $w_i$, which does not change the power diagram (Bourne et al., 2014).

The relation to centroidal Voronoi tessellations is direct. Centroidal Voronoi tessellations are tessellations by Voronoi cells such that each generator coincides with the centroid of its cell; CPDs generalize them by replacing Voronoi cells with power cells and by adapting the weights through the mass–cost relation $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$0 (Bourne et al., 2014). In the geodesic setting, geodesic centroidal Voronoi tessellations are recovered when all weights are equal, so geodesic centroidal power diagrams generalize geodesic centroidal Voronoi tessellations in the same sense (Camilli et al., 2021).

2. Variational and optimal-transport formulations

The principal variational formulation in the Euclidean theory considers the energy

$$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$1

where the second term is the squared $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$2-Wasserstein distance between the density $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$3 and the atomic measure $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$4 (Bourne et al., 2014). The same work states that the target energy can be written either in terms of atomic measures and the Wasserstein distance or in terms of weighted points and power diagrams, and that the latter formulation is more suitable for computation (Bourne et al., 2014).

Under suitable conditions, the Wasserstein term admits the geometric representation

$$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$5

where the partition $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$6 is a power diagram such that $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$7 (Bourne et al., 2014). This leads to the reformulated energy

$$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$8

with each $$P_i=\left\{\mathbf{x}\in\Omega:\ |\mathbf{x}-\mathbf{x}_i|^2-w_i\leq |\mathbf{x}-\mathbf{x}_j|^2-w_j,\ \forall j\right\},$$9 defined by the power-diagram inequalities (Bourne et al., 2014).

The optimal-transport interpretation is made explicit again in the adaptive-binning context. There the problem is posed as a semi-discrete optimal transport or quantization problem: partition the domain into bins $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$0 minimizing

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$1

subject to

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$2

where $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$3 is the desired capacity for bin $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$4 (Cappellari, 8 Sep 2025). The solution is a power diagram with properly chosen weights, adjusted until each cell’s capacity equals its target (Cappellari, 8 Sep 2025). A dual functional is given by

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$5

with gradients

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$6

so that a CPD is achieved when capacity is matched for all bins and generators coincide with the centroids of their cells (Cappellari, 8 Sep 2025).

This suggests a unifying interpretation: in both the optimal-location formulation and the adaptive-binning formulation, the CPD is the geometric structure induced by an energy or transport problem in which both assignment and mass distribution are optimized.

3. Critical-point conditions and generalized Lloyd iteration

The Euclidean theory identifies CPDs as critical points of the energy $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$7. The paper states that, at a critical point, the “CPD conditions” are

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$8

again up to adding the same constant to all $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\,\rho(\mathbf{x})\,d\mathbf{x}, \qquad m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}.$$9 (Bourne et al., 2014). The first condition is the centroid condition; the second links the weight to the cell mass through the derivative of the cost function.

These conditions motivate a generalization of Lloyd’s algorithm. Given current generators $d_C$0, the procedure is: compute the power diagram $d_C$1, compute the mass $d_C$2, compute the centroid

$d_C$3

and then update

$d_C$4

optionally removing empty cells or generators (Bourne et al., 2014). In the case $d_C$5 and fixed weights, this reduces to classical Lloyd’s algorithm for centroidal Voronoi tessellations (Bourne et al., 2014).

The convergence properties reported for this generalized Lloyd algorithm are explicit. The algorithm is energy decreasing,

$d_C$6

with equality if and only if a fixed point, hence a critical point, is reached (Bourne et al., 2014). Under mild assumptions—fixed $d_C$7, no cell elimination in the limit, finitely many critical points per energy level—convergence to a critical point is proven (Bourne et al., 2014). Numerical results indicate linear convergence (Bourne et al., 2014).

The same work also records a structural limitation: the algorithm can remove but not create new generators, so overestimating $d_C$8 in initialization is recommended (Bourne et al., 2014). That property is often important in practice because the effective number of nonempty cells may decrease during optimization.

4. Geodesic extension and Hamilton–Jacobi characterization

The geodesic theory replaces Euclidean distance by a geodesic distance $d_C$9 defined in terms of a family $|\mathbf{x}-\mathbf{x}_i|^2-w_i$0 of convex, compact sets in $|\mathbf{x}-\mathbf{x}_i|^2-w_i$1, with trajectories satisfying the inclusion $|\mathbf{x}-\mathbf{x}_i|^2-w_i$2 (Camilli et al., 2021). Within this framework, the geodesic power diagram associated with centers $|\mathbf{x}-\mathbf{x}_i|^2-w_i$3 and weights $|\mathbf{x}-\mathbf{x}_i|^2-w_i$4 is

$|\mathbf{x}-\mathbf{x}_i|^2-w_i$5

which reduces to the geodesic Voronoi diagram if all $|\mathbf{x}-\mathbf{x}_i|^2-w_i$6 are equal (Camilli et al., 2021).

A geodesic centroidal power diagram is then a power diagram, using geodesic distances and weights, where each generator coincides with the centroid of its cell with respect to $|\mathbf{x}-\mathbf{x}_i|^2-w_i$7: $|\mathbf{x}-\mathbf{x}_i|^2-w_i$8 The formulation also allows capacity constraints

$|\mathbf{x}-\mathbf{x}_i|^2-w_i$9

with prescribed capacities $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$0 (Camilli et al., 2021).

A central result is the characterization of geodesic centroidal tessellations by a system of Hamilton–Jacobi equations. For fixed $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$1, the geodesic distance $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$2 solves

$$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$3

For geodesic centroidal power diagrams, the paper gives the system

$$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$4

Its interpretation is also given explicitly: $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$5, the sets $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$6 define the power cells, the centroids $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$7 minimize the sum of distances over each cell, and the weights $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$8 are adjusted to meet the capacity constraints (Camilli et al., 2021).

This extension broadens CPDs from Euclidean quadratic-cost geometry to geodesic, possibly anisotropic or Riemannian metrics (Camilli et al., 2021). A plausible implication is that the centroidal-power-diagram paradigm is not tied to Euclidean convex polytopes alone, but to a more general coupling of assignment geometry, centroid conditions, and cell-capacity control.

5. Numerical methods and computational properties

In the Euclidean setting, power diagrams can be computed efficiently using lifting methods, with $$\mathcal{V}_{\rm pow}(\mathbf{g}_j)=\big\{\mathbf{x}\mid \|\mathbf{x}-\mathbf{g}_j\|^2-w_j\leq \|\mathbf{x}-\mathbf{g}_k\|^2-w_k,\ \forall k\neq j\big\},$$9 complexity in $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$0D (Bourne et al., 2014). The same source states that centroids and masses for each cell can be computed either analytically for constant $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$1 or numerically otherwise, and that the generalized Lloyd algorithm is practical in $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$2D and $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$3D (Bourne et al., 2014).

In the geodesic setting, numerical computation proceeds through a discretized Hamilton–Jacobi solve. The reported steps are: domain discretization by triangulating the domain into a mesh $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$4 with centroids $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$5; initialization with generators $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$6 and weights $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$7; solution of the discrete Hamilton–Jacobi equation

$$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$8

using a semi-Lagrangian scheme or a Fast Marching method adapted to unstructured grids; cell assignment by the smallest value of $$\mathbf{x}_i=\frac{1}{m_i}\int_{P_i}\mathbf{x}\rho(\mathbf{x})\,d\mathbf{x},$$9; centroid update; weight update for capacity-constrained power diagrams, for example by maximizing a concave Lagrangian functional; and iteration until convergence (Camilli et al., 2021). The Fast Marching method is described as updating each node once in the “correct” upwind order, making it highly efficient, especially on large or unstructured grids (Camilli et al., 2021).

The adaptive-binning literature introduces a distinct computational perspective. There, formal CPD solvers are described as unstable with real data, particularly because optimal-transport theory assumes continuous densities whereas real data are pixelated, and because non-additive capacities break the assumptions behind the analytic convex dual functional and gradient-based approaches (Cappellari, 8 Sep 2025). The paper reports the empirical result that rigorous gradient-based optimization as described in prior works works only for additive capacities and large bins, while for non-additive signals and small bins, “the method fails catastrophically” (Cappellari, 8 Sep 2025).

To address this, PowerBin introduces a heuristic based on a soap-bubble analogy. Each bin is represented by a disk of radius $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$0, with weight $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$1, and the area is approximated by $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$2 (Cappellari, 8 Sep 2025). The update rule is

$m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$3

or equivalently

$m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$4

followed by updating each generator to be the geometric centroid of its pixels (Cappellari, 8 Sep 2025). The iterative process alternates pixel assignment, bin-size update, and generator-position update until convergence (Cappellari, 8 Sep 2025).

The computational complexity emphasized in that work is $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$5: power-diagram assignment is $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$6 via KD-Tree and geometric lifting, generator and radius updates are $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$7 per iteration, and the combined bin-accretion and CPD-regularization stages scale as $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$8 (Cappellari, 8 Sep 2025). The same paper states that the resulting method is about two orders of magnitude faster than previous methods on million-pixel datasets and gives the benchmark statement that PowerBin is $m_i=\int_{P_i}\rho(\mathbf{x})\,d\mathbf{x}$9 faster than VorBin for million-pixel datasets (Cappellari, 8 Sep 2025).

6. Applications, comparisons, and recurring misconceptions

The range of applications associated with CPDs is broad in the supplied literature. In the Euclidean optimal-location formulation, applications include optimal location problems, quantization and data compression, approximation of a density by atomic measures in Wasserstein metric, and pattern formation in material science (Bourne et al., 2014). For block copolymers, the paper considers the energy

$w_i=-f'(m_i),$0

and reports that, for $w_i=-f'(m_i),$1, minimizers in $w_i=-f'(m_i),$2D give hexagonal or triangular patterns, while in $w_i=-f'(m_i),$3D a BCC lattice is numerically verified in the paper (Bourne et al., 2014).

The geodesic literature provides numerical experiments for Euclidean and non-Euclidean metrics, including Minkowski and Riemannian metrics, isotropic and anisotropic metrics, constant and variable metrics, regular domains such as the disk, square, and L-shape, as well as complex domains such as text shapes and a rabbit silhouette (Camilli et al., 2021). It also reports demonstrations for $w_i=-f'(m_i),$4 on unstructured grids, and observations that centroidal diagrams adapt well to non-uniform densities while capacity-constrained power diagrams produce cells of prescribed sizes, even in complicated domains (Camilli et al., 2021).

The adaptive-binning application focuses on astronomy and integral-field spectroscopy. There, CPDs are presented as the solution class for adaptive binning with prescribed capacity, such as target signal-to-noise ratio, while guaranteeing convex bins (Cappellari, 8 Sep 2025). The paper reports performance on mock galaxies, correlated-noise data, large mosaics with multiple sources and backgrounds, a classic real-world IFS dataset, and general image binning such as stippling for computer graphics (Cappellari, 8 Sep 2025).

Several misconceptions are addressed implicitly by the literature.

First, a CPD is not merely a Voronoi diagram with renamed parameters. The Euclidean definition uses power cells with variable weights, and the critical-point condition includes both centroid matching and the weight law $w_i=-f'(m_i),$5 (Bourne et al., 2014).

Second, the relation between CPDs and CVTs is one of strict specialization, not equivalence. CVTs are the special case of CPDs with constant weights in the Euclidean formulation, and geodesic centroidal Voronoi tessellations are the special case of geodesic centroidal power diagrams when all weights are equal [(Bourne et al., 2014); (Camilli et al., 2021)].

Third, formal CPD optimal-transport solvers are not universally robust on discrete real data. The adaptive-binning literature explicitly states that, for non-additive signals and small bins, formal optimization may fail catastrophically (Cappellari, 8 Sep 2025). This does not negate the theoretical role of CPDs; rather, it distinguishes the ideal geometric object from the numerical strategy used to approximate it in data analysis.

Aspect	CVT	CPD
Tessellation	Voronoi (equal weights)	Power diagram (variable weights)
Generator condition	Must be centroids of their cells	Must be centroids of their power cells
Weight condition	Fixed/ignored	Adapted via mass and cost $w_i=-f'(m_i),$6

Across these formulations, CPDs function as a common geometric language for partition problems in which location, mass, and transport cost are jointly constrained. The Euclidean theory emphasizes energy minimization and optimal location (Bourne et al., 2014), the geodesic theory emphasizes Hamilton–Jacobi characterization and non-Euclidean metrics (Camilli et al., 2021), and recent large-scale data analysis emphasizes capacity control, convexity, and $w_i=-f'(m_i),$7 computational scaling (Cappellari, 8 Sep 2025).