Reflectors, power cells and a sphere

Scientific reference

Pedro Machado Manhães de Castro, Quentin Mérigot and Boris Thibert. Far-field reflector problem and intersection of paraboloids. Numerische Mathematik 134(2), 389–411, 2016. DOI: 10.1007/s00211-015-0780-z.

Demo by Pedro M. M. de Castro. Derived from: P. M. M. de Castro, Q. Mérigot, and B. Thibert, “Far-field reflector problem and intersection of paraboloids,” Numerische Mathematik, 134(2), 389–411 (2016). https://doi.org/10.1007/s00211-015-0780-z.

We would like to thank the paper’s co-authors Quentin Mérigot and Boris Thibert. Tools for computational exploration and implementation: OpenAI Codex and Anthropic Claude. These acknowledgements do not imply implementation participation or endorsement by the co-authors.

Equation, proposition, theorem, section and figure numbers below are those of this journal version. An earlier preprint by the same authors (arXiv:1403.0062v1, 17 pages, titled Intersection of paraboloids and application to Minkowski-type problems) and the SoCG 2014 conference paper share the transformation and the structure of the algorithm with different numbering; they are supporting context only.

From reflectors to power cells

A solid paraboloid of revolution with focus at the origin, unit axis y and parameter λ > 0 is P(y, λ) = {x : ‖x‖ − x·y ≤ λ}; λ/2 is the usual focus–vertex distance and the vertex is −(λ/2)y. Its boundary is the radial graph of r(u) = λ/(1 − y·u) over the unit sphere (Eq. 2.1). For a family (y_i, λ_i), the reflector is the boundary of K = ∩_i P(y_i, λ_i), the radial graph of r(u) = min_i r_i(u) (Eq. 2.2), and the cell of i is the set of directions where i attains this minimum (Definition 2.1). Ties belong to every cell that attains the minimum. A ray leaving the focus along u in the cell of i is reflected into the direction y_i: with n = (u − y_i)/‖u − y_i‖, u − 2(u·n)n = y_i. The direction u = y_i itself is singular (r_i = +∞); with at least two distinct directions the solid K is bounded.

Write q_i = 1/λ_i. Proposition 1 uses the additive power D_i(x) = ‖x − p_i‖² + ω_i with

On the sphere this gives D_i(u) = 1 − a_i(u) with a_i(u) = q_i(1 − y_i·u) = 1/r_i(u): the minimal radius, the maximal a_i and the minimal additive power select the same indices, ties included, and r(u) = 1/max_i a_i(u). The site p_i points away from the outgoing direction y_i. CGAL’s Regular_triangulation_3 subtracts the weight, ‖x − p‖² − w, so the same diagram is obtained with w = −ω_i = q_i + q_i²/4; the demo labels both conventions wherever weights appear.

The §4 algorithm as implemented

The demo follows Section 4 (journal pp. 401–406). Its output is the oriented boundary of every intersection cell S² ∩ Pow(p_i): vertices on the sphere, oriented arcs and oriented cycles, with the cell on the right of every arc.

1. Regular triangulation. The weighted points are lifted to (p_i, ‖p_i‖² + ω_i) and the projection of the lower convex hull is built incrementally (Bowyer–Watson) with an infinite vertex, as in CGAL. If all sites are coplanar or collinear, the triangulation is built in their affine hull: its ridges are then lines (coplanar sites) or absent (collinear sites, whose facets are whole planes). A site whose lifted point lies above the lower hull is hidden: its power cell is empty. A site whose lifted point lies exactly on the lower hull without being a vertex has a power cell of zero volume, reported with the face on which it ties; identical weighted points are reported as duplicates.
2. Facets. Every edge (i, j) of the triangulation is dual to a facet of Pow(p_i) in the radical plane n·x = b, n = 2(p_j − p_i), b = ‖p_j‖² + ω_j − ‖p_i‖² − ω_i; the cell is on the side n·x ≤ b. plane_crossing_sphere compares b² with ‖n‖². A bounded facet whose vertices are all on or inside the sphere carries no arc (Fig. 6a); a vertex exactly on S² is recorded as a point contact.
3. Ridges. The tetrahedra around the edge are visited in the cyclic order given by their positive orientation; consecutive tetrahedra share a triangle whose dual ridge is a segment between power centres, a ray towards the outward hull normal, or a line. Infinite tetrahedra are handled as in CGAL; no facet is decided from its finite vertices alone.
4. Sources and targets. number_of_intersections counts the roots of the ridge’s supporting line with the sphere of radius 1 + ε (ε → 0) inside its parameter range: an endpoint (power centre) lying exactly on S² counts as inside the ball, so a root there is counted only if the ridge leaves the ball at that point, and a tangent ridge counts twice; walking the facet, an exit from the ball is a source and an entry is a target. Matching each source with the next target in the same cyclic sequence gives the arcs (Fig. 6b). If no ridge meets the sphere, the circle is enclosed exactly when its centre, the projection of the origin, lies in the facet (Fig. 6c).
5. Cycles. A vertex on the sphere is identified by its ridge (a triangle of the triangulation) and its root, never by rounded coordinates. It is the target of one arc of the cell and the source of the next, so the arcs close into cycles: the oriented connected components of the boundary graph.

