From 04b1ac3ef427f768b09a2a0b9eb757ce476cfdea Mon Sep 17 00:00:00 2001 From: Devguru Date: Sun, 1 Mar 2026 16:40:25 +0530 Subject: [PATCH 1/3] feat: add WLS IVIM fitting algorithm (Feature #110) 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). --- src/original/DT_IIITN/__init__.py | 1 + src/original/DT_IIITN/wls_ivim_fitting.py | 139 ++++++++++++++++++++ src/standardized/DT_IIITN_WLS.py | 87 ++++++++++++ tests/IVIMmodels/unit_tests/algorithms.json | 18 ++- 4 files changed, 239 insertions(+), 6 deletions(-) create mode 100644 src/original/DT_IIITN/__init__.py create mode 100644 src/original/DT_IIITN/wls_ivim_fitting.py create mode 100644 src/standardized/DT_IIITN_WLS.py diff --git a/src/original/DT_IIITN/__init__.py b/src/original/DT_IIITN/__init__.py new file mode 100644 index 00000000..70146a93 --- /dev/null +++ b/src/original/DT_IIITN/__init__.py @@ -0,0 +1 @@ +# WLS IVIM fitting by Devguru Tiwari, IIIT Nagpur diff --git a/src/original/DT_IIITN/wls_ivim_fitting.py b/src/original/DT_IIITN/wls_ivim_fitting.py new file mode 100644 index 00000000..87cfce5d --- /dev/null +++ b/src/original/DT_IIITN/wls_ivim_fitting.py @@ -0,0 +1,139 @@ +""" +Weighted Least Squares (WLS) IVIM fitting. + +Author: Devguru Tiwari, IIIT Nagpur +Date: 2026-03-01 + +Implements a segmented approach for IVIM parameter estimation: +1. Estimate D from high b-values using weighted linear regression on log-signal +2. Estimate f from the intercept of the Step 1 fit +3. Estimate D* from residuals at low b-values using weighted linear regression + +Weighting follows Veraart et al. (2013): weights = signal^2 to account +for heteroscedasticity introduced by the log-transform. + +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 + +Requirements: + numpy +""" + +import numpy as np +import warnings + + +def _weighted_linreg(x, y, weights): + """Fast weighted linear regression: y = a + b*x. + + Args: + x: 1D array, independent variable. + y: 1D array, dependent variable. + weights: 1D array, weights for each observation. + + Returns: + (intercept, slope) tuple. + """ + W = np.diag(weights) + X = np.column_stack([np.ones_like(x), x]) + # Weighted normal equations: (X^T W X) beta = X^T W y + XtW = X.T @ W + beta = np.linalg.solve(XtW @ X, XtW @ y) + return beta[0], beta[1] # intercept, slope + + +def wls_ivim_fit(bvalues, signal, cutoff=200): + """ + Weighted Least Squares IVIM fit (segmented approach). + + Step 1: Fit D from high b-values using WLS on log-signal. + Weights = S(b)^2 (Veraart et al. 2013). + Step 2: Fit D* from residuals at low b-values using WLS. + + Args: + bvalues (array-like): 1D array of b-values (s/mm²). + signal (array-like): 1D array of signal intensities (will be normalized). + cutoff (float): b-value threshold separating D from D* fitting. + Default: 200 s/mm². + + Returns: + tuple: (D, f, Dp) where + D (float): True diffusion coefficient (mm²/s). + f (float): Perfusion fraction (0-1). + Dp (float): Pseudo-diffusion coefficient (mm²/s). + """ + bvalues = np.array(bvalues, dtype=float) + signal = np.array(signal, dtype=float) + + # Normalize signal to S(b=0) + s0_vals = signal[bvalues == 0] + if len(s0_vals) == 0 or np.mean(s0_vals) <= 0: + return 0.0, 0.0, 0.0 + s0 = np.mean(s0_vals) + signal = signal / s0 + + try: + # ── Step 1: Estimate D from high b-values ───────────────────── + # At high b, perfusion component ≈ 0, so: + # S(b) ≈ (1 - f) * exp(-b * D) + # ln(S(b)) = ln(1 - f) - b * D + # Weighted linear fit: weights = S(b)^2 (Veraart correction) + + high_mask = bvalues >= cutoff + b_high = bvalues[high_mask] + s_high = signal[high_mask] + + # Guard against zero/negative signal values + s_high = np.maximum(s_high, 1e-8) + log_s = np.log(s_high) + + # Veraart weights: w = S^2 (corrects for noise amplification in log-domain) + weights_high = s_high ** 2 + + # WLS: ln(S) = intercept + slope * (-b) ⟹ slope = D + intercept, D = _weighted_linreg(-b_high, log_s, weights_high) + + # Extract f from intercept: intercept = ln(1 - f) + f = 1.0 - np.exp(intercept) + + # Clamp to physically meaningful ranges + D = np.clip(D, 0, 0.005) + f = np.clip(f, 0, 1) + + # ── Step 2: Estimate D* from low b-value residuals ──────────── + # Subtract the diffusion component: + # residual(b) = S(b) - (1 - f) * exp(-b * D) + # ≈ f * exp(-b * D*) + # ln(residual) = ln(f) - b * D* + + residual = signal - (1 - f) * np.exp(-bvalues * D) + + low_mask = (bvalues < cutoff) & (bvalues > 0) + b_low = bvalues[low_mask] + r_low = residual[low_mask] + + # Guard against zero/negative residuals + r_low = np.maximum(r_low, 1e-8) + log_r = np.log(r_low) + + weights_low = r_low ** 2 + + if len(b_low) >= 2: + _, Dp = _weighted_linreg(-b_low, log_r, weights_low) + Dp = np.clip(Dp, 0.005, 0.2) + else: + Dp = 0.01 # fallback + + # Ensure D* > D (by convention) + if Dp < D: + D, Dp = Dp, D + f = 1 - f + + return D, f, Dp + + except Exception: + # If fit fails, return zeros (consistent with other algorithms) + return 0.0, 0.0, 0.0 diff --git a/src/standardized/DT_IIITN_WLS.py b/src/standardized/DT_IIITN_WLS.py new file mode 100644 index 00000000..b89858fb --- /dev/null +++ b/src/standardized/DT_IIITN_WLS.py @@ -0,0 +1,87 @@ +from src.wrappers.OsipiBase import OsipiBase +from src.original.DT_IIITN.wls_ivim_fitting import wls_ivim_fit +import numpy as np + + +class DT_IIITN_WLS(OsipiBase): + """ + Weighted Least Squares IVIM fitting using statsmodels Robust Linear Model. + + Segmented approach: + 1. Estimate D from high b-values using robust linear regression on log-signal + 2. Estimate D* from residuals at low b-values using robust linear regression + + Author: Devguru Tiwari, IIIT Nagpur + + 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 + """ + + # Algorithm identification + id_author = "Devguru Tiwari, IIIT Nagpur" + id_algorithm_type = "Weighted least squares segmented fit" + id_return_parameters = "f, D*, D" + id_units = "seconds per milli metre squared or milliseconds per micro metre squared" + id_ref = "https://doi.org/10.1016/j.neuroimage.2013.05.028" + + # Algorithm requirements + required_bvalues = 4 + required_thresholds = [0, 0] + required_bounds = False + required_bounds_optional = True + required_initial_guess = False + required_initial_guess_optional = True + + # Supported inputs + supported_bounds = False + supported_initial_guess = False + supported_thresholds = True + supported_dimensions = 1 + supported_priors = False + + def __init__(self, bvalues=None, thresholds=None, + bounds=None, initial_guess=None): + """ + Initialize the WLS IVIM fitting algorithm. + + Args: + bvalues (array-like, optional): b-values for the fitted signals. + thresholds (array-like, optional): Threshold b-value for segmented + fitting. The first value is used as the cutoff between high + and low b-values. Default: 200 s/mm². + bounds (dict, optional): Not used by this algorithm. + initial_guess (dict, optional): Not used by this algorithm. + """ + super(DT_IIITN_WLS, self).__init__( + bvalues=bvalues, bounds=bounds, + initial_guess=initial_guess, thresholds=thresholds + ) + + def ivim_fit(self, signals, bvalues, **kwargs): + """Perform the IVIM fit using WLS. + + Args: + signals (array-like): Signal intensities at each b-value. + bvalues (array-like, optional): b-values for the signals. + + Returns: + dict: Dictionary with keys "D", "f", "Dp". + """ + bvalues = np.array(bvalues) + + # Use threshold as cutoff if available + cutoff = 200 + if self.thresholds is not None and len(self.thresholds) > 0: + cutoff = self.thresholds[0] + + D, f, Dp = wls_ivim_fit(bvalues, signals, cutoff=cutoff) + + results = {} + results["D"] = D + results["f"] = f + results["Dp"] = Dp + + return results diff --git a/tests/IVIMmodels/unit_tests/algorithms.json b/tests/IVIMmodels/unit_tests/algorithms.json index 7da8258e..83ad20e1 100644 --- a/tests/IVIMmodels/unit_tests/algorithms.json +++ b/tests/IVIMmodels/unit_tests/algorithms.json @@ -25,14 +25,18 @@ "OJ_GU_seg", "OJ_GU_segMATLAB", "OJ_GU_bayesMATLAB", - "TF_reference_IVIMfit" + "TF_reference_IVIMfit", + "DT_IIITN_WLS" ], "TCML_TechnionIIT_lsqBOBYQA": { - "xfail_names": {"pericardium": false, "bounds": false} + "xfail_names": { + "pericardium": false, + "bounds": false + } }, - "IVIM_NEToptim": { + "IVIM_NEToptim": { "deep_learning": true, - "n":3000000 + "n": 3000000 }, "Super_IVIM_DC": { "deep_learning": true @@ -51,7 +55,9 @@ }, "ETP_SRI_LinearFitting": { "options": { - "thresholds": [500] + "thresholds": [ + 500 + ] }, "tolerances": { "rtol": { @@ -81,4 +87,4 @@ }, "test_bounds": false } -} +} \ No newline at end of file From 1fbd37e2428fb6299ba318108d8299b60a85e520 Mon Sep 17 00:00:00 2001 From: Devguru Date: Thu, 12 Mar 2026 16:34:36 +0530 Subject: [PATCH 2/3] feat: add RLM method toggle to WLS IVIM fitting --- src/original/DT_IIITN/wls_ivim_fitting.py | 73 ++++++++++++++++------- src/standardized/DT_IIITN_WLS.py | 23 ++++--- 2 files changed, 67 insertions(+), 29 deletions(-) diff --git a/src/original/DT_IIITN/wls_ivim_fitting.py b/src/original/DT_IIITN/wls_ivim_fitting.py index 87cfce5d..d30e88ba 100644 --- a/src/original/DT_IIITN/wls_ivim_fitting.py +++ b/src/original/DT_IIITN/wls_ivim_fitting.py @@ -1,16 +1,17 @@ """ -Weighted Least Squares (WLS) IVIM fitting. +Weighted Least Squares (WLS) / Robust Linear Model (RLM) IVIM fitting. Author: Devguru Tiwari, IIIT Nagpur Date: 2026-03-01 Implements a segmented approach for IVIM parameter estimation: -1. Estimate D from high b-values using weighted linear regression on log-signal +1. Estimate D from high b-values using weighted/robust linear regression on log-signal 2. Estimate f from the intercept of the Step 1 fit -3. Estimate D* from residuals at low b-values using weighted linear regression +3. Estimate D* from residuals at low b-values using weighted/robust linear regression -Weighting follows Veraart et al. (2013): weights = signal^2 to account -for heteroscedasticity introduced by the log-transform. +Two regression methods are available: +- WLS: Weighted Linear Least Squares with Veraart weights (w = S^2) +- RLM: Robust Linear Model using Huber's T norm (statsmodels) Reference: Veraart, J. et al. (2013). "Weighted linear least squares estimation of @@ -20,6 +21,7 @@ Requirements: numpy + statsmodels (only for method="RLM") """ import numpy as np @@ -29,6 +31,8 @@ def _weighted_linreg(x, y, weights): """Fast weighted linear regression: y = a + b*x. + Uses Veraart et al. (2013) approach with weights = S^2. + Args: x: 1D array, independent variable. y: 1D array, dependent variable. @@ -45,19 +49,41 @@ def _weighted_linreg(x, y, weights): return beta[0], beta[1] # intercept, slope -def wls_ivim_fit(bvalues, signal, cutoff=200): +def _rlm_linreg(x, y): + """Robust linear regression using statsmodels RLM with Huber's T norm. + + RLM down-weights outlier observations via iteratively reweighted least + squares (IRLS), making the fit resistant to corrupted/noisy voxels. + + Args: + x: 1D array, independent variable. + y: 1D array, dependent variable. + + Returns: + (intercept, slope) tuple. + """ + import statsmodels.api as sm + X = sm.add_constant(x) + model = sm.RLM(y, X, M=sm.robust.norms.HuberT()) + result = model.fit() + return result.params[0], result.params[1] # intercept, slope + + +def wls_ivim_fit(bvalues, signal, cutoff=200, method="WLS"): """ - Weighted Least Squares IVIM fit (segmented approach). + IVIM fit using WLS or RLM (segmented approach). - Step 1: Fit D from high b-values using WLS on log-signal. - Weights = S(b)^2 (Veraart et al. 2013). - Step 2: Fit D* from residuals at low b-values using WLS. + Step 1: Fit D from high b-values on log-signal. + Step 2: Fit D* from residuals at low b-values. Args: bvalues (array-like): 1D array of b-values (s/mm²). signal (array-like): 1D array of signal intensities (will be normalized). cutoff (float): b-value threshold separating D from D* fitting. Default: 200 s/mm². + method (str): Regression method to use. + - "WLS": Weighted Least Squares with Veraart S² weights (default). + - "RLM": Robust Linear Model with Huber's T norm (statsmodels). Returns: tuple: (D, f, Dp) where @@ -65,6 +91,10 @@ def wls_ivim_fit(bvalues, signal, cutoff=200): f (float): Perfusion fraction (0-1). Dp (float): Pseudo-diffusion coefficient (mm²/s). """ + method = method.upper() + if method not in ("WLS", "RLM"): + raise ValueError(f"Unknown method '{method}'. Use 'WLS' or 'RLM'.") + bvalues = np.array(bvalues, dtype=float) signal = np.array(signal, dtype=float) @@ -80,8 +110,6 @@ def wls_ivim_fit(bvalues, signal, cutoff=200): # At high b, perfusion component ≈ 0, so: # S(b) ≈ (1 - f) * exp(-b * D) # ln(S(b)) = ln(1 - f) - b * D - # Weighted linear fit: weights = S(b)^2 (Veraart correction) - high_mask = bvalues >= cutoff b_high = bvalues[high_mask] s_high = signal[high_mask] @@ -90,11 +118,13 @@ def wls_ivim_fit(bvalues, signal, cutoff=200): s_high = np.maximum(s_high, 1e-8) log_s = np.log(s_high) - # Veraart weights: w = S^2 (corrects for noise amplification in log-domain) - weights_high = s_high ** 2 - - # WLS: ln(S) = intercept + slope * (-b) ⟹ slope = D - intercept, D = _weighted_linreg(-b_high, log_s, weights_high) + if method == "WLS": + # Veraart weights: w = S^2 (corrects for noise in log-domain) + weights_high = s_high ** 2 + intercept, D = _weighted_linreg(-b_high, log_s, weights_high) + else: + # RLM: robust regression, no explicit weights needed + intercept, D = _rlm_linreg(-b_high, log_s) # Extract f from intercept: intercept = ln(1 - f) f = 1.0 - np.exp(intercept) @@ -108,7 +138,6 @@ def wls_ivim_fit(bvalues, signal, cutoff=200): # residual(b) = S(b) - (1 - f) * exp(-b * D) # ≈ f * exp(-b * D*) # ln(residual) = ln(f) - b * D* - residual = signal - (1 - f) * np.exp(-bvalues * D) low_mask = (bvalues < cutoff) & (bvalues > 0) @@ -119,10 +148,12 @@ def wls_ivim_fit(bvalues, signal, cutoff=200): r_low = np.maximum(r_low, 1e-8) log_r = np.log(r_low) - weights_low = r_low ** 2 - if len(b_low) >= 2: - _, Dp = _weighted_linreg(-b_low, log_r, weights_low) + if method == "WLS": + weights_low = r_low ** 2 + _, Dp = _weighted_linreg(-b_low, log_r, weights_low) + else: + _, Dp = _rlm_linreg(-b_low, log_r) Dp = np.clip(Dp, 0.005, 0.2) else: Dp = 0.01 # fallback diff --git a/src/standardized/DT_IIITN_WLS.py b/src/standardized/DT_IIITN_WLS.py index b89858fb..a79c8a95 100644 --- a/src/standardized/DT_IIITN_WLS.py +++ b/src/standardized/DT_IIITN_WLS.py @@ -5,11 +5,15 @@ class DT_IIITN_WLS(OsipiBase): """ - Weighted Least Squares IVIM fitting using statsmodels Robust Linear Model. + Segmented IVIM fitting with selectable regression method. + + Two methods are available: + - WLS: Weighted Least Squares with Veraart S² weights (default) + - RLM: Robust Linear Model with Huber's T norm (statsmodels) Segmented approach: - 1. Estimate D from high b-values using robust linear regression on log-signal - 2. Estimate D* from residuals at low b-values using robust linear regression + 1. Estimate D from high b-values using linear regression on log-signal + 2. Estimate D* from residuals at low b-values using linear regression Author: Devguru Tiwari, IIIT Nagpur @@ -22,7 +26,7 @@ class DT_IIITN_WLS(OsipiBase): # Algorithm identification id_author = "Devguru Tiwari, IIIT Nagpur" - id_algorithm_type = "Weighted least squares segmented fit" + id_algorithm_type = "Weighted least squares / robust linear model segmented fit" id_return_parameters = "f, D*, D" id_units = "seconds per milli metre squared or milliseconds per micro metre squared" id_ref = "https://doi.org/10.1016/j.neuroimage.2013.05.028" @@ -43,9 +47,9 @@ class DT_IIITN_WLS(OsipiBase): supported_priors = False def __init__(self, bvalues=None, thresholds=None, - bounds=None, initial_guess=None): + bounds=None, initial_guess=None, method="WLS"): """ - Initialize the WLS IVIM fitting algorithm. + Initialize the IVIM fitting algorithm. Args: bvalues (array-like, optional): b-values for the fitted signals. @@ -54,14 +58,16 @@ def __init__(self, bvalues=None, thresholds=None, and low b-values. Default: 200 s/mm². bounds (dict, optional): Not used by this algorithm. initial_guess (dict, optional): Not used by this algorithm. + method (str): Regression method — "WLS" (default) or "RLM". """ super(DT_IIITN_WLS, self).__init__( bvalues=bvalues, bounds=bounds, initial_guess=initial_guess, thresholds=thresholds ) + self.method = method.upper() def ivim_fit(self, signals, bvalues, **kwargs): - """Perform the IVIM fit using WLS. + """Perform the IVIM fit using the selected method (WLS or RLM). Args: signals (array-like): Signal intensities at each b-value. @@ -77,7 +83,8 @@ def ivim_fit(self, signals, bvalues, **kwargs): if self.thresholds is not None and len(self.thresholds) > 0: cutoff = self.thresholds[0] - D, f, Dp = wls_ivim_fit(bvalues, signals, cutoff=cutoff) + D, f, Dp = wls_ivim_fit(bvalues, signals, cutoff=cutoff, + method=self.method) results = {} results["D"] = D From edd53526bf172a0f1e147ae46b708f2292ec4195 Mon Sep 17 00:00:00 2001 From: Devguru Date: Sat, 14 Mar 2026 01:40:44 +0530 Subject: [PATCH 3/3] set required_bounds_optional and required_initial_guess_optional to False per review --- src/standardized/DT_IIITN_WLS.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/standardized/DT_IIITN_WLS.py b/src/standardized/DT_IIITN_WLS.py index a79c8a95..4925dbed 100644 --- a/src/standardized/DT_IIITN_WLS.py +++ b/src/standardized/DT_IIITN_WLS.py @@ -35,9 +35,9 @@ class DT_IIITN_WLS(OsipiBase): required_bvalues = 4 required_thresholds = [0, 0] required_bounds = False - required_bounds_optional = True + required_bounds_optional = False required_initial_guess = False - required_initial_guess_optional = True + required_initial_guess_optional = False # Supported inputs supported_bounds = False