Adaptive Tiling for Least-Squares Phase Unwrapping

Let $\psi$ denote the wrapped observation and $\mathcal{W}$ reduce an angle to $(-\pi,\pi]$. A cell consists of four neighboring pixels arranged in a square. Compute the principal differences $\mathcal{W}(\psi_{q}-\psi_{p})$ along its four edges and sum them with signs following a consistent orientation around the cell. A nonzero sum, measured in multiples of $2\pi$, is a residue: the four observed differences cannot all belong to one consistent phase field. For the main runtime comparison, the tile score is the sum of absolute residue charges over cells wholly inside it, excluding cells crossed by its boundary. Alternative scores are examined in Sec. 4.

The quadtree starts with the image rectangle. A split bisects both sides and replaces one terminal tile, or leaf, with four children (Fig. 3). Both sides must be at least $2s$, where $s$ is the minimum permitted tile side. With a tile budget, the largest-scoring eligible leaf is split first. Refinement stops when the budget is reached or no eligible positive-score leaf remains. The final runtime comparison removes the budget restriction: every eligible tile containing a residue is split. In that setting, depth-first traversal produces the same leaves without maintaining their priority in a queue.

With a consistent oriented-edge convention, a residue-free rectangle admits an exact fit to its observed differences in exact arithmetic: zero circulation around every cell makes the integral between two pixels independent of the path. This is a local consistency statement, not a guarantee of the true phase. For example, slopes of $3\pi/2$ and $-\pi/2$ per pixel give the same wrapped observations and no residues. Recovering the physical phase requires additional sampling assumptions [9]. Inconsistencies can also lie across tile boundaries and therefore escape each tile’s internal test.

A kd-tree splits a rectangle into two children with one horizontal or vertical cut. It can refine one direction without refining the other, which may suit elongated structures. The cut is eligible only if both children retain at least $s$ pixels along the divided axis. We test KD midpoint, which bisects the longest eligible side, and KD dyadic, which places the cut at the largest power of two no greater than half that side, clamped to respect the minimum size. The latter favors dimensions supported by fast transforms without guaranteeing that every resulting tile has power-of-two dimensions.

Both kd-trees use the same residue stopping test as the quadtree, so their comparison measures the effect of cut geometry and the resulting tile shapes. We also examine a residue-balanced rule, which places the cut near the median of the residue distribution instead of dividing side lengths. Its measured cost and accuracy are reported with the other geometric comparisons in Sec. 3.3.

Let $g_{T}$ collect the observed principal differences inside tile $T$, and let $D_{T}$ compute neighboring differences of a candidate phase $u$. The local least-squares problem is

When the observed differences are inconsistent, this fit finds a compromise between them; a small fitting residual does not by itself guarantee agreement with the reference phase. The normal equations form a discrete Poisson system with no connections outside the tile. A discrete cosine transform (DCT) diagonalizes this system: transform the right-hand side, divide by the nonzero eigenvalues, and apply the inverse transform [4]. The constant coefficient is set to zero because local differences do not determine an additive offset.

The complete-runtime comparison in Sec. 3 uses BRiDCT [10], our high-performance single-precision DCT implementation, for both the grid and adaptive partitions. Based on the Shao–Johnson factorization [11], it supports square and rectangular tiles with power-of-two side lengths from 8 to 1024 pixels. This range is particularly useful for quadtree reconstruction, where repeated bisection produces tiles of widely varying sizes. Other shapes use cosine transforms from the FFTW library. Transform plans, which store reusable setup information, and working buffers are cached by shape.

After each local solve, phase values are rounded to multiples of $\delta=2\pi/256$ radians, using nearest-value rounding with ties to even. The remaining additive offsets must then be reconciled across shared boundaries. We follow the fixed-grid boundary-difference construction of Moevus and Mignotte [6], and account for unequal seam lengths in the median/MAD reliability rule below.

