diff --git a/cpp/src/common/maths.cpp b/cpp/src/common/maths.cpp index 21d0b00..438bde8 100644 --- a/cpp/src/common/maths.cpp +++ b/cpp/src/common/maths.cpp @@ -1,9 +1,21 @@ #include "common/maths.hpp" #include +#include namespace larmor { +namespace { + +// y below which exp(-y) I0(y) sums its power series, and above which its asymptotic series. +constexpr double BESSEL_SERIES_LIMIT = 25.0; +// Relative size of the last term at which a series stops. +constexpr double SERIES_TOLERANCE = 1e-17; +// Maximum number of terms of a series. +constexpr int MAX_SERIES_TERMS = 200; + +} // namespace + double norm3(const Vector3& vector) { return std::sqrt(vector.x() * vector.x() + vector.y() * vector.y() + vector.z() * vector.z()); } @@ -29,4 +41,26 @@ Vector3 rotate_to(const Vector3& reference, const Vector3& direction) { }; } +double scaled_bessel_i0(double y) { + if (y < BESSEL_SERIES_LIMIT) { + const double quarter_squared = 0.25 * y * y; + double term = 1.0; + double sum = 1.0; + for (int order = 1; order < MAX_SERIES_TERMS && term > SERIES_TOLERANCE * sum; ++order) { + term *= quarter_squared / (static_cast(order) * order); + sum += term; + } + return std::exp(-y) * sum; + } + // 1 / sqrt(2 pi y) sum_k c_k / y^k with c_k = c_{k-1} (2k - 1)^2 / (8 k); the terms fall + // until k is about 2 y, far beyond the precision needed. + double term = 1.0; + double sum = 1.0; + for (int order = 1; order < MAX_SERIES_TERMS && term > SERIES_TOLERANCE * sum; ++order) { + term *= (2.0 * order - 1.0) * (2.0 * order - 1.0) / (8.0 * order * y); + sum += term; + } + return sum / std::sqrt(2.0 * std::numbers::pi * y); +} + } // namespace larmor diff --git a/cpp/src/common/maths.hpp b/cpp/src/common/maths.hpp index f733537..0d9cee9 100644 --- a/cpp/src/common/maths.hpp +++ b/cpp/src/common/maths.hpp @@ -8,6 +8,9 @@ namespace larmor { /// Returns the Euclidean norm of a 3-vector. double norm3(const Vector3& vector); +/// Returns exp(-y) I0(y) for y >= 0, with I0 the modified Bessel function of the first kind. +double scaled_bessel_i0(double y); + /// Rotates `direction` by the rotation that takes +z to `reference`, a unit vector. /// /// When `reference` is -z, `direction` is negated. diff --git a/cpp/src/dwi/CMakeLists.txt b/cpp/src/dwi/CMakeLists.txt index 6490329..f9711af 100644 --- a/cpp/src/dwi/CMakeLists.txt +++ b/cpp/src/dwi/CMakeLists.txt @@ -1,6 +1,8 @@ -# Diffusion-weighted images: Gibbs-ringing removal, MP-PCA denoising and the log-signal fits. +# Diffusion-weighted images: Gibbs-ringing removal, MP-PCA denoising, the log-signal fits and +# the NODDI fit. # Lists this folder's sources, relative to it; the folder's tests go in tests/dwi/CMakeLists.txt. target_sources(larmor_numerics PRIVATE gibbs.cpp log_fit.cpp - mppca.cpp) + mppca.cpp + noddi.cpp) diff --git a/cpp/src/dwi/noddi.cpp b/cpp/src/dwi/noddi.cpp new file mode 100644 index 0000000..fbc3751 --- /dev/null +++ b/cpp/src/dwi/noddi.cpp @@ -0,0 +1,767 @@ +// The model, noise level and fit follow the NODDI Matlab Toolbox 1.05 (Zhang et al. 2012); the +// start follows AMICO (Daducci et al. 2015). +#include "dwi/noddi.hpp" + +#include +#include +#include +#include +#include +#include + +#include "common/maths.hpp" +#include "sd/nnls.hpp" + +namespace larmor::dwi { + +namespace { + +constexpr double PI = std::numbers::pi; +constexpr double SQRT_PI = 1.0 / std::numbers::inv_sqrtpi; + +// x from which erfi(x) is taken as exp(x^2) / (x sqrt(pi)). +constexpr double ERFI_ASYMPTOTIC = 5.7; +// Relative size of the last term at which the erfi series stops. +constexpr double SERIES_TOLERANCE = 1e-17; +// Maximum number of terms of the erfi series. +constexpr int MAX_SERIES_TERMS = 400; +// kappa above which the Watson coefficients come from fits in ln(kappa / 30). +constexpr double WATSON_FIT_KAPPA = 30.0; +// kappa up to which the Watson coefficients come from their series. +constexpr double WATSON_SERIES_KAPPA = 0.1; +// x = b d up to which the Legendre-Gaussian integrals come from their series. +constexpr double INTEGRAL_SERIES_X = 0.05; +// kappa below which the Watson-averaged diffusivities come from their series. +constexpr double HINDERED_SERIES_KAPPA = 1e-5; + +// The noise level sigma is the spread of the b0s, at least MIN_SNR times their mean, divided by +// SIGMA_SCALING. +constexpr double MIN_SNR = 0.02; +constexpr double SIGMA_SCALING = 100.0; +// The largest kappa of the fit. +constexpr double KAPPA_MAX = 64.0; +// The fit's parameter is kappa divided by this. +constexpr double KAPPA_SCALE = 10.0; +// Bounds of the encoded fractions and angles. +constexpr double FRACTION_UPPER = 1e5; +constexpr double ANGLE_LIMIT = 100.0; +// Argument above which log I0(x) is taken as x - log(2 pi x) / 2. +constexpr double LOG_BESSEL_ASYMPTOTIC = 700.0; +// Most iterations of each stage. +constexpr int MAX_ITERATIONS = 1000; +// Largest change of an encoded parameter at which a stage stops. +constexpr double STEP_TOLERANCE = 1e-6; +// Relative change of the objective below which a stage stops. +constexpr double OBJECTIVE_TOLERANCE = 1e-12; +// Line search: sufficient-decrease constant, and the number of halvings. +constexpr double ARMIJO = 1e-4; +constexpr int MAX_HALVINGS = 40; +// Largest change of an encoded parameter in the first step of a stage. +constexpr double FIRST_STEP = 0.1; +// How far inside [0, 1] the fit starts the fractions. +constexpr double START_MARGIN = 0.01; +// Fewest positive diffusion-weighted measurements the fit needs. +constexpr Index MIN_DIFFUSION_WEIGHTED = 6; + +// Returns erfi(x) for x >= 0: its power series below ERFI_ASYMPTOTIC, the leading asymptotic term +// above. +double erfi(double x) { + if (x >= ERFI_ASYMPTOTIC) { + return std::exp(x * x) / x / SQRT_PI; + } + // erfi(x) = 2 / sqrt(pi) sum_n x^(2n + 1) / (n! (2n + 1)); every term is positive. + const double x2 = x * x; + double power = x; + double sum = x; + for (int index = 1; index < MAX_SERIES_TERMS; ++index) { + power *= x2 / index; + const double term = power / (2.0 * index + 1.0); + sum += term; + if (term <= SERIES_TOLERANCE * sum) { + break; + } + } + return 2.0 / SQRT_PI * sum; +} + +// Returns P_0(x), P_2(x), ..., P_12(x), the even Legendre polynomials, by their recurrence. +NoddiTerms even_legendre(double x) { + NoddiTerms even{}; + double previous = 1.0; + double current = x; + even[0] = 1.0; + for (int degree = 1; degree < 2 * (NODDI_TERMS - 1); ++degree) { + const double next = ((2.0 * degree + 1.0) * x * current - degree * previous) / (degree + 1); + previous = current; + current = next; + if ((degree + 1) % 2 == 0) { + even.at(static_cast((degree + 1) / 2)) = current; + } + } + return even; +} + +// Returns log I0(x) for x >= 0. +double log_bessel_i0(double x) { + if (x >= LOG_BESSEL_ASYMPTOTIC) { + return x - 0.5 * std::log(2.0 * PI * x); + } + return std::log(scaled_bessel_i0(x)) + x; +} + +// The coefficients of the fits of the Watson coefficients of degrees 2 to 12 above +// kappa = 30, in powers 0 to 6 of ln(kappa / 30). +constexpr std::array, NODDI_TERMS - 1> WATSON_FITS{{ + {7.52308, 0.411538, -0.214588, 0.0784091, -0.023981, 0.00731537, -0.0026467}, + {8.93718, 1.62147, -0.733421, 0.191568, -0.0202906, -0.00779095, 0.00574847}, + {8.87905, 3.35689, -1.15935, 0.0673053, 0.121857, -0.066642, 0.0180215}, + {7.84352, 5.03178, -1.0193, -0.426362, 0.328816, -0.0688176, -0.0229398}, + {6.30113, 6.09914, -0.16088, -1.05578, 0.338069, 0.0937157, -0.106935}, + {4.65678, 6.30069, 1.13754, -1.38393, -0.0134758, 0.331686, -0.105954}, +}}; + +// The Watson coefficients above kappa = 30, from the fits. +NoddiTerms watson_fitted(double kappa) { + NoddiTerms coefficients{}; + coefficients[0] = 2.0 * SQRT_PI; + const double log_ratio = std::log(kappa) - std::log(WATSON_FIT_KAPPA); + for (std::size_t term = 1; term < coefficients.size(); ++term) { + double power = 1.0; + double value = 0.0; + for (const double coefficient : WATSON_FITS.at(term - 1)) { + value += coefficient * power; + power *= log_ratio; + } + coefficients.at(term) = value; + } + return coefficients; +} + +// The Watson coefficients for kappa up to 0.1, from their series. +NoddiTerms watson_series(double kappa) { + const double k2 = kappa * kappa; + const double k3 = k2 * kappa; + const double k4 = k3 * kappa; + const double k5 = k4 * kappa; + const double k6 = k5 * kappa; + return {2.0 * SQRT_PI, + (4.0 / 3.0 * kappa + 8.0 / 63.0 * k2) * std::sqrt(PI / 5.0), + (8.0 / 21.0 * k2 + 32.0 / 693.0 * k3) * (SQRT_PI * 0.2), + (16.0 / 693.0 * k3 + 32.0 / 10395.0 * k4) * std::sqrt(PI / 13.0), + 32.0 / 19305.0 * k4 * std::sqrt(PI / 17.0), + 64.0 * std::sqrt(PI / 21.0) * k5 / 692835.0, + 128.0 * SQRT_PI * k6 / 152108775.0}; +} + +// The Watson coefficients for 0.1 < kappa <= 30, in closed form. +NoddiTerms watson_closed_form(double kappa) { + const double sk = std::sqrt(kappa); + const double sk2 = sk * kappa; + const double sk3 = sk2 * kappa; + const double sk4 = sk3 * kappa; + const double sk5 = sk4 * kappa; + const double sk6 = sk5 * kappa; + const double k2 = kappa * kappa; + const double k3 = k2 * kappa; + const double k4 = k3 * kappa; + const double k5 = k4 * kappa; + const double k6 = k5 * kappa; + const double erfik = erfi(sk); + const double ierfik = 1.0 / erfik; + const double ek = std::exp(kappa); + const double dawsonk = 0.5 * SQRT_PI * erfik / ek; + + NoddiTerms coefficients{}; + coefficients[0] = 2.0 * SQRT_PI; + coefficients[1] = + std::sqrt(5.0) * (3.0 * sk - (3.0 + 2.0 * kappa) * dawsonk) * ek * ierfik / kappa; + coefficients[2] = 0.375 * + ((105.0 + 60.0 * kappa + 12.0 * k2) * dawsonk - 105.0 * sk + 10.0 * sk2) * + ek / k2 * ierfik; + coefficients[3] = ((-3465.0 - 1890.0 * kappa - 420.0 * k2 - 40.0 * k3) * dawsonk + 3465.0 * sk - + 420.0 * sk2 + 84.0 * sk3) * + std::sqrt(13.0 * PI) / 64.0 / k3 / dawsonk; + coefficients[4] = + std::sqrt(17.0) * + ((675675.0 + 360360.0 * kappa + 83160.0 * k2 + 10080.0 * k3 + 560.0 * k4) * dawsonk - + 675675.0 * sk + 90090.0 * sk2 - 23100.0 * sk3 + 744.0 * sk4) * + ek / 512.0 / k4 * ierfik; + coefficients[5] = + std::sqrt(21.0 * PI) * + ((-43648605.0 - 22972950.0 * kappa - 5405400.0 * k2 - 720720.0 * k3 - 55440.0 * k4 - + 2016.0 * k5) * + dawsonk + + 43648605.0 * sk - 6126120.0 * sk2 + 1729728.0 * sk3 - 82368.0 * sk4 + 5104.0 * sk5) / + 4096.0 / k5 / dawsonk; + coefficients[6] = 5.0 * + ((7027425405.0 + 3666482820.0 * kappa + 872972100.0 * k2 + 122522400.0 * k3 + + 10810800.0 * k4 + 576576.0 * k5 + 14784.0 * k6) * + dawsonk - + 7027425405.0 * sk + 1018467450.0 * sk2 - 302630328.0 * sk3 + + 17153136.0 * sk4 - 1553552.0 * sk5 + 25376.0 * sk6) * + ek / 16384.0 / k6 * ierfik; + return coefficients; +} + +// The parts of the signal that depend on the direction, for every volume: g . mu, and +// sqrt((2 l + 1) / (4 pi)) P_l(g . mu) of each even degree l. +struct DirectionTerms { + Vector cosines; + Matrix harmonics; +}; + +DirectionTerms direction_terms(const NoddiProtocol& protocol, const Vector3& mu) { + const Index n_volumes = protocol.bvals.size(); + DirectionTerms terms{.cosines = Vector(n_volumes), .harmonics = Matrix(n_volumes, NODDI_TERMS)}; + for (Index volume = 0; volume < n_volumes; ++volume) { + const double cosine = std::clamp(protocol.bvecs.row(volume).dot(mu), -1.0, 1.0); + terms.cosines(volume) = cosine; + const NoddiTerms legendre = even_legendre(cosine); + for (int term = 0; term < NODDI_TERMS; ++term) { + terms.harmonics(volume, term) = + std::sqrt((term + 0.25) / PI) * legendre.at(static_cast(term)); + } + } + return terms; +} + +// Returns the Watson sticks' signal of every volume. A value that is not positive, from the +// series' cut at degree 12, becomes 0.1 times the smallest positive one. +Vector stick_signal(const NoddiProtocol& protocol, const DirectionTerms& terms, double kappa) { + const NoddiTerms coefficients = watson_sh_coefficients(kappa); + const Index n_volumes = protocol.bvals.size(); + Vector signal(n_volumes); + for (Index volume = 0; volume < n_volumes; ++volume) { + double sum = 0.0; + for (int term = 0; term < NODDI_TERMS; ++term) { + sum += protocol.integrals(volume, term) * + coefficients.at(static_cast(term)) * terms.harmonics(volume, term); + } + signal(volume) = 0.5 * sum; + } + double smallest = std::numeric_limits::infinity(); + for (Index volume = 0; volume < n_volumes; ++volume) { + if (signal(volume) > 0.0) { + smallest = std::min(smallest, signal(volume)); + } + } + for (Index volume = 0; volume < n_volumes; ++volume) { + if (!(signal(volume) > 0.0)) { + signal(volume) = 0.1 * smallest; + } + } + return signal; +} + +// Returns the zeppelin's signal of every volume, with the tortuous perpendicular diffusivity +// d (1 - v_ic) averaged over the Watson density. +Vector zeppelin_signal(const NoddiProtocol& protocol, const NoddiSettings& settings, + const DirectionTerms& terms, double v_ic, double kappa) { + const auto [along, across] = + watson_hindered_diffusivities(settings.d_par, settings.d_par * (1.0 - v_ic), kappa); + const Index n_volumes = protocol.bvals.size(); + Vector signal(n_volumes); + for (Index volume = 0; volume < n_volumes; ++volume) { + const double c2 = terms.cosines(volume) * terms.cosines(volume); + signal(volume) = std::exp(-protocol.bvals(volume) * (across + (along - across) * c2)); + } + return signal; +} + +// Returns the dictionary about mu: one column per (v_ic, kappa) of the grids, v_ic outer, and +// free water last, (N, n_v_ic n_kappa + 1). +Matrix noddi_dictionary(const NoddiProtocol& protocol, const NoddiSettings& settings, + const DirectionTerms& terms) { + const Index n_volumes = protocol.bvals.size(); + const auto n_kappa = static_cast(settings.kappa_grid.size()); + const auto n_atoms = static_cast(settings.v_ic_grid.size()) * n_kappa; + Matrix dictionary(n_volumes, n_atoms + 1); + std::vector sticks; + sticks.reserve(settings.kappa_grid.size()); + for (const double kappa : settings.kappa_grid) { + sticks.push_back(stick_signal(protocol, terms, kappa)); + } + Index column = 0; + for (const double v_ic : settings.v_ic_grid) { + for (std::size_t index = 0; index < sticks.size(); ++index) { + const double kappa = settings.kappa_grid.at(index); + dictionary.col(column) = + v_ic * sticks.at(index) + + (1.0 - v_ic) * zeppelin_signal(protocol, settings, terms, v_ic, kappa); + ++column; + } + } + dictionary.col(n_atoms) = protocol.free_water; + return dictionary; +} + +// The start of the fit: the dictionary estimate, or no tissue. +struct NoddiStart { + double v_ic = 0.0; + double kappa = 0.0; + double v_iso = 0.0; + bool has_tissue = false; +}; + +// Solves the dictionary: free water from a plain NNLS, the tissue atoms' support +// from a non-negative L1 + L2 solve of the rest, their weights from a plain NNLS on the support, +// and v_ic and kappa as the weighted means of the atoms' values. +NoddiStart dictionary_start(const NoddiSettings& settings, const Matrix& dictionary, + const Vector& normalised) { + const Index n_atoms = dictionary.cols() - 1; + const auto solve = [](const Matrix& design, const Vector& data, double ridge, double l1) { + const sd::RidgeDesign prepared = sd::make_ridge_design(design, ridge); + return sd::solve_ridge_nnls(prepared, data, l1, sd::default_gradient_tolerance(design)); + }; + NoddiStart start; + start.v_iso = solve(dictionary, normalised, 0.0, 0.0)(n_atoms); + const Matrix tissue = dictionary.leftCols(n_atoms); + const Vector without_water = normalised - start.v_iso * dictionary.col(n_atoms); + Vector weights = solve(tissue, without_water, settings.ridge, settings.l1); + + std::vector support; + for (Index atom = 0; atom < n_atoms; ++atom) { + if (weights(atom) > 0.0) { + support.push_back(atom); + } + } + if (support.empty()) { + return start; + } + const Matrix kept = tissue(Eigen::placeholders::all, support); + const Vector refitted = solve(kept, without_water, 0.0, 0.0); + double total = 0.0; + const auto n_kappa = static_cast(settings.kappa_grid.size()); + for (std::size_t index = 0; index < support.size(); ++index) { + const double weight = refitted(static_cast(index)); + const Index atom = support.at(index); + total += weight; + start.v_ic += weight * settings.v_ic_grid.at(static_cast(atom / n_kappa)); + start.kappa += weight * settings.kappa_grid.at(static_cast(atom % n_kappa)); + } + if (!(total > 0.0)) { + start.v_ic = 0.0; + start.kappa = 0.0; + return start; + } + start.v_ic /= total; + start.kappa /= total; + start.has_tissue = true; + return start; +} + +// The five fitted parameters, encoded: asin(sqrt(v_ic)), sqrt(kappa / 10), +// asin(sqrt(v_iso)), and the polar and azimuthal angles of mu. +using Encoded = Eigen::Matrix; + +Encoded encode(double v_ic, double kappa, double v_iso, const Vector3& mu) { + Encoded x; + x << std::asin(std::sqrt(std::clamp(v_ic, 0.0, 1.0))), + std::sqrt(std::min(kappa, KAPPA_MAX) / KAPPA_SCALE), + std::asin(std::sqrt(std::clamp(v_iso, 0.0, 1.0))), std::acos(std::clamp(mu.z(), -1.0, 1.0)), + std::atan2(mu.y(), mu.x()); + return x; +} + +Vector3 direction_of(const Encoded& x) { + return {std::cos(x(4)) * std::sin(x(3)), std::sin(x(4)) * std::sin(x(3)), std::cos(x(3))}; +} + +// The measurements a voxel's fit uses, its positive ones. +struct Measurements { + NoddiProtocol protocol; + Vector signal; + double s0 = 0.0; + double sigma = 0.0; +}; + +// Returns the voxel's positive measurements with their S0 and noise level, or false when there is +// no positive b0 or too few positive diffusion-weighted measurements. +bool positive_measurements(const NoddiProtocol& protocol, const VectorView& signal, + Measurements& out) { + std::vector kept; + std::vector b0s; + Index n_weighted = 0; + for (Index volume = 0; volume < signal.size(); ++volume) { + if (signal(volume) > 0.0) { + kept.push_back(volume); + if (protocol.is_b0.at(static_cast(volume))) { + b0s.push_back(signal(volume)); + } else { + ++n_weighted; + } + } + } + if (b0s.empty() || n_weighted < MIN_DIFFUSION_WEIGHTED) { + return false; + } + out.protocol.bvals = protocol.bvals(kept); + out.protocol.bvecs = protocol.bvecs(kept, Eigen::placeholders::all); + out.protocol.integrals = protocol.integrals(kept, Eigen::placeholders::all); + out.protocol.free_water = protocol.free_water(kept); + out.signal = signal(kept); + double mean = 0.0; + for (const double b0 : b0s) { + mean += b0; + } + mean /= static_cast(b0s.size()); + // The spread about the mean, over n rather than n - 1. + double spread = 0.0; + for (const double b0 : b0s) { + spread += (b0 - mean) * (b0 - mean); + } + spread = std::sqrt(spread / static_cast(b0s.size())); + out.s0 = mean; + out.sigma = std::max(spread, MIN_SNR * mean) / SIGMA_SCALING; + return true; +} + +// Returns the negative Rician log-likelihood of the measurements for the encoded parameters. Like +// the toolbox, it adds log E where the Rician density has log m. +double negative_log_likelihood(const Measurements& data, const NoddiSettings& settings, + const Encoded& x) { + const double v_ic = std::sin(x(0)) * std::sin(x(0)); + const double kappa = KAPPA_SCALE * x(1) * x(1); + const double v_iso = std::sin(x(2)) * std::sin(x(2)); + const DirectionTerms terms = direction_terms(data.protocol, direction_of(x)); + const Vector sticks = stick_signal(data.protocol, terms, kappa); + const Vector zeppelin = zeppelin_signal(data.protocol, settings, terms, v_ic, kappa); + const double variance = data.sigma * data.sigma; + double total = 0.0; + for (Index volume = 0; volume < data.signal.size(); ++volume) { + const double tissue = v_ic * sticks(volume) + (1.0 - v_ic) * zeppelin(volume); + const double model = + data.s0 * ((1.0 - v_iso) * tissue + v_iso * data.protocol.free_water(volume)); + const double measured = data.signal(volume); + total += -2.0 * std::log(data.sigma) - + (model * model + measured * measured) / (2.0 * variance) + std::log(model) + + log_bessel_i0(measured * model / variance); + } + return std::isfinite(total) ? -total : std::numeric_limits::infinity(); +} + +// Returns the lower bounds of the encoded parameters. +Encoded lower_bounds() { + Encoded lower; + lower << 0.0, 0.0, 0.0, -ANGLE_LIMIT, -ANGLE_LIMIT; + return lower; +} + +// Returns the upper bounds of the encoded parameters. +Encoded upper_bounds() { + Encoded upper; + upper << FRACTION_UPPER, std::sqrt(KAPPA_MAX / KAPPA_SCALE), FRACTION_UPPER, ANGLE_LIMIT, + ANGLE_LIMIT; + return upper; +} + +using Hessian = Eigen::Matrix; +using Objective = std::function; +// Whether each encoded parameter is fitted or held. +using FreeParameters = std::array; + +// Returns the forward-difference gradient of the objective at point, where it has the given value; +// the parameters that are not free get 0. +Encoded forward_gradient(const Objective& objective, const FreeParameters& free, + const Encoded& point, double value) { + const Encoded upper = upper_bounds(); + Encoded gradient = Encoded::Zero(); + for (Index param = 0; param < 5; ++param) { + if (!free.at(static_cast(param))) { + continue; + } + double step = std::sqrt(std::numeric_limits::epsilon()) * + std::max(std::abs(point(param)), 1.0); + if (point(param) + step > upper(param)) { + step = -step; + } + Encoded moved = point; + moved(param) += step; + gradient(param) = (objective(moved) - value) / step; + } + return gradient; +} + +// Returns the quasi-Newton direction from x, holding the parameters that are not free and those at +// a bound that the gradient pushes out of. When that is no descent direction, the inverse Hessian +// is reset, restarted is set, and the direction is the steepest descent of the free parameters. +Encoded descent_direction(Hessian& inverse_hessian, const Encoded& gradient, const Encoded& x, + const FreeParameters& free, bool& restarted) { + const Encoded lower = lower_bounds(); + const Encoded upper = upper_bounds(); + Encoded direction = -inverse_hessian * gradient; + for (Index param = 0; param < 5; ++param) { + const bool held = !free.at(static_cast(param)) || + (x(param) <= lower(param) && gradient(param) > 0.0) || + (x(param) >= upper(param) && gradient(param) < 0.0); + if (held) { + direction(param) = 0.0; + } + } + if (direction.dot(gradient) < 0.0) { + return direction; + } + inverse_hessian.setIdentity(); + restarted = true; + direction = -gradient; + for (Index param = 0; param < 5; ++param) { + if (!free.at(static_cast(param))) { + direction(param) = 0.0; + } + } + return direction; +} + +// A point found by the line search, the objective there, and the fraction of the step taken. +struct LineStep { + Encoded point; + double value = 0.0; + double length = 1.0; + bool accepted = false; +}; + +// Halves the step along direction from x, projected onto the bounds, until the objective falls +// enough below value (the Armijo condition). +LineStep line_search(const Objective& objective, const Encoded& x, double value, + const Encoded& gradient, const Encoded& direction) { + const Encoded lower = lower_bounds(); + const Encoded upper = upper_bounds(); + LineStep step{.point = x, .value = value, .length = 1.0, .accepted = false}; + for (int halving = 0; halving < MAX_HALVINGS; ++halving) { + step.point = (x + step.length * direction).cwiseMax(lower).cwiseMin(upper); + step.value = objective(step.point); + if (step.value <= value + ARMIJO * gradient.dot(step.point - x)) { + step.accepted = true; + return step; + } + step.length *= 0.5; + } + return step; +} + +// Updates the inverse Hessian by BFGS from a step and the change of the gradient over it, scaled +// first when it is still the identity. Returns false, changing nothing, when the curvature along +// the step is not positive. +bool update_inverse_hessian(Hessian& inverse_hessian, const Encoded& step, + const Encoded& difference, bool scale) { + const double curvature = step.dot(difference); + if (!(curvature > 0.0)) { + return false; + } + if (scale) { + inverse_hessian *= curvature / difference.squaredNorm(); + } + const double rho = 1.0 / curvature; + const Hessian identity = Hessian::Identity(); + inverse_hessian = (identity - rho * step * difference.transpose()) * inverse_hessian * + (identity - rho * difference * step.transpose()) + + rho * step * step.transpose(); + return true; +} + +// Minimises the objective over the free encoded parameters within the bounds, from x, by a +// quasi-Newton method (BFGS) with forward-difference gradients and a backtracking line search. +// Returns the objective at the result. +double minimise(const Objective& objective, const FreeParameters& free, Encoded& x) { + double value = objective(x); + if (!std::isfinite(value)) { + return value; + } + Encoded gradient = forward_gradient(objective, free, x, value); + Hessian inverse_hessian = Hessian::Identity(); + bool first = true; + for (int iteration = 0; iteration < MAX_ITERATIONS; ++iteration) { + Encoded direction = descent_direction(inverse_hessian, gradient, x, free, first); + const double largest = direction.cwiseAbs().maxCoeff(); + if (largest == 0.0) { + break; + } + if (first) { + direction *= FIRST_STEP / largest; + } + const LineStep next = line_search(objective, x, value, gradient, direction); + if (!next.accepted) { + break; + } + const Encoded step = next.point - x; + const double change = value - next.value; + x = next.point; + value = next.value; + // A step shortened by the line search is no sign of convergence; only a full one is. + const bool full_step = next.length == 1.0 && !first; + if ((full_step && step.cwiseAbs().maxCoeff() < STEP_TOLERANCE) || + change <= OBJECTIVE_TOLERANCE * std::abs(value)) { + break; + } + const Encoded next_gradient = forward_gradient(objective, free, x, value); + if (update_inverse_hessian(inverse_hessian, step, next_gradient - gradient, first)) { + first = false; + } + gradient = next_gradient; + } + return value; +} + +} // namespace + +NoddiTerms watson_sh_coefficients(double kappa) { + if (kappa > WATSON_FIT_KAPPA) { + return watson_fitted(kappa); + } + if (kappa > WATSON_SERIES_KAPPA) { + return watson_closed_form(kappa); + } + return watson_series(kappa); +} + +NoddiTerms legendre_gaussian_integrals(double x) { + if (x <= INTEGRAL_SERIES_X) { + const double x2 = x * x; + const double x3 = x2 * x; + const double x4 = x3 * x; + const double x5 = x4 * x; + const double x6 = x5 * x; + return {2.0 - 2.0 * x / 3.0 + x2 / 5.0 - x3 / 21.0 + x4 / 108.0, + -4.0 * x / 15.0 + 4.0 * x2 / 35.0 - 2.0 * x3 / 63.0 + 2.0 * x4 / 297.0, + 8.0 * x2 / 315.0 - 8.0 * x3 / 693.0 + 4.0 * x4 / 1287.0, + -16.0 * x3 / 9009.0 + 16.0 * x4 / 19305.0, + 32.0 * x4 / 328185.0, + -64.0 * x5 / 14549535.0, + 128.0 * x6 / 760543875.0}; + } + // I_i = integral of u^(2i) exp(-x u^2) over [-1, 1], by recurrence. + std::array moments{}; + const double root = std::sqrt(x); + moments[0] = SQRT_PI * std::erf(root) / root; + const double minus_exp = -std::exp(-x); + for (std::size_t power = 1; power < moments.size(); ++power) { + moments.at(power) = + (minus_exp + (static_cast(power) - 0.5) * moments.at(power - 1)) / x; + } + const auto& mom = moments; + return {mom[0], + -0.5 * mom[0] + 1.5 * mom[1], + 0.375 * mom[0] - 3.75 * mom[1] + 4.375 * mom[2], + -0.3125 * mom[0] + 6.5625 * mom[1] - 19.6875 * mom[2] + 14.4375 * mom[3], + 0.2734375 * mom[0] - 9.84375 * mom[1] + 54.140625 * mom[2] - 93.84375 * mom[3] + + 50.2734375 * mom[4], + (-63.0 * mom[0] + 3465.0 * mom[1] - 30030.0 * mom[2] + 90090.0 * mom[3] - + 109395.0 * mom[4] + 46189.0 * mom[5]) / + 256.0, + (231.0 * mom[0] - 18018.0 * mom[1] + 225225.0 * mom[2] - 1021020.0 * mom[3] + + 2078505.0 * mom[4] - 1939938.0 * mom[5] + 676039.0 * mom[6]) / + 1024.0}; +} + +std::array watson_hindered_diffusivities(double d_par, double d_perp, double kappa) { + const double difference = d_par - d_perp; + if (kappa < HINDERED_SERIES_KAPPA) { + const double mean_sum = d_par + 2.0 * d_perp; + const double k2 = kappa * kappa; + return {mean_sum / 3.0 + 4.0 * difference * kappa / 45.0 + 8.0 * difference * k2 / 945.0, + mean_sum / 3.0 - 2.0 * difference * kappa / 45.0 - 4.0 * difference * k2 / 945.0}; + } + const double root = std::sqrt(kappa); + const double dawson = 0.5 * std::exp(-kappa) * SQRT_PI * erfi(root); + const double factor = root / dawson; + return {(-difference + 2.0 * d_perp * kappa + difference * factor) / (2.0 * kappa), + (difference + 2.0 * (d_par + d_perp) * kappa - difference * factor) / (4.0 * kappa)}; +} + +NoddiProtocol make_noddi_protocol(const VectorView& bvals, const MatrixView& bvecs, + const std::vector& is_b0, const NoddiSettings& settings) { + const Index n_volumes = bvals.size(); + NoddiProtocol protocol{.bvals = bvals, + .bvecs = bvecs, + .is_b0 = is_b0, + .integrals = Matrix(n_volumes, NODDI_TERMS), + .free_water = Vector(n_volumes)}; + for (Index volume = 0; volume < n_volumes; ++volume) { + const NoddiTerms integrals = legendre_gaussian_integrals(bvals(volume) * settings.d_par); + for (int term = 0; term < NODDI_TERMS; ++term) { + protocol.integrals(volume, term) = integrals.at(static_cast(term)); + } + protocol.free_water(volume) = std::exp(-bvals(volume) * settings.d_iso); + } + return protocol; +} + +Vector noddi_signal(const NoddiProtocol& protocol, const NoddiSettings& settings, double v_ic, + double kappa, double v_iso, const Vector3& mu) { + const DirectionTerms terms = direction_terms(protocol, mu); + const Vector tissue = v_ic * stick_signal(protocol, terms, kappa) + + (1.0 - v_ic) * zeppelin_signal(protocol, settings, terms, v_ic, kappa); + return (1.0 - v_iso) * tissue + v_iso * protocol.free_water; +} + +NoddiVoxelFit fit_noddi_voxel(const NoddiProtocol& protocol, const NoddiSettings& settings, + const VectorView& signal, const Vector3& mu) { + NoddiVoxelFit fit; + double b0_sum = 0.0; + Index n_b0 = 0; + for (Index volume = 0; volume < signal.size(); ++volume) { + if (protocol.is_b0.at(static_cast(volume))) { + b0_sum += signal(volume); + ++n_b0; + } + } + const double b0_mean = n_b0 > 0 ? b0_sum / static_cast(n_b0) : 0.0; + Measurements data; + const double mu_norm = mu.norm(); + if (!(b0_mean > 0.0) || !signal.allFinite() || !std::isfinite(mu_norm) || mu_norm == 0.0 || + !positive_measurements(protocol, signal, data)) { + fit.fitted = false; + return fit; + } + const Vector3 start_mu = mu / mu_norm; + + const Matrix dictionary = + noddi_dictionary(protocol, settings, direction_terms(protocol, start_mu)); + const NoddiStart start = dictionary_start(settings, dictionary, signal / b0_mean); + fit.v_iso = start.v_iso; + fit.mu = start_mu; + if (!start.has_tissue) { + return fit; + } + + const auto objective = [&](const Encoded& x) { + return negative_log_likelihood(data, settings, x); + }; + // v = sin^2(u) has no slope at v = 0 and 1, where a fit could not move, so the start keeps + // the fractions just inside. + const Encoded initial = + encode(std::clamp(start.v_ic, START_MARGIN, 1.0 - START_MARGIN), start.kappa, + std::clamp(start.v_iso, START_MARGIN, 1.0 - START_MARGIN), start_mu); + const double initial_value = objective(initial); + // Two stages: the fractions with kappa and mu held, then everything. + Encoded x = initial; + minimise(objective, {true, false, true, false, false}, x); + double value = minimise(objective, {true, true, true, true, true}, x); + // A fraction that ends at 0 or 1 cannot move again, so the second stage runs once more with it + // just inside, and the better of the two fits is kept. + Encoded moved = x; + for (const Index fraction : {Index{0}, Index{2}}) { + const double share = std::sin(x(fraction)) * std::sin(x(fraction)); + if (share < START_MARGIN || share > 1.0 - START_MARGIN) { + moved(fraction) = + std::asin(std::sqrt(std::clamp(share, START_MARGIN, 1.0 - START_MARGIN))); + } + } + if (moved != x) { + const double moved_value = minimise(objective, {true, true, true, true, true}, moved); + if (moved_value < value) { + x = moved; + value = moved_value; + } + } + if (!(value < initial_value)) { + x = initial; + } + fit.v_ic = std::sin(x(0)) * std::sin(x(0)); + fit.kappa = KAPPA_SCALE * x(1) * x(1); + fit.v_iso = std::sin(x(2)) * std::sin(x(2)); + fit.mu = direction_of(x); + return fit; +} + +} // namespace larmor::dwi diff --git a/cpp/src/dwi/noddi.hpp b/cpp/src/dwi/noddi.hpp new file mode 100644 index 0000000..8338aa6 --- /dev/null +++ b/cpp/src/dwi/noddi.hpp @@ -0,0 +1,97 @@ +// Computes the NODDI signal and fits it to each voxel, from a start found in a dictionary of NODDI +// signals, by maximum likelihood with Rician noise. +// +// References: Zhang et al., NeuroImage 61 (2012) 1000-1016; NODDI Matlab Toolbox 1.05; Daducci et +// al., NeuroImage 105 (2015) 32-44. +#pragma once + +#include +#include + +#include "common/types.hpp" + +namespace larmor::dwi { + +/// The number of even Legendre terms in the toolbox's Watson stick, degrees 0 to 12. +constexpr int NODDI_TERMS = 7; + +/// One value per even Legendre degree 0, 2, ..., 12. +using NoddiTerms = std::array; + +/// Returns the coefficients of the Watson density of concentration kappa >= 0 in the even +/// spherical harmonics: by series up to kappa = 0.1, in closed form up to kappa = 30, and by fits +/// in ln(kappa) above. +NoddiTerms watson_sh_coefficients(double kappa); + +/// Returns the integrals over u in [-1, 1] of exp(-x u^2) P_2i(u), i = 0 to 6, for x = b d >= 0. +NoddiTerms legendre_gaussian_integrals(double x); + +/// Returns the parallel and perpendicular diffusivities of a cylindrically symmetric tensor with +/// diffusivities d_par and d_perp, averaged over a Watson density of concentration kappa. +std::array watson_hindered_diffusivities(double d_par, double d_perp, double kappa); + +/// The fixed values of the model and of the fit. +struct NoddiSettings { + /// The diffusivity along the neurites, mm^2/s. + double d_par = 1.7e-3; + /// The diffusivity of free water, mm^2/s. + double d_iso = 3.0e-3; + /// The neurite densities of the dictionary's atoms. + std::vector v_ic_grid; + /// The Watson concentrations of the dictionary's atoms. + std::vector kappa_grid; + /// The ridge of the dictionary's sparse solve, the square root of its lambda. + double ridge = 0.0; + /// The L1 weight of the dictionary's sparse solve. + double l1 = 0.0; +}; + +/// A gradient table prepared once for the fits of every voxel. +struct NoddiProtocol { + /// The b-value of each volume, s/mm^2, (N,). + Vector bvals; + /// The unit gradient direction of each volume, (N, 3). + Matrix bvecs; + /// Whether each volume is a b0. + std::vector is_b0; + /// legendre_gaussian_integrals(b d_par) of each volume, (N, NODDI_TERMS). + Matrix integrals; + /// exp(-b d_iso) of each volume, (N,). + Vector free_water; +}; + +/// Prepares a gradient table: b-values (N,), unit directions (N, 3) and the b0 flags. +NoddiProtocol make_noddi_protocol(const VectorView& bvals, const MatrixView& bvecs, + const std::vector& is_b0, const NoddiSettings& settings); + +/// Returns the normalised NODDI signal of every volume for the neurite density v_ic, Watson +/// concentration kappa and free-water fraction v_iso about the unit direction mu, (N,). +Vector noddi_signal(const NoddiProtocol& protocol, const NoddiSettings& settings, double v_ic, + double kappa, double v_iso, const Vector3& mu); + +/// The fit of one voxel. +struct NoddiVoxelFit { + /// The neurite density. + double v_ic = 0.0; + /// The Watson concentration. + double kappa = 0.0; + /// The free-water fraction. + double v_iso = 0.0; + /// The fitted mean direction, a unit vector, or zero when the voxel is not fitted. + Vector3 mu = Vector3::Zero(); + /// False for a voxel that is not fitted and gets zeros: one with a signal that is not finite, + /// without a positive mean b0, with fewer than six positive diffusion-weighted measurements, + /// or with a starting direction that is zero or not finite. + bool fitted = true; +}; + +/// Fits one voxel's signal (N,) from the starting direction mu, which is normalised. +/// +/// The start comes from the dictionary of NODDI signals about mu; the parameters and the direction +/// are then fitted by maximum likelihood with Rician noise, first the fractions and then all of +/// them. A voxel whose dictionary solve finds no tissue gets v_ic = kappa = 0 and the +/// dictionary's v_iso, without the second step. +NoddiVoxelFit fit_noddi_voxel(const NoddiProtocol& protocol, const NoddiSettings& settings, + const VectorView& signal, const Vector3& mu); + +} // namespace larmor::dwi diff --git a/cpp/src/python/dwi.cpp b/cpp/src/python/dwi.cpp index 0a1dc61..68dbebd 100644 --- a/cpp/src/python/dwi.cpp +++ b/cpp/src/python/dwi.cpp @@ -1,4 +1,5 @@ -// Bindings of the log-signal fits of DTI and DKI, MP-PCA denoising and Gibbs-ringing removal. +// Bindings of the log-signal fits of DTI and DKI, the NODDI fit, MP-PCA denoising and Gibbs-ringing +// removal. #include #include @@ -6,6 +7,7 @@ #include #include #include +#include #include #include @@ -18,6 +20,7 @@ #include "dwi/gibbs.hpp" #include "dwi/log_fit.hpp" #include "dwi/mppca.hpp" +#include "dwi/noddi.hpp" #include "python/arrays.hpp" #include "python/bindings.hpp" @@ -31,6 +34,8 @@ namespace { constexpr Index BLOCK_SIZE = 64; // Voxels per block of the parallel loop of the denoising. constexpr Index DENOISE_BLOCK_SIZE = 8; +// Voxels per block of the parallel loop of the NODDI fit. +constexpr Index NODDI_BLOCK_SIZE = 16; // Fits the log signal of every voxel and returns the parameters and the constrained flags. std::pair, NumpyArray> @@ -76,6 +81,101 @@ fit_log_signal(const InArray& design_array, const InArray& return {to_numpy(std::move(parameters)), to_numpy(std::move(constrained))}; } +// Returns "(rows, columns)", for error messages. +std::string shape_tuple(Index rows, Index columns) { + return "(" + std::to_string(rows) + ", " + std::to_string(columns) + ")"; +} + +// Fits NODDI to every voxel and returns the parameters and which voxels were fitted. +std::pair, NumpyArray> +fit_noddi(const InArray& signals_array, const InArray& directions_array, + const InArray& bvals_array, const InArray& bvecs_array, + const InArray& is_b0_array, const InArray& v_ic_grid_array, + const InArray& kappa_grid_array, double d_par, double d_iso, double ridge, + double l1, int n_threads) { + const auto signals = matrix_view(signals_array); // checks contiguity + const auto directions = matrix_view(directions_array); + const auto bvals = flat_view(bvals_array); + const auto bvecs = matrix_view(bvecs_array); + const auto is_b0 = flat_view(is_b0_array); + const auto v_ic_grid = flat_view(v_ic_grid_array); + const auto kappa_grid = flat_view(kappa_grid_array); + const Index n_volumes = bvals.size(); + if (signals.cols() != n_volumes || bvecs.rows() != n_volumes || bvecs.cols() != 3 || + is_b0.size() != n_volumes || directions.rows() != signals.rows() || + directions.cols() != 3) { + throw nb::value_error(("signals must be (v, n_volumes), mu (v, 3), bvals (n_volumes,), " + "bvecs (n_volumes, 3) and is_b0 (n_volumes,): got " + "signals " + + shape_tuple(signals.rows(), signals.cols()) + ", mu " + + shape_tuple(directions.rows(), directions.cols()) + ", bvals (" + + std::to_string(n_volumes) + ",), bvecs " + + shape_tuple(bvecs.rows(), bvecs.cols()) + ", is_b0 (" + + std::to_string(is_b0.size()) + ",)") + .c_str()); + } + if (v_ic_grid.size() == 0 || kappa_grid.size() == 0) { + throw nb::value_error("the v_ic and kappa grids must not be empty"); + } + if (!(ridge > 0.0) || !(l1 >= 0.0)) { + throw nb::value_error("ridge must be positive and l1 at least 0"); + } + require_threads(n_threads); + + const dwi::NoddiSettings settings{ + .d_par = d_par, + .d_iso = d_iso, + .v_ic_grid = std::vector(v_ic_grid.begin(), v_ic_grid.end()), + .kappa_grid = std::vector(kappa_grid.begin(), kappa_grid.end()), + .ridge = ridge, + .l1 = l1}; + const Index n_voxels = signals.rows(); + Matrix parameters = Matrix::Zero(n_voxels, 6); + dwi::Flags fitted(n_voxels); + { + const nb::gil_scoped_release release; + const dwi::NoddiProtocol protocol = dwi::make_noddi_protocol( + Vector(bvals), Matrix(bvecs), std::vector(is_b0.begin(), is_b0.end()), settings); + parallel_for_blocks( + n_voxels, NODDI_BLOCK_SIZE, n_threads, [] { return 0; }, + [&](Index begin, Index end, int& /*scratch*/) { + for (Index voxel = begin; voxel < end; ++voxel) { + const Vector signal = signals.row(voxel).transpose(); + const Vector3 mu = directions.row(voxel).transpose(); + const dwi::NoddiVoxelFit fit = + dwi::fit_noddi_voxel(protocol, settings, signal, mu); + parameters.row(voxel) << fit.v_ic, fit.kappa, fit.v_iso, fit.mu.x(), fit.mu.y(), + fit.mu.z(); + fitted(voxel) = fit.fitted; + } + }); + } + return {to_numpy(std::move(parameters)), to_numpy(std::move(fitted))}; +} + +// Returns the normalised NODDI signal of every volume. +NumpyArray noddi_signal(const InArray& bvals_array, + const InArray& bvecs_array, double v_ic, double kappa, + double v_iso, const std::array& mu, double d_par, + double d_iso) { + const auto bvals = flat_view(bvals_array); // checks contiguity + const auto bvecs = matrix_view(bvecs_array); + if (bvecs.rows() != bvals.size() || bvecs.cols() != 3) { + throw nb::value_error(("bvals must be (n,) and bvecs (n, 3): got (" + + std::to_string(bvals.size()) + ",) and " + + shape_tuple(bvecs.rows(), bvecs.cols())) + .c_str()); + } + const dwi::NoddiSettings settings{ + .d_par = d_par, .d_iso = d_iso, .v_ic_grid = {}, .kappa_grid = {}, .ridge = 0.0, .l1 = 0.0}; + const dwi::NoddiProtocol protocol = dwi::make_noddi_protocol( + Vector(bvals), Matrix(bvecs), std::vector(static_cast(bvals.size())), + settings); + Vector signal = dwi::noddi_signal(protocol, settings, v_ic, kappa, v_iso, + Vector3(mu[0], mu[1], mu[2]).normalized()); + return to_numpy(std::move(signal)); +} + // The denoised signals, the noise levels and the numbers of signal components. using DenoisedArrays = std::tuple, NumpyArray, NumpyArray>; @@ -229,6 +329,35 @@ void register_dwi(nb::module_& module) { "Raises:\n" " ValueError: If an array is not C-contiguous, the shapes do not match,\n" " max_iter is negative or n_threads is 0.\n"); + module.def("fit_noddi", &fit_noddi, nb::arg("signals").noconvert(), nb::arg("mu").noconvert(), + nb::arg("bvals").noconvert(), nb::arg("bvecs").noconvert(), + nb::arg("is_b0").noconvert(), nb::arg("v_ic_grid").noconvert(), + nb::arg("kappa_grid").noconvert(), nb::arg("d_par"), nb::arg("d_iso"), + nb::arg("ridge"), nb::arg("l1"), nb::arg("n_threads"), + "Fits NODDI to each voxel's signal, a row of signals (v, n), from the direction\n" + "mu (v, 3).\n" + "\n" + "The start comes from a dictionary of NODDI signals over v_ic_grid and\n" + "kappa_grid, solved with the ridge and l1 weights; the fit is then by maximum\n" + "likelihood with Rician noise.\n" + "d_par and d_iso are in mm^2/s and bvals in s/mm^2.\n" + "\n" + "Returns v_ic, kappa, v_iso and the fitted direction of each voxel (v, 6), and\n" + "whether each voxel was fitted (v,).\n" + "\n" + "Raises:\n" + " ValueError: If an array is not C-contiguous, the shapes do not match, a\n" + " grid is empty, ridge is not positive, l1 is negative or n_threads is 0.\n"); + module.def("noddi_signal", &noddi_signal, nb::arg("bvals").noconvert(), + nb::arg("bvecs").noconvert(), nb::arg("v_ic"), nb::arg("kappa"), nb::arg("v_iso"), + nb::arg("mu"), nb::arg("d_par"), nb::arg("d_iso"), + "Returns the normalised NODDI signal of each volume for the neurite density\n" + "v_ic, Watson concentration kappa and free-water fraction v_iso about the\n" + "direction mu. bvals (n,) are in s/mm^2, bvecs (n, 3) unit directions, and\n" + "d_par and d_iso in mm^2/s.\n" + "\n" + "Raises:\n" + " ValueError: If an array is not C-contiguous or the shapes do not match.\n"); module.def("denoise_mppca", &denoise_mppca, nb::arg("image").noconvert(), nb::arg("voxels").noconvert(), nb::arg("patch"), nb::arg("n_threads"), "Denoises by MP-PCA the voxels at flat indices voxels (k,) int64 of a\n" diff --git a/cpp/src/python/sd.cpp b/cpp/src/python/sd.cpp index 2b24847..413aece 100644 --- a/cpp/src/python/sd.cpp +++ b/cpp/src/python/sd.cpp @@ -1,5 +1,6 @@ // Bindings of the spherical-deconvolution code. +#include #include #include #include @@ -226,6 +227,36 @@ NumpyArray lsqnonneg(const InArray& design_array, return to_numpy(std::move(solution)); } +// Solves one non-negative least-squares problem with a ridge and an L1 penalty. +NumpyArray ridge_l1_nnls(const InArray& design_array, + const InArray& data_array, double ridge, double l1, + std::optional gradient_tolerance) { + const auto design = matrix_view(design_array); // checks contiguity + const Vector data = copy_vector(data_array); + if (data.size() != design.rows()) { + throw nb::value_error(("c must be (m, n) and d (m,): got " + + shape_tuple(design.rows(), design.cols()) + " and " + + std::to_string(data.size())) + .c_str()); + } + // Written so that NaN fails too. + if (!std::isfinite(ridge) || ridge < 0.0 || !std::isfinite(l1) || l1 < 0.0) { + throw nb::value_error("ridge and l1 must be finite and >= 0."); + } + if (l1 > 0.0 && ridge == 0.0) { + throw nb::value_error("l1 > 0 needs ridge > 0."); + } + Vector solution; + { + const nb::gil_scoped_release release; + const double tolerance = + gradient_tolerance ? *gradient_tolerance : sd::default_gradient_tolerance(design); + const sd::RidgeDesign ridge_design = sd::make_ridge_design(design, ridge); + solution = sd::solve_ridge_nnls(ridge_design, data, l1, tolerance); + } + return to_numpy(std::move(solution)); +} + // Constructs a TissueKernel in place from copies of the kernel matrices. void construct_kernel(sd::TissueKernel* kernel, const InArray& anisotropic_array, const InArray& isotropic_array, @@ -447,6 +478,19 @@ void register_sd(nb::module_& module) { "Raises:\n" " ValueError: If c is not C-contiguous or the shapes do not match.\n"); + module.def("ridge_l1_nnls", &ridge_l1_nnls, nb::arg("c").noconvert(), nb::arg("d").noconvert(), + nb::arg("ridge"), nb::arg("l1"), nb::arg("tol_x") = nb::none(), + "Solves min ||c x - d||^2 + ridge^2 ||x||^2 + 2 l1 sum(x) with x >= 0.\n" + "\n" + "c is (m, n) and d (m,). Since x >= 0, the last term is an L1 penalty, which\n" + "sets the coefficients the data do not support to zero; ridge = l1 = 0 is\n" + "plain NNLS. A coefficient enters the solution only when its gradient\n" + "exceeds tol_x; None uses 10 * eps * norm(c, 1) * max(m, n).\n" + "\n" + "Raises:\n" + " ValueError: If c is not C-contiguous, the shapes do not match, ridge or\n" + " l1 is negative or not finite, or l1 > 0 with ridge = 0.\n"); + nb::class_(module, "DeconvKernel", "One resolution of a multi-tissue kernel, copied into C++.") .def("__init__", &construct_kernel, nb::arg("k_aniso").noconvert(), diff --git a/cpp/src/python/simulation.cpp b/cpp/src/python/simulation.cpp index 9656de5..e3f67f6 100644 --- a/cpp/src/python/simulation.cpp +++ b/cpp/src/python/simulation.cpp @@ -78,11 +78,14 @@ NumpyArray compartment_attenuation(const InArray& mo const simulation::PulseTiming timing{.small_delta_ms = small_delta_ms, .big_delta_ms = big_delta_ms}; Matrix result(static_cast(compartments.size()), bvals.size()); - for (Index row = 0; row < result.rows(); ++row) { - for (Index column = 0; column < result.cols(); ++column) { - result(row, column) = - simulation::attenuation(compartments[static_cast(row)], bvals(column), - cosines(column), timing); + { + const nb::gil_scoped_release release; + for (Index row = 0; row < result.rows(); ++row) { + for (Index column = 0; column < result.cols(); ++column) { + result(row, column) = + simulation::attenuation(compartments[static_cast(row)], + bvals(column), cosines(column), timing); + } } } return to_numpy(std::move(result)); diff --git a/cpp/src/sd/nnls.cpp b/cpp/src/sd/nnls.cpp index 4f3264e..67012a6 100644 --- a/cpp/src/sd/nnls.cpp +++ b/cpp/src/sd/nnls.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -290,14 +291,17 @@ RidgeDesign make_ridge_design(const MatrixView& kernel, double ridge) { .normal_equations = normal_equations}; } -Vector solve_ridge_nnls(const RidgeDesign& design, const VectorView& signal, +Vector solve_ridge_nnls(const RidgeDesign& design, const VectorView& signal, double l1, double gradient_tolerance) { + if (l1 > 0 && design.ridge == 0) { + throw std::invalid_argument("l1 > 0 needs ridge > 0."); + } const Index n_measurements = design.kernel.rows(); const Index n_columns = design.kernel.cols(); - // C^T d, with C = [K; ridge I] and d = [s; 0]. - const Vector correlations = design.kernel.transpose() * signal; + // C^T d = K^T s - l1, with C = [K; ridge I] and d = [s; -l1 / ridge]. + const Vector correlations = (design.kernel.transpose() * signal).array() - l1; - // C^T (d - C x) = K^T s - (K^T K + ridge^2 I) x, over the non-zero coefficients only. + // C^T (d - C x) = K^T s - l1 - (K^T K + ridge^2 I) x, over the non-zero coefficients only. const auto gradient_at = [&](const Vector& solution) { Vector gradients = correlations; for (Index column = 0; column < n_columns; ++column) { @@ -323,6 +327,7 @@ Vector solve_ridge_nnls(const RidgeDesign& design, const VectorView& signal, data.head(n_measurements) = signal; if (has_ridge) { active_design.bottomRows(n_active).diagonal().setConstant(design.ridge); + data.tail(n_active).setConstant(-l1 / design.ridge); } return scatter(least_squares(active_design, data), indices, n_columns); }; diff --git a/cpp/src/sd/nnls.hpp b/cpp/src/sd/nnls.hpp index d169c26..0bd5ad1 100644 --- a/cpp/src/sd/nnls.hpp +++ b/cpp/src/sd/nnls.hpp @@ -57,11 +57,12 @@ void remove_cholesky_column(Matrix& factor, Index n_columns, Index position); /// Prepares a kernel for solve_ridge_nnls. RidgeDesign make_ridge_design(const MatrixView& kernel, double ridge); -/// Solves min ||K x - s||^2 + ridge^2 ||x||^2 subject to x >= 0. +/// Solves min ||K x - s||^2 + ridge^2 ||x||^2 + 2 l1 sum(x) subject to x >= 0. /// -/// The steps and result are those of solve_nnls on the stacked matrix [K; ridge I] and data -/// [s; 0]. -Vector solve_ridge_nnls(const RidgeDesign& design, const VectorView& signal, +/// Since x >= 0, the last term is an L1 penalty. The steps and result are those of solve_nnls on +/// the stacked matrix [K; ridge I] and data [s; -l1 / ridge]. l1 > 0 needs ridge > 0, and throws +/// std::invalid_argument otherwise. +Vector solve_ridge_nnls(const RidgeDesign& design, const VectorView& signal, double l1, double gradient_tolerance); } // namespace larmor::sd diff --git a/cpp/src/sd/solver.cpp b/cpp/src/sd/solver.cpp index fbcae42..b5dbea5 100644 --- a/cpp/src/sd/solver.cpp +++ b/cpp/src/sd/solver.cpp @@ -15,7 +15,8 @@ void L2NnlsSolver::solve_block(const MatrixView& block_signals, Eigen::Ref block_fods) const { for (Index voxel = 0; voxel < block_signals.rows(); ++voxel) { const Vector signal = block_signals.row(voxel).transpose(); - block_fods.row(voxel) = solve_ridge_nnls(design, signal, gradient_tolerance).transpose(); + block_fods.row(voxel) = + solve_ridge_nnls(design, signal, 0.0, gradient_tolerance).transpose(); } } diff --git a/cpp/src/simulation/compartments/compartments.cpp b/cpp/src/simulation/compartments/compartments.cpp index 846a4f1..e4ac69d 100644 --- a/cpp/src/simulation/compartments/compartments.cpp +++ b/cpp/src/simulation/compartments/compartments.cpp @@ -5,6 +5,8 @@ #include #include +#include "common/maths.hpp" + namespace larmor::simulation { namespace { @@ -31,13 +33,6 @@ constexpr int N_NEWTON_ROOTS = 200; // Number of Newton steps per cylinder root. constexpr int NEWTON_STEPS = 8; -// y below which exp(-y) I0(y) sums its power series, and above which its asymptotic series. -constexpr double BESSEL_SERIES_LIMIT = 25.0; -// Relative size of the last term at which a series stops. -constexpr double SERIES_TOLERANCE = 1e-17; -// Maximum number of terms of a series. -constexpr int MAX_SERIES_TERMS = 200; - // The Watson integrals in u = 1 - t stop at u = WATSON_CUTOFF / rate, where the integrand is // below exp(-WATSON_CUTOFF). constexpr double WATSON_CUTOFF = 40.0; @@ -289,28 +284,6 @@ double sphere(double b_value, double diffusivity, double radius_um, const PulseT return std::exp(-2.0 * gamma_g_squared * sum); } -double scaled_bessel_i0(double y) { - if (y < BESSEL_SERIES_LIMIT) { - const double quarter_squared = 0.25 * y * y; - double term = 1.0; - double sum = 1.0; - for (int order = 1; order < MAX_SERIES_TERMS && term > SERIES_TOLERANCE * sum; ++order) { - term *= quarter_squared / (static_cast(order) * order); - sum += term; - } - return std::exp(-y) * sum; - } - // 1 / sqrt(2 pi y) sum_k c_k / y^k with c_k = c_{k-1} (2k - 1)^2 / (8 k); the terms fall - // until k is about 2 y, far beyond the precision needed. - double term = 1.0; - double sum = 1.0; - for (int order = 1; order < MAX_SERIES_TERMS && term > SERIES_TOLERANCE * sum; ++order) { - term *= (2.0 * order - 1.0) * (2.0 * order - 1.0) / (8.0 * order * y); - sum += term; - } - return sum / std::sqrt(2.0 * std::numbers::pi * y); -} - double watson_stick(double b_value, double diffusivity, double kappa, double cosine) { const double x = std::max(b_value * diffusivity, 0.0); const double c2 = std::min(cosine * cosine, 1.0); diff --git a/cpp/src/simulation/compartments/compartments.hpp b/cpp/src/simulation/compartments/compartments.hpp index 115920a..456997c 100644 --- a/cpp/src/simulation/compartments/compartments.hpp +++ b/cpp/src/simulation/compartments/compartments.hpp @@ -89,9 +89,6 @@ const std::array& sphere_roots(); /// D = 0. double sphere(double b_value, double diffusivity, double radius_um, const PulseTiming& timing); -/// Returns exp(-y) I0(y) for y >= 0, with I0 the modified Bessel function of the first kind. -double scaled_bessel_i0(double y); - /// Returns the attenuation of sticks with diffusivity D dispersed about the fibre direction mu by /// a Watson density (Zhang et al. 2012, Eqs. 2 and 3): /// diff --git a/cpp/tests/dwi/CMakeLists.txt b/cpp/tests/dwi/CMakeLists.txt index 44ac529..8048930 100644 --- a/cpp/tests/dwi/CMakeLists.txt +++ b/cpp/tests/dwi/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(larmor_tests PRIVATE test_gibbs.cpp test_log_fit.cpp - test_mppca.cpp) + test_mppca.cpp + test_noddi.cpp) diff --git a/cpp/tests/dwi/test_noddi.cpp b/cpp/tests/dwi/test_noddi.cpp new file mode 100644 index 0000000..051a8b1 --- /dev/null +++ b/cpp/tests/dwi/test_noddi.cpp @@ -0,0 +1,207 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include "common/types.hpp" +#include "dwi/noddi.hpp" +#include "simulation/compartments/compartments.hpp" + +namespace { + +using larmor::Index; +using larmor::Matrix; +using larmor::Vector; +using larmor::Vector3; + +constexpr double D_PAR = 1.7e-3; + +// Returns n unit directions spread over the half sphere, by the golden spiral. +Matrix spiral_directions(Index count) { + Matrix directions(count, 3); + const double golden = std::numbers::pi * (3.0 - std::sqrt(5.0)); + for (Index index = 0; index < count; ++index) { + const double z = 1.0 - (static_cast(index) + 0.5) / static_cast(count); + const double radius = std::sqrt(1.0 - z * z); + const double angle = golden * static_cast(index); + directions.row(index) << radius * std::cos(angle), radius * std::sin(angle), z; + } + return directions; +} + +// A two-shell table: 6 b0s, 30 directions at b = 1000 and 60 at b = 2000 s/mm^2. +larmor::dwi::NoddiProtocol two_shells(const larmor::dwi::NoddiSettings& settings) { + const Index n_volumes = 96; + Vector bvals = Vector::Zero(n_volumes); + Matrix bvecs = Matrix::Zero(n_volumes, 3); + std::vector is_b0(static_cast(n_volumes), false); + for (std::size_t volume = 0; volume < 6; ++volume) { + is_b0.at(volume) = true; + } + bvals.segment(6, 30).setConstant(1000.0); + bvals.segment(36, 60).setConstant(2000.0); + bvecs.middleRows(6, 30) = spiral_directions(30); + bvecs.middleRows(36, 60) = spiral_directions(60); + return larmor::dwi::make_noddi_protocol(bvals, bvecs, is_b0, settings); +} + +larmor::dwi::NoddiSettings fit_settings() { + larmor::dwi::NoddiSettings settings; + for (int atom = 0; atom < 12; ++atom) { + settings.v_ic_grid.push_back(0.1 + 0.9 * atom / 11.0); + settings.kappa_grid.push_back(20.0 * atom / 11.0); + } + settings.ridge = std::sqrt(0.001); + settings.l1 = 0.5; + return settings; +} + +} // namespace + +TEST_CASE("noddi: the toolbox's Watson stick matches the exact integral") { + // The series stops at degree 12, so the gap grows with kappa and b; these bounds sit above + // the measured ones (3e-4 at kappa 64 and b = 3000). + const auto settings = fit_settings(); + const Matrix directions = spiral_directions(40); + for (const double bval : {1000.0, 3000.0}) { + Vector bvals = Vector::Constant(40, bval); + const auto protocol = larmor::dwi::make_noddi_protocol( + bvals, directions, std::vector(40, false), settings); + for (const double kappa : {0.0, 0.05, 0.1001, 0.3, 2.0, 16.0, 29.9, 30.1, 64.0}) { + const Vector series = + larmor::dwi::noddi_signal(protocol, settings, 1.0, kappa, 0.0, Vector3::UnitZ()); + for (Index volume = 0; volume < 40; ++volume) { + const double exact = + larmor::simulation::watson_stick(bval, D_PAR, kappa, directions(volume, 2)); + CAPTURE(bval); + CAPTURE(kappa); + CHECK(std::abs(series(volume) - exact) < (kappa > 16.0 ? 1e-3 : 1e-4)); + } + } + } +} + +TEST_CASE("noddi: the whole signal matches Zhang et al.'s model") { + // Eq. 1 with the exact Watson stick, the zeppelin of Eqs. 4 to 6 and tau1 of Eq. 7, the last + // integrated by the midpoint rule. + const auto settings = fit_settings(); + const auto protocol = two_shells(settings); + const Vector3 mu = Vector3(0.3, -0.5, 0.8).normalized(); + for (const double kappa : {0.5, 3.0, 12.0}) { + double weight_sum = 0.0; + double moment = 0.0; + const int n_steps = 100000; + for (int step = 0; step < n_steps; ++step) { + const double cosine = (step + 0.5) / n_steps; + const double weight = std::exp(kappa * cosine * cosine); + weight_sum += weight; + moment += weight * cosine * cosine; + } + const double tau = moment / weight_sum; + const double v_ic = 0.6; + const double v_iso = 0.15; + const double along = D_PAR - D_PAR * v_ic * (1.0 - tau); + const double across = D_PAR - D_PAR * v_ic * (1.0 + tau) / 2.0; + const Vector signal = larmor::dwi::noddi_signal(protocol, settings, v_ic, kappa, v_iso, mu); + for (Index volume = 0; volume < protocol.bvals.size(); ++volume) { + const double bval = protocol.bvals(volume); + const double cosine = protocol.bvecs.row(volume).dot(mu); + const double sticks = larmor::simulation::watson_stick(bval, D_PAR, kappa, cosine); + const double zeppelin = std::exp(-bval * (across + (along - across) * cosine * cosine)); + const double expected = (1.0 - v_iso) * (v_ic * sticks + (1.0 - v_ic) * zeppelin) + + v_iso * std::exp(-bval * 3.0e-3); + CAPTURE(kappa); + CHECK(std::abs(signal(volume) - expected) < 1e-4); + } + } +} + +TEST_CASE("noddi: the Watson coefficients join up where the toolbox changes formula") { + // Just above kappa = 0.1 the closed form of degrees 10 and 12 loses digits to cancellation, as + // in the toolbox; the stick test above bounds what that does to the signal. + for (const auto& [edge, n_terms] : + {std::pair{0.1, std::size_t{5}}, std::pair{30.0, std::size_t{7}}}) { + const auto below = larmor::dwi::watson_sh_coefficients(edge * (1.0 - 1e-9)); + const auto above = larmor::dwi::watson_sh_coefficients(edge * (1.0 + 1e-9)); + for (std::size_t term = 0; term < n_terms; ++term) { + CAPTURE(edge); + CAPTURE(term); + // The fits above kappa = 30 are the toolbox's, good to about 1e-3 of the value. + CHECK(std::abs(below.at(term) - above.at(term)) <= + 2e-3 * std::abs(below.at(term)) + 1e-7); + } + } +} + +TEST_CASE("noddi: the averaged zeppelin follows the tortuosity model") { + // Zhang et al. 2012, Eqs. 5 and 6: d_par' = d - d v (1 - tau1) and d_perp' = d - d v (1 + + // tau1) / 2, with tau1 the mean squared cosine of the Watson density (Eq. 7), integrated here + // by the midpoint rule. + for (const double kappa : {0.0, 1e-6, 0.3, 4.0, 25.0}) { + double weight_sum = 0.0; + double moment = 0.0; + const int n_steps = 200000; + for (int step = 0; step < n_steps; ++step) { + const double cosine = (step + 0.5) / n_steps; + const double weight = std::exp(kappa * cosine * cosine); + weight_sum += weight; + moment += weight * cosine * cosine; + } + const double tau = moment / weight_sum; + for (const double v_ic : {0.2, 0.7}) { + const auto [along, across] = + larmor::dwi::watson_hindered_diffusivities(D_PAR, D_PAR * (1.0 - v_ic), kappa); + CAPTURE(kappa); + CHECK(along == doctest::Approx(D_PAR - D_PAR * v_ic * (1.0 - tau)).epsilon(1e-8)); + CHECK(across == + doctest::Approx(D_PAR - D_PAR * v_ic * (1.0 + tau) / 2.0).epsilon(1e-8)); + } + } +} + +TEST_CASE("noddi: a noise-free voxel's parameters and direction are recovered") { + const auto settings = fit_settings(); + const auto protocol = two_shells(settings); + const Vector3 mu = Vector3(0.3, -0.5, 0.8).normalized(); + // The start's direction is off by a few degrees, as a DTI direction can be. + const Vector3 start = Vector3(0.35, -0.45, 0.82).normalized(); + for (const auto& [v_ic, kappa, v_iso] : + std::vector>{{0.5, 3.0, 0.2}, {0.8, 12.0, 0.05}, {0.3, 0.8, 0.4}}) { + const Vector signal = + 900.0 * larmor::dwi::noddi_signal(protocol, settings, v_ic, kappa, v_iso, mu); + const auto fit = larmor::dwi::fit_noddi_voxel(protocol, settings, signal, start); + CAPTURE(v_ic); + CHECK(fit.fitted); + CHECK(fit.v_ic == doctest::Approx(v_ic).epsilon(1e-3)); + CHECK(fit.kappa == doctest::Approx(kappa).epsilon(1e-2)); + CHECK(fit.v_iso == doctest::Approx(v_iso).epsilon(1e-3)); + CHECK(std::abs(fit.mu.dot(mu)) > std::cos(0.1 * std::numbers::pi / 180.0)); + } +} + +TEST_CASE("noddi: a voxel without a positive b0 is not fitted") { + const auto settings = fit_settings(); + const auto protocol = two_shells(settings); + Vector signal = + 900.0 * larmor::dwi::noddi_signal(protocol, settings, 0.5, 3.0, 0.2, Vector3::UnitZ()); + signal.head(6).setConstant(-1.0); + const auto fit = larmor::dwi::fit_noddi_voxel(protocol, settings, signal, Vector3::UnitZ()); + CHECK_FALSE(fit.fitted); + CHECK(fit.v_ic == 0.0); + CHECK(fit.v_iso == 0.0); +} + +TEST_CASE("noddi: a starting direction that is zero or not finite is not fitted") { + const auto settings = fit_settings(); + const auto protocol = two_shells(settings); + const Vector signal = + 900.0 * larmor::dwi::noddi_signal(protocol, settings, 0.5, 3.0, 0.2, Vector3::UnitZ()); + for (const Vector3& start : {Vector3::Zero().eval(), Vector3(std::nan(""), 0.0, 1.0), + Vector3(std::numeric_limits::infinity(), 0.0, 0.0)}) { + CHECK_FALSE(larmor::dwi::fit_noddi_voxel(protocol, settings, signal, start).fitted); + } +} diff --git a/cpp/tests/sd/test_nnls.cpp b/cpp/tests/sd/test_nnls.cpp index fb4e9ac..cbef692 100644 --- a/cpp/tests/sd/test_nnls.cpp +++ b/cpp/tests/sd/test_nnls.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -225,20 +226,23 @@ TEST_CASE("nnls: the ridge solve matches plain NNLS on the stacked matrix") { } } for (const double ridge : {0.1, 1e-7}) { - const larmor::sd::RidgeDesign design = larmor::sd::make_ridge_design(kernel, ridge); - CHECK(design.normal_equations == (ridge == 0.1)); - Matrix stacked = Matrix::Zero(50, 40); - stacked.topRows(10) = kernel; - stacked.bottomRows(40).diagonal().setConstant(ridge); - for (int voxel = 0; voxel < 5; ++voxel) { - Vector signal(10); - for (Index row = 0; row < 10; ++row) { - signal(row) = 2.0 * random.next() - 0.5; + for (const double l1 : {0.0, 0.01}) { + const larmor::sd::RidgeDesign design = larmor::sd::make_ridge_design(kernel, ridge); + CHECK(design.normal_equations == (ridge == 0.1)); + Matrix stacked = Matrix::Zero(50, 40); + stacked.topRows(10) = kernel; + stacked.bottomRows(40).diagonal().setConstant(ridge); + for (int voxel = 0; voxel < 5; ++voxel) { + Vector signal(10); + for (Index row = 0; row < 10; ++row) { + signal(row) = 2.0 * random.next() - 0.5; + } + Vector data = Vector::Zero(50); + data.head(10) = signal; + data.tail(40).setConstant(-l1 / ridge); + check_close(larmor::sd::solve_ridge_nnls(design, signal, l1, 1e-12), + larmor::sd::solve_nnls(stacked, data, 1e-12), 1e-9); } - Vector data = Vector::Zero(50); - data.head(10) = signal; - check_close(larmor::sd::solve_ridge_nnls(design, signal, 1e-12), - larmor::sd::solve_nnls(stacked, data, 1e-12), 1e-9); } } } @@ -259,7 +263,7 @@ TEST_CASE("nnls: with no ridge, the normal equations give plain NNLS") { for (Index row = 0; row < 30; ++row) { signal(row) = 2.0 * random.next() - 0.5; } - check_close(larmor::sd::solve_ridge_nnls(design, signal, 1e-12), + check_close(larmor::sd::solve_ridge_nnls(design, signal, 0.0, 1e-12), larmor::sd::solve_nnls(kernel, signal, 1e-12), 1e-9); } } @@ -286,3 +290,20 @@ TEST_CASE("nnls: removing a column keeps the Cholesky factor of the others") { CHECK((lower.diagonal().array() > 0.0).all()); } } + +TEST_CASE("nnls: exception thrown if l1 > 0 and ridge == 0") { + PseudoRandom random; + Matrix kernel(10, 40); + for (Index row = 0; row < 10; ++row) { + for (Index column = 0; column < 40; ++column) { + kernel(row, column) = random.next(); + } + } + Vector signal(10); + for (Index row = 0; row < 10; ++row) { + signal(row) = 2.0 * random.next() - 0.5; + } + CHECK_THROWS_AS(larmor::sd::solve_ridge_nnls(larmor::sd::make_ridge_design(kernel, 0.0), signal, + 0.5, 1e-12), + std::invalid_argument); +} diff --git a/cpp/tests/simulation/test_compartments.cpp b/cpp/tests/simulation/test_compartments.cpp index 09bd636..afff0af 100644 --- a/cpp/tests/simulation/test_compartments.cpp +++ b/cpp/tests/simulation/test_compartments.cpp @@ -4,6 +4,7 @@ #include +#include "common/maths.hpp" #include "simulation/compartments/compartments.hpp" namespace { @@ -71,15 +72,11 @@ TEST_CASE("compartments: the cylinder roots solve J1'(x) = 0 and satisfy the sum TEST_CASE("compartments: exp(-y) I0(y) has its reference values on both sides of y = 25") { // The references are scipy.special.i0e. - CHECK(larmor::simulation::scaled_bessel_i0(0.0) == 1.0); - CHECK(larmor::simulation::scaled_bessel_i0(1.0) == - doctest::Approx(0.4657596075936404).epsilon(1e-14)); - CHECK(larmor::simulation::scaled_bessel_i0(24.999) == - doctest::Approx(0.08019839425680447).epsilon(1e-14)); - CHECK(larmor::simulation::scaled_bessel_i0(25.0) == - doctest::Approx(0.08019677354743669).epsilon(1e-14)); - CHECK(larmor::simulation::scaled_bessel_i0(40.0) == - doctest::Approx(0.06327827987523532).epsilon(1e-14)); + CHECK(larmor::scaled_bessel_i0(0.0) == 1.0); + CHECK(larmor::scaled_bessel_i0(1.0) == doctest::Approx(0.4657596075936404).epsilon(1e-14)); + CHECK(larmor::scaled_bessel_i0(24.999) == doctest::Approx(0.08019839425680447).epsilon(1e-14)); + CHECK(larmor::scaled_bessel_i0(25.0) == doctest::Approx(0.08019677354743669).epsilon(1e-14)); + CHECK(larmor::scaled_bessel_i0(40.0) == doctest::Approx(0.06327827987523532).epsilon(1e-14)); } TEST_CASE("compartments: the dot is 1 and the Watson sticks span the powder stick to the stick") { diff --git a/docs/_extra/llms.txt b/docs/_extra/llms.txt index 92010f3..ca19274 100644 --- a/docs/_extra/llms.txt +++ b/docs/_extra/llms.txt @@ -19,6 +19,7 @@ Larmor simulates and processes diffusion MRI from Python and from its command li - [Gradient tables and shells](https://larmor.org/guides/gradient-tables/): bvals and bvecs, MRtrix gradient files, the FSL and world frames, and shells - [DTI](https://larmor.org/guides/dti/): fitting the diffusion tensor, and the FA, MD, AD, RD and direction maps - [DKI](https://larmor.org/guides/dki/): fitting diffusion kurtosis, constrained and weighted, and the MK, AK, RK and KA maps +- [NODDI](https://larmor.org/guides/noddi/): fitting NODDI, and the neurite density, orientation dispersion and free-water maps - [FODs](https://larmor.org/guides/fods/): response functions, the GRL, mFOD, CSD and MSMT-CSD methods and their solvers, and SH bases - [Peaks](https://larmor.org/guides/peaks/): the peaks of FODs, and the peaks image - [Tractography](https://larmor.org/guides/tractography/): deterministic and probabilistic tracking, seeding, stopping and output formats diff --git a/docs/about/citing.md b/docs/about/citing.md index bc0fa4f..8eb9391 100644 --- a/docs/about/citing.md +++ b/docs/about/citing.md @@ -51,6 +51,13 @@ The same information is in the repository's `CITATION.cff` file, which GitHub an - **Weighted fits of DTI and DKI.** Veraart J, Sijbers J, Sunaert S, Leemans A, Jeurissen B. Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls. NeuroImage 81 (2013) 335-346. +- **NODDI.** Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite + orientation dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016. +- **The start of the NODDI fit,** a dictionary solved by convex optimisation. Daducci A, Canales-Rodriguez EJ, + Zhang H, Dyrby TB, Alexander DC, Thiran JP. Accelerated Microstructure Imaging via Convex Optimization (AMICO) + from diffusion MRI data. NeuroImage 105 (2015) 32-44. +- **The NODDI model and fit as Larmor computes them** follow the NODDI Matlab Toolbox 1.05, + [nitrc.org/projects/noddi_toolbox](https://www.nitrc.org/projects/noddi_toolbox). ### Spherical deconvolution diff --git a/docs/cli.md b/docs/cli.md index f78d89a..751e0bf 100644 --- a/docs/cli.md +++ b/docs/cli.md @@ -19,6 +19,7 @@ $ lar [OPTIONS] COMMAND [ARGS]... * `dti`: Fit the diffusion tensor (DTI) in every... * `dki`: Fit diffusion kurtosis imaging (DKI) to... +* `noddi`: Fit NODDI (neurite density and orientation... * `denoise`: Denoise a DWI volume with MP-PCA and... * `mask`: Compute a binary brain mask from the mean... * `degibbs`: Remove Gibbs ringing from an image by... @@ -95,6 +96,34 @@ $ lar dki [OPTIONS] {volume} {output} * `--debug`: Show debug messages and full tracebacks. * `--help`: Show this message and exit. +## `lar noddi` + +Fit NODDI (neurite density and orientation dispersion) and write its maps. + +**Usage**: + +```console +$ lar noddi [OPTIONS] {volume} {output} +``` + +**Arguments**: + +* `volume`: Path to the diffusion volume (.nii, .nii.gz, .mif, .mif.gz or .mih). [required] +* `output`: Prefix of the output maps. 'out' writes the neurite density to out_v_ic.nii.gz, the orientation dispersion index to out_od, the Watson concentration to out_kappa, the free-water fraction to out_v_iso and the mean neurite direction to out_mu, all as .nii.gz. The direction is in world coordinates (RAS), not in the voxel frame of the image. [required] + +**Options**: + +* `--bvals `: Path to the b-values file. It is not needed for a .mif file with a gradient table in its header; given with such a file, --bvals and --bvecs replace that table. +* `--bvecs `: Path to the b-vectors file. It is not needed for a .mif file with a gradient table in its header; given with such a file, --bvals and --bvecs replace that table. +* `--bvec-frame `: Coordinate frame of the b-vectors. 'fsl' means the frame of FSL bvecs files, along the image axes with x flipped when the affine determinant is positive, and 'world' means scanner coordinates (RAS). [default: fsl] +* `-m, --mask `: Brain mask; voxels outside it are not fitted. Without it, a mask is computed from the mean b = 0 image. +* `--force`: Overwrite existing output files. Without it, the command stops before any work when an output exists. +* `--nt `: Number of threads; -1 uses all CPU cores. [default: -1] +* `--quiet`: Show only errors. +* `--info`: Show progress and informational messages. +* `--debug`: Show debug messages and full tracebacks. +* `--help`: Show this message and exit. + ## `lar denoise` Denoise a DWI volume with MP-PCA and estimate its noise level. diff --git a/docs/guides/index.md b/docs/guides/index.md index 7a7d633..1ebf91d 100644 --- a/docs/guides/index.md +++ b/docs/guides/index.md @@ -19,6 +19,7 @@ The task pages show each step both from Python and from the command line. - [Gradient tables and shells](gradient-tables.md) - [DTI](dti.md) - [DKI](dki.md) +- [NODDI](noddi.md) - [FODs](fods.md) - [Peaks](peaks.md) - [Tractography](tractography.md) @@ -37,6 +38,7 @@ Loading and saving data Gradient tables and shells DTI DKI +NODDI FODs Peaks Tractography diff --git a/docs/guides/noddi.md b/docs/guides/noddi.md new file mode 100644 index 0000000..97704c2 --- /dev/null +++ b/docs/guides/noddi.md @@ -0,0 +1,91 @@ +--- +myst: + html_meta: + description: "How to fit NODDI with Larmor, from Python and with lar noddi, to get maps of neurite density, orientation dispersion and free water." +--- + +# NODDI + +Neurite orientation dispersion and density imaging (NODDI) separates two things that both lower the fractional +anisotropy of DTI: how many neurites a voxel holds, and how widely their directions spread. It models the signal +of each voxel with three kinds of water (Zhang et al. 2012). Water inside axons and dendrites moves only along +them, as sticks whose directions follow a Watson distribution about a mean direction; water between them is +hindered, more along that direction than across it; and free water, such as cerebrospinal fluid, diffuses +equally in every direction. The fit gives the neurite density v_ic, the orientation dispersion index OD, from +0 for aligned neurites to 1 for neurites in every direction, and the free-water fraction v_iso. NODDI needs at +least two non-zero shells. [Citing](../about/citing.md#models) lists the papers. + +```{include} ../_snippets/synthetic-dwi.md +``` + +## Fitting the model + +::::{tab-set} +:sync-group: interface + +:::{tab-item} Python +:sync: python +```python +from larmor.dwi import NoddiModel +from larmor.io import load_dwi, load_volume, save_volume + +dwi, gtab = load_dwi("dwi.nii.gz", "dwi.bval", "dwi.bvec") +fit = NoddiModel().fit(dwi, gtab, mask=load_volume("mask.nii.gz")) + +save_volume(dwi.with_data(fit.od), "od.nii.gz") +print(f"OD in a bundle: {fit.od[0, 5, 0]:.2f}, where they cross: {fit.od[5, 5, 0]:.2f}") +``` + +The result, a {py:class}`~larmor.dwi.NoddiFit`, holds the neurite density `v_ic`, the orientation dispersion +index `od`, the Watson concentration `kappa` it comes from, the free-water fraction `v_iso` and the fitted mean +direction `mu`. All maps are zero outside the mask. +::: + +:::{tab-item} Command line +:sync: cli +```bash +lar noddi dwi.nii.gz noddi --bvals dwi.bval --bvecs dwi.bvec -m mask.nii.gz +``` + +The command writes `noddi_v_ic.nii.gz`, `noddi_od.nii.gz`, `noddi_kappa.nii.gz`, `noddi_v_iso.nii.gz` and the +mean direction `noddi_mu.nii.gz`. +::: +:::: + +## How it is fitted + +Each voxel is fitted in two steps. The first matches its signal with a non-negative mix of NODDI signals computed +in advance for a grid of neurite densities and concentrations about the voxel's DTI direction, which gives a +starting point quickly and without a search (Daducci et al. 2015). The second fits v_ic, kappa, v_iso and the +direction from there by maximum likelihood with Rician noise, with the model, the noise level from the b=0 +images, and the two stages of the NODDI Matlab Toolbox: first the fractions alone, then all the parameters. The +start and the optimiser differ from the toolbox's, so a few voxels can reach another optimum. The diffusivities are fixed at the toolbox's in vivo values, 1.7e-3 mm^2/s along the neurites and 3.0e-3 mm^2/s for +free water. + +## Options and caveats + +- **Signal fractions.** v_ic and v_iso are fractions of the signal, not of the volume. They differ when the + compartments have different T2 or proton density, as free water does. +- **Voxels of free water.** Where v_iso is high, little signal is left to measure the neurites, and v_ic is not + reliable there; it often reaches 1. Read v_ic and OD where v_iso is low. +- **Voxels that are not fitted.** A voxel with a signal that is not finite (NaN or infinite), without a + positive b=0 signal, or with fewer than six positive diffusion-weighted measurements gets zeros, and a warning + gives their number. A voxel whose signal the first step explains with free water alone + gets v_ic = 0 and OD = 0. +- **Starting direction.** By default the fit starts from the principal eigenvector of a DTI fit. From Python, + `fit(..., mu=directions)` starts from directions of your own, one vector per voxel; they are normalised, and + one that is zero or not finite in the mask is refused. +- **Mask and threads.** As for [DTI](dti.md#options). + +## See also + +- API: {py:class}`~larmor.dwi.NoddiModel`, {py:class}`~larmor.dwi.NoddiFit`. +- Command: [`lar noddi`](../cli.md#lar-noddi). + +## References + +- Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation + dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016. +- Daducci A, Canales-Rodriguez EJ, Zhang H, Dyrby TB, Alexander DC, Thiran JP. Accelerated Microstructure Imaging + via Convex Optimization (AMICO) from diffusion MRI data. NeuroImage 105 (2015) 32-44. +- NODDI Matlab Toolbox 1.05, [nitrc.org/projects/noddi_toolbox](https://www.nitrc.org/projects/noddi_toolbox). diff --git a/docs/guides/what-larmor-does.md b/docs/guides/what-larmor-does.md index 11f534b..89dacef 100644 --- a/docs/guides/what-larmor-does.md +++ b/docs/guides/what-larmor-does.md @@ -22,6 +22,7 @@ under [Citing](../about/citing.md#methods-and-papers). ## Models - DTI and DKI, fitted by ordinary or weighted least squares (`lar dti`, `lar dki`). +- NODDI, for neurite density, orientation dispersion and free water (`lar noddi`). - Fibre orientation distributions (FODs) by spherical deconvolution (`lar fod`), with generalized Richardson-Lucy (GRL), with mFOD, which fits one FOD for each anisotropic tissue, or with constrained spherical deconvolution (CSD). GRL and mFOD can be solved by damped Richardson-Lucy (dRL), non-negative least squares or diff --git a/docs/reference/api/dwi.rst b/docs/reference/api/dwi.rst index 0b76d59..22d9e94 100644 --- a/docs/reference/api/dwi.rst +++ b/docs/reference/api/dwi.rst @@ -1,5 +1,5 @@ .. meta:: - :description: Diffusion-weighted images in Larmor: brain masks, MP-PCA denoising, Gibbs ringing removal, and the DTI and DKI models. + :description: Diffusion-weighted images in Larmor: brain masks, MP-PCA denoising, Gibbs ringing removal, and the DTI, DKI and NODDI models. larmor.dwi ========== diff --git a/pyproject.toml b/pyproject.toml index da5fb9a..d1383e2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -270,6 +270,7 @@ layers = ["main", "cli", "viewer", "plot", "simulation", "tractography", "sd", " # `from larmor import _core` reads as an import of the package: the compiled module is invisible here. ignore_imports = [ "larmor.dwi.fitting -> larmor", + "larmor.dwi.noddi -> larmor", "larmor.dwi.denoise -> larmor", "larmor.dwi.gibbs -> larmor", "larmor.simulation.instant -> larmor", diff --git a/src/larmor/_core.pyi b/src/larmor/_core.pyi index 6aee759..a6e0205 100644 --- a/src/larmor/_core.pyi +++ b/src/larmor/_core.pyi @@ -164,6 +164,20 @@ def lsqnonneg(c: Annotated[NDArray[numpy.float64], dict(shape=(None, None), devi ValueError: If c is not C-contiguous or the shapes do not match. """ +def ridge_l1_nnls(c: Annotated[NDArray[numpy.float64], dict(shape=(None, None), device='cpu', writable=False)], d: Annotated[NDArray[numpy.float64], dict(shape=(None,), device='cpu', writable=False)], ridge: float, l1: float, tol_x: float | None = None) -> Annotated[NDArray[numpy.float64], dict(shape=(None,))]: + """ + Solves min ||c x - d||^2 + ridge^2 ||x||^2 + 2 l1 sum(x) with x >= 0. + + c is (m, n) and d (m,). Since x >= 0, the last term is an L1 penalty, which + sets the coefficients the data do not support to zero; ridge = l1 = 0 is + plain NNLS. A coefficient enters the solution only when its gradient + exceeds tol_x; None uses 10 * eps * norm(c, 1) * max(m, n). + + Raises: + ValueError: If c is not C-contiguous, the shapes do not match, ridge or + l1 is negative or not finite, or l1 > 0 with ridge = 0. + """ + class DeconvKernel: """One resolution of a multi-tissue kernel, copied into C++.""" @@ -235,6 +249,35 @@ def fit_log_signal(x: Annotated[NDArray[numpy.float64], dict(shape=(None, None), max_iter is negative or n_threads is 0. """ +def fit_noddi(signals: Annotated[NDArray[numpy.float64], dict(shape=(None, None), device='cpu', writable=False)], mu: Annotated[NDArray[numpy.float64], dict(shape=(None, None), device='cpu', writable=False)], bvals: Annotated[NDArray[numpy.float64], dict(shape=(None,), device='cpu', writable=False)], bvecs: Annotated[NDArray[numpy.float64], dict(shape=(None, None), device='cpu', writable=False)], is_b0: Annotated[NDArray[numpy.bool_], dict(shape=(None,), device='cpu', writable=False)], v_ic_grid: Annotated[NDArray[numpy.float64], dict(shape=(None,), device='cpu', writable=False)], kappa_grid: Annotated[NDArray[numpy.float64], dict(shape=(None,), device='cpu', writable=False)], d_par: float, d_iso: float, ridge: float, l1: float, n_threads: int) -> tuple[Annotated[NDArray[numpy.float64], dict(shape=(None, None))], Annotated[NDArray[numpy.bool_], dict(shape=(None,))]]: + """ + Fits NODDI to each voxel's signal, a row of signals (v, n), from the direction + mu (v, 3). + + The start comes from a dictionary of NODDI signals over v_ic_grid and + kappa_grid, solved with the ridge and l1 weights; the fit is then by maximum + likelihood with Rician noise. + d_par and d_iso are in mm^2/s and bvals in s/mm^2. + + Returns v_ic, kappa, v_iso and the fitted direction of each voxel (v, 6), and + whether each voxel was fitted (v,). + + Raises: + ValueError: If an array is not C-contiguous, the shapes do not match, a + grid is empty, ridge is not positive, l1 is negative or n_threads is 0. + """ + +def noddi_signal(bvals: Annotated[NDArray[numpy.float64], dict(shape=(None,), device='cpu', writable=False)], bvecs: Annotated[NDArray[numpy.float64], dict(shape=(None, None), device='cpu', writable=False)], v_ic: float, kappa: float, v_iso: float, mu: Sequence[float], d_par: float, d_iso: float) -> Annotated[NDArray[numpy.float64], dict(shape=(None,))]: + """ + Returns the normalised NODDI signal of each volume for the neurite density + v_ic, Watson concentration kappa and free-water fraction v_iso about the + direction mu. bvals (n,) are in s/mm^2, bvecs (n, 3) unit directions, and + d_par and d_iso in mm^2/s. + + Raises: + ValueError: If an array is not C-contiguous or the shapes do not match. + """ + def denoise_mppca(image: Annotated[NDArray[numpy.float32], dict(shape=(None, None, None, None), device='cpu', writable=False)], voxels: Annotated[NDArray[numpy.int64], dict(shape=(None,), device='cpu', writable=False)], patch: Sequence[int], n_threads: int) -> tuple[Annotated[NDArray[numpy.float64], dict(shape=(None, None))], Annotated[NDArray[numpy.float64], dict(shape=(None,))], Annotated[NDArray[numpy.int64], dict(shape=(None,))]]: """ Denoises by MP-PCA the voxels at flat indices voxels (k,) int64 of a diff --git a/src/larmor/cli/noddi.py b/src/larmor/cli/noddi.py new file mode 100644 index 0000000..60298fd --- /dev/null +++ b/src/larmor/cli/noddi.py @@ -0,0 +1,77 @@ +"""The `lar noddi` command, which fits NODDI to multi-shell data and writes its maps.""" + +import logging +import typing + +import typer + +from larmor._log import LarmorLogger +from larmor.cli.options import ( + BValsOpt, + BvecFrame, + BvecFrameOpt, + BVecsOpt, + DiffusionArg, + ForceOpt, + ThreadsOpt, +) +from larmor.cli.outputs import check_outputs, output_prefix, with_suffix +from larmor.io import load_dwi, load_volume, save_volume + +MaskArg = typing.Annotated[ + str | None, + typer.Option( + "--mask", + "-m", + help="Brain mask; voxels outside it are not fitted. Without it, a mask is " + "computed from the mean b = 0 image.", + ), +] + +OutputArg = typing.Annotated[ + str, + typer.Argument( + help="Prefix of the output maps. 'out' writes the neurite density to " + "out_v_ic.nii.gz, the orientation dispersion index to out_od, the Watson " + "concentration to out_kappa, the free-water fraction to out_v_iso and the " + "mean neurite direction to out_mu, all as .nii.gz. The direction is in " + "world coordinates (RAS), not in the voxel frame of the image.", + ), +] + + +# The NoddiFit attributes saved, each to _.nii.gz. +METRICS = ("v_ic", "od", "kappa", "v_iso", "mu") + + +def noddi( + volume: DiffusionArg, + output: OutputArg, + bvals: BValsOpt = None, + bvecs: BVecsOpt = None, + bvec_frame: BvecFrameOpt = BvecFrame.fsl, + mask: MaskArg = None, + force: ForceOpt = False, + jobs: ThreadsOpt = -1, +) -> None: + """Fit NODDI (neurite density and orientation dispersion) and write its maps.""" + # Imported here, so that starting lar and Tab completion do not load scipy. + from larmor.dwi import NoddiModel + + prefix = output_prefix(output) + outputs = {name: with_suffix(prefix, f"_{name}.nii.gz") for name in METRICS} + check_outputs(list(outputs.values()), force) + + image, gtab = load_dwi(volume, bvals=bvals, bvecs=bvecs, frame=bvec_frame.value) + mask_volume = load_volume(mask) if mask is not None else None + noddi_fit = NoddiModel().fit( + image, + gtab, + mask_volume, + n_jobs=jobs, + progress=LarmorLogger.isEnabledFor(logging.WARNING), + ) + + for name, path in outputs.items(): + save_volume(image.with_data(getattr(noddi_fit, name)), path) + LarmorLogger.info(f"Saved {name.upper()} to {path}") diff --git a/src/larmor/dwi/__init__.py b/src/larmor/dwi/__init__.py index a1bea7c..d52a9e9 100644 --- a/src/larmor/dwi/__init__.py +++ b/src/larmor/dwi/__init__.py @@ -1,4 +1,4 @@ -"""Preparing diffusion-weighted images (DWIs) and fitting the DTI and DKI models.""" +"""Preparing diffusion-weighted images (DWIs) and fitting DTI, DKI and NODDI.""" from larmor.dwi.denoise import DenoisedDWI, denoise_mppca from larmor.dwi.dki import DKIFit, DKIModel @@ -15,6 +15,7 @@ largest_component, threshold_mask, ) +from larmor.dwi.noddi import NoddiFit, NoddiModel from larmor.dwi.shells import shell_mask __all__ = [ @@ -23,6 +24,8 @@ "DTIFit", "DTIModel", "DenoisedDWI", + "NoddiFit", + "NoddiModel", "combine_masks", "create_brain_mask", "denoise_mppca", diff --git a/src/larmor/dwi/dki.py b/src/larmor/dwi/dki.py index 0b044c2..3301d74 100644 --- a/src/larmor/dwi/dki.py +++ b/src/larmor/dwi/dki.py @@ -15,7 +15,9 @@ check_dwi, estimate_s0, floor_bad_samples, + largest_finite_signal, mean_b0, + signal_chunks, ) from larmor.dwi.mask import create_brain_mask, unmask from larmor.types import GradientTable, Volume @@ -216,7 +218,6 @@ def fit( else: mask = as_mask(volume, mask) - s = volume.as_float()[mask].astype(np.float64) fitted = mask.copy() is_b0 = gtab.b0_mask s0: np.ndarray | None = None @@ -230,17 +231,26 @@ def fit( int((~usable).sum()), ) fitted[mask] = usable - s, s0 = s[usable], s0[usable] - own = s[:, is_b0].mean(axis=1) - repaired = ~np.isclose(s0, own, rtol=1e-9, atol=0.0) - s[np.ix_(repaired, is_b0)] = s0[repaired, None] - s[:, ~is_b0] = np.minimum(s[:, ~is_b0], s0[:, None]) - floor_bad_samples(s, s0) + s0 = s0[usable] + voxels = np.nonzero(fitted) + # Without b0, every voxel's floor comes from the largest signal of all the + # voxels, not of a chunk. + largest = largest_finite_signal(volume, voxels) if s0 is None else None constraints = _constraint_matrix(gtab) if self.constrained else None - params, broke = fit_log_signal( - design, s, self.fit_method == "wls", constraints, threads - ) + params = np.empty((design.shape[1], len(voxels[0]))) + broke = np.empty(len(voxels[0]), dtype=bool) + for rows, s in signal_chunks(volume, voxels): + chunk_s0 = None if s0 is None else s0[rows] + if chunk_s0 is not None: + own = s[:, is_b0].mean(axis=1) + repaired = ~np.isclose(chunk_s0, own, rtol=1e-9, atol=0.0) + s[np.ix_(repaired, is_b0)] = chunk_s0[repaired, None] + s[:, ~is_b0] = np.minimum(s[:, ~is_b0], chunk_s0[:, None]) + floor_bad_samples(s, chunk_s0, largest=largest) + params[:, rows], broke[rows] = fit_log_signal( + design, s, self.fit_method == "wls", constraints, threads + ) if self.constrained: LarmorLogger.info( "%d of %d voxels break a constraint.", diff --git a/src/larmor/dwi/dti.py b/src/larmor/dwi/dti.py index 6bb23f5..5f4ad9c 100644 --- a/src/larmor/dwi/dti.py +++ b/src/larmor/dwi/dti.py @@ -12,7 +12,9 @@ check_dwi, estimate_s0, floor_bad_samples, + largest_finite_signal, mean_b0, + signal_chunks, ) from larmor.dwi.mask import create_brain_mask, unmask from larmor.types import GradientTable, Volume @@ -125,7 +127,6 @@ def fit( """ threads = n_threads(n_jobs) check_dwi(volume, gtab) - s = volume.as_float() b_mat = gtab.b_matrix # The fit itself needs no b0 (log S0 is its intercept); only the automatic @@ -135,11 +136,15 @@ def fit( else: mask = as_mask(volume, mask) - s_masked = s[mask].astype(np.float64) s0 = estimate_s0(volume, gtab, mask)[mask] if gtab.b0_mask.any() else None - floor_bad_samples(s_masked, s0) - - N = s.shape[-1] + voxels = np.nonzero(mask) + # The floor of voxels without a usable S0 comes from the largest signal of + # all the voxels, not of a chunk. + largest = None + if s0 is None or not (np.isfinite(s0) & (s0 > 0)).all(): + largest = largest_finite_signal(volume, voxels) + + N = volume.data.shape[-1] design_matrix = np.column_stack([np.ones(N), -b_mat]) # A rank-deficient system has no unique fit (e.g. one shell without b0). if np.linalg.matrix_rank(design_matrix) < design_matrix.shape[1]: @@ -147,9 +152,12 @@ def fit( "The b-values and b-vectors cannot determine a tensor: DTI needs at " "least two distinct b-values and six non-coplanar directions." ) - x, _ = fit_log_signal( - design_matrix, s_masked, self.fit_method == "wls", n_threads=threads - ) + x = np.empty((design_matrix.shape[1], len(voxels[0]))) + for rows, s in signal_chunks(volume, voxels): + floor_bad_samples(s, None if s0 is None else s0[rows], largest=largest) + x[:, rows], _ = fit_log_signal( + design_matrix, s, self.fit_method == "wls", n_threads=threads + ) d_params = x[1:7, :] dt, fa, md, ad, rd, pe, eigen_vals, npd = compute_dti_metrics(d_params, mask) diff --git a/src/larmor/dwi/inputs.py b/src/larmor/dwi/inputs.py index b282a2d..a312d0c 100644 --- a/src/larmor/dwi/inputs.py +++ b/src/larmor/dwi/inputs.py @@ -1,4 +1,6 @@ -"""Checks the inputs of the DWI fits, and computes the mean b0 and S0.""" +"""Checks the inputs of the DWI fits, reads their voxels, and computes b0 and S0.""" + +from collections.abc import Iterator import numpy as np from scipy.ndimage import uniform_filter @@ -11,6 +13,8 @@ S0_NEIGHBOURHOODS = (3, 5) # Fraction of S0 that replaces a measurement that is not finite and positive. SIGNAL_FLOOR = 1e-3 +# Voxels read at a time by `signal_chunks`. +READ_CHUNK = 2**15 def check_dwi(volume: Volume, gtab: GradientTable) -> None: @@ -130,21 +134,117 @@ def estimate_s0(volume: Volume, gtab: GradientTable, mask: np.ndarray) -> np.nda return s0 -def floor_bad_samples(s: np.ndarray, s0: np.ndarray | None) -> None: +def signal_rows( + volume: Volume, + voxels: tuple[np.ndarray, ...], + rows: slice, + columns: np.ndarray | None = None, +) -> np.ndarray: + """Returns the scaled signal of some voxels as float64, one row per voxel. + + The values are scaled in float32, as `Volume.as_float` does, then made + float64, so a fit reading them in chunks sees the same numbers as one reading + the whole image. + + Args: + volume: The DWI, shape (X, Y, Z, N). + voxels: The voxel indices, as `np.nonzero` of a mask gives them. + rows: Which of those voxels to read. + columns: The volumes to keep, a boolean mask of length N. Defaults to + None, every volume. + + Returns: + The signals, shape (n_rows, N), or (n_rows, n_columns) with `columns`. + """ + index = tuple(axis[rows] for axis in voxels) + signal = volume.data[index] + if columns is not None: + signal = signal[:, columns] + # The steps of Volume.as_float, in float32. + scaled = signal.astype(np.float32) + if volume.slope != 1.0 or volume.intercept != 0.0: + scaled *= volume.slope + scaled += volume.intercept + result: np.ndarray = scaled.astype(np.float64) + return result + + +def signal_chunks( + volume: Volume, + voxels: tuple[np.ndarray, ...], + columns: np.ndarray | None = None, + size: int | None = None, +) -> Iterator[tuple[slice, np.ndarray]]: + """Yields the voxels' signals chunk by chunk, with the rows each chunk covers. + + Args: + volume: The DWI, shape (X, Y, Z, N). + voxels: The voxel indices, as `np.nonzero` of a mask gives them. + columns: The volumes to keep; see `signal_rows`. Defaults to None. + size: Voxels per chunk. Defaults to None, READ_CHUNK. + + Yields: + The rows of the chunk, and their signals as `signal_rows` returns them. + """ + size = READ_CHUNK if size is None else size + n_voxels = len(voxels[0]) + for start in range(0, n_voxels, size): + rows = slice(start, min(start + size, n_voxels)) + yield rows, signal_rows(volume, voxels, rows, columns) + + +def finite_rows( + volume: Volume, + voxels: tuple[np.ndarray, ...], + columns: np.ndarray | None = None, +) -> np.ndarray: + """Returns whether each voxel's signal is finite in every kept volume. + + Args: + volume: The DWI, shape (X, Y, Z, N). + voxels: The voxel indices, as `np.nonzero` of a mask gives them. + columns: The volumes to check; see `signal_rows`. Defaults to None. + + Returns: + One flag per voxel, bool. + """ + flags = [ + np.isfinite(s).all(axis=1) for _, s in signal_chunks(volume, voxels, columns) + ] + return np.concatenate(flags) if flags else np.zeros(0, dtype=bool) + + +def largest_finite_signal(volume: Volume, voxels: tuple[np.ndarray, ...]) -> float: + """Returns the largest finite signal of the voxels, or 1 when none is positive.""" + largest = -np.inf + for _, s in signal_chunks(volume, voxels): + finite = s[np.isfinite(s)] + if finite.size: + largest = max(largest, float(finite.max())) + return largest if largest > 0 else 1.0 + + +def floor_bad_samples( + s: np.ndarray, s0: np.ndarray | None, *, largest: float | None = None +) -> None: """Sets the measurements that are not finite and positive to 1e-3 S0, in place. A voxel without a finite, positive S0 (or every voxel when `s0` is None, as - for a table without b0) uses 1e-3 times the largest finite value of `s`. + for a table without b0) uses 1e-3 times the largest finite value of `s`, or + of `largest` when given, for a fit that floors its voxels in chunks. Args: s: The measurements, shape (n_voxels, N). s0: The S0 of each voxel, shape (n_voxels,), or None. + largest: The largest finite signal of all the voxels being fitted. + Defaults to None, the largest of `s`. """ bad = ~(np.isfinite(s) & (s > 0)) if not bad.any(): return - finite = s[np.isfinite(s)] - largest = float(finite.max()) if finite.size and finite.max() > 0 else 1.0 + if largest is None: + finite = s[np.isfinite(s)] + largest = float(finite.max()) if finite.size and finite.max() > 0 else 1.0 floor = np.full(len(s), SIGNAL_FLOOR * largest) if s0 is not None: usable = np.isfinite(s0) & (s0 > 0) diff --git a/src/larmor/dwi/noddi.py b/src/larmor/dwi/noddi.py new file mode 100644 index 0000000..7fd4254 --- /dev/null +++ b/src/larmor/dwi/noddi.py @@ -0,0 +1,182 @@ +"""The NODDI model, started from a signal dictionary and fitted by Rician likelihood.""" + +from dataclasses import dataclass +from math import sqrt + +import numpy as np +from tqdm import tqdm + +from larmor import _core +from larmor._log import LarmorLogger +from larmor._threads import n_threads +from larmor.dwi.dti import DTIModel +from larmor.dwi.inputs import as_mask, check_dwi, mean_b0, signal_rows +from larmor.dwi.mask import create_brain_mask +from larmor.types import GradientTable, Volume + +# The neurite densities and Watson concentrations of the dictionary's atoms. +V_IC = np.linspace(0.1, 1, 12) +KAPPA = np.linspace(0, 20, 12) +# The diffusivities NODDI fixes, along the neurites and of free water, mm^2/s. +D_PAR = 1.7e-3 +D_ISO = 3.0e-3 +# Weights of the ridge (L2) and L1 penalties of the dictionary's sparse solve. +LAMBDA = 0.001 +GAMMA = 0.5 +# Voxels per call to the compiled fit. +CHUNK_VOXELS = 4096 + + +@dataclass(frozen=True) +class NoddiFit: + """Holds the maps of a NODDI fit, on the image grid and zero outside the mask. + + Each map has the spatial shape (X, Y, Z) of the DWI, with a last axis of 3 + for mu, and is float32. The fractions are signal fractions, which differ from + volume fractions when the compartments have different T2 or proton density. + + Attributes: + v_ic: Neurite density, the neurites' share of the tissue outside free + water. + od: Orientation dispersion index, (2 / pi) arctan(1 / kappa), from 0 + (aligned neurites) to 1 (isotropic); 0 where v_ic is 0. + kappa: Watson concentration of the neurites about mu. + v_iso: Free-water fraction. + mu: Fitted mean neurite direction in world coordinates (RAS+); its sign + is arbitrary. + """ + + v_ic: np.ndarray + od: np.ndarray + kappa: np.ndarray + v_iso: np.ndarray + mu: np.ndarray + + +class NoddiModel: + """NODDI model, which fits neurite density, dispersion and free water per voxel. + + `fit` returns a `NoddiFit`. The fit is that of the NODDI Matlab Toolbox, by + maximum likelihood with Rician noise, but it starts from a dictionary solve + (Daducci et al. 2015) rather than a grid search. + + References: + Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical + in vivo neurite orientation dispersion and density imaging of the human + brain. NeuroImage 61 (2012) 1000-1016. + Daducci A, Canales-Rodriguez EJ, Zhang H, Dyrby TB, Alexander DC, Thiran + JP. Accelerated Microstructure Imaging via Convex Optimization (AMICO) + from diffusion MRI data. NeuroImage 105 (2015) 32-44. + NODDI Matlab Toolbox 1.05, https://www.nitrc.org/projects/noddi_toolbox + """ + + def fit( + self, + volume: Volume, + gtab: GradientTable, + mask: Volume | np.ndarray | None = None, + mu: np.ndarray | None = None, + *, + n_jobs: int = -1, + progress: bool = False, + ) -> NoddiFit: + """Fits NODDI in every voxel of the mask. + + Each voxel is fitted on its own, so the result does not depend on n_jobs. + A voxel with a signal that is not finite, without a positive mean b0, or + with fewer than six positive diffusion-weighted measurements is not fitted + and gets zeros, and a warning gives the number of such voxels. + + Args: + volume: Diffusion-weighted volume, shape (X, Y, Z, N), with b0 volumes. + gtab: Its gradient table, with at least two non-zero shells. + mask: Voxels to fit (nonzero inside), on the volume's grid. Defaults + to None, which computes a brain mask from the mean b0 volume. + mu: The starting direction of each voxel, (X, Y, Z, 3), in world + coordinates; it need not be unit length. Defaults to None, which + takes the principal eigenvector of a DTI fit. + n_jobs: Number of threads for the fit, at most the number of cores; + -1 uses every core. Defaults to -1. + progress: Whether to show a progress bar. Defaults to False. + + Returns: + NoddiFit: The maps, zero outside the mask. + + Raises: + ValueError: If the table or the mask does not match the volume, the + table has fewer than two non-zero shells, mu does not have the + shape (X, Y, Z, 3) or is zero or not finite in a mask voxel, or + n_jobs is 0 or below -1. + """ + threads = n_threads(n_jobs) + check_dwi(volume, gtab) + if gtab.count_shells() < 2: + raise ValueError( + "NODDI needs at least two non-zero b-value shells; the gradient table " + f"has {gtab.count_shells()}." + ) + shape = volume.data.shape[:3] + if mu is not None and mu.shape != (*shape, 3): + raise ValueError( + f"mu must have the shape {(*shape, 3)} of the volume, got {mu.shape}." + ) + mask = ( + create_brain_mask(mean_b0(volume, gtab)) + if mask is None + else as_mask(volume, mask) + ) + if mu is None: + mu = DTIModel().fit(volume, gtab, mask, n_jobs=n_jobs).pe + else: + norms = np.linalg.norm(mu[mask], axis=-1) + if not (np.isfinite(norms) & (norms > 0)).all(): + raise ValueError("mu must be finite and nonzero in every mask voxel.") + voxels = np.nonzero(mask) + n_voxels = len(voxels[0]) + values = np.zeros((n_voxels, 6)) + fitted = np.zeros(n_voxels, dtype=bool) + bvals = np.ascontiguousarray(gtab.bvals, dtype=np.float64) + bvecs = np.ascontiguousarray(gtab.bvecs, dtype=np.float64) + is_b0 = np.ascontiguousarray(gtab.b0_mask) + + with tqdm(total=n_voxels, desc="Fitting voxels", disable=not progress) as bar: + for start in range(0, n_voxels, CHUNK_VOXELS): + # Each chunk is read and made float64 on its own, so the whole + # masked image is never copied. + chunk = slice(start, min(start + CHUNK_VOXELS, n_voxels)) + index = tuple(axis[chunk] for axis in voxels) + values[chunk], fitted[chunk] = _core.fit_noddi( + np.ascontiguousarray(signal_rows(volume, voxels, chunk)), + np.ascontiguousarray(mu[index], dtype=np.float64), + bvals, + bvecs, + is_b0, + V_IC, + KAPPA, + D_PAR, + D_ISO, + sqrt(LAMBDA), + GAMMA, + threads, + ) + bar.update(chunk.stop - chunk.start) + if not fitted.all(): + LarmorLogger.warning( + "%d mask voxels have a signal that is not finite, no positive b0 or " + "too few positive measurements; they are not fitted and come out as 0.", + int((~fitted).sum()), + ) + + v_ic = np.zeros(shape, dtype=np.float32) + kappa = np.zeros(shape, dtype=np.float32) + v_iso = np.zeros(shape, dtype=np.float32) + fitted_mu = np.zeros((*shape, 3), dtype=np.float32) + v_ic[mask], kappa[mask], v_iso[mask] = values[:, :3].T + fitted_mu[mask] = values[:, 3:] + + # kappa = 0 is fully dispersed, OD = 1; without neurites OD stays 0. + od = np.zeros(shape, dtype=np.float32) + neurites = mask & (v_ic > 0) + with np.errstate(divide="ignore", over="ignore"): + od[neurites] = (2 / np.pi) * np.arctan(1 / kappa[neurites]) + return NoddiFit(v_ic=v_ic, od=od, kappa=kappa, v_iso=v_iso, mu=fitted_mu) diff --git a/src/larmor/main.py b/src/larmor/main.py index 6236da7..590c862 100644 --- a/src/larmor/main.py +++ b/src/larmor/main.py @@ -16,6 +16,7 @@ fod, info, mask, + noddi, peaks, sh_convert, simulate, @@ -34,6 +35,7 @@ ) app.command(name="dti")(larmor_command(dti.dti)) app.command(name="dki")(larmor_command(dki.dki)) +app.command(name="noddi")(larmor_command(noddi.noddi)) app.command(name="denoise")(larmor_command(denoise.denoise)) app.command(name="mask")(larmor_command(mask.mask)) app.command(name="degibbs")(larmor_command(degibbs.degibbs)) diff --git a/src/larmor/sd/csd.py b/src/larmor/sd/csd.py index c74fcff..70027cf 100644 --- a/src/larmor/sd/csd.py +++ b/src/larmor/sd/csd.py @@ -4,10 +4,16 @@ from larmor._log import LarmorLogger from larmor._threads import n_threads -from larmor.dwi.inputs import as_mask, check_dwi +from larmor.dwi.inputs import as_mask, check_dwi, finite_rows, signal_chunks from larmor.dwi.mask import unmask from larmor.sd.csd_rf import estimate_csd_rf -from larmor.sd.csd_solver import CSDSolver, csd_shell, fit_csd +from larmor.sd.csd_solver import ( + CSDSolver, + csd_problem, + csd_shell, + solve_csd, + warn_unconverged, +) from larmor.sd.model import SDFit from larmor.sd.progress import RESPONSE_STAGE, report_stage, report_voxels from larmor.sd.rf import ZonalSHRF @@ -123,8 +129,8 @@ def fit( f"CSD fits one shell: using b = {bval:g} s/mm^2 " f"({int(on_shell.sum())} volumes) and leaving out the others." ) - signals = volume.as_float()[mask][:, on_shell].astype(np.float64) - finite = np.isfinite(signals).all(axis=1) + voxels = np.nonzero(mask) + finite = finite_rows(volume, voxels, on_shell) if not finite.all(): LarmorLogger.warning( f"{int((~finite).sum())} voxels have a NaN or infinite signal; their " @@ -136,17 +142,17 @@ def fit( ) n_fitted = int(finite.sum()) report_voxels(0, n_fitted) - fitted = fit_csd( - signals[finite], - gtab.bvecs[on_shell], - zonal, - self.lmax, - self.solver, - threads, - ) + problem = csd_problem(gtab.bvecs[on_shell], zonal, self.lmax, self.solver) + coeffs = np.full((len(finite), problem.forward.shape[1]), np.nan) + rows_of_fitted = np.flatnonzero(finite) + unconverged = 0 + fitted_voxels = tuple(axis[finite] for axis in voxels) + for rows, signals in signal_chunks(volume, fitted_voxels, on_shell): + fods, converged = solve_csd(signals, problem, threads) + coeffs[rows_of_fitted[rows]] = fods + unconverged += int((~converged).sum()) + warn_unconverged(unconverged, self.solver.max_iter) report_voxels(n_fitted, n_fitted) - coeffs = np.full((len(signals), fitted.shape[1]), np.nan) - coeffs[finite] = fitted if self.basis != "tournier": coeffs = convert_sh_basis(coeffs, "tournier", self.basis) return SDFit( diff --git a/src/larmor/sd/csd_rf.py b/src/larmor/sd/csd_rf.py index 8ad6afe..12e4211 100644 --- a/src/larmor/sd/csd_rf.py +++ b/src/larmor/sd/csd_rf.py @@ -14,8 +14,20 @@ from larmor import _core from larmor._log import LarmorLogger from larmor._threads import n_threads -from larmor.dwi.inputs import as_mask, check_dwi -from larmor.sd.csd_solver import CSDSolver, csd_shell, fit_csd +from larmor.dwi.inputs import ( + as_mask, + check_dwi, + finite_rows, + signal_chunks, + signal_rows, +) +from larmor.sd.csd_solver import ( + CSDSolver, + csd_problem, + csd_shell, + solve_csd, + warn_unconverged, +) from larmor.sd.rf import ZonalSHRF from larmor.sd.sh import build_basis_matrix from larmor.sd.sphere import get_theta_phi @@ -383,8 +395,8 @@ def estimate_csd_rf( bval = csd_shell(gtab) on_shell = gtab.shell_index == len(gtab.shells) - 1 bvecs = gtab.bvecs[on_shell] - data = volume.as_float()[..., on_shell] - usable = np.isfinite(data).all(axis=-1) + usable = np.zeros(brain.shape, dtype=bool) + usable[brain] = finite_rows(volume, np.nonzero(brain), on_shell) response = INITIAL_RESPONSE candidates = brain & usable @@ -393,8 +405,9 @@ def estimate_csd_rf( converged = False for iteration in range(max_iter): fod_lmax = min(2 * (len(response) - 1), MAX_FOD_LMAX) - signals = data[candidates].astype(np.float64) - fods = fit_csd(signals, bvecs, response, fod_lmax, CSDSolver(), threads) + fods = _candidate_fods( + volume, candidates, on_shell, bvecs, response, fod_lmax, threads + ) peaks = _two_largest_peaks(fods, fod_lmax, threads) cost = np.full(brain.shape, np.nan) with np.errstate(divide="ignore", invalid="ignore"): @@ -403,15 +416,15 @@ def estimate_csd_rf( ) selected = _top(cost, num_voxels) # The cost is finite only on candidates, and both masks list voxels in the - # same order, so these are the selected voxels' FODs in the order of - # data[selected]. + # same order, so these are the selected voxels' FODs in the order of their + # signals. fibres = np.array( [_first_lobe_direction(f, fod_lmax) for f in fods[selected[candidates]]] ) response_lmax = INITIAL_LMAX if iteration == 0 else lmax - response = _fit_zonal_response( - data[selected].astype(np.float64), bvecs, fibres, response_lmax - ) + selected_voxels = np.nonzero(selected) + signals = signal_rows(volume, selected_voxels, slice(None), on_shell) + response = _fit_zonal_response(signals, bvecs, fibres, response_lmax) LarmorLogger.info( f"Response iteration {iteration + 1}: {int(selected.sum())} voxels, " f"coefficients {np.array2string(response, precision=3)}" @@ -435,3 +448,38 @@ def estimate_csd_rf( num_iter=iteration + 1, converged=converged, ) + + +def _candidate_fods( + volume: Volume, + candidates: np.ndarray, + on_shell: np.ndarray, + bvecs: np.ndarray, + response: np.ndarray, + lmax: int, + threads: int, +) -> np.ndarray: + """Returns the CSD FODs of the candidate voxels on the shell, read in chunks. + + Args: + volume: The DWI. + candidates: The voxels to fit, (X, Y, Z) bool. + on_shell: The volumes of the shell, bool. + bvecs: Their gradient directions, (n_meas, 3). + response: The response's zonal SH coefficients on the shell. + lmax: Highest degree of the FODs. + threads: Number of threads. + + Returns: + The FODs in the tournier basis, one row per candidate in mask order. + """ + solver = CSDSolver() + problem = csd_problem(bvecs, response, lmax, solver) + voxels = np.nonzero(candidates) + fods = np.empty((len(voxels[0]), problem.forward.shape[1])) + unconverged = 0 + for rows, signals in signal_chunks(volume, voxels, on_shell): + fods[rows], converged = solve_csd(signals, problem, threads) + unconverged += int((~converged).sum()) + warn_unconverged(unconverged, solver.max_iter) + return fods diff --git a/src/larmor/sd/csd_solver.py b/src/larmor/sd/csd_solver.py index f4a49b0..1c61442 100644 --- a/src/larmor/sd/csd_solver.py +++ b/src/larmor/sd/csd_solver.py @@ -159,35 +159,46 @@ def _initial_matrix( return initial -def fit_csd( - signals: npt.NDArray[np.float64], - bvecs: np.ndarray, - zonal: np.ndarray, - lmax: int, - solver: CSDSolver, - n_threads: int, -) -> npt.NDArray[np.float64]: - """Fits the FOD of many voxels of one shell by CSD, in C++. +@dataclass(frozen=True) +class CSDProblem: + """Holds the matrices and settings of the CSD fit of one shell, for `solve_csd`. + + Attributes: + initial: Maps a voxel's signal to its starting FOD. + forward: Maps FOD coefficients to the signal on the shell. + constraint: The FOD amplitudes on the constraint directions, weighted. + norm_weight: Weight of the norm penalty. + threshold: Amplitude below which a direction is constrained. + max_iter: Iteration limit of each voxel's fit. + """ + + initial: npt.NDArray[np.float64] + forward: npt.NDArray[np.float64] + constraint: npt.NDArray[np.float64] + norm_weight: float + threshold: float + max_iter: int + + +def csd_problem( + bvecs: np.ndarray, zonal: np.ndarray, lmax: int, solver: CSDSolver +) -> CSDProblem: + """Checks the settings of a CSD fit of one shell and builds its matrices. Args: - signals: Signal of each voxel on the shell, one per row, shape - (n_voxels, n_meas), in the units of the response; negative values are - fitted as 0. bvecs: Unit gradient direction of each measurement, shape (n_meas, 3). zonal: The response's zonal SH coefficients on this shell, l = 0, 2, ... lmax: Highest degree of the FOD. solver: The CSD settings. - n_threads: Number of threads. Returns: - The FODs in the tournier basis, shape (n_voxels, C). + The problem `solve_csd` fits voxels with. Raises: ValueError: If the shell has fewer than 6 measurements, the response's l = 0 coefficient is not positive or it has no anisotropic degree, or norm_lambda is 0 while the data or the response do not determine every - coefficient (more coefficients than measurements, or lmax above the - response's degrees), which leaves the solve singular. + coefficient, which leaves the solve singular. """ if len(bvecs) < 6: raise ValueError( @@ -220,20 +231,84 @@ def fit_csd( constraint = _basis(directions, lmax) * weight norm_weight = solver.norm_lambda * NORM_SCALE * float(forward[:, 0] @ forward[:, 0]) initial = _initial_matrix(bvecs, zonal, lmax, solver.init_lmax) - fods, converged = _core.fit_csd( - np.ascontiguousarray(np.maximum(signals, 0.0), dtype=np.float64), - np.ascontiguousarray(initial), - np.ascontiguousarray(forward), - np.ascontiguousarray(constraint), + return CSDProblem( + initial=np.ascontiguousarray(initial), + forward=np.ascontiguousarray(forward), + constraint=np.ascontiguousarray(constraint), norm_weight=norm_weight, threshold=solver.threshold * weight, max_iter=solver.max_iter, + ) + + +def solve_csd( + signals: npt.NDArray[np.float64], problem: CSDProblem, n_threads: int +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.bool_]]: + """Fits the FODs of some voxels, and says which converged. + + Args: + signals: Signal of each voxel on the shell, one per row, shape + (n_voxels, n_meas); negative values are fitted as 0. + problem: The matrices from `csd_problem`. + n_threads: Number of threads. + + Returns: + The FODs in the tournier basis, shape (n_voxels, C), and whether each + voxel's fit converged, shape (n_voxels,). + """ + fods, converged = _core.fit_csd( + np.ascontiguousarray(np.maximum(signals, 0.0), dtype=np.float64), + problem.initial, + problem.forward, + problem.constraint, + norm_weight=problem.norm_weight, + threshold=problem.threshold, + max_iter=problem.max_iter, n_threads=n_threads, ) - if not converged.all(): + return fods, converged + + +def warn_unconverged(n_unconverged: int, max_iter: int) -> None: + """Warns that some voxels' CSD fits stopped at the iteration limit.""" + if n_unconverged: LarmorLogger.warning( - f"{int((~converged).sum())} voxels did not converge within " - f"{solver.max_iter} iterations; their FODs are those of the last one." + f"{n_unconverged} voxels did not converge within {max_iter} " + "iterations; their FODs are those of the last one." ) - result: npt.NDArray[np.float64] = fods - return result + + +def fit_csd( + signals: npt.NDArray[np.float64], + bvecs: np.ndarray, + zonal: np.ndarray, + lmax: int, + solver: CSDSolver, + n_threads: int, +) -> npt.NDArray[np.float64]: + """Fits the FOD of many voxels of one shell by CSD, in C++. + + Args: + signals: Signal of each voxel on the shell, one per row, shape + (n_voxels, n_meas), in the units of the response; negative values are + fitted as 0. + bvecs: Unit gradient direction of each measurement, shape (n_meas, 3). + zonal: The response's zonal SH coefficients on this shell, l = 0, 2, ... + lmax: Highest degree of the FOD. + solver: The CSD settings. + n_threads: Number of threads. + + Returns: + The FODs in the tournier basis, shape (n_voxels, C). + + Raises: + ValueError: If the shell has fewer than 6 measurements, the response's l = 0 + coefficient is not positive or it has no anisotropic degree, or + norm_lambda is 0 while the data or the response do not determine every + coefficient (more coefficients than measurements, or lmax above the + response's degrees), which leaves the solve singular. + """ + problem = csd_problem(bvecs, zonal, lmax, solver) + fods, converged = solve_csd(signals, problem, n_threads) + warn_unconverged(int((~converged).sum()), solver.max_iter) + return fods diff --git a/src/larmor/sd/model.py b/src/larmor/sd/model.py index 5295061..e91e84c 100644 --- a/src/larmor/sd/model.py +++ b/src/larmor/sd/model.py @@ -11,7 +11,7 @@ from larmor._log import LarmorLogger from larmor._threads import n_threads from larmor.dwi import estimate_s0, shell_mask -from larmor.dwi.inputs import as_mask, check_dwi +from larmor.dwi.inputs import as_mask, check_dwi, signal_chunks from larmor.dwi.mask import unmask from larmor.sd.kernel import Kernel, order_rfs from larmor.sd.method import SDMethod @@ -212,16 +212,15 @@ def fit( else: mask = as_mask(volume, mask) - # Scaled in float32 as the image is read, then float64 for the solvers. s0 = estimate_s0(volume, gtab, mask)[mask] - signals = volume.as_float()[mask].astype(np.float64) / s0[:, np.newaxis] + voxels = np.nonzero(mask) # The solvers would turn a NaN signal into a small, plausible FOD, so such # voxels are left out and given NaN. Checked before the cap, which would # turn +inf into 1. - finite = np.isfinite(signals).all(axis=1) - # Diffusion weighting only lowers the signal, so a value above S0 is noise. - np.minimum(signals, 1.0, out=signals) + finite = np.zeros(len(s0), dtype=bool) + for rows, signals in signal_chunks(volume, voxels): + finite[rows] = np.isfinite(signals / s0[rows, np.newaxis]).all(axis=1) if not finite.all(): LarmorLogger.warning( f"{int((~finite).sum())} voxels have a NaN or infinite signal or " @@ -232,28 +231,24 @@ def fit( f"Processing {int(finite.sum())} voxels across {threads} cores " f"using {self.method} with {solver}..." ) - fitted_fods, fitted_fractions = self._run_chunks( - signals[finite], hr_kernel.k_aniso.shape[1], fit_chunk + fitted_sh, fitted_fractions = self._run_chunks( + volume, + tuple(axis[finite] for axis in voxels), + s0[finite], + hr_kernel.aniso_blocks, + fit_chunk, + lmax, + exact_integral=isinstance(solver, (L2NNLS, Lasso)), ) - amplitudes = np.full((len(signals), fitted_fods.shape[1]), np.nan) - amplitudes[finite] = fitted_fods - exact_integral = isinstance(solver, (L2NNLS, Lasso)) - fods = [ - unmask( - _to_sh( - amplitudes[:, block], - lmax, - self.basis, - exact_integral=exact_integral, - ), - mask, - ) - for block in hr_kernel.aniso_blocks - ] + fods = [] + for sh in fitted_sh: + coefficients = np.full((len(finite), sh.shape[1]), np.nan) + coefficients[finite] = sh + fods.append(unmask(coefficients, mask)) fractions = None if fitted_fractions is not None: - fractions = np.full((len(signals), fitted_fractions.shape[1]), np.nan) + fractions = np.full((len(finite), fitted_fractions.shape[1]), np.nan) fractions[finite] = fitted_fractions fractions = unmask(fractions, mask) return SDFit( @@ -326,38 +321,74 @@ def fit_multi( return hr, fit_multi def _run_chunks( - self, signals: np.ndarray, n_cols: int, fit_chunk: _ChunkFit - ) -> tuple[np.ndarray, np.ndarray | None]: + self, + volume: Volume, + voxels: tuple[np.ndarray, ...], + s0: np.ndarray, + blocks: list[slice], + fit_chunk: _ChunkFit, + lmax: int, + *, + exact_integral: bool, + ) -> tuple[list[np.ndarray], np.ndarray | None]: """Fits the voxels in chunks of FIT_CHUNK, with a progress bar. - Warns when some solves stopped at the solver's iteration limit. + Each chunk is read, divided by S0 and capped at 1, fitted, and its FOD + amplitudes turned into SH coefficients, so only one chunk of signals and + amplitudes is held at a time. Warns when some solves stopped at the + solver's iteration limit. + + Args: + volume: The DWI. + voxels: The voxels to fit, as `np.nonzero` gives them. + s0: Their S0. + blocks: The columns of each anisotropic tissue's amplitudes. + fit_chunk: Fits one chunk of normalised signals. + lmax: Maximum SH degree of the FODs. + exact_integral: Passed to `_to_sh`. Returns: - The FOD amplitudes, and the fractions of a multi-tissue fit or None. + The SH coefficients of each anisotropic tissue, one row per voxel, and + the fractions of a multi-tissue fit or None. """ - fods = np.empty((len(signals), n_cols)) + n_voxels = len(s0) + n_coefficients = (lmax + 1) * (lmax + 2) // 2 + sh = [np.empty((n_voxels, n_coefficients)) for _ in blocks] fractions = ( - np.empty((len(signals), len(self.rfs))) if self.is_multi_tissue else None + np.empty((n_voxels, len(self.rfs))) if self.is_multi_tissue else None ) unconverged = 0 - starts = range(0, len(signals), FIT_CHUNK) - report_voxels(0, len(signals)) - for start in tqdm(starts, desc="Fitting voxels", disable=not self.progress): - stop = min(start + FIT_CHUNK, len(signals)) - fods[start:stop], chunk_fractions, chunk_unconverged = fit_chunk( - signals[start:stop] - ) + chunks = signal_chunks(volume, voxels, size=FIT_CHUNK) + report_voxels(0, n_voxels) + for rows, signals in tqdm( + chunks, + desc="Fitting voxels", + total=-(-n_voxels // FIT_CHUNK), + disable=not self.progress, + ): + signals /= s0[rows, np.newaxis] + # Diffusion weighting only lowers the signal, so a value above S0 is + # noise. + np.minimum(signals, 1.0, out=signals) + amplitudes, chunk_fractions, chunk_unconverged = fit_chunk(signals) + for tissue, block in zip(sh, blocks, strict=True): + tissue[rows] = _to_sh( + amplitudes[:, block], + lmax, + self.basis, + exact_integral=exact_integral, + ) unconverged += chunk_unconverged if fractions is not None and chunk_fractions is not None: - fractions[start:stop] = chunk_fractions - report_voxels(stop, len(signals)) + fractions[rows] = chunk_fractions + report_voxels(rows.stop, n_voxels) # Only the elastic net counts solves that stop at their limit. if unconverged: LarmorLogger.warning( f"{unconverged} elastic net solves did not converge within max_iter " "sweeps; their FODs are those of the last sweep." ) - return fods, fractions + return sh, fractions def low_resolution_dirs(rfs: list[ResponseFunction], n_meas: int) -> list[np.ndarray]: diff --git a/src/larmor/sd/msmt_csd.py b/src/larmor/sd/msmt_csd.py index 758182d..19615f9 100644 --- a/src/larmor/sd/msmt_csd.py +++ b/src/larmor/sd/msmt_csd.py @@ -5,7 +5,7 @@ from larmor._log import LarmorLogger from larmor._threads import n_threads -from larmor.dwi.inputs import as_mask, check_dwi +from larmor.dwi.inputs import as_mask, check_dwi, finite_rows, signal_chunks from larmor.dwi.mask import unmask from larmor.sd.model import FIT_CHUNK, SDFit from larmor.sd.msmt import fit_msmt, prepare_msmt @@ -148,8 +148,8 @@ def fit( starts = np.concatenate([[0], np.cumsum(sizes)[:-1]]).astype(int) design = prepare_msmt(rfs, lmaxes, gtab.bvals, gtab.bvecs, self.directions) - signals = volume.as_float()[mask].astype(np.float64) - finite = np.isfinite(signals).all(axis=1) + voxels = np.nonzero(mask) + finite = finite_rows(volume, voxels) if not finite.all(): LarmorLogger.warning( f"{int((~finite).sum())} voxels have a NaN or infinite signal; their " @@ -159,23 +159,28 @@ def fit( f"Processing {int(finite.sum())} voxels across {threads} cores using " "MSMT-CSD..." ) - fitted = signals[finite] - fitted_coefficients = np.empty((len(fitted), int(sum(sizes)))) - converged = np.empty(len(fitted), dtype=bool) - chunks = range(0, len(fitted), FIT_CHUNK) - report_voxels(0, len(fitted)) - for start in tqdm(chunks, desc="Fitting voxels", disable=not self.progress): - stop = min(start + FIT_CHUNK, len(fitted)) - fitted_coefficients[start:stop], converged[start:stop] = fit_msmt( - fitted[start:stop], design, threads + fitted_voxels = tuple(axis[finite] for axis in voxels) + n_fitted = int(finite.sum()) + fitted_coefficients = np.empty((n_fitted, int(sum(sizes)))) + converged = np.empty(n_fitted, dtype=bool) + chunks = signal_chunks(volume, fitted_voxels, size=FIT_CHUNK) + report_voxels(0, n_fitted) + for rows, signals in tqdm( + chunks, + desc="Fitting voxels", + total=-(-n_fitted // FIT_CHUNK), + disable=not self.progress, + ): + fitted_coefficients[rows], converged[rows] = fit_msmt( + signals, design, threads ) - report_voxels(stop, len(fitted)) + report_voxels(rows.stop, n_fitted) if not converged.all(): LarmorLogger.warning( f"The constrained fit of {int((~converged).sum())} voxels stopped at " "its iteration limit; their FODs may be slightly negative." ) - coefficients = np.full((len(signals), fitted_coefficients.shape[1]), np.nan) + coefficients = np.full((len(finite), fitted_coefficients.shape[1]), np.nan) coefficients[finite] = fitted_coefficients fods, fractions = _fods_and_fractions( diff --git a/src/larmor/sd/rf.py b/src/larmor/sd/rf.py index ed45e5b..34dc76b 100644 --- a/src/larmor/sd/rf.py +++ b/src/larmor/sd/rf.py @@ -11,7 +11,7 @@ from larmor._log import LarmorLogger from larmor.dwi.fitting import fit_log_signal -from larmor.dwi.inputs import as_mask +from larmor.dwi.inputs import as_mask, signal_chunks from larmor.dwi.shells import SHELL_TOLERANCE from larmor.sd.sphere import get_theta_phi from larmor.types import GradientTable, Volume @@ -199,9 +199,13 @@ def estimate_wm_rf( columns.append((gtab.bvals * 1e-3) ** 2 / 6) design = np.column_stack(columns) - s = volume.as_float()[mask].astype(np.float64) - usable = np.isfinite(s).all(axis=1) & (s > 0).all(axis=1) - params, _ = fit_log_signal(design, s[usable], weighted=True) + fitted = [] + for _, s in signal_chunks(volume, np.nonzero(mask)): + usable = np.isfinite(s).all(axis=1) & (s > 0).all(axis=1) + fitted.append(fit_log_signal(design, s[usable], weighted=True)[0]) + params = ( + np.concatenate(fitted, axis=1) if fitted else np.empty((design.shape[1], 0)) + ) e_vals = _tensor_eigenvalues(params[1:7]) fa = _fractional_anisotropy(e_vals) diff --git a/src/larmor/sd/tissue_responses.py b/src/larmor/sd/tissue_responses.py index a370bb4..4e09f19 100644 --- a/src/larmor/sd/tissue_responses.py +++ b/src/larmor/sd/tissue_responses.py @@ -10,7 +10,7 @@ from larmor import _core from larmor._log import LarmorLogger from larmor._threads import n_threads -from larmor.dwi.inputs import as_mask, check_dwi +from larmor.dwi.inputs import as_mask, check_dwi, signal_chunks, signal_rows from larmor.sd.msmt import fit_msmt, prepare_msmt from larmor.sd.peaks import find_peaks from larmor.sd.rf import ZonalSHRF @@ -182,15 +182,13 @@ def estimate_tissue_responses( rows = np.where(gtab.b0_mask, 0, gtab.shell_index + 1) eroded = _erode(brain, erode) - signal = volume.as_float()[eroded].astype(np.float64) + eroded_voxels = np.nonzero(eroded) - sdm, valid = _signal_decay_metric(signal, rows, len(shell_bvals)) - fa, fibre_dirs = _tensor_fa_and_direction(signal[valid], gtab) - # The metric is kept as the float32 values an image would hold, compared in - # float64. - sdm = sdm[valid].astype(np.float64) - signal = signal[valid] + valid, sdm, fa, fibre_dirs = _voxel_metrics( + volume, eroded_voxels, gtab, rows, len(shell_bvals) + ) safe = _flat_to_mask(eroded, valid) + safe_voxels = tuple(axis[valid] for axis in eroded_voxels) LarmorLogger.info( f"Mask: {int(brain.sum())} -> {int(eroded.sum())} -> {int(safe.sum())} voxels" ) @@ -223,21 +221,29 @@ def estimate_tissue_responses( n_shells = len(shell_bvals) csf = ZonalSHRF( bvals=shell_bvals, - coefficients=_isotropic_response(signal[csf_voxels], rows, n_shells), + coefficients=_isotropic_response( + _signals_of(volume, safe_voxels, csf_voxels), rows, n_shells + ), ) gm = ZonalSHRF( bvals=shell_bvals, - coefficients=_isotropic_response(signal[gm_voxels], rows, n_shells), + coefficients=_isotropic_response( + _signals_of(volume, safe_voxels, gm_voxels), rows, n_shells + ), ) n_wm = round(refined_wm.sum() * wm_percent / 100) wm_voxels = _single_fibre_voxels( - signal, refined_wm, csf, gtab, shell_bvals, n_wm, threads + volume, safe_voxels, refined_wm, csf, gtab, shell_bvals, n_wm, threads ) wm = ZonalSHRF( bvals=shell_bvals, coefficients=_anisotropic_response( - signal[wm_voxels], fibre_dirs[wm_voxels], gtab, rows, n_shells + _signals_of(volume, safe_voxels, wm_voxels), + fibre_dirs[wm_voxels], + gtab, + rows, + n_shells, ), ) LarmorLogger.info( @@ -279,6 +285,47 @@ def _wm_outliers(sdm: np.ndarray, crude_wm: np.ndarray) -> np.ndarray: return outliers +def _voxel_metrics( + volume: Volume, + voxels: tuple[np.ndarray, ...], + gtab: GradientTable, + rows: np.ndarray, + n_shells: int, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Returns the numbers the voxel classes come from, reading the voxels in chunks. + + Returns: + Which voxels have a valid signal decay metric, (n_voxels,) bool, and for + those voxels the metric as float64, the FA and the fibre direction. + """ + valid_parts, sdm_parts, fa_parts, direction_parts = [], [], [], [] + for _, chunk in signal_chunks(volume, voxels): + chunk_sdm, chunk_valid = _signal_decay_metric(chunk, rows, n_shells) + chunk_fa, chunk_dirs = _tensor_fa_and_direction(chunk[chunk_valid], gtab) + valid_parts.append(chunk_valid) + sdm_parts.append(chunk_sdm[chunk_valid]) + fa_parts.append(chunk_fa) + direction_parts.append(chunk_dirs) + if not valid_parts: + return np.zeros(0, bool), np.zeros(0), np.zeros(0, np.float32), np.zeros((0, 3)) + # The metric is kept as the float32 values an image would hold, compared in + # float64. + return ( + np.concatenate(valid_parts), + np.concatenate(sdm_parts).astype(np.float64), + np.concatenate(fa_parts), + np.concatenate(direction_parts), + ) + + +def _signals_of( + volume: Volume, voxels: tuple[np.ndarray, ...], selected: np.ndarray +) -> np.ndarray: + """Returns the float64 signals of the voxels that `selected` keeps, one per row.""" + kept = tuple(axis[selected] for axis in voxels) + return signal_rows(volume, kept, slice(None)) + + def _flat_to_mask(mask: np.ndarray, selected: np.ndarray) -> np.ndarray: """Returns the voxels of mask that `selected`, one entry per True voxel, keeps.""" out = np.zeros(mask.shape, dtype=bool) @@ -597,7 +644,8 @@ def _anisotropic_response( def _single_fibre_voxels( - signal: np.ndarray, + volume: Volume, + voxels: tuple[np.ndarray, ...], refined_wm: np.ndarray, csf: ZonalSHRF, gtab: GradientTable, @@ -612,11 +660,16 @@ def _single_fibre_voxels( largest peak amplitude over the summed l = 0 coefficients; a fit at lmax 6 over those keeps the count best. + The signals of `voxels` are read chunk by chunk from `volume`; `refined_wm` + has one entry per voxel. + Returns: - Which voxels (entries of signal) are kept, (n_voxels,) bool. + Which voxels are kept, (n_voxels,) bool. """ + wm_voxels = tuple(axis[refined_wm] for axis in voxels) + means = [chunk.mean(axis=1) for _, chunk in signal_chunks(volume, wm_voxels)] level = np.sqrt(4 * np.pi) * _median( - signal[refined_wm].mean(axis=1).astype(np.float32) + (np.concatenate(means) if means else np.zeros(0)).astype(np.float32) ) generic = np.tile([level, -level, level, -level], (len(shell_bvals), 1)) generic[0, 1:] = 0.0 @@ -628,11 +681,17 @@ def _single_fibre_voxels( ) # Selection only ranks the voxels, so one that stopped at the iteration limit # still counts with its last fit. - coefficients, _ = fit_msmt(signal[selected], design, threads) + selected_voxels = tuple(axis[selected] for axis in voxels) + coefficients = np.concatenate( + [ + fit_msmt(chunk, design, threads)[0] + for _, chunk in signal_chunks(volume, selected_voxels) + ] + ) n_fod = (lmax + 1) * (lmax + 2) // 2 peaks = find_peaks(coefficients[:, :n_fod], num_peaks=1, n_jobs=threads) total = coefficients[:, 0] + coefficients[:, n_fod] - metric = np.full(len(signal), 0.0, dtype=np.float32) + metric = np.full(len(refined_wm), 0.0, dtype=np.float32) metric[selected] = (peaks.amplitudes[:, 0] / total).astype(np.float32) selected = selected & _top(metric, n_keep) return selected diff --git a/tests/cli/test_cli.py b/tests/cli/test_cli.py index d217376..fc42655 100644 --- a/tests/cli/test_cli.py +++ b/tests/cli/test_cli.py @@ -20,7 +20,7 @@ COMMANDS = ["dti", "mask", "fod"] # Commands that need a multi-shell DWI, which the single-shell fixtures here are # not; only their checks before any work are tested in this file. -MULTI_SHELL_COMMANDS = ["dki"] +MULTI_SHELL_COMMANDS = ["dki", "noddi"] # Gradient options for runs whose inputs do not exist. NO_GRADIENTS = ("--bvals", "a", "--bvecs", "b") @@ -72,7 +72,7 @@ def test_every_command_has_the_shared_flags(command: str) -> None: assert flag in result.output -@pytest.mark.parametrize("command", ["fod", "dti", "dki"]) +@pytest.mark.parametrize("command", ["fod", "dti", "dki", "noddi"]) def test_fitting_commands_have_the_threads_flag(command: str) -> None: assert "--nt" in _run(command, "--help").output diff --git a/tests/cli/test_noddi_command.py b/tests/cli/test_noddi_command.py new file mode 100644 index 0000000..0b8844f --- /dev/null +++ b/tests/cli/test_noddi_command.py @@ -0,0 +1,31 @@ +from pathlib import Path + +import numpy as np +from typer.testing import CliRunner, Result + +from larmor.dwi import NoddiModel +from larmor.io import load_volume, read_bvals_bvecs +from larmor.main import app + +NODDI_DIR = Path(__file__).parents[1] / "data" / "noddi" +VOLUME, MASK = NODDI_DIR / "dwi.nii.gz", NODDI_DIR / "mask.nii.gz" +GRADIENTS = ("--bvals", str(NODDI_DIR / "bvals"), "--bvecs", str(NODDI_DIR / "bvecs")) +METRICS = ("v_ic", "od", "kappa", "v_iso", "mu") + + +def _run(*args: str | Path) -> Result: + return CliRunner().invoke(app, ["noddi", *map(str, args)], env={"TERM": "dumb"}) + + +def test_noddi_writes_the_maps_of_the_python_fit(tmp_path: Path) -> None: + result = _run(VOLUME, tmp_path / "out", *GRADIENTS, "--mask", MASK) + + assert result.exit_code == 0, result.output + assert sorted(p.name for p in tmp_path.iterdir()) == sorted( + f"out_{name}.nii.gz" for name in METRICS + ) + gtab = read_bvals_bvecs(NODDI_DIR / "bvals", NODDI_DIR / "bvecs", ref_path=VOLUME) + expected = NoddiModel().fit(load_volume(VOLUME), gtab, load_volume(MASK)) + for name in METRICS: + written = load_volume(tmp_path / f"out_{name}.nii.gz").as_float() + np.testing.assert_array_equal(written, getattr(expected, name), err_msg=name) diff --git a/tests/data/README.md b/tests/data/README.md index 61df45a..79deb5d 100644 --- a/tests/data/README.md +++ b/tests/data/README.md @@ -10,7 +10,7 @@ Every file here is either synthetic or made from one of two openly licensed data - **Synthetic**: made by the scripts named in the table from numbers alone (numpy, or MATLAB for `matlab/`), and the tools' outputs on them: `sh_tournier_ref.npz`, `grl_kurtosis_kernel.mat`, `mif/`, `tck/`, `vtk/`, `tracking/`, `registration/`, `matlab/`, `gibbs/phantom.nii.gz`, `editing/edge_cases*.tck`, `editing/weights.txt`, - `tsf/edges.*`, `connectome/template.nii.gz` and `connectome/weights.txt`. These are part of Larmor, under its + `tsf/edges.*`, `connectome/template.nii.gz`, `connectome/weights.txt` and `noddi/`. These are part of Larmor, under its licence. - **CFIN multi-shell** (Hansen B, Jespersen SN. Data for evaluation of fast kurtosis strategies, b-value optimization and exploration of diffusion MRI contrast. Scientific Data 3 (2016) 160072, @@ -71,6 +71,8 @@ Every file here is either synthetic or made from one of two openly licensed data | `mfod/mfod_ref.mat` | the same run (`PerformDeconv` with 0 workers) | On the 512 mask voxels: their indices in MRIToolkit's voxel order and in NIfTI order, their signals, the FODs (600 columns, single), the fractions (DKI, NODDI, CSF) and the RSS; the mask | | `mfod/mfod_merge_ref.mat` | `../matlab_scripts/make_mfod_merge_ref.m`, with MATLAB R2025a (25.1.0.2943329), MRIToolkit 8a2c9b8 and the NODDI toolbox 1.05, the same setup as `mfod_kernel.mat` | The DKI and NODDI FODs (45 SH coefficients) and the three fractions as `MRTTrack.SaveOutputToNii` writes them, and the FOD `EDTI_Library.WholeBrainTracking_mDRL_fast_exe` merges from them with its `linear` (mFOD-WS) and `majority` interpolation, on the 8x8x8 grid in MRIToolkit's voxel order | | `lasso/lasso_problems.mat` | `../matlab_scripts/make_lasso_ref.m`, with MATLAB R2025a (25.1.0.2943329), from `mfod/mfod_kernel.mat`, `mfod/mfod_ref.mat` and `rng(0)` | Six elastic-net problems (`X` single, `y`): kernel columns of the mFOD setup against b0-normalised voxel signals (60, 120, 135 and 600 columns, 133 rows) and two random ones (100 and 250 columns, 198 rows); for each, `lambda_max` and MATLAB's `lasso` with `Alpha` 0.5 and `Lambda` 0.1, 0.01, 0.001: with the default `RelTol` 1e-4, with `RelTol` 1e-12 (`MaxIter` 1e7), and in single precision, coefficients and intercepts, with the largest violation of the optimality conditions. The tests use the problems only, since Larmor's elastic net has no intercept | +| `noddi/dwi.nii.gz`, `noddi/dwi_noisy.nii.gz`, `noddi/mask.nii.gz`, `noddi/bvals`, `noddi/bvecs`, `noddi/truth.npz` | `make_noddi_phantom.py` (Larmor) | A 9x6x3 grid of NODDI voxels (`noddi_tissue`, no relaxation, S0 = 1000) with v_ic, OD and v_iso along x, y and z and a random mu each; 6 b0 volumes, 30 at b = 1000 and 60 at b = 2000; noise-free and with Rician noise at SNR 30; the true parameters, mu in world coordinates | +| `noddi/noddi_ref.mat` | `../matlab_scripts/make_noddi_ref.m`, with MATLAB R2025a and the NODDI toolbox 1.05 (`WatsonSHStickTortIsoV_B0`, `batch_fitting_single`) | The toolbox's fit of both `noddi/` signals: `ficvf`, `kappa` (divided by 10, as the toolbox keeps it), `odi`, `fiso`, `b0`, `error_code` and `fibredirs` (bvecs frame) | | `tracking/fod.nii.gz`, `tracking/mask.nii.gz`, `tracking/seeds.nii.gz` | `make_tckgen_det_ref.py` (numpy) | A 16x16x5 synthetic FOD (2 mm voxels, x flipped, lmax 8, tournier, float32): a bundle along +x through a crossing, around a quarter circle to +y, fading in its outer half, with a 40 degree kink near the mask's edge; its tracking mask and a 4-voxel seed mask | | `tracking/det_ref.tck`, `tracking/det_ref_seeds.txt` | the same script: `tckgen fod.nii.gz det_ref.tck -algorithm SD_Stream -seed_grid_per_voxel seeds.nii.gz 2 -seed_direction 1,0,0 -mask mask.nii.gz -step 0.2 -angle 30 -minlength 5 -maxlength 80 -cutoff 0.3 -nthreads 0 -output_seeds det_ref_seeds.txt` (MRtrix 3.0.8; without `-noprecomputed`, which crashes SD_Stream in 3.0.8) | The 32 streamlines from 8 grid seeds per voxel, and the index of each one's seed point | | `tracking/det_ref_maxlen.tck`, `tracking/det_ref_maxlen_seeds.txt` | the same command with `-maxlength 10` | The same seeds stopped by the maximum length | diff --git a/tests/data/make_noddi_phantom.py b/tests/data/make_noddi_phantom.py new file mode 100644 index 0000000..ced4eba --- /dev/null +++ b/tests/data/make_noddi_phantom.py @@ -0,0 +1,113 @@ +"""Generates the NODDI phantom in tests/data/noddi/, with known parameters per voxel. + +Each voxel holds one NODDI tissue (`noddi_tissue`) with its own v_ic, OD and +v_iso, about a random mean orientation mu. The grid axes are the parameters: +x runs over V_IC, y over OD and z over V_ISO. The signal is the sum of the +compartment attenuations weighted by their volume fractions, without +relaxation, so the signal fractions a fit recovers equal the volume fractions. + +The protocol has 6 b = 0 volumes and two shells, b = 1000 s/mm^2 with 30 +directions and b = 2000 s/mm^2 with 60. The files are: + +- dwi.nii.gz: the noise-free signal, with S0 = 1000 in every voxel. +- dwi_noisy.nii.gz: the same with Rician noise of sigma = S0 / SNR. +- mask.nii.gz: every voxel. +- bvals, bvecs: the gradients, FSL format. +- truth.npz: v_ic, od, kappa and v_iso (X, Y, Z), and mu (X, Y, Z, 3) in world + coordinates. + +The NODDI Matlab Toolbox is run on these files by +tests/matlab_scripts/make_noddi_ref.m. + +Run from the repository root: + + uv run python tests/data/make_noddi_phantom.py +""" + +import pathlib + +import numpy as np + +from larmor import _core +from larmor.io import save_volume, write_bvals_bvecs +from larmor.simulation import add_noise, create_gradient_table, noddi_tissue +from larmor.types import GradientTable, Volume + +OUT = pathlib.Path("tests/data/noddi") +SEED = 0 +S0 = 1000.0 +SNR = 30.0 +VOXEL_MM = 2.0 + +# Values between the dictionary's grid points too, so a fit is not tested only +# where it is exact. +V_IC = np.linspace(0.1, 0.9, 9) +OD = np.array([0.05, 0.1, 0.2, 0.3, 0.5, 0.7]) +V_ISO = np.array([0.0, 0.1, 0.3]) + +# The NODDI compartments do not depend on the pulse timings. +SMALL_DELTA_MS = 20.0 +BIG_DELTA_MS = 40.0 + + +def kappa_of(od: np.ndarray) -> np.ndarray: + """Returns the Watson kappa of an orientation dispersion index.""" + return np.asarray(1 / np.tan(od * np.pi / 2)) + + +def random_directions(count: int, rng: np.random.Generator) -> np.ndarray: + """Returns unit vectors drawn evenly on the sphere, (count, 3).""" + v = rng.standard_normal((count, 3)) + return v / np.linalg.norm(v, axis=1, keepdims=True) + + +def voxel_signal( + v_ic: float, kappa: float, v_iso: float, mu: np.ndarray, gradients: GradientTable +) -> np.ndarray: + """Returns the b0-normalised NODDI signal of one voxel, one value per volume.""" + tissue = noddi_tissue(intra_fraction=v_ic, free_water_fraction=v_iso, kappa=kappa) + compartments = tissue.compartments + attenuation = _core.compartment_attenuation( + np.array([c.model.code for c in compartments], dtype=np.int64), + np.array([c.engine_parameters for c in compartments], dtype=np.float64), + np.ascontiguousarray(gradients.bvals, dtype=np.float64), + np.ascontiguousarray(gradients.bvecs @ mu, dtype=np.float64), + SMALL_DELTA_MS, + BIG_DELTA_MS, + ) + fractions = np.array([c.fraction for c in compartments]) + return np.asarray(fractions @ attenuation) + + +rng = np.random.default_rng(SEED) +gradients = create_gradient_table([(1000.0, 30), (2000.0, 60)], num_b0=6, seed=SEED) +shape = (len(V_IC), len(OD), len(V_ISO)) +num_volumes = len(gradients.bvals) + +v_ic, od, v_iso = np.meshgrid(V_IC, OD, V_ISO, indexing="ij") +kappa = kappa_of(od) +mu = random_directions(v_ic.size, rng).reshape(*shape, 3) + +signal = np.empty((*shape, num_volumes)) +for index in np.ndindex(shape): + signal[index] = S0 * voxel_signal( + float(v_ic[index]), + float(kappa[index]), + float(v_iso[index]), + mu[index], + gradients, + ) + +affine = np.diag([VOXEL_MM, VOXEL_MM, VOXEL_MM, 1.0]) +noise_free = Volume(signal.astype(np.float32), affine) +(noisy,) = add_noise([noise_free], S0 / SNR, seed=SEED) + +OUT.mkdir(parents=True, exist_ok=True) +save_volume(noise_free, OUT / "dwi.nii.gz") +save_volume(noisy, OUT / "dwi_noisy.nii.gz") +save_volume(Volume(np.ones(shape, dtype=np.uint8), affine), OUT / "mask.nii.gz") +write_bvals_bvecs(gradients, OUT / "bvals", OUT / "bvecs", ref_path=OUT / "dwi.nii.gz") +np.savez_compressed( + OUT / "truth.npz", v_ic=v_ic, od=od, kappa=kappa, v_iso=v_iso, mu=mu +) +print(f"Wrote {OUT}: {shape} voxels, {num_volumes} volumes.") # noqa: T201 (a script run by hand) diff --git a/tests/data/noddi/bvals b/tests/data/noddi/bvals new file mode 100644 index 0000000..d25b4fc --- /dev/null +++ b/tests/data/noddi/bvals @@ -0,0 +1 @@ +0 0 0 0 0 0 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 diff --git a/tests/data/noddi/bvecs b/tests/data/noddi/bvecs new file mode 100644 index 0000000..0ca4f2b --- /dev/null +++ b/tests/data/noddi/bvecs @@ -0,0 +1,3 @@ +0 0 0 0 0 0 -0.10526488494871286 0.1741124754938399 0.61225369954022824 -0.85827374744700757 -0.83607396356365327 -0.37419338632304522 0.26117375921073327 -0.92223969506586545 0.31891766945634775 0.86729205601633141 -0.99057535489632964 -0.53195867001673447 0.33154906776037929 -0.29533860680182777 -0.7681553552671142 0.17330472480764425 -0.5526408251817525 0.70476888621769584 0.77197656125920111 0.39281965394022478 -0.00048336509251634076 -0.59701004176848693 -0.074963952419659322 0.93356316167575326 -0.44667550199905132 0.52291047489037557 -0.1794401434057355 -0.12623422432698922 0.61122443266024107 -0.68799471625480879 -0.11276820733882392 0.6614594893900404 0.34613394216518267 0.85750842039096864 0.27439126129378139 0.74329303636171162 0.99143044192126983 -0.55458179631042015 -0.24357278848578248 -0.38632482632123405 -0.77270081211963193 -0.84260873939241476 -0.9631187878600701 0.68049289675697977 0.37263071715992163 -0.50547226704920523 0.83351537891734095 -0.71380120004089942 0.0605987943818087 0.09504351729327519 -0.41052847231391948 -0.27606942706087945 0.26545144933840092 0.9707050822718325 0.46156209196748932 -0.09304394389455789 -0.9198086580574284 0.055467612365367437 0.55928313447448086 0.5993469328672103 0.22084728368714501 0.53522925647887865 0.38520205695432891 -0.18472681168986937 0.895231608948163 0.66860030313207564 -0.36348780907413303 0.47900365744065054 -0.88128015665919435 0.7132706831741279 0.15511312822452375 -0.67981607432320679 -0.041203099240123608 0.8011311759287344 -0.43373908354159763 0.88867024070799394 -0.25838649957093734 0.43702013276170781 -0.58147700346622033 -0.059809326910878584 -0.82211369484928043 0.25779732886499535 -0.15707268810313146 0.951448377459853 -0.66767525498254943 -0.70492897169804913 -0.45326037119248297 -0.54509628637608787 -0.045324383840348363 0.15134757523553347 +0 0 0 0 0 0 -0.026676429932943418 -0.81779959605984187 -0.72912998332573276 0.5125707089165521 -0.044213377316755391 0.68188850824506697 -0.95964612667736671 -0.3548950574602473 0.29088505467180198 -0.29878862131070855 0.09535677650332551 0.19235927723079338 -0.17411762541424533 -0.75536649356277163 0.42140822422717172 0.72036030528785999 0.8028665768996941 0.021358967946468774 0.59847685012036822 0.87590503569967848 -0.5066658362353782 -0.78261931792466399 0.94041394903954012 0.1647368995558024 -0.29428000484007577 -0.49795250299106442 -0.97139335990485776 0.4280147558827716 0.49775261516798514 -0.49465105436246198 -0.85205590443838819 0.31747641344670258 0.32712223843465743 0.33025011049968056 -0.55083928830939977 -0.56472077724334535 -0.1291550137569929 -0.38089016473925197 -0.38252096941740743 0.50916014651841113 0.58214458760592769 0.11474459064697048 -0.069962867721885091 -0.72892223948706392 -0.92735917933641054 0.6894425958448831 0.53871880106400816 0.43283069667185925 0.23208358438777668 -0.33067603351441499 -0.081071854977031704 0.23555521527142959 -0.031252611319390416 0.23620446512709223 -0.28337090611081162 -0.067611308830606801 0.28983183994157413 0.89978518686743258 0.038441680710718608 0.77883054878463098 -0.79249232505583278 -0.58832416109175745 0.80425137912522138 0.74921931928147112 -0.27727445967735742 0.59234627303104515 -0.64294873260098884 -0.80808492650312136 -0.41923782900292378 -0.29905430594899679 -0.94712886254879836 -0.075386557160480663 0.50886680797931882 0.013781917485963052 -0.80787282880138389 -0.45580981629475459 0.90608933231525846 0.5663281972624914 0.2353795288649384 0.99520938511180634 -0.26480137861887526 0.9509158821987661 -0.96889032943144937 0.035523279704920564 -0.57010765749129921 -0.6906520607954123 -0.88737100623936283 0.80967288508015844 -0.6358863122350995 0.6964966574178435 +0 0 0 0 0 0 0.9940863504156805 0.54853319549665047 0.30580201898404885 0.02524762967198671 0.54683224549875353 0.62849611928428761 0.10424768128290184 0.15337354085730562 0.90204068925803704 0.39815807080380877 0.098323707450978653 0.82463196752082923 0.92723150733100534 0.58497569841407837 0.48202952059822585 0.67159995750813428 0.2235919901065318 0.70911537249430034 0.21419068312235603 0.28014833198544481 0.86214238774699647 0.17631226060238464 0.33166550964711466 0.31831050420417017 0.844914359367398 0.69181496082034533 0.1555511982187241 0.89491244787139657 0.6153430157339872 0.53102128471801358 0.51115953197144681 0.67947043408884644 0.87930787281969069 0.39447962364808459 0.78821666702779247 0.35861665583999192 0.019612782903671615 0.7398389781607978 0.89126315118781119 0.76909627080438936 0.25306349810034967 0.5261598532924997 0.25982185745444786 0.074845081631713578 0.033930828598517086 0.51881277382000379 0.12261389190642075 0.55058644629043041 0.97080636379104623 0.93894626613012133 0.90823660339703483 0.93182584853646777 0.96361756020210709 0.044035143983629899 0.84063152738196345 0.99336374778968006 0.2644797480312519 0.43279320867179155 0.82808490668320667 0.18494926423886709 0.5684913297629175 0.6061388656774509 0.45254733950530923 0.63603963450268164 0.34882551562452113 0.44955481031563571 0.67416136043329333 0.34286184931678348 0.21813969884379425 0.63388599496887188 0.28086800312499505 0.72949775331475397 0.85985863742242408 0.59833009092427702 0.3990136587799285 0.050026139660561499 0.33502020641144514 0.6987744818940913 0.7787688821667601 0.07733772817873498 0.50399335573279858 0.17117103204752826 0.19125820292545367 0.30575166659892322 0.47873480420141917 0.16147716798291015 0.084425903561605065 0.21748530466440841 0.77045071103967355 0.70141736339792427 diff --git a/tests/data/noddi/dwi.nii.gz b/tests/data/noddi/dwi.nii.gz new file mode 100644 index 0000000..4a9dc0a Binary files /dev/null and b/tests/data/noddi/dwi.nii.gz differ diff --git a/tests/data/noddi/dwi_noisy.nii.gz b/tests/data/noddi/dwi_noisy.nii.gz new file mode 100644 index 0000000..eeeca5e Binary files /dev/null and b/tests/data/noddi/dwi_noisy.nii.gz differ diff --git a/tests/data/noddi/mask.nii.gz b/tests/data/noddi/mask.nii.gz new file mode 100644 index 0000000..7bd0f83 Binary files /dev/null and b/tests/data/noddi/mask.nii.gz differ diff --git a/tests/data/noddi/noddi_ref.mat b/tests/data/noddi/noddi_ref.mat new file mode 100644 index 0000000..7dba179 Binary files /dev/null and b/tests/data/noddi/noddi_ref.mat differ diff --git a/tests/data/noddi/truth.npz b/tests/data/noddi/truth.npz new file mode 100644 index 0000000..308a54a Binary files /dev/null and b/tests/data/noddi/truth.npz differ diff --git a/tests/dwi/test_chunked_reading.py b/tests/dwi/test_chunked_reading.py new file mode 100644 index 0000000..96fcb50 --- /dev/null +++ b/tests/dwi/test_chunked_reading.py @@ -0,0 +1,72 @@ +from collections.abc import Callable + +import numpy as np +import pytest + +import larmor.sd.model +import larmor.sd.msmt_csd +from larmor.dwi import DKIModel, DTIModel, inputs +from larmor.sd import CSDModel, estimate_csd_rf +from larmor.sd.sphere import get_dirs +from larmor.types import GradientTable, Volume + +SHAPE = (7, 6, 5) + + +def _gradients(shells: list[float]) -> GradientTable: + dirs = get_dirs(30) + bvals = np.concatenate([np.zeros(4), *[np.full(30, b) for b in shells]]) + bvecs = [np.zeros((4, 3))] + [np.roll(dirs, k, axis=1) for k in range(len(shells))] + return GradientTable(bvals, np.vstack(bvecs)) + + +def _damaged(gtab: GradientTable, scaled: bool) -> Volume: + """One fibre per voxel in a random direction, with a few broken voxels.""" + rng = np.random.default_rng(3) + fibres = rng.standard_normal((*SHAPE, 3)) + fibres /= np.linalg.norm(fibres, axis=-1, keepdims=True) + cos2 = (fibres @ gtab.bvecs.T) ** 2 + signal = 1000 * np.exp(-gtab.bvals * (0.3e-3 + 1.4e-3 * cos2)) + data = np.abs(signal + 10 * rng.standard_normal(signal.shape)).astype(np.float32) + data[1, 1, 1, gtab.b0_mask] = 0 + data[2, 2, 2] = 0 + data[3, 1, 2, 20:25] = -5 + if scaled: + return Volume( + np.round(data * 4).astype(np.int16), np.eye(4), slope=0.25, intercept=1.0 + ) + data[0, 1, 2, 10] = np.nan + data[4, 3, 1, gtab.b0_mask] = np.inf + return Volume(data, np.eye(4)) + + +def _csd(volume: Volume, gtab: GradientTable, mask: np.ndarray) -> np.ndarray: + rf = estimate_csd_rf(volume, gtab, mask, num_voxels=30).rf + return CSDModel(rf, lmax=6).fit(volume, gtab, mask=mask).fods[0] + + +FITS: dict[str, Callable[[Volume, GradientTable, np.ndarray], np.ndarray]] = { + "dti": lambda volume, gtab, mask: DTIModel().fit(volume, gtab, mask).fa, + "dki": lambda volume, gtab, mask: DKIModel().fit(volume, gtab, mask).mk, + "csd": _csd, +} + + +@pytest.mark.parametrize("scaled", [False, True]) +@pytest.mark.parametrize("name", FITS) +def test_fits_do_not_depend_on_the_chunk_size( + monkeypatch: pytest.MonkeyPatch, name: str, scaled: bool +) -> None: + # Chunks of a few voxels split the image at many places, and between the broken + # voxels; only rounding of the batched solves may differ. + gtab = _gradients([2500.0] if name == "csd" else [1000.0, 2500.0]) + volume = _damaged(gtab, scaled) + mask = np.ones(SHAPE, dtype=bool) + whole = FITS[name](volume, gtab, mask) + + monkeypatch.setattr(inputs, "READ_CHUNK", 7) + monkeypatch.setattr(larmor.sd.model, "FIT_CHUNK", 5) + monkeypatch.setattr(larmor.sd.msmt_csd, "FIT_CHUNK", 5) + chunked = FITS[name](volume, gtab, mask) + + np.testing.assert_allclose(chunked, whole, rtol=1e-6, atol=1e-7) diff --git a/tests/dwi/test_noddi_fit.py b/tests/dwi/test_noddi_fit.py new file mode 100644 index 0000000..31f0d17 --- /dev/null +++ b/tests/dwi/test_noddi_fit.py @@ -0,0 +1,227 @@ +import logging +from pathlib import Path +from typing import Any, cast + +import numpy as np +import pytest + +from larmor.dwi import DTIModel, NoddiFit, NoddiModel, noddi +from larmor.io import load_volume, read_bvals_bvecs, read_matlab +from larmor.types import GradientTable, Volume + +NODDI_DIR = Path(__file__).parents[1] / "data" / "noddi" +MAPS = ("v_ic", "od", "v_iso") +# The toolbox's names of the same maps. +TOOLBOX_MAPS = {"v_ic": "ficvf", "od": "odi", "v_iso": "fiso"} +# Bounds on the mean difference to the toolbox, about twice the measured ones. +TOOLBOX_MEAN_DIFFERENCE = {"v_ic": 0.0015, "od": 0.0005, "v_iso": 0.0015} + + +def _mask() -> np.ndarray: + return load_volume(NODDI_DIR / "mask.nii.gz").data > 0 + + +def _gradients() -> GradientTable: + return read_bvals_bvecs( + NODDI_DIR / "bvals", NODDI_DIR / "bvecs", ref_path=NODDI_DIR / "dwi.nii.gz" + ) + + +@pytest.fixture(scope="module") +def truth() -> dict[str, np.ndarray]: + return dict(np.load(NODDI_DIR / "truth.npz")) + + +@pytest.fixture(scope="module") +def toolbox() -> dict[str, dict[str, np.ndarray]]: + ref = cast(dict[str, Any], read_matlab(NODDI_DIR / "noddi_ref.mat")["ref"]) + return { + name: { + ours: np.asarray(ref[name][theirs]) for ours, theirs in TOOLBOX_MAPS.items() + } + for name in ("noise_free", "noisy") + } + + +@pytest.fixture(scope="module") +def fits() -> dict[str, NoddiFit]: + gtab = _gradients() + return { + name: NoddiModel().fit(load_volume(NODDI_DIR / f"{file}.nii.gz"), gtab, _mask()) + for name, file in (("noise_free", "dwi"), ("noisy", "dwi_noisy")) + } + + +def _along_z(volume: Volume) -> np.ndarray: + mu = np.zeros((*volume.data.shape[:3], 3)) + mu[..., 2] = 1 + return mu + + +def _mean_error( + maps: NoddiFit | dict[str, np.ndarray], truth: dict[str, np.ndarray], name: str +) -> float: + value = getattr(maps, name) if isinstance(maps, NoddiFit) else maps[name] + return float(np.abs(value - truth[name]).mean()) + + +@pytest.mark.parametrize("name", MAPS) +def test_noise_free_fit_recovers_the_phantom( + fits: dict[str, NoddiFit], truth: dict[str, np.ndarray], name: str +) -> None: + error = np.abs(getattr(fits["noise_free"], name) - truth[name]) + + assert error.max() < 1e-3 + + +@pytest.mark.parametrize("name", MAPS) +def test_noisy_fit_errors_are_those_of_the_toolbox( + fits: dict[str, NoddiFit], + truth: dict[str, np.ndarray], + toolbox: dict[str, dict[str, np.ndarray]], + name: str, +) -> None: + ours = _mean_error(fits["noisy"], truth, name) + theirs = _mean_error(toolbox["noisy"], truth, name) + + assert ours < 1.1 * theirs, (ours, theirs) + + +@pytest.mark.parametrize("name", MAPS) +def test_noisy_fit_matches_the_toolbox_voxel_by_voxel( + fits: dict[str, NoddiFit], toolbox: dict[str, dict[str, np.ndarray]], name: str +) -> None: + # Same model and likelihood, another start and optimiser: a few voxels reach + # another optimum. + difference = np.abs(getattr(fits["noisy"], name) - toolbox["noisy"][name]) + assert difference.mean() < TOOLBOX_MEAN_DIFFERENCE[name] + assert (difference > 0.02).sum() <= 2 + + +def test_toolbox_reference_recovers_the_noise_free_phantom( + toolbox: dict[str, dict[str, np.ndarray]], truth: dict[str, np.ndarray] +) -> None: + # Checks the reference itself, and that phantom and toolbox share conventions. + for name in MAPS: + assert _mean_error(toolbox["noise_free"], truth, name) < 0.001, name + + +def test_od_follows_from_kappa(fits: dict[str, NoddiFit]) -> None: + fit = fits["noise_free"] + + with np.errstate(divide="ignore"): + expected = (2 / np.pi) * np.arctan(1 / fit.kappa.astype(np.float64)) + expected[fit.v_ic == 0] = 0 + np.testing.assert_allclose(fit.od, expected, atol=1e-6) + + +def test_fit_is_zero_outside_the_mask() -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + mask = np.zeros(volume.data.shape[:3], dtype=bool) + mask[4, 2, :] = True + + fit = NoddiModel().fit(volume, _gradients(), mask) + + for name in ("v_ic", "od", "kappa", "v_iso", "mu"): + value = getattr(fit, name) + assert np.all(value[~mask] == 0), name + assert np.any(value[mask] != 0), name + + +def test_free_water_voxel_has_no_neurites() -> None: + gtab = _gradients() + signal = 1000 * np.exp(-gtab.bvals * 3.0e-3) + volume = Volume(np.tile(signal, (1, 1, 1, 1)).astype(np.float32), np.eye(4)) + mu = np.zeros((1, 1, 1, 3)) + mu[..., 2] = 1 + + fit = NoddiModel().fit(volume, gtab, np.ones((1, 1, 1), dtype=bool), mu=mu) + + assert fit.v_iso[0, 0, 0] == pytest.approx(1, abs=1e-6) + assert fit.v_ic[0, 0, 0] == 0 + assert fit.od[0, 0, 0] == 0 + + +@pytest.mark.parametrize("bad", ["no_b0", "nan"]) +def test_voxel_that_cannot_be_fitted_gets_zeros_and_a_warning( + bad: str, caplog: pytest.LogCaptureFixture +) -> None: + gtab = _gradients() + signal = 1000 * np.exp(-gtab.bvals * 1.0e-3) + if bad == "no_b0": + signal[gtab.b0_mask] = 0 + else: + signal[10] = np.nan + volume = Volume(np.tile(signal, (1, 2, 1, 1)).astype(np.float32), np.eye(4)) + mask = np.array([[[True], [False]]]) + + with caplog.at_level(logging.WARNING, logger="larmor"): + fit = NoddiModel().fit(volume, gtab, mask, mu=_along_z(volume)) + + for name in ("v_ic", "od", "kappa", "v_iso", "mu"): + assert np.all(getattr(fit, name) == 0), name + assert "1 mask voxels have a signal that is not finite" in caplog.text + + +def test_one_shell_is_refused() -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + gtab = _gradients() + keep = gtab.b0_mask | (gtab.bvals == gtab.bvals.max()) + one_shell = GradientTable(gtab.bvals[keep], gtab.bvecs[keep]) + + with pytest.raises(ValueError, match="at least two non-zero b-value shells"): + NoddiModel().fit(volume.with_data(volume.data[..., keep]), one_shell, _mask()) + + +@pytest.mark.parametrize("bad", [0.0, np.nan]) +def test_mu_that_is_zero_or_not_finite_is_refused(bad: float) -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + mu = _along_z(volume) + mu[tuple(np.argwhere(_mask())[0])] = bad + + with pytest.raises(ValueError, match="mu must be finite and nonzero"): + NoddiModel().fit(volume, _gradients(), _mask(), mu=mu) + + +def test_mu_need_not_be_unit_length() -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + mask = np.zeros(volume.data.shape[:3], dtype=bool) + mask[4, 2, :] = True + mu = DTIModel().fit(volume, _gradients(), mask).pe.astype(np.float64) + + unit = NoddiModel().fit(volume, _gradients(), mask, mu=mu) + scaled = NoddiModel().fit(volume, _gradients(), mask, mu=3 * mu) + + for name in ("v_ic", "kappa", "v_iso", "mu"): + np.testing.assert_allclose( + getattr(scaled, name), getattr(unit, name), atol=1e-6 + ) + + +def test_mu_of_the_wrong_shape_is_refused() -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + + with pytest.raises(ValueError, match="mu must have the shape"): + NoddiModel().fit(volume, _gradients(), mu=np.zeros((2, 2, 2, 3))) + + +def test_fit_does_not_depend_on_the_number_of_threads( + monkeypatch: pytest.MonkeyPatch, +) -> None: + # Chunks small enough that the phantom's voxels are shared among the threads. + monkeypatch.setattr(noddi, "CHUNK_VOXELS", 16) + volume = load_volume(NODDI_DIR / "dwi_noisy.nii.gz") + gtab = _gradients() + + one = NoddiModel().fit(volume, gtab, _mask(), n_jobs=1) + many = NoddiModel().fit(volume, gtab, _mask(), n_jobs=4, progress=True) + + for name in ("v_ic", "od", "kappa", "v_iso", "mu"): + np.testing.assert_array_equal(getattr(one, name), getattr(many, name), name) + + +def test_n_jobs_of_zero_is_refused() -> None: + volume = load_volume(NODDI_DIR / "dwi.nii.gz") + + with pytest.raises(ValueError, match="n_jobs"): + NoddiModel().fit(volume, _gradients(), _mask(), n_jobs=0) diff --git a/tests/matlab_scripts/make_noddi_ref.m b/tests/matlab_scripts/make_noddi_ref.m new file mode 100644 index 0000000..d4d1303 --- /dev/null +++ b/tests/matlab_scripts/make_noddi_ref.m @@ -0,0 +1,73 @@ +% Makes ../data/noddi/noddi_ref.mat: the NODDI Matlab Toolbox fit of the +% phantom that ../data/make_noddi_phantom.py writes, for the noise-free and the +% noisy signal. +% +% The model is WatsonSHStickTortIsoV_B0, fitted with batch_fitting_single. +% For each of noise_free and noisy the file holds the maps (X, Y, Z) ficvf +% (v_ic), kappa, odi, fiso (v_iso), b0 and error_code, and fibredirs +% (X, Y, Z, 3). The toolbox keeps kappa divided by 10: odi is computed from +% 10 * kappa, as SaveParamsAsNIfTI does. fibredirs is in the frame of the +% bvecs file, so along the voxel axes with x flipped, since the phantom's +% affine has a positive determinant. +% +% The images are read with niftiread and the ROI file is written here in the +% layout of CreateROI, since the toolbox's own reader needs SPM's compiled +% files. +% +% Run from inside tests/matlab_scripts, with the NODDI Matlab Toolbox (v1.05) +% and the Image Processing Toolbox on the MATLAB path (MATLAB R2025a). + +phantom_dir = fullfile(pwd, '../data/noddi'); +work = tempname; +mkdir(work); +cleanup = onCleanup(@() rmdir(work, 's')); + +% load reads a file without an extension as a .mat file. +copyfile(fullfile(phantom_dir, 'bvals'), fullfile(work, 'bvals.txt')); +copyfile(fullfile(phantom_dir, 'bvecs'), fullfile(work, 'bvecs.txt')); +protocol = FSL2Protocol(fullfile(work, 'bvals.txt'), fullfile(work, 'bvecs.txt')); +model = MakeModel('WatsonSHStickTortIsoV_B0'); +column = @(name) find(strcmp(model.paramsStr, name)); + +in_mask = niftiread(fullfile(phantom_dir, 'mask.nii.gz')) > 0; +[xsize, ysize, zsize] = size(in_mask); + +names = {'noise_free', 'noisy'}; +inputs = {'dwi.nii.gz', 'dwi_noisy.nii.gz'}; +for n = 1:numel(names) + dwi = double(niftiread(fullfile(phantom_dir, inputs{n}))); + + % CreateROI's layout: voxels in mask, x slowest, z fastest. + [k, j, i] = ndgrid(1:zsize, 1:ysize, 1:xsize); + idx = [i(:) j(:) k(:)]; + idx = idx(in_mask(sub2ind(size(in_mask), idx(:, 1), idx(:, 2), idx(:, 3))), :); + flat = reshape(dwi, [], size(dwi, 4)); + roi = flat(sub2ind(size(in_mask), idx(:, 1), idx(:, 2), idx(:, 3)), :); + mask = zeros(size(in_mask)); + mask(sub2ind(size(in_mask), idx(:, 1), idx(:, 2), idx(:, 3))) = 1:size(idx, 1); + roi_file = fullfile(work, [names{n} '_roi.mat']); + save(roi_file, 'roi', 'mask', 'idx'); + + params_file = fullfile(work, [names{n} '_params.mat']); + batch_fitting_single(roi_file, protocol, model, params_file); + params = load(params_file); + mlps = params.mlps; + + voxels = sub2ind(size(in_mask), idx(:, 1), idx(:, 2), idx(:, 3)); + as_map = @(values) accumarray(voxels, values, [numel(in_mask) 1]); + fit = struct(); + for name = {'ficvf', 'kappa', 'fiso', 'b0'} + fit.(name{1}) = reshape(as_map(mlps(:, column(name{1}))), size(in_mask)); + end + odi = atan2(1, mlps(:, column('kappa')) * 10) * 2 / pi; + fit.odi = reshape(as_map(odi), size(in_mask)); + fit.error_code = reshape(as_map(double(params.error_code(:))), size(in_mask)); + dirs = GetFibreOrientation(model.name, mlps)'; + fit.fibredirs = zeros(xsize, ysize, zsize, 3); + for a = 1:3 + fit.fibredirs(:, :, :, a) = reshape(as_map(dirs(:, a)), size(in_mask)); + end + ref.(names{n}) = fit; +end + +save(fullfile(phantom_dir, 'noddi_ref.mat'), 'ref', '-v7'); diff --git a/tests/sd/test_matlab_match.py b/tests/sd/test_matlab_match.py index b71bde9..366ccab 100644 --- a/tests/sd/test_matlab_match.py +++ b/tests/sd/test_matlab_match.py @@ -152,7 +152,7 @@ def test_match_mfod_fit() -> None: hr, fit_chunk = model._prepare( model.method.default_solver(), bvals, bvecs, n_threads=2 ) - fods, fractions = model._run_chunks(signals, hr.k_aniso.shape[1], fit_chunk) + fods, fractions, _ = fit_chunk(signals) # The reference is single precision and the NNLS active set can flip at # rounding level, so the bounds leave headroom over the observed 5e-5 and 2e-6. diff --git a/tests/test_docs.py b/tests/test_docs.py index e87b9f9..d56b4d0 100644 --- a/tests/test_docs.py +++ b/tests/test_docs.py @@ -79,6 +79,7 @@ def test_examples_that_download_data_carry_the_marker() -> None: "gradient-tables", "dti", "dki", + "noddi", "fods", "peaks", "tractography",