Medical Imaging Interaction Toolkit  2018.4.99-c7ee88da
Medical Imaging Interaction Toolkit
mitkNumericTwoCompartmentExchangeModel.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
15 #include "mitkTimeGridHelper.h"
17 #include <vnl/algo/vnl_fft_1d.h>
18 #include <boost/numeric/odeint.hpp>
19 #include <fstream>
20 
22  "Numeric Two Compartment Exchange Model";
23 
28 
29 const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_F = "ml/min/100ml";
30 const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_PS = "ml/min/100ml";
33 
38 
40 
42 
43 
45 {
46  return MODEL_DISPLAY_NAME;
47 };
48 
50 {
51  return "Perfusion.MR";
52 };
53 
54 
56 {
57 
58 }
59 
61 {
62 
63 }
64 
66 {
67  ParameterNamesType result;
68 
69  result.push_back(NAME_STATIC_PARAMETER_AIF);
70  result.push_back(NAME_STATIC_PARAMETER_AIFTimeGrid);
71  result.push_back(NAME_STATIC_PARAMETER_ODEINTStepSize);
72 
73  return result;
74 }
75 
77 const
78 {
79  return 3;
80 }
81 
82 
84  const StaticParameterValuesType& values)
85 {
86  if (name == NAME_STATIC_PARAMETER_AIF)
87  {
89 
91  }
92 
94  {
96 
98  }
99 
101  {
102  SetODEINTStepSize(values[0]);
103  }
104 };
105 
107  const ParameterNameType& name) const
108 {
110 
111  if (name == NAME_STATIC_PARAMETER_AIF)
112  {
114  }
115 
117  {
119  }
121  {
122  result.push_back(GetODEINTStepSize());
123  }
124 
125  return result;
126 };
127 
128 
131 {
132  ParameterNamesType result;
133 
134  result.push_back(NAME_PARAMETER_F);
135  result.push_back(NAME_PARAMETER_PS);
136  result.push_back(NAME_PARAMETER_ve);
137  result.push_back(NAME_PARAMETER_vp);
138 
139  return result;
140 }
141 
144 {
145  return NUMBER_OF_PARAMETERS;
146 }
147 
148 
151 {
152  ParamterUnitMapType result;
153 
154  result.insert(std::make_pair(NAME_PARAMETER_F, UNIT_PARAMETER_F));
155  result.insert(std::make_pair(NAME_PARAMETER_PS, UNIT_PARAMETER_PS));
156  result.insert(std::make_pair(NAME_PARAMETER_vp, UNIT_PARAMETER_vp));
157  result.insert(std::make_pair(NAME_PARAMETER_ve, UNIT_PARAMETER_ve));
158 
159  return result;
160 };
161 
164 const
165 {
166  typedef itk::Array<double> ConcentrationCurveType;
167  typedef std::vector<double> ConcentrationVectorType;
168 
169  if (this->m_TimeGrid.GetSize() == 0)
170  {
171  itkExceptionMacro("No Time Grid Set! Cannot Calculate Signal");
172  }
173 
174  AterialInputFunctionType aterialInputFunction;
175  aterialInputFunction = GetAterialInputFunction(this->m_TimeGrid);
176 
177  unsigned int timeSteps = this->m_TimeGrid.GetSize();
178 
182  aterialInputFunction);
185 
187  aifODE.push_back(aif[timeSteps - 1]);
189  gridODE.push_back(grid[timeSteps - 1] + (grid[timeSteps - 1] - grid[timeSteps - 2]));
190 
191 
192 
193  //Model Parameters
194  double F = (double) parameters[POSITION_PARAMETER_F] / 6000.0;
195  double PS = (double) parameters[POSITION_PARAMETER_PS] / 6000.0;
196  double ve = (double) parameters[POSITION_PARAMETER_ve];
197  double vp = (double) parameters[POSITION_PARAMETER_vp];
198 
199 
202  ode.initialize(F, PS, ve, vp);
203  ode.setAIF(aifODE);
204  ode.setAIFTimeGrid(gridODE);
205 
206  state_type x(2);
207  x[0] = 0.0;
208  x[1] = 0.0;
209  typedef boost::numeric::odeint::runge_kutta_cash_karp54<state_type> error_stepper_type;
210  //typedef boost::numeric::odeint::runge_kutta4< state_type > stepper_type;
211 
213  ConcentrationVectorType Cp;
214  ConcentrationVectorType Ce;
215  ConcentrationVectorType odeTimeGrid;
216 
217  error_stepper_type stepper;
218 
219 
221 // const double dt = 0.05;
222  const double dt = this->m_ODEINTStepSize;
223 
224 
226  for (double t = 0.0; t < this->m_TimeGrid(this->m_TimeGrid.GetSize() - 1) - 2*dt; t += dt)
227  {
228  stepper.do_step(ode, x, t, dt);
229  Cp.push_back(x[0]);
230  Ce.push_back(x[1]);
231  odeTimeGrid.push_back(t);
232  }
233 
235  ConcentrationCurveType plasmaConcentration = mitk::convertParameterToArray(Cp);
236  ConcentrationCurveType EESConcentration = mitk::convertParameterToArray(Ce);
237  ConcentrationCurveType rungeKuttaTimeGrid = mitk::convertParameterToArray(odeTimeGrid);
238 
239  mitk::ModelBase::ModelResultType C_Plasma = mitk::InterpolateSignalToNewTimeGrid(plasmaConcentration, rungeKuttaTimeGrid, m_TimeGrid);
240  mitk::ModelBase::ModelResultType C_EES = mitk::InterpolateSignalToNewTimeGrid(EESConcentration, rungeKuttaTimeGrid, m_TimeGrid);
241 
242 
243  //Signal that will be returned by ComputeModelFunction
244  mitk::ModelBase::ModelResultType signal(timeSteps);
245  signal.fill(0.0);
246 
247  mitk::ModelBase::ModelResultType::iterator signalPos = signal.begin();
248  mitk::ModelBase::ModelResultType::const_iterator CePos = C_EES.begin();
249 
250  mitk::ModelBase::ModelResultType::const_iterator t = this->m_TimeGrid.begin();
251  mitk::ModelBase::ModelResultType::const_iterator Cin = aterialInputFunction.begin();
252 
253 
254 
255  for (mitk::ModelBase::ModelResultType::const_iterator CpPos = C_Plasma.begin();
256  CpPos != C_Plasma.end(); ++CpPos, ++CePos, ++signalPos, ++t, ++Cin)
257  {
258  *signalPos = vp * (*CpPos) + ve * (*CePos);
259 
260  }
261 
262  return signal;
263 
264 }
265 
266 
267 
268 
270 {
272 
273  newClone->SetTimeGrid(this->m_TimeGrid);
274 
275  return newClone.GetPointer();
276 }
277 
279  ::itk::Indent indent) const
280 {
281  Superclass::PrintSelf(os, indent);
282 
283 
284 }
285 
286 
ModelTraitsInterface::ParametersType ParametersType
Definition: mitkModelBase.h:59
TimeGridType m_TimeGrid
virtual const double & GetODEINTStepSize()
static const std::string NAME_STATIC_PARAMETER_AIFTimeGrid
ModelResultType ComputeModelfunction(const ParametersType &parameters) const override
void PrintSelf(std::ostream &os, ::itk::Indent indent) const override
itk::Array< double > TimeGridType
Definition: mitkModelBase.h:62
Helper Class for NumericTwoCompartmentExchangeModel: Defines the differential equations (Mass Balance...
const AterialInputFunctionType GetAterialInputFunction(TimeGridType currentTimeGrid) const
MITKPHARMACOKINETICS_EXPORT ModelBase::StaticParameterValuesType convertArrayToParameter(itk::Array< double > array)
itk::LightObject::Pointer InternalClone() const override
ModelTraitsInterface::ParametersSizeType ParametersSizeType
Definition: mitkModelBase.h:65
AterialInputFunctionType m_AterialInputFunctionValues
std::vector< StaticParameterValueType > StaticParameterValuesType
Definition: mitkModelBase.h:71
ModelTraitsInterface::ModelResultType ModelResultType
Definition: mitkModelBase.h:55
virtual void SetODEINTStepSize(double _arg)
MITKMODELFIT_EXPORT ModelBase::ModelResultType InterpolateSignalToNewTimeGrid(const ModelBase::ModelResultType &inputSignal, const ModelBase::TimeGridType &inputGrid, const ModelBase::TimeGridType &outputGrid)
MITKPHARMACOKINETICS_EXPORT itk::Array< double > convertParameterToArray(ModelBase::StaticParameterValuesType)
itk::Array< double > AterialInputFunctionType
void PrintSelf(std::ostream &os, ::itk::Indent indent) const override
ModelTraitsInterface::ParameterNameType ParameterNameType
Definition: mitkModelBase.h:63
void SetStaticParameter(const ParameterNameType &name, const StaticParameterValuesType &values) override
virtual void SetAterialInputFunctionTimeGrid(TimeGridType _arg)
static const std::string NAME_STATIC_PARAMETER_AIF
virtual void SetAterialInputFunctionValues(AterialInputFunctionType _arg)
ModelTraitsInterface::ParameterNamesType ParameterNamesType
Definition: mitkModelBase.h:64
void initialize(double Fp, double ps, double fi, double fp)
Initialize class with parameters F/Vp, PS/Vp, fi and fp that are free fit parameters.
TimeGridType m_AterialInputFunctionTimeGrid
std::map< ParameterNameType, std::string > ParamterUnitMapType
StaticParameterValuesType GetStaticParameterValue(const ParameterNameType &name) const override