Wavefront Fitting over Arbitrary Freeform Apertures via CSF-Guided Progressive Quasi-Conformal Mapping

1. Introduction

Advanced optical manufacturing increasingly demands nanometer surface accuracy for freeform optical elements [1,2]. Wavefront testing is therefore a critical component of the metrology pipeline, and reconstruction accuracy directly affects the reliability of closed-loop fabrication. Conventional wavefront analysis commonly relies on Zernike polynomials defined on the unit disk, where orthogonality and physical interpretability are advantageous [3,4,5]. However, for arbitrary freeform apertures (e.g., non-convex regions, high-aspect-ratio pupils, or apertures with obscurations), a mismatch between the physical support and the unit disk induces significant modal crosstalk and coefficient coupling [6,7,8].

Two broad strategies have been used to mitigate these limitations.

First, customized orthogonal bases can be constructed by Gram–Schmidt orthogonalization or singular value decomposition (SVD) on sampled apertures[9,10]. Recent developments in this direction have further optimized orthogonal fitting algorithms for aberration removal on arbitrary shaped apertures[11] and data-driven approaches utilizing deep neural networks [12] have emerged. While effective, numerical bases are aperture-specific and complicate standardization, whereas deep learning models often lack physical interpretability and require extensive training datasets.

Second, conformal maps, such as the Fornberg algorithm [13], variants for slender regions [14], or the Schwarz–Christoffel transformation [15], can map a non-circular region to a canonical domain before Zernike fitting. Notably, the Schwarz–Christoffel mapping has been recently adapted for modal wavefront reconstruction on noncircular pupils[16], demonstrating the enduring relevance of mapping-based approaches. However, standard conformal maps may be numerically fragile on non-smooth boundaries due to derivative singularities (the crowding phenomenon)[17], producing severe grid distortion.

Although advances in computational geometry, such as Ricci flow [18,19] and optimal mass transport [20], have improved mesh quality for complex topologies, for doubly connected apertures (e.g., annular pupils), mapping must also respect topological invariants; particularly, preserving the conformal modulus is necessary to avoid excessive shear.

Quasi-conformal mapping (QCM) provides additional flexibility by allowing for bounded angular distortion, which is quantitatively characterized by the Beltrami coefficient [21]. When the distortion magnitude is bounded (${\parallel \mu \parallel }_{\infty }<1$), the mapping remains a homeomorphism; the conformal case is recovered when $\mu \equiv 0$. However, leveraging this flexibility to automatically construct a low-distortion, topology-consistent, and boundary-robust quasi-conformal parameterization for arbitrary freeform apertures remains a significant challenge.

This paper introduces a CSF-guided progressive quasi-conformal mapping framework (CSF-QCM). By coupling boundary smoothing through geometric evolution with topology-aware canonicalization and interior harmonic relaxation, the method yields a unified low-distortion parameterization for simply and doubly connected apertures. The resulting diffeomorphic map improves sampling uniformity on the canonical domain, reduces design-matrix ill-conditioning, suppresses boundary ringing artifacts, and improves numerical stability in wavefront fitting.

2. Methodology

2.1. Notation and Problem Setup

Let ${\Omega }_{\mathrm{phy}}\subset {\mathbb{R}}^2$ denote the physical aperture domain, sampled on a triangular mesh with vertices $\mathcal{V}$ and faces $\mathcal{F}$ using standard discrete geometric processing techniques. The boundary may be simply connected (only an outer boundary $\partial {\Omega }_{\mathrm{out}}$) or doubly connected (an outer boundary $\partial {\Omega }_{\mathrm{out}}$ and an inner boundary $\partial {\Omega }_{\mathrm{in}}$ representing an obscuration).

We seek a diffeomorphic parameterization $\Phi :{\Omega }_{\mathrm{phy}}\to {\Omega }_{\mathrm{can}}$, where ${\Omega }_{\mathrm{can}}$ is a canonical domain:

