Adaptive Tiling for Least-Squares Phase Unwrapping: Runtime and Accuracy

Antoine Moevus; Max Mignotte

eess.IVSep 22, 2026

Adaptive Tiling for Least-Squares Phase Unwrapping: Runtime and Accuracy

Authors: Antoine Moevus, Max Mignotte

Organizations: Department of Computer Science and Operations Research Université de Montréal

Abstract

Phase unwrapping estimates the missing multiples of $2π$ in measured phase images. For large images, tiling limits the size of local reconstruction problems and enables parallel processing. Adaptive tiling could further reduce the number of local problems and boundaries by retaining large tiles where little refinement is needed. We investigate whether this reduction makes reconstruction faster. We compare complete reconstruction time and accuracy for a regular grid, quadtree, and kd-tree partitions. We also evaluate nine criteria for deciding where quadtree tiles should be subdivided, including residue count, fringe density, and measures of phase variation, at different tile sizes and budgets. In single-threaded experiments on a heterogeneous image dataset, optimized adaptive partitions use fewer tiles but remain slower than the optimized grid, and some reconstructions lose substantial accuracy. Stage measurements explain why: constructing the partition and solving larger retained tiles outweigh the savings at tile boundaries. The criterion comparison also shows that more refinement does not consistently improve accuracy. These results motivate evaluating adaptive partitions by the complete time needed to reach a chosen reconstruction accuracy, including whether limited refinement can provide a faster approximate result.

Figures & tables

Figure 1: Regular and adaptive partitions on three $256\times 256$ phase fields: A, a generated interferometric synthetic aperture radar (InSAR) field; B, a field derived from tree-photograph intensities; C, a generated vortex field. Columns show the known generating phase (ground truth, GT), wrapped phase, signed cell residues (red: $+1$ , blue: $-1$ , white: zero), a grid of $8\times 8$ -pixel tiles, and a residue-driven quadtree. Each GT panel has its own color scale in radians; wrapped panels share the cyclic scale $[-\pi,\pi]$ . White lines mark tile boundaries. The quadtree has minimum side $s=8$ pixels and stopping budget $B=512$ leaves. It retains 268, 514, and four tiles in A–C, versus 1024 grid tiles in each case. A four-way split adds three leaves, allowing the budget to be exceeded by two. These unquantized examples illustrate allocation, not reconstruction accuracy or runtime.

Figure 2: Common pipeline for grid, quadtree, and kd-tree reconstruction. The partition sets the tile boundaries; local Poisson solves use the BRiDCT discrete cosine transform (DCT), with FFTW for unsupported shapes. Seam discrepancies are summarized by their median and median absolute deviation (MAD). A separate reliability tree connects the reconstructed tiles and determines their relative phase offsets. The complete-call timer includes preparation, all reconstruction stages, and temporary-buffer cleanup. Transform plans are reused in the main timing comparison; first-call latency is evaluated separately.

Figure 3: Priority-based quadtree refinement on a

32\times 32

image, with minimum tile side

s=8

pixels and stopping budget

B=7

leaves.

C

is an illustrative nonnegative tile score, and

L

is the leaf count. (a) The root covers the image. (b) Bisecting both sides creates four

16\times 16

children; the orange outline selects the highest-scoring eligible child. (c) Splitting that child creates four

8\times 8

tiles and brings the total to seven leaves. Refinement stops at the budget even though another tile remains eligible with positive score. Zero-score tiles are retained.

Family	Entries	Source and reference
Photographs	45	Image intensities mapped to known phase fields, including USC-SIPI images.
InSAR	69	19 generated fields; 50 InSAR-DLPU entries, including 30 with terrain-model references [ 12 ] .
Synthetic surfaces	18	Analytic or generated fields with known phase.
Published tests	4	Ghiglia–Pritt images and supplied reference arrays [ 5 ] .
Holography	40	20 synthetic and 20 experimental entries; the latter use supplied algorithmic wrap-count references [ 13 ] .
MRI simulations	32	Magnetic resonance imaging (MRI) simulations from the quantitative susceptibility mapping (QSM) challenge, with known phase and brain masks [ 14 ] .

