Medical Imaging Interaction Toolkit  2018.4.99-a3d2e8fb
Medical Imaging Interaction Toolkit
mitkNumericTwoTissueCompartmentModel.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 Tissue Compartment Model";
23 
29 
35 
41 
43 
44 
46 {
47  return MODEL_DISPLAY_NAME;
48 };
49 
51 {
52  return "Dynamic.PET";
53 };
54 
56 {
57 
58 }
59 
61 {
62 
63 }
64 
67 {
68  ParameterNamesType result;
69 
70  result.push_back(NAME_PARAMETER_K1);
71  result.push_back(NAME_PARAMETER_k2);
72  result.push_back(NAME_PARAMETER_k3);
73  result.push_back(NAME_PARAMETER_k4);
74  result.push_back(NAME_PARAMETER_VB);
75 
76  return result;
77 }
78 
81 {
82  return NUMBER_OF_PARAMETERS;
83 }
84 
87 {
88  ParamterUnitMapType result;
89 
90  result.insert(std::make_pair(NAME_PARAMETER_K1, UNIT_PARAMETER_K1));
91  result.insert(std::make_pair(NAME_PARAMETER_k2, UNIT_PARAMETER_k2));
92  result.insert(std::make_pair(NAME_PARAMETER_k3, UNIT_PARAMETER_k3));
93  result.insert(std::make_pair(NAME_PARAMETER_k4, UNIT_PARAMETER_k4));
94  result.insert(std::make_pair(NAME_PARAMETER_VB, UNIT_PARAMETER_VB));
95 
96  return result;
97 };
98 
99 
102 {
103  typedef itk::Array<double> ConcentrationCurveType;
104  typedef std::vector<double> ConcentrationVectorType;
105 
106  if (this->m_TimeGrid.GetSize() == 0)
107  {
108  itkExceptionMacro("No Time Grid Set! Cannot Calculate Signal");
109  }
110 
111  AterialInputFunctionType aterialInputFunction;
112  aterialInputFunction = GetAterialInputFunction(this->m_TimeGrid);
113 
114  unsigned int timeSteps = this->m_TimeGrid.GetSize();
115 
119  aterialInputFunction);
121  m_TimeGrid);
122 
124  aifODE.push_back(aif[timeSteps - 1]);
126  gridODE.push_back(grid[timeSteps - 1] + (grid[timeSteps - 1] - grid[timeSteps - 2]));
127 
128  //Model Parameters
129  double K1 = (double)parameters[POSITION_PARAMETER_K1] / 60.0;
130  double k2 = (double)parameters[POSITION_PARAMETER_k2] / 60.0;
131  double k3 = (double)parameters[POSITION_PARAMETER_k3] / 60.0;
132  double k4 = (double)parameters[POSITION_PARAMETER_k4] / 60.0;
133  double VB = parameters[POSITION_PARAMETER_VB];
134 
135 
138  ode.initialize(K1, k2, k3, k4);
139  ode.setAIF(aifODE);
140  ode.setAIFTimeGrid(gridODE);
141 
142  state_type x(2);
143  x[0] = 0.0;
144  x[1] = 0.0;
145  typedef boost::numeric::odeint::runge_kutta_cash_karp54<state_type> error_stepper_type;
146  //typedef boost::numeric::odeint::runge_kutta4< state_type > stepper_type;
147 
149  ConcentrationVectorType C1;
150  ConcentrationVectorType C2;
151  ConcentrationVectorType odeTimeGrid;
152 
153  error_stepper_type stepper;
155  const double dt = 0.1;
158  double T = this->m_TimeGrid(timeSteps - 1) + (grid[timeSteps - 1] - grid[timeSteps - 2]);
159 
160  for (double t = 0.0; t < T; t += dt)
161  {
162  stepper.do_step(ode, x, t, dt);
163  C1.push_back(x[0]);
164  C2.push_back(x[1]);
165  odeTimeGrid.push_back(t);
166 
167  }
168 
170  ConcentrationCurveType ConcentrationCompartment1 = mitk::convertParameterToArray(C1);
171  ConcentrationCurveType ConcentrationCompartment2 = mitk::convertParameterToArray(C2);
172  ConcentrationCurveType rungeKuttaTimeGrid = mitk::convertParameterToArray(odeTimeGrid);
173 
175  ConcentrationCompartment1, rungeKuttaTimeGrid, m_TimeGrid);
177  ConcentrationCompartment2, rungeKuttaTimeGrid, m_TimeGrid);
178 
179 
180  //Signal that will be returned by ComputeModelFunction
181  mitk::ModelBase::ModelResultType signal(timeSteps);
182  signal.fill(0.0);
183  mitk::ModelBase::ModelResultType::iterator signalPos = signal.begin();
184  mitk::ModelBase::ModelResultType::const_iterator C1Pos = C_1.begin();
185  mitk::ModelBase::ModelResultType::const_iterator C2Pos = C_2.begin();
186 
187 
188  for (AterialInputFunctionType::const_iterator aifpos = aterialInputFunction.begin();
189  aifpos != aterialInputFunction.end(); ++aifpos, ++C1Pos, ++C2Pos, ++signalPos)
190  {
191  *signalPos = VB * (*aifpos) + (1 - VB) * (*C1Pos + *C2Pos);
192  }
193 
194  return signal;
195 
196 }
197 
199 {
201 
202  newClone->SetTimeGrid(this->m_TimeGrid);
203 
204  return newClone.GetPointer();
205 }
206 
207 void mitk::NumericTwoTissueCompartmentModel::PrintSelf(std::ostream& os, ::itk::Indent indent) const
208 {
209  Superclass::PrintSelf(os, indent);
210 
211 
212 }
213 
214 
215 
216 
void initialize(double k_1, double k_2, double k_3, double k_4)
Initialize class with parameters K1, k2, k3 and k4 that are free fit parameters.
ModelTraitsInterface::ParametersType ParametersType
Definition: mitkModelBase.h:59
TimeGridType m_TimeGrid
void PrintSelf(std::ostream &os, ::itk::Indent indent) const override
void PrintSelf(std::ostream &os, ::itk::Indent indent) const override
const AterialInputFunctionType GetAterialInputFunction(TimeGridType currentTimeGrid) const
MITKPHARMACOKINETICS_EXPORT ModelBase::StaticParameterValuesType convertArrayToParameter(itk::Array< double > array)
ModelTraitsInterface::ParametersSizeType ParametersSizeType
Definition: mitkModelBase.h:65
ModelTraitsInterface::ModelResultType ModelResultType
Definition: mitkModelBase.h:55
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
ModelResultType ComputeModelfunction(const ParametersType &parameters) const override
ModelTraitsInterface::ParameterNamesType ParameterNamesType
Definition: mitkModelBase.h:64
Helper Class for NumericTwoTissueCompartment Model: Defines the differential equations (Mass Balance ...
std::map< ParameterNameType, std::string > ParamterUnitMapType
itk::LightObject::Pointer InternalClone() const override