diff --git a/Modules/ModelFit/src/Common/mitkTimeGridHelper.cpp b/Modules/ModelFit/src/Common/mitkTimeGridHelper.cpp index f422b9ccb7..685cef5d67 100644 --- a/Modules/ModelFit/src/Common/mitkTimeGridHelper.cpp +++ b/Modules/ModelFit/src/Common/mitkTimeGridHelper.cpp @@ -1,94 +1,93 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTimeGridHelper.h" #include "itkExceptionObject.h" bool mitk::TimeGridIsMonotonIncreasing(const mitk::ModelBase::TimeGridType timeGrid) { - mitk::ModelBase::TimeGridType::ValueType lastTime = itk::NumericTraits::NonpositiveMin(); - const auto beginPos = timeGrid.begin(); - for(mitk::ModelBase::TimeGridType::const_iterator posTime = beginPos; posTime != timeGrid.end(); ++posTime) + const auto endPos = timeGrid.end(); + for(mitk::ModelBase::TimeGridType::const_iterator posTime = beginPos; posTime != endPos; ++posTime) { if (posTime != beginPos && *(posTime-1)<*posTime) return false; } return true; }; mitk::ModelBase::ModelResultType mitk::InterpolateSignalToNewTimeGrid(const ModelBase::ModelResultType& inputSignal, const ModelBase::TimeGridType& inputGrid, const ModelBase::TimeGridType& outputGrid) { mitk::ModelBase::ModelResultType result(outputGrid.GetSize()); if (! inputSignal.GetSize()) { return result; } if (inputSignal.GetSize() != inputGrid.GetSize()) { itkGenericExceptionMacro("Input signal and input time grid have not the same size."); } mitk::ModelBase::ModelResultType::ValueType lastValue = inputSignal[0]; mitk::ModelBase::TimeGridType::ValueType lastTime = itk::NumericTraits::NonpositiveMin(); mitk::ModelBase::TimeGridType::const_iterator posITime = inputGrid.begin(); mitk::ModelBase::ModelResultType::const_iterator posValue = inputSignal.begin(); mitk::ModelBase::ModelResultType::iterator posResult = result.begin(); for(mitk::ModelBase::TimeGridType::const_iterator posOTime = outputGrid.begin(); posOTime != outputGrid.end(); ++posResult, ++posOTime) { while(*posOTime > *posITime && posITime!=inputGrid.end()) { //forward in the input grid until the current output point //is between last and the current input point. lastValue = *posValue; lastTime = *posITime; ++posValue; ++posITime; } double weightLast = 1 - (*posOTime - lastTime)/(*posITime - lastTime); double weightNext = 1 - (*posITime - *posOTime)/(*posITime - lastTime); *posResult = weightLast * lastValue + weightNext * (*posValue); } return result; }; mitk::ModelBase::TimeGridType mitk::GenerateSupersampledTimeGrid(const mitk::ModelBase::TimeGridType& grid, const unsigned int samplingRate) { unsigned int origGridSize = grid.size(); mitk::ModelBase::TimeGridType interpolatedTimeGrid(((origGridSize - 1) * samplingRate) + 1); for (unsigned int t = 0; t < origGridSize - 1; ++t) { double delta = (grid[t + 1] - grid[t]) / samplingRate; for (unsigned int i = 0; i < samplingRate; ++i) { interpolatedTimeGrid[(t * samplingRate) + i] = grid[t] + i * delta; } } interpolatedTimeGrid[interpolatedTimeGrid.size() - 1] = grid[grid.size() - 1]; return interpolatedTimeGrid; -}; \ No newline at end of file +}; diff --git a/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentFDGModel.cpp b/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentFDGModel.cpp index 550b3f17a7..066b087bf9 100644 --- a/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentFDGModel.cpp +++ b/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentFDGModel.cpp @@ -1,162 +1,160 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTwoTissueCompartmentFDGModel.h" #include "mitkConvolutionHelper.h" #include const std::string mitk::TwoTissueCompartmentFDGModel::MODEL_DISPLAY_NAME = "Two Tissue Compartment Model for FDG (Sokoloff Model)"; const std::string mitk::TwoTissueCompartmentFDGModel::NAME_PARAMETER_K1 = "K1"; const std::string mitk::TwoTissueCompartmentFDGModel::NAME_PARAMETER_k2 = "k2"; const std::string mitk::TwoTissueCompartmentFDGModel::NAME_PARAMETER_k3 = "k3"; const std::string mitk::TwoTissueCompartmentFDGModel::NAME_PARAMETER_VB = "V_B"; const std::string mitk::TwoTissueCompartmentFDGModel::UNIT_PARAMETER_K1 = "1/min"; const std::string mitk::TwoTissueCompartmentFDGModel::UNIT_PARAMETER_k2 = "1/min"; const std::string mitk::TwoTissueCompartmentFDGModel::UNIT_PARAMETER_k3 = "1/min"; const std::string mitk::TwoTissueCompartmentFDGModel::UNIT_PARAMETER_VB = "ml/ml"; const unsigned int mitk::TwoTissueCompartmentFDGModel::POSITION_PARAMETER_K1 = 0; const unsigned int mitk::TwoTissueCompartmentFDGModel::POSITION_PARAMETER_k2 = 1; const unsigned int mitk::TwoTissueCompartmentFDGModel::POSITION_PARAMETER_k3 = 2; const unsigned int mitk::TwoTissueCompartmentFDGModel::POSITION_PARAMETER_VB = 3; const unsigned int mitk::TwoTissueCompartmentFDGModel::NUMBER_OF_PARAMETERS = 4; std::string mitk::TwoTissueCompartmentFDGModel::GetModelDisplayName() const { return MODEL_DISPLAY_NAME; }; std::string mitk::TwoTissueCompartmentFDGModel::GetModelType() const { return "Dynamic.PET"; }; mitk::TwoTissueCompartmentFDGModel::TwoTissueCompartmentFDGModel() { } mitk::TwoTissueCompartmentFDGModel::~TwoTissueCompartmentFDGModel() { } mitk::TwoTissueCompartmentFDGModel::ParameterNamesType mitk::TwoTissueCompartmentFDGModel::GetParameterNames() const { ParameterNamesType result; result.push_back(NAME_PARAMETER_K1); result.push_back(NAME_PARAMETER_k2); result.push_back(NAME_PARAMETER_k3); result.push_back(NAME_PARAMETER_VB); return result; } mitk::TwoTissueCompartmentFDGModel::ParametersSizeType mitk::TwoTissueCompartmentFDGModel::GetNumberOfParameters() const { return NUMBER_OF_PARAMETERS; } mitk::TwoTissueCompartmentFDGModel::ParamterUnitMapType mitk::TwoTissueCompartmentFDGModel::GetParameterUnits() const { ParamterUnitMapType result; result.insert(std::make_pair(NAME_PARAMETER_K1, UNIT_PARAMETER_K1)); result.insert(std::make_pair(NAME_PARAMETER_k2, UNIT_PARAMETER_k2)); result.insert(std::make_pair(NAME_PARAMETER_k3, UNIT_PARAMETER_k3)); result.insert(std::make_pair(NAME_PARAMETER_VB, UNIT_PARAMETER_VB)); return result; }; mitk::TwoTissueCompartmentFDGModel::ModelResultType mitk::TwoTissueCompartmentFDGModel::ComputeModelfunction(const ParametersType& parameters) const { - typedef itk::Array ResidueFunctionType; - if (this->m_TimeGrid.GetSize() == 0) { itkExceptionMacro("No Time Grid Set! Cannot Calculate Signal"); } AterialInputFunctionType aterialInputFunction; aterialInputFunction = GetAterialInputFunction(this->m_TimeGrid); unsigned int timeSteps = this->m_TimeGrid.GetSize(); //Model Parameters double k1 = (double)parameters[POSITION_PARAMETER_K1] / 60.0; double k2 = (double)parameters[POSITION_PARAMETER_k2] / 60.0; double k3 = (double)parameters[POSITION_PARAMETER_k3] / 60.0; double VB = parameters[POSITION_PARAMETER_VB]; double lambda = k2+k3; //double lambda2 = -alpha2; mitk::ModelBase::ModelResultType exp = mitk::convoluteAIFWithExponential(this->m_TimeGrid, aterialInputFunction, lambda); mitk::ModelBase::ModelResultType CA = mitk::convoluteAIFWithConstant(this->m_TimeGrid, aterialInputFunction, k3); //Signal that will be returned by ComputeModelFunction mitk::ModelBase::ModelResultType signal(timeSteps); signal.fill(0.0); mitk::ModelBase::ModelResultType::const_iterator expPos = exp.begin(); mitk::ModelBase::ModelResultType::const_iterator CAPos = CA.begin(); AterialInputFunctionType::const_iterator aifPos = aterialInputFunction.begin(); for (mitk::ModelBase::ModelResultType::iterator signalPos = signal.begin(); signalPos != signal.end(); ++expPos, ++signalPos, ++aifPos) { double Ci = k1 * k2 /lambda *(*expPos) + k1*k3/lambda*(*CAPos); *signalPos = VB * (*aifPos) + (1 - VB) * Ci; } return signal; } itk::LightObject::Pointer mitk::TwoTissueCompartmentFDGModel::InternalClone() const { TwoTissueCompartmentFDGModel::Pointer newClone = TwoTissueCompartmentFDGModel::New(); newClone->SetTimeGrid(this->m_TimeGrid); return newClone.GetPointer(); } void mitk::TwoTissueCompartmentFDGModel::PrintSelf(std::ostream& os, ::itk::Indent indent) const { Superclass::PrintSelf(os, indent); } diff --git a/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentModel.cpp b/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentModel.cpp index c314fac87c..ade10da373 100644 --- a/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentModel.cpp +++ b/Modules/Pharmacokinetics/src/Models/mitkTwoTissueCompartmentModel.cpp @@ -1,178 +1,176 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTwoTissueCompartmentModel.h" #include "mitkConvolutionHelper.h" #include const std::string mitk::TwoTissueCompartmentModel::MODEL_DISPLAY_NAME = "Two Tissue Compartment Model"; const std::string mitk::TwoTissueCompartmentModel::NAME_PARAMETER_K1 = "K1"; const std::string mitk::TwoTissueCompartmentModel::NAME_PARAMETER_k2 = "k2"; const std::string mitk::TwoTissueCompartmentModel::NAME_PARAMETER_k3 = "k3"; const std::string mitk::TwoTissueCompartmentModel::NAME_PARAMETER_k4 = "k4"; const std::string mitk::TwoTissueCompartmentModel::NAME_PARAMETER_VB = "V_B"; const std::string mitk::TwoTissueCompartmentModel::UNIT_PARAMETER_K1 = "1/min"; const std::string mitk::TwoTissueCompartmentModel::UNIT_PARAMETER_k2 = "1/min"; const std::string mitk::TwoTissueCompartmentModel::UNIT_PARAMETER_k3 = "1/min"; const std::string mitk::TwoTissueCompartmentModel::UNIT_PARAMETER_k4 = "1/min"; const std::string mitk::TwoTissueCompartmentModel::UNIT_PARAMETER_VB = "ml/ml"; const unsigned int mitk::TwoTissueCompartmentModel::POSITION_PARAMETER_K1 = 0; const unsigned int mitk::TwoTissueCompartmentModel::POSITION_PARAMETER_k2 = 1; const unsigned int mitk::TwoTissueCompartmentModel::POSITION_PARAMETER_k3 = 2; const unsigned int mitk::TwoTissueCompartmentModel::POSITION_PARAMETER_k4 = 3; const unsigned int mitk::TwoTissueCompartmentModel::POSITION_PARAMETER_VB = 4; const unsigned int mitk::TwoTissueCompartmentModel::NUMBER_OF_PARAMETERS = 5; inline double square(double a) { return a * a; } std::string mitk::TwoTissueCompartmentModel::GetModelDisplayName() const { return MODEL_DISPLAY_NAME; }; std::string mitk::TwoTissueCompartmentModel::GetModelType() const { return "Dynamic.PET"; }; mitk::TwoTissueCompartmentModel::TwoTissueCompartmentModel() { } mitk::TwoTissueCompartmentModel::~TwoTissueCompartmentModel() { } mitk::TwoTissueCompartmentModel::ParameterNamesType mitk::TwoTissueCompartmentModel::GetParameterNames() const { ParameterNamesType result; result.push_back(NAME_PARAMETER_K1); result.push_back(NAME_PARAMETER_k2); result.push_back(NAME_PARAMETER_k3); result.push_back(NAME_PARAMETER_k4); result.push_back(NAME_PARAMETER_VB); return result; } mitk::TwoTissueCompartmentModel::ParametersSizeType mitk::TwoTissueCompartmentModel::GetNumberOfParameters() const { return NUMBER_OF_PARAMETERS; } mitk::TwoTissueCompartmentModel::ParamterUnitMapType mitk::TwoTissueCompartmentModel::GetParameterUnits() const { ParamterUnitMapType result; result.insert(std::make_pair(NAME_PARAMETER_K1, UNIT_PARAMETER_K1)); result.insert(std::make_pair(NAME_PARAMETER_k2, UNIT_PARAMETER_k2)); result.insert(std::make_pair(NAME_PARAMETER_k3, UNIT_PARAMETER_k3)); result.insert(std::make_pair(NAME_PARAMETER_k4, UNIT_PARAMETER_k4)); result.insert(std::make_pair(NAME_PARAMETER_VB, UNIT_PARAMETER_VB)); return result; }; mitk::TwoTissueCompartmentModel::ModelResultType mitk::TwoTissueCompartmentModel::ComputeModelfunction(const ParametersType& parameters) const { - typedef itk::Array ResidueFunctionType; - if (this->m_TimeGrid.GetSize() == 0) { itkExceptionMacro("No Time Grid Set! Cannot Calculate Signal"); } AterialInputFunctionType aterialInputFunction; aterialInputFunction = GetAterialInputFunction(this->m_TimeGrid); unsigned int timeSteps = this->m_TimeGrid.GetSize(); //Model Parameters double k1 = (double)parameters[POSITION_PARAMETER_K1] / 60.0; double k2 = (double)parameters[POSITION_PARAMETER_k2] / 60.0; double k3 = (double)parameters[POSITION_PARAMETER_k3] / 60.0; double k4 = (double)parameters[POSITION_PARAMETER_k4] / 60.0; double VB = parameters[POSITION_PARAMETER_VB]; double alpha1 = 0.5 * ((k2 + k3 + k4) - sqrt(square(k2 + k3 + k4) - 4 * k2 * k4)); double alpha2 = 0.5 * ((k2 + k3 + k4) + sqrt(square(k2 + k3 + k4) - 4 * k2 * k4)); //double lambda1 = -alpha1; //double lambda2 = -alpha2; mitk::ModelBase::ModelResultType exp1 = mitk::convoluteAIFWithExponential(this->m_TimeGrid, aterialInputFunction, alpha1); mitk::ModelBase::ModelResultType exp2 = mitk::convoluteAIFWithExponential(this->m_TimeGrid, aterialInputFunction, alpha2); //Signal that will be returned by ComputeModelFunction mitk::ModelBase::ModelResultType signal(timeSteps); signal.fill(0.0); mitk::ModelBase::ModelResultType::const_iterator exp1Pos = exp1.begin(); mitk::ModelBase::ModelResultType::const_iterator exp2Pos = exp2.begin(); AterialInputFunctionType::const_iterator aifPos = aterialInputFunction.begin(); for (mitk::ModelBase::ModelResultType::iterator signalPos = signal.begin(); signalPos != signal.end(); ++exp1Pos, ++exp2Pos, ++signalPos, ++aifPos) { double Ci = k1 / (alpha2 - alpha1) * ((k4 - alpha1 + k3) * (*exp1Pos) + (alpha2 - k4 - k3) * (*exp2Pos)); *signalPos = VB * (*aifPos) + (1 - VB) * Ci; } return signal; } itk::LightObject::Pointer mitk::TwoTissueCompartmentModel::InternalClone() const { TwoTissueCompartmentModel::Pointer newClone = TwoTissueCompartmentModel::New(); newClone->SetTimeGrid(this->m_TimeGrid); return newClone.GetPointer(); } void mitk::TwoTissueCompartmentModel::PrintSelf(std::ostream& os, ::itk::Indent indent) const { Superclass::PrintSelf(os, indent); }