Table 1: Composition of the 208-entry dataset. Counts denote test inputs, not independent acquisitions. Ten entries guide implementation choices; the other 198 provide the main comparison. Reference phases comprise 158 known generating fields, 30 terrain-model references, and 20 algorithmic references. Agreement with the last category does not establish agreement with independent ground truth.

$s$ (px)	Partition	Median $t/t_{G}$	$\sum t/\sum t_{G}$	$\Delta S_{\pi}$ (pp)	$\Delta$ RMSE (rad)
8	Quadtree	1.244	1.371	-0.20	+0.105
8	KD midpoint	1.228	1.381	-0.04	+0.133
8	KD dyadic	1.229	1.294	-0.06	+0.145
16	Quadtree	1.532	1.777	-1.02	+0.184
16	KD midpoint	1.578	1.796	-0.99	+0.213
16	KD dyadic	1.570	1.644	-1.01	+0.225

Table 2: Complete reconstruction time and reference agreement on 198 entries.

s

is the minimum adaptive tile side in pixels; the grid uses nominal

s\times s

tiles. For each entry,

t

and

t_{G}

are the median times over seven repeats for the adaptive method and its grid control. The median of

t/t_{G}

describes a typical entry, while

\sum t/\sum t_{G}

compares aggregate time across entries. Ratios above one mean slower reconstruction.

\Delta S_{\pi}

and

\Delta

RMSE are mean paired changes, adaptive minus grid, in percentage points (pp) and radians.

S_{\pi}

is the fraction of pixels with absolute error below

\pi

after removing a common mean phase offset. Thus positive

\Delta S_{\pi}

and negative

\Delta

RMSE indicate better reference agreement. All methods share BRiDCT, FFTW fallback, and reconciliation.

Figure 4: Complete reconstruction time relative to the grid and its stage breakdown. (a) Each adaptive method uses minimum tile side

s

and is paired with a grid of nominal side

s

on the same input. Per-entry times are medians over seven repeats. Dots show the median ratio across 198 entries; whiskers show the 10th–90th percentiles across entries, not confidence intervals. The dashed line at one marks equal time. (b) Stage times are summed across entries using their per-entry medians and divided by the summed median complete-call times. Partitioning includes residue computation and tree construction; local solves include data preparation, the Poisson right-hand side, and transforms. Boundary statistics compute seam medians and MAD; joining propagates offsets. The remainder includes cleanup, uninstrumented overhead, and differences caused by summarizing each stage separately.

Figure 5: Effect of the subdivision criterion on a fixed observation. All panels show the same unquantized $256\times 256$ tree-photograph phase field as row B of Fig. 1 , with white quadtree boundaries over the wrapped phase. Panels (a)–(i) use criteria A–I from Sec. 4 ; $L$ is the resulting number of tiles. Minimum side $s=8$ pixels and stopping budget $B=512$ leaves are identical throughout. A four-way split can exceed this budget by two leaves. The shared cyclic color scale gives phase in radians. These partitions illustrate how scores allocate a limited tile budget; they are not the unrestricted partitions used in the main runtime comparison.

	Full budget $B_{1}$			Half budget $B_{1/2}$
Criterion	$L$	$\Delta S_{\pi}$ (pp)	$\Delta$ RMSE (rad)	$L$	$\Delta S_{\pi}$ (pp)	$\Delta$ RMSE (rad)
Minimum side $s=8$ pixels
A. Residue count	302.5	+0.00	+0.000	302.5	+0.00	+0.000
B. Local net charge	419.5	-0.25	+0.266	419.5	+0.50	+0.236
C. Paired cut length	362.5	-0.42	+0.026	362.5	-0.18	+0.070
D. Fringe density	748	+0.07	+0.021	514	+0.37	-0.226

