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

Change the integration method #39

Closed
wants to merge 18 commits into from
Closed

Change the integration method #39

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
f010ecd
Changed integration method to create good performance on DRO
AnnaBeers Jan 13, 2017
d284f5b
Removed debugging comments
AnnaBeers Jan 13, 2017
79816a1
Indentation Edits
AnnaBeers Jan 17, 2017
e4b8ab4
Second Indentation Edits..
AnnaBeers Jan 18, 2017
3ee8cdc
(temporary?) Solution for NaN Streak from warm fitting
AnnaBeers Jan 19, 2017
cfd0488
Slight fixes to possibly avoid NaN values
AnnaBeers Feb 2, 2017
6466688
New Integration Method - Simplex
AnnaBeers Feb 20, 2017
fbf24df
Added a radio button under "Advanced Options" to select a preferred
AnnaBeers Mar 10, 2017
bbd1d61
Changing the fitting method means switching between two classes. To
AnnaBeers Mar 22, 2017
08fae37
Made much more sensible changes to allow switching fitting algorithm,
AnnaBeers Mar 23, 2017
4cc4203
Potentially found a solution to the multiple fitting algorithm problem,
AnnaBeers Mar 27, 2017
59753e1
These options finally work. Now it's just time to clean up the code.
AnnaBeers Mar 28, 2017
c714194
Some updates to consolidate code. More updates to come.
AnnaBeers Apr 5, 2017
4d8e0c6
Cleaned the code up. It's more or less ready for production. The clock
AnnaBeers Apr 12, 2017
5edf8c0
Changed the name of the Algorithm button to be shorter, take up less
AnnaBeers Apr 14, 2017
8ea4721
Merge pull request #1 from QTIM-Lab/ChooseFittingMethod
AnnaBeers Apr 14, 2017
e585ab0
Reverted accidental parameter changes to Prostate tests.
AnnaBeers Apr 20, 2017
1cc1199
Merge pull request #2 from QTIM-Lab/ChooseFittingMethod
AnnaBeers Apr 20, 2017
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
22 changes: 22 additions & 0 deletions Data/Tofts_DROs/dynamic - as a 1321 frames MultiVolume by TriggerTime.nhdr
View file
Open in desktop

Large diffs are not rendered by default.

Binary file added Data/Tofts_DROs/dynamic - as a 1321 frames MultiVolume by TriggerTime.raw.gz
View file
Open in desktop
Binary file not shown.
Binary file added Data/Tofts_DROs/tofts_v6_label_AIF.nrrd
View file
Open in desktop
Binary file not shown.
Binary file added Data/Tofts_DROs/tofts_v6_label_ROI.nrrd
View file
Open in desktop
Binary file not shown.
7 changes: 6 additions & 1 deletion PkSolver/PkSolver.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,11 @@
#include <itkGradientMagnitudeImageFilter.h>
#include <itkImageRegionIterator.h>
#include <itkLevenbergMarquardtOptimizer.h>
// #include <itkAmoebaOptimizer.h>
#include "PkSolver.h"
#include "itkTimeProbesCollectorBase.h"
#include <string>

// hello
namespace itk
{
//
Expand Down Expand Up @@ -58,6 +59,10 @@ bool pk_solver (int signalSize, const float* timeAxis,

// Levenberg Marquardt optimizer
itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New();

// Simplex optimizer
// itk::AmoebaOptimizer::Pointer optimizer = itk::AmoebaOptimizer::New();

LMCostFunction::Pointer costFunction = LMCostFunction::New();

LMCostFunction::ParametersType initialValue;
Expand Down
120 changes: 91 additions & 29 deletions PkSolver/PkSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#define PkSolver_h_

#include "itkLevenbergMarquardtOptimizer.h"
// #include "itkAmoebaOptimizer.h"
#include <math.h>
#include <vnl/algo/vnl_convolve.h>
#include "itkArray.h"
Expand Down Expand Up @@ -83,7 +84,6 @@ class LMCostFunction: public itk::MultipleValuedCostFunction
typedef Superclass::MeasureType MeasureType, ArrayType;
typedef Superclass::ParametersValueType ValueType;


float m_Hematocrit;

int m_ModelType;
Expand Down Expand Up @@ -133,37 +133,76 @@ class LMCostFunction: public itk::MultipleValuedCostFunction
MeasureType GetValue( const ParametersType & parameters) const
{
MeasureType measure(RangeDimension);
MeasureType measure2(RangeDimension);
MeasureType measure3(RangeDimension);

ValueType Ktrans = parameters[0];
ValueType Ve = parameters[1];
//Ktrans = .2;
//Ve = .1;

ArrayType VeTerm;
VeTerm = -Ktrans/Ve*Time;
ValueType deltaT = Time(1) - Time(0);
ArrayType VeTerm;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AndrewBeers can you fix the indentation to be consistent with the rest of the code?

VeTerm = -Ktrans / Ve*Time;
ValueType deltaT = Time(1) - Time(0);

if( m_ModelType == TOFTS_3_PARAMETER)
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));
}

}
else if(m_ModelType == TOFTS_2_PARAMETER)
{
measure = Cv - (1/(1.0-m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb,Exponential(VeTerm))));
}

