diff --git a/CLI/PkModeling.cxx b/CLI/PkModeling.cxx index 7e010ec..ab20554 100644 --- a/CLI/PkModeling.cxx +++ b/CLI/PkModeling.cxx @@ -372,6 +372,7 @@ type Get##name(itk::MetaDataDictionary& dictionary) \ converter->SetconstantBAT(ConstantBAT); converter->SetRGD_relaxivity(RelaxivityValue); converter->SetS0GradThresh(S0GradValue); + converter->SetToftsIntegrationMethod(ToftsIntegrationMethod); if (T1MapFileName != "") { @@ -425,6 +426,7 @@ type Get##name(itk::MetaDataDictionary& dictionary) \ quantifier->Sethematocrit(Hematocrit); quantifier->SetconstantBAT(ConstantBAT); quantifier->SetBATCalculationMode(BATCalculationMode); + quantifier->SetToftsIntegrationMethod(ToftsIntegrationMethod); if (ROIMaskFileName != "") { quantifier->SetROIMask(roiMaskVolume); diff --git a/CLI/PkModeling.xml b/CLI/PkModeling.xml index 85b153e..1cddf0f 100644 --- a/CLI/PkModeling.xml +++ b/CLI/PkModeling.xml @@ -7,7 +7,7 @@ 0.1.0.$Revision: 19276 $(alpha) http://wiki.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modules/PkModeling - Yingxuan Zhu (while at GE), Andrey Fedorov (SPL), John Evans (MGH), Jim Miller (GE) + Yingxuan Zhu (while at GE), Andrey Fedorov (SPL), John Evans (MGH), Jim Miller (GE), Andrew Beers (MGH) @@ -212,6 +212,15 @@ input 1 + + ToftsIntegrationMethod + ToftsIntegrationMethod + + Determine which integration method will be used when calculating Mean Squared Error in the fitting process. + Recursive + Recursive + Convolutional + OutputRSquaredFileName outputRSquared diff --git a/CLI/itkConcentrationToQuantitativeImageFilter.h b/CLI/itkConcentrationToQuantitativeImageFilter.h index 7bce656..d94114f 100644 --- a/CLI/itkConcentrationToQuantitativeImageFilter.h +++ b/CLI/itkConcentrationToQuantitativeImageFilter.h @@ -135,6 +135,8 @@ namespace itk itkSetMacro(constantBAT, int); itkGetMacro(BATCalculationMode, std::string); itkSetMacro(BATCalculationMode, std::string); + itkGetMacro(ToftsIntegrationMethod, std::string); + itkSetMacro(ToftsIntegrationMethod, std::string); void SetTiming(const std::vector& inputTiming); const std::vector& GetTiming(); @@ -246,6 +248,7 @@ namespace itk int m_ModelType; bool m_MaskByRSquared; int m_constantBAT; + std::string m_ToftsIntegrationMethod; std::string m_BATCalculationMode; std::vector m_Timing; diff --git a/CLI/itkConcentrationToQuantitativeImageFilter.hxx b/CLI/itkConcentrationToQuantitativeImageFilter.hxx index 4fc10af..ceb0964 100644 --- a/CLI/itkConcentrationToQuantitativeImageFilter.hxx +++ b/CLI/itkConcentrationToQuantitativeImageFilter.hxx @@ -37,6 +37,7 @@ namespace itk m_UsePrescribedAIF = false; m_MaskByRSquared = true; m_ModelType = itk::LMCostFunction::TOFTS_2_PARAMETER; + m_ToftsIntegrationMethod = "Recursive"; m_constantBAT = 0; m_BATCalculationMode = "PeakGradient"; this->Superclass::SetNumberOfRequiredInputs(1); @@ -438,6 +439,7 @@ namespace itk optimizerErrorCode = pk_solver(timeSize, &timeMinute[0], const_cast(shiftedVectorVoxel.GetDataPointer()), &m_AIF[0], + m_ToftsIntegrationMethod, tempKtrans, tempVe, tempFpv, m_fTol, m_gTol, m_xTol, m_epsilon, m_maxIter, m_hematocrit, @@ -862,6 +864,7 @@ namespace itk ::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os, indent); + os << indent << "Integration Method: " << m_ToftsIntegrationMethod << std::endl; os << indent << "Function tolerance: " << m_fTol << std::endl; os << indent << "Gradient tolerance: " << m_gTol << std::endl; os << indent << "Parameter tolerance: " << m_xTol << std::endl; diff --git a/CLI/itkSignalIntensityToConcentrationImageFilter.h b/CLI/itkSignalIntensityToConcentrationImageFilter.h index 872f15b..262568c 100644 --- a/CLI/itkSignalIntensityToConcentrationImageFilter.h +++ b/CLI/itkSignalIntensityToConcentrationImageFilter.h @@ -105,6 +105,8 @@ namespace itk itkSetMacro(S0GradThresh, float); itkGetMacro(BATCalculationMode, std::string); itkSetMacro(BATCalculationMode, std::string); + itkGetMacro(ToftsIntegrationMethod, std::string); + itkSetMacro(ToftsIntegrationMethod, std::string); itkGetMacro(constantBAT, int); itkSetMacro(constantBAT, int); @@ -173,6 +175,7 @@ namespace itk float m_RGD_relaxivity; float m_S0GradThresh; std::string m_BATCalculationMode; + std::string m_ToftsIntegrationMethod; int m_constantBAT; }; diff --git a/PkSolver/PkSolver.cxx b/PkSolver/PkSolver.cxx index 075ee96..49a4dff 100644 --- a/PkSolver/PkSolver.cxx +++ b/PkSolver/PkSolver.cxx @@ -29,6 +29,7 @@ namespace itk int m_ConstantBAT; std::string m_BATCalculationMode; + std::string m_ToftsIntegrationMethod; // // Implementation of the PkSolver API @@ -38,6 +39,7 @@ namespace itk bool pk_solver(int signalSize, const float* timeAxis, const float* PixelConcentrationCurve, const float* BloodConcentrationCurve, + const std::string ToftsIntegrationMethod, float& Ktrans, float& Ve, float& Fpv, float fTol, float gTol, float xTol, float epsilon, int maxIter, @@ -55,6 +57,7 @@ namespace itk m_BATCalculationMode = BATCalculationMode; m_ConstantBAT = constantBAT; + m_ToftsIntegrationMethod = ToftsIntegrationMethod; // Levenberg Marquardt optimizer itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New(); @@ -79,6 +82,7 @@ namespace itk costFunction->SetCv(PixelConcentrationCurve, signalSize); //Signal Y costFunction->SetTime(timeAxis, signalSize); //Signal X costFunction->SetHematocrit(hematocrit); + costFunction->SetIntegrationType(m_ToftsIntegrationMethod); costFunction->GetValue(initialValue); costFunction->SetModelType(modelType); @@ -147,6 +151,7 @@ namespace itk unsigned pk_solver(int signalSize, const float* timeAxis, const float* PixelConcentrationCurve, const float* BloodConcentrationCurve, + const std::string ToftsIntegrationMethod, float& Ktrans, float& Ve, float& Fpv, float fTol, float gTol, float xTol, float epsilon, int maxIter, @@ -170,6 +175,7 @@ namespace itk m_BATCalculationMode = BATCalculationMode; m_ConstantBAT = constantBAT; + m_ToftsIntegrationMethod = ToftsIntegrationMethod; // Levenberg Marquardt optimizer @@ -194,6 +200,7 @@ namespace itk costFunction->SetCv(PixelConcentrationCurve, signalSize); //Signal Y costFunction->SetTime(timeAxis, signalSize); //Signal X costFunction->SetHematocrit(hematocrit); + costFunction->SetIntegrationType(m_ToftsIntegrationMethod); costFunction->GetValue(initialValue); //... costFunction->SetModelType(modelType); diff --git a/PkSolver/PkSolver.h b/PkSolver/PkSolver.h index 8277326..3c097af 100644 --- a/PkSolver/PkSolver.h +++ b/PkSolver/PkSolver.h @@ -87,7 +87,7 @@ namespace itk float m_Hematocrit; - + std::string m_IntegrationType; int m_ModelType; LMCostFunction() @@ -99,6 +99,11 @@ namespace itk m_Hematocrit = hematocrit; } + void SetIntegrationType(std::string ToftsIntegrationMethod) + { + m_IntegrationType = ToftsIntegrationMethod; + } + void SetModelType(int model) { m_ModelType = model; @@ -137,25 +142,50 @@ namespace itk MeasureType GetValue(const ParametersType & parameters) const { MeasureType measure(RangeDimension); + MeasureType cost_metric(RangeDimension); ValueType Ktrans = parameters[0]; ValueType Ve = parameters[1]; - ArrayType VeTerm; + Array VeTerm; VeTerm = -Ktrans / Ve*Time; ValueType deltaT = Time(1) - Time(0); + // Precalculated function chunks to make the recursive integration method more intelligble. + ValueType log_e = (-Ktrans / Ve)*deltaT; + ValueType capital_E = exp(log_e); + ValueType log_e_2 = pow(log_e, 2); + ValueType block_A = capital_E - log_e - 1; + ValueType block_B = capital_E - (capital_E * log_e) - 1; + ValueType block_ktrans = Ktrans * deltaT / log_e_2; + if (m_ModelType == TOFTS_3_PARAMETER) { ValueType f_pv = parameters[2]; measure = Cv - (1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm)) + f_pv*Cb)); + cost_metric = measure; } else if (m_ModelType == TOFTS_2_PARAMETER) { - measure = Cv - (1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm)))); + if (m_IntegrationType == "Convolutional") + { + measure = Cv - (1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm)))); + cost_metric = measure; + } + + else if (m_IntegrationType == "Recursive") + { + measure[0] = 0; + cost_metric[0] = 0; + for (unsigned int t = 1; t < RangeDimension; ++t) + { + measure[t] = measure[t - 1] * capital_E + (1 / (1.0 - m_Hematocrit)) * block_ktrans * (Cb[t] * block_A - Cb[t - 1] * block_B); + cost_metric[t] = Cv[t] - measure[t]; + } + } } - return measure; + return cost_metric; } MeasureType GetFittedFunction(const ParametersType & parameters) const @@ -169,6 +199,14 @@ namespace itk VeTerm = -Ktrans / Ve*Time; ValueType deltaT = Time(1) - Time(0); + // Precalculated function chunks to make the recursive integration method more intelligble. + ValueType log_e = (-Ktrans / Ve)*deltaT; + ValueType capital_E = exp(log_e); + ValueType log_e_2 = pow(log_e, 2); + ValueType block_A = capital_E - log_e - 1; + ValueType block_B = capital_E - (capital_E * log_e) - 1; + ValueType block_ktrans = Ktrans * deltaT / log_e_2; + if (m_ModelType == TOFTS_3_PARAMETER) { ValueType f_pv = parameters[2]; @@ -176,7 +214,19 @@ namespace itk } else if (m_ModelType == TOFTS_2_PARAMETER) { - measure = 1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm))); + if (m_IntegrationType == "Convolutional") + { + measure = 1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm))); + } + + else if (m_IntegrationType == "Recursive") + { + measure[0] = 0; + for (unsigned int t = 1; t < RangeDimension; ++t) + { + measure[t] = measure[t - 1] * capital_E + (1 / (1.0 - m_Hematocrit)) * block_ktrans * (Cb[t] * block_A - Cb[t - 1] * block_B); + } + } } return measure; @@ -288,6 +338,7 @@ namespace itk bool pk_solver(int signalSize, const float* timeAxis, const float* PixelConcentrationCurve, const float* BloodConcentrationCurve, + const std::string ToftsIntegrationMethod, float& Ktrans, float& Ve, float& Fpv, float fTol = 1e-4f, float gTol = 1e-4f, @@ -303,10 +354,16 @@ namespace itk // as defined by OptimizerDiagnosticCodes, and masked to indicate // wheather Ktrans or Ve were clamped. unsigned pk_solver(int signalSize, const float* timeAxis, - const float* PixelConcentrationCurve, const float* BloodConcentrationCurve, + const float* PixelConcentrationCurve, + const float* BloodConcentrationCurve, + const std::string ToftsIntegrationMethod, float& Ktrans, float& Ve, float& Fpv, - float fTol, float gTol, float xTol, - float epsilon, int maxIter, float hematocrit, + float fTol, + float gTol, + float xTol, + float epsilon, + int maxIter, + float hematocrit, itk::LevenbergMarquardtOptimizer* optimizer, LMCostFunction* costFunction, int modelType = itk::LMCostFunction::TOFTS_2_PARAMETER,