Each arc of cell i on facet (i, j) is computed independently from both sides; the reversed twin in cell j is used as a consistency check. Two cells may share several disjoint arcs of the same circle (Fig. 5, right); they are kept as separate arcs.

What is exact

Inputs are rationals: decimals such as 0.1 are read as the decimal they denote. Directions must be exact unit vectors; a picked or typed vector that is not one is replaced, visibly, by a nearby rational point (2s, 2t, ±(1 − s² − t²))/(1 + s² + t²), and the inspector reports the angular offset. Every predicate — orientation, power test, plane and ridge against the sphere, membership of a circle centre in a facet, the order of the roots along a ridge, equality of two vertex positions, ownership at a vertex and whether an arc turns by more than π — is an exact sign of a rational expression or of a + b√r or a + b√r₁ + c√r₂ + d√(r₁r₂). For the costliest of these signs (a power centre against the sphere, ownership at a vertex, equality of two vertex positions and the turn of an arc) a filter answers first, in the spirit of the dynamic filtering of the journal implementation: the expression is evaluated in BigInt fixed point (scale 2⁻¹⁶⁰) with a proven bound on its error, and its sign is accepted only when the value exceeds that bound; otherwise it is computed with exact rationals. A filter never changes a decision, only its cost; all other signs are computed with exact rationals directly. Candidate pairs of coincident vertices are proposed by a sweep within 10⁻⁶ on coordinates obtained from the same fixed-point evaluation (accurate to about 10⁻¹⁶), so no exact coincidence is missed, and each candidate is then decided exactly. No input value is perturbed, rounded or deduplicated by tolerance.

Contact conventions. As in §4.1, a point on the sphere counts as inside, and a ridge tangent to the sphere has two coincident intersections (the limit of a sphere of radius 1 + ε). As in §4.3 (p. 405), a radical plane tangent to the sphere creates no arc; its tangency point, when it lies in the facet, is recorded as a point contact of the two cells. Only when a ridge passes through that tangency point do the two coincident roots of the ridge produce arcs, all of zero length; they are kept so that every cycle stays closed, and flagged. A power vertex exactly on S² of a facet inside the ball, and a point that a cell owns without any arc of its own reaching it (found from the exact owners of every vertex point), are also recorded as point contacts. A cell whose only contacts with the sphere are such points has zero area. Zero-size features are kept and flagged rather than merged; the equality of two vertex positions is decided exactly through the canonical supporting lines of their ridges.

The certificate pane re-checks the triangulation exhaustively (adjacency, combinatorial orientation, local regularity of every facet, convexity and lower-dimensional regularity of the hull). The test suite compares the output with an independent brute-force oracle written in Python; this is evidence, not a formal proof of the program.

What is approximate

• binary64 Coordinates of vertices and arc angles (square roots), and the polylines drawn on the sphere.
• binary64 Cell areas and masses for the uniform density 1/4π. With a pole S chosen off every boundary circle, α = (x dy − y dx)/(1 + z) in a frame whose south pole is S satisfies dα = area form, so Area(C) = −Σ_arcs ∫ α + 4π [S ∈ C] for cycles with the cell on the right; membership of S is decided exactly. The integrals use adaptive 8-point Gauss–Legendre quadrature: panels of at most π/8, each bisected until it agrees with its two halves within its share of an absolute tolerance of 10⁻¹⁴ (the share is halved at each level), at most 14 levels deep and never below a rounding floor of 8·10⁻¹⁶ times the magnitude of the panel; the result is compared with a run at tolerance 10⁻⁹ and the difference is shown. Small circles are handled without cancellation: the centre and the squared radius 1 − b²/‖n‖² are exact rationals converted once.
• display Figures (WebGL, three.js): on the sphere, the owner of each pixel is the minimal power evaluated on the graphics card in single precision (float32), and the certified arcs are drawn over it. The reflector is a mesh of about 39,000 triangles of directions moved radially to r(u) = 1/max_i a_i(u); the owner of each pixel is recomputed from its direction, the surface is transparent and drawn back faces first (K is convex), and the selected paraboloid is drawn whole up to a stated distance from the focus, because a paraboloid is infinite. The Section mode cuts every surface by the plane through the focus that contains u and the reflected direction. A click on a figure is resolved by binary64 ray casting. Without WebGL 2 the page uses the Canvas 2D fallback renderer (per-pixel binary64 evaluation and ray casting).