• Simply connected: ${\Omega }_{\mathrm{can}}=\mathbb{D}=\{w\in \mathbb{C}:|w|\le 1\}$ (unit disk).
• Doubly connected: ${\Omega }_{\mathrm{can}}=\mathbb{A}=\{w\in \mathbb{C}:R_{\mathrm{in}}\le \left|w\right|\le 1\}$ (concentric annulus).

Wavefront measurements are given as samples $W(x_i,y_i)$ on ${\Omega }_{\mathrm{phy}}$. After mapping to ${\Omega }_{\mathrm{can}}$, we fit W using an orthogonal basis on ${\Omega }_{\mathrm{can}}$.

2.2. Framework overview

The CSF-QCM pipeline (Figure 1) consists of four distinct steps:

• Boundary geometric regularization (CSF).We evolve the physical boundary $\partial \mathcal{M}$ under discrete Curve Shortening Flow. This iterative process smooths high-curvature features (e.g., sharp corners of a star-shaped aperture) and generates a sequence of regularized boundary frames that converge to a canonical circle, minimizing grid crowding effects.
• Topology-aware canonical domain specification. The target domain ${\Omega }_{\mathrm{can}}$ is determined by the aperture’s topology. For simply connected domains, ${\Omega }_{\mathrm{can}}$ is set to the unit disk $\mathbb{D}$. For doubly connected domains (e.g., apertures with central obscurations like a triangular hole), we solve the harmonic equation to compute the Dirichlet energy and derive the conformal modulus. This uniquely defines the inner radius $R_{\mathrm{in}}$ of the target annulus $\mathbb{A}$, ensuring a conformal bijection.
• Interior mesh optimization via Laplacian construction. We construct the discrete Laplacian matrix $\mathbf{L}$ using cotangent weights to approximate the harmonic energy. The interior parameterization is obtained by solving the sparse linear system $\mathbf{L}\mathbf{u}=-{\mathbf{L}}_{\mathrm{bd}}{\mathbf{u}}_{\mathrm{bd}}$, where the boundary conditions ${\mathbf{u}}_{\mathrm{bd}}$ are updated progressively using the CSF frames. This step relaxes the internal mesh vertices to minimize angular distortion.
• Wavefront resampling and orthogonal fitting. Using the computed bijection $\Phi$, we map the wavefront data into the canonical domain. Depending on the topology determined in step 2, the wavefront is fitted using standard Zernike polynomials (for disk $\mathbb{D}$) or Annular Zernike polynomials (for annulus $\mathbb{A}$), allowing for high-precision reconstruction over arbitrary free-form apertures.

2.3. Boundary Regularization via Curve Shortening Flow

Conformal maps can be numerically fragile on non-smooth boundaries because corners and sharp curvature variations induce large derivatives and crowding. We therefore smooth the boundary geometry prior to parameterization.

Let $\mathcal{C}(u,t):S^1\times [0,T]\to {\mathbb{R}}^2$ denote a parameterized boundary curve, where t is evolution time and $\mathbf{N}(u,t)$ is the inward unit normal. Under curve shortening flow (CSF), the boundary evolves by[22]

$${\displaystyle \frac{\partial \mathcal{C}(u,t)}{\partial t}}=\kappa (u,t)\mathbf{N}(u,t),$$

(1)

where $\kappa (u,t)$ is curvature. Intuitively, high-curvature regions move faster, smoothing the curve.

Stopping criterion.

Because CSF ultimately shrinks embedded plane curves to a round point, the evolution must be stopped early. We monitor a circularity deviation metric $E_{\mathrm{circ}}^{\left(k\right)}$ (defined on the discrete boundary at iteration k) and terminate when the relative change is sufficiently small:

$$\left|{\displaystyle \frac{E_{\mathrm{circ}}^{(k+1)}-E_{\mathrm{circ}}^{\left(k\right)}}{E_{\mathrm{circ}}^{\left(k\right)}+\epsilon }}\right|<{ϵ}_{\mathrm{tol}},$$

(2)