Table 3: Criterion comparison on 208 entries at minimum side

s

pixels. For original dimensions

H\times W

, the full budget is

B_{1}=\max(1,\lfloor H/s\rfloor\lfloor W/s\rfloor)

and the half budget is

B_{1/2}=\max(1,\lfloor B_{1}/2\rfloor)

. For a

256\times 256

image at

s=8

, these targets are 1024 and 512 tiles; a split may exceed the target by two.

L

is the median tile count.

\Delta S_{\pi}

(percentage points, pp) and

\Delta

RMSE (radians, rad) are mean paired changes relative to residue count at the same

s

and budget. Positive

\Delta S_{\pi}

and negative

\Delta

RMSE improve agreement. The zero control rows do not imply equal accuracy between budgets: absolute control means are given below. Budgets limit tiles, not error, and this table does not compare runtimes.

Family	Entries	Dimensions	Source and reference phase
Photograph-derived	45	$256\times 256$ , $512\times 512$ , $1024\times 1024$ ; $250\times 288$	Grayscale test images, including USC-SIPI images, scaled into phase fields; the generating field is the reference. Clean and noisy variants are included.
InSAR	69	$256\times 256$	19 locally generated interferograms with known phase; 20 simulated and 30 real TanDEM-X entries from InSAR-DLPU [ 12 ] . The real-data references derive from a digital elevation model.
Synthetic surfaces	18	$256\times 256$ , $257\times 257$ ; $458\times 152$	Analytic and generated phase fields, with the generating field retained as the reference.
Published test images	4	$257\times 257$ ; $458\times 157$ ; $458\times 152$	Ghiglia–Pritt test images [ 5 ] , using their supplied reference arrays.
Holography	40	$256\times 256$	20 synthetic fields with known phase and 20 experimental inputs from Gontarz et al. [ 13 ] . Experimental references are reconstructed from the supplied quality-guided phase-unwrapping (QGPU) wrap counts.
MRI simulations	32	$92$ – $167$ rows; $101$ – $155$ columns	Cropped slices from four QSM challenge simulation/noise combinations [ 14 ] , two echo times per view. Reference phase is $2\pi f\,\mathrm{TE}$ from the supplied frequency $f$ and echo time $\mathrm{TE}$ .

Table S1: Dataset families, original image dimensions, and reference-phase construction. Dimensions are height

\times

width in pixels before square padding; ranges cover the selected entries, not the entire source collections. Known generating references include both analytic fields and phases derived from photograph intensities. Terrain-model and algorithmic references are identified separately in the last column.

ID	Input	Size	Origin and role in the comparison
F1	Tree photograph	$256\times 256$	Source recorded as USC-SIPI; intensity-derived phase with fine texture.
F2	Portrait	$1024\times 1024$	Source recorded as USC-SIPI; larger intensity-derived phase image.
F3	Zernike field	$256\times 256$	Locally generated polynomial phase field; curved fringes.
F4	Peaks field	$256\times 256$	Locally generated smooth test surface; no residues in the input.
F5	Noisy baboon	$512\times 512$	Locally prepared noisy photograph-derived field; dense subdivision.
F6	Lena	$512\times 512$	Photograph-derived field retained in the local collection; intermediate subdivision.

Table S2: Ten inputs used for implementation choices and first-call measurements. Dimensions are original height

\times

width in pixels, before padding. F1–F10 identify entries within this report. The last column states their origin and the image property relevant to these tests. Source-collection citations identify provenance; the exact derived benchmark arrays are not publicly available.

Figure S1: Cell-score maps underlying the tree-image partitions in the main report’s Fig. 5 . Panels (a)–(i) correspond to criteria A–I. Color represents $\sqrt{a_{c}/a_{\max}}$ , where $a_{c}$ is the cell score and $a_{\max}$ is the maximum within that panel; a zero map would be displayed as zero. Independent normalization reveals each map’s spatial support but prevents comparison of absolute scores between criteria. Panel labels report the proportion or count of nonzero cells. The uniform control is one everywhere; its inset shows the resulting partition, not a second score map.

