NODDI fitting

cpp/src/common/maths.cpp

@@ -1,9 +1,21 @@
 #include "common/maths.hpp"
 
 #include <cmath>
+#include <numbers>
 
 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<double>(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

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.

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)