mitkTwoTissueCompartmentModelDifferentialEquations.h

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.

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

#ifndef MITKTWOTISSUECOMPARTMENTMODELDIFFERENTIALEQUATIONS_H
#define MITKTWOTISSUECOMPARTMENTMODELDIFFERENTIALEQUATIONS_H
#include "mitkNumericTwoTissueCompartmentModel.h"

namespace mitk{
class TwoTissueCompartmentModelDifferentialEquations
{
public:

    typedef std::vector< double > AIFType;

    void operator() (const mitk::NumericTwoTissueCompartmentModel::state_type &x, mitk::NumericTwoTissueCompartmentModel::state_type &dxdt, const double t)
    {
        double Ca_t = InterpolateAIFToCurrentTimeStep(t);

        dxdt[0] = this->K1*Ca_t-(this->k2+this->k3)*x[0] + this->k4*x[1];
        dxdt[1] = this->k3*x[0] - this->k4*x[1];
    }

    TwoTissueCompartmentModelDifferentialEquations() : K1(0), k2(0), k3(0), k4(0), m_AIF(0), m_AIFTimeGrid(0)
    {
    }

    void initialize(double k_1, double k_2, double k_3, double k_4)
    {
        this->K1 = k_1;
        this->k2 = k_2;
        this->k3 = k_3;
        this->k4 = k_4;
    }


    void setAIF(AIFType &aif)
    {
        this->m_AIF = aif;
    }


    void setAIFTimeGrid(AIFType &grid)
    {
        this->m_AIFTimeGrid = grid;
    }

private:

    double K1;
    double k2;
    double k3;
    double k4;

    AIFType m_AIF;
    AIFType m_AIFTimeGrid;


    double InterpolateAIFToCurrentTimeStep(double t)
    {
        double lastValue = m_AIF[0];
        double lastTime = std::numeric_limits<double>::min();

        AIFType::const_iterator posITime = m_AIFTimeGrid.begin();
        AIFType::const_iterator posValue = m_AIF.begin();

        while(t > *posITime)
        {
            lastValue = *posValue;
            lastTime = *posITime;
            ++posValue;
            ++posITime;
        }
        double weightLast = 1 - (t - lastTime)/(*posITime - lastTime);
        double weightNext = 1- (*posITime - t)/(*posITime - lastTime);
        double result = weightLast * lastValue + weightNext * (*posValue);

        return result;
    }

};
}

#endif // MITKTWOTISSUECOMPARTMENTMODELDIFFERENTIALEQUATIONS_H