From eb597bfe2bd6b34c5d214eb505953e85a373ad4c Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Mon, 13 Oct 2025 16:32:52 +0200 Subject: [PATCH 01/12] bash: line 1: q: command not found --- applications/pctfdk/pctfdk.py | 162 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_pct__.py | 1 + wrapping/pctDDParkerShortScanImageFilter.wrap | 4 +- ...TProjectionsConvolutionImageFilterPCT.wrap | 5 + wrapping/rtkFFTRampImageFilterPCT.wrap | 5 + 6 files changed, 177 insertions(+), 1 deletion(-) create mode 100644 applications/pctfdk/pctfdk.py create mode 100644 wrapping/rtkFFTProjectionsConvolutionImageFilterPCT.wrap create mode 100644 wrapping/rtkFFTRampImageFilterPCT.wrap diff --git a/applications/pctfdk/pctfdk.py b/applications/pctfdk/pctfdk.py new file mode 100644 index 0000000..b6e6b15 --- /dev/null +++ b/applications/pctfdk/pctfdk.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python +import argparse +import numpy as np +import itk +from itk import PCT as pct + + +def build_parser(): + parser = pct.PCTArgumentParser( + description="Reconstruct a 3D volume from a sequence of projections [Feldkamp, David, Kress, 1984]." + ) + # General + parser.add_argument( + "--verbose", "-v", help="Verbose execution", action="store_true" + ) + parser.add_argument( + "--geometry", "-g", help="XML geometry file name", type=str, required=True + ) + parser.add_argument( + "--path", "-p", help="Path containing projections", type=str, required=True + ) + parser.add_argument( + "--regexp", + "-r", + help="Regular expression to select projection files in path", + type=str, + required=True, + ) + parser.add_argument( + "--output", "-o", help="Output file name", type=str, required=True + ) + parser.add_argument( + "--lowmem", + "-l", + help="Load only one projection per thread in memory", + action="store_true", + ) + parser.add_argument( + "--wpc", + help="Water precorrection coefficients (default is no correction)", + type=float, + nargs="+", + ) + + # Ramp filter + parser.add_argument( + "--pad", + help="Data padding parameter to correct for truncation", + type=float, + default=0.0, + ) + parser.add_argument( + "--hann", + help="Cut frequency for hann window in ]0,1] (0.0 disables it)", + type=float, + default=0.0, + ) + parser.add_argument( + "--hannY", + help="Cut frequency for hann window in ]0,1] (0.0 disables it)", + type=float, + default=0.0, + ) + + # Volume properties + parser.add_argument( + "--origin", help="Origin (default=centered)", type=float, nargs="+" + ) + parser.add_argument( + "--dimension", help="Deprecated. Use --size instead.", type=int, nargs="+" + ) + parser.add_argument("--size", help="Size", type=int, nargs="+", default=[256]) + parser.add_argument( + "--spacing", help="Spacing", type=float, nargs="+", default=[1.0] + ) + parser.add_argument( + "--direction", help="Direction (3x3 row-major)", type=float, nargs="+" + ) + parser.add_argument( + "--like", + help="Copy info from image (origin, size, spacing, direction)", + type=str, + ) + return parser + + +def process(args_info: argparse.Namespace): + from itk import RTK as rtk + + # Generate file names + names = itk.RegularExpressionSeriesFileNames.New() + names.SetDirectory(args_info.path) + names.SetNumericSort(False) + names.SetRegularExpression(args_info.regexp) + names.SetSubMatch(0) + if args_info.verbose: + print(f"Regular expression matches {len(names.GetFileNames())} file(s)...") + + # Projections reader + OutputPixelType = itk.F + ProjectionImageType = itk.Image[OutputPixelType, 4] + ReaderType = rtk.ProjectionsReader[ProjectionImageType] + reader = ReaderType.New() + reader.SetFileNames(names.GetFileNames()) + if args_info.wpc: + reader.SetWaterPrecorrectionCoefficients([float(c) for c in args_info.wpc]) + + # Geometry + if args_info.verbose: + print(f"Reading geometry information from {args_info.geometry}...") + geometryReader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New() + geometryReader.SetFilename(args_info.geometry) + geometryReader.GenerateOutputInformation() + + # Short scan image filter + pssf = pct.DDParkerShortScanImageFilter[ProjectionImageType].New() + pssf.SetInput(reader.GetOutput()) + pssf.SetGeometry(geometryReader.GetOutputObject()) + pssf.InPlaceOff() + + # Create reconstructed image + OutputImageType = itk.Image[OutputPixelType, 3] + constantImageSource = rtk.ConstantImageSource[OutputImageType].New() + rtk.SetConstantImageSourceFromArgParse(constantImageSource, args_info) + + # FDK reconstruction filtering + FeldkampType = pct.FDKDDConeBeamReconstructionFilter[OutputImageType] + feldkamp = FeldkampType.New() + feldkamp.SetInput(0, constantImageSource.GetOutput()) + feldkamp.SetProjectionStack(pssf.GetOutput()) + feldkamp.SetGeometry(geometryReader.GetOutputObject()) + feldkamp.GetRampFilter().SetTruncationCorrection(args_info.pad) + feldkamp.GetRampFilter().SetHannCutFrequency(args_info.hann) + feldkamp.GetRampFilter().SetHannCutFrequencyY(args_info.hannY) + + # Write + WriterType = itk.ImageFileWriter[OutputImageType] + writer = WriterType.New() + writer.SetFileName(args_info.output) + writer.SetInput(feldkamp.GetOutput()) + + if args_info.verbose: + print("Reconstructing and writing... ", end="", flush=True) + + writerProbe = itk.TimeProbe() + writerProbe.Start() + writer.Update() + writerProbe.Stop() + + if args_info.verbose: + print(f"It took {writerProbe.GetMean()} {writerProbe.GetUnit()}") + feldkamp.PrintTiming(itk.ostream) + + +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 8e46967..01b2e48 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ ] [project.scripts] +pctfdk = "itk.pctfdk:main" pctfdktwodweights = "itk.pctfdktwodweights:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 19b90d5..514320f 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -18,6 +18,7 @@ # Application modules _app_modules = [ + "pctfdk", "pctfdktwodweights", "pctpairprotons", "pctweplfit", diff --git a/wrapping/pctDDParkerShortScanImageFilter.wrap b/wrapping/pctDDParkerShortScanImageFilter.wrap index 7b69f49..febee28 100644 --- a/wrapping/pctDDParkerShortScanImageFilter.wrap +++ b/wrapping/pctDDParkerShortScanImageFilter.wrap @@ -1,3 +1,5 @@ itk_wrap_class("pct::DDParkerShortScanImageFilter" POINTER) - itk_wrap_template("IF4IF4" "itk::Image<${ITKT_F}, 4>, itk::Image<${ITKT_F}, 4>") + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("I${ITKM_${t}}4" "itk::Image<${ITKT_${t}}, 4>") + endforeach() itk_end_wrap_class() diff --git a/wrapping/rtkFFTProjectionsConvolutionImageFilterPCT.wrap b/wrapping/rtkFFTProjectionsConvolutionImageFilterPCT.wrap new file mode 100644 index 0000000..b9a3d12 --- /dev/null +++ b/wrapping/rtkFFTProjectionsConvolutionImageFilterPCT.wrap @@ -0,0 +1,5 @@ +itk_wrap_class("rtk::FFTProjectionsConvolutionImageFilter" POINTER) + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("I${ITKM_${t}}4I${ITKM_${t}}4${ITKM_D}" "itk::Image<${ITKT_${t}}, 4>, itk::Image<${ITKT_${t}}, 4>, ${ITKT_D}") + endforeach() +itk_end_wrap_class() diff --git a/wrapping/rtkFFTRampImageFilterPCT.wrap b/wrapping/rtkFFTRampImageFilterPCT.wrap new file mode 100644 index 0000000..230675a --- /dev/null +++ b/wrapping/rtkFFTRampImageFilterPCT.wrap @@ -0,0 +1,5 @@ +itk_wrap_class("rtk::FFTRampImageFilter" POINTER) + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("I${ITKM_${t}}4I${ITKM_${t}}4${ITKM_D}" "itk::Image<${ITKT_${t}}, 4>, itk::Image<${ITKT_${t}}, 4>, ${ITKT_D}") + endforeach() +itk_end_wrap_class() From b8987d8db5cf1afbb5433a009fbd9f9db448ad0a Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 14 Oct 2025 16:03:25 +0200 Subject: [PATCH 02/12] BUG: Use Get/SetObjectMacro instead of Get/SetMacro for pointers --- include/pctDDParkerShortScanImageFilter.h | 4 +-- include/pctFDKDDBackProjectionImageFilter.h | 4 +-- .../pctFDKDDConeBeamReconstructionFilter.h | 13 +++++----- .../pctFDKDDConeBeamReconstructionFilter.hxx | 25 +++---------------- include/pctFDKDDWeightProjectionFilter.h | 4 +-- include/pctFDKDDWeightProjectionFilter.hxx | 2 +- include/pctProtonPairsToBackProjection.h | 16 ++++++------ ...pctProtonPairsToDistanceDrivenProjection.h | 14 +++++------ 8 files changed, 32 insertions(+), 50 deletions(-) diff --git a/include/pctDDParkerShortScanImageFilter.h b/include/pctDDParkerShortScanImageFilter.h index a5e0f8a..696aa84 100644 --- a/include/pctDDParkerShortScanImageFilter.h +++ b/include/pctDDParkerShortScanImageFilter.h @@ -47,8 +47,8 @@ class ITK_TEMPLATE_EXPORT DDParkerShortScanImageFilter : public itk::InPlaceImag itkOverrideGetNameOfClassMacro(DDParkerShortScanImageFilter); /** Get / Set the object pointer to projection geometry */ - itkGetMacro(Geometry, GeometryPointer); - itkSetMacro(Geometry, GeometryPointer); + itkGetConstObjectMacro(Geometry, GeometryType); + itkSetObjectMacro(Geometry, GeometryType); protected: DDParkerShortScanImageFilter() { this->SetInPlace(true); } diff --git a/include/pctFDKDDBackProjectionImageFilter.h b/include/pctFDKDDBackProjectionImageFilter.h index f37a9fc..2d1c4e1 100644 --- a/include/pctFDKDDBackProjectionImageFilter.h +++ b/include/pctFDKDDBackProjectionImageFilter.h @@ -36,8 +36,8 @@ class ITK_TEMPLATE_EXPORT FDKDDBackProjectionImageFilter GetDDProjection(const unsigned int iProj); /** Get / Set the stack of projection images */ - itkGetMacro(ProjectionStack, ProjectionStackPointer); - itkSetMacro(ProjectionStack, ProjectionStackPointer); + itkGetModifiableObjectMacro(ProjectionStack, ProjectionStackType); + itkSetObjectMacro(ProjectionStack, ProjectionStackType); /** Creates the #iProj index to index projection matrix with current inputs instead of the physical point to physical point projection matrix provided by Geometry */ diff --git a/include/pctFDKDDConeBeamReconstructionFilter.h b/include/pctFDKDDConeBeamReconstructionFilter.h index e9699a2..6008fc0 100644 --- a/include/pctFDKDDConeBeamReconstructionFilter.h +++ b/include/pctFDKDDConeBeamReconstructionFilter.h @@ -46,10 +46,8 @@ class ITK_TEMPLATE_EXPORT FDKDDConeBeamReconstructionFilter : public itk::ImageT itkOverrideGetNameOfClassMacro(FDKDDConeBeamReconstructionFilter); /** Get / Set the object pointer to projection geometry */ - virtual rtk::ThreeDCircularProjectionGeometry::Pointer - GetGeometry(); - virtual void - SetGeometry(const rtk::ThreeDCircularProjectionGeometry::Pointer _arg); + itkGetConstObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); + itkSetObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); /** Get pointer to the ramp filter used by the feldkamp reconstruction */ typename RampFilterType::Pointer @@ -62,8 +60,8 @@ class ITK_TEMPLATE_EXPORT FDKDDConeBeamReconstructionFilter : public itk::ImageT PrintTiming(std::ostream & os) const; /** Get / Set the stack of projection images */ - itkGetMacro(ProjectionStack, ProjectionStackPointer); - itkSetMacro(ProjectionStack, ProjectionStackPointer); + itkGetModifiableObjectMacro(ProjectionStack, ProjectionStackType); + itkSetObjectMacro(ProjectionStack, ProjectionStackType); protected: FDKDDConeBeamReconstructionFilter(); @@ -101,7 +99,8 @@ class ITK_TEMPLATE_EXPORT FDKDDConeBeamReconstructionFilter : public itk::ImageT itk::TimeProbe m_FilterProbe; itk::TimeProbe m_BackProjectionProbe; - ProjectionStackPointer m_ProjectionStack; + ProjectionStackPointer m_ProjectionStack; + rtk::ThreeDCircularProjectionGeometry::Pointer m_Geometry; }; // end of class } // end namespace pct diff --git a/include/pctFDKDDConeBeamReconstructionFilter.hxx b/include/pctFDKDDConeBeamReconstructionFilter.hxx index 489a29b..500f866 100644 --- a/include/pctFDKDDConeBeamReconstructionFilter.hxx +++ b/include/pctFDKDDConeBeamReconstructionFilter.hxx @@ -45,6 +45,10 @@ FDKDDConeBeamReconstructionFilter::Gen { const unsigned int Dimension = this->InputImageDimension; + // Propagate geometry to subfilters + m_WeightFilter->SetGeometry(m_Geometry); + m_BackProjectionFilter->SetGeometry(m_Geometry); + // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function typename ExtractFilterType::InputImageRegionType projRegion; @@ -110,27 +114,6 @@ FDKDDConeBeamReconstructionFilter::Gen Superclass::GraftOutput(m_BackProjectionFilter->GetOutput()); } -template -rtk::ThreeDCircularProjectionGeometry::Pointer -FDKDDConeBeamReconstructionFilter::GetGeometry() -{ - return this->m_WeightFilter->GetGeometry(); -} - -template -void -FDKDDConeBeamReconstructionFilter::SetGeometry( - const rtk::ThreeDCircularProjectionGeometry::Pointer _arg) -{ - itkDebugMacro("setting GeometryPointer to " << _arg); - if (this->GetGeometry() != _arg) - { - m_WeightFilter->SetGeometry(_arg); - m_BackProjectionFilter->SetGeometry(_arg.GetPointer()); - this->Modified(); - } -} - template void FDKDDConeBeamReconstructionFilter::PrintTiming(std::ostream & os) const diff --git a/include/pctFDKDDWeightProjectionFilter.h b/include/pctFDKDDWeightProjectionFilter.h index cc0ffb5..5d4ad2d 100644 --- a/include/pctFDKDDWeightProjectionFilter.h +++ b/include/pctFDKDDWeightProjectionFilter.h @@ -43,8 +43,8 @@ class ITK_TEMPLATE_EXPORT FDKDDWeightProjectionFilter : public itk::InPlaceImage itkOverrideGetNameOfClassMacro(FDKDDWeightProjectionFilter); /** Get/ Set geometry structure */ - itkGetMacro(Geometry, rtk::ThreeDCircularProjectionGeometry::Pointer); - itkSetMacro(Geometry, rtk::ThreeDCircularProjectionGeometry::Pointer); + itkGetConstObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); + itkSetObjectMacro(Geometry, rtk::ThreeDCircularProjectionGeometry); protected: FDKDDWeightProjectionFilter() {} diff --git a/include/pctFDKDDWeightProjectionFilter.hxx b/include/pctFDKDDWeightProjectionFilter.hxx index 58a7da3..e609814 100644 --- a/include/pctFDKDDWeightProjectionFilter.hxx +++ b/include/pctFDKDDWeightProjectionFilter.hxx @@ -12,7 +12,7 @@ void FDKDDWeightProjectionFilter::BeforeThreadedGenerateData() { // Get angular weights from geometry - m_AngularWeightsAndRampFactor = this->GetGeometry()->GetAngularGaps(m_Geometry->GetGantryAngles()); + m_AngularWeightsAndRampFactor = m_Geometry->GetAngularGaps(m_Geometry->GetGantryAngles()); for (unsigned int k = 0; k < m_AngularWeightsAndRampFactor.size(); k++) { diff --git a/include/pctProtonPairsToBackProjection.h b/include/pctProtonPairsToBackProjection.h index b5bf403..e5bda73 100644 --- a/include/pctProtonPairsToBackProjection.h +++ b/include/pctProtonPairsToBackProjection.h @@ -71,22 +71,22 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToBackProjection : public itk::InPlaceImage itkSetMacro(MostLikelyPathPolynomialDegree, int); /** Get/Set the boundaries of the object. */ - itkGetMacro(QuadricIn, RQIType::Pointer); - itkSetMacro(QuadricIn, RQIType::Pointer); - itkGetMacro(QuadricOut, RQIType::Pointer); - itkSetMacro(QuadricOut, RQIType::Pointer); + itkGetConstObjectMacro(QuadricIn, RQIType); + itkSetObjectMacro(QuadricIn, RQIType); + itkGetConstObjectMacro(QuadricOut, RQIType); + itkSetObjectMacro(QuadricOut, RQIType); /** Get/Set the count of proton pairs per pixel. */ - itkGetMacro(Counts, CountImagePointer); - itkSetMacro(Counts, CountImagePointer); + itkGetModifiableObjectMacro(Counts, CountImageType); + itkSetObjectMacro(Counts, CountImageType); /** Get/Set the ionization potential used in the Bethe-Bloch equation. */ itkGetMacro(IonizationPotential, double); itkSetMacro(IonizationPotential, double); /** Get / Set the object pointer to projection geometry */ - itkGetMacro(Geometry, GeometryPointer); - itkSetMacro(Geometry, GeometryPointer); + itkGetConstObjectMacro(Geometry, GeometryType); + itkSetObjectMacro(Geometry, GeometryType); /** Get / Set the boolean desactivating rotation to bin in coordinate orientation. Default is off. */ itkGetMacro(DisableRotation, bool); diff --git a/include/pctProtonPairsToDistanceDrivenProjection.h b/include/pctProtonPairsToDistanceDrivenProjection.h index af7fe00..ec02726 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.h +++ b/include/pctProtonPairsToDistanceDrivenProjection.h @@ -70,19 +70,19 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection itkSetMacro(MaterialBudget, double); /** Get/Set the boundaries of the object. */ - itkGetMacro(QuadricIn, RQIType::Pointer); - itkSetMacro(QuadricIn, RQIType::Pointer); - itkGetMacro(QuadricOut, RQIType::Pointer); - itkSetMacro(QuadricOut, RQIType::Pointer); + itkGetConstObjectMacro(QuadricIn, RQIType); + itkSetObjectMacro(QuadricIn, RQIType); + itkGetConstObjectMacro(QuadricOut, RQIType); + itkSetObjectMacro(QuadricOut, RQIType); /** Get/Set the count of proton pairs per pixel. */ - itkGetMacro(Count, CountImagePointer); + itkGetModifiableObjectMacro(Count, CountImageType); /** Get/Set the angle of proton pairs per pixel. */ - itkGetMacro(Angle, AngleImagePointer); + itkGetModifiableObjectMacro(Angle, AngleImageType); /** Get/Set the angle of proton pairs per pixel. */ - itkGetMacro(SquaredOutput, OutputImagePointer); + itkGetModifiableObjectMacro(SquaredOutput, OutputImageType); /** Get/Set the ionization potential used in the Bethe-Bloch equation. */ itkGetMacro(IonizationPotential, double); From 893697b355bb8536ba25795ac98c45cced2f0751 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 15 Oct 2025 16:35:50 +0200 Subject: [PATCH 03/12] ENH: Refactor imports in pctweplfit.py and pctpairprotons.py --- applications/pctpairprotons/pctpairprotons.py | 4 ++-- applications/pctweplfit/pctweplfit.py | 20 ++++++++----------- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/applications/pctpairprotons/pctpairprotons.py b/applications/pctpairprotons/pctpairprotons.py index 0bb42c1..147f401 100644 --- a/applications/pctpairprotons/pctpairprotons.py +++ b/applications/pctpairprotons/pctpairprotons.py @@ -1,10 +1,8 @@ #!/usr/bin/env python import argparse import json -import uproot import itk from itk import PCT as pct - import numpy as np import numpy.lib.recfunctions as rfn @@ -78,6 +76,8 @@ def build_parser(): def process(args_info: argparse.Namespace): + import uproot + if args_info.verbose: def verbose(message): diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index 8ed9dd9..f5b77e9 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -4,23 +4,12 @@ import sys from multiprocessing import Pool, Manager -import matplotlib.pyplot as plt -import opengate as gate -import uproot import numpy as np import itk from itk import PCT as pct epsilon_mm = 1e-5 -# Units -nm = gate.g4_units.nm -mm = gate.g4_units.mm -cm = gate.g4_units.cm -m = gate.g4_units.m -sec = gate.g4_units.second -MeV = gate.g4_units.MeV - def pv(verbose, *args, **kwargs): if verbose: @@ -39,6 +28,10 @@ def tof_fit_mc( visu, verbose, ): + import opengate as gate + + u = gate.g4_units + nm, mm, cm, m, sec, MeV = u.nm, u.mm, u.cm, u.m, u.second, u.MeV pv(verbose, "Starting simulation with following parameters: " + str(locals())) @@ -96,7 +89,6 @@ def tof_fit_mc( sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" # Phase spaces - def add_detector(name, translation, attach_to_phantom=False): plane = sim.add_volume("Box", "PlanePhaseSpace" + name) plane.size = [phantom_width_cm * cm, phantom_width_cm * cm, 1 * nm] @@ -149,6 +141,8 @@ def process_phantom_length( elosses, verbose, ): + import uproot + tofs_phantom_length = [] wepls_phantom_length = [] elosses_phantom_length = [] @@ -304,6 +298,8 @@ def fit(xs, ys, xlabel, ylabel): f.write(f"{xpoint} {ypoint}\n") if display or savefig: + import matplotlib.pyplot as plt + plt.figure() plt.plot(xmedians, ymedians, "+", label="Medians") From 17d9199580b8b7cf455c1361f581bbd7a1fe68cb Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 15 Oct 2025 21:51:13 +0100 Subject: [PATCH 04/12] STYLE: Correct wrapping file name of pct::ZengWeightedBackProjectionImageFilter --- ...eFilter.wrap => pctZengWeightedBackProjectionImageFilter.wrap} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename wrapping/{pctZengBackProjectionImageFilter.wrap => pctZengWeightedBackProjectionImageFilter.wrap} (100%) diff --git a/wrapping/pctZengBackProjectionImageFilter.wrap b/wrapping/pctZengWeightedBackProjectionImageFilter.wrap similarity index 100% rename from wrapping/pctZengBackProjectionImageFilter.wrap rename to wrapping/pctZengWeightedBackProjectionImageFilter.wrap From 8e343b7c681bfcd1962e7663948cbb8663068238 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 15 Oct 2025 22:05:53 +0100 Subject: [PATCH 05/12] ENH: Remove output stream parameter from PrintTiming and use std::cout --- applications/pctfdk/pctfdk.cxx | 2 +- applications/pctfdk/pctfdk.py | 2 +- include/pctEnergyAdaptiveMLPFunction.h | 2 +- include/pctFDKDDConeBeamReconstructionFilter.h | 2 +- include/pctFDKDDConeBeamReconstructionFilter.hxx | 12 +++++++----- include/pctMostLikelyPathFunction.h | 2 +- include/pctPolynomialMLPFunction.h | 2 +- include/pctProtonPairsToBackProjection.hxx | 2 +- include/pctProtonPairsToDistanceDrivenProjection.hxx | 2 +- include/pctSchulteMLPFunction.h | 2 +- src/pctEnergyAdaptiveMLPFunction.cxx | 6 +++--- src/pctPolynomialMLPFunction.cxx | 6 +++--- src/pctSchulteMLPFunction.cxx | 8 ++++---- 13 files changed, 26 insertions(+), 24 deletions(-) diff --git a/applications/pctfdk/pctfdk.cxx b/applications/pctfdk/pctfdk.cxx index 6ce02bb..bd4b2e7 100644 --- a/applications/pctfdk/pctfdk.cxx +++ b/applications/pctfdk/pctfdk.cxx @@ -86,7 +86,7 @@ main(int argc, char * argv[]) if (args_info.verbose_flag) { std::cout << "It took " << writerProbe.GetMean() << ' ' << writerProbe.GetUnit() << std::endl; - feldkamp->PrintTiming(std::cout); + feldkamp->PrintTiming(); } return EXIT_SUCCESS; diff --git a/applications/pctfdk/pctfdk.py b/applications/pctfdk/pctfdk.py index b6e6b15..181b73d 100644 --- a/applications/pctfdk/pctfdk.py +++ b/applications/pctfdk/pctfdk.py @@ -149,7 +149,7 @@ def process(args_info: argparse.Namespace): if args_info.verbose: print(f"It took {writerProbe.GetMean()} {writerProbe.GetUnit()}") - feldkamp.PrintTiming(itk.ostream) + feldkamp.PrintTiming() def main(argv=None): diff --git a/include/pctEnergyAdaptiveMLPFunction.h b/include/pctEnergyAdaptiveMLPFunction.h index 262158f..173ac7c 100644 --- a/include/pctEnergyAdaptiveMLPFunction.h +++ b/include/pctEnergyAdaptiveMLPFunction.h @@ -135,7 +135,7 @@ class PCT_EXPORT EnergyAdaptiveMLPFunction : public MostLikelyPathFunction::Gen template void -FDKDDConeBeamReconstructionFilter::PrintTiming(std::ostream & os) const +FDKDDConeBeamReconstructionFilter::PrintTiming() const { - os << "FDKDDConeBeamReconstructionFilter timing:" << std::endl; - os << " Prefilter operations: " << m_PreFilterProbe.GetTotal() << ' ' << m_PreFilterProbe.GetUnit() << std::endl; - os << " Ramp filter: " << m_FilterProbe.GetTotal() << ' ' << m_FilterProbe.GetUnit() << std::endl; - os << " Backprojection: " << m_BackProjectionProbe.GetTotal() << ' ' << m_BackProjectionProbe.GetUnit() << std::endl; + std::cout << "FDKDDConeBeamReconstructionFilter timing:" << std::endl; + std::cout << " Prefilter operations: " << m_PreFilterProbe.GetTotal() << ' ' << m_PreFilterProbe.GetUnit() + << std::endl; + std::cout << " Ramp filter: " << m_FilterProbe.GetTotal() << ' ' << m_FilterProbe.GetUnit() << std::endl; + std::cout << " Backprojection: " << m_BackProjectionProbe.GetTotal() << ' ' << m_BackProjectionProbe.GetUnit() + << std::endl; } } // end namespace pct diff --git a/include/pctMostLikelyPathFunction.h b/include/pctMostLikelyPathFunction.h index 3b35532..7d56e15 100644 --- a/include/pctMostLikelyPathFunction.h +++ b/include/pctMostLikelyPathFunction.h @@ -85,7 +85,7 @@ class ITK_TEMPLATE_EXPORT MostLikelyPathFunction : public itk::LightObject #ifdef MLP_TIMING /** Print timing information */ virtual void - PrintTiming(std::ostream & os) + PrintTiming() {} #endif diff --git a/include/pctPolynomialMLPFunction.h b/include/pctPolynomialMLPFunction.h index 4fa5524..8b554f1 100644 --- a/include/pctPolynomialMLPFunction.h +++ b/include/pctPolynomialMLPFunction.h @@ -169,7 +169,7 @@ class PCT_EXPORT PolynomialMLPFunction : public MostLikelyPathFunction #ifdef MLP_TIMING /** Print timing information */ virtual void - PrintTiming(std::ostream & os) override; + PrintTiming() override; #endif protected: diff --git a/include/pctProtonPairsToBackProjection.hxx b/include/pctProtonPairsToBackProjection.hxx index 2628d33..b42812d 100644 --- a/include/pctProtonPairsToBackProjection.hxx +++ b/include/pctProtonPairsToBackProjection.hxx @@ -307,7 +307,7 @@ ProtonPairsToBackProjection::GenerateData() << m_ProtonPairs->GetLargestPossibleRegion().GetSize(1) << " pairs of protons processed (100%) in all threads." << std::endl; #ifdef MLP_TIMING - mlp->PrintTiming(std::cout); + mlp->PrintTiming(); #endif } this->AfterThreadedGenerateData(); diff --git a/include/pctProtonPairsToDistanceDrivenProjection.hxx b/include/pctProtonPairsToDistanceDrivenProjection.hxx index 37db9f6..879b2d9 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.hxx +++ b/include/pctProtonPairsToDistanceDrivenProjection.hxx @@ -457,7 +457,7 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera { std::cout << '\r' << region.GetSize(1) << " pairs of protons processed (100%) in thread 1" << std::endl; #ifdef MLP_TIMING - mlp->PrintTiming(std::cout); + mlp->PrintTiming(); #endif } } diff --git a/include/pctSchulteMLPFunction.h b/include/pctSchulteMLPFunction.h index f46a873..d5d153e 100644 --- a/include/pctSchulteMLPFunction.h +++ b/include/pctSchulteMLPFunction.h @@ -179,7 +179,7 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction #ifdef MLP_TIMING /** Print timing information */ virtual void - PrintTiming(std::ostream & os) override; + PrintTiming() override; #endif protected: diff --git a/src/pctEnergyAdaptiveMLPFunction.cxx b/src/pctEnergyAdaptiveMLPFunction.cxx index 0c48415..89146d7 100644 --- a/src/pctEnergyAdaptiveMLPFunction.cxx +++ b/src/pctEnergyAdaptiveMLPFunction.cxx @@ -138,10 +138,10 @@ EnergyAdaptiveMLPFunction ::EvaluateError(const double u, itk::Matrix #ifdef MLP_TIMING void -PolynomialMLPFunction ::PrintTiming(std::ostream & os) +PolynomialMLPFunction ::PrintTiming() { - os << "PolynomialMLPFunction timing:" << std::endl; - os << " EvaluateProbe1: " << m_EvaluateProbe1.GetTotal() << ' ' << m_EvaluateProbe1.GetUnit() << std::endl; + std::cout << "PolynomialMLPFunction timing:" << std::endl; + std::cout << " EvaluateProbe1: " << m_EvaluateProbe1.GetTotal() << ' ' << m_EvaluateProbe1.GetUnit() << std::endl; // os << " EvaluateProbe2: " << m_EvaluateProbe2.GetTotal() // << ' ' << m_EvaluateProbe2.GetUnit() << std::endl; } diff --git a/src/pctSchulteMLPFunction.cxx b/src/pctSchulteMLPFunction.cxx index a9d1146..6f8e8d4 100644 --- a/src/pctSchulteMLPFunction.cxx +++ b/src/pctSchulteMLPFunction.cxx @@ -222,11 +222,11 @@ SchulteMLPFunction ::EvaluateError(const double u, itk::Matrix & e #ifdef MLP_TIMING void -SchulteMLPFunction ::PrintTiming(std::ostream & os) +SchulteMLPFunction ::PrintTiming() { - os << "SchulteMLPFunction timing:" << std::endl; - os << " EvaluateProbe1: " << m_EvaluateProbe1.GetTotal() << ' ' << m_EvaluateProbe1.GetUnit() << std::endl; - os << " EvaluateProbe2: " << m_EvaluateProbe2.GetTotal() << ' ' << m_EvaluateProbe2.GetUnit() << std::endl; + std::cout << "SchulteMLPFunction timing:" << std::endl; + std::cout << " EvaluateProbe1: " << m_EvaluateProbe1.GetTotal() << ' ' << m_EvaluateProbe1.GetUnit() << std::endl; + std::cout << " EvaluateProbe2: " << m_EvaluateProbe2.GetTotal() << ' ' << m_EvaluateProbe2.GetUnit() << std::endl; } #endif From bf3ba285e40a744140a9bd3d5d3fd727d375d9f2 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Thu, 16 Oct 2025 09:06:45 +0100 Subject: [PATCH 06/12] BUG: Remove pctfdktwodweights from the wheel --- pyproject.toml | 1 - wrapping/__init_pct__.py | 1 - 2 files changed, 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 01b2e48..7e8869b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,7 +40,6 @@ dependencies = [ [project.scripts] pctfdk = "itk.pctfdk:main" -pctfdktwodweights = "itk.pctfdktwodweights:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 514320f..d0f6eaa 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -19,7 +19,6 @@ # Application modules _app_modules = [ "pctfdk", - "pctfdktwodweights", "pctpairprotons", "pctweplfit", ] From eaa35a0fdbceea9e79027ecb728c90ccf8d6d965 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 17 Oct 2025 09:53:43 +0200 Subject: [PATCH 07/12] ENH: Add pctbinning python application --- applications/pctbinning/pctbinning.py | 237 ++++++++++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_pct__.py | 1 + 3 files changed, 239 insertions(+) create mode 100644 applications/pctbinning/pctbinning.py diff --git a/applications/pctbinning/pctbinning.py b/applications/pctbinning/pctbinning.py new file mode 100644 index 0000000..5c009a0 --- /dev/null +++ b/applications/pctbinning/pctbinning.py @@ -0,0 +1,237 @@ +#!/usr/bin/env python +import argparse +import sys +from typing import Optional +import itk +from itk import PCT as pct + + +# Build a parser in the same style as pctpairprotons +def build_parser(): + parser = pct.PCTArgumentParser( + description="Create the 3D sequence (2D + distance) of proton radiographies" + ) + parser.add_argument( + "-i", + "--input", + help="Input file name containing the proton pairs", + required=True, + ) + parser.add_argument("-o", "--output", help="Output file name") + parser.add_argument("--elosswepl", help="Output file name (alias for --output)") + parser.add_argument( + "-c", "--count", help="Image of count of proton pairs per pixel" + ) + parser.add_argument( + "--scatwepl", help="Image of scattering WEPL of proton pairs per pixel" + ) + parser.add_argument("--noise", help="Image of WEPL variance per pixel") + parser.add_argument( + "-r", + "--robust", + help="Use robust estimation of scattering using 19.1 percentile.", + action="store_true", + default=False, + ) + parser.add_argument( + "-v", "--verbose", help="Verbose execution", action="store_true", default=False + ) + parser.add_argument("--config", help="Config file") + + parser.add_argument( + "-s", "--source", help="Source position", type=float, default=0.0 + ) + parser.add_argument( + "--quadricIn", + help=( + "Parameters of the entrance quadric support function, " + "see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" + ), + type=float, + nargs="+", + ) + parser.add_argument( + "--quadricOut", + help=( + "Parameters of the exit quadric support function, " + "see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" + ), + type=float, + nargs="+", + ) + parser.add_argument( + "--mlptype", + help="Type of most likely path (schulte, polynomial, or krah)", + choices=["schulte", "polynomial", "krah"], + default="schulte", + ) + parser.add_argument( + "--mlptrackeruncert", + help="Consider tracker uncertainties in MLP [Krah 2018, PMB]", + action="store_true", + default=False, + ) + parser.add_argument( + "--mlppolydeg", + help="Degree of the polynomial to approximate 1/beta^2p^2", + type=int, + default=5, + ) + parser.add_argument( + "--ionpot", + help="Ionization potential used in the reconstruction in eV", + type=float, + default=68.9984, + ) + parser.add_argument( + "--fill", + help="Fill holes, i.e. pixels that were not hit by protons", + action="store_true", + default=False, + ) + parser.add_argument( + "--trackerresolution", help="Tracker resolution in mm", type=float + ) + parser.add_argument( + "--trackerspacing", help="Tracker pair spacing in mm", type=float + ) + parser.add_argument( + "--materialbudget", help="Material budget x/X0 of tracker", type=float + ) + + parser.add_argument( + "--origin", help="Origin (default=centered)", type=float, nargs="+" + ) + parser.add_argument( + "--dimension", + help="Dimension(Deprecated) Use --size instead.", + type=int, + nargs="+", + default=[256], + ) + parser.add_argument("--size", help="Size", type=int, nargs="+", default=[256]) + parser.add_argument( + "--spacing", help="Spacing", type=float, nargs="+", default=[1.0] + ) + parser.add_argument("--direction", help="Direction", type=float, nargs="+") + parser.add_argument( + "--like", + help="Copy information from this image (origin, size, spacing, direction)", + ) + + return parser + + +def process(args_info: argparse.Namespace): + from itk import RTK as rtk + + if args_info.elosswepl and args_info.output: + print("Only --output or --elosswepl should be provided", file=sys.stderr) + sys.exit(1) + + OutputPixelType = itk.F + Dimension = 3 + OutputImageType = itk.Image[OutputPixelType, Dimension] + + max_threads = int(min(8, itk.MultiThreaderBase.GetGlobalMaximumNumberOfThreads())) + itk.MultiThreaderBase.SetGlobalMaximumNumberOfThreads(max_threads) + + # Create a stack of empty projection images + constantImageSource = rtk.ConstantImageSource[OutputImageType].New() + rtk.SetConstantImageSourceFromArgParse(constantImageSource, args_info) + + # Projection filter + projection = pct.ProtonPairsToDistanceDrivenProjection[ + OutputImageType, OutputImageType + ].New() + projection.SetInput(constantImageSource.GetOutput()) + projection.SetProtonPairsFileName(args_info.input) + projection.SetSourceDistance(args_info.source) + projection.SetMostLikelyPathType(args_info.mlptype) + projection.SetMostLikelyPathPolynomialDegree(args_info.mlppolydeg) + projection.SetMostLikelyPathTrackerUncertainties(args_info.mlptrackeruncert) + if args_info.trackerresolution is not None: + projection.SetTrackerResolution(args_info.trackerresolution) + if args_info.trackerspacing is not None: + projection.SetTrackerPairSpacing(args_info.trackerspacing) + if args_info.materialbudget is not None: + projection.SetMaterialBudget(args_info.materialbudget) + projection.SetIonizationPotential(args_info.ionpot * 1e-6) + projection.SetRobust(args_info.robust) + projection.SetComputeScattering(bool(args_info.scatwepl)) + projection.SetComputeNoise(bool(args_info.noise)) + + if args_info.quadricIn: + # quadric = object surface + qIn = rtk.QuadricShape.New() + qIn.SetA(args_info.quadricIn[0]) + qIn.SetB(args_info.quadricIn[1]) + qIn.SetC(args_info.quadricIn[2]) + qIn.SetD(args_info.quadricIn[3]) + qIn.SetE(args_info.quadricIn[4]) + qIn.SetF(args_info.quadricIn[5]) + qIn.SetG(args_info.quadricIn[6]) + qIn.SetH(args_info.quadricIn[7]) + qIn.SetI(args_info.quadricIn[8]) + qIn.SetJ(args_info.quadricIn[9]) + projection.SetQuadricIn(qIn) + if args_info.quadricOut: + qOut = rtk.QuadricShape.New() + qOut.SetA(args_info.quadricOut[0]) + qOut.SetB(args_info.quadricOut[1]) + qOut.SetC(args_info.quadricOut[2]) + qOut.SetD(args_info.quadricOut[3]) + qOut.SetE(args_info.quadricOut[4]) + qOut.SetF(args_info.quadricOut[5]) + qOut.SetG(args_info.quadricOut[6]) + qOut.SetH(args_info.quadricOut[7]) + qOut.SetI(args_info.quadricOut[8]) + qOut.SetJ(args_info.quadricOut[9]) + projection.SetQuadricOut(qOut) + + projection.Update() + + filler = pct.SmallHoleFiller[OutputImageType]() + if args_info.fill: + filler.SetImage(projection.GetOutput()) + filler.SetHolePixel(0.0) + filler.Fill() + + cii = itk.ChangeInformationImageFilter[OutputImageType].New() + if args_info.fill: + cii.SetInput(filler.GetOutput()) + else: + cii.SetInput(projection.GetOutput()) + cii.ChangeOriginOn() + cii.ChangeDirectionOn() + cii.ChangeSpacingOn() + cii.SetOutputDirection(projection.GetOutput().GetDirection()) + cii.SetOutputOrigin(projection.GetOutput().GetOrigin()) + cii.SetOutputSpacing(projection.GetOutput().GetSpacing()) + + if args_info.elosswepl or args_info.output: + # Write + out_name = args_info.elosswepl or args_info.output + itk.imwrite(cii.GetOutput(), out_name) + + if args_info.count: + # Write + itk.imwrite(projection.GetCount(), args_info.count) + + if args_info.scatwepl: + # Write + itk.imwrite(projection.GetAngle(), args_info.scatwepl) + + if args_info.noise: + # Write + itk.imwrite(projection.GetSquaredOutput(), args_info.noise) + + +def main(argv: Optional[list[str]] = 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 7e8869b..26ac8f0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ ] [project.scripts] +pctbinning = "itk.pctbinning:main" pctfdk = "itk.pctfdk:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index d0f6eaa..f695de2 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -18,6 +18,7 @@ # Application modules _app_modules = [ + "pctbinning", "pctfdk", "pctpairprotons", "pctweplfit", From 3cd8684e1d2a4e9ad24e78778a41c653228ee20d Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 17 Oct 2025 09:55:06 +0200 Subject: [PATCH 08/12] ENH: Use const pointers for image parameters --- include/SmallHoleFiller.h | 8 ++++---- include/SmallHoleFiller.hxx | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/SmallHoleFiller.h b/include/SmallHoleFiller.h index e2e15b9..6a4e35f 100644 --- a/include/SmallHoleFiller.h +++ b/include/SmallHoleFiller.h @@ -10,12 +10,12 @@ class SmallHoleFiller // Inputs void - SetImage(typename TImage::Pointer image); + SetImage(const TImage * image); void SetHolePixel(typename TImage::PixelType pixel); // Outputs - typename TImage::Pointer + TImage * GetOutput(); // This is the main loop. It simply calls Iterate() until complete. @@ -33,7 +33,7 @@ class SmallHoleFiller private: // The input image. - typename TImage::Pointer Image; + typename TImage::ConstPointer Image; // The intermediate and eventually output image. typename TImage::Pointer Output; @@ -45,7 +45,7 @@ class SmallHoleFiller // This function copies the data from 'input' to 'output' template void -DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output); +DeepCopy(const TImage * input, TImage * output); #ifndef ITK_MANUAL_INSTANTIATION # include "SmallHoleFiller.hxx" diff --git a/include/SmallHoleFiller.hxx b/include/SmallHoleFiller.hxx index 5de6d7e..6fdfc99 100644 --- a/include/SmallHoleFiller.hxx +++ b/include/SmallHoleFiller.hxx @@ -18,7 +18,7 @@ SmallHoleFiller::SmallHoleFiller() template void -SmallHoleFiller::SetImage(typename TImage::Pointer image) +SmallHoleFiller::SetImage(const TImage * image) { this->Image = image; this->Output = TImage::New(); @@ -32,10 +32,10 @@ SmallHoleFiller::SetHolePixel(typename TImage::PixelType pixel) } template -typename TImage::Pointer +TImage * SmallHoleFiller::GetOutput() { - return this->Output; + return this->Output.GetPointer(); } template @@ -153,7 +153,7 @@ SmallHoleFiller::HasEmptyPixels() /** Copy the input to the output*/ template void -DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output) +DeepCopy(const TImage * input, TImage * output) { output->SetRegions(input->GetLargestPossibleRegion()); output->Allocate(); From d3f9f3f23bfb5f08f89b4fbf22f2c4d6c06dde1d Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Fri, 17 Oct 2025 15:14:02 +0100 Subject: [PATCH 09/12] STYLE: Removed unused config file option --- .../pctbackprojectionbinning/pctbackprojectionbinning.ggo | 1 - applications/pctbinning/pctbinning.ggo | 1 - applications/pctfdk/pctfdk.ggo | 1 - applications/pctpaircuts/pctpaircuts.ggo | 1 - applications/pctpairgeometry/pctpairgeometry.ggo | 1 - applications/pctprojections/pctprojections.ggo | 1 - applications/pctzengbackprojections/pctzengbackprojections.ggo | 1 - 7 files changed, 7 deletions(-) diff --git a/applications/pctbackprojectionbinning/pctbackprojectionbinning.ggo b/applications/pctbackprojectionbinning/pctbackprojectionbinning.ggo index 084afdf..75f253b 100644 --- a/applications/pctbackprojectionbinning/pctbackprojectionbinning.ggo +++ b/applications/pctbackprojectionbinning/pctbackprojectionbinning.ggo @@ -2,7 +2,6 @@ package "pct" version "Backprojects the proton pairs according to mlp before applying mlp" option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "path" p "Path containing pair files" string yes option "regexp" r "Regular expression to select pair files in path" string yes option "output" o "Output file name" string yes diff --git a/applications/pctbinning/pctbinning.ggo b/applications/pctbinning/pctbinning.ggo index 86a4a06..9b37507 100644 --- a/applications/pctbinning/pctbinning.ggo +++ b/applications/pctbinning/pctbinning.ggo @@ -2,7 +2,6 @@ package "pct" version "Create the 3D sequence (2D + distance) of proton radiographies" option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "input" i "Input file name containing the proton pairs" string yes option "output" o "Output file name" string no option "elosswepl" - "Output file name (alias for --output)" string no diff --git a/applications/pctfdk/pctfdk.ggo b/applications/pctfdk/pctfdk.ggo index 7c51685..d4b6bb1 100644 --- a/applications/pctfdk/pctfdk.ggo +++ b/applications/pctfdk/pctfdk.ggo @@ -2,7 +2,6 @@ package "pct" version "Reconstruct a 3D volume from a sequence of projections [Feldkamp, David, Kress, 1984]." option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "geometry" g "XML geometry file name" string yes option "path" p "Path containing projections" string yes option "regexp" r "Regular expression to select projection files in path" string yes diff --git a/applications/pctpaircuts/pctpaircuts.ggo b/applications/pctpaircuts/pctpaircuts.ggo index 2d33b3f..f281d15 100644 --- a/applications/pctpaircuts/pctpaircuts.ggo +++ b/applications/pctpaircuts/pctpaircuts.ggo @@ -2,7 +2,6 @@ package "pct" version "Select proton pairs according to relative exit angle and energy [Schulte, MedPhys, 2008]." option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "input" i "Input file name containing the proton pairs" string yes option "output" o "Output file name" string yes option "source" s "Source position" double no default="0." diff --git a/applications/pctpairgeometry/pctpairgeometry.ggo b/applications/pctpairgeometry/pctpairgeometry.ggo index 834b9ae..0c376df 100644 --- a/applications/pctpairgeometry/pctpairgeometry.ggo +++ b/applications/pctpairgeometry/pctpairgeometry.ggo @@ -2,7 +2,6 @@ package "pct" version "Try to find source-to-detector1 and source-to-detector2 from list-mode information." option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "input" i "Input file name containing the proton pairs" string yes option "weplmax" - "WEPL threshold for protons in air." double no default="5." option "mindist" - "Minimum projected distance from central line" double no default="5." diff --git a/applications/pctprojections/pctprojections.ggo b/applications/pctprojections/pctprojections.ggo index bf43f2c..fbd0a6e 100644 --- a/applications/pctprojections/pctprojections.ggo +++ b/applications/pctprojections/pctprojections.ggo @@ -2,7 +2,6 @@ package "pct" version "Stack projections" option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "path" p "Path containing projections" string yes option "regexp" r "Regular expression to select projection files in path" string yes option "output" o "Output file name" string yes diff --git a/applications/pctzengbackprojections/pctzengbackprojections.ggo b/applications/pctzengbackprojections/pctzengbackprojections.ggo index ef60bca..11b47bc 100644 --- a/applications/pctzengbackprojections/pctzengbackprojections.ggo +++ b/applications/pctzengbackprojections/pctzengbackprojections.ggo @@ -2,7 +2,6 @@ package "pct" version "Computes the weighted backprojections described in [Zeng, Med Phys, 2007] from backprojections binned by direction" option "verbose" v "Verbose execution" flag off -option "config" - "Config file" string no option "input" i "Input file name" string yes option "outputc" c "Output file name for the bp weighted by the cosinus" string yes option "outputs" s "Output file name for the bp weighted by the sinus" string yes From 8ed496e1313220bdc57cf612ec40ed74e74787f5 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 17 Oct 2025 14:30:58 +0200 Subject: [PATCH 10/12] ENH: Centralize runtime output location for application targets in CMake --- applications/CMakeLists.txt | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 53f75e4..09b5308 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -24,3 +24,27 @@ add_subdirectory(pctpaircuts) add_subdirectory(pctpairgeometry) add_subdirectory(pctprojections) add_subdirectory(pctzengbackprojections) + +# Centralize runtime output location for application targets so we don't rely +# on the global CMAKE_RUNTIME_OUTPUT_DIRECTORY which can be overridden by ITK +set( + PCT_APPS + pctbackprojectionbinning + pctbinning + pctfdk + pctpaircuts + pctpairgeometry + pctprojections + pctzengbackprojections +) + +foreach(APP IN LISTS PCT_APPS) + if(TARGET ${APP}) + set_target_properties( + ${APP} + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" + ) + endif() +endforeach() From 62a07c5d19e652c8a24333628107110344ab4799 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 22 Oct 2025 13:59:59 +0200 Subject: [PATCH 11/12] ENH: Set the runtime output location inside each CMakeLists --- applications/CMakeLists.txt | 24 ------------------- .../pctbackprojectionbinning/CMakeLists.txt | 8 +++++++ applications/pctbinning/CMakeLists.txt | 8 +++++++ applications/pctfdk/CMakeLists.txt | 8 +++++++ applications/pctpaircuts/CMakeLists.txt | 8 +++++++ applications/pctpairgeometry/CMakeLists.txt | 8 +++++++ applications/pctprojections/CMakeLists.txt | 8 +++++++ .../pctzengbackprojections/CMakeLists.txt | 8 +++++++ 8 files changed, 56 insertions(+), 24 deletions(-) diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 09b5308..53f75e4 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -24,27 +24,3 @@ add_subdirectory(pctpaircuts) add_subdirectory(pctpairgeometry) add_subdirectory(pctprojections) add_subdirectory(pctzengbackprojections) - -# Centralize runtime output location for application targets so we don't rely -# on the global CMAKE_RUNTIME_OUTPUT_DIRECTORY which can be overridden by ITK -set( - PCT_APPS - pctbackprojectionbinning - pctbinning - pctfdk - pctpaircuts - pctpairgeometry - pctprojections - pctzengbackprojections -) - -foreach(APP IN LISTS PCT_APPS) - if(TARGET ${APP}) - set_target_properties( - ${APP} - PROPERTIES - RUNTIME_OUTPUT_DIRECTORY - "${PCT_BINARY_DIR}/bin" - ) - endif() -endforeach() diff --git a/applications/pctbackprojectionbinning/CMakeLists.txt b/applications/pctbackprojectionbinning/CMakeLists.txt index a6dfa13..d8ee376 100644 --- a/applications/pctbackprojectionbinning/CMakeLists.txt +++ b/applications/pctbackprojectionbinning/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctbackprojectionbinning PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctbackprojectionbinning + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctbinning/CMakeLists.txt b/applications/pctbinning/CMakeLists.txt index c141555..de0b4bf 100644 --- a/applications/pctbinning/CMakeLists.txt +++ b/applications/pctbinning/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctbinning PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctbinning + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctfdk/CMakeLists.txt b/applications/pctfdk/CMakeLists.txt index 6045c1d..1d2cd80 100644 --- a/applications/pctfdk/CMakeLists.txt +++ b/applications/pctfdk/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctfdk PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctfdk + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctpaircuts/CMakeLists.txt b/applications/pctpaircuts/CMakeLists.txt index d02d0cf..5e8958f 100644 --- a/applications/pctpaircuts/CMakeLists.txt +++ b/applications/pctpaircuts/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctpaircuts PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctpaircuts + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctpairgeometry/CMakeLists.txt b/applications/pctpairgeometry/CMakeLists.txt index 3b32131..265352f 100644 --- a/applications/pctpairgeometry/CMakeLists.txt +++ b/applications/pctpairgeometry/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctpairgeometry PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctpairgeometry + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctprojections/CMakeLists.txt b/applications/pctprojections/CMakeLists.txt index 5d20001..d02e0cc 100644 --- a/applications/pctprojections/CMakeLists.txt +++ b/applications/pctprojections/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctprojections PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctprojections + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS diff --git a/applications/pctzengbackprojections/CMakeLists.txt b/applications/pctzengbackprojections/CMakeLists.txt index e2ef6e1..3d15205 100644 --- a/applications/pctzengbackprojections/CMakeLists.txt +++ b/applications/pctzengbackprojections/CMakeLists.txt @@ -6,6 +6,14 @@ add_executable( ) target_link_libraries(pctzengbackprojections PCT) +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctzengbackprojections + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + # Installation code install( TARGETS From 9786ea95b8bb824d6f607a3997a355c9fe1c78e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 24 Oct 2025 20:05:35 +0200 Subject: [PATCH 12/12] ENH: Fit WEPL to angle standard deviation in pctweplfit --- applications/pctweplfit/pctweplfit.py | 55 ++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 10 deletions(-) diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index f5b77e9..18b58af 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -107,6 +107,7 @@ def add_detector(name, translation, attach_to_phantom=False): "EventID", "TrackID", "Position", + "Direction", "PreGlobalTime", "KineticEnergy", ] @@ -139,6 +140,7 @@ def process_phantom_length( tofs, wepls, elosses, + angles, verbose, ): import uproot @@ -146,6 +148,7 @@ def process_phantom_length( tofs_phantom_length = [] wepls_phantom_length = [] elosses_phantom_length = [] + angles_phantom_length = [] data = uproot.concatenate(f"{output}/l{int(phantom_length)}_*.root", library="np") pv( @@ -177,15 +180,20 @@ def process_phantom_length( ): continue + us = data["Position_X"][event_mask] + vs = data["Position_Y"][event_mask] + ws = data["Position_Z"][event_mask] + dus = data["Direction_X"][event_mask] + dvs = data["Direction_Y"][event_mask] + dws = data["Direction_Z"][event_mask] + times = data["PreGlobalTime"][event_mask] + tof = times[-1] - times[0] if phantom_length == 0.0: wepl = 0.0 else: - us = data["Position_X"][event_mask] - vs = data["Position_Y"][event_mask] - ws = data["Position_Z"][event_mask] if path_type == "phantom_length": # Length of the phantom wepl = phantom_length @@ -215,13 +223,38 @@ def process_phantom_length( data["KineticEnergy"][event_mask][0] - data["KineticEnergy"][event_mask][-1] ) + d_in_uw = np.array([dus[0], dws[0]]) + d_in_vw = np.array([dvs[0], dws[0]]) + d_out_uw = np.array([dus[-1], dws[-1]]) + d_out_vw = np.array([dvs[-1], dws[-1]]) + angle = [ + np.arccos( + np.minimum( + 1.0, + d_in_uw + @ d_out_uw + / (np.sqrt(d_in_uw @ d_in_uw) * np.sqrt(d_out_uw @ d_out_uw)), + ) + ), + np.arccos( + np.minimum( + 1.0, + d_in_vw + @ d_out_vw + / (np.sqrt(d_in_vw @ d_in_vw) * np.sqrt(d_out_vw @ d_out_vw)), + ) + ), + ] + tofs_phantom_length.append(tof) wepls_phantom_length.append(wepl) elosses_phantom_length.append(eloss) + angles_phantom_length += angle tofs[phantom_length] = tofs_phantom_length wepls[phantom_length] = wepls_phantom_length elosses[phantom_length] = elosses_phantom_length + angles[phantom_length] = angles_phantom_length def tof_fit( @@ -240,6 +273,7 @@ def tof_fit( tofs = manager.dict() wepls = manager.dict() elosses = manager.dict() + angles = manager.dict() results = [] with Pool() as pool: @@ -254,6 +288,7 @@ def tof_fit( tofs, wepls, elosses, + angles, verbose, ), ) @@ -263,9 +298,11 @@ def tof_fit( for result in results: result.get() - def fit(xs, ys, xlabel, ylabel): + def fit(xs, ys, xlabel, ylabel, xreduction=np.median): - xmedians = [np.median(xs[phantom_length]) for phantom_length in phantom_lengths] + xmedians = [ + xreduction(xs[phantom_length]) for phantom_length in phantom_lengths + ] ymedians = [np.median(ys[phantom_length]) for phantom_length in phantom_lengths] xpercentile25 = [ np.percentile(xs[phantom_length], 25) for phantom_length in phantom_lengths @@ -301,17 +338,14 @@ def fit(xs, ys, xlabel, ylabel): import matplotlib.pyplot as plt plt.figure() - plt.plot(xmedians, ymedians, "+", label="Medians") + plt.plot(xmedians, ymedians, "+") tof_xs = np.linspace(np.min(xmedians), np.max(xmedians), 100) for d, p in ps.items(): - plt.plot( - tof_xs, np.polyval(p, tof_xs), label=f"Polynomial fit (degree {d})" - ) + plt.plot(tof_xs, np.polyval(p, tof_xs)) plt.xlabel(xlabel) plt.ylabel(ylabel) - plt.legend() if savefig: plt.savefig(f"{output}/{xlabel}_to_{ylabel}_fit.pdf") @@ -326,6 +360,7 @@ def fit(xs, ys, xlabel, ylabel): fit(tofs, wepls, "tof", "wepl") fit(elosses, wepls, "eloss", "wepl") + fit(angles, wepls, xlabel="angle", ylabel="wepl", xreduction=np.std) def pctweplfit(