Skip to content

Commit ee2b7f5

Browse files
spinicistspinicist
committed
Swap ASE to full integral equation
1 parent d2ae0d8 commit ee2b7f5

File tree

3 files changed

+64
-45
lines changed

3 files changed

+64
-45
lines changed

‎Python/Tests/test_perfusion.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,8 @@ def test_oef(self):
6262
noise=noise, verbose=vb).run()
6363
diff_DBV = Diff(in_file='ASE_DBV.nii.gz', baseline='DBV.nii.gz',
6464
noise=noise, verbose=vb).run()
65-
self.assertLessEqual(diff_R2p.outputs.out_diff, 12)
66-
self.assertLessEqual(diff_DBV.outputs.out_diff, 50)
65+
self.assertLessEqual(diff_R2p.outputs.out_diff, 20)
66+
self.assertLessEqual(diff_DBV.outputs.out_diff, 75)
6767

6868
def test_oef_fixed_dbv(self):
6969
# Use MultiEchoFlex as a proxy for ASE

‎Source/Core/FitScaledAuto.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,8 @@ template <typename ModelType_, typename FlagType_ = int> struct ScaledAutoDiffFi
3636
using AutoCost = ceres::AutoDiffCostFunction<Cost, ceres::DYNAMIC, ModelType::NV>;
3737
auto *cost = new Cost{this->model, fixed, data};
3838
auto *auto_cost = new AutoCost(cost, this->model.sequence.size());
39-
problem.AddResidualBlock(auto_cost, NULL, varying.data());
39+
auto *loss = new ceres::HuberLoss(1.0);
40+
problem.AddResidualBlock(auto_cost, loss, varying.data());
4041
for (int i = 0; i < ModelType::NV; i++) {
4142
problem.SetParameterLowerBound(varying.data(), i, this->model.bounds_lo[i]);
4243
problem.SetParameterUpperBound(varying.data(), i, this->model.bounds_hi[i]);

‎Source/Perfusion/qi_ase_oef.cpp‎

Lines changed: 60 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -11,9 +11,12 @@
1111

1212
#include <Eigen/Dense>
1313

14+
#include "NumericalIntegration.h"
15+
1416
// #define QI_DEBUG_BUILD
1517

1618
#include "Args.h"
19+
1720
#include "FitFunction.h"
1821
#include "FitScaledAuto.h"
1922
#include "ImageIO.h"
@@ -28,45 +31,52 @@ using namespace std::literals;
2831
constexpr double kappa = 0.03; // Conversion factor
2932
constexpr double gyro_gamma = 2 * M_PI * 42.577e6; // Gyromagnetic Ratio
3033
constexpr double delta_X0 = 0.264e-6; // Susc diff oxy and fully de-oxy blood
31-
constexpr double Hct = 0.34;
32-
constexpr double Hb = Hct / kappa;
34+
35+
template <typename T> struct fcFunctor {
36+
const T &dw, tau;
37+
T operator()(const T u) const {
38+
const auto Jterm = ceres::BesselJ0(1.5 * dw * tau * u);
39+
return (2. + u) * sqrt(1. - u) * (1. - Jterm) / (u * u);
40+
}
41+
};
3342

3443
struct ASEModel : QI::Model<double, double, 4, 0, 1, 3> {
3544
using SequenceType = QI::MultiEchoSequence;
3645
const SequenceType &sequence;
37-
const double B0;
46+
const double B0, Hct;
3847

3948
const std::array<const std::string, NV> varying_names{"S0"s, "dT"s, "R2p"s, "DBV"s};
4049
const std::array<const std::string, ND> derived_names{"Tc"s, "OEF"s, "dHb"s};
4150

42-
const VaryingArray start{0.98, 0., 2.0, 0.02};
51+
const VaryingArray start{0.98, 0., 5.0, 0.025};
4352
const VaryingArray bounds_lo{0.1, -0.1, 0.25, 0.001};
44-
const VaryingArray bounds_hi{2., 0.1, 20.0, 0.05};
53+
const VaryingArray bounds_hi{2., 0.1, 50.0, 0.5};
4554

4655
int input_size(const int /* Unused */) const { return sequence.size(); }
4756

4857
template <typename Derived>
4958
auto signal(const Eigen::ArrayBase<Derived> &varying, const FixedArray & /* Unused */) const
5059
-> QI_ARRAY(typename Derived::Scalar) {
51-
using T = typename Derived::Scalar;
52-
const T &S0 = varying[0];
53-
const T &dT = varying[1];
54-
const T &R2p = varying[2];
55-
const T &DBV = varying[3];
56-
57-
const auto dw = R2p / DBV;
58-
const auto Tc = 1.5 * (1. / dw);
60+
using T = typename Derived::Scalar;
61+
const T & S0 = varying[0];
62+
const T & dT = varying[1];
63+
const T & R2p = varying[2];
64+
const T & DBV = varying[3];
65+
const auto dw = R2p / DBV;
5966

60-
const auto aTE = (sequence.TE + dT).abs();
61-
QI_ARRAY(T) S(sequence.size());
67+
Eigen::Integrator<T> integrator(5);
68+
const auto quad_rule = Eigen::Integrator<T>::GaussKronrod15;
69+
T abs_error{1.e-9};
70+
T rel_error{1.e-9};
71+
const auto aTE = (sequence.TE + dT).abs();
72+
QI_ARRAY(T) fc(sequence.size());
6273
for (int i = 0; i < sequence.size(); i++) {
63-
const auto tau = aTE(i);
64-
if (tau < Tc) {
65-
S(i) = S0 * exp(-(0.3) * (R2p * tau) * (R2p * tau) / DBV);
66-
} else {
67-
S(i) = S0 * exp(-tau * R2p + DBV);
68-
}
74+
const auto tau = aTE(i);
75+
fcFunctor<T> functor{dw, tau};
76+
fc[i] = (1. / 3.) * integrator.quadratureAdaptive(
77+
functor, T{0.0}, T{1.0}, abs_error, rel_error, quad_rule);
6978
}
79+
QI_ARRAY(T) S = S0 * exp(-DBV * fc);
7080
return S;
7181
}
7282

@@ -76,9 +86,11 @@ struct ASEModel : QI::Model<double, double, 4, 0, 1, 3> {
7686
const auto &R2p = varying[2];
7787
const auto &DBV = varying[3];
7888

79-
const auto dw = R2p / DBV;
80-
const auto Tc = 1.5 * (1. / dw);
81-
const auto OEF = dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct);
89+
const auto dw = R2p / DBV;
90+
const auto Tc = 1.5 * (1. / dw);
91+
const auto OEF =
92+
std::clamp(dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct), 0., 1.);
93+
const auto Hb = Hct / kappa;
8294
const auto dHb = OEF * Hb;
8395

8496
derived[0] = Tc;
@@ -91,14 +103,14 @@ using ASEFit = QI::ScaledAutoDiffFit<ASEModel>;
91103
struct ASEFixDBVModel : QI::Model<double, double, 3, 0, 1, 3> {
92104
using SequenceType = QI::MultiEchoSequence;
93105
const SequenceType &sequence;
94-
const double B0, DBV;
106+
const double B0, Hct, DBV;
95107

96108
const std::array<const std::string, NV> varying_names{"S0"s, "dT"s, "R2p"s};
97109
const std::array<const std::string, ND> derived_names{"Tc"s, "OEF"s, "dHb"s};
98110

99-
const VaryingArray start{0.98, 0., 1.0};
111+
const VaryingArray start{0.98, 0., 5.0};
100112
const VaryingArray bounds_lo{0.1, -0.1, 0.25};
101-
const VaryingArray bounds_hi{2., 0.1, 20.0};
113+
const VaryingArray bounds_hi{2., 0.1, 50.0};
102114

103115
int input_size(const int /* Unused */) const { return sequence.size(); }
104116

@@ -110,18 +122,21 @@ struct ASEFixDBVModel : QI::Model<double, double, 3, 0, 1, 3> {
110122
const T & dT = varying[1];
111123
const T & R2p = varying[2];
112124
const auto dw = R2p / DBV;
113-
const auto Tc = 1.5 * (1. / dw);
114125

115-
const auto aTE = (sequence.TE + dT).abs();
116-
QI_ARRAY(T) S(sequence.size());
126+
Eigen::Integrator<T> integrator(5);
127+
const auto quad_rule = Eigen::Integrator<T>::GaussKronrod15;
128+
T abs_error{1.e-9};
129+
T rel_error{1.e-9};
130+
const auto aTE = (sequence.TE + dT).abs();
131+
QI_ARRAY(T) fc(sequence.size());
117132
for (int i = 0; i < sequence.size(); i++) {
118-
const auto tau = aTE(i);
119-
if (tau < Tc) {
120-
S(i) = S0 * exp(-(0.3) * (R2p * tau) * (R2p * tau) / DBV);
121-
} else {
122-
S(i) = S0 * exp(-tau * R2p + DBV);
123-
}
133+
const auto tau = aTE(i);
134+
fcFunctor<T> functor{dw, tau};
135+
fc[i] = (1. / 3.) * integrator.quadratureAdaptive(
136+
functor, T{0.0}, T{1.0}, abs_error, rel_error, quad_rule);
124137
}
138+
QI_ARRAY(T) S = S0 * exp(-DBV * fc);
139+
125140
return S;
126141
}
127142

@@ -131,8 +146,10 @@ struct ASEFixDBVModel : QI::Model<double, double, 3, 0, 1, 3> {
131146
const auto &R2p = varying[2];
132147
const auto dw = R2p / DBV;
133148
const auto Tc = 1.5 * (1. / dw);
134-
const auto OEF = dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct);
135-
const auto dHb = OEF * Hb;
149+
const auto OEF =
150+
std::clamp(dw / ((4. * M_PI / 3.) * gyro_gamma * B0 * delta_X0 * Hct), 0., 1.);
151+
const auto Hb = Hct / kappa;
152+
const auto dHb = OEF * Hb;
136153

137154
derived[0] = Tc;
138155
derived[1] = OEF;
@@ -149,6 +166,7 @@ int ase_oef_main(args::Subparser &parser) {
149166

150167
QI_COMMON_ARGS;
151168
args::ValueFlag<double> B0(parser, "B0", "Field-strength (Tesla), default 3", {'B', "B0"}, 3.0);
169+
args::ValueFlag<double> Hct(parser, "HCT", "Hematocrit (default 0.34)", {'h', "Hct"}, 0.34);
152170
args::ValueFlag<double> DBV(parser, "DBV", "Fix DBV and only fit R2'", {'d', "DBV"}, 0.0);
153171

154172
parser.Parse();
@@ -157,7 +175,7 @@ int ase_oef_main(args::Subparser &parser) {
157175

158176
if (simulate) {
159177
if (DBV) {
160-
ASEFixDBVModel model{{}, sequence, B0.Get(), DBV.Get()};
178+
ASEFixDBVModel model{{}, sequence, B0.Get(), Hct.Get(), DBV.Get()};
161179
QI::SimulateModel<ASEFixDBVModel, false>(input,
162180
model,
163181
{},
@@ -168,7 +186,7 @@ int ase_oef_main(args::Subparser &parser) {
168186
threads.Get(),
169187
subregion.Get());
170188
} else {
171-
ASEModel model{{}, sequence, B0.Get()};
189+
ASEModel model{{}, sequence, B0.Get(), Hct.Get()};
172190
QI::SimulateModel<ASEModel, false>(input,
173191
model,
174192
{},
@@ -188,11 +206,11 @@ int ase_oef_main(args::Subparser &parser) {
188206
fit_filter->WriteOutputs(prefix.Get() + "ASE_");
189207
};
190208
if (DBV) {
191-
ASEFixDBVModel model{{}, sequence, B0.Get(), DBV.Get()};
209+
ASEFixDBVModel model{{}, sequence, B0.Get(), Hct.Get(), DBV.Get()};
192210
ASEFixDBVFit fit{model};
193211
process(fit);
194212
} else {
195-
ASEModel model{{}, sequence, B0.Get()};
213+
ASEModel model{{}, sequence, Hct.Get(), B0.Get()};
196214
ASEFit fit{model};
197215
process(fit);
198216
}

0 commit comments

Comments
 (0)