From daac113bfac63fc9627e2a4df959cf4b7414536d Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 9 Apr 2026 16:27:56 +0200 Subject: [PATCH 1/3] ENH: Add pctinputprojections_group module --- applications/pctfdk/pctfdk.py | 32 +--- applications/pctinputprojections_group.py | 183 ++++++++++++++++++++++ wrapping/__init_pct__.py | 1 + 3 files changed, 188 insertions(+), 28 deletions(-) create mode 100644 applications/pctinputprojections_group.py diff --git a/applications/pctfdk/pctfdk.py b/applications/pctfdk/pctfdk.py index 23a3075e..53306aef 100644 --- a/applications/pctfdk/pctfdk.py +++ b/applications/pctfdk/pctfdk.py @@ -16,16 +16,6 @@ def build_parser(): 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 ) @@ -35,12 +25,6 @@ def build_parser(): 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( @@ -81,29 +65,21 @@ def build_parser(): help="Copy info from image (origin, size, spacing, direction)", type=str, ) + + pct.add_pctinputprojections_group(parser) + 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]) + pct.SetProjectionsReaderFromArgParse(reader, args_info) # Geometry if args_info.verbose: diff --git a/applications/pctinputprojections_group.py b/applications/pctinputprojections_group.py new file mode 100644 index 00000000..977ba0ea --- /dev/null +++ b/applications/pctinputprojections_group.py @@ -0,0 +1,183 @@ +import itk +from itk import PCT as pct +import numpy as np + +__all__ = [ + "add_pctinputprojections_group", + "GetProjectionsFileNamesFromArgParse", +] + + +# Mimicks pctinputprojections_section.ggo +def add_pctinputprojections_group(parser): + pctinputprojections_group = parser.add_argument_group( + "Input projections and their pre-processing" + ) + pctinputprojections_group.add_argument( + "--path", "-p", help="Path containing projections", required=True + ) + pctinputprojections_group.add_argument( + "--regexp", + "-r", + help="Regular expression to select projection files in path", + required=True, + ) + pctinputprojections_group.add_argument( + "--nsort", + help="Numeric sort for regular expression matches", + action="store_true", + ) + pctinputprojections_group.add_argument( + "--submatch", + help="Index of the submatch that will be used to sort matches", + type=int, + default=0, + ) + pctinputprojections_group.add_argument( + "--newdirection", + help="New value of input projections (before pre-processing)", + type=float, + nargs="+", + ) + pctinputprojections_group.add_argument( + "--neworigin", + help="New origin of input projections (before pre-processing)", + type=float, + nargs="+", + ) + pctinputprojections_group.add_argument( + "--newspacing", + help="New spacing of input projections (before pre-processing)", + type=float, + nargs="+", + ) + pctinputprojections_group.add_argument( + "--lowercrop", + help="Lower boundary crop size", + type=int, + nargs="+", + default=[0], + ) + pctinputprojections_group.add_argument( + "--uppercrop", + help="Upper boundary crop size", + type=int, + nargs="+", + default=[0], + ) + pctinputprojections_group.add_argument( + "--binning", + help="Shrink / Binning factos in each direction", + type=int, + nargs="+", + default=[1], + ) + pctinputprojections_group.add_argument( + "--wpc", + help="Water precorrection coefficients (default is no correction)", + type=float, + nargs="+", + ) + pctinputprojections_group.add_argument( + "--radius", + help="Radius of neighborhood for conditional median filtering", + type=int, + nargs="+", + default=[0], + ) + pctinputprojections_group.add_argument( + "--multiplier", + help="Threshold multiplier for conditional median filtering", + type=float, + default=0, + ) + + +# Mimicks GetProjectionsFileNamesFromGgo +def GetProjectionsFileNamesFromArgParse(args_info): + # Generate file names + names = itk.RegularExpressionSeriesFileNames.New() + names.SetDirectory(args_info.path) + names.SetNumericSort(args_info.nsort) + names.SetRegularExpression(args_info.regexp) + names.SetSubMatch(args_info.submatch) + + if args_info.verbose: + print(f"Regular expression matches {len(names.GetFileNames())} file(s)...") + + fileNames = [] + for fn in names.GetFileNames(): + imageio = itk.ImageIOFactory.CreateImageIO( + fn, itk.CommonEnums.IOFileMode_ReadMode + ) + if imageio is None: + print(f"Ignoring file: {fn}") + continue + fileNames.append(fn) + + return fileNames + + +def SetProjectionsReaderFromArgParse(reader, args_info): + fileNames = GetProjectionsFileNamesFromArgParse(args_info) + + # Vector component extraction (not in PCT ggo) + + # Change image information + Dimension = reader.GetOutput().GetImageDimension() + if args_info.newdirection is not None: + direction = [args_info.newdirection[0]] * 9 + for i in range(min(9, len(args_info.newdirection))): + direction[i] = args_info.newdirection[i] + direction = np.array(direction).reshape((3, 3)) + reader.SetDirection(itk.matrix_from_array(direction)) + + if args_info.newspacing is not None: + spacing = itk.Vector[itk.D, Dimension]() + spacing.Fill(args_info.newspacing[0]) + for i in range(len(args_info.newspacing)): + spacing[i] = args_info.newspacing[i] + reader.SetSpacing(spacing) + + if args_info.neworigin is not None: + origin = itk.Point[itk.D, Dimension]() + origin.Fill(args_info.neworigin[0]) + for i in range(len(args_info.neworigin)): + origin[i] = args_info.neworigin[i] + reader.SetOrigin(origin) + + # Crop boundaries + upperCrop = [0] * Dimension + lowerCrop = [0] * Dimension + if args_info.lowercrop is not None: + for i in range(len(args_info.lowercrop)): + lowerCrop[i] = args_info.lowercrop[i] + reader.SetLowerBoundaryCropSize(lowerCrop) + if args_info.uppercrop is not None: + for i in range(len(args_info.uppercrop)): + upperCrop[i] = args_info.uppercrop[i] + reader.SetUpperBoundaryCropSize(upperCrop) + + # Conditional median + medianRadius = reader.GetMedianRadius() + if args_info.radius is not None: + for i in range(len(args_info.radius)): + medianRadius[i] = args_info.radius[i] + reader.SetMedianRadius(medianRadius) + if args_info.multiplier is not None: + reader.SetConditionalMedianThresholdMultiplier(args_info.multiplier) + + # Shrink / Binning + binFactors = reader.GetShrinkFactors() + if args_info.binning is not None: + for i in range(len(args_info.binning)): + binFactors[i] = args_info.binning[i] + reader.SetShrinkFactors(binFactors) + + # Water precorrection + if args_info.wpc is not None: + reader.SetWaterPrecorrectionCoefficients(args_info.wpc) + + # Pass list to projections reader and update information + reader.SetFileNames(fileNames) + reader.UpdateOutputInformation() diff --git a/wrapping/__init_pct__.py b/wrapping/__init_pct__.py index 3466b211..df780597 100644 --- a/wrapping/__init_pct__.py +++ b/wrapping/__init_pct__.py @@ -12,6 +12,7 @@ pct_submodules = [ "itk.pctargumentparser", "itk.pctExtras", + "itk.pctinputprojections_group", ] for mod_name in pct_submodules: mod = importlib.import_module(mod_name) From f3aad8c6f2a1a7fab3294b02f7d184c1f75ff786 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 10 Apr 2026 10:24:29 +0200 Subject: [PATCH 2/3] ENH: Automates python application modules discovery --- wrapping/CMakeLists.txt | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/wrapping/CMakeLists.txt b/wrapping/CMakeLists.txt index 303fda63..d2c68198 100644 --- a/wrapping/CMakeLists.txt +++ b/wrapping/CMakeLists.txt @@ -8,24 +8,22 @@ configure_file( ) file( - GLOB PCT_PYTHON_APP + GLOB PCT_GROUP_SCRIPTS + CONFIGURE_DEPENDS + "${PCT_SOURCE_DIR}/applications/pct*_group.py" +) +file( + GLOB PCT_APP_SCRIPTS CONFIGURE_DEPENDS "${PCT_SOURCE_DIR}/applications/pct*/pct*.py" ) wrap_itk_python_bindings_install(/itk "PCT" __init_pct__.py - ${PCT_PYTHON_APP} + ${PCT_GROUP_SCRIPTS} + ${PCT_APP_SCRIPTS} + ${PCT_VERSION_SCRIPT} ${PCT_SOURCE_DIR}/wrapping/pctExtras.py ${PCT_SOURCE_DIR}/applications/pctargumentparser.py ${PCT_BINARY_DIR}/Wrapping/Generators/Python/pctConfig.py ) - -# Copy python applications to the ITK wrapping directory to ensure they can be imported in tests. -# This directory is added to the PYTHONPATH by itk_python_add_test. -if(ITK_DIR) - set(itk_wrap_python_binary_dir "${ITK_DIR}/Wrapping/Generators/Python") -else() - set(itk_wrap_python_binary_dir "${ITK_BINARY_DIR}/Wrapping/Generators/Python") -endif() -file(COPY ${PCT_PYTHON_APP} DESTINATION "${itk_wrap_python_binary_dir}") From c3a45206f78d7bcf88a451676b81eb72a7bd561a Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 16 Jul 2026 13:44:54 +0200 Subject: [PATCH 3/3] ENH: Add pctinputprojections_section module --- applications/pctfdk/CMakeLists.txt | 2 +- applications/pctfdk/pctfdk.cxx | 20 +-- applications/pctfdk/pctfdk.ggo | 3 - applications/pctinputprojections_section.ggo | 14 ++ include/pctGgoFunctions.h | 177 +++++++++++++++++++ 5 files changed, 194 insertions(+), 22 deletions(-) create mode 100644 applications/pctinputprojections_section.ggo create mode 100644 include/pctGgoFunctions.h diff --git a/applications/pctfdk/CMakeLists.txt b/applications/pctfdk/CMakeLists.txt index 48c90913..ab3a5e72 100644 --- a/applications/pctfdk/CMakeLists.txt +++ b/applications/pctfdk/CMakeLists.txt @@ -1,4 +1,4 @@ -wrap_ggo(pctfdk_GGO_C pctfdk.ggo) +wrap_ggo(pctfdk_GGO_C pctfdk.ggo ../pctinputprojections_section.ggo) add_executable( pctfdk pctfdk.cxx diff --git a/applications/pctfdk/pctfdk.cxx b/applications/pctfdk/pctfdk.cxx index 809673fe..0209c07a 100644 --- a/applications/pctfdk/pctfdk.cxx +++ b/applications/pctfdk/pctfdk.cxx @@ -1,4 +1,5 @@ #include "pctfdk_ggo.h" +#include "pctGgoFunctions.h" #include "rtkGgoFunctions.h" #include "rtkThreeDCircularProjectionGeometryXMLFile.h" @@ -7,7 +8,6 @@ #include "pctFDKDDConeBeamReconstructionFilter.h" #include "pctFDKDDConeBeamVarianceReconstructionFilter.h" -#include #include int @@ -23,26 +23,10 @@ main(int argc, char * argv[]) itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads( std::min(8, itk::MultiThreaderBase::GetGlobalMaximumNumberOfThreads())); - // Generate file names - itk::RegularExpressionSeriesFileNames::Pointer names = itk::RegularExpressionSeriesFileNames::New(); - names->SetDirectory(args_info.path_arg); - names->SetNumericSort(false); - names->SetRegularExpression(args_info.regexp_arg); - names->SetSubMatch(0); - - if (args_info.verbose_flag) - std::cout << "Regular expression matches " << names->GetFileNames().size() << " file(s)..." << std::endl; - // Projections reader using ProjectionImageType = itk::Image; auto reader = rtk::ProjectionsReader::New(); - reader->SetFileNames(names->GetFileNames()); - if (args_info.wpc_given) - { - std::vector coeffs; - coeffs.assign(args_info.wpc_arg, args_info.wpc_arg + args_info.wpc_given); - reader->SetWaterPrecorrectionCoefficients(coeffs); - } + pct::SetProjectionsReaderFromGgo, args_info_pctfdk>(reader, args_info); // Geometry if (args_info.verbose_flag) diff --git a/applications/pctfdk/pctfdk.ggo b/applications/pctfdk/pctfdk.ggo index 8d963e6e..3d2bafee 100644 --- a/applications/pctfdk/pctfdk.ggo +++ b/applications/pctfdk/pctfdk.ggo @@ -3,11 +3,8 @@ version "Reconstruct a 3D volume from a sequence of projections [Feldkamp, David option "verbose" v "Verbose execution" flag off 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 option "output" o "Output file name" string yes option "lowmem" l "Load only one projection per thread in memory" flag off -option "wpc" - "Water precorrection coefficients (default is no correction)" double multiple no section "Ramp filter" option "pad" - "Data padding parameter to correct for truncation" double no default="0.0" diff --git a/applications/pctinputprojections_section.ggo b/applications/pctinputprojections_section.ggo new file mode 100644 index 00000000..53bad7ef --- /dev/null +++ b/applications/pctinputprojections_section.ggo @@ -0,0 +1,14 @@ +section "Input projections and their pre-processing" +option "path" p "Path containing projections" string yes +option "regexp" r "Regular expression to select projection files in path" string yes +option "nsort" - "Numeric sort for regular expression matches" flag off +option "submatch" - "Index of the submatch that will be used to sort matches" int no default="0" +option "newdirection" - "New value of input projections (before pre-processing)" double multiple no +option "neworigin" - "New origin of input projections (before pre-processing)" double multiple no +option "newspacing" - "New spacing of input projections (before pre-processing)" double multiple no +option "lowercrop" - "Lower boundary crop size" int multiple no default="0" +option "uppercrop" - "Upper boundary crop size" int multiple no default="0" +option "binning" - "Shrink / Binning factos in each direction" int multiple no default="1" +option "wpc" - "Water precorrection coefficients (default is no correction)" double multiple no +option "radius" - "Radius of neighborhood for conditional median filtering" int multiple no default="0" +option "multiplier" - "Threshold multiplier for conditional median filtering" double no default="0" diff --git a/include/pctGgoFunctions.h b/include/pctGgoFunctions.h new file mode 100644 index 00000000..8bf014d3 --- /dev/null +++ b/include/pctGgoFunctions.h @@ -0,0 +1,177 @@ +/*========================================================================= + * + * Copyright PCT Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * 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. + * + *=========================================================================*/ + +#ifndef pctGgoFunctions_h +#define pctGgoFunctions_h + +#include "rtkMacro.h" +#include "rtkConstantImageSource.h" +#include "rtkIOFactories.h" +#include "rtkProjectionsReader.h" +#include +#include + +namespace pct +{ + +/** \brief Read a stack of 2D projections from gengetopt specifications. + * + * This function returns the file names of a projection series from command + * line options stored in ggo struct. + * The required options in the ggo struct are: + * - verbose + * - path: path containing projections + * - regexp: regular expression to select projection files in path + * - nsort: boolean to (des-)activate the numeric sort for expression matches + * - submatch: index of the submatch that will be used to sort matches + * + * \author Simon Rit + * + * \ingroup PCT + */ +template +std::vector +GetProjectionsFileNamesFromGgo(const TArgsInfo & args_info) +{ + auto names = itk::RegularExpressionSeriesFileNames::New(); + names->SetDirectory(args_info.path_arg); + names->SetNumericSort(args_info.nsort_flag); + names->SetRegularExpression(args_info.regexp_arg); + names->SetSubMatch(args_info.submatch_arg); + + if (args_info.verbose_flag) + std::cout << "Regular expression matches " << names->GetFileNames().size() << " file(s)..." << std::endl; + + if (args_info.submatch_given) + { + itksys::RegularExpression reg; + if (!reg.compile(args_info.regexp_arg)) + { + itkGenericExceptionMacro(<< "Error compiling regular expression " << args_info.regexp_arg); + } + + for (const std::string & name : names->GetFileNames()) + { + reg.find(name); + if (reg.match(args_info.submatch_arg) == std::string("")) + { + itkGenericExceptionMacro(<< "Cannot find submatch " << args_info.submatch_arg << " in " << name + << " from regular expression " << args_info.regexp_arg); + } + } + } + + std::vector fileNames = names->GetFileNames(); + rtk::RegisterIOFactories(); + std::vector idxtopop; + size_t i = 0; + for (const auto & fn : fileNames) + { + itk::ImageIOBase::Pointer imageio = + itk::ImageIOFactory::CreateImageIO(fn.c_str(), itk::ImageIOFactory::IOFileModeEnum::ReadMode); + + if (imageio.IsNull()) + { + std::cerr << "Ignoring file: " << fn << "\n"; + idxtopop.push_back(i); + } + i++; + } + std::reverse(idxtopop.begin(), idxtopop.end()); + for (const auto & id : idxtopop) + { + fileNames.erase(fileNames.begin() + id); + } + + return fileNames; +} + +template +void +SetProjectionsReaderFromGgo(TProjectionsReaderType * reader, const TArgsInfo & args_info) +{ + const std::vector fileNames = GetProjectionsFileNamesFromGgo(args_info); + + const unsigned int Dimension = TProjectionsReaderType::OutputImageType::GetImageDimension(); + + typename TProjectionsReaderType::OutputImageDirectionType direction; + if (args_info.newdirection_given) + { + direction.Fill(args_info.newdirection_arg[0]); + for (unsigned int i = 0; i < args_info.newdirection_given; i++) + direction[i / Dimension][i % Dimension] = args_info.newdirection_arg[i]; + reader->SetDirection(direction); + } + + typename TProjectionsReaderType::OutputImageSpacingType spacing; + if (args_info.newspacing_given) + { + spacing.Fill(args_info.newspacing_arg[0]); + for (unsigned int i = 0; i < args_info.newspacing_given; i++) + spacing[i] = args_info.newspacing_arg[i]; + reader->SetSpacing(spacing); + } + + typename TProjectionsReaderType::OutputImagePointType origin; + if (args_info.neworigin_given) + { + origin.Fill(args_info.neworigin_arg[0]); + for (unsigned int i = 0; i < args_info.neworigin_given; i++) + origin[i] = args_info.neworigin_arg[i]; + reader->SetOrigin(origin); + } + + auto lowerCrop = itk::MakeFilled(0); + for (unsigned int i = 0; i < args_info.lowercrop_given; i++) + lowerCrop[i] = args_info.lowercrop_arg[i]; + if (args_info.lowercrop_given) + reader->SetLowerBoundaryCropSize(lowerCrop); + + auto upperCrop = itk::MakeFilled(0); + for (unsigned int i = 0; i < args_info.uppercrop_given; i++) + upperCrop[i] = args_info.uppercrop_arg[i]; + if (args_info.uppercrop_given) + reader->SetUpperBoundaryCropSize(upperCrop); + + typename TProjectionsReaderType::MedianRadiusType medianRadius{}; + for (unsigned int i = 0; i < args_info.radius_given; i++) + medianRadius[i] = args_info.radius_arg[i]; + reader->SetMedianRadius(medianRadius); + if (args_info.multiplier_given) + reader->SetConditionalMedianThresholdMultiplier(args_info.multiplier_arg); + + typename TProjectionsReaderType::ShrinkFactorsType binFactors; + binFactors.Fill(1); + for (unsigned int i = 0; i < args_info.binning_given; i++) + binFactors[i] = args_info.binning_arg[i]; + reader->SetShrinkFactors(binFactors); + + if (args_info.wpc_given) + { + std::vector coeffs; + coeffs.assign(args_info.wpc_arg, args_info.wpc_arg + args_info.wpc_given); + reader->SetWaterPrecorrectionCoefficients(coeffs); + } + + reader->SetFileNames(fileNames); + TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->UpdateOutputInformation()); +} + +} // namespace pct + +#endif // pctGgoFunctions_h