diff --git a/Python/Tests/test_perfusion.py b/Python/Tests/test_perfusion.py index 0079d40..db8d0e9 100644 --- a/Python/Tests/test_perfusion.py +++ b/Python/Tests/test_perfusion.py @@ -62,8 +62,8 @@ def test_oef(self): noise=noise, verbose=vb).run() diff_DBV = Diff(in_file='ASE_DBV.nii.gz', baseline='DBV.nii.gz', noise=noise, verbose=vb).run() - self.assertLessEqual(diff_R2p.outputs.out_diff, 12) - self.assertLessEqual(diff_DBV.outputs.out_diff, 50) + self.assertLessEqual(diff_R2p.outputs.out_diff, 20) + self.assertLessEqual(diff_DBV.outputs.out_diff, 75) def test_oef_fixed_dbv(self): # Use MultiEchoFlex as a proxy for ASE diff --git a/Source/Core/FitScaledAuto.h b/Source/Core/FitScaledAuto.h index 7f1fb6f..84ba6c8 100644 --- a/Source/Core/FitScaledAuto.h +++ b/Source/Core/FitScaledAuto.h @@ -36,7 +36,8 @@ template struct ScaledAutoDiffFi using AutoCost = ceres::AutoDiffCostFunction; auto *cost = new Cost{this->model, fixed, data}; auto *auto_cost = new AutoCost(cost, this->model.sequence.size()); - problem.AddResidualBlock(auto_cost, NULL, varying.data()); + auto *loss = new ceres::HuberLoss(1.0); + problem.AddResidualBlock(auto_cost, loss, varying.data()); for (int i = 0; i < ModelType::NV; i++) { problem.SetParameterLowerBound(varying.data(), i, this->model.bounds_lo[i]); problem.SetParameterUpperBound(varying.data(), i, this->model.bounds_hi[i]); diff --git a/Source/Perfusion/qi_ase_oef.cpp b/Source/Perfusion/qi_ase_oef.cpp index b761d49..12b14de 100644 --- a/Source/Perfusion/qi_ase_oef.cpp +++ b/Source/Perfusion/qi_ase_oef.cpp @@ -11,9 +11,12 @@ #include +#include "NumericalIntegration.h" + // #define QI_DEBUG_BUILD #include "Args.h" + #include "FitFunction.h" #include "FitScaledAuto.h" #include "ImageIO.h" @@ -28,45 +31,52 @@ using namespace std::literals; constexpr double kappa = 0.03; // Conversion factor constexpr double gyro_gamma = 2 * M_PI * 42.577e6; // Gyromagnetic Ratio constexpr double delta_X0 = 0.264e-6; // Susc diff oxy and fully de-oxy blood -constexpr double Hct = 0.34; -constexpr double Hb = Hct / kappa; + +template struct fcFunctor { + const T &dw, tau; + T operator()(const T u) const { + const auto Jterm = ceres::BesselJ0(1.5 * dw * tau * u); + return (2. + u) * sqrt(1. - u) * (1. - Jterm) / (u * u); + } +}; struct ASEModel : QI::Model { using SequenceType = QI::MultiEchoSequence; const SequenceType &sequence; - const double B0; + const double B0, Hct; const std::array varying_names{"S0"s, "dT"s, "R2p"s, "DBV"s}; const std::array derived_names{"Tc"s, "OEF"s, "dHb"s}; - const VaryingArray start{0.98, 0., 2.0, 0.02}; + const VaryingArray start{0.98, 0., 5.0, 0.025}; const VaryingArray bounds_lo{0.1, -0.1, 0.25, 0.001}; - const VaryingArray bounds_hi{2., 0.1, 20.0, 0.05}; + const VaryingArray bounds_hi{2., 0.1, 50.0, 0.5}; int input_size(const int /* Unused */) const { return sequence.size(); } template auto signal(const Eigen::ArrayBase &varying, const FixedArray & /* Unused */) const -> QI_ARRAY(typename Derived::Scalar) { - using T = typename Derived::Scalar; - const T &S0 = varying[0]; - const T &dT = varying[1]; - const T &R2p = varying[2]; - const T &DBV = varying[3]; - - const auto dw = R2p / DBV; - const auto Tc = 1.5 * (1. / dw); + using T = typename Derived::Scalar; + const T & S0 = varying[0]; + const T & dT = varying[1]; + const T & R2p = varying[2]; + const T & DBV = varying[3]; + const auto dw = R2p / DBV; - const auto aTE = (sequence.TE + dT).abs(); - QI_ARRAY(T) S(sequence.size()); + Eigen::Integrator integrator(5); + const auto quad_rule = Eigen::Integrator::GaussKronrod15; + T abs_error{1.e-9}; + T rel_error{1.e-9}; + const auto aTE = (sequence.TE + dT).abs(); + QI_ARRAY(T) fc(sequence.size()); for (int i = 0; i < sequence.size(); i++) { - const auto tau = aTE(i); - if (tau < Tc) { - S(i) = S0 * exp(-(0.3) * (R2p * tau) * (R2p * tau) / DBV); - } else { - S(i) = S0 * exp(-tau * R2p + DBV); - } + const auto tau = aTE(i); + fcFunctor functor{dw, tau}; + fc[i] = (1. / 3.) * integrator.quadratureAdaptive( + functor, T{0.0}, T{1.0}, abs_error, rel_error, quad_rule); } + QI_ARRAY(T) S = S0 * exp(-DBV * fc); return S; } @@ -76,9 +86,11 @@ struct ASEModel : QI::Model { const auto &R2p = varying[2]; const auto &DBV = varying[3]; - const auto dw = R2p / DBV; - const auto Tc = 1.5 * (1. / dw); - const auto OEF = dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct); + const auto dw = R2p / DBV; + const auto Tc = 1.5 * (1. / dw); + const auto OEF = + std::clamp(dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct), 0., 1.); + const auto Hb = Hct / kappa; const auto dHb = OEF * Hb; derived[0] = Tc; @@ -91,14 +103,14 @@ using ASEFit = QI::ScaledAutoDiffFit; struct ASEFixDBVModel : QI::Model { using SequenceType = QI::MultiEchoSequence; const SequenceType &sequence; - const double B0, DBV; + const double B0, Hct, DBV; const std::array varying_names{"S0"s, "dT"s, "R2p"s}; const std::array derived_names{"Tc"s, "OEF"s, "dHb"s}; - const VaryingArray start{0.98, 0., 1.0}; + const VaryingArray start{0.98, 0., 5.0}; const VaryingArray bounds_lo{0.1, -0.1, 0.25}; - const VaryingArray bounds_hi{2., 0.1, 20.0}; + const VaryingArray bounds_hi{2., 0.1, 50.0}; int input_size(const int /* Unused */) const { return sequence.size(); } @@ -110,18 +122,21 @@ struct ASEFixDBVModel : QI::Model { const T & dT = varying[1]; const T & R2p = varying[2]; const auto dw = R2p / DBV; - const auto Tc = 1.5 * (1. / dw); - const auto aTE = (sequence.TE + dT).abs(); - QI_ARRAY(T) S(sequence.size()); + Eigen::Integrator integrator(5); + const auto quad_rule = Eigen::Integrator::GaussKronrod15; + T abs_error{1.e-9}; + T rel_error{1.e-9}; + const auto aTE = (sequence.TE + dT).abs(); + QI_ARRAY(T) fc(sequence.size()); for (int i = 0; i < sequence.size(); i++) { - const auto tau = aTE(i); - if (tau < Tc) { - S(i) = S0 * exp(-(0.3) * (R2p * tau) * (R2p * tau) / DBV); - } else { - S(i) = S0 * exp(-tau * R2p + DBV); - } + const auto tau = aTE(i); + fcFunctor functor{dw, tau}; + fc[i] = (1. / 3.) * integrator.quadratureAdaptive( + functor, T{0.0}, T{1.0}, abs_error, rel_error, quad_rule); } + QI_ARRAY(T) S = S0 * exp(-DBV * fc); + return S; } @@ -131,8 +146,10 @@ struct ASEFixDBVModel : QI::Model { const auto &R2p = varying[2]; const auto dw = R2p / DBV; const auto Tc = 1.5 * (1. / dw); - const auto OEF = dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct); - const auto dHb = OEF * Hb; + const auto OEF = + std::clamp(dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct), 0., 1.); + const auto Hb = Hct / kappa; + const auto dHb = OEF * Hb; derived[0] = Tc; derived[1] = OEF; @@ -149,6 +166,7 @@ int ase_oef_main(args::Subparser &parser) { QI_COMMON_ARGS; args::ValueFlag B0(parser, "B0", "Field-strength (Tesla), default 3", {'B', "B0"}, 3.0); + args::ValueFlag Hct(parser, "HCT", "Hematocrit (default 0.34)", {'h', "Hct"}, 0.34); args::ValueFlag DBV(parser, "DBV", "Fix DBV and only fit R2'", {'d', "DBV"}, 0.0); parser.Parse(); @@ -157,7 +175,7 @@ int ase_oef_main(args::Subparser &parser) { if (simulate) { if (DBV) { - ASEFixDBVModel model{{}, sequence, B0.Get(), DBV.Get()}; + ASEFixDBVModel model{{}, sequence, B0.Get(), Hct.Get(), DBV.Get()}; QI::SimulateModel(input, model, {}, @@ -168,7 +186,7 @@ int ase_oef_main(args::Subparser &parser) { threads.Get(), subregion.Get()); } else { - ASEModel model{{}, sequence, B0.Get()}; + ASEModel model{{}, sequence, B0.Get(), Hct.Get()}; QI::SimulateModel(input, model, {}, @@ -188,11 +206,11 @@ int ase_oef_main(args::Subparser &parser) { fit_filter->WriteOutputs(prefix.Get() + "ASE_"); }; if (DBV) { - ASEFixDBVModel model{{}, sequence, B0.Get(), DBV.Get()}; + ASEFixDBVModel model{{}, sequence, B0.Get(), Hct.Get(), DBV.Get()}; ASEFixDBVFit fit{model}; process(fit); } else { - ASEModel model{{}, sequence, B0.Get()}; + ASEModel model{{}, sequence, Hct.Get(), B0.Get()}; ASEFit fit{model}; process(fit); }