Let $u_{a},u_{b}$ be the local solutions on neighboring tiles, with unknown additive offsets $o_{a},o_{b}$. A seam is their shared boundary. For adjacent pixel pairs $(p_{k},q_{k})$ crossing it from tile $a$ to tile $b$, form

Its median $m_{ab}=\operatorname{median}_{k}d_{k}$ estimates $o_{b}-o_{a}$. The median absolute deviation (MAD), $D_{ab}=\operatorname{median}_{k}|d_{k}-m_{ab}|$, measures disagreement along the seam rather than the size of its offset. With seam length $n_{ab}$ (the number of pixel pairs), define

where $\eta=\delta$ is one phase-quantization step and keeps the weight finite when the seam discrepancies agree exactly. The factor 1.4826 expresses MAD on a Gaussian standard-deviation scale. The seam-length factor is heuristic: correlated boundary samples do not justify treating this weight as calibrated uncertainty.

Reliability-tree joining treats tiles as graph vertices and shared boundaries as graph edges. A maximum-weight spanning tree connects every tile without cycles while maximizing the sum of selected reliability weights. Median offsets are propagated from a fixed root along this tree. This joining tree is distinct from the spatial subdivision tree: it determines relative phase offsets, not tile boundaries.

Algorithms 2.4 and 2.4 separate the decisions that are varied from the operations held fixed. The score-map builder $A$ selects where refinement is useful, while the split rule $G$ determines its geometry. Section 4 defines the candidate scores; Sec. 2.2 defines the binary alternatives.

Algorithm 2.4 uses a summed-area table to obtain each rectangular score from four cumulative entries. It splits the highest-scoring eligible tile while retaining zero-score and ineligible leaves. Checking the budget before a split permits an overshoot of two leaves for quadtrees, but none for binary trees. With no budget limit, depth-first traversal gives the same leaves without a priority queue.

Algorithm 2.4 accepts either this adaptive partition or a regular grid. It solves the local problems, estimates seam offsets, and connects the tiles in decreasing order of reliability without creating cycles. A fixed root supplies the additive reference; assembly, minimum subtraction, and rounding complete the output. Ground truth is used only for subsequent accuracy evaluation.

The dataset contains 208 entries from six image families (Table 1). Ten inputs covering smooth fields, dense texture, and non-square dimensions guide implementation choices, while the other 198 evaluate the selected configurations. Because the dataset includes correlated cases and had already been used in exploratory work, this larger comparison evaluates fixed configurations rather than providing an independent test of generalization.

The exploratory subset comprises five photograph-derived inputs (trees, a portrait, noisy baboon, Lena, and a Barbara crop), two generated surfaces (Zernike and peaks), one synthetic InSAR field, one experimental holographic input, and one MRI simulation. Their dimensions range from $126\times 134$ to $1024\times 1024$ pixels. The complete dataset also includes narrow rectangular test fields and non-power-of-two dimensions, so the comparison exercises both the BRiDCT and fallback paths.

The reconstruction timings use one thread on an Apple M3 Max. For these experiments, observations are encoded in 256 intensity levels per cycle, and local and assembled solutions are quantized to the same phase step. Non-square inputs are padded with zero intensity to a square of side $\max(H,W)$ for every method, where $H$ and $W$ are the original height and width. The whole padded domain is reconstructed; output is cropped before evaluation, and MRI accuracy uses the supplied brain mask. These choices are shared controls, but they limit direct extrapolation to floating-point observations or reconstruction only within a mask.

The timer surrounds the complete C++ reconstruction call, including partition construction, local solves, seam processing, output assembly, and temporary-buffer destruction. Input loading, encoding, process startup, and accuracy evaluation are excluded. Transform plans are prepared during an untimed warm-up. Seven randomly ordered repetitions are recorded for every image and configuration at minimum sides $s=8$ and $s=16$, without an adaptive leaf-budget limit. Grids use nominal side $s$, with partial tiles at image edges. First-call latency is evaluated separately.

