fix(dsc): restore circulant deconvolution for cSVD and oSVD

[MohamedNasser8]

Fixes #166

The bug

_build_circulant_matrix_xp built the full 2n x 2n block-circulant matrix and then returned only its upper-left n x n block:

idx = (arange(2n)[:, None] - arange(2n)[None, :]) % (2n)
A = aif_padded[idx]
return A[:n, :n]      # <- the problem

Because the padded half of the AIF is zero, that block is exactly the lower-triangular Toeplitz matrix that sSVD uses. Verified directly: circulant[:n, :n] == toeplitz is True, bit for bit.

So cSVD and oSVD were mathematically identical to sSVD. The wraparound terms in the upper-right of the matrix were thrown away before the SVD ever ran. That's why the reporter saw identical residue functions and CBF maps from all three methods.

The fix

Move the truncation from the matrix to the solution:

• return the matrix in full (2n x 2n)
• zero-pad the concentration curve to 2n before solving
• keep only the first n points of the recovered residue

The wraparound terms now in the SVD instead of being discarded before it.

The same padding is applied everywhere the matrix is used.

Second bug: oSVD threshold selection

While verifying the above I found the oSVD threshold search was inverted:

improved = (oi < target_oi) & (oi < best_oi)   # keeps the LOWEST OI

The oscillation index falls monotonically as truncation increases, so this never stopped, it ran to the top of the candidate range and over-smoothed the residue, flattening the peak and underestimating CBF.

Wu et al. intend the minimal regularization that suppresses ringing, so the search now stops at the first threshold meeting the target.

Verification

Measured on the OSIPI ground-truth cases with an injected bolus delay.

CBF stability under delay (mean change across all 14 cases, lower is better):

method	mean change
sSVD	4.14
cSVD	1.69 (2.4x better)
oSVD	0.99 (4.2x better)

The gap widens with flow, as expected, at CBF=70, sSVD's estimate moves by 12.3 mL/100mL/min from bolus timing alone, versus 2.5 for oSVD.

Delay recovery: cSVD recovers the injected delay exactly in 14/14 cases. sSVD is consistently off by one frame.

Note on the compliance suite

cSVD and oSVD now fail 4/14 cases on the delay-free compliance data and are marked xfail.

Worth flagging: all 14 compliance cases are delay0_dispersion0, so the suite only measures the regime where block-circulant is expected to be slightly worse, and never tests the delay-insensitivity these methods exist for. Adding delayed cases to the fixture would be a good follow-up.

The existing cSVD xfail reason was also stale, it said "13/14 cases pass", which described the pre-fix behaviour where cSVD was silently inheriting sSVD's results.

[julienpoublanc-uhn]

Thank you for fixing that so quickly. I will try it again.
Talk to you later.

[MohamedNasser8]

Thanks for the tests @julienpoublanc-uhn, that was a bug.

deconvolve() only read a keyword literally named params, so SSVDConfig /
CSVDConfig / OSVDConfig were silently ignored and the default 0.2 always
applied. Both forms work now:

d.deconvolve(..., threshold=0.1)
d.deconvolve(..., params=SSVDConfig(method="sSVD", threshold=0.1))

I also widened oSVD's threshold search from 0.01–0.5 to 0.01–0.95. On shorter
acquisitions the target OI wasn't reachable in the old range, so every voxel
fell back to 0.2, which is why the two looked identical.

Note: for oSVD use oscillation_index; default_threshold is only the fallback.

If there are more concerns about oSVD deconvolution we can look for possible refinements in another PR

[julienpoublanc-uhn]

Thank you for fixing that! It looks good.
I will test it and let you know.
Talk to you later.

[julienpoublanc-uhn]

Both forms of threshold setting you mentioned above are working now. Thanks!

[MohamedNasser8]

Hey @ltorres6 I think this PR is ready, whenever you have time.

[julienpoublanc-uhn]

Hi Mohamed, Luis,

A while ago, I developped a DSC analysis program with scripts and a graphical interface, including visual, automated AIF selection, deconvoltion ...

Anyway, that also why I wanted to compare my residue function with your. I am glad because I obtained exactly the same results for both sSVD and cSVD as OSIPI (I do not have an implementation of oSVD).

However, I noticed that OSIPI sets negative values of the residue function R(t) to zero. I was a little surprised by that choice. The tail of R(t) typically oscillates around zero, with both negative and positive values. If only the positive values are retained, it seems to me that this could introduce a bias.

Actually, I think MTT calculation might reveals the potential problem. Apart from using the central volume theorem, MTT can also be estimated from the area under the residue function. However, if negative values are set to zero, the estimated MTT can become heavily dependent on the chosen time window used for the bolus. For example, the MTT could continue to increase depending on whether a 30, 45, or 60-second window is selected.

On the other hand, if both the negative and positive oscillations are retained, they tend to cancel each other out as the time window increases. This seems more consistent with what we would expect from the true residue function, which should decay toward zero, like the ideal exponential residue function.

Any thought on this?

Thank you.

[ltorres6]

Hi Julien - Thanks for taking a look and for your valuable feedback.

I agree, we need to retain the negative values and this is a bug. Could you open a separate issue for this so we can track and resolve independently from this threshold wiring issue?