Raw File
OxCalculatorT2.hxx
/*!
 * \file OxCalculatorT2.hxx
 * \author Konrad Werys
 * \date 2019/11/05
 */

#ifndef Tomato_OXCALCULATORT2_HXX
#define Tomato_OXCALCULATORT2_HXX

namespace Ox {

    template< typename MeasureType >
    int
    CalculatorT2<MeasureType>
    ::calculate(){

        // some initialisation
        this->_Results = std::map <std::string, MeasureType>();

        // get ready, continue if prepareToCalculate EXIT_SUCCESS
        if (this->prepareToCalculate() == 1) {
            return 1; // EXIT_FAILURE
        }

        // calculate if higher than the cutoff
        if (KWUtil::calcMeanArray(this->getNSamples(), this->getSigMag()) < this->getMeanCutOff()) {
            return 0; // EXIT_SUCCESS
        }

        MeasureType tempParameters[3];
        KWUtil::copyArrayToArray(this->_nDims, tempParameters, this->_StartPoint); // start from the starting point

//        MeasureType tempSig[4];
//        KWUtil::copyArrayToArray(3, tempSig, this->getSigMag());
//        tempSig[3] = 0;
//
//        MeasureType tempEcho[4];
//        KWUtil::copyArrayToArray(3, tempEcho, this->getEchoTimes());
//        tempEcho[3] = 10000;
//
//        this->getModel()->setNSamples(4);
//        this->getModel()->setSignal(tempSig);
//        this->getModel()->setEchoTimes(tempEcho);

        this->getModel()->setNSamples(this->getNSamples());
        this->getModel()->setSignal(this->getSigMag());
        this->getModel()->setEchoTimes(this->getEchoTimes());

        // configure Fitter
        this->getFitter()->setParameters(tempParameters);
        this->getFitter()->setModel(this->getModel());

        // fit
        this->getFitter()->performFitting();

        MeasureType tempResults[3];
        KWUtil::copyArrayToArray(this->_nDims, tempResults, this->getFitter()->getParameters());
        if (tempResults[0] != 0) {
            if (this->_nDims == 2) {
                this->_Results["A"]  = tempResults[0];
                this->_Results["T2"] = tempResults[1];
                this->_Results["R2"] = calculateR2AbsFromModel(this->getNSamples(), this->getEchoTimes(),this->getSigMag(), tempResults);
            } else if (this->_nDims == 3){
                this->_Results["A"]  = tempResults[0];
                this->_Results["B"]  = tempResults[1];
                this->_Results["T2"] = tempResults[2];
                this->_Results["R2"] = calculateR2AbsFromModel(this->getNSamples(), this->getEchoTimes(), this->getSigMag(), tempResults);
            }
        }

        return 0; // EXIT_SUCCESS
    }

    template< typename MeasureType >
    int
    CalculatorT2<MeasureType>
    ::prepareToCalculate(){

        // if fitter does not have to iterate, do not calculate
        if (this->getFitter()->getMaxFunctionEvals() == 0){
            return 1; // EXIT_FAILURE
        }

        // verify echoTimes are sorted
        for (int i = 0; i < this->getNSamples()-1; i++){
            if (this->getEchoTimes()[i] > this->getEchoTimes()[i+1]){
                throw std::runtime_error("Echo times have to be sorted!");
            }
        }

        // calculate start point
        MeasureType const defaultStartPoint2[] = {500, 100};
        MeasureType const defaultStartPoint3[] = {5, 100, 50};

        if (!this->getStartPointCalculator()){
            // no start point calculator
            if (this->_nDims == 2) {
                KWUtil::copyArrayToArray(this->_nDims, this->_StartPoint, defaultStartPoint2);
            } else if (this->_nDims == 3) {
                KWUtil::copyArrayToArray(this->_nDims, this->_StartPoint, defaultStartPoint3);
            } else {
                throw std::runtime_error("CalculatorT2: No default start point for this number of dimensions, "
                                         "consider using StartPointCalculator ");
            }
        } else {
            // start point calculator
            this->getStartPointCalculator()->setNDims(this->_nDims);
            if (!this->getStartPointCalculator()->getInputStartPoint()) {
                if (this->_nDims == 2) {
                    this->getStartPointCalculator()->setInputStartPoint(defaultStartPoint2);
                } else if (this->_nDims == 3) {
                    this->getStartPointCalculator()->setInputStartPoint(defaultStartPoint3);
                } else {
                    throw std::runtime_error("CalculatorT2: Set InputStartPoint in StartPointCalculator");
                }
            }
            this->getStartPointCalculator()->setNSamples(this->getNSamples());
            this->getStartPointCalculator()->setEchoTimes(this->getEchoTimes());
            this->getStartPointCalculator()->setSigMag(this->getSigMag());
            this->getStartPointCalculator()->setSigns(this->getSigns());
            this->getStartPointCalculator()->setCalculatedStartPoint(this->_StartPoint);

            this->getStartPointCalculator()->calculateStartPoint();
        }

        return 0; // EXIT_SUCCESS
    }

    template< typename MeasureType >
    MeasureType
    CalculatorT2<MeasureType>
    ::calculateR2AbsFromModel(int nSamples, const MeasureType* times, const MeasureType* signal, const MeasureType* parameters) {

        MeasureType *absFitted  = new MeasureType[nSamples];
        MeasureType *absYsignal = new MeasureType[nSamples];

        for (int i = 0; i < nSamples; i++){
            MeasureType fitted;
            fitted = this->_Model->calcModelValue(parameters, times[i]);
            absFitted[i] = fabs(fitted);
            absYsignal[i] = fabs(signal[i]);
        }

        double result = KWUtil::calcR2cor(nSamples, absFitted, absYsignal);

        delete[] absFitted;
        delete[] absYsignal;
        return result;
    }

    template< typename MeasureType >
    const MeasureType *
    CalculatorT2<MeasureType>
    ::getEchoTimes() const {
        if (!this->_EchoTimes) {
            throw std::runtime_error("_EchoTimes equals 0. Set _EchoTimes");
        };
        return this->_EchoTimes;
    }

} //namespace Ox

#endif //Tomato_OXCALCULATORT2_HXX
back to top