Merge pull request #148 from Devguru-codes/feature-wls-ivim-fitting

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

Author: oliverchampion

src/original/DT_IIITN/__init__.py

@@ -0,0 +1 @@
+# WLS IVIM fitting by Devguru Tiwari, IIIT Nagpur

src/original/DT_IIITN/wls_ivim_fitting.py

@@ -0,0 +1,170 @@
+"""
+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/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/robust linear regression
+
+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
+    diffusion MRI parameters: strengths, limitations, and pitfalls."
+    NeuroImage, 81, 335-346.
+    DOI: 10.1016/j.neuroimage.2013.05.028
+
+Requirements:
+    numpy
+    statsmodels (only for method="RLM")
+"""
+
+import numpy as np
+import warnings
+
+
+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.
+        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 _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"):
+    """
+    IVIM fit using WLS or RLM (segmented approach).
+
+    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
+            D (float): True diffusion coefficient (mm²/s).
+            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)
+
+    # 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
+        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)
+
+        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)
+
+        # 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)
+
+        if len(b_low) >= 2:
+            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
+
+        # 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

src/standardized/DT_IIITN_WLS.py

@@ -0,0 +1,94 @@
+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):
+    """
+    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 linear regression on log-signal
+    2. Estimate D* from residuals at low b-values using 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 / 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"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0, 0]
+    required_bounds = False
+    required_bounds_optional = False
+    required_initial_guess = False
+    required_initial_guess_optional = False
+
+    # 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, method="WLS"):
+        """
+        Initialize the 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.
+            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 the selected method (WLS or RLM).
+
+        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,
+                                method=self.method)
+
+        results = {}
+        results["D"] = D
+        results["f"] = f
+        results["Dp"] = Dp
+
+        return results

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
     }
-}
+}