unwrapdelay

Description:

Resolves sidelobe ambiguity in a rapidtide delay map by phase unwrapping.

The similarity function of an LFO signal has sidelobes spaced by the dominant period of the waveform (rapidtide writes its estimate to XXX_autocorr_sidelobetime_passN.txt). When a sidelobe exceeds the main lobe in a voxel, peak picking assigns a delay wrong by nearly one full period. These errors are large, and because neighboring voxels see the same waveform they fail together, in coherent patches rather than isolated speckles - which is what a median filter based despeckler struggles with, since the neighbors it votes with are wrong too.

The framing is phase unwrapping: the absolute delay is ambiguous modulo the sidelobe period, while the local gradient is not. Each voxel offers candidate delays (the local maxima of its similarity function), and quality guided region growing assigns each voxel the candidate closest to the value predicted from its already assigned neighbors, starting from the least ambiguous voxels.

Before you use this, check two things:

Measure your sidelobe. It depends on the LFO spectrum of the particular acquisition and is not a constant. rapidtide writes its estimate to XXX_autocorr_sidelobetime_passN.txt and records acsidelobelag_passN in the runoptions, leaving it None when it finds none; the autocorrelation itself is saved as XXX_desc-autocorr_timeseries. Two real examples, measured from that file within the search range:

dataset

strongest sidelobe

normalized amplitude

sub-RAPIDTIDETEST

13.25 s

0.026

a CVR study subject

12.00 s

0.001 (i.e. none)

With no sidelobe there is no periodic ambiguity to resolve, and large delay changes are being driven by something else - usually noise repicking the peak more or less at random inside the search range. This program will still “repair” those, but it is then smoothing, not unwrapping. The run prints a diagnostic: a genuine alias moves delays in ONE direction by close to one period (the test dataset gives 0% positive changes, median 14.01 s, IQR 5.24 s), while noise repicking moves them BOTH ways over a broad range (the CVR subject gives 13% positive, median 12.06 s, IQR 5.91 s, and triggers a warning).

Long delays are not automatically wrong. Vascular pathology can produce genuine delays of tens of seconds. What separates real long transit from a repicking error is not magnitude but spatial behaviour - real delay varies smoothly from surrounding tissue, an error jumps. The smoothness prior used here cannot tell the difference at a discontinuity, so on a patient with focal delay pathology this program may erase the finding of interest. Inspect the unwrapchanged mask before trusting the output.

Measured performance:

On the test dataset, processed without --despecklepasses and without --robustdelayfit, scored against the despeckled run:

method

sidelobe errors

rapidtide raw peak picking

380 (1.45%)

unwrapdelay

89 (0.34%)

a 77% reduction, in under a second. Under controlled injection - corrupting the similarity function so the sidelobe outranks the main lobe - it recovers about 97% of injected errors even when they are placed in the lowest maxcorrsq quartile, where real errors actually live.

Error propagation, and what fixed it:

Region growing propagates its own mistakes - a voxel assigned a wrapped value drags its correct neighbors to match - so the method creates new wrapped voxels while fixing old ones. Scoring only the errors present in the input hides this and flatters any aggressive method.

Measured on HCP data with a strong measured sidelobe (13.2 s, amplitude 0.25), using an alias-specific detector (a voxel is wrapped if it sits one sidelobe period from its local median, which does not penalise genuine smooth long delays):

Predicting from the median over every already assigned neighbor, rather than from whichever single neighbor popped off the heap, cut new wraps by 44% and moved the method from a tie with despeckling to 37% fewer residual wrapped voxels. It also raised cross-stream consistency from 98.93% to 99.37%.

--minconfidence was tried for the same purpose and does not work: net totals of 1685, 2449, 1982, 1706, 1685 for floors of 0.0, 0.25, 0.5, 0.75, 0.9. Every nonzero setting is neutral or worse, because excluding sources shrinks the consensus and makes the median less robust. It defaults to off and is retained only to document the negative result.

Multiple passes:

Simply rerunning the region grow does nothing - it is deterministic and takes no delay map as input, so a literal second pass is a verified no-op. Iteration only means something as feedback: smooth the current solution, then re-snap every voxel to whichever candidate lies nearest that smoothed field. That is iterated conditional modes on a smoothness regularised labelling problem, controlled by --numpasses (default 3).

passes

wrapped

non-alias jumps

raw

9591

27052

1 (region grow)

1685

24949

2

1463

24914

3

1272

24941

5

1184

25013

9

1129

25013

21

1000

25227

