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.

Explore similar work

Translation-Invariant Tile-Based Phase Unwrapping with Residual-Weighted Multipath Averaging
Sep 11, 2026Antoine Moevus, Max MignotteImage Inverse Problems Inverse Problems
D´epartement d’informatique et de recherche opérationnelle, Universit´e de Montr´eal, Montr´eal, QC, Canada
Deeply unrolled alternating projections for phase retrieval
Oct 12, 2026Jiayi Yao, Vasileios Charisopoulos, Rebecca WillettImage Inverse Problems Deep Learning Optimization
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 Reconstruction from Phase with Untrained Neural Priors
Sep 25, 2026Ene Meco, Ahmet Enis CetinImage Reconstruction Image Inverse Problems
Dept. of Electrical and Computer Engineering, University of Illinois at Chicago, IL, 60607, USA

Sep 11, 2026J/K move · Enter open · S save

Translation-Invariant Tile-Based Phase Unwrapping with Residual-Weighted Multipath Averaging

Antoine Moevus, Max Mignotte

D´epartement d’informatique et de recherche opérationnelle, Universit´e de Montr´eal, Montr´eal, QC, Canada

Image Inverse Problems Inverse Problems Image Reconstruction

Open paper AskPDF arXiv

Phase unwrapping is a key step in interferometric and coherent imaging, where the physical quantity of interest is carried by a phase that the instrument delivers only modulo 2*pi. The difficulty in two dimensions is to separate the jumps caused by wrapping from those produced by noise, by true discontinuities, by under-sampling or by decorrelation. Spatial-domain and frequency-domain methods have both been studied extensively, each with advantages the other lacks; hybrid schemes combining the two remain scarce. We propose a semi-global tile-based strategy in which every tile is unwrapped in the frequency domain, through the Discrete Cosine Transform (DCT) and the least squares (LS) formalism of Ghiglia et al., the tiles being merged spatially. Unwrapping and noise filtering are performed jointly, which regularizes an otherwise ill-posed inverse problem. An error stays confined to the tile in which it arose; the tiling artifacts are removed by averaging over every shift of the grid and over the symmetries of the square; and that average is weighted by the Poisson residual each pass leaves behind, so that a pass whose tile boundaries fell on a discontinuity does not impose its seam on the result. Experiments on synthetic and real data, against four reference algorithms from four distinct families and under six complementary metrics, show that the proposed method matches or improves on the state of the art. One of these metrics, a corrected cyclic re-wrap residual introduced here, needs no ground truth and therefore remains available on real acquisitions.

Abstract

Explore similar work

Translation-Invariant Tile-Based Phase Unwrapping with Residual-Weighted Multipath Averaging

Deeply unrolled alternating projections for phase retrieval

Image Reconstruction from Phase with Untrained Neural Priors