feat: add small flip angle approximation as an option

gdevenyi · Mar 26, 2024 · 5a693c3 · 5a693c3
1 parent c68a269
commit 5a693c3
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 7 deletions.
diff --git a/Source/MT/MTSatModel.cpp b/Source/MT/MTSatModel.cpp
@@ -8,15 +8,22 @@ auto MTSatModel::fit(QI_ARRAYN(DataType, NI) const &     in,
     auto const  S_mt = in[2];
     auto        B1   = fixed[0];
     auto const &s    = sequence;
-
+
+    auto al_t1_b1corr = B1 * s.al_t1;
+    auto al_pd_b1corr = B1 * s.al_pd;
+    if (!smallangle) {
+        al_t1_b1corr = 2 * tan(al_t1_b1corr / 2);
+        al_pd_b1corr = 2 * tan(al_pd_b1corr / 2);
+    }
+
     auto const R1 =
-        std::clamp(2 * (S_t1 * tan(B1 * s.al_t1 / 2) / s.TR_t1 - S_pd * tan(B1 * s.al_pd / 2) / s.TR_pd) /
-                       (S_pd / tan(B1 * s.al_pd / 2) - S_t1 / tan(B1 * s.al_t1 / 2)),
+        std::clamp(0.5 * (S_t1 * al_t1_b1corr / s.TR_t1 - S_pd * al_pd_b1corr / s.TR_pd) /
+                       (S_pd / al_pd_b1corr - S_t1 / al_t1_b1corr),
                    0.,
                    10.);
     auto const A =
-        std::max(0.5 * S_pd * S_t1 * (s.TR_pd * tan(B1 * s.al_t1 / 2) / tan(B1 * s.al_pd / 2) - s.TR_t1 * tan(B1 * s.al_pd / 2) / tan(B1 * s.al_t1 / 2)) /
-                     (S_t1 * s.TR_pd * tan(B1 * s.al_t1 / 2) - S_pd * s.TR_t1 * tan(B1 * s.al_pd / 2)),
+        std::max((S_pd * S_t1) * (s.TR_pd * al_t1_b1corr / al_pd_b1corr - s.TR_t1 * al_pd_b1corr / al_t1_b1corr) /
+                     (S_t1 * s.TR_pd * al_t1_b1corr - S_pd * s.TR_t1 * al_pd_b1corr),
                  0.);
     auto const d           = (A * s.al_mt / S_mt - 1.0) * R1 * s.TR_mt - s.al_mt * s.al_mt / 2;
     auto const d_corrected = std::clamp(d * (1.0 - C) / (1.0 - C * B1), 0., 0.1);

diff --git a/Source/MT/MTSatModel.h b/Source/MT/MTSatModel.h
@@ -7,6 +7,7 @@ using namespace std::literals;
 struct MTSatModel : QI::Model<double, double, 3, 1, 3, 0> {
     QI::MTSatSequence const &sequence;
     double const             C;
+    bool const               smallangle;
 
     std::array<const std::string, NV> const varying_names{{"PD"s, "R1"s, "delta"s}};
     std::array<const std::string, 3> const  derived_names{};

diff --git a/Source/MT/qi_mtsat.cpp b/Source/MT/qi_mtsat.cpp
@@ -29,13 +29,14 @@ int mtsat_main(args::Subparser &parser) {
     args::ValueFlag<std::string> b1_path(parser, "B1", "Path to B1 map", {'b', "B1"});
     args::ValueFlag<double>      C(
         parser, "C", "Correction factor for delta (default 0.4)", {'C', "C"}, 0.4);
+    args::Flag smallangle(parser, "smallangle", "Use small flip angle approx for R1 and PD calculation", {'s', "smallangle"});
     parser.Parse();
-
+    
     QI::Log(verbose, "Reading sequence parameters");
     json              input = json_file ? QI::ReadJSON(json_file.Get()) : QI::ReadJSON(std::cin);
     QI::MTSatSequence seq(input["MTSat"]);
 
-    MTSatModel model{{}, seq, C.Get()};
+    MTSatModel model{{}, seq, C.Get(), smallangle};
     if (simulate) {
         QI::SimulateModel<MTSatModel, true>(
             input,