feat: add RLM method toggle to WLS IVIM fitting

Author: Devguru-codes

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,26 +49,52 @@ 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
             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)
 
@@ -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

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