An exact combinatorial structure does not make these numbers exact; the quadrature, the drawing and any optimisation built on them are numerical.

Scaling λ

Under the coupled change λ_i → kλ_i (k > 0), q_i → q_i/k, so D_i(x) − D_j(x) = q_i(x·y_i − 1) − q_j(x·y_j − 1) is multiplied by 1/k at every point x ∈ ℝ³. Hence the full power cells in ℝ³ and their traces on the sphere are unchanged, while the representation changes: the sites move to p_i/k, the weights become −q_i/k − q_i²/(4k²), and every radius scales, r(u) → k r(u), so the reflector is scaled by k about the focus. The article remarks (p. 394) that the intersection of the power cells with the sphere does not change under such a scaling; with the parametrisation of Proposition 1 the cells in ℝ³ are invariant as well, and what changes is the weighted-point data. The demo’s scaling control recomputes the diagram and compares the certified boundaries before and after; the tests check both statements exactly.

Power diagrams and the inverse map

The “Power diagram” representation accepts arbitrary weighted points. Such a diagram can have cells with holes (Fig. 5, middle: a cube-shaped cell whose trace has six circular holes) and disconnected cells (a rod-shaped cell that crosses the sphere twice), unlike a paraboloid diagram, whose non-empty cells are connected (Prop. 3). Conversion back to paraboloids requires the canonical form: p_i ≠ 0, y_i = −p_i/‖p_i‖, λ_i = 1/(2‖p_i‖) and ω_i + ‖p_i‖² + 2‖p_i‖ = c for one common constant c (adding c to every weight does not change the cells). The test is exact even when ‖p_i‖ is irrational. It certifies the canonical representation only: it does not decide whether some other choice of weighted points with the same cells could be realised by paraboloids, and Lemma 9 (which concerns ellipsoids) is not an inverse for paraboloids. Incompatible input is rejected with the failing condition; no weight is dropped or rewritten.

Direct and inverse problems

The demo solves the direct problem: given λ it computes the cells. The far-field reflector problem asks for λ such that the light of a density ρ received by each cell matches prescribed amounts α_i (Section 5.1). Theorem 12 shows that γ_i = log λ_i solves it exactly when γ maximises the concave function Φ(γ) = ∫ min_i[c(u, y_i) + γ_i] ρ(u) du − Σ γ_i α_i with c(u, y) = −log(1 − u·y), whose gradient is M_i − α_i (Eq. 5.2). The authors maximise Φ with a quasi-Newton method, recomputing the cells and integrating with Gaussian quadrature (§5.2). The optional solver in this demo is described in its own panel, with its update rule, gauge, quadrature and stopping rule; it is not the authors’ solver.

Complexity and limits

For confocal paraboloids the diagram on the sphere has O(N) vertices, edges and faces (Theorem 2), while power diagrams and ellipsoid diagrams can have Ω(N²) edges on the sphere (Theorem 8; the “Quadratic complexity” data set follows the construction of its proof). Theorem 10 gives an O(N log N + C) algorithm, where C is the complexity of the 3D power diagram; whether C is O(N) for the power diagrams of Proposition 1 is Open Question 11. This demo uses its own incremental construction in input order with a stochastic walk and exact BigInt arithmetic, and makes no complexity claim. Limits: at most 64 paraboloids or weighted points; each input value at most 96 bits of numerator plus denominator; a bounded number of triangulation steps; intermediate exact values of at most 6000 bits; and 20 s of construction, after which it stops with a message and the previous diagram stays on screen. The construction yields to the page between cells and can be cancelled. The data sets build in tens to hundreds of milliseconds on a laptop; 64 sites whose coordinates have unrelated 30- to 47-bit denominators take about 2–4 s.

Data sets

The ten reference data sets provide fixed examples with expected geometric results. Other data sets reproduce configurations of the article (Fig. 1, Fig. 5, Prop. 4, the proof of Theorem 8) or isolate degeneracies: four cells meeting at both poles, a power centre exactly on the sphere, tangent planes and ridges, identical and dominated paraboloids, parameters across twelve orders of magnitude and directions about 2·10⁻⁶ rad apart. Spread directions are exact rational points obtained from a Fibonacci sphere by the inverse stereographic map.

Typography, assets and licences

Latin Modern is used for headings, explanations and mathematical symbols, Inter for controls and labels, and Latin Modern Mono for exact values. The fonts and the CIn/UFPE mark are local copies of the public assets of the academic site; the GUST Font License and the Inter licence are included. Interface pictograms are local SVG shapes.

The 3D figures use three.js r186 (npm package three@0.186.1, MIT licence, licence text), vendored locally in vendor/three/; only its fat-line add-ons were changed, to import three.js by a relative path. The geometric core does not depend on it. Details: licenses/README.md.