Skip to content

feature: Add Weighted Least Squares (WLS) IVIM fitting algorithm — DT_IIITN (Feature #110)#148

Open
Devguru-codes wants to merge 1 commit intoOSIPI:mainfrom
Devguru-codes:feature-wls-ivim-fitting
Open

feature: Add Weighted Least Squares (WLS) IVIM fitting algorithm — DT_IIITN (Feature #110)#148
Devguru-codes wants to merge 1 commit intoOSIPI:mainfrom
Devguru-codes:feature-wls-ivim-fitting

Conversation

@Devguru-codes
Copy link

@Devguru-codes Devguru-codes commented Mar 1, 2026

Hello @oliverchampion. The man who was assigned this issue didnt update his PR and had issues in it so I worked according to "how to commit" given in readme. I have also taken 1 decision explained in this PR ( statsmodels.robust.robust_linear_model.RLM vs numpy WSL) . If you suggest that we should use the former than latter - please let me know and I will update this PR.

What this PR does

Closes #110. Implements a Weighted Least Squares (WLS) segmented IVIM fitting algorithm in accordance with the repository's contribution structure, following the approach recommended in Veraart et al. (2013).

This PR addresses the rejected PR #136, which implemented the same feature but placed code in the wrong location. This PR follows the correct structure:

  • Fit code → src/original/Initials_Institution/
  • Standardized wrapper → src/standardized/
  • Test registration → algorithms.json (no custom test file needed)

Algorithm

Segmented two-step WLS approach:

Step 1 — Estimate D (high b-values, b ≥ 200 s/mm²):
At high b-values the perfusion component decays to ~0, leaving:

ln(S(b)) ≈ ln(1 - f) - b·D

Weighted linear regression on log-signal with weights w = S(b)² (Veraart correction for heteroscedasticity introduced by the log-transform). Intercept yields f.

Step 2 — Estimate D* (low b-values):
Subtract the diffusion component to isolate the perfusion signal:

residual(b) = S(b) - (1-f)·exp(-b·D)
ln(residual(b)) ≈ ln(f) - b·D*

Weighted linear regression with weights w = residual² gives D*.

Implementation: Pure numpy (no extra dependencies). Uses the closed-form weighted normal equations:

β = (X^T W X)^{-1} X^T W y

Files Changed

Action File Description
NEW src/original/DT_IIITN/__init__.py Package marker
NEW src/original/DT_IIITN/wls_ivim_fitting.py Raw WLS algorithm (numpy only)
NEW src/standardized/DT_IIITN_WLS.py OsipiBase standardized wrapper
MODIFY tests/IVIMmodels/unit_tests/algorithms.json Registered DT_IIITN_WLS

API Usage

from src.wrappers.OsipiBase import OsipiBase
import numpy as np

bvals = np.array([0, 10, 20, 50, 100, 200, 400, 800])

# Basic usage
fit = OsipiBase(algorithm="DT_IIITN_WLS", bvalues=bvals)
results = fit.ivim_fit(signal, bvals)
# → {"D": 0.0010, "f": 0.08, "Dp": 0.013}

# Custom segmentation threshold (default: 200 s/mm²)
fit = OsipiBase(algorithm="DT_IIITN_WLS", bvalues=bvals, thresholds=[150])

# Direct class import
from src.standardized.DT_IIITN_WLS import DT_IIITN_WLS
fit = DT_IIITN_WLS(bvalues=bvals)

Design Decision: numpy WLS vs statsmodels RLM

The issue description suggests using statsmodels.robust.robust_linear_model.RLM. After evaluating both approaches:

statsmodels RLM This PR (numpy WLS)
Weighting Iterative reweighting (IRLS) via M-estimators Closed-form S² weighting (Veraart 2013)
Test suite runtime ~4+ hours (per-voxel IRLS is very slow) ~17 minutes
Extra dependency statsmodels None (numpy only)
Statistical basis Robust to outliers via IRLS Optimal WLS for Rician-noise log-signal

The numpy WLS approach implements the same weighting strategy recommended by Veraart et al. (w = S²) but using the closed-form normal equations — which is exactly what the paper validates. The statsmodels.RLM iterative approach adds robustness to outliers but at a ~15× runtime cost that is impractical for voxel-wise fitting on large datasets.

If you (anyone who will see my pr) prefers the full RLM approach, it can be easily switched in wls_ivim_fitting.py by replacing the _weighted_linreg call with sm.RLM(...).fit().


Reference

Veraart, J. et al. (2013). "Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls." NeuroImage, 81, 335–346.
DOI: 10.1016/j.neuroimage.2013.05.028


Testing

Testing is handled automatically by the existing framework via algorithms.json registration.

Metric Result
Total passed 1184 (57 new tests for DT_IIITN_WLS)
Regressions 0
test_ivim_fit_saved 593 passed, exit code 0 ✅
test_volume / test_parallel errors Pre-existing on main (Windows FileNotFoundError in subprocess, unrelated to this PR — confirmed by running same tests on main before this branch)

Sample fit result on synthetic signal (f=0.10, D=0.0010, D*=0.010):

Fitted: f=0.080, D=0.0010, Dp=0.013  ✅ (expected for linearized WLS)

Checklist

  • Self-review of changed code
  • Follows src/original/Initials_Institution/ + src/standardized/ contribution structure
  • No custom test file — registered in algorithms.json
  • No new dependencies (numpy only — statsmodels evaluated but removed for performance; see Design Decision above)
  • Tested against automated test suite — zero regressions

Implement Weighted Least Squares (WLS) segmented IVIM fitting following
Veraart et al. (2013) NeuroImage 81:335-346.

Algorithm:
- Step 1: Fit D from high b-values via WLS on log-signal (w=S^2)
- Step 2: Fit D* from residuals at low b-values via WLS

New files:
- src/original/DT_IIITN/wls_ivim_fitting.py (raw algorithm, numpy only)
- src/standardized/DT_IIITN_WLS.py (OsipiBase standardized wrapper)

Modified files:
- tests/IVIMmodels/unit_tests/algorithms.json (register for automated tests)

Follows repository contribution structure:
- Fit code in src/original/Initials_Institution/
- Standardized wrapper in src/standardized/
- Registered in algorithms.json (no custom tests needed)

Test results: 1184 passed, 167 skipped, 22 xfailed, 6 xpassed.
27 test_volume errors are pre-existing (FileNotFoundError on Windows).
required_bvalues = 4
required_thresholds = [0, 0]
required_bounds = False
required_bounds_optional = True
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@IvanARashid , I'm a bit confused what the difference is between
requiered bounds
requiered bounds optional
supported bounds

Does it make sense for the second to be true and the others false?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[Feature] Weighted Least Squares algorithm

2 participants