Returns diminish sharply after two or three passes, and the count of large non-periodic jumps starts creeping up - the signature of a smoothness prior beginning to invent structure rather than repair it. Three is a reasonable stopping point; running to convergence trades a little more wrap reduction against slowly accumulating damage elsewhere.

A caveat on all the despeckling comparisons here:

rapidtide has a --despeckle-patch-detection option that detects large connected patches of shifted delay and flags them for refitting - which is exactly the failure mode this program was written to address. It was not enabled in any of the runs compared above. A fair head to head needs that option turned on, and it has not been done.

On the choice of prior:

--prior smooth (the default) predicts that a voxel resembles its neighbors. --prior flow predicts using grad(tau) = v/|v|**2 from the corrflow optical flow field.

The flow prior was the original motivation, on the argument that a sidelobe of a travelling wave moves at the same velocity as the main lobe, so optical flow is immune to the ambiguity that defeats peak picking. That argument is correct, but it turns out not to matter much here: the sidelobe period is about 13.8 s while the median delay difference between adjacent voxels is 0.24 s, a ratio of nearly 60, so the ambiguity is enormous compared to the local gradient and almost any sensible local prior picks the right candidate. Measured head to head, flow gave 84-100 residual errors against 89-99 for smooth across several injection conditions - no reliable advantage, at the cost of computing the flow field.

Use --prior flow if your sidelobe spacing is small relative to your delay gradients, which is where it should begin to pay.

Inputs:

The similarity function, XXX_desc-corrout_info.nii.gz. The mask defaults to XXX_desc-corrfit_mask.nii.gz if present.

Outputs:

maxtimeunwrapped: the corrected delay map.

maxtimenaive: what naive peak picking would have given, for comparison.

unwrapchanged: mask of voxels assigned something other than their strongest peak.

unwrapconfidence: the ambiguity gap between the best and second best peak, which is what orders the region growing.

Usage:

Resolve sidelobe ambiguity in a rapidtide delay map by unwrapping it against the optical flow velocity field from corrflow.

usage: unwrapdelay [-h] [--maskfile FILE] [--prior {smooth,flow}]
                   [--maxcandidates N] [--numpasses N]
                   [--minconfidence QUANTILE] [--maxdeltatau SECONDS]
                   [--fitradius RADIUS] [--lagoversamp N] [--noprogressbar]
                   [--debug]
                   corrfile outputroot

Positional Arguments

corrfile

The rapidtide similarity function, XXX_desc-corrout_info.nii.gz.

outputroot

The root name of the output files.

Named Arguments

--maskfile

The analysis mask. Defaults to XXX_desc-corrfit_mask.nii.gz.

--prior

Possible choices: smooth, flow

Which local prior to predict each voxel’s delay with. “smooth” predicts that a voxel resembles its neighbors and needs no optical flow, so it is much faster. “flow” predicts using grad(tau) = v/|v|**2 from corrflow. The default is smooth because on the test data the flow prior was not measurably better - the sidelobe period dwarfs the local delay gradient, so the ambiguity is easy to resolve either way. Try flow if your sidelobe spacing is small relative to your delay gradients.

Default: 'smooth'

--maxcandidates

Consider at most this many candidate peaks per voxel, strongest first. Default is 6.

Default: 6

--numpasses

Number of passes. Pass 1 is the region grow; later passes re-snap every voxel to the candidate nearest a smoothed version of the current solution. Note that simply rerunning the region grow does nothing - it is deterministic and takes no delay map as input - so iteration only means anything in this feedback form. Returns diminish sharply and the smoothness prior slowly starts inventing structure, so more is not better. Default is 3.

Default: 3

--minconfidence

Confidence floor, as a quantile (0-1) of the per voxel peak ambiguity gap. Voxels below it are still assigned, but are not trusted to predict for their neighbors. Region growing otherwise propagates its own mistakes - a voxel assigned a wrapped value drags its correct neighbors along - which is where newly wrapped voxels come from. Default is 0.0 (off).

Default: 0.0

--maxdeltatau

Clip the predicted delay change between adjacent voxels to this, in seconds. This keeps a wild velocity estimate from dragging the unwrap off course. Default is 3.0.

Default: 3.0

--fitradius

Radius for the optical flow fit, in mm. Default is 6.0.

Default: 6.0

--lagoversamp

Interpolate the lag axis N times finer before estimating flow. Default is 1.

Default: 1

--noprogressbar

Disable progress bars.

Default: True

--debug

Turn on debugging information.

Default: False