For each image and method we take the median repeated time, then form the adaptive/grid ratio. Table 2 reports the median of those paired ratios and the ratio of summed case times, which gives greater weight to expensive images. A time ratio of 1.25 therefore means 25% slower reconstruction.

Accuracy is evaluated on the original image, or the supplied MRI mask. We subtract the mean reconstruction-minus-reference difference over those pixels to remove one common additive offset before computing the fraction of evaluated pixels with absolute error below $\pi$, denoted $S_{\pi}$, and the root-mean-square error (RMSE). Mean paired changes in these measures show whether a difference in speed is accompanied by a difference in reference agreement. For example, a success change of $-1$ percentage point means a smaller fraction of pixels within tolerance.

We reduce the cost of constructing the tree, processing its interfaces, and solving its tiles. Tree construction is specific to the adaptive methods; improvements to local solves and reconciliation are also given to the grid.

For unrestricted residue refinement, depth-first traversal replaces the priority queue because processing order no longer affects the final leaves. Residue computation and cumulative summation are combined in one pass, avoiding a separate residue-map pass. Since the input is encoded on 256 levels, residue evaluation uses integer phase differences. A lookup table handles differences of exactly half a cycle consistently with the reference wrapping convention.

For reconciliation, we collect each shared boundary as a contiguous interface. Finding its median and MAD requires selecting central values rather than sorting every sample. A disjoint-set structure tracks the relative offsets of tile groups as they are joined, avoiding a separate propagation traversal. Finally, FFTW replaces the slow direct-transform fallback on shapes unsupported by BRiDCT.

To measure the benefit of these changes, we compare with a quadtree implementation that uses a priority queue, ordered-map seam collection, and direct rectangular fallbacks. Both implementations already use BRiDCT on supported shapes. The changes reduce median paired runtime by about 51% at $s=8$ and 42% at $s=16$. The comparison with the grid below therefore evaluates adaptivity after these avoidable costs have been reduced.

Two further shortcuts also fail to reduce paired runtime on the ten exploratory inputs. Direct integration with an explicit consistency check is about 1–4% slower than DCT solves, while searching each tile only until its first residue is found is about 3–4% slower than building the fused residue table. The retained pipeline therefore keeps DCT solves and the full table.

The optimized grid remains faster (Table 2). Quadtree reconstruction takes 24% longer in median paired time at $s=8$, and 53% longer at $s=16$. Both kd-tree variants also remain slower than the grid. Across the 198 entries, no tested adaptive configuration has a lower per-image median time than the optimized grid at either minimum size. Figure 4 shows the spread across inputs and the contribution of each stage.

The two runtime summaries reveal different effects of binary subdivision. At $s=8$, both binary variants are close to the quadtree in median runtime; at $s=16$, neither improves that statistic. Dyadic cuts reduce summed time relative to midpoint cuts more than they change the median ratio, suggesting that their benefit is concentrated among expensive inputs. This advantage over midpoint cuts remains insufficient to beat the grid. Binary subdivision also has a construction cost: for the same number of leaves $L$, it requires $L-1$ splits, compared with $(L-1)/3$ for a full quadtree. Greater directional freedom therefore comes with additional node-processing work.

Balancing residues across a cut does not improve this result. On the ten exploratory inputs, the residue-balanced rule takes about 1.53 and 1.85 times the midpoint kd-tree runtime at $s=8$ and $s=16$, respectively, in median paired time with otherwise matched reconstruction. It also reduces mean success at both sizes. The broader comparison therefore retains the simpler midpoint and dyadic rules.