Side (px)	SJ baseline	FFTW	Accelerate	Ooura general	Ooura specialized
8	0.048	0.131	—	0.530	0.101
16	0.189	1.189	0.956	1.712	0.489
32	0.882	4.374	3.350	6.451	—
64	6.141	17.668	11.417	25.823	—
128	67.924	94.720	44.160	104.776	—
256	296.262	468.333	186.678	473.060	—

Table S3: Poisson spectral-solve time in

\mu

s per square tile on an Apple M3 Max, using one thread and single precision. Each value is the median of nine randomly ordered timing batches cycling through 16 inputs. The operation comprises a forward DCT, division by nonzero Poisson eigenvalues, and an inverse DCT, including normalization and data rearrangement; setup, allocation, and right-hand-side construction are excluded. SJ baseline is the Shao–Johnson implementation preceding BRiDCT, not a BRiDCT measurement. Bold marks the fastest tested implementation at each size. Dashes indicate sizes unsupported by the tested interface, not by the entire algorithm family. Smaller values are faster; these kernel timings are not complete reconstruction times.

Variant	Time ratio		$\Delta S_{\pi}$ (pp)		$\Delta$ RMSE (rad)
	$s=8$	$s=16$	$s=8$	$s=16$	$s=8$	$s=16$
Density $0.001$	1.027	1.016	-5.93	-0.78	+0.240	+0.090
Density $0.01$	1.061	1.124	-1.52	-1.30	+0.556	+0.016
Low-charge cutoff	0.991	1.023	-1.44	+0.46	+0.873	+0.106
Verified integration	1.039	1.012	-0.04	0.00	+0.002	0.000
Residue-balanced kd-tree	1.525	1.847	-0.70	-3.07	+0.420	-0.190

Table S4: Paired effects of implementation variants on ten inputs at minimum side

s

pixels. Time ratios compare each variant with the control specified below, using per-input median timings; the table reports the median ratio across inputs. Values above one are slower.

\Delta S_{\pi}

and

\Delta

RMSE are mean paired changes, variant minus control, in percentage points (pp) and radians.

S_{\pi}

is the fraction of pixels with absolute reference error below

\pi

after mean-offset removal. Positive

\Delta S_{\pi}

and negative

\Delta

RMSE indicate improvement.

N

is the processed pixel count. Rows from different timing sessions do not establish a ranking between variants.

$s$ (px)	Partition	First/warm	First/grid
8	Grid	2.245	1.000
8	Quadtree	2.151	1.219
8	KD dyadic	2.244	1.244
16	Grid	2.938	1.000
16	Quadtree	2.391	1.378
16	KD dyadic	2.459	1.373

Table S5: First-call overhead and adaptive/grid comparison on the ten exploratory inputs.

s

is the minimum adaptive side, or nominal grid side, in pixels. Each entry/configuration uses three fresh processes, each with one first call and seven warm calls. First-call time is the median over the three processes; warm time is the median of their seven-call medians. Columns report medians of the resulting per-entry ratios. First/warm compares the same method before and after initialization; first/grid compares its first call with the grid’s first call at the same

s

. Ratios above one mean longer reconstruction. Process startup, input loading, and encoding are excluded.

Image Reconstruction

Image reconstruction recovers images, volumes, or scenes from incomplete, noisy, indirect, or compressed measurements using physical models and learned priors.

See all 274 papers

Sep 11, 2026

Tile-Based Phase Unwrapping with Reliability-Ordered Joining and Residual-Weighted Multipath Averaging

