mitkNumericTwoCompartmentExchangeModel.cpp

Go to the documentation of this file.

/*============================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.

Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.

============================================================================*/

#include "mitkNumericTwoCompartmentExchangeModel.h"
#include "mitkAIFParametrizerHelper.h"
#include "mitkTimeGridHelper.h"
#include "mitkTwoCompartmentExchangeModelDifferentialEquations.h"
#include <vnl/algo/vnl_fft_1d.h>
#include <boost/numeric/odeint.hpp>
#include <fstream>

const std::string mitk::NumericTwoCompartmentExchangeModel::MODEL_DISPLAY_NAME =
  "Numeric Two Compartment Exchange Model";

const std::string mitk::NumericTwoCompartmentExchangeModel::NAME_PARAMETER_F = "F";
const std::string mitk::NumericTwoCompartmentExchangeModel::NAME_PARAMETER_PS = "PS";
const std::string mitk::NumericTwoCompartmentExchangeModel::NAME_PARAMETER_ve = "ve";
const std::string mitk::NumericTwoCompartmentExchangeModel::NAME_PARAMETER_vp = "vp";

const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_F = "ml/min/100ml";
const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_PS = "ml/min/100ml";
const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_ve = "ml/ml";
const std::string mitk::NumericTwoCompartmentExchangeModel::UNIT_PARAMETER_vp = "ml/ml";

const unsigned int mitk::NumericTwoCompartmentExchangeModel::POSITION_PARAMETER_F = 0;
const unsigned int mitk::NumericTwoCompartmentExchangeModel::POSITION_PARAMETER_PS = 1;
const unsigned int mitk::NumericTwoCompartmentExchangeModel::POSITION_PARAMETER_ve = 2;
const unsigned int mitk::NumericTwoCompartmentExchangeModel::POSITION_PARAMETER_vp = 3;

const unsigned int mitk::NumericTwoCompartmentExchangeModel::NUMBER_OF_PARAMETERS = 4;

const std::string mitk::NumericTwoCompartmentExchangeModel::NAME_STATIC_PARAMETER_ODEINTStepSize = "ODEIntStepSize";


std::string mitk::NumericTwoCompartmentExchangeModel::GetModelDisplayName() const
{
  return MODEL_DISPLAY_NAME;
};

std::string mitk::NumericTwoCompartmentExchangeModel::GetModelType() const
{
  return "Perfusion.MR";
};


mitk::NumericTwoCompartmentExchangeModel::NumericTwoCompartmentExchangeModel()
{

}

mitk::NumericTwoCompartmentExchangeModel::~NumericTwoCompartmentExchangeModel()
{

}

mitk::NumericTwoCompartmentExchangeModel::ParameterNamesType mitk::NumericTwoCompartmentExchangeModel::GetStaticParameterNames() const
{
  ParameterNamesType result;

  result.push_back(NAME_STATIC_PARAMETER_AIF);
  result.push_back(NAME_STATIC_PARAMETER_AIFTimeGrid);
  result.push_back(NAME_STATIC_PARAMETER_ODEINTStepSize);

  return result;
}

mitk::NumericTwoCompartmentExchangeModel::ParametersSizeType  mitk::NumericTwoCompartmentExchangeModel::GetNumberOfStaticParameters()
const
{
  return 3;
}


void mitk::NumericTwoCompartmentExchangeModel::SetStaticParameter(const ParameterNameType& name,
    const StaticParameterValuesType& values)
{
  if (name == NAME_STATIC_PARAMETER_AIF)
  {
    AterialInputFunctionType aif = mitk::convertParameterToArray(values);

    SetAterialInputFunctionValues(aif);
  }

  if (name == NAME_STATIC_PARAMETER_AIFTimeGrid)
  {
    TimeGridType timegrid = mitk::convertParameterToArray(values);

    SetAterialInputFunctionTimeGrid(timegrid);
  }

  if (name == NAME_STATIC_PARAMETER_ODEINTStepSize)
  {
      SetODEINTStepSize(values[0]);
  }
};

mitk::NumericTwoCompartmentExchangeModel::StaticParameterValuesType mitk::NumericTwoCompartmentExchangeModel::GetStaticParameterValue(
  const ParameterNameType& name) const
{
  StaticParameterValuesType result;

  if (name == NAME_STATIC_PARAMETER_AIF)
  {
    result = mitk::convertArrayToParameter(this->m_AterialInputFunctionValues);
  }

  if (name == NAME_STATIC_PARAMETER_AIFTimeGrid)
  {
    result = mitk::convertArrayToParameter(this->m_AterialInputFunctionTimeGrid);
  }
  if (name == NAME_STATIC_PARAMETER_ODEINTStepSize)
  {
    result.push_back(GetODEINTStepSize());
  }

  return result;
};


mitk::NumericTwoCompartmentExchangeModel::ParameterNamesType
mitk::NumericTwoCompartmentExchangeModel::GetParameterNames() const
{
  ParameterNamesType result;

  result.push_back(NAME_PARAMETER_F);
  result.push_back(NAME_PARAMETER_PS);
  result.push_back(NAME_PARAMETER_ve);
  result.push_back(NAME_PARAMETER_vp);

  return result;
}

mitk::NumericTwoCompartmentExchangeModel::ParametersSizeType
mitk::NumericTwoCompartmentExchangeModel::GetNumberOfParameters() const
{
  return NUMBER_OF_PARAMETERS;
}