The additional runtime is not accompanied by a consistent improvement in reference agreement. Mean success changes are modest, but averages hide severe failures: two related shear-field test images fall from $S_{\pi}=100\%$ with the grid to $S_{\pi}=0\%$ with each adaptive method. Thus, after removing the common phase offset, none of their adaptive pixels lies within $\pi$ of the reference. Mean RMSE increases for all three adaptive configurations. The shear failures also occur with the quadtree before the construction optimizations, so they are not introduced by the faster construction routine. These observations do not prove that every adaptive reconstruction is worse; they rule out a claim that the tested allocation preserves accuracy across the dataset.

Partition construction accounts for about 24–25% of summed quadtree time, and local reconstruction for 47–54% (Fig. 4). A regular grid avoids image-dependent partitioning. Adaptive tiling reduces the number of interfaces, but every pixel still contributes to a local solve, and retaining larger tiles can increase the work per pixel. For a fast transform on $n_{t}$ pixels, the work scales as $O(n_{t}\log n_{t})$. Combining small tiles increases the transform size while leaving the total number of pixels unchanged, so fewer local problems need not mean less transform work.

The complexity follows the same separation of work. Let $N$ denote the processed pixel count, $L$ the number of leaves, $P$ the number of pixel pairs crossing tile boundaries, and $E$ the number of adjacent tile pairs. Residue preprocessing costs $O(N)$, after which the unrestricted traversal creates $O(L)$ nodes and sorts leaves into a deterministic order in $O(L\log L)$. The transform work over all tiles is $O(\sum_{t}n_{t}\log n_{t})$. Replacing a direct separable transform on an $h\times w$ tile, whose cost is $O(hw(h+w))$, by a fast DCT avoids an additional penalty on unsupported rectangular shapes.

Fewer tiles mainly help the boundary stages: seam processing is linear in $P$ on average for the selection routine used, and reliability sorting costs $O(E\log E)$. Pixel preparation and output assembly still require $O(N)$ work, while storage is $O(N+E+L)$ apart from cached transform plans. Thus reducing interfaces leaves both the initial residue pass and the work over all pixels in place.

This distribution of work limits what further DCT acceleration can achieve on its own. The local-solve stage includes input extraction and right-hand-side construction as well as the transforms themselves. A faster DCT changes only part of that stage, and the same improvement is available to the grid. A useful optimization must therefore be judged by the complete paired runtime, rather than a kernel speedup in isolation.

For a single reconstruction, initialization adds a cost absent from the warm timings. On the ten exploratory inputs, median first-call/warm-call ratios are about 2.2–2.9. Quadtree and dyadic kd-tree remain slower than the grid in paired first-call time; process startup, loading, and encoding are excluded.

A criterion assigns a nonnegative score to each cell, and a tile sums scores over its interior cells. The resulting score serves two purposes: zero stops refinement, while positive values determine priority when a leaf budget limits subdivision. Consequently, a sparse score can leave large zero-score regions untouched, whereas a score that is positive everywhere continues refining even where the local fit is already consistent.

With fixed split locations and no budget limit, the stopping test determines the leaves; positive priorities affect only processing order. A restricted budget makes those priorities matter, as Fig. 5 illustrates. Criteria A–C use residues, D–F describe phase variation, G–H test phase jumps, and I ignores image content. Their neighborhoods and thresholds are fixed across entries, without a claim that these settings are optimal.

For each entry, we cross the nine criteria with minimum sides $s\in\{8,16\}$ pixels and the full and half budgets defined in Table 3. This gives 7488 reconstructions.

For this allocation and accuracy study, every criterion uses the same double-precision DCT solver, local and final quantization, and median/MAD reliability-tree reconciliation. Scores and local right-hand sides use the original, unquantized observations without square padding. All nine criteria are compared under these same conditions. Their error differences describe the effect of allocation within this study and should not be combined numerically with the encoded-input BRiDCT results of Sec. 3.

The two error measures in Table 3 need not improve together: more pixels can fall below the $\pi$ threshold while the largest errors increase. Nor do equal budgets imply equal tile counts, because a sparse criterion may stop before using its budget.