// This is the original integration method.
// measure = Cv - (1/(1.0-m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb,Exponential(VeTerm))));

return measure;
// This is the new method.
measure2[0] = Cv[0];
for (unsigned int t = 1; t < RangeDimension; ++t) {
measure2[t] = measure2[t - 1] * capital_E + (1 / (1.0 - m_Hematocrit)) * block_ktrans * (Cb[t] * block_A - Cb[t - 1] * block_B);
}

// I fully understand that this is a strange way to subtract from observed.
for (unsigned int t = 1; t < RangeDimension; ++t) {
measure3[t] = Cv[t] - measure2[t];
}

}

measure = measure3;

return measure;
}

MeasureType GetFittedFunction( const ParametersType & parameters) const
MeasureType GetFittedFunction(const ParametersType & parameters) const
{
MeasureType measure(RangeDimension);
MeasureType measure(RangeDimension);
MeasureType measure2(RangeDimension);

ValueType Ktrans = parameters[0];
ValueType Ve = parameters[1];
ValueType Ktrans = parameters[0];
ValueType Ve = parameters[1];

ArrayType VeTerm;
VeTerm = -Ktrans / Ve*Time;
ValueType deltaT = Time(1) - Time(0);

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;

ArrayType VeTerm;
VeTerm = -Ktrans/Ve*Time;
ValueType deltaT = Time(1) - Time(0);

if( m_ModelType == TOFTS_3_PARAMETER)
{
Expand All @@ -172,9 +211,20 @@ class LMCostFunction: public itk::MultipleValuedCostFunction
}
else if(m_ModelType == TOFTS_2_PARAMETER)
{
measure = 1/(1.0-m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb,Exponential(VeTerm)));

// original method
// measure = 1/(1.0-m_Hematocrit)*(Ktrans*deltaT*Convolution(Cb,Exponential(VeTerm)));

// new method
measure2[0] = 0;
for (unsigned int t = 1; t < RangeDimension; ++t) {
measure2[t] = (measure2[t - 1] * capital_E + (1 / (1.0 - m_Hematocrit)) * block_ktrans * (Cb[t] * block_A - Cb[t - 1] * block_B));
}

}

measure = measure2;

return measure;
}

Expand Down Expand Up @@ -261,6 +311,7 @@ class CommandIterationUpdateLevenbergMarquardt : public itk::Command

void Execute(const itk::Object * object, const itk::EventObject & event)
{
//std::cout << m_IterationNumber++ << std::endl;
//std::cout << "Observer::Execute() " << std::endl;
OptimizerPointer optimizer =
dynamic_cast< OptimizerPointer >( object );
Expand All @@ -282,10 +333,13 @@ class CommandIterationUpdateLevenbergMarquardt : public itk::Command
itk::FunctionEvaluationIterationEvent m_FunctionEvent;
itk::GradientEvaluationIterationEvent m_GradientEvent;
};
bool pk_solver(int signalSize, const float* timeAxis,
bool pk_solver(int signalSize,
const float* timeAxis,
const float* PixelConcentrationCurve,
const float* BloodConcentrationCurve,
float& Ktrans, float& Ve, float& Fpv,
float& Ktrans,
float& Ve,
float& Fpv,
float fTol = 1e-4f,
float gTol = 1e-4f,
float xTol = 1e-5f,
Expand All @@ -299,16 +353,24 @@ bool pk_solver(int signalSize, const float* timeAxis,
// returns diagnostic error code from the VNL optimizer,
// 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,
float& Ktrans, float& Ve, float& Fpv,
float fTol, float gTol,float xTol,
float epsilon, int maxIter, float hematocrit,
itk::LevenbergMarquardtOptimizer* optimizer,
LMCostFunction* costFunction,
int modelType = itk::LMCostFunction::TOFTS_2_PARAMETER,
int constantBAT = 0,
const std::string BATCalculationMode = "PeakGradient");
unsigned pk_solver(int signalSize,
const float* timeAxis,
const float* PixelConcentrationCurve,
const float* BloodConcentrationCurve,
float& Ktrans,
float& Ve,
float& Fpv,
float fTol,
float gTol,
float xTol,
float epsilon,
int maxIter,
float hematocrit,
itk::LevenbergMarquardtOptimizer* optimizer,
LMCostFunction* costFunction,
int modelType = itk::LMCostFunction::TOFTS_2_PARAMETER,
int constantBAT = 0,
const std::string BATCalculationMode = "PeakGradient");

void pk_report();
void pk_clear();
Expand Down