mitk::NumericTwoCompartmentExchangeModel::ParamterUnitMapType
mitk::NumericTwoCompartmentExchangeModel::GetParameterUnits() const
{
  ParamterUnitMapType result;

  result.insert(std::make_pair(NAME_PARAMETER_F, UNIT_PARAMETER_F));
  result.insert(std::make_pair(NAME_PARAMETER_PS, UNIT_PARAMETER_PS));
  result.insert(std::make_pair(NAME_PARAMETER_vp, UNIT_PARAMETER_vp));
  result.insert(std::make_pair(NAME_PARAMETER_ve, UNIT_PARAMETER_ve));

  return result;
};

mitk::NumericTwoCompartmentExchangeModel::ModelResultType
mitk::NumericTwoCompartmentExchangeModel::ComputeModelfunction(const ParametersType& parameters)
const
{
  typedef itk::Array<double> ConcentrationCurveType;
  typedef std::vector<double> ConcentrationVectorType;

  if (this->m_TimeGrid.GetSize() == 0)
  {
    itkExceptionMacro("No Time Grid Set! Cannot Calculate Signal");
  }

  AterialInputFunctionType aterialInputFunction;
  aterialInputFunction = GetAterialInputFunction(this->m_TimeGrid);

  unsigned int timeSteps = this->m_TimeGrid.GetSize();

  mitk::TwoCompartmentExchangeModelDifferentialEquations::AIFType aif = mitk::convertArrayToParameter(
        aterialInputFunction);
  mitk::TwoCompartmentExchangeModelDifferentialEquations::AIFType grid =
    mitk::convertArrayToParameter(m_TimeGrid);

  mitk::TwoCompartmentExchangeModelDifferentialEquations::AIFType aifODE = aif;
  aifODE.push_back(aif[timeSteps - 1]);
  mitk::TwoCompartmentExchangeModelDifferentialEquations::AIFType gridODE = grid;
  gridODE.push_back(grid[timeSteps - 1] + (grid[timeSteps - 1] - grid[timeSteps - 2]));



  //Model Parameters
  double F = (double) parameters[POSITION_PARAMETER_F] / 6000.0;
  double PS  = (double) parameters[POSITION_PARAMETER_PS] / 6000.0;
  double ve = (double) parameters[POSITION_PARAMETER_ve];
  double vp = (double) parameters[POSITION_PARAMETER_vp];


  mitk::TwoCompartmentExchangeModelDifferentialEquations ode;
  ode.initialize(F, PS, ve, vp);
  ode.setAIF(aifODE);
  ode.setAIFTimeGrid(gridODE);

  state_type x(2);
  x[0] = 0.0;
  x[1] = 0.0;
  typedef boost::numeric::odeint::runge_kutta_cash_karp54<state_type> error_stepper_type;
  //typedef boost::numeric::odeint::runge_kutta4< state_type > stepper_type;

  ConcentrationVectorType Cp;
  ConcentrationVectorType Ce;
  ConcentrationVectorType odeTimeGrid;

  error_stepper_type stepper;


//  const double dt = 0.05;
   const double dt = this->m_ODEINTStepSize;


  for (double t = 0.0; t < this->m_TimeGrid(this->m_TimeGrid.GetSize() - 1) - 2*dt; t += dt)
  {
    stepper.do_step(ode, x, t, dt);
    Cp.push_back(x[0]);
    Ce.push_back(x[1]);
    odeTimeGrid.push_back(t);
  }

  ConcentrationCurveType plasmaConcentration = mitk::convertParameterToArray(Cp);
  ConcentrationCurveType EESConcentration = mitk::convertParameterToArray(Ce);
  ConcentrationCurveType rungeKuttaTimeGrid = mitk::convertParameterToArray(odeTimeGrid);

  mitk::ModelBase::ModelResultType C_Plasma = mitk::InterpolateSignalToNewTimeGrid(plasmaConcentration, rungeKuttaTimeGrid, m_TimeGrid);
  mitk::ModelBase::ModelResultType C_EES = mitk::InterpolateSignalToNewTimeGrid(EESConcentration, rungeKuttaTimeGrid, m_TimeGrid);


  //Signal that will be returned by ComputeModelFunction
  mitk::ModelBase::ModelResultType signal(timeSteps);
  signal.fill(0.0);

  mitk::ModelBase::ModelResultType::iterator signalPos = signal.begin();
  mitk::ModelBase::ModelResultType::const_iterator CePos = C_EES.begin();

  mitk::ModelBase::ModelResultType::const_iterator t = this->m_TimeGrid.begin();
  mitk::ModelBase::ModelResultType::const_iterator Cin = aterialInputFunction.begin();



  for (mitk::ModelBase::ModelResultType::const_iterator CpPos = C_Plasma.begin();
       CpPos != C_Plasma.end(); ++CpPos, ++CePos, ++signalPos, ++t, ++Cin)
  {
    *signalPos = vp * (*CpPos) + ve * (*CePos);

  }

  return signal;

}




itk::LightObject::Pointer mitk::NumericTwoCompartmentExchangeModel::InternalClone() const
{
  NumericTwoCompartmentExchangeModel::Pointer newClone = NumericTwoCompartmentExchangeModel::New();

  newClone->SetTimeGrid(this->m_TimeGrid);

  return newClone.GetPointer();
}

void mitk::NumericTwoCompartmentExchangeModel::PrintSelf(std::ostream& os,
    ::itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);


}