where ${ϵ}_{\mathrm{tol}}={10}^{-6}$ and $\epsilon$ is a small constant to avoid division by zero. This criterion indicates that high-frequency boundary irregularities have been sufficiently attenuated while preserving topology.

For doubly connected apertures, the same evolution is applied to both $\partial {\Omega }_{\mathrm{out}}$ and $\partial {\Omega }_{\mathrm{in}}$ (with appropriate inward normals defined with respect to the domain).

2.4. Progressive Quasi-Conformal Mapping and Distortion

Let $z=x+\mathrm{i}y$ denote the complex coordinate in the physical plane and $w=\Phi \left(z\right)$ denote the canonical coordinate. The Beltrami coefficient of $\Phi$ is defined by[23]

$${\mu }_{\Phi }\left(z\right)={\displaystyle \frac{\partial \Phi /\partial \overline{z}}{\partial \Phi /\partial z}},{\parallel {\mu }_{\Phi }\parallel }_{\infty }<1,$$

(3)

where ${\mu }_{\Phi }\equiv 0$ corresponds to conformality.

Progressive strategy.

To avoid instability under large deformations, we decompose the mapping into a sequence of small transitions. The intermediate boundary states are provided by the CSF frames. At each frame t, we impose Dirichlet boundary positions on $\partial {\Omega }_{\mathrm{can}}$ and compute interior vertex positions by harmonic relaxation. To ensure mesh robustness and prevent triangle flips, we employ the graph Laplacian[24], denoted by $\mathbf{L}$,which is used here purely as a numerical smoothing operator rather than a geometric Laplace–Beltrami discretization.

The interior coordinates are solved via the linear system:

$${\mathbf{L}}_{II}{\mathbf{P}}_I^{\left(t\right)}=-{\mathbf{L}}_{IB}{\mathbf{P}}_B^{\left(t\right)},$$

(4)

where ${\mathbf{P}}_B^{\left(t\right)}$ stores the fixed boundary vertex coordinates on $\partial {\Omega }_{\mathrm{can}}$, and ${\mathbf{P}}_I^{\left(t\right)}$ stores the unknown interior coordinates. The matrices ${\mathbf{L}}_{II}$ and ${\mathbf{L}}_{IB}$ are the sub-blocks of $\mathbf{L}$ corresponding to the interior (I) and boundary (B) vertex indices, respectively.

The rationale for employing the Laplacian operator lies in its connection to Dirichlet energy minimization. The solution to the harmonic equation $\Delta u=0$ corresponds to the configuration that minimizes the stretching energy of the mapping, analogous to the equilibrium state of a stretched rubber sheet or a spring network.

In our framework, once the boundary geometric tension is released by the CSF evolution, the Laplacian solver naturally propagates this relaxation into the domain interior. By minimizing the local Dirichlet energy, the Laplacian operator promotes smoothness and uniformity in the coordinate field. This intrinsic smoothing property prevents the formation of singular gradients, thereby ensuring that the resulting Beltrami coefficient ${\mu }_{\Phi }$—which depends on the coordinate derivatives—remains low and spatially continuous.

Boundary correspondence.

Rather than prescribing point-to-point matching manually, boundary correspondence is induced by CSF vertex trajectories with arc-length redistribution. This allocates parameter-space resolution to regions that are geometrically difficult (high curvature or near concavities), mitigating the crowding problem typical of static conformal maps.

Area uniformity and monitoring.

We monitor distortion using both $|{\mu }_{\Phi }|$ (angular distortion) and an area-scaling statistic derived from the Jacobian determinant $J_{\Phi }\left(z\right)=|detD\Phi \left(z\right)|$. In discrete form, $J_{\Phi }$ is evaluated per triangle as the area ratio between mapped and physical triangles. Ideally, the mapping should preserve the area element to maintain the orthogonality of the polynomial basis defined on the physical domain. Significant variations in the Jacobian determinant $J_{\Phi }$ would effectively introduce a non-uniform weight function in the inner product integral, destroying the orthogonality and necessitating computationally expensive re-orthogonalization processes like the Gram-Schmidt procedure[9]. Therefore, maintaining a quasi-uniform $J_{\Phi }$ allows us to use standard bases with improved discrete orthogonality under uniform sampling.

