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() 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/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/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.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/pctfdk/pctfdk.py b/applications/pctfdk/pctfdk.py new file mode 100644 index 0000000..181b73d --- /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() + + +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.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/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/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/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") 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 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(); 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/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 { 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,35 +114,16 @@ 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 +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/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/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..5ac59b3 100644 --- a/include/pctPolynomialMLPFunction.h +++ b/include/pctPolynomialMLPFunction.h @@ -166,10 +166,16 @@ class PCT_EXPORT PolynomialMLPFunction : public MostLikelyPathFunction void SetPolynomialDegree(const int polydeg); + void + SetCoefficients(const std::vector & coeffs); + + std::vector + GetCoefficients() const; + #ifdef MLP_TIMING /** Print timing information */ virtual void - PrintTiming(std::ostream & os) override; + PrintTiming() override; #endif protected: 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/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.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); 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/pyproject.toml b/pyproject.toml index 8e46967..26ac8f0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,8 @@ dependencies = [ ] [project.scripts] -pctfdktwodweights = "itk.pctfdktwodweights:main" +pctbinning = "itk.pctbinning:main" +pctfdk = "itk.pctfdk:main" pctpairprotons = "itk.pctpairprotons:main" pctweplfit = "itk.pctweplfit:main" 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 & coeffs) +{ + m_bm = coeffs; + m_PolynomialDegree = m_bm.size() - 1; + + m_PolynomialDegreePlusThree = m_PolynomialDegree + 3; +} + +std::vector +PolynomialMLPFunction ::GetCoefficients() const +{ + return m_bm; +} + void PolynomialMLPFunction ::Init(const VectorType posIn, const VectorType posOut, const VectorType dirIn, const VectorType dirOut) { + // Ensure we have coefficient values: if the user supplied custom + // coefficients via SetCoefficients(), m_bm will be non-empty and we'll + // keep them. Otherwise populate built-in coefficients for the current + // polynomial degree. + if (m_bm.empty()) + { + SetPolynomialDegree(m_PolynomialDegree); + } + m_uOrigin = posIn[2]; m_u2 = posOut[2] - m_uOrigin; @@ -187,10 +213,10 @@ PolynomialMLPFunction ::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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f9ffca5..464fe80 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,13 +3,20 @@ itk_module_test() #----------------------------------------------------------------------------- # CXX tests -set(PCTTests pctProtonPairsToDistanceDrivenProjectionTest.cxx) +set( + PCTTests + pctProtonPairsToDistanceDrivenProjectionTest.cxx + pctPolynomialMLPFunctionTest.cxx +) createtestdriver(PCT "${PCT-Test_LIBRARIES}" "${PCTTests}") itk_add_test(NAME pctProtonPairsToDistanceDrivenProjectionTest COMMAND PCTTestDriver pctProtonPairsToDistanceDrivenProjectionTest ) +itk_add_test(NAME pctPolynomialMLPFunctionTest + COMMAND PCTTestDriver pctPolynomialMLPFunctionTest +) #----------------------------------------------------------------------------- # Python tests diff --git a/test/pctPolynomialMLPFunctionTest.cxx b/test/pctPolynomialMLPFunctionTest.cxx new file mode 100644 index 0000000..2b28ab6 --- /dev/null +++ b/test/pctPolynomialMLPFunctionTest.cxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright Centre National de la Recherche Scientifique + * + * 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 + * + * https://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. + * + *=========================================================================*/ + +#include "pctPolynomialMLPFunction.h" + +#include "itkTestingMacros.h" + +#include +#include + +int +pctPolynomialMLPFunctionTest(int, char **) +{ + auto mlp = pct::PolynomialMLPFunction::New(); + + // Simple geometry for Init() + using VectorType = pct::MostLikelyPathFunction::VectorType; + VectorType posIn = itk::MakeVector(0.0, 0.0, 0.0); + VectorType posOut = itk::MakeVector(0.0, 0.0, 100.0); + VectorType dirIn = itk::MakeVector(0.0, 0.0, 1.0); + VectorType dirOut = itk::MakeVector(0.0, 0.0, 1.0); + + // Populate built-in coefficients + mlp->Init(posIn, posOut, dirIn, dirOut); + + // Default degree must be 5 + std::vector coeffs = mlp->GetCoefficients(); + if (coeffs.size() != 6) + { + std::cerr << "Default coefficients size != 6 (degree 5). Got " << coeffs.size() << std::endl; + return EXIT_FAILURE; + } + + // Set degree to 2 and re-Init + mlp->SetPolynomialDegree(2); + mlp->Init(posIn, posOut, dirIn, dirOut); + + std::vector u = { -50.0, 0.0, 50.0 }; + std::vector x, y; + + // Test MostLikelyPathFunction::Evaluate + mlp->Evaluate(u, x, y); + + // Basic sanity checks on Evaluate output: correct sizes + if (x.size() != u.size() || y.size() != u.size()) + { + std::cerr << "Evaluate returned vectors of wrong size." << std::endl; + return EXIT_FAILURE; + } + + // Custom coefficients (degree inferred from vector length). Use degree = 3 (4 coefficients). + std::vector custom = { 1.0e-06, 2.0e-07, 3.0e-08, 4.0e-09 }; + mlp->SetCoefficients(custom); + mlp->Init(posIn, posOut, dirIn, dirOut); + // Verify coefficients round-trip through GetCoefficients + std::vector got = mlp->GetCoefficients(); + if (got.size() != custom.size()) + { + std::cerr << "GetCoefficients size mismatch: " << got.size() << " != " << custom.size() << std::endl; + return EXIT_FAILURE; + } + for (unsigned int i = 0; i < got.size(); ++i) + { + if (got[i] != custom[i]) + { + std::cerr << "Coefficient mismatch at " << i << ": " << got[i] << " != " << custom[i] << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "pctPolynomialMLPFunctionTest PASSED" << std::endl; + return EXIT_SUCCESS; +} diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 19b90d5..f695de2 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -18,7 +18,8 @@ # Application modules _app_modules = [ - "pctfdktwodweights", + "pctbinning", + "pctfdk", "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/pctZengBackProjectionImageFilter.wrap b/wrapping/pctZengWeightedBackProjectionImageFilter.wrap similarity index 100% rename from wrapping/pctZengBackProjectionImageFilter.wrap rename to wrapping/pctZengWeightedBackProjectionImageFilter.wrap 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()