From 8dbba6b8470de20a8be8ca473cfaf52bcc6eb31d Mon Sep 17 00:00:00 2001 From: "Goetz.Stefanie" Date: Wed, 8 Oct 2025 15:30:12 +0100 Subject: [PATCH 01/11] ENH: Implemented helium coefficients Amended by Simon Rit for merging and avoiding inclusion of _ggo.h files --- applications/pctbinning/pctbinning.cxx | 1 + applications/pctbinning/pctbinning.ggo | 1 + include/pctBetheBlochFunctor.h | 18 +++-- ...pctProtonPairsToDistanceDrivenProjection.h | 7 ++ ...tProtonPairsToDistanceDrivenProjection.hxx | 11 ++- include/pctSchulteMLPFunction.h | 68 +++++++++++++++- src/pctSchulteMLPFunction.cxx | 81 ++++++++++++++----- 7 files changed, 157 insertions(+), 30 deletions(-) diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index 5921b6f..b9a29ea 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/applications/pctbinning/pctbinning.cxx @@ -52,6 +52,7 @@ main(int argc, char * argv[]) projection->SetRobust(args_info.robust_flag); projection->SetComputeScattering(args_info.scatwepl_given); projection->SetComputeNoise(args_info.noise_given); + projection->SetParticle(args_info.particle_arg); if (args_info.quadricIn_given) { diff --git a/applications/pctbinning/pctbinning.ggo b/applications/pctbinning/pctbinning.ggo index 9b37507..e5c4c0b 100644 --- a/applications/pctbinning/pctbinning.ggo +++ b/applications/pctbinning/pctbinning.ggo @@ -9,6 +9,7 @@ option "count" c "Image of count of proton pairs per pixel" option "scatwepl" - "Image of scattering WEPL of proton pairs per pixel" string no option "noise" - "Image of WEPL variance per pixel" string no option "robust" r "Use robust estimation of scattering using 19.1 %ile." flag off +option "particle" p "Particle used for imaging (proton or alpha)" string no default="proton" option "source" s "Source position" double no default="0." option "quadricIn" - "Parameters of the entrance quadric support function, see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" double multiple no diff --git a/include/pctBetheBlochFunctor.h b/include/pctBetheBlochFunctor.h index 3631a55..d9d51d7 100644 --- a/include/pctBetheBlochFunctor.h +++ b/include/pctBetheBlochFunctor.h @@ -29,8 +29,13 @@ class BetheBlochProtonStoppingPower ~BetheBlochProtonStoppingPower() {} TOutput - GetValue(const TInput e, const double I) const + GetValue(const TInput e, const double I, const std::string particle) const { + if (particle != "proton") + { + // only implemented for protons + throw std::invalid_argument("Invalid particle type in Bethe Block Functor."); + } /** Physical constants */ static const double K = 4. * CLHEP::pi * CLHEP::classic_electr_radius * CLHEP::classic_electr_radius * CLHEP::electron_mass_c2 * 3.343e+23 / CLHEP::cm3; @@ -55,10 +60,12 @@ template class IntegratedBetheBlochProtonStoppingPowerInverse { public: - IntegratedBetheBlochProtonStoppingPowerInverse(const double I, - const double maxEnergy, - const double binSize = 1. * CLHEP::keV) + IntegratedBetheBlochProtonStoppingPowerInverse(const double I, + const double maxEnergy, + const double binSize = 1. * CLHEP::keV, + const std::string particle = "proton") : m_BinSize(binSize) + , m_Particle(particle) { unsigned int lowBinLimit, numberOfBins; lowBinLimit = itk::Math::Ceil(m_S.GetLowEnergyLimit() / m_BinSize); @@ -71,7 +78,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse } for (unsigned int i = lowBinLimit; i < numberOfBins; i++) { - m_LUT[i] = m_LUT[i - 1] + binSize / m_S.GetValue(TOutput(i) * binSize, I); + m_LUT[i] = m_LUT[i - 1] + binSize / m_S.GetValue(TOutput(i) * binSize, I, m_Particle); // Create inverse lut, i.e., get energy from length in water for (unsigned j = m_Length.size(); j < unsigned(m_LUT[i] * CLHEP::mm) + 1; j++) { @@ -118,6 +125,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse std::vector m_Length; double m_BinSize; + std::string m_Particle; std::vector m_LUT; }; } // end namespace Functor diff --git a/include/pctProtonPairsToDistanceDrivenProjection.h b/include/pctProtonPairsToDistanceDrivenProjection.h index ec02726..03630f4 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.h +++ b/include/pctProtonPairsToDistanceDrivenProjection.h @@ -110,6 +110,10 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection itkGetConstMacro(ComputeNoise, bool); itkBooleanMacro(ComputeNoise); + /** Get/Set the particle type. */ + itkGetMacro(Particle, std::string); + itkSetMacro(Particle, std::string); + protected: ProtonPairsToDistanceDrivenProjection(); virtual ~ProtonPairsToDistanceDrivenProjection() {} @@ -182,6 +186,9 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection bool m_Robust; bool m_ComputeScattering; bool m_ComputeNoise; + + /** The particle type, enum comes from gengetopt header. */ + std::string m_Particle; }; } // end namespace pct diff --git a/include/pctProtonPairsToDistanceDrivenProjection.hxx b/include/pctProtonPairsToDistanceDrivenProjection.hxx index 472fd21..f09caae 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.hxx +++ b/include/pctProtonPairsToDistanceDrivenProjection.hxx @@ -40,8 +40,9 @@ ProtonPairsToDistanceDrivenProjection::BeforeThreaded if (m_QuadricOut.GetPointer() == NULL) m_QuadricOut = m_QuadricIn; + const int zsq = (m_Particle == "alpha") ? 4 : 1; // helium energy is 800 MeV instead of 200 MeV for protons m_ConvFunc = new Functor::IntegratedBetheBlochProtonStoppingPowerInverse( - m_IonizationPotential, 600. * CLHEP::MeV, 0.1 * CLHEP::keV); + m_IonizationPotential, zsq * 600. * CLHEP::MeV, 0.1 * CLHEP::keV, m_Particle); // Read pairs using ReaderType = itk::ImageFileReader; @@ -57,12 +58,16 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera const OutputImageRegionType & itkNotUsed(outputRegionForThread), rtk::ThreadIdType threadId) { + if (m_Particle != "proton" && m_MostLikelyPathType != "schulte") + itkGenericExceptionMacro("Only Schulte's MLP can be used with particles which are not proton"); + // Create MLP depending on type pct::MostLikelyPathFunction::Pointer mlp; if (m_MostLikelyPathType == "polynomial") mlp = pct::ThirdOrderPolynomialMLPFunction::New(); else if (m_MostLikelyPathType == "krah") { + mlp = pct::ThirdOrderPolynomialMLPFunction::New(); pct::PolynomialMLPFunction::Pointer mlp_poly; mlp_poly = pct::PolynomialMLPFunction::New(); mlp_poly->SetPolynomialDegree(m_MostLikelyPathPolynomialDegree); @@ -74,7 +79,9 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera } else if (m_MostLikelyPathType == "schulte") { - mlp = pct::SchulteMLPFunction::New(); + pct::SchulteMLPFunction::Pointer mlpSchulte = pct::SchulteMLPFunction::New(); + mlpSchulte->SetParticle(m_Particle); + mlp = mlpSchulte; } else { diff --git a/include/pctSchulteMLPFunction.h b/include/pctSchulteMLPFunction.h index 93cf6f2..ad429c8 100644 --- a/include/pctSchulteMLPFunction.h +++ b/include/pctSchulteMLPFunction.h @@ -37,6 +37,12 @@ static const double a2 = -9.986645e-08 * aunit / (CLHEP::cm2); static const double a3 = 2.026409e-08 * aunit / (CLHEP::cm3); static const double a4 = -1.420501e-09 * aunit / (CLHEP::cm3 * CLHEP::cm); static const double a5 = 3.899100e-11 * aunit / (CLHEP::cm3 * CLHEP::cm2); +static const double a0He = 4.665553448715765888e-07 * aunit; +static const double a1He = 2.817219469097221029e-08 * aunit / (CLHEP::cm); +static const double a2He = -3.534701631484195519e-09 * aunit / (CLHEP::cm2); +static const double a3He = 7.963860554062791474e-10 * aunit / (CLHEP::cm3); +static const double a4He = -5.635546522420130709e-11 * aunit / (CLHEP::cm3 * CLHEP::cm); +static const double a5He = 1.635484142530505741e-12 * aunit / (CLHEP::cm3 * CLHEP::cm2); #elif defined(NICO_COEFF) static const double a0 = 3.599e-06 * aunit; static const double a1 = 7.303e-08 * aunit / (CLHEP::cm); @@ -69,6 +75,11 @@ class ConstantPartOfIntegrals double v = 1. + 0.038 * std::log((diffU)*invX0); return c * v * v; } + static double + GetValueHe(const double ux, const double uy) + { + return GetValue(ux, uy); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 8 @@ -88,6 +99,20 @@ class IntegralForSigmaSqTheta static const double a5Theta = a5 / 6.; return u * (a0Theta + u * (a1Theta + u * (a2Theta + u * (a3Theta + u * (a4Theta + u * a5Theta))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0Theta = a0He; + static const double a1Theta = a1He / 2.; + static const double a2Theta = a2He / 3.; + static const double a3Theta = a3He / 4.; + static const double a4Theta = a4He / 5.; + static const double a5Theta = a5He / 6.; + return u * (a0Theta + u * (a1Theta + u * (a2Theta + u * (a3Theta + u * (a4Theta + u * a5Theta))))); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 9 @@ -107,6 +132,20 @@ class IntegralForSigmaSqTTheta static const double a5TTheta = a5 / 7.; return u * u * (a0TTheta + u * (a1TTheta + u * (a2TTheta + u * (a3TTheta + u * (a4TTheta + u * a5TTheta))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0TTheta = a0He / 2.; + static const double a1TTheta = a1He / 3.; + static const double a2TTheta = a2He / 4.; + static const double a3TTheta = a3He / 5.; + static const double a4TTheta = a4He / 6.; + static const double a5TTheta = a5He / 7.; + return u * u * (a0TTheta + u * (a1TTheta + u * (a2TTheta + u * (a3TTheta + u * (a4TTheta + u * a5TTheta))))); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 7 @@ -126,6 +165,20 @@ class IntegralForSigmaSqT static const double a5T = a5 / 8.; return u * u * u * (a0T + u * (a1T + u * (a2T + u * (a3T + u * (a4T + u * a5T))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0T = a0He / 3.; + static const double a1T = a1He / 4.; + static const double a2T = a2He / 5.; + static const double a3T = a3He / 6.; + static const double a4T = a4He / 7.; + static const double a5T = a5He / 8.; + return u * u * u * (a0T + u * (a1T + u * (a2T + u * (a3T + u * (a4T + u * a5T))))); + } }; } // end namespace SchulteMLP @@ -177,6 +230,18 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction void EvaluateError(const double u1, itk::Matrix & error); + /** Get/Set the particle type. */ + void + SetParticle(const std::string particle) + { + m_Particle = particle; + } + std::string + GetParticle() + { + return m_Particle; + } + #ifdef MLP_TIMING /** Print timing information */ virtual void @@ -239,6 +304,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction double m_IntForSigmaSqTTheta2; double m_IntForSigmaSqT2; + std::string m_Particle{ "proton" }; + #ifdef MLP_TIMING itk::TimeProbe m_EvaluateProbe1; itk::TimeProbe m_EvaluateProbe2; @@ -247,5 +314,4 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction } // end namespace pct - #endif diff --git a/src/pctSchulteMLPFunction.cxx b/src/pctSchulteMLPFunction.cxx index 7522d7f..e575a56 100644 --- a/src/pctSchulteMLPFunction.cxx +++ b/src/pctSchulteMLPFunction.cxx @@ -105,9 +105,22 @@ SchulteMLPFunction ::Init(const PointType posIn, m_uOrigin = posIn[2]; m_u2 = posOut[2] - m_uOrigin; - m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(m_u2); - m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(m_u2); - m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(m_u2); + if (m_Particle == "proton") + { + m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(m_u2); + m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(m_u2); + m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(m_u2); + } + else if (m_Particle == "alpha") + { + m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValueHe(m_u2); + m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValueHe(m_u2); + m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValueHe(m_u2); + } + else + { + throw std::invalid_argument("Invalid particle in Schulte MLP function."); + } // Parameters vectors m_x0[0] = posIn[0]; @@ -141,27 +154,51 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d InverseMatrix(R1_Inv); // Constants used in both integrals - const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(u1); - const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(u1); - const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(u1); - - // Construct Sigma1 (equations 6-9) - m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; - m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; - m_Sigma1(1, 0) = m_Sigma1(0, 1); - m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; - m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(m_u0, u1); - - double sigma1 = std::sqrt(m_Sigma2(1, 1)); - - // Construct Sigma2 (equations 15-18) - m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; - m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; - m_Sigma2(1, 0) = m_Sigma2(0, 1); - m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; - m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1, m_u2); + if (m_Particle == "proton") + { + const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(u1); + const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(u1); + const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(u1); + + // Construct Sigma1 (equations 6-9) + m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; + m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; + m_Sigma1(1, 0) = m_Sigma1(0, 1); + m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; + m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(m_u0, u1); + + double sigma1 = std::sqrt(m_Sigma2(1, 1)); + + // Construct Sigma2 (equations 15-18) + m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; + m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; + m_Sigma2(1, 0) = m_Sigma2(0, 1); + m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; + m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1, m_u2); + } + else if (m_Particle == "alpha") + { + const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValueHe(u1); + const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValueHe(u1); + const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValueHe(u1); + + // Construct Sigma1 (equations 6-9) + m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; + m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; + m_Sigma1(1, 0) = m_Sigma1(0, 1); + m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; + m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValueHe(m_u0, u1); + + // Construct Sigma2 (equations 15-18) + m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; + m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; + m_Sigma2(1, 0) = m_Sigma2(0, 1); + m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; + m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValueHe(u1, m_u2); + } #ifdef MLP_TIMING + m_EvaluateProbe1.Stop(); m_EvaluateProbe2.Start(); #endif From 18dcd5c0ba6e7228384ef2e61a028d5dd42056e5 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 8 Oct 2025 20:48:04 +0100 Subject: [PATCH 02/11] ENH: Add LMU variance reconstruction code --- applications/pctbinning/pctbinning.cxx | 33 ++ applications/pctbinning/pctbinning.ggo | 2 + applications/pctfdk/pctfdk.cxx | 13 +- applications/pctfdk/pctfdk.ggo | 1 + ...DKDDConeBeamVarianceReconstructionFilter.h | 68 ++++ ...DDConeBeamVarianceReconstructionFilter.hxx | 33 ++ include/pctFDKDDWeightProjectionFilter.h | 6 + include/pctFDKDDWeightProjectionFilter.hxx | 5 + include/pctFFTVarianceImageFilter.h | 96 ++++++ include/pctFFTVarianceImageFilter.hxx | 311 ++++++++++++++++++ ...pctProtonPairsToDistanceDrivenProjection.h | 19 ++ ...tProtonPairsToDistanceDrivenProjection.hxx | 98 +++++- 12 files changed, 683 insertions(+), 2 deletions(-) create mode 100644 include/pctFDKDDConeBeamVarianceReconstructionFilter.h create mode 100644 include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx create mode 100644 include/pctFFTVarianceImageFilter.h create mode 100644 include/pctFFTVarianceImageFilter.hxx diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index b9a29ea..126542f 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/applications/pctbinning/pctbinning.cxx @@ -52,6 +52,7 @@ main(int argc, char * argv[]) projection->SetRobust(args_info.robust_flag); projection->SetComputeScattering(args_info.scatwepl_given); projection->SetComputeNoise(args_info.noise_given); + projection->SetComputeVariance(args_info.variance_given); projection->SetParticle(args_info.particle_arg); if (args_info.quadricIn_given) @@ -129,6 +130,38 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(cwriter->Update()) } + if (args_info.variance_given) + { + + SmallHoleFiller fillerVar; + if (args_info.variance_given && args_info.fillvariance_flag) + { + fillerVar.SetImage(projection->GetVariance()); + fillerVar.SetHolePixel(0.); + fillerVar.Fill(); + } + + typedef itk::ChangeInformationImageFilter CIIVarType; + CIIVarType::Pointer ciiVar = CIIVarType::New(); + if (args_info.fillvariance_flag) + ciiVar->SetInput(fillerVar.GetOutput()); + else + ciiVar->SetInput(projection->GetVariance()); + ciiVar->ChangeOriginOn(); + ciiVar->ChangeDirectionOn(); + ciiVar->ChangeSpacingOn(); + ciiVar->SetOutputDirection(projection->GetOutput()->GetDirection()); + ciiVar->SetOutputOrigin(projection->GetOutput()->GetOrigin()); + ciiVar->SetOutputSpacing(projection->GetOutput()->GetSpacing()); + + // Write + typedef itk::ImageFileWriter VarianceWriterType; + VarianceWriterType::Pointer vwriter = VarianceWriterType::New(); + vwriter->SetFileName(args_info.variance_arg); + vwriter->SetInput(ciiVar->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(vwriter->Update()) + } + if (args_info.scatwepl_given) { // Write diff --git a/applications/pctbinning/pctbinning.ggo b/applications/pctbinning/pctbinning.ggo index e5c4c0b..7685f72 100644 --- a/applications/pctbinning/pctbinning.ggo +++ b/applications/pctbinning/pctbinning.ggo @@ -6,6 +6,7 @@ option "input" i "Input file name containing the proton pairs" option "output" o "Output file name" string no option "elosswepl" - "Output file name (alias for --output)" string no option "count" c "Image of count of proton pairs per pixel" string no +option "variance" - "Image of variance of WEPL values per pixel" string no option "scatwepl" - "Image of scattering WEPL of proton pairs per pixel" string no option "noise" - "Image of WEPL variance per pixel" string no option "robust" r "Use robust estimation of scattering using 19.1 %ile." flag off @@ -19,6 +20,7 @@ option "mlptrackeruncert" - "Consider tracker uncertainties in MLP [Krah 2018 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 "fillvariance" - "Fill holes, i.e. pixels that were not hit by protons (for variance projections)" 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 diff --git a/applications/pctfdk/pctfdk.cxx b/applications/pctfdk/pctfdk.cxx index bd4b2e7..809673f 100644 --- a/applications/pctfdk/pctfdk.cxx +++ b/applications/pctfdk/pctfdk.cxx @@ -5,6 +5,7 @@ #include "rtkProjectionsReader.h" #include "pctDDParkerShortScanImageFilter.h" #include "pctFDKDDConeBeamReconstructionFilter.h" +#include "pctFDKDDConeBeamVarianceReconstructionFilter.h" #include #include @@ -62,7 +63,17 @@ main(int argc, char * argv[]) rtk::SetConstantImageSourceFromGgo(constantImageSource, args_info); // FDK reconstruction filtering - auto feldkamp = pct::FDKDDConeBeamReconstructionFilter::New(); + using FDKCPUType = pct::FDKDDConeBeamReconstructionFilter; + FDKCPUType::Pointer feldkamp = nullptr; + if (args_info.variance_flag) + { + feldkamp = pct::FDKDDConeBeamVarianceReconstructionFilter::New(); + } + else + { + feldkamp = FDKCPUType::New(); + } + feldkamp->SetInput(0, constantImageSource->GetOutput()); feldkamp->SetProjectionStack(pssf->GetOutput()); feldkamp->SetGeometry(geometryReader->GetOutputObject()); diff --git a/applications/pctfdk/pctfdk.ggo b/applications/pctfdk/pctfdk.ggo index d4b6bb1..8d963e6 100644 --- a/applications/pctfdk/pctfdk.ggo +++ b/applications/pctfdk/pctfdk.ggo @@ -13,6 +13,7 @@ section "Ramp filter" option "pad" - "Data padding parameter to correct for truncation" double no default="0.0" option "hann" - "Cut frequency for hann window in ]0,1] (0.0 disables it)" double no default="0.0" option "hannY" - "Cut frequency for hann window in ]0,1] (0.0 disables it)" double no default="0.0" +option "variance" - "Use variance reconstruction kernel" flag off section "Volume properties" option "origin" - "Origin (default=centered)" double multiple no diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h new file mode 100644 index 0000000..5b8d14a --- /dev/null +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h @@ -0,0 +1,68 @@ +#ifndef __pctFDKDDConeBeamVarianceReconstructionFilter_h +#define __pctFDKDDConeBeamVarianceReconstructionFilter_h + +#include "pctFFTVarianceImageFilter.h" +#include "pctFDKDDConeBeamReconstructionFilter.h" +#include "pctFFTVarianceImageFilter.h" + +#include +#include + +/** \class FDKDDConeBeamReconstructionFilter + * TODO + * + * \author Jannis Dickmann + */ +namespace pct +{ + +template +class ITK_EXPORT FDKDDConeBeamVarianceReconstructionFilter + : public pct::FDKDDConeBeamReconstructionFilter +{ +public: + typedef pct::FDKDDConeBeamReconstructionFilter Baseclass; + + /** Standard class typedefs. */ + typedef FDKDDConeBeamVarianceReconstructionFilter Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Some convenient typedefs. */ + typedef typename Baseclass::InputImageType InputImageType; + typedef typename Baseclass::OutputImageType OutputImageType; + + typedef typename Baseclass::ProjectionStackType ProjectionStackType; + typedef typename Baseclass::ProjectionStackPointer ProjectionStackPointer; + + /** Typedefs of each subfilter of this composite filter */ + typedef typename Baseclass::ExtractFilterType ExtractFilterType; + typedef typename Baseclass::WeightFilterType WeightFilterType; + typedef pct::FFTVarianceImageFilter VarianceFilterType; + typedef typename Baseclass::BackProjectionFilterType BackProjectionFilterType; + + /** Standard New method. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkOverrideGetNameOfClassMacro(FDKDDConeBeamVarianceReconstructionFilter); + +protected: + FDKDDConeBeamVarianceReconstructionFilter(); + ~FDKDDConeBeamVarianceReconstructionFilter() {} + +private: + // purposely not implemented + FDKDDConeBeamVarianceReconstructionFilter(const Self &); + void + operator=(const Self &); +}; // end of class + +} // end namespace pct + +#ifndef ITK_MANUAL_INSTANTIATION +# include "pctFDKDDConeBeamVarianceReconstructionFilter.hxx" +#endif + +#endif diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx b/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx new file mode 100644 index 0000000..0daa2cc --- /dev/null +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx @@ -0,0 +1,33 @@ +namespace pct +{ + +template +FDKDDConeBeamVarianceReconstructionFilter:: + FDKDDConeBeamVarianceReconstructionFilter() +{ + this->SetNumberOfRequiredInputs(1); + + // Create each filter of the composite filter + this->m_ExtractFilter = ExtractFilterType::New(); + this->m_WeightFilter = WeightFilterType::New(); + this->m_RampFilter = VarianceFilterType::New(); + this->m_BackProjectionFilter = BackProjectionFilterType::New(); + + // Permanent internal connections + this->m_WeightFilter->SetInput(this->m_ExtractFilter->GetOutput()); + this->m_RampFilter->SetInput(this->m_WeightFilter->GetOutput()); + this->m_BackProjectionFilter->SetProjectionStack(this->m_RampFilter->GetOutput()); + + // Default parameters +#if ITK_VERSION_MAJOR >= 4 + this->m_ExtractFilter->SetDirectionCollapseToSubmatrix(); +#endif + this->m_WeightFilter->InPlaceOn(); + this->m_BackProjectionFilter->InPlaceOn(); + this->m_BackProjectionFilter->SetTranspose(true); + + this->m_WeightFilter->SetVarianceReconstruction(true); +} + + +} // end namespace pct diff --git a/include/pctFDKDDWeightProjectionFilter.h b/include/pctFDKDDWeightProjectionFilter.h index 5d4ad2d..4869f91 100644 --- a/include/pctFDKDDWeightProjectionFilter.h +++ b/include/pctFDKDDWeightProjectionFilter.h @@ -46,6 +46,10 @@ class ITK_TEMPLATE_EXPORT FDKDDWeightProjectionFilter : public itk::InPlaceImage itkGetConstObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); itkSetObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); + /** Get/ Set variance mode */ + itkGetMacro(VarianceReconstruction, bool); + itkSetMacro(VarianceReconstruction, bool); + protected: FDKDDWeightProjectionFilter() {} ~FDKDDWeightProjectionFilter() {} @@ -55,6 +59,8 @@ class ITK_TEMPLATE_EXPORT FDKDDWeightProjectionFilter : public itk::InPlaceImage virtual void DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + bool m_VarianceReconstruction = false; + private: FDKDDWeightProjectionFilter(const Self &); // purposely not implemented void diff --git a/include/pctFDKDDWeightProjectionFilter.hxx b/include/pctFDKDDWeightProjectionFilter.hxx index e609814..2d01fee 100644 --- a/include/pctFDKDDWeightProjectionFilter.hxx +++ b/include/pctFDKDDWeightProjectionFilter.hxx @@ -27,6 +27,11 @@ FDKDDWeightProjectionFilter::BeforeThreadedGenerateDa const double rampFactor = sdd / (2. * sid); m_AngularWeightsAndRampFactor[k] *= rampFactor; } + if (this->m_VarianceReconstruction) + { + // square if in variance mode + m_AngularWeightsAndRampFactor[k] = m_AngularWeightsAndRampFactor[k] * m_AngularWeightsAndRampFactor[k]; + } } } diff --git a/include/pctFFTVarianceImageFilter.h b/include/pctFFTVarianceImageFilter.h new file mode 100644 index 0000000..dc6918b --- /dev/null +++ b/include/pctFFTVarianceImageFilter.h @@ -0,0 +1,96 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef pctFFTVarianceImageFilter_h +#define pctFFTVarianceImageFilter_h + +#include "rtkFFTRampImageFilter.h" +#include "rtkFFTProjectionsConvolutionImageFilter.h" +#include "itkHalfHermitianToRealInverseFFTImageFilter.h" +#include "itkImageFileWriter.h" + +#include +#include + +namespace pct +{ + +/** \class FFTVarianceImageFilter + * \brief Implements the variance image filter of the filtered backprojection algorithm. + * + * The filter code is based on FFTConvolutionImageFilter by Gaetan Lehmann + * (see http://hdl.handle.net/10380/3154) + * + * \test rtkvariancefiltertest.cxx + * + * \author Simon Rit + * + * \ingroup ImageToImageFilter + */ + +template +class ITK_EXPORT FFTVarianceImageFilter : public rtk::FFTRampImageFilter +{ +public: + typedef rtk::FFTRampImageFilter Baseclass; + + /** Standard class typedefs. */ + typedef FFTVarianceImageFilter Self; + typedef rtk::FFTRampImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Some convenient typedefs. */ + typedef typename Baseclass::InputImageType InputImageType; + typedef typename Baseclass::OutputImageType OutputImageType; + typedef typename Baseclass::FFTPrecisionType FFTPrecisionType; + typedef typename Baseclass::IndexType IndexType; + typedef typename Baseclass::SizeType SizeType; + + typedef typename Baseclass::FFTInputImageType FFTInputImageType; + typedef typename Baseclass::FFTInputImagePointer FFTInputImagePointer; + typedef typename Baseclass::FFTOutputImageType FFTOutputImageType; + typedef typename Baseclass::FFTOutputImagePointer FFTOutputImagePointer; + + /** Standard New method. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkOverrideGetNameOfClassMacro(FFTVarianceImageFilter); + +protected: + FFTVarianceImageFilter(); + ~FFTVarianceImageFilter() {} + + /** Creates and return a pointer to one line of the variance kernel in Fourier space. + * Used in generate data functions. */ + void + UpdateFFTProjectionsConvolutionKernel(const SizeType size) override; + + FFTVarianceImageFilter(const Self &); // purposely not implemented + void + operator=(const Self &); // purposely not implemented +}; // end of class + +} // namespace pct + +#ifndef ITK_MANUAL_INSTANTIATION +# include "pctFFTVarianceImageFilter.hxx" +#endif + +#endif diff --git a/include/pctFFTVarianceImageFilter.hxx b/include/pctFFTVarianceImageFilter.hxx new file mode 100644 index 0000000..19be1b1 --- /dev/null +++ b/include/pctFFTVarianceImageFilter.hxx @@ -0,0 +1,311 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef pctFFTVarianceImageFilter_hxx +#define pctFFTVarianceImageFilter_hxx + +namespace pct +{ + +template +FFTVarianceImageFilter::FFTVarianceImageFilter() +{ + // this->m_HannCutFrequency = 0.; + // this->m_CosineCutFrequency = 0.; + // this->m_HammingFrequency = 0.; + // this->m_HannCutFrequency = 0.; + // this->m_HannCutFrequencyY = 0.; + // this->m_RamLakCutFrequency = 0.; + // this->m_SheppLoganCutFrequency = 0.; +} + +template +void +FFTVarianceImageFilter::UpdateFFTProjectionsConvolutionKernel( + const SizeType s) +{ + // if(this->m_KernelFFT.GetPointer() != ITK_NULLPTR && s == this->m_PreviousKernelUpdateSize) + // { + // return; + // } + // this->m_PreviousKernelUpdateSize = s; + + const int width = s[0]; + const int height = s[1]; + + // Allocate kernel + SizeType size; + size.Fill(1); + size[0] = width; + FFTInputImagePointer kernel = FFTInputImageType::New(); + kernel->SetRegions(size); + kernel->Allocate(); + kernel->FillBuffer(0.); + + // Compute kernel in space domain (see Kak & Slaney, chapter 3 equation 61 + // page 72) although spacing is not squared according to equation 69 page 75 + double spacing = this->GetInput()->GetSpacing()[0]; + IndexType ix, jx; + ix.Fill(0); + jx.Fill(0); + kernel->SetPixel(ix, 1. / (4. * spacing)); + for (ix[0] = 1, jx[0] = size[0] - 1; ix[0] < typename IndexType::IndexValueType(size[0] / 2); ix[0] += 2, jx[0] -= 2) + { + double v = ix[0] * vnl_math::pi; + v = -1. / (v * v * spacing); + kernel->SetPixel(ix, v); + kernel->SetPixel(jx, v); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // typedef itk::ImageRegionIteratorWithIndex InverseIteratorType2; + // InverseIteratorType2 itiK2(kernel, kernel->GetLargestPossibleRegion()); + // std::ofstream out0("kernel_real.csv"); + // itiK2.GoToBegin(); + // for(; !itiK2.IsAtEnd(); ++itiK2) + // { + // out0 << itiK2.Get() << std::endl; + // } + // out0.close(); + // // DEBUG + + // FFT kernel + typedef itk::RealToHalfHermitianForwardFFTImageFilter FFTType; + typename FFTType::Pointer fftK = FFTType::New(); + fftK->SetInput(kernel); + fftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK->Update(); + this->m_KernelFFT = fftK->GetOutput(); + + // Windowing (if enabled) + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType itK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion()); + + unsigned int n = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); + + itK.GoToBegin(); + if (this->GetHannCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetHannCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = 0.5 * (1 + std::cos(itk::Math::pi * i / ncut)); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetCosineCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetCosineCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = std::cos(0.5 * itk::Math::pi * i / ncut); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetHammingFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetHammingFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = 0.54 + 0.46 * (std::cos(itk::Math::pi * i / ncut)); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetRamLakCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetRamLakCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + } + } + else if (this->GetSheppLoganCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetSheppLoganCutFrequency())); + // sinc(0) --> is 1 + ++itK; + for (unsigned int i = 1; i < ncut; i++, ++itK) + { + double x = 0.5 * vnl_math::pi * i / ncut; + double k = std::sin(x) / x; + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else + { + itK.GoToReverseBegin(); + ++itK; + } + + for (; !itK.IsAtEnd(); ++itK) + { + itK.Set(itK.Get() * TFFTPrecision(0.)); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out1("kernel_fourier.csv"); + // itK.GoToBegin(); + // for(; !itK.IsAtEnd(); ++itK) + // { + // out1 << itK.Get().real() << " " << itK.Get().imag() << std::endl; + // } + // out1.close(); + // // DEBUG + + // iFFT kernel + typedef itk::HalfHermitianToRealInverseFFTImageFilter IFFTType; + typename IFFTType::Pointer ifftK = IFFTType::New(); + ifftK->SetInput(this->m_KernelFFT); + ifftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + ifftK->Update(); + typename FFTType::InputImageType::Pointer KernelIFFT = ifftK->GetOutput(); + + // calculate ratio g_C/g^2 + typedef itk::ImageRegionIteratorWithIndex InverseIteratorType; + InverseIteratorType itiK(KernelIFFT, KernelIFFT->GetLargestPossibleRegion()); + typename FFTType::InputImageType::PixelType sumgc = 0.; + typename FFTType::InputImageType::PixelType sumg2 = 0.; + typename FFTType::InputImageType::IndexType idxshifted; + const unsigned int maxidx = KernelIFFT->GetLargestPossibleRegion().GetSize()[0]; + for (; !itiK.IsAtEnd(); ++itiK) + { + idxshifted = itiK.GetIndex(); + if (idxshifted[0] == 0) + idxshifted[0] = maxidx - 1; + else + idxshifted[0] -= 1; + + sumgc += itiK.Get() * KernelIFFT->GetPixel(idxshifted); + sumg2 += itiK.Get() * itiK.Get(); + } + const typename FFTType::InputImageType::PixelType ratiogcg2 = sumgc / sumg2; + + // numerical integration to calculate f_interp + const double aprecision = 0.00001; + double finterp = 0.; + for (unsigned int i = 0; i < int(1. / aprecision); i++) + { + const double a = double(i) * aprecision; + finterp += (1 - a) * (1 - a) + 2 * ratiogcg2 * (1 - a) * a + a * a; + } + finterp *= aprecision; + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out9("kernel_window_real.csv"); + // itiK.GoToBegin(); + // for(; !itiK.IsAtEnd(); ++itiK) + // { + // out9 << itiK.Get() << std::endl; + // } + // out9.close(); + // // DEBUG + + // square kernel and multiply with finterp + itiK.GoToBegin(); + for (; !itiK.IsAtEnd(); ++itiK) + { + itiK.Set(itiK.Get() * itiK.Get() * finterp); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out2("kernel_sq_real.csv"); + // itiK.GoToBegin(); + // for(; !itiK.IsAtEnd(); ++itiK) + // { + // out2 << itiK.Get() << std::endl; + // } + // out2.close(); + // // DEBUG + + // FFT kernel + typename FFTType::Pointer fftK2 = FFTType::New(); + fftK2->SetInput(ifftK->GetOutput()); + fftK2->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK2->Update(); + this->m_KernelFFT = fftK2->GetOutput(); + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out3("kernel_sq_fourier.csv"); + // IteratorType itK2(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion() ); + // itK2.GoToBegin(); + // for(; !itK2.IsAtEnd(); ++itK2) + // { + // out3 << itK2.Get().real() << " " << itK2.Get().imag() << std::endl; + // } + // out3.close(); + // // DEBUG + + // Replicate and window if required + if (this->GetHannCutFrequencyY() > 0.) + { + std::cerr << "HannY not implemented for variance image filter." << std::endl; + + size.Fill(1); + size[0] = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); + size[1] = height; + + const unsigned int ncut = itk::Math::Round((height / 2 + 1) * std::min(1.0, this->GetHannCutFrequencyY())); + + this->m_KernelFFT = FFTOutputImageType::New(); + this->m_KernelFFT->SetRegions(size); + this->m_KernelFFT->Allocate(); + this->m_KernelFFT->FillBuffer(0.); + + IteratorType itTwoDK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion()); + for (unsigned int j = 0; j < ncut; j++) + { + itK.GoToBegin(); + const TFFTPrecision win(0.5 * (1 + std::cos(itk::Math::pi * j / ncut))); + for (unsigned int i = 0; i < size[0]; ++itK, ++itTwoDK, i++) + { + itTwoDK.Set(win * itK.Get()); + } + } + itTwoDK.GoToReverseBegin(); + for (unsigned int j = 1; j < ncut; j++) + { + itK.GoToReverseBegin(); + const TFFTPrecision win(0.5 * (1 + std::cos(itk::Math::pi * j / ncut))); + for (unsigned int i = 0; i < size[0]; --itK, --itTwoDK, i++) + { + itTwoDK.Set(win * itK.Get()); + } + } + } + + this->m_KernelFFT->DisconnectPipeline(); + + // ////////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out4("finterp.csv"); + // out4 << finterp << std::endl; + // out4 << ratiogcg2 << std::endl; + // out4.close(); + // std::cerr << std::endl << "finterp = " << finterp << std::endl; + // std::cerr << "ratiogcg2 = " << ratiogcg2 << std::endl; + // std::cerr << "Aborting here since filters were alrady written out." << std::endl; + // exit(9); + // // DEBUG +} + +} // namespace pct +#endif diff --git a/include/pctProtonPairsToDistanceDrivenProjection.h b/include/pctProtonPairsToDistanceDrivenProjection.h index 03630f4..08c0a23 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.h +++ b/include/pctProtonPairsToDistanceDrivenProjection.h @@ -29,6 +29,9 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection using CountImageType = itk::Image; using CountImagePointer = CountImageType::Pointer; + typedef itk::Image VarianceImageType; + typedef VarianceImageType::Pointer VarianceImagePointer; + using AngleImageType = itk::Image; using AngleImagePointer = AngleImageType::Pointer; @@ -78,6 +81,9 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection /** Get/Set the count of proton pairs per pixel. */ itkGetModifiableObjectMacro(Count, CountImageType); + /** Get/Set the variance of proton pairs per pixel. */ + itkGetMacro(Variance, VarianceImagePointer); + /** Get/Set the angle of proton pairs per pixel. */ itkGetModifiableObjectMacro(Angle, AngleImageType); @@ -110,6 +116,12 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection itkGetConstMacro(ComputeNoise, bool); itkBooleanMacro(ComputeNoise); + /** Do the variance projections + ** Default is off. */ + itkSetMacro(ComputeVariance, bool); + itkGetConstMacro(ComputeVariance, bool); + itkBooleanMacro(ComputeVariance); + /** Get/Set the particle type. */ itkGetMacro(Particle, std::string); itkSetMacro(Particle, std::string); @@ -154,6 +166,12 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection CountImagePointer m_Count; std::vector m_Counts; + /** Calculate variance in each thread */ + VarianceImagePointer m_Variance; + std::vector m_Variances; + VarianceImagePointer m_VarianceSq; + std::vector m_VariancesSq; + AngleImagePointer m_Angle; std::vector m_Angles; std::vector> m_AnglesVectors; @@ -185,6 +203,7 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection ProtonPairsImageType::Pointer m_ProtonPairs; bool m_Robust; bool m_ComputeScattering; + bool m_ComputeVariance; bool m_ComputeNoise; /** The particle type, enum comes from gengetopt header. */ diff --git a/include/pctProtonPairsToDistanceDrivenProjection.hxx b/include/pctProtonPairsToDistanceDrivenProjection.hxx index f09caae..0d23463 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.hxx +++ b/include/pctProtonPairsToDistanceDrivenProjection.hxx @@ -15,6 +15,7 @@ ProtonPairsToDistanceDrivenProjection::ProtonPairsToD : m_Robust(false) , m_ComputeScattering(false) , m_ComputeNoise(false) + , m_ComputeVariance(false) { this->DynamicMultiThreadingOff(); this->SetNumberOfWorkUnits(itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads()); @@ -26,6 +27,13 @@ ProtonPairsToDistanceDrivenProjection::BeforeThreaded { m_Outputs.resize(this->GetNumberOfWorkUnits()); m_Counts.resize(this->GetNumberOfWorkUnits()); + + if (m_ComputeVariance) + { + m_Variances.resize(this->GetNumberOfWorkUnits()); + m_VariancesSq.resize(this->GetNumberOfWorkUnits()); + } + if (m_ComputeScattering) { m_Angles.resize(this->GetNumberOfWorkUnits()); @@ -99,6 +107,19 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera m_Counts[threadId]->Allocate(); m_Counts[threadId]->FillBuffer(0); + if (m_ComputeVariance) + { + m_Variances[threadId] = VarianceImageType::New(); + m_Variances[threadId]->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_Variances[threadId]->Allocate(); + m_Variances[threadId]->FillBuffer(0); + + m_VariancesSq[threadId] = VarianceImageType::New(); + m_VariancesSq[threadId]->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_VariancesSq[threadId]->Allocate(); + m_VariancesSq[threadId]->FillBuffer(0); + } + if (m_ComputeScattering && (!m_Robust || threadId == 0)) // Note NK: is this condition correct? Should it not be !(m_Robust || threadId==0) ? { @@ -125,6 +146,11 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera { m_Outputs[0] = this->GetOutput(); m_Count = m_Counts[0]; + if (m_ComputeVariance) + { + m_Variance = m_Variances[0]; + m_VarianceSq = m_VariancesSq[0]; + } if (m_ComputeScattering) { m_Angle = m_Angles[0]; @@ -159,7 +185,15 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera typename OutputImageType::PixelType * imgData = m_Outputs[threadId]->GetBufferPointer(); typename OutputImageType::PixelType * imgSquaredData = NULL; unsigned int * imgCountData = m_Counts[threadId]->GetBufferPointer(); - float * imgAngleData = NULL, *imgAngleSqData = NULL; + + float *imgVarianceData = NULL, *imgVarianceSqData = NULL; + if (m_ComputeVariance) + { + imgVarianceData = m_Variances[threadId]->GetBufferPointer(); + imgVarianceSqData = m_VariancesSq[threadId]->GetBufferPointer(); + } + + float *imgAngleData = NULL, *imgAngleSqData = NULL; if (m_ComputeScattering && !m_Robust) { imgAngleData = m_Angles[threadId]->GetBufferPointer(); @@ -438,6 +472,11 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera { imgSquaredData[idx] += value * value; } + if (m_ComputeVariance) + { + imgVarianceData[idx] += value; + imgVarianceSqData[idx] += value * value; + } imgCountData[idx]++; if (m_ComputeScattering) { @@ -552,6 +591,61 @@ ProtonPairsToDistanceDrivenProjection::AfterThreadedG } } + if (m_ComputeVariance) + { + typedef itk::ImageRegionIterator ImageVarianceIteratorType; + ImageVarianceIteratorType itVarianceOut(m_Variances[0], m_Outputs[0]->GetLargestPossibleRegion()); + + typedef itk::ImageRegionIterator ImageVarianceSqIteratorType; + ImageVarianceSqIteratorType itVarianceSqOut(m_VariancesSq[0], m_Outputs[0]->GetLargestPossibleRegion()); + + for (unsigned int i = 1; i < this->GetNumberOfWorkUnits(); i++) + { + if (m_Outputs[i].GetPointer() == NULL) + continue; + ImageVarianceIteratorType itVarianceOutThread(m_Variances[i], m_Outputs[i]->GetLargestPossibleRegion()); + ImageVarianceSqIteratorType itVarianceSqOutThread(m_VariancesSq[i], m_Outputs[i]->GetLargestPossibleRegion()); + + while (!itVarianceOut.IsAtEnd()) + { + itVarianceOut.Set(itVarianceOut.Get() + itVarianceOutThread.Get()); + ++itVarianceOutThread; + ++itVarianceOut; + + itVarianceSqOut.Set(itVarianceSqOut.Get() + itVarianceSqOutThread.Get()); + ++itVarianceSqOutThread; + ++itVarianceSqOut; + } + itVarianceOut.GoToBegin(); + itVarianceSqOut.GoToBegin(); + } + + // Set variance image information + m_Variance->SetSpacing(this->GetOutput()->GetSpacing()); + m_Variance->SetOrigin(this->GetOutput()->GetOrigin()); + + itCOut.GoToBegin(); + while (!itCOut.IsAtEnd()) + { + if (itCOut.Get() >= 2) + { + const float vsum = itVarianceOut.Get(); + const float vsqsum = itVarianceSqOut.Get(); + const float n = itCOut.Get(); + const float variance_wepl = (vsqsum - vsum * vsum / n) / (n - 1); + const float variance_proj = variance_wepl / n; + itVarianceOut.Set(variance_proj); + } + else + { + itVarianceOut.Set(0.); + } + ++itCOut; + ++itVarianceOut; + ++itVarianceSqOut; + } + } + if (m_ComputeScattering) { using ImageAngleIteratorType = itk::ImageRegionIterator; @@ -632,6 +726,8 @@ ProtonPairsToDistanceDrivenProjection::AfterThreadedG m_Outputs.resize(0); m_SquaredOutputs.resize(0); m_Counts.resize(0); + m_Variances.resize(0); + m_VariancesSq.resize(0); m_Angles.resize(0); m_AnglesSq.resize(0); m_AnglesVectors.resize(0); From 5ba3c828865295292522f45c40ddd2e93010c499 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 8 Oct 2025 20:51:54 +0100 Subject: [PATCH 03/11] ENH: Adapt pctbinning to use HoleFillingImageFilter --- applications/pctbinning/pctbinning.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index 126542f..38b6180 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/applications/pctbinning/pctbinning.cxx @@ -133,18 +133,17 @@ main(int argc, char * argv[]) if (args_info.variance_given) { - SmallHoleFiller fillerVar; + auto varianceHoleFilter = pct::HoleFillingImageFilter::New(); if (args_info.variance_given && args_info.fillvariance_flag) { - fillerVar.SetImage(projection->GetVariance()); - fillerVar.SetHolePixel(0.); - fillerVar.Fill(); + varianceHoleFilter->SetInput(projection->GetVariance()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(varianceHoleFilter->Update()); } typedef itk::ChangeInformationImageFilter CIIVarType; CIIVarType::Pointer ciiVar = CIIVarType::New(); if (args_info.fillvariance_flag) - ciiVar->SetInput(fillerVar.GetOutput()); + ciiVar->SetInput(varianceHoleFilter->GetOutput()); else ciiVar->SetInput(projection->GetVariance()); ciiVar->ChangeOriginOn(); From 9fd28a8fb2be84ab46e2a3750a9c25c8fade610e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 3 Apr 2026 15:55:15 +0200 Subject: [PATCH 04/11] ENH: Convert pctpairprotonsLomaLinda to Python application --- applications/pctlomalinda/pctlomalinda.py | 129 ++++++++++++++++++++++ pyproject.toml | 2 +- wrapping/__init_pct__.py | 1 + 3 files changed, 131 insertions(+), 1 deletion(-) create mode 100644 applications/pctlomalinda/pctlomalinda.py diff --git a/applications/pctlomalinda/pctlomalinda.py b/applications/pctlomalinda/pctlomalinda.py new file mode 100644 index 0000000..fbb679d --- /dev/null +++ b/applications/pctlomalinda/pctlomalinda.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python +import argparse +import uproot +import itk +from itk import PCT as pct + +import numpy as np + + +def build_parser(): + + parser = pct.PCTArgumentParser(description="Convert Loma Linda data to PCT data") + parser.add_argument( + "-i", "--input", help="Root phase space file of particles", required=True + ) + parser.add_argument("-o", "--output", help="Output file name", required=True) + parser.add_argument( + "--plane-in", + help="Plane position of incoming protons", + required=True, + type=float, + ) + parser.add_argument( + "--plane-out", + help="Plane position of outgoing protons", + required=True, + type=float, + ) + parser.add_argument( + "--min-run", help="Minimum run (inclusive)", default=0, type=int + ) + parser.add_argument( + "--max-run", help="Maximum run (exclusive)", default=1e6, type=int + ) + parser.add_argument( + "--verbose", "-v", help="Verbose execution", default=False, action="store_true" + ) + parser.add_argument( + "--ps", help="Name of tree in input phase space", default="PhaseSpace" + ) + + return parser + + +def process(args_info: argparse.Namespace): + + if args_info.verbose: + + def verbose(message): + print(message) + + else: + + def verbose(message): + pass + + tree = uproot.open(args_info.input)[args_info.ps] + pairs = tree.arrays(library="np") + verbose("Read input phase space:\n" + str(pairs)) + + # should match what it in the file, no checks is performed + pairs["u_hit1"] = np.full_like(pairs["u_hit1"], args_info.plane_in) + pairs["u_hit2"] = np.full_like(pairs["u_hit2"], args_info.plane_out) + + verbose("Filter RunIDs…") + # "projection_angle" acts as the run ID here + interception = np.logical_and( + pairs["projection_angle"] < args_info.max_run, + pairs["projection_angle"] >= args_info.min_run, + ) + for k, v in pairs.items(): + pairs[k] = v[interception] + + verbose("Calculating directions…") + dir_in_t = pairs["t_hit1"] - pairs["t_hit0"] + dir_in_v = pairs["v_hit1"] - pairs["v_hit0"] + dir_in_u = pairs["u_hit1"] - pairs["u_hit0"] + norm_in = np.sqrt(dir_in_t**2 + dir_in_v**2 + dir_in_u**2) + dir_out_t = pairs["t_hit3"] - pairs["t_hit2"] + dir_out_v = pairs["v_hit3"] - pairs["v_hit2"] + dir_out_u = pairs["u_hit3"] - pairs["u_hit2"] + norm_out = np.sqrt(dir_out_t**2 + dir_out_v**2 + dir_out_u**2) + + ComponentType = itk.ctype("float") + PixelType = itk.Vector[ComponentType, 3] + ImageType = itk.Image[PixelType, 2] + + runs = np.unique(pairs["projection_angle"]) + number_of_runs = len(runs) + verbose("Identified number of runs: " + str(number_of_runs)) + run_range = range(args_info.min_run, min(number_of_runs, args_info.max_run)) + for r in run_range: + ps_run = pairs["projection_angle"] == runs[r] + len_ps_run = ps_run.sum() + if len_ps_run == 0: + continue + + ps_np = np.empty(shape=(len_ps_run, 5, 3), dtype=np.float32) + ps_np[:, 0, 0] = pairs["t_hit1"][ps_run] + ps_np[:, 0, 1] = pairs["v_hit1"][ps_run] + ps_np[:, 0, 2] = pairs["u_hit1"][ps_run] + ps_np[:, 1, 0] = pairs["t_hit2"][ps_run] + ps_np[:, 1, 1] = pairs["v_hit2"][ps_run] + ps_np[:, 1, 2] = pairs["u_hit2"][ps_run] + ps_np[:, 2, 0] = dir_in_t[ps_run] / norm_in[ps_run] + ps_np[:, 2, 1] = dir_in_v[ps_run] / norm_in[ps_run] + ps_np[:, 2, 2] = dir_in_u[ps_run] / norm_in[ps_run] + ps_np[:, 3, 0] = dir_in_t[ps_run] / norm_out[ps_run] + ps_np[:, 3, 1] = dir_in_v[ps_run] / norm_out[ps_run] + ps_np[:, 3, 2] = dir_in_u[ps_run] / norm_out[ps_run] + ps_np[:, 4, 0] = 0.0 # WEPL is already present in the ROOT file + ps_np[:, 4, 1] = pairs["calculated_WEPL"][ps_run] + ps_np[:, 4, 2] = pairs["timestamp"][ps_run] + + df_itk = itk.GetImageFromArray(ps_np, ttype=ImageType) + + output_file = args_info.output.replace(".", f"{r:04d}.") + itk.imwrite(df_itk, output_file) + verbose(f"Wrote file {output_file}.") + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 086950d..bffa5e4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,7 +45,7 @@ pctbinning = "itk.pctbinning:main" pctfdk = "itk.pctfdk:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" - +pctlomalinda = "itk.pctlomalinda:main" [project.urls] Download = "https://github.com/RTKConsortium/PCT" diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 57e5d85..3466b21 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -25,6 +25,7 @@ "pctfdk", "pctpairprotons", "pctweplfit", + "pctlomalinda", ] # Dynamically access make_application_func from pctExtras From c442c9ea185d8a5376c0bd2a828181a4c79d4c60 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Thu, 9 Apr 2026 16:57:44 +0200 Subject: [PATCH 05/11] ENH: Add pctfilterprotons application --- .../pctfilterprotons/pctfilterprotons.py | 103 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_pct__.py | 1 + 3 files changed, 105 insertions(+) create mode 100644 applications/pctfilterprotons/pctfilterprotons.py diff --git a/applications/pctfilterprotons/pctfilterprotons.py b/applications/pctfilterprotons/pctfilterprotons.py new file mode 100644 index 0000000..8cc0649 --- /dev/null +++ b/applications/pctfilterprotons/pctfilterprotons.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python +import argparse +import itk +from itk import PCT as pct +import numpy as np + + +def build_parser(): + + parser = pct.PCTArgumentParser( + description="Filter protons according to various criteria" + ) + parser.add_argument("-i", "--input", help="Input list-mode file", required=True) + parser.add_argument("-o", "--output", help="Output list-mode file", required=True) + parser.add_argument("--roi-radius", help="Radius of the ROI in mm", type=float) + parser.add_argument( + "--fluence", help="Fluence level as fraction of full fluence", type=float + ) + parser.add_argument("--min-wepl", help="Minimum WEPL", type=float) + parser.add_argument("--max-wepl", help="Maximum WEPL", type=float) + parser.add_argument( + "--verbose", "-v", help="Verbose execution", default=False, action="store_true" + ) + + return parser + + +def process(args_info: argparse.Namespace): + + if args_info.verbose: + + def verbose(message): + print(message) + + else: + + def verbose(message): + pass + + pairs_itk = itk.imread(args_info.input) + pairs = itk.GetArrayFromImage(pairs_itk) + + u_in = np.s_[:, 0, 0] + w_in = np.s_[:, 0, 2] + u_out = np.s_[:, 1, 0] + w_out = np.s_[:, 1, 2] + e_in = np.s_[:, 4, 0] + wepl = np.s_[:, 4, 1] + + if np.any(pairs[e_in] != 0.0) and ( + args_info.min_wepl is not None or args_info.max_wepl is not None + ): + raise ValueError( + "Cannot filter WEPLs because input file does not contain WEPL (e_in != 0)! Aborting." + ) + + interception = np.full_like(pairs[wepl], True) + + if args_info.roi_radius is not None: + # Circular ROI intersection + verbose("Filtering ions outside of the ROI…") + radius = args_info.roi_radius + du = pairs[u_out] - pairs[u_in] + dt = pairs[w_out] - pairs[w_in] + dr_sq = (du * du) + (dt * dt) + d = (pairs[u_in] * pairs[w_out]) - (pairs[u_out] * pairs[w_in]) + delta = (radius * radius * dr_sq) - (d * d) + interception = np.logical_and(interception, delta >= 0.0) + + if args_info.fluence is not None: + # Fluence filtering + verbose("Applying fluence level…") + rng = np.random.default_rng() + fluence_filter = rng.random(interception.shape) + interception = np.logical_and(interception, fluence_filter < args_info.fluence) + + if args_info.min_wepl is not None: + verbose("Filtering low WEPLs…") + interception = np.logical_and(interception, pairs[wepl] >= args_info.min_wepl) + + if args_info.max_wepl is not None: + verbose("Filtering high WEPLs…") + interception = np.logical_and(interception, pairs[wepl] <= args_info.max_wepl) + + verbose("Applying interception filter…") + + ComponentType = itk.ctype("float") + PixelType = itk.Vector[ComponentType, 3] + ImageType = itk.Image[PixelType, 2] + + output_itk = itk.GetImageFromArray(pairs[interception], ttype=ImageType) + itk.imwrite(output_itk, args_info.output) + verbose(f"Wrote file {args_info.output}.") + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index bffa5e4..2e9f903 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,6 +46,7 @@ pctfdk = "itk.pctfdk:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" pctlomalinda = "itk.pctlomalinda:main" +pctfilterprotons = "itk.pctfilterprotons:main" [project.urls] Download = "https://github.com/RTKConsortium/PCT" diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 3466b21..4f6d9df 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -26,6 +26,7 @@ "pctpairprotons", "pctweplfit", "pctlomalinda", + "pctfilterprotons", ] # Dynamically access make_application_func from pctExtras From afe447287d6adbf1c9a7acd4458322d953551a1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Wed, 22 Apr 2026 17:30:08 +0200 Subject: [PATCH 06/11] DOC: Fix PCTInDoxygenGroup test --- include/pctFDKDDConeBeamVarianceReconstructionFilter.h | 4 +++- include/pctFFTVarianceImageFilter.h | 4 +--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h index 5b8d14a..f1ba98a 100644 --- a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h @@ -8,10 +8,12 @@ #include #include -/** \class FDKDDConeBeamReconstructionFilter +/** \class FDKDDConeBeamVarianceReconstructionFilter * TODO * * \author Jannis Dickmann + * + * \ingroup PCT */ namespace pct { diff --git a/include/pctFFTVarianceImageFilter.h b/include/pctFFTVarianceImageFilter.h index dc6918b..cfee56a 100644 --- a/include/pctFFTVarianceImageFilter.h +++ b/include/pctFFTVarianceImageFilter.h @@ -36,11 +36,9 @@ namespace pct * The filter code is based on FFTConvolutionImageFilter by Gaetan Lehmann * (see http://hdl.handle.net/10380/3154) * - * \test rtkvariancefiltertest.cxx - * * \author Simon Rit * - * \ingroup ImageToImageFilter + * \ingroup PCT */ template From bc0ed2bd5f1c3c97d807073b78e1994f06ce8479 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Wed, 22 Apr 2026 18:16:51 +0200 Subject: [PATCH 07/11] ENH: Create small tests for pctlomalinda and pctfilterprotons --- .../pctfilterprotons/pctfilterprotons.py | 3 +- test/pct_application_test.py | 48 +++++++++++++++++++ 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/applications/pctfilterprotons/pctfilterprotons.py b/applications/pctfilterprotons/pctfilterprotons.py index 8cc0649..f2cfed1 100644 --- a/applications/pctfilterprotons/pctfilterprotons.py +++ b/applications/pctfilterprotons/pctfilterprotons.py @@ -18,6 +18,7 @@ def build_parser(): ) parser.add_argument("--min-wepl", help="Minimum WEPL", type=float) parser.add_argument("--max-wepl", help="Maximum WEPL", type=float) + parser.add_argument("--seed", help="Random seed (for reproducibility)", type=int) parser.add_argument( "--verbose", "-v", help="Verbose execution", default=False, action="store_true" ) @@ -70,7 +71,7 @@ def verbose(message): if args_info.fluence is not None: # Fluence filtering verbose("Applying fluence level…") - rng = np.random.default_rng() + rng = np.random.default_rng(args_info.seed) fluence_filter = rng.random(interception.shape) interception = np.logical_and(interception, fluence_filter < args_info.fluence) diff --git a/test/pct_application_test.py b/test/pct_application_test.py index 750ea8f..821e6ca 100644 --- a/test/pct_application_test.py +++ b/test/pct_application_test.py @@ -74,3 +74,51 @@ def test_weplfit_application(tmp_path): ] ) assert np.allclose(eloss_to_wepl_fit, reference) + + +lomalinda_data = download_file_fixture( + "69e21803ed08a1c077afd077", "projection_045.root" +) +baseline_lomalinda_mhd = download_file_fixture( + "69e89639ed08a1c077afd0d9", "baseline_lomalinda0000.mhd" +) +baseline_lomalinda_raw = download_file_fixture( + "69e8963eed08a1c077afd0dc", "baseline_lomalinda0000.raw" +) + + +@pytest.fixture(scope="session") +def test_lomalinda_application( + tmp_path_factory, lomalinda_data, baseline_lomalinda_mhd, baseline_lomalinda_raw +): + output = tmp_path_factory.getbasetemp() / "lomalinda.mhd" + pct.pctlomalinda( + f"-i {lomalinda_data} -o {output} --plane-in -167.2 --plane-out 167.2 --ps recoENTRY -v" + ) + output0000 = str(output).replace(".", "0000.") + test_lomalinda = itk.array_from_image(itk.imread(output0000)) + reference_lomalinda = itk.array_from_image(itk.imread(baseline_lomalinda_mhd)) + assert np.array_equal(test_lomalinda, reference_lomalinda) + return output0000 + + +baseline_filterprotons_mhd = download_file_fixture( + "69e89be4ed08a1c077afd0e5", "baseline_filterprotons.mhd") +baseline_filterprotons_raw = download_file_fixture( + "69e89be5ed08a1c077afd0e8", "baseline_filterprotons.raw") + + +def test_filterprotons_application( + tmp_path, + test_lomalinda_application, + baseline_filterprotons_mhd, + baseline_filterprotons_raw, +): + output = tmp_path / "filterprotons.mhd" + pct.pctfilterprotons( + f"-i {test_lomalinda_application} -o {output} --roi-radius 200 --fluence .5 -v --min-wepl 70 --max-wepl 180 --seed 1234" + ) + test_filterprotons = itk.array_from_image(itk.imread(output)) + reference_filterprotons = itk.array_from_image( + itk.imread(baseline_filterprotons_mhd)) + assert np.array_equal(test_filterprotons, reference_filterprotons) From d1d09daf7ce2ac32e6d4495c2d9766814c189998 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Tue, 19 May 2026 16:25:41 +0200 Subject: [PATCH 08/11] ENH: Move pctfilterprotons features to pctpaircuts --- .../pctfilterprotons/pctfilterprotons.py | 104 ------------------ applications/pctpaircuts/pctpaircuts.cxx | 44 ++++++++ applications/pctpaircuts/pctpaircuts.ggo | 5 + pyproject.toml | 1 - test/Baseline/baseline_paircuts.mhd.sha512 | 1 + test/Baseline/baseline_paircuts.raw.sha512 | 1 + test/CMakeLists.txt | 8 ++ test/Input/input_paircuts.mhd.sha512 | 1 + test/Input/input_paircuts.raw.sha512 | 1 + test/pct_application_test.py | 22 ---- wrapping/__init_pct__.py | 1 - 11 files changed, 61 insertions(+), 128 deletions(-) delete mode 100644 applications/pctfilterprotons/pctfilterprotons.py create mode 100644 test/Baseline/baseline_paircuts.mhd.sha512 create mode 100644 test/Baseline/baseline_paircuts.raw.sha512 create mode 100644 test/Input/input_paircuts.mhd.sha512 create mode 100644 test/Input/input_paircuts.raw.sha512 diff --git a/applications/pctfilterprotons/pctfilterprotons.py b/applications/pctfilterprotons/pctfilterprotons.py deleted file mode 100644 index f2cfed1..0000000 --- a/applications/pctfilterprotons/pctfilterprotons.py +++ /dev/null @@ -1,104 +0,0 @@ -#!/usr/bin/env python -import argparse -import itk -from itk import PCT as pct -import numpy as np - - -def build_parser(): - - parser = pct.PCTArgumentParser( - description="Filter protons according to various criteria" - ) - parser.add_argument("-i", "--input", help="Input list-mode file", required=True) - parser.add_argument("-o", "--output", help="Output list-mode file", required=True) - parser.add_argument("--roi-radius", help="Radius of the ROI in mm", type=float) - parser.add_argument( - "--fluence", help="Fluence level as fraction of full fluence", type=float - ) - parser.add_argument("--min-wepl", help="Minimum WEPL", type=float) - parser.add_argument("--max-wepl", help="Maximum WEPL", type=float) - parser.add_argument("--seed", help="Random seed (for reproducibility)", type=int) - parser.add_argument( - "--verbose", "-v", help="Verbose execution", default=False, action="store_true" - ) - - return parser - - -def process(args_info: argparse.Namespace): - - if args_info.verbose: - - def verbose(message): - print(message) - - else: - - def verbose(message): - pass - - pairs_itk = itk.imread(args_info.input) - pairs = itk.GetArrayFromImage(pairs_itk) - - u_in = np.s_[:, 0, 0] - w_in = np.s_[:, 0, 2] - u_out = np.s_[:, 1, 0] - w_out = np.s_[:, 1, 2] - e_in = np.s_[:, 4, 0] - wepl = np.s_[:, 4, 1] - - if np.any(pairs[e_in] != 0.0) and ( - args_info.min_wepl is not None or args_info.max_wepl is not None - ): - raise ValueError( - "Cannot filter WEPLs because input file does not contain WEPL (e_in != 0)! Aborting." - ) - - interception = np.full_like(pairs[wepl], True) - - if args_info.roi_radius is not None: - # Circular ROI intersection - verbose("Filtering ions outside of the ROI…") - radius = args_info.roi_radius - du = pairs[u_out] - pairs[u_in] - dt = pairs[w_out] - pairs[w_in] - dr_sq = (du * du) + (dt * dt) - d = (pairs[u_in] * pairs[w_out]) - (pairs[u_out] * pairs[w_in]) - delta = (radius * radius * dr_sq) - (d * d) - interception = np.logical_and(interception, delta >= 0.0) - - if args_info.fluence is not None: - # Fluence filtering - verbose("Applying fluence level…") - rng = np.random.default_rng(args_info.seed) - fluence_filter = rng.random(interception.shape) - interception = np.logical_and(interception, fluence_filter < args_info.fluence) - - if args_info.min_wepl is not None: - verbose("Filtering low WEPLs…") - interception = np.logical_and(interception, pairs[wepl] >= args_info.min_wepl) - - if args_info.max_wepl is not None: - verbose("Filtering high WEPLs…") - interception = np.logical_and(interception, pairs[wepl] <= args_info.max_wepl) - - verbose("Applying interception filter…") - - ComponentType = itk.ctype("float") - PixelType = itk.Vector[ComponentType, 3] - ImageType = itk.Image[PixelType, 2] - - output_itk = itk.GetImageFromArray(pairs[interception], ttype=ImageType) - itk.imwrite(output_itk, args_info.output) - verbose(f"Wrote file {args_info.output}.") - - -def main(argv=None): - parser = build_parser() - args_info = parser.parse_args(argv) - process(args_info) - - -if __name__ == "__main__": - main() diff --git a/applications/pctpaircuts/pctpaircuts.cxx b/applications/pctpaircuts/pctpaircuts.cxx index 8cb1554..89d1f4b 100644 --- a/applications/pctpaircuts/pctpaircuts.cxx +++ b/applications/pctpaircuts/pctpaircuts.cxx @@ -8,6 +8,7 @@ #include #include +#include #include #define PAIRS_IN_RAM 1000000 @@ -47,6 +48,14 @@ main(int argc, char * argv[]) std::vector> energies(npixels); std::vector> angles(npixels); + // RNG for fluence + using GeneratorType = itk::Statistics::MersenneTwisterRandomVariateGenerator; + auto generator = GeneratorType::New(); + if (args_info.seed_given) + generator->Initialize(args_info.seed_arg); + else + generator->Initialize(); + // Read pairs using VectorType = itk::Vector; using PairsImageType = itk::Image; @@ -275,6 +284,41 @@ main(int argc, char * argv[]) ++it; } + // Fluence filtering + if (args_info.fluence_given) + { + if (generator->GetVariate() >= args_info.fluence_arg) + continue; + } + + // Check ROI intersection + if (args_info.roiradius_given) + { + float radius = args_info.roiradius_arg; + float dx = pOut[0] - pIn[0]; + float dz = pOut[2] - pIn[2]; + float dr_sq = (dx * dx) + (dz * dz); + float d = (pIn[0] * pOut[2]) - (pOut[0] * pIn[2]); + float delta = (radius * radius * dr_sq) - (d * d); + if (delta < 0.) + continue; + } + + // Filter low or high WEPLs + if (data[0] == 0.) + { + if ((args_info.minwepl_given && (data[1] < args_info.minwepl_arg)) || + (args_info.maxwepl_given && (args_info.maxwepl_arg < data[1]))) + continue; + } + else + { + if (args_info.minwepl_given || args_info.maxwepl_given) + std::cerr << "Cannot specify WEPL bounds because the input file does not contain WEPL data. Aborting." + << std::endl; + return EXIT_FAILURE; + } + static double mag = (args_info.source_arg == 0.) ? 1. : (args_info.source_arg - pOut[2]) / (args_info.source_arg - pIn[2]); diff --git a/applications/pctpaircuts/pctpaircuts.ggo b/applications/pctpaircuts/pctpaircuts.ggo index 2f20103..65e6ee6 100644 --- a/applications/pctpaircuts/pctpaircuts.ggo +++ b/applications/pctpaircuts/pctpaircuts.ggo @@ -11,6 +11,11 @@ option "robust" r "Use robust estimation using 50/19.1 %ile." flag of option "robustopt" - "Use newer options for robust cut." int no default="0" option "primaries" - "Consider only primary protons" flag off option "nonuclear" - "Consider only primary protons without nuclear interactions" flag off +option "roiradius" - "Radius of the ROI in mm" double no +option "minwepl" - "Minimum WEPL" double no +option "maxwepl" - "Maximum WEPL" double no +option "fluence" - "Fluence level as fraction of full fluence" double no +option "seed" - "Random seed (for reproducibility)" int no section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no diff --git a/pyproject.toml b/pyproject.toml index 2e9f903..bffa5e4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,6 @@ pctfdk = "itk.pctfdk:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" pctlomalinda = "itk.pctlomalinda:main" -pctfilterprotons = "itk.pctfilterprotons:main" [project.urls] Download = "https://github.com/RTKConsortium/PCT" diff --git a/test/Baseline/baseline_paircuts.mhd.sha512 b/test/Baseline/baseline_paircuts.mhd.sha512 new file mode 100644 index 0000000..db4ca3e --- /dev/null +++ b/test/Baseline/baseline_paircuts.mhd.sha512 @@ -0,0 +1 @@ +fb6e1fec89f16fb5b85a494e406e0f343e33eef7546720d62c15031a43cc0b1018f6de0446642cc527fce11dd71e931b4fde8df85f21db858c0bba9b93738967 diff --git a/test/Baseline/baseline_paircuts.raw.sha512 b/test/Baseline/baseline_paircuts.raw.sha512 new file mode 100644 index 0000000..4017571 --- /dev/null +++ b/test/Baseline/baseline_paircuts.raw.sha512 @@ -0,0 +1 @@ +49e739d0cc26ae1a40567dbc1aca0d9fd3654bba160704b36e3e0d446e4bb007667629717c029cc18f82ac8eaa35f1d3b40b7bd77e33686ddfa3875e016bc7a8 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6aea818..88d856a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,4 +38,12 @@ if(PCT_BUILD_APPLICATIONS) ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/pctbinning -i DATA{Input/pairs0000.mhd,pairs0000.raw} -o ${ITK_TEST_OUTPUT_DIR}/projections.mhd --source 1000 --size=400,20,220 ) + itk_add_test(NAME pctapppaircutstest + COMMAND itkTestDriver + --compare + DATA{Baseline/baseline_paircuts.mhd,baseline_paircuts.raw} + ${ITK_TEST_OUTPUT_DIR}/paircuts.mhd + ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/pctpaircuts + -i DATA{Input/input_paircuts.mhd,input_paircuts.raw} -o ${ITK_TEST_OUTPUT_DIR}/paircuts.mhd --roiradius 200 --minwepl 70 --maxwepl 180 --fluence .5 --seed 1234 + ) endif() diff --git a/test/Input/input_paircuts.mhd.sha512 b/test/Input/input_paircuts.mhd.sha512 new file mode 100644 index 0000000..85fd2c5 --- /dev/null +++ b/test/Input/input_paircuts.mhd.sha512 @@ -0,0 +1 @@ +37f417e870ab37e307c3c5414fefd7c04e67fe5821ed802ebc573e050ea575d2b76901090baf15d2880cc73877603777cfc7a6c5aa521efa67478d8b1c80796a diff --git a/test/Input/input_paircuts.raw.sha512 b/test/Input/input_paircuts.raw.sha512 new file mode 100644 index 0000000..8ccad1b --- /dev/null +++ b/test/Input/input_paircuts.raw.sha512 @@ -0,0 +1 @@ +bffdb9ff142eec5586509d575dc2c66b83689473bfb36a5c57299f91a5a46c2b032588e62647f528f550c4a13ec311f569582677a084f5c3f333412f25997a41 diff --git a/test/pct_application_test.py b/test/pct_application_test.py index 821e6ca..2e87177 100644 --- a/test/pct_application_test.py +++ b/test/pct_application_test.py @@ -100,25 +100,3 @@ def test_lomalinda_application( reference_lomalinda = itk.array_from_image(itk.imread(baseline_lomalinda_mhd)) assert np.array_equal(test_lomalinda, reference_lomalinda) return output0000 - - -baseline_filterprotons_mhd = download_file_fixture( - "69e89be4ed08a1c077afd0e5", "baseline_filterprotons.mhd") -baseline_filterprotons_raw = download_file_fixture( - "69e89be5ed08a1c077afd0e8", "baseline_filterprotons.raw") - - -def test_filterprotons_application( - tmp_path, - test_lomalinda_application, - baseline_filterprotons_mhd, - baseline_filterprotons_raw, -): - output = tmp_path / "filterprotons.mhd" - pct.pctfilterprotons( - f"-i {test_lomalinda_application} -o {output} --roi-radius 200 --fluence .5 -v --min-wepl 70 --max-wepl 180 --seed 1234" - ) - test_filterprotons = itk.array_from_image(itk.imread(output)) - reference_filterprotons = itk.array_from_image( - itk.imread(baseline_filterprotons_mhd)) - assert np.array_equal(test_filterprotons, reference_filterprotons) diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 4f6d9df..3466b21 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -26,7 +26,6 @@ "pctpairprotons", "pctweplfit", "pctlomalinda", - "pctfilterprotons", ] # Dynamically access make_application_func from pctExtras From 481c0413c858abd57b8b0b65c658b2e927ad1dae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Thu, 21 May 2026 10:03:41 +0200 Subject: [PATCH 09/11] STYLE: Clean debug comments in FFTVarianceImageFilter --- include/pctFFTVarianceImageFilter.hxx | 86 +-------------------------- 1 file changed, 1 insertion(+), 85 deletions(-) diff --git a/include/pctFFTVarianceImageFilter.hxx b/include/pctFFTVarianceImageFilter.hxx index 19be1b1..ac6314e 100644 --- a/include/pctFFTVarianceImageFilter.hxx +++ b/include/pctFFTVarianceImageFilter.hxx @@ -24,27 +24,13 @@ namespace pct template FFTVarianceImageFilter::FFTVarianceImageFilter() -{ - // this->m_HannCutFrequency = 0.; - // this->m_CosineCutFrequency = 0.; - // this->m_HammingFrequency = 0.; - // this->m_HannCutFrequency = 0.; - // this->m_HannCutFrequencyY = 0.; - // this->m_RamLakCutFrequency = 0.; - // this->m_SheppLoganCutFrequency = 0.; -} +{} template void FFTVarianceImageFilter::UpdateFFTProjectionsConvolutionKernel( const SizeType s) { - // if(this->m_KernelFFT.GetPointer() != ITK_NULLPTR && s == this->m_PreviousKernelUpdateSize) - // { - // return; - // } - // this->m_PreviousKernelUpdateSize = s; - const int width = s[0]; const int height = s[1]; @@ -72,19 +58,6 @@ FFTVarianceImageFilter::UpdateFFTProje kernel->SetPixel(jx, v); } - // ///////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // typedef itk::ImageRegionIteratorWithIndex InverseIteratorType2; - // InverseIteratorType2 itiK2(kernel, kernel->GetLargestPossibleRegion()); - // std::ofstream out0("kernel_real.csv"); - // itiK2.GoToBegin(); - // for(; !itiK2.IsAtEnd(); ++itiK2) - // { - // out0 << itiK2.Get() << std::endl; - // } - // out0.close(); - // // DEBUG - // FFT kernel typedef itk::RealToHalfHermitianForwardFFTImageFilter FFTType; typename FFTType::Pointer fftK = FFTType::New(); @@ -157,17 +130,6 @@ FFTVarianceImageFilter::UpdateFFTProje itK.Set(itK.Get() * TFFTPrecision(0.)); } - // ///////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // std::ofstream out1("kernel_fourier.csv"); - // itK.GoToBegin(); - // for(; !itK.IsAtEnd(); ++itK) - // { - // out1 << itK.Get().real() << " " << itK.Get().imag() << std::endl; - // } - // out1.close(); - // // DEBUG - // iFFT kernel typedef itk::HalfHermitianToRealInverseFFTImageFilter IFFTType; typename IFFTType::Pointer ifftK = IFFTType::New(); @@ -206,17 +168,6 @@ FFTVarianceImageFilter::UpdateFFTProje } finterp *= aprecision; - // ///////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // std::ofstream out9("kernel_window_real.csv"); - // itiK.GoToBegin(); - // for(; !itiK.IsAtEnd(); ++itiK) - // { - // out9 << itiK.Get() << std::endl; - // } - // out9.close(); - // // DEBUG - // square kernel and multiply with finterp itiK.GoToBegin(); for (; !itiK.IsAtEnd(); ++itiK) @@ -224,17 +175,6 @@ FFTVarianceImageFilter::UpdateFFTProje itiK.Set(itiK.Get() * itiK.Get() * finterp); } - // ///////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // std::ofstream out2("kernel_sq_real.csv"); - // itiK.GoToBegin(); - // for(; !itiK.IsAtEnd(); ++itiK) - // { - // out2 << itiK.Get() << std::endl; - // } - // out2.close(); - // // DEBUG - // FFT kernel typename FFTType::Pointer fftK2 = FFTType::New(); fftK2->SetInput(ifftK->GetOutput()); @@ -242,18 +182,6 @@ FFTVarianceImageFilter::UpdateFFTProje fftK2->Update(); this->m_KernelFFT = fftK2->GetOutput(); - // ///////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // std::ofstream out3("kernel_sq_fourier.csv"); - // IteratorType itK2(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion() ); - // itK2.GoToBegin(); - // for(; !itK2.IsAtEnd(); ++itK2) - // { - // out3 << itK2.Get().real() << " " << itK2.Get().imag() << std::endl; - // } - // out3.close(); - // // DEBUG - // Replicate and window if required if (this->GetHannCutFrequencyY() > 0.) { @@ -293,18 +221,6 @@ FFTVarianceImageFilter::UpdateFFTProje } this->m_KernelFFT->DisconnectPipeline(); - - // ////////////////////////////////////////////////////////////////////////////////////////////// - // // DEBUG - // std::ofstream out4("finterp.csv"); - // out4 << finterp << std::endl; - // out4 << ratiogcg2 << std::endl; - // out4.close(); - // std::cerr << std::endl << "finterp = " << finterp << std::endl; - // std::cerr << "ratiogcg2 = " << ratiogcg2 << std::endl; - // std::cerr << "Aborting here since filters were alrady written out." << std::endl; - // exit(9); - // // DEBUG } } // namespace pct From d08b91b3b913a5d18b669547d029ebc03fb80038 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 12 Jun 2026 10:24:30 +0200 Subject: [PATCH 10/11] DOC: Add short description to FDKDDConeBeamVarianceReconstructionFilter --- include/pctFDKDDConeBeamVarianceReconstructionFilter.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h index f1ba98a..48d8c42 100644 --- a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h @@ -9,7 +9,8 @@ #include /** \class FDKDDConeBeamVarianceReconstructionFilter - * TODO + * Implements variance reconstruction with distance-driven filtered backprojection reconstruction + * (doi:10.1118/1.4789589) from variance projections. * * \author Jannis Dickmann * From dae94d324da70a561c9a97a6c6372f9c840fd512 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 12 Jun 2026 10:47:43 +0200 Subject: [PATCH 11/11] PERF: Use =default constructor in FFTVarianceImageFilter --- include/pctFFTVarianceImageFilter.h | 2 +- include/pctFFTVarianceImageFilter.hxx | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/include/pctFFTVarianceImageFilter.h b/include/pctFFTVarianceImageFilter.h index cfee56a..86b9d51 100644 --- a/include/pctFFTVarianceImageFilter.h +++ b/include/pctFFTVarianceImageFilter.h @@ -72,7 +72,7 @@ class ITK_EXPORT FFTVarianceImageFilter : public rtk::FFTRampImageFilter -FFTVarianceImageFilter::FFTVarianceImageFilter() -{} - template void FFTVarianceImageFilter::UpdateFFTProjectionsConvolutionKernel(