2.5. Topology-Aware Canonical Annulus for Doubly Connected Domains

For an aperture with an internal obscuration, mapping to the unit disk is topologically invalid[25]. We therefore map to a concentric annulus

$$\mathbb{A}=\{w\in \mathbb{C}:R_{\mathrm{in}}\le |w|\le 1\}.$$

The inner radius $R_{\mathrm{in}}$ is determined by the conformal modulus of the physical domain. Assigning an arbitrary radius would introduce significant shear distortion[26].

We compute this conformal invariant by solving a harmonic Dirichlet problem on the physical domain ${\Omega }_{\mathrm{phy}}$[27,28]:

$$\left\{\begin{array}{cc}\Delta u=0,\hfill & \mathrm{in}{\Omega }_{\mathrm{phy}},\hfill \\ u=0,\hfill & \mathrm{on}\partial {\Omega }_{\mathrm{out}},\hfill \\ u=1,\hfill & \mathrm{on}\partial {\Omega }_{\mathrm{in}}.\hfill \end{array}\right.$$

(5)

On a triangular mesh, we discretize the Laplace-Beltrami operator $\Delta$ using the cotangent weights[29], leading to the stiffness matrix $\mathbf{W}$. We then solve the resulting sparse linear system associated with $\mathrm{W}$ under Dirichlet boundary constraints. The corresponding solution minimizes the Dirichlet energy, defined as follows[30]:

$$E={\int }_{{\Omega }_{\mathrm{phy}}}{\parallel \nabla u\parallel }^2\mathrm{d}\sigma \approx {\mathbf{u}}^{\mathsf{T}}\mathbf{W}\mathbf{u}.$$

(6)

For the canonical annulus, the analytic energy is given by $E=2\pi /ln(1/R_{\mathrm{in}})$. By equating the discrete and analytic energies, we obtain the unique radius:

$$R_{\mathrm{in}}=exp\left(-{\displaystyle \frac{2\pi }{E}}\right).$$

(7)

This ensures the canonical annulus preserves the conformal modulus of the original domain.

2.6. Wavefront Reconstruction and Numerical Stability

After obtaining $\Phi :{\Omega }_{\mathrm{phy}}\to {\Omega }_{\mathrm{can}}$, wavefront samples are mapped as[31]

$$w_i=\Phi (x_i,y_i),W_i=W(x_i,y_i).$$

We then fit W on the canonical domain using an orthogonal basis.

Simply connected (disk).

On $\mathbb{D}$, we use the standard Zernike basis ${\left\{Z_j\right\}}_{j=1}^M$ (Noll indexing)[4]. The expansion is

$$W\approx \sum _{j=1}^M{c}_j{Z}_j\left(w\right),w\in \mathbb{D}.$$

(8)

Doubly connected (annulus).

On $\mathbb{A}$, we use annular Zernike polynomials ${\left\{Z_j^{\mathrm{ann}}(·;R_{\mathrm{in}})\right\}}_{j=1}^M$ , which reduce to standard Zernike polynomials as $R_{\mathrm{in}}\to 0$[32]. The expansion is

$$W\approx \sum _{j=1}^M{c}_j{Z}_j^{\mathrm{ann}}(w;R_{\mathrm{in}}),w\in \mathbb{A}.$$

(9)

Discrete least squares.

With N samples, define the design matrix $\mathbf{H}\in {\mathbb{R}}^{N\times M}$ by

$$H_{ij}={\Psi }_j\left(w_i\right),$$

where ${\Psi }_j$ denotes either $Z_j$ (disk) or $Z_j^{\mathrm{ann}}(·;R_{\mathrm{in}})$ (annulus). The coefficients are obtained by a numerically stable least-squares solve (QR or SVD):

$$\mathbf{c}=\underset{\mathbf{c}\in {\mathbb{R}}^M}{arg\; min}{\parallel \mathbf{H}\mathbf{c}-\mathbf{W}\parallel }_2^2.$$

(10)

Conditioning and discrete orthogonality.

In exact arithmetic with continuous uniform sampling, these bases are orthogonal on their canonical domains. In discrete computations, stability is governed by the conditioning of $\mathbf{H}$ (or the Gram matrix $\mathbf{G}={\mathbf{H}}^{\mathsf{T}}\mathbf{H}$). Severe area distortion (high variance of $J_{\Phi }$) effectively introduces non-uniform sampling weights on ${\Omega }_{\mathrm{can}}$, degrading discrete orthogonality and increasing $\kappa \left(\mathbf{H}\right)$[6]. CSF-QCM leverages bounded quasi-conformal relaxation to mitigate extreme area distortion, thereby improving conditioning and preventing error amplification.

The overall procedure is summarized in Algorithm 1.

3. Results

We evaluate three representative freeform aperture types exhibiting pronounced geometric challenges, corresponding to specific optical metrology scenarios:

• Type I: a non-convex butterfly-shaped aperture, typical of the irregular interference regions encountered in speckle metrology,as shown in Figure 2;
• Type II: a high-aspect-ratio rounded rectangle, representing the geometry of primary or secondary mirrors in wide-field-of-view off-axis three-mirror anastigmat (TMA) systems,as shown in Figure 2;
• Type III: a highly eccentric doubly connected annulus, modeling the pupil in fundus aberration interferometry where the central macular region creates an off-center obscuration,as shown in Figure 2.

3.1. Mesh Distribution and Sampling Uniformity

In mapping-based wavefront fitting, local cell areas determine discrete quadrature weights and influence numerical stability. Excessive area compression or stretching leads to oversampling or undersampling and degrades least-squares conditioning.

3.1.1. Visual Assessment of Parameterization Meshes

Figure 3 illustrates that strict conformality can impose geometric rigidity, producing highly non-uniform meshes on irregular boundaries:

• Type I (Butterfly): Fornberg-type conformal mapping exhibits strong crowding near concave regions, yielding redundant sampling.
• Type II (Rounded rectangle): Schwarz–Christoffel (SC) mapping compresses the grid near the ends of the long axis due to crowding, producing disproportionate sampling.
• Type III (Annulus): Discrete Ricci flow preserves angles but can introduce severe area distortion in narrow eccentric gaps.

By contrast, CSF-QCM smooths high-frequency boundary features and redistributes arc length during CSF evolution, resulting in quasi-uniform meshes across all cases while preserving topology.

Algorithm 1: CSF-Guided Progressive Quasi-Conformal Mapping (CSF-QCM)

Require: Physical domain mesh ${\Omega }_{\mathrm{phy}}$ with vertices $\mathcal{V}$ and faces $\mathcal{F}$; Wavefront samples vector $\mathbf{W}$; Boundary definition $\partial {\Omega }_{\mathrm{out}}$ (and $\partial {\Omega }_{\mathrm{in}}$ if doubly connected); Tolerance ${ϵ}_{\mathrm{tol}}$. Ensure: Canonical mapping $\Phi$, Fitted coefficients $\mathbf{c}$.  Step 1: Topology-Aware Target Specification         ▹ See Sec. 2.5 1: if domain is simply connected then 2:     Set canonical domain ${\Omega }_{\mathrm{can}}\leftarrow \mathbb{D}$ (Unit Disk) 3: else    ▹ Doubly connected case 4:     Solve harmonic Eq. (5): $\Delta u=0$ s.t. ${u|}_{\partial {\Omega }_{\mathrm{out}}}{=0,u|}_{\partial {\Omega }_{\mathrm{in}}}=1$ 5:     Compute Dirichlet energy $E\approx {\mathbf{u}}^{\mathsf{T}}\mathbf{W}\mathbf{u}$           ▹Eq. (6) 6:     Compute conformal modulus: $R_{\mathrm{in}}\leftarrow exp(-2\pi /E)$     ▹Eq. (7) 7:     Set canonical domain ${\Omega }_{\mathrm{can}}\leftarrow \mathbb{A}$ (Annulus with radii $R_{\mathrm{in}},1$) 8: end if  Step 2: Boundary Geometric Regularization (CSF)        ▹See Sec. 2.3 9: Initialize $k\leftarrow 0$, boundary configuration ${\mathcal{C}}^{\left(0\right)}$ 10: repeat 11:     Compute curvature $\kappa$ and normals $\mathbf{N}$ 12:     Evolve boundary: ${\displaystyle \frac{\partial \mathcal{C}}{\partial t}}\leftarrow \kappa \mathbf{N}$              ▹Eq. (1) 13:     Update circularity metric $E_{\mathrm{circ}}^{(k+1)}$ 14:     $k\leftarrow k+1$ 15: until Stopping criterion Eq. (2) is met: $|\Delta E_{\mathrm{circ}}|<{ϵ}_{\mathrm{tol}}$ 16: Store sequence of boundary frames ${\left\{{\mathcal{C}}^{\left(t\right)}\right\}}_{t=0}^k$  Step 3: Progressive Interior Optimization           ▹ See Sec. 2.4 17: Construct Graph Laplacian $\mathbf{L}$ 18: for each time step $t=0\dots k$ do 19:     Update boundary positions ${\mathbf{P}}_B^{\left(t\right)}$ mapping ${\mathcal{C}}^{\left(t\right)}$ to $\partial {\Omega }_{\mathrm{can}}$ 20:     Solve linear system for interior positions ${\mathbf{P}}_I^{\left(t\right)}$: 21:        ${\mathbf{L}}_{II}{\mathbf{P}}_I^{\left(t\right)}=-{\mathbf{L}}_{IB}{\mathbf{P}}_B^{\left(t\right)}$                 ▹ Eq. (4) 22:     Compute Beltrami coefficient ${\mu }_{\Phi }$ for validation       ▹ Eq. (3) 23: end for 24: Construct final map $\Phi$ from ${\mathbf{P}}^{\left(k\right)}$  Step 4: Wavefront Fitting                ▹ See Sec. 2.6 25: Map sample points: $w_i\leftarrow \Phi (x_i,y_i)$ 26: if${\Omega }_{\mathrm{can}}=\mathbb{D}$then 27:     Construct $\mathbf{H}$ using Standard Zernike $\left\{Z_j\right\}$         ▹ Eq. (8) 28: else 29:     Construct $\mathbf{H}$ using Annular Zernike $\left\{Z_j^{\mathrm{ann}}(·;R_{\mathrm{in}})\right\}$       ▹ Eq. (9) 30: end if 31: Solve $\mathbf{c}={arg\; min\parallel \mathbf{H}\mathbf{c}-\mathbf{W}\parallel }_2^2$ via QR/SVD 32: return$\Phi ,\mathbf{c}$

3.1.2. Quantitative Statistics of Area Distortion

For each triangle $\tau$, define an area scaling factor

$$s_{\tau }={\displaystyle \frac{\mathrm{Area}(\Phi (\tau \left)\right)}{\mathrm{Area}\left(\tau \right)}}\approx J_{\Phi }.$$

Uniform sampling on ${\Omega }_{\mathrm{can}}$ is promoted when the distribution of $s_{\tau }$ avoids extreme tails. We summarize the distribution using the empirical cumulative distribution function (ECDF). Curves closer to a steep transition around a constant scale indicate more uniform area scaling.

Figure 4. ECDF comparison of the local area scaling factor $s_{\tau }$. Curves with reduced tail spread indicate improved area uniformity.

Figure 4. ECDF comparison of the local area scaling factor $s_{\tau }$. Curves with reduced tail spread indicate improved area uniformity.

3.2. Characterization of Quasi-Conformal Distortion and Discrete Orthogonality

The key mechanism of CSF-QCM is to implicitly accommodate local angular distortion to improve global area uniformity. This section examines the Beltrami coefficient distribution as a post-mapping indicator and the resulting recovery of discrete orthogonality.

3.2.1. Beltrami Coefficient Distribution

Strict conformality enforces ${\mu }_{\Phi }\equiv 0$, which can be overly restrictive on irregular boundaries. CSF-QCM relaxes this constraint via progressive quasi-conformal deformation. In geometrically constrained regions (e.g., concave turns in Type I or long-edge endpoints in Type II), moderate $|{\mu }_{\Phi }|$ values arise to relieve boundary-induced tension and reduce area distortion, while maintaining $|{\mu }_{\Phi }|<1$ to preserve diffeomorphism. Figure 5 illustrates the spatial distributions (a–c) and statistical histograms (d–f) of the Beltrami coefficients for all three aperture types.

3.2.2. Recovery of Discrete Orthogonality

When the mapping $\Phi$ induces highly non-uniform area scaling, the discrete orthogonality of the canonical basis is compromised. This degradation is visualized by the structure of the Gram matrix $\mathbf{G}={\mathbf{H}}^{\mathsf{T}}\mathbf{H}$. In an ideal scenario with uniform sampling, $\mathbf{G}$ should be an identity matrix (up to normalization). However, as shown in the top row of Figure 6, baseline conformal mappings produce Gram matrices with significant off-diagonal components (modal crosstalk), particularly for the high-aspect-ratio Type II and eccentric Type III apertures. This structure indicates that the basis functions have become numerically linearly dependent on the sampled grid, leading to an ill-conditioned inverse problem.

By contrast, CSF-QCM explicitly optimizes for sampling uniformity. As illustrated in the bottom row of Figure 6, the resulting Gram matrices are strictly diagonally dominant with suppressed off-diagonal terms, confirming the recovery of discrete orthogonality. We quantify this improvement using the condition number $\kappa \left(\mathbf{G}\right)$ for the first 36 Zernike modes (Table 1). While baseline methods yield condition numbers ranging from ∼162 to over 900, implying amplification of measurement noise, CSF-QCM consistently reduces $\kappa \left(\mathbf{G}\right)$ to single digits ($<\mathbf{8}$) across all aperture types. This near-optimal conditioning ensures that the wavefront fitting remains numerically stable and minimizes the propagation of potential measurement errors.

3.3. Wavefront Fitting Accuracy and Computational Efficiency

To validate the reconstruction fidelity, we generated a synthetic wavefront composed of the first 36 Zernike modes (using Noll indexing for simply connected domains and the corresponding annular basis for Type III), with a peak-to-valley (PV) amplitude of $1/20\lambda$. This setup mimics typical high-precision testing scenarios where minimizing residual fitting error is critical.

Figure 7 presents a comprehensive visual comparison of the reconstructed surfaces and their associated error maps. The baseline methods (Figure 7(b, c)) exhibit distinct boundary ringing and structural artifacts, particularly near high-curvature regions and narrow gaps. These errors are consistent with the ill-conditioning analysis in Section 3.2: the local sampling sparsity caused by conformal crowding prevents the orthogonal basis from resolving edge behaviors accurately. In contrast, CSF-QCM (Figure 7(d, e)) yields a highly accurate reconstruction with spatially homogeneous residuals. The absence of localized error spikes confirms that the quasi-uniform sampling strategy effectively stabilizes the numerical fit.

Table 2 quantitatively summarizes the Peak-to-Valley (PV) and Root-Mean-Square (RMS) residuals. Across all three aperture types, CSF-QCM significantly outperforms the baseline approaches. For visible wavelengths (e.g., $\lambda \approx 633$ nm), the method achieves nanometer RMS accuracy. Notably, for the most challenging Type II (rounded rectangle) and Type III (eccentric annulus) apertures, the RMS fitting error is reduced by over 90%, demonstrating the method’s robustness against complex boundary topologies.

Beyond accuracy, practical metrology demands computational efficiency. We evaluated the runtime on a standard PC (Intel i5 CPU, 16 GB RAM) using meshes with approximately 15,000 vertices. Table 3 details the computational cost broken down by processing stage. Although CSF-QCM introduces an iterative boundary evolution step, the total runtime remains competitive (∼3 seconds). This is achieved because the subsequent interior mapping relies on solving sparse linear systems, which is computationally inexpensive compared to the nonlinear optimization required by Ricci flow. The proposed framework thus offers a favorable trade-off, providing high-precision reconstruction with a computational cost suitable for routine laboratory testing.

4. Discussion

CSF-QCM improves accuracy and numerical stability for wavefront fitting on simply and doubly connected apertures by effectively managing the trade-off between conformality (angle preservation) and area uniformity (sampling regularity). Several aspects merit further investigation.

4.1. Extension to Multiply Connected Apertures

The current framework uses a capacity-based invariant to determine the canonical annulus for doubly connected domains. For more complex topologies (e.g., triply connected apertures, segmented mirrors, or pupils with multiple obscurations), the canonical domain is no longer a simple annulus. Possible extensions include circle-domain parameterizations leveraging fast boundary integral equation methods[33] or analytical approaches based on the Schottky–Klein prime function[34]. These modern numerical tools offer superior convergence rates compared to classical Koebe-type iterative constructions for general n-connected planar domains. CSF-based boundary smoothing remains applicable and can be combined with these multi-boundary canonicalization techniques to handle complex segmented pupil geometries.

4.2. Interaction Between Quasi-Conformal Distortion and Aberration Estimation

The present optimization emphasizes geometric uniformity, while mapping distortion may interact with finite sampling and specific aberration modes. A useful next step is to quantify the sensitivity of estimated coefficients to local quasi-conformal distortion, e.g., through a distortion–aberration sensitivity matrix. Such a model could enable distortion allocation strategies that prioritize regions most influential to the targeted aberration terms, conceptually analogous to the weighted Zernike decomposition strategies employed in high-contrast imaging.[35]

4.3. Acceleration Toward Real-Time Metrology

Although CSF-QCM achieves second-level runtime on typical meshes, online metrology may require remeshing and remapping when the valid aperture changes dynamically. The dominant cost arises from large sparse linear solves (Eq. (4)). While GPU-accelerated sparse factorizations can reduce latency, a more transformative direction is the adoption of physics-informed machine learning. Specifically, Physics-Informed Neural Networks (PINNs)[36] or Deep Operator Networks (DeepONet)[37] could serve as real-time surrogate models to predict the mapping functions directly by learning the underlying Laplace operator, potentially bypassing the iterative mesh generation process entirely.

5. Conclusion

We presented a CSF-guided progressive quasi-conformal mapping framework (CSF-QCM) to address the loss of orthogonality and numerical instability in wavefront fitting on non-circular freeform apertures. By introducing boundary evolution as a geometric preprocessing step, the framework constructs low-distortion parameterizations from irregular physical apertures—including non-convex, high-aspect-ratio, and doubly connected domains—to canonical computational domains.

The quantitative validations presented in this study highlight three key advancements. Firstly, the numerical stability is drastically improved; by managing the conformality–uniformity trade-off, CSF-QCM reduces the condition number of the Gram matrix from severe levels (e.g., $>900$ for butterfly apertures) to single digits ($<8$) across all tested geometries, effectively resolving the ill-conditioning caused by crowding. Secondly, this stability translates into superior reconstruction accuracy. The method eliminates boundary ringing artifacts and achieves nanometer-level precision (RMS $\approx 0.003\lambda$), reducing fitting errors by over 80% (up to 92.5% for high-aspect-ratio shapes) compared to baseline Fornberg and Ricci flow algorithms. Thirdly, the framework maintains computational efficiency, processing typical meshes (∼15k vertices) in approximately 3 seconds, which is faster than iterative Ricci flow and suitable for routine laboratory testing.

In summary, CSF-QCM provides a unified, robust, and fast preprocessing route for wavefront fitting over arbitrary apertures, addressing critical challenges in high-precision freeform metrology and the alignment of complex off-axis optical systems.