Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH: Add ability to choose integration method [Issue43] #47

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

ENH: Add ability to choose integration method [Issue43] #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CLI/PkModeling.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 != "")
{
Expand Down Expand Up @@ -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);
Expand Down
11 changes: 10 additions & 1 deletion CLI/PkModeling.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
<version>0.1.0.$Revision: 19276 $(alpha)</version>
<documentation-url>http://wiki.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modules/PkModeling</documentation-url>
<license/>
<contributor>Yingxuan Zhu (while at GE), Andrey Fedorov (SPL), John Evans (MGH), Jim Miller (GE)</contributor>
<contributor>Yingxuan Zhu (while at GE), Andrey Fedorov (SPL), John Evans (MGH), Jim Miller (GE), Andrew Beers (MGH)</contributor>
<acknowledgements><![CDATA[This work is part of the National Alliance for Medical Image Computing (NA-MIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, and by National Cancer Institute as part of the Quantitative Imaging Network initiative (U01CA151261) and QIICR (U24CA180918).]]></acknowledgements>
<parameters>
<label>PkModeling Parameters</label>
Expand Down Expand Up @@ -212,6 +212,15 @@
<channel>input</channel>
<default>1</default>
</integer>
<string-enumeration>
<name>ToftsIntegrationMethod</name>
<longflag>ToftsIntegrationMethod</longflag>
<label>Fitting Integration Method</label>
<description>Determine which integration method will be used when calculating Mean Squared Error in the fitting process.</description>
<default>Recursive</default>
<element>Recursive</element>
<element>Convolutional</element>
</string-enumeration>
<image>
<name>OutputRSquaredFileName</name>
<longflag>outputRSquared</longflag>
Expand Down
3 changes: 3 additions & 0 deletions CLI/itkConcentrationToQuantitativeImageFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<float>& inputTiming);
const std::vector<float>& GetTiming();
Expand Down Expand Up @@ -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<float> m_Timing;
Expand Down
3 changes: 3 additions & 0 deletions CLI/itkConcentrationToQuantitativeImageFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -438,6 +439,7 @@ namespace itk
optimizerErrorCode = pk_solver(timeSize, &timeMinute[0],
const_cast<float *>(shiftedVectorVoxel.GetDataPointer()),
&m_AIF[0],
m_ToftsIntegrationMethod,
tempKtrans, tempVe, tempFpv,
m_fTol, m_gTol, m_xTol,
m_epsilon, m_maxIter, m_hematocrit,
Expand Down Expand Up @@ -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;
Expand Down
3 changes: 3 additions & 0 deletions CLI/itkSignalIntensityToConcentrationImageFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -173,6 +175,7 @@ namespace itk
float m_RGD_relaxivity;
float m_S0GradThresh;
std::string m_BATCalculationMode;
std::string m_ToftsIntegrationMethod;
int m_constantBAT;
};

Expand Down
7 changes: 7 additions & 0 deletions PkSolver/PkSolver.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ namespace itk

int m_ConstantBAT;
std::string m_BATCalculationMode;
std::string m_ToftsIntegrationMethod;

//
// Implementation of the PkSolver API
Expand All @@ -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,
Expand All @@ -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();
Expand All @@ -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);

Expand Down Expand Up @@ -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,
Expand All @@ -170,6 +175,7 @@ namespace itk

m_BATCalculationMode = BATCalculationMode;
m_ConstantBAT = constantBAT;
m_ToftsIntegrationMethod = ToftsIntegrationMethod;

// Levenberg Marquardt optimizer

Expand All @@ -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);

Expand Down
73 changes: 65 additions & 8 deletions PkSolver/PkSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ namespace itk


float m_Hematocrit;

std::string m_IntegrationType;
int m_ModelType;

LMCostFunction()
Expand All @@ -99,6 +99,11 @@ namespace itk
m_Hematocrit = hematocrit;
}

void SetIntegrationType(std::string ToftsIntegrationMethod)
{
m_IntegrationType = ToftsIntegrationMethod;
}

void SetModelType(int model)
{
m_ModelType = model;
Expand Down Expand Up @@ -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 <double> 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
Expand All @@ -169,14 +199,34 @@ 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];
measure = 1 / (1.0 - m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb, Exponential(VeTerm)) + f_pv*Cb);
}
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;
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down