|
| 1 | +""" |
| 2 | +Weighted Least Squares (WLS) / Robust Linear Model (RLM) IVIM fitting. |
| 3 | +
|
| 4 | +Author: Devguru Tiwari, IIIT Nagpur |
| 5 | +Date: 2026-03-01 |
| 6 | +
|
| 7 | +Implements a segmented approach for IVIM parameter estimation: |
| 8 | +1. Estimate D from high b-values using weighted/robust linear regression on log-signal |
| 9 | +2. Estimate f from the intercept of the Step 1 fit |
| 10 | +3. Estimate D* from residuals at low b-values using weighted/robust linear regression |
| 11 | +
|
| 12 | +Two regression methods are available: |
| 13 | +- WLS: Weighted Linear Least Squares with Veraart weights (w = S^2) |
| 14 | +- RLM: Robust Linear Model using Huber's T norm (statsmodels) |
| 15 | +
|
| 16 | +Reference: |
| 17 | + Veraart, J. et al. (2013). "Weighted linear least squares estimation of |
| 18 | + diffusion MRI parameters: strengths, limitations, and pitfalls." |
| 19 | + NeuroImage, 81, 335-346. |
| 20 | + DOI: 10.1016/j.neuroimage.2013.05.028 |
| 21 | +
|
| 22 | +Requirements: |
| 23 | + numpy |
| 24 | + statsmodels (only for method="RLM") |
| 25 | +""" |
| 26 | + |
| 27 | +import numpy as np |
| 28 | +import warnings |
| 29 | + |
| 30 | + |
| 31 | +def _weighted_linreg(x, y, weights): |
| 32 | + """Fast weighted linear regression: y = a + b*x. |
| 33 | +
|
| 34 | + Uses Veraart et al. (2013) approach with weights = S^2. |
| 35 | +
|
| 36 | + Args: |
| 37 | + x: 1D array, independent variable. |
| 38 | + y: 1D array, dependent variable. |
| 39 | + weights: 1D array, weights for each observation. |
| 40 | +
|
| 41 | + Returns: |
| 42 | + (intercept, slope) tuple. |
| 43 | + """ |
| 44 | + W = np.diag(weights) |
| 45 | + X = np.column_stack([np.ones_like(x), x]) |
| 46 | + # Weighted normal equations: (X^T W X) beta = X^T W y |
| 47 | + XtW = X.T @ W |
| 48 | + beta = np.linalg.solve(XtW @ X, XtW @ y) |
| 49 | + return beta[0], beta[1] # intercept, slope |
| 50 | + |
| 51 | + |
| 52 | +def _rlm_linreg(x, y): |
| 53 | + """Robust linear regression using statsmodels RLM with Huber's T norm. |
| 54 | +
|
| 55 | + RLM down-weights outlier observations via iteratively reweighted least |
| 56 | + squares (IRLS), making the fit resistant to corrupted/noisy voxels. |
| 57 | +
|
| 58 | + Args: |
| 59 | + x: 1D array, independent variable. |
| 60 | + y: 1D array, dependent variable. |
| 61 | +
|
| 62 | + Returns: |
| 63 | + (intercept, slope) tuple. |
| 64 | + """ |
| 65 | + import statsmodels.api as sm |
| 66 | + X = sm.add_constant(x) |
| 67 | + model = sm.RLM(y, X, M=sm.robust.norms.HuberT()) |
| 68 | + result = model.fit() |
| 69 | + return result.params[0], result.params[1] # intercept, slope |
| 70 | + |
| 71 | + |
| 72 | +def wls_ivim_fit(bvalues, signal, cutoff=200, method="WLS"): |
| 73 | + """ |
| 74 | + IVIM fit using WLS or RLM (segmented approach). |
| 75 | +
|
| 76 | + Step 1: Fit D from high b-values on log-signal. |
| 77 | + Step 2: Fit D* from residuals at low b-values. |
| 78 | +
|
| 79 | + Args: |
| 80 | + bvalues (array-like): 1D array of b-values (s/mm²). |
| 81 | + signal (array-like): 1D array of signal intensities (will be normalized). |
| 82 | + cutoff (float): b-value threshold separating D from D* fitting. |
| 83 | + Default: 200 s/mm². |
| 84 | + method (str): Regression method to use. |
| 85 | + - "WLS": Weighted Least Squares with Veraart S² weights (default). |
| 86 | + - "RLM": Robust Linear Model with Huber's T norm (statsmodels). |
| 87 | +
|
| 88 | + Returns: |
| 89 | + tuple: (D, f, Dp) where |
| 90 | + D (float): True diffusion coefficient (mm²/s). |
| 91 | + f (float): Perfusion fraction (0-1). |
| 92 | + Dp (float): Pseudo-diffusion coefficient (mm²/s). |
| 93 | + """ |
| 94 | + method = method.upper() |
| 95 | + if method not in ("WLS", "RLM"): |
| 96 | + raise ValueError(f"Unknown method '{method}'. Use 'WLS' or 'RLM'.") |
| 97 | + |
| 98 | + bvalues = np.array(bvalues, dtype=float) |
| 99 | + signal = np.array(signal, dtype=float) |
| 100 | + |
| 101 | + # Normalize signal to S(b=0) |
| 102 | + s0_vals = signal[bvalues == 0] |
| 103 | + if len(s0_vals) == 0 or np.mean(s0_vals) <= 0: |
| 104 | + return 0.0, 0.0, 0.0 |
| 105 | + s0 = np.mean(s0_vals) |
| 106 | + signal = signal / s0 |
| 107 | + |
| 108 | + try: |
| 109 | + # ── Step 1: Estimate D from high b-values ───────────────────── |
| 110 | + # At high b, perfusion component ≈ 0, so: |
| 111 | + # S(b) ≈ (1 - f) * exp(-b * D) |
| 112 | + # ln(S(b)) = ln(1 - f) - b * D |
| 113 | + high_mask = bvalues >= cutoff |
| 114 | + b_high = bvalues[high_mask] |
| 115 | + s_high = signal[high_mask] |
| 116 | + |
| 117 | + # Guard against zero/negative signal values |
| 118 | + s_high = np.maximum(s_high, 1e-8) |
| 119 | + log_s = np.log(s_high) |
| 120 | + |
| 121 | + if method == "WLS": |
| 122 | + # Veraart weights: w = S^2 (corrects for noise in log-domain) |
| 123 | + weights_high = s_high ** 2 |
| 124 | + intercept, D = _weighted_linreg(-b_high, log_s, weights_high) |
| 125 | + else: |
| 126 | + # RLM: robust regression, no explicit weights needed |
| 127 | + intercept, D = _rlm_linreg(-b_high, log_s) |
| 128 | + |
| 129 | + # Extract f from intercept: intercept = ln(1 - f) |
| 130 | + f = 1.0 - np.exp(intercept) |
| 131 | + |
| 132 | + # Clamp to physically meaningful ranges |
| 133 | + D = np.clip(D, 0, 0.005) |
| 134 | + f = np.clip(f, 0, 1) |
| 135 | + |
| 136 | + # ── Step 2: Estimate D* from low b-value residuals ──────────── |
| 137 | + # Subtract the diffusion component: |
| 138 | + # residual(b) = S(b) - (1 - f) * exp(-b * D) |
| 139 | + # ≈ f * exp(-b * D*) |
| 140 | + # ln(residual) = ln(f) - b * D* |
| 141 | + residual = signal - (1 - f) * np.exp(-bvalues * D) |
| 142 | + |
| 143 | + low_mask = (bvalues < cutoff) & (bvalues > 0) |
| 144 | + b_low = bvalues[low_mask] |
| 145 | + r_low = residual[low_mask] |
| 146 | + |
| 147 | + # Guard against zero/negative residuals |
| 148 | + r_low = np.maximum(r_low, 1e-8) |
| 149 | + log_r = np.log(r_low) |
| 150 | + |
| 151 | + if len(b_low) >= 2: |
| 152 | + if method == "WLS": |
| 153 | + weights_low = r_low ** 2 |
| 154 | + _, Dp = _weighted_linreg(-b_low, log_r, weights_low) |
| 155 | + else: |
| 156 | + _, Dp = _rlm_linreg(-b_low, log_r) |
| 157 | + Dp = np.clip(Dp, 0.005, 0.2) |
| 158 | + else: |
| 159 | + Dp = 0.01 # fallback |
| 160 | + |
| 161 | + # Ensure D* > D (by convention) |
| 162 | + if Dp < D: |
| 163 | + D, Dp = Dp, D |
| 164 | + f = 1 - f |
| 165 | + |
| 166 | + return D, f, Dp |
| 167 | + |
| 168 | + except Exception: |
| 169 | + # If fit fails, return zeros (consistent with other algorithms) |
| 170 | + return 0.0, 0.0, 0.0 |
0 commit comments