mitkConvertToConcentrationViaT1Functor.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 mitkConvertToConcentrationViaT1Functor_h
#define mitkConvertToConcentrationViaT1Functor_h
 
#include "itkMath.h"
 
#include "MitkPharmacokineticsExports.h"
 
namespace mitk {
 
    template <class TInputPixel1, class TInputPixel2, class TInputPixel3, class TOutputpixel>
    class MITKPHARMACOKINETICS_EXPORT ConvertToConcentrationViaT1CalcFunctor
    {
 
    public:
        ConvertToConcentrationViaT1CalcFunctor(): m_relaxivity(0.0), m_TR(0.0),  m_flipangle(0.0) {};
                ~ConvertToConcentrationViaT1CalcFunctor() {};
 
        void initialize(double relaxivity, double TR, double flipangle)
        {
 
                        m_relaxivity = relaxivity;
                        m_TR = TR;
            m_flipangle = flipangle;
        }
 
        bool operator!=( const ConvertToConcentrationViaT1CalcFunctor & other) const
        {
            return !(*this == other);
 
        }
       bool operator==( const ConvertToConcentrationViaT1CalcFunctor & other) const
        {
           return (this->m_relaxivity == other.m_relaxivity) && (this->m_TR == other.m_TR) && (this->m_flipangle == other.m_flipangle);
        }
 
 
        inline TOutputpixel operator()( const TInputPixel1 & value, const TInputPixel2 & baseline, const TInputPixel3 & nativeT1)
        {
            TOutputpixel concentration(0);
            double R10, R1;
            if (baseline !=0 && nativeT1 != 0 && value != 0)
            {
                double s =  (double) value/baseline;
                R10 = (double) 1/nativeT1;
                double tmp1 = log(1-s + s*exp(-R10*m_TR) - exp(-R10*m_TR)* cos(m_flipangle));
                double tmp2 = (1-s*cos(m_flipangle) + s*exp(-R10*m_TR)*cos(m_flipangle) - exp(-R10*m_TR)* cos(m_flipangle));
 
                R1 = (double) -1/m_TR * tmp1/tmp2 ;
 
                concentration = (double) (R1 - R10)/ m_relaxivity;
            }
            else
            {
                concentration = 0;
            }
 
                return concentration;
        }
    private:
                double m_relaxivity;
                double m_TR;
        double m_flipangle;
 
    };
 
}
#endif