Phase unwrapping recovers a continuous phase map from measurements known only modulo 2π, in the presence of noise and discontinuities. We propose a tile-based method that combines local spectral reconstruction, reliability-ordered tile joining, and residual-weighted multipath averaging. Each tile is reconstructed by a discrete cosine transform (DCT) least-squares solve, and robust measurements along tile borders guide their joining along a spanning tree. We repeat this reconstruction on shifted tile grids and on rotated and reflected copies of the image, which changes tile borders and joining paths. A residual-weighted average combines the resulting reconstructions, which are computed in parallel and averaged on the fly without being stored. We also introduce the corrected circular rewrapping error, a congruence measure: after removing a global phase offset, it measures how far the rewrapped reconstruction is from the observed phase. It applies to least-squares and path-based methods, even without an unwrapped reference. Used with reference-based errors, it separates agreement with the measurement from recovery of the correct number of cycles. Experiments on radar-interferometric, magnetic-resonance, holographic, and constructed phase images compare our method with seven reference algorithms. Among the methods that complete every test image, ours has the lowest mean fraction of pixels with an error above π. On the images all methods complete, it lies on several speed-accuracy Pareto fronts. Error maps also reveal isolated cycle errors that average scores or congruence alone can hide.

Antoine Moevus, Max Mignotte

The authors are with the D´epartement d’informatique et de recherche opérationnelle, Université de Montréal, Montréal, QC, Canada.

Image Inverse Problems Inverse Problems Image Reconstruction

Ask PDF

Oct 12, 2026

Deeply unrolled alternating projections for phase retrieval

Deep unrolling methods for inverse problems augment traditional reconstruction algorithms with learned priors to improve their performance. These methods have been developed under the implicit assumption that the geometry of the measurement operator does not change between training and deployment; when this is no longer the case, the quality of obtained reconstructions can degrade significantly. In this paper, we argue that this is not an inherent limitation of the deep unrolling paradigm and can be mitigated by a careful choice of the underlying iterative algorithm. We propose a deep unrolling framework for the phase retrieval problem that augments the alternating projections method with a learned denoising step. Our method obtains high-quality reconstructions without the need for retraining in settings where conventional model-based alternatives struggle, such as when the dynamic range of the measurement operator changes significantly at test time.

Jiayi Yao, Vasileios Charisopoulos, Rebecca Willett

Electrical & Computer Engineering, University of Washington · NSF-Simons National Institute for Theory and Mathematics in Biology · Department of Computer Science, University of Chicago +3

Image Inverse Problems Deep Learning Optimization

Ask PDF

Sep 25, 2026

Image Reconstruction from Phase with Untrained Neural Priors

Fourier phase encodes important spatial image structure, but recovering an image without measured spectral magnitude requires additional constraints and leaves absolute intensity ambiguous. We propose a projection-based two-stage framework that combines Fourier-phase and spatial-support constraints with an image-specific neural prior. The first stage alternates constraint enforcement with regularized neural-prior updates, while the second performs phase/support refinement alone with guaranteed convergence. We evaluate two neural-prior implementations on the same 77 microscopy images and compare them with a constraint-only baseline. After 500 final refinement passes, the best-performing variant achieves 31.41 dB pooled PSNR, 35.75 dB mean PSNR, and 0.9531 mean SSIM, improving pooled PSNR by 1.51~dB and reducing pooled MSE by 29.3% relative to the baseline. The results demonstrate the benefit of combining neural guidance with explicit constraint refinement at the evaluated iteration budget, while showing that lower phase residual alone does not guarantee greater reconstruction accuracy.

Ene Meco, Ahmet Enis Cetin

Dept. of Electrical and Computer Engineering, University of Illinois at Chicago, IL, 60607, USA

Image Reconstruction Image Inverse Problems

Ask PDF

Abstract

Image Reconstruction

Explore similar work

Tile-Based Phase Unwrapping with Reliability-Ordered Joining and Residual-Weighted Multipath Averaging

Deeply unrolled alternating projections for phase retrieval

Image Reconstruction from Phase with Untrained Neural Priors