Skip to content

ENH: Refactor MLP with uncertainty to pass 2D uncertainty matrices #31

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

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

ENH: Refactor MLP with uncertainty to pass 2D uncertainty matrices #31

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
31 changes: 27 additions & 4 deletions applications/pctbinning/pctbinning.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,15 +44,38 @@ main(int argc, char * argv[])
projection->SetSourceDistance(args_info.source_arg);
projection->SetMostLikelyPathType(args_info.mlptype_arg);
projection->SetMostLikelyPathPolynomialDegree(args_info.mlppolydeg_arg);
projection->SetMostLikelyPathTrackerUncertainties(args_info.mlptrackeruncert_flag);
projection->SetTrackerResolution(args_info.trackerresolution_arg);
projection->SetTrackerPairSpacing(args_info.trackerspacing_arg);
projection->SetMaterialBudget(args_info.materialbudget_arg);
projection->SetIonizationPotential(args_info.ionpot_arg * CLHEP::eV);
projection->SetRobust(args_info.robust_flag);
projection->SetComputeScattering(args_info.scatwepl_given);
projection->SetComputeNoise(args_info.noise_given);

// Uncertainties
ProjectionFilter::TwoDMatrixType sigmaIn, sigmaOut;
if (args_info.sigmaIn_given)
{
if (args_info.sigmaIn_given != 4)
{
std::cerr << "--sigmaIn needs exactly 4 values" << std::endl;
exit(EXIT_FAILURE);
}
for (unsigned int i = 0; i < 2; i++)
for (unsigned int j = 0; j < 2; j++)
sigmaIn[i][j] = args_info.sigmaIn_arg[i * 2 + j];
projection->SetSigmaIn(sigmaIn);
}
if (args_info.sigmaOut_given)
{
if (args_info.sigmaOut_given != 4)
{
std::cerr << "--sigmaOut needs exactly 4 values" << std::endl;
exit(EXIT_FAILURE);
}
for (unsigned int i = 0; i < 2; i++)
for (unsigned int j = 0; j < 2; j++)
sigmaOut[i][j] = args_info.sigmaOut_arg[i * 2 + j];
projection->SetSigmaOut(sigmaOut);
}

if (args_info.quadricIn_given)
{
// quadric = object surface
Expand Down
9 changes: 5 additions & 4 deletions applications/pctbinning/pctbinning.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,14 @@ option "source" s "Source position"
option "quadricIn" - "Parameters of the entrance quadric support function, see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" double multiple no
option "quadricOut" - "Parameters of the exit quadric support function, see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" double multiple no
option "mlptype" - "Type of most likely path (schulte, polynomial, or krah)" string no default="schulte"
option "mlptrackeruncert" - "Consider tracker uncertainties in MLP [Krah 2018, PMB]" flag off
option "mlppolydeg" - "Degree of the polynomial to approximate 1/beta^2p^2" int no default="5"
option "ionpot" - "Ionization potential used in the reconstruction in eV" double no default="68.9984"
option "fill" - "Fill holes, i.e. pixels that were not hit by protons" flag off
option "trackerresolution" - "Tracker resolution in mm" double no
option "trackerspacing" - "Tracker pair spacing in mm" double no
option "materialbudget" - "Material budget x/X0 of tracker" double no

section "MLP uncertainty [Krah et al, 2018, PMB]"
option "sigmaIn" - "Entrance uncertainty matrix (equation 16)" double multiple no
option "sigmaOut" - "Exit uncertainty matrix (equation 17)" double multiple no


section "Projections parameters"
option "origin" - "Origin (default=centered)" double multiple no
Expand Down
4 changes: 3 additions & 1 deletion include/pctEnergyAdaptiveMLPFunction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,8 @@ class PCT_EXPORT EnergyAdaptiveMLPFunction : public MostLikelyPathFunction<doubl
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;

using TwoDMatrixType = Superclass::TwoDMatrixType;

/** Method for creation through the object factory. */
itkNewMacro(Self);

Expand All @@ -130,7 +132,7 @@ class PCT_EXPORT EnergyAdaptiveMLPFunction : public MostLikelyPathFunction<doubl

/** Evaluate the error (x,y) (equation 27) at depth z. */
void
EvaluateError(const double u1, itk::Matrix<double, 2, 2> & error);
EvaluateError(const double u1, TwoDMatrixType & error);

#ifdef MLP_TIMING
/** Print timing information */
Expand Down
11 changes: 2 additions & 9 deletions include/pctMostLikelyPathFunction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ class ITK_TEMPLATE_EXPORT MostLikelyPathFunction : public itk::LightObject

/** Useful defines. */
using VectorType = itk::Vector<TCoordRep, 3>;
using TwoDMatrixType = typename itk::Matrix<TCoordRep, 2, 2>;

/** Init the mlp parameters from the input and output directions and positions. */
virtual void
Expand All @@ -52,15 +53,7 @@ class ITK_TEMPLATE_EXPORT MostLikelyPathFunction : public itk::LightObject

/** Init with additional parameters to consider tracker uncertainties */
virtual void
InitUncertain(const VectorType posIn,
const VectorType posOut,
const VectorType dirIn,
const VectorType dirOut,
double dEntry,
double dExit,
double m_TrackerResolution,
double m_TrackerPairSpacing,
double m_MaterialBudget)
InitUncertain(double dEntry, double dExit, TwoDMatrixType sigmaIn, TwoDMatrixType sigmaOut)
{
itkGenericExceptionMacro("Not implemented in the derived class.");
}
Expand Down
33 changes: 16 additions & 17 deletions include/pctProtonPairsToDistanceDrivenProjection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,15 @@
#ifndef __pctProtonPairsToDistanceDrivenProjection_h
#define __pctProtonPairsToDistanceDrivenProjection_h

#include "rtkConfiguration.h"
#include "pctBetheBlochFunctor.h"
#include <mutex>

#include <rtkQuadricShape.h>
#include <itkInPlaceImageFilter.h>
#include <mutex>

#include <rtkConfiguration.h>
#include <rtkQuadricShape.h>

#include "pctBetheBlochFunctor.h"
#include "pctMostLikelyPathFunction.h"

namespace pct
{
Expand Down Expand Up @@ -36,6 +39,7 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection : public itk::In
using OutputImageRegionType = typename OutputImageType::RegionType;

using RQIType = rtk::QuadricShape;
using TwoDMatrixType = typename MostLikelyPathFunction<double>::TwoDMatrixType;

/** Method for creation through the object factory. */
itkNewMacro(Self);
Expand All @@ -55,19 +59,16 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection : public itk::In
itkGetMacro(MostLikelyPathType, std::string);
itkSetMacro(MostLikelyPathType, std::string);

itkGetMacro(MostLikelyPathTrackerUncertainties, bool);
itkSetMacro(MostLikelyPathTrackerUncertainties, bool);
/** Tracker uncertainty parameters ([Krah et al, 2018, PMB] section 2.5) */
itkGetMacro(SigmaIn, TwoDMatrixType);
itkSetMacro(SigmaIn, TwoDMatrixType);

itkGetMacro(SigmaOut, TwoDMatrixType);
itkSetMacro(SigmaOut, TwoDMatrixType);

itkGetMacro(MostLikelyPathPolynomialDegree, int);
itkSetMacro(MostLikelyPathPolynomialDegree, int);

itkGetMacro(TrackerResolution, double);
itkSetMacro(TrackerResolution, double);
itkGetMacro(TrackerPairSpacing, double);
itkSetMacro(TrackerPairSpacing, double);
itkGetMacro(MaterialBudget, double);
itkSetMacro(MaterialBudget, double);

/** Get/Set the boundaries of the object. */
itkGetMacro(QuadricIn, RQIType::Pointer);
itkSetMacro(QuadricIn, RQIType::Pointer);
Expand Down Expand Up @@ -140,10 +141,8 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection : public itk::In
bool m_VariableBeamEnergy = false;

// MLP considering tracker uncertainties
bool m_MostLikelyPathTrackerUncertainties;
double m_TrackerResolution;
double m_TrackerPairSpacing;
double m_MaterialBudget;
TwoDMatrixType m_SigmaIn;
TwoDMatrixType m_SigmaOut;

/** Count event in each thread */
CountImagePointer m_Count;
Expand Down
23 changes: 11 additions & 12 deletions include/pctProtonPairsToDistanceDrivenProjection.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,10 @@ ProtonPairsToDistanceDrivenProjection<TInputImage, TOutputImage>::ThreadedGenera
itkGenericExceptionMacro("MLP must either be schulte, polynomial, krah, or adaptive, not [" << m_MostLikelyPathType
<< ']');
}
if (m_MostLikelyPathTrackerUncertainties && m_MostLikelyPathType != "schulte")

const TwoDMatrixType id{ TwoDMatrixType::GetIdentity() };
bool bMLPTrackerUncertainties = !(m_SigmaIn==id && m_SigmaOut==id);
if (bMLPTrackerUncertainties && m_MostLikelyPathType != "schulte")
{
itkGenericExceptionMacro("Tracker uncertainties can currently only be considered with MLP type 'Schulte'.");
}
Expand Down Expand Up @@ -298,17 +301,13 @@ ProtonPairsToDistanceDrivenProjection<TInputImage, TOutputImage>::ThreadedGenera
// Init MLP before mm to voxel conversion
double xIn, xOut, yIn, yOut;
double dxIn, dxOut, dyIn, dyOut;
if (m_MostLikelyPathTrackerUncertainties)
if (bMLPTrackerUncertainties)
{
mlp->InitUncertain(pSIn,
pSOut,
dIn,
dOut,
distanceEntry,
mlp->InitUncertain(distanceEntry,
distanceExit,
m_TrackerResolution,
m_TrackerPairSpacing,
m_MaterialBudget);
m_SigmaIn,
m_SigmaOut);
mlp->Init(pSIn, pSOut, dIn, dOut);
mlp->Evaluate(pSIn[2], xIn, yIn, dxIn, dyIn); // get entrance and exit position according to MLP
mlp->Evaluate(pSOut[2], xOut, yOut, dxOut, dyOut);
}
Expand Down Expand Up @@ -337,7 +336,7 @@ ProtonPairsToDistanceDrivenProjection<TInputImage, TOutputImage>::ThreadedGenera
double dxDummy, dyDummy;

double dInMLP[2];
if (m_MostLikelyPathTrackerUncertainties && QuadricIntersected)
if (bMLPTrackerUncertainties && QuadricIntersected)
{
dInMLP[0] = dxIn;
dInMLP[1] = dyIn;
Expand All @@ -349,7 +348,7 @@ ProtonPairsToDistanceDrivenProjection<TInputImage, TOutputImage>::ThreadedGenera
}

double dOutMLP[2];
if (m_MostLikelyPathTrackerUncertainties && QuadricIntersected)
if (bMLPTrackerUncertainties && QuadricIntersected)
{
dOutMLP[0] = dxOut;
dOutMLP[1] = dyOut;
Expand Down
58 changes: 26 additions & 32 deletions include/pctSchulteMLPFunction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction<double>
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;

using TwoDMatrixType = Superclass::TwoDMatrixType;

/** Method for creation through the object factory. */
itkNewMacro(Self);

Expand All @@ -158,23 +160,15 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction<double>

/** Init with additional parameters to consider tracker uncertainties */
virtual void
InitUncertain(const VectorType posIn,
const VectorType posOut,
const VectorType dirIn,
const VectorType dirOut,
double dEntry,
double dExit,
double TrackerResolution,
double TrackerPairSpacing,
double MaterialBudget) override;
InitUncertain(double dEntry, double dExit, TwoDMatrixType sigmaIn, TwoDMatrixType sigmaOut) override;

/** Evaluate the coordinates (x,y) at depth u1. */
virtual void
Evaluate(const double u1, double & x, double & y, double & dx, double & dy) override;

/** Evaluate the error (x,y) (equation 27) at depth z. */
void
EvaluateError(const double u1, itk::Matrix<double, 2, 2> & error);
EvaluateError(const double u1, TwoDMatrixType & error);

#ifdef MLP_TIMING
/** Print timing information */
Expand All @@ -185,10 +179,10 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction<double>
protected:
/// Implementation of 2x2 matrix inversion, faster than itk/vnl inversion
void
InverseMatrix(itk::Matrix<double, 2, 2> & mat);
InverseMatrix(TwoDMatrixType & mat);

/// Constructor
SchulteMLPFunction();
SchulteMLPFunction() = default;

/// Destructor
~SchulteMLPFunction() {}
Expand All @@ -200,8 +194,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction<double>

// Depth position at entrance and exit, only u1 is variable
double m_uOrigin;
double m_u0;
double m_u2;
double m_u0{0.};
double m_u2{0.};

// Entrance and exit parameters (equation 1)
itk::Vector<double, 2> m_x0;
Expand All @@ -210,26 +204,26 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction<double>
itk::Vector<double, 2> m_y2;

// Part of the rotation matrices which is constant for the trajectory
itk::Matrix<double, 2, 2> m_R0;
itk::Matrix<double, 2, 2> m_R1;
itk::Matrix<double, 2, 2> m_R0T;
itk::Matrix<double, 2, 2> m_R1T;
TwoDMatrixType m_R0{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_R1{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_R0T{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_R1T{TwoDMatrixType::GetIdentity()};

// Scattering matrices
itk::Matrix<double, 2, 2> m_Sigma1;
itk::Matrix<double, 2, 2> m_Sigma2;

bool m_considerTrackerUncertainties;

// Addtional matrices needed for MLP with tracker uncertainties
itk::Matrix<double, 2, 2> m_SigmaIn;
itk::Matrix<double, 2, 2> m_SigmaOut;
itk::Matrix<double, 2, 2> m_Sin;
itk::Matrix<double, 2, 2> m_SinT;
itk::Matrix<double, 2, 2> m_Sout;
itk::Matrix<double, 2, 2> m_SoutT;
itk::Matrix<double, 2, 2> m_Sout_Inv;
itk::Matrix<double, 2, 2> m_SoutT_Inv;
TwoDMatrixType m_Sigma1;
TwoDMatrixType m_Sigma2;

// Members for MLP with tracker uncertainties
bool m_ConsiderTrackerUncertainties;
TwoDMatrixType m_SigmaIn;
TwoDMatrixType m_SigmaOut;
TwoDMatrixType m_SIn{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_SInT{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_SOut{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_SOutT{TwoDMatrixType::GetIdentity()};
TwoDMatrixType m_SOut_Inv;
TwoDMatrixType m_SOutT_Inv;

// Part common to all positions along the trajectory
// double m_IntForSigmaSqTheta0; //Always 0. because m_u0=0.
// double m_IntForSigmaSqTTheta0; //Always 0. because m_u0=0.
Expand Down
Loading
Loading