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