The absolute control values show that changing the budget can affect the two error measures differently. At $s=16$, halving the budget lowers mean success but also lowers mean RMSE. A reduction in one error statistic therefore does not establish that accuracy is preserved.

Dense scores often spend more tiles without improving reconstruction. At $s=8$ and budget $B_{1}$, residue count produces a median of 302.5 tiles, compared with 1024 for gradient variance, Laplacian magnitude, and the uniform control. Their mean success is about 1.50–1.53 percentage points lower, and mean RMSE is about 0.36 radians higher. Conversely, stopping almost everywhere is not sufficient: the near-$\pi$ criterion has a median of one tile but loses 5.50 percentage points of success. Tile reduction must therefore be assessed together with error.

Residue count is not the most accurate criterion in every setting. At $s=8$ and budget $B_{1/2}$, fringe density improves mean success by 0.37 percentage points and reduces mean RMSE by 0.23 radians, while its median tile count rises from 302.5 to 514. At $s=16$ with the same budget fraction, fringe density also improves both error measures, although both median counts are 130. Local net charge and gradient variance sometimes improve success while increasing RMSE. These findings support residue count as an interpretable baseline that often uses fewer tiles, rather than a universal accuracy optimum.

The criterion comparison shows how allocation affects error, but does not establish whether coarser partitions save time. We therefore test relaxed residue stopping on the ten exploratory inputs, measuring whether allowing a few residues inside a tile reduces reconstruction time without losing reference agreement. The density rules stop a tile when its residue count per interior cell is at most $0.001$ or $0.01$. A second rule stops when the count is below $0.25(h+w)$ for tile dimensions $h,w$, provided the larger side is at most $4s$. Each is paired with zero-residue stopping at $s=8$ and $s=16$, with five repeated timings per input. Both arms use the same joining rule and permit direct integration in place of a DCT solve only after checking consistency against every internal edge. These timings exclude temporary-buffer destruction and are interpreted only within this comparison.

These tests do not yield a consistent time–accuracy improvement. Across the two minimum sizes, the density rules take 1.6–12.4% longer in median paired time and reduce mean success by 0.78–5.93 percentage points. The perimeter-based cutoff is about 0.9% faster at $s=8$ and 2.3% slower at $s=16$, while mean RMSE increases at both sizes. Allowing additional local inconsistency therefore does not by itself produce a useful speedup in these tests. A broader study should instead measure the fastest reconstruction available at each acceptable error level.

Let $\phi$ be a phase field and $\psi=\mathcal{W}(\phi)$ its wrapped observation, where $\mathcal{W}$ reduces angles to $(-\pi,\pi]$. An oriented edge $e=(p,q)$ joins two neighboring pixels. On tile $T$, collect the observed principal differences $g_{e}=\mathcal{W}(\psi_{q}-\psi_{p})$ into $g_{T}$. The operator $D_{T}$ computes the corresponding differences of a candidate phase, so $(D_{T}u)_{e}=u_{q}-u_{p}$. The local problem is

Its normal equations, $D_{T}^{\mathsf{T}}D_{T}u=D_{T}^{\mathsf{T}}g_{T}$, form a discrete Poisson system with the natural boundary conditions of the local fit. No edge outside the tile is included. The DCT diagonalizes this system [4]. Its zero-frequency coefficient is fixed to zero because neighboring differences do not determine an additive constant. Local solutions are rounded to a phase step $\delta=2\pi/256$ radians, with ties to even. After offset assembly, the output minimum is subtracted and the values are rounded to the same step. Input encoding is used in the timing comparison, whereas the criterion study retains the original floating-point observations.

For a four-pixel cell $c$, let $r_{c}=(2\pi)^{-1}\sum_{e\in\partial c}\epsilon_{ce}g_{e}$, with signs following the oriented cell boundary. The residue criterion is

using only cells wholly inside the tile. In exact arithmetic on a hole-free rectangle, $C_{r}(T)=0$ exactly when these differences admit a consistent potential: zero circulation on all cells makes path integration independent of the path. The unrounded least-squares residual is then zero, and the solution is unique up to a constant.

This concerns observed differences, not the true phase. Slopes $3\pi/2$ and $-\pi/2$ per pixel have identical wrapped observations and zero residues. Recovering the physical phase needs additional sampling assumptions [9]. Cells crossed by tile boundaries are also absent from the local count, leaving inconsistencies for reconciliation.

The reconstruction and median/MAD reconciliation are defined in the main report, Sec. 2.3 and Algorithms 2.4–2.4. In the criterion comparison, the size-based budget is $B_{1}=\max(1,\lfloor H/s\rfloor\lfloor W/s\rfloor)$; the restricted budget is $B_{1/2}=\max(1,\lfloor B_{1}/2\rfloor)$. A regular grid includes partial edge tiles and has $\lceil H/s\rceil\lceil W/s\rceil$ tiles. On non-dyadic images, bisection need not reproduce that grid.

Let $N$ be the processed pixel count, including padding, and $L$ the final tile count. Residue evaluation and a summed-area table cost $O(N)$ and give constant-time tile-score queries. Priority refinement costs $O(L\log L)$ in queue operations. The final unrestricted traversal instead creates $O(L)$ nodes, then sorts leaves into a deterministic order in $O(L\log L)$; it removes queue overhead without changing this overall bound. Fixed-neighborhood residue, fringe, variance, and divergence maps are linear in $N$, whereas paired-path criteria also require residue searches and path construction.

Fast-transform work on $n_{t}$ pixels scales as $O(n_{t}\log n_{t})$, yielding $O(\sum_{t}n_{t}\log n_{t})$ local work. A direct separable fallback on an $h_{t}\times w_{t}$ rectangle costs $O(h_{t}w_{t}(h_{t}+w_{t}))$; replacing it can dominate aggregate speedups on unusual shapes. Input extraction, right-hand-side construction, pixel ownership, and output assembly additionally require linear pixel work regardless of tile count.

Let $P$ count pixel pairs across tile boundaries, $E$ adjacent tile pairs, and $p_{e}$ the length of seam $e$. The final collector walks contiguous interfaces in $O(P)$ and computes median/MAD by selection, with average linear work in $p_{e}$ for the implementation used. Interfaces are first sorted by their tile identifiers to fix tie order, then stably sorted by decreasing reliability. These two sorts cost $O(E\log E)$. Weighted union–find propagates offsets while joining groups in $O(E\alpha(L))$ amortized work, where $\alpha$ is the inverse Ackermann function. The spatial and reconciliation storage is $O(N+E+L)$, excluding reusable transform plans. An ordered map over boundary samples would instead add $O(P\log(E+1))$ grouping work.

For the same leaf count $L$, a binary tree performs $L-1$ splits, compared with $(L-1)/3$ for a full quadtree. The tested residue-balanced cut uses binary searches of summed-area queries, adding $O(\log\ell_{v})$ work at split $v$ along an axis of length $\ell_{v}$. The rejected early-exit search caps inspected cells at a fixed multiple of $N$ before building a prefix table, bounding repeated scanning but adding work before fallback. These bounds describe costs; the stage measurements determine their practical balance.

Tile scores sum cell scores $a_{c}$ over cells wholly inside the tile. These definitions instantiate the criteria in the main report, Sec. 4. Table 3 in the main report compares their allocation and accuracy. The complete-runtime comparison retains residue count (A), whose additional stopping variants are defined in Sec. S4. Positive global rescaling changes neither priority order nor the zero test.

Figure S1 shows the scores underlying the main report’s partition comparison (Fig. 5). Positive rescaling leaves the partition rule unchanged, so each map is displayed with its own normalization. This visualization makes the difference between sparse and nearly everywhere-positive criteria visible without claiming a runtime or accuracy ranking.