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