diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index f37b969..dea21c2 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/applications/pctbinning/pctbinning.cxx @@ -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 diff --git a/applications/pctbinning/pctbinning.ggo b/applications/pctbinning/pctbinning.ggo index 86a4a06..3f1ef4c 100644 --- a/applications/pctbinning/pctbinning.ggo +++ b/applications/pctbinning/pctbinning.ggo @@ -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 diff --git a/include/pctEnergyAdaptiveMLPFunction.h b/include/pctEnergyAdaptiveMLPFunction.h index 262158f..2b6507f 100644 --- a/include/pctEnergyAdaptiveMLPFunction.h +++ b/include/pctEnergyAdaptiveMLPFunction.h @@ -109,6 +109,8 @@ class PCT_EXPORT EnergyAdaptiveMLPFunction : public MostLikelyPathFunction; using ConstPointer = itk::SmartPointer; + using TwoDMatrixType = Superclass::TwoDMatrixType; + /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -130,7 +132,7 @@ class PCT_EXPORT EnergyAdaptiveMLPFunction : public MostLikelyPathFunction & error); + EvaluateError(const double u1, TwoDMatrixType & error); #ifdef MLP_TIMING /** Print timing information */ diff --git a/include/pctMostLikelyPathFunction.h b/include/pctMostLikelyPathFunction.h index 3b35532..01a9d6c 100644 --- a/include/pctMostLikelyPathFunction.h +++ b/include/pctMostLikelyPathFunction.h @@ -30,6 +30,7 @@ class ITK_TEMPLATE_EXPORT MostLikelyPathFunction : public itk::LightObject /** Useful defines. */ using VectorType = itk::Vector; + using TwoDMatrixType = typename itk::Matrix; /** Init the mlp parameters from the input and output directions and positions. */ virtual void @@ -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."); } diff --git a/include/pctProtonPairsToDistanceDrivenProjection.h b/include/pctProtonPairsToDistanceDrivenProjection.h index c4f17d1..109de9c 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.h +++ b/include/pctProtonPairsToDistanceDrivenProjection.h @@ -1,12 +1,15 @@ #ifndef __pctProtonPairsToDistanceDrivenProjection_h #define __pctProtonPairsToDistanceDrivenProjection_h -#include "rtkConfiguration.h" -#include "pctBetheBlochFunctor.h" +#include -#include #include -#include + +#include +#include + +#include "pctBetheBlochFunctor.h" +#include "pctMostLikelyPathFunction.h" namespace pct { @@ -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::TwoDMatrixType; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -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); @@ -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; diff --git a/include/pctProtonPairsToDistanceDrivenProjection.hxx b/include/pctProtonPairsToDistanceDrivenProjection.hxx index e32b505..9779a4f 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.hxx +++ b/include/pctProtonPairsToDistanceDrivenProjection.hxx @@ -81,7 +81,10 @@ ProtonPairsToDistanceDrivenProjection::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'."); } @@ -298,17 +301,13 @@ ProtonPairsToDistanceDrivenProjection::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); } @@ -337,7 +336,7 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera double dxDummy, dyDummy; double dInMLP[2]; - if (m_MostLikelyPathTrackerUncertainties && QuadricIntersected) + if (bMLPTrackerUncertainties && QuadricIntersected) { dInMLP[0] = dxIn; dInMLP[1] = dyIn; @@ -349,7 +348,7 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera } double dOutMLP[2]; - if (m_MostLikelyPathTrackerUncertainties && QuadricIntersected) + if (bMLPTrackerUncertainties && QuadricIntersected) { dOutMLP[0] = dxOut; dOutMLP[1] = dyOut; diff --git a/include/pctSchulteMLPFunction.h b/include/pctSchulteMLPFunction.h index f46a873..512209c 100644 --- a/include/pctSchulteMLPFunction.h +++ b/include/pctSchulteMLPFunction.h @@ -146,6 +146,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; + using TwoDMatrixType = Superclass::TwoDMatrixType; + /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -158,15 +160,7 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction /** 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 @@ -174,7 +168,7 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction /** Evaluate the error (x,y) (equation 27) at depth z. */ void - EvaluateError(const double u1, itk::Matrix & error); + EvaluateError(const double u1, TwoDMatrixType & error); #ifdef MLP_TIMING /** Print timing information */ @@ -185,10 +179,10 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction protected: /// Implementation of 2x2 matrix inversion, faster than itk/vnl inversion void - InverseMatrix(itk::Matrix & mat); + InverseMatrix(TwoDMatrixType & mat); /// Constructor - SchulteMLPFunction(); + SchulteMLPFunction() = default; /// Destructor ~SchulteMLPFunction() {} @@ -200,8 +194,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction // 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 m_x0; @@ -210,26 +204,26 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction itk::Vector m_y2; // Part of the rotation matrices which is constant for the trajectory - itk::Matrix m_R0; - itk::Matrix m_R1; - itk::Matrix m_R0T; - itk::Matrix 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 m_Sigma1; - itk::Matrix m_Sigma2; - - bool m_considerTrackerUncertainties; - - // Addtional matrices needed for MLP with tracker uncertainties - itk::Matrix m_SigmaIn; - itk::Matrix m_SigmaOut; - itk::Matrix m_Sin; - itk::Matrix m_SinT; - itk::Matrix m_Sout; - itk::Matrix m_SoutT; - itk::Matrix m_Sout_Inv; - itk::Matrix 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. diff --git a/src/pctSchulteMLPFunction.cxx b/src/pctSchulteMLPFunction.cxx index 8d9b0f9..14f12a3 100644 --- a/src/pctSchulteMLPFunction.cxx +++ b/src/pctSchulteMLPFunction.cxx @@ -21,72 +21,25 @@ namespace pct { -SchulteMLPFunction ::SchulteMLPFunction() -{ - // We apply a change of origin, u0 is always 0 - m_u0 = 0.; - - // Construct the constant part of R0 and R1 (equations 11 and 14) - m_R0(0, 0) = 1.; - m_R0(1, 0) = 0.; - m_R0(1, 1) = 1.; - m_R1 = m_R0; - - // Transpose - m_R0T = m_R0.GetTranspose(); - m_R1T = m_R1.GetTranspose(); - - // Construct the constant part of Sin and Sout (Eq. 14 & 15 in Krah 2018, PMB) - // Needed only when tracker uncertainties are considered - m_Sin(0, 0) = 1.; - m_Sin(1, 0) = 0.; - m_Sin(1, 1) = 1.; - m_Sout = m_Sin; -} - // Initialize terms needed to include tracker uncertainties void SchulteMLPFunction -::InitUncertain(const VectorType posIn, const VectorType posOut, const VectorType dirIn, const VectorType dirOut, - double dEntry, double dExit, double TrackerResolution, double TrackerPairSpacing, double MaterialBudget) +::InitUncertain(const double dEntry, const double dExit, const TwoDMatrixType sigmaIn, const TwoDMatrixType sigmaOut) { - m_considerTrackerUncertainties = true; // NK: maybe this should actually go into constructor - m_u2 = posOut[2] - posIn[2]; - const double sigmaPSq = TrackerResolution * TrackerResolution; - // Finish constructing Sin and Sout matrices (Eq. 14 & 15 in Krah 2018, PMB) - m_Sin(0, 1) = dEntry; - m_Sout(0, 1) = dExit; - - // Transpose - m_SinT = m_Sin.GetTranspose(); - m_SoutT = m_Sout.GetTranspose(); - - m_Sout_Inv = m_Sout; - InverseMatrix(m_Sout_Inv); - m_SoutT_Inv = m_SoutT; - InverseMatrix(m_SoutT_Inv); - - m_SigmaIn(0, 0) = 1; - m_SigmaIn(0, 1) = 1 / TrackerPairSpacing; - m_SigmaIn(1, 0) = m_SigmaIn(0, 1); - m_SigmaIn(1, 1) = 2 / TrackerPairSpacing / TrackerPairSpacing; - m_SigmaIn *= sigmaPSq; - - m_SigmaOut(0, 0) = m_SigmaIn(0, 0); - m_SigmaOut(0, 1) = -m_SigmaIn(0, 1); - m_SigmaOut(1, 0) = -m_SigmaIn(1, 0); - m_SigmaOut(1, 1) = m_SigmaIn(1, 1); - m_SigmaOut *= sigmaPSq; - - const double c = 13.6 * CLHEP::MeV * 13.6 * CLHEP::MeV / (36.1 * CLHEP::cm); - const double trackerThickness = MaterialBudget * 36.1 * CLHEP::cm; - m_SigmaIn(1, 1) += Functor::SchulteMLP::IntegralForSigmaSqTheta::GetValue(trackerThickness) * c; - // m_SigmaOut(1,1) += Functor::SchulteMLP::IntegralForSigmaSqTheta::GetValue(trackerThickness) * c; - m_SigmaOut(1, 1) += (Functor::SchulteMLP::IntegralForSigmaSqTheta::GetValue(m_u2 + trackerThickness) - - Functor::SchulteMLP::IntegralForSigmaSqTheta::GetValue(m_u2)) * - c; - - SchulteMLPFunction::Init(posIn, posOut, dirIn, dirOut); + m_ConsiderTrackerUncertainties = true; + + // Finish constructing Sin and Sout matrices + transpose + inverse (Eq. 14 & 15 in Krah 2018, PMB) + m_SIn(0, 1) = m_SInT(1, 0) = dEntry; + m_SOut(0, 1) = m_SOutT(1, 0) = dExit; + + m_SOut_Inv = m_SOut; + InverseMatrix(m_SOut_Inv); + m_SOutT_Inv = m_SOutT; + InverseMatrix(m_SOutT_Inv); + + // Init uncertainty matrices + m_SigmaIn = sigmaIn; + m_SigmaOut = sigmaOut; } // standard part of the Initialization @@ -129,9 +82,9 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d m_R0T(1, 0) = m_R0(0, 1); m_R1T(1, 0) = m_R1(0, 1); - itk::Matrix R1T_Inv(m_R1T); + TwoDMatrixType R1T_Inv(m_R1T); InverseMatrix(R1T_Inv); - itk::Matrix R1_Inv(m_R1); + TwoDMatrixType R1_Inv(m_R1); InverseMatrix(R1_Inv); // Constants used in both integrals @@ -162,15 +115,15 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d itk::Vector xMLP; itk::Vector yMLP; - if (m_considerTrackerUncertainties) + if (m_ConsiderTrackerUncertainties) { - itk::Matrix C1 = m_R0 * m_Sin * m_SigmaIn * m_SinT * m_R0T + m_Sigma1; - itk::Matrix C2 = - R1_Inv * m_Sout_Inv * m_SigmaOut * m_SoutT_Inv * R1T_Inv + R1_Inv * m_Sigma2 * R1T_Inv; - itk::Matrix C1plusC2(C1 + C2); + TwoDMatrixType C1 = m_R0 * m_SIn * m_SigmaIn * m_SInT * m_R0T + m_Sigma1; + TwoDMatrixType C2 = + R1_Inv * m_SOut_Inv * m_SigmaOut * m_SOutT_Inv * R1T_Inv + R1_Inv * m_Sigma2 * R1T_Inv; + TwoDMatrixType C1plusC2(C1 + C2); InverseMatrix(C1plusC2); - itk::Matrix factorIn(C2 * C1plusC2 * m_R0); - itk::Matrix factorOut(C1 * C1plusC2 * R1_Inv); + TwoDMatrixType factorIn(C2 * C1plusC2 * m_R0); + TwoDMatrixType factorOut(C1 * C1plusC2 * R1_Inv); xMLP = factorIn * m_x0 + factorOut * m_x2; yMLP = factorIn * m_y0 + factorOut * m_y2; @@ -180,13 +133,13 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d // This version here is better than the previously implemented one // because it avoids inverting the matrices Sigma. // See comment in [Krah 2018, PMB] - itk::Matrix sum1(R1_Inv * m_Sigma2 + m_Sigma1 * m_R1T); + TwoDMatrixType sum1(R1_Inv * m_Sigma2 + m_Sigma1 * m_R1T); InverseMatrix(sum1); - itk::Matrix sum2(m_R1 * m_Sigma1 + m_Sigma2 * R1T_Inv); + TwoDMatrixType sum2(m_R1 * m_Sigma1 + m_Sigma2 * R1T_Inv); InverseMatrix(sum2); - itk::Matrix part1(R1_Inv * m_Sigma2 * sum1 * m_R0); - itk::Matrix part2(m_Sigma1 * sum2); + TwoDMatrixType part1(R1_Inv * m_Sigma2 * sum1 * m_R0); + TwoDMatrixType part2(m_Sigma1 * sum2); xMLP = part1 * m_x0 + part2 * m_x2; yMLP = part1 * m_y0 + part2 * m_y2; @@ -204,7 +157,7 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d } void -SchulteMLPFunction ::EvaluateError(const double u, itk::Matrix & error) +SchulteMLPFunction ::EvaluateError(const double u, TwoDMatrixType & error) { double x, y; double dx, dy; @@ -225,7 +178,7 @@ SchulteMLPFunction ::PrintTiming(std::ostream & os) #endif void -SchulteMLPFunction ::InverseMatrix(itk::Matrix & mat) +SchulteMLPFunction ::InverseMatrix(TwoDMatrixType & mat) { double det = 1. / (mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0)); std::swap(mat(0, 0), mat(1, 1));