mitkTwoCompartmentExchangeModelDifferentialEquations.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 MITKTWOCOMPARTMENTEXCHANGEMODELDIFFERENTIALEQUATIONS_H
#define MITKTWOCOMPARTMENTEXCHANGEMODELDIFFERENTIALEQUATIONS_H

#include "mitkNumericTwoCompartmentExchangeModel.h"

namespace mitk{
class TwoCompartmentExchangeModelDifferentialEquations
{
public:

    typedef std::vector< double > AIFType;

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

//        dxdt[0] = -(this->FVP + this->PSVP)*x[0] - this->PSVP*x[1]+this->FVP*Ca_t;
        dxdt[0] = (1/this->vp) * ( this->F*(Ca_t - x[0]) - this->PS*(x[0] - x[1]) );
        dxdt[1] = (1/this->ve) * this->PS * (x[0] - x[1]);

    }

    TwoCompartmentExchangeModelDifferentialEquations() : F(0), PS(0), ve(0), vp(0), m_AIF(0), m_AIFTimeGrid(0)
    {
    }

    void initialize(double Fp, double ps, double fi, double fp)
    {
        this->F = Fp;
        this->PS = ps;
        this->ve = fi;
        this->vp = fp;
    }


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


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

private:

    double F;
    double PS;
    double ve;
    double vp;

    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 // MITKTWOCOMPARTMENTEXCHANGEMODELDIFFERENTIALEQUATIONS_H