ENH: Add rtksubselect python application

Author: axel-grc

applications/README.md

@@ -21,6 +21,7 @@ The following are examples using RTK applications:
 ./rtkshowgeometry/README.md
 ./rtkosem/README.md
 ./rtksart/README.md
+./rtksubselect/README.md
 ```
 
 In [applications/rtktutorialapplication/](https://github.com/RTKConsortium/RTK/blob/main/applications/rtktutorialapplication), you will find a very basic RTK application that can be used as a starting point for building your own new application.

applications/rtksubselect/README.md

@@ -0,0 +1,11 @@
+# Projection selection
+
+This example illustrates the `rtksubselect` application, which extracts a subset of projections from an input stack and geometry and writes a reduced stack and matching geometry.
+
+![sin](../../documentation/docs/ExternalData/SheppLogan-Sinogram-3D.png){w=200px alt="SheppLogan sinogram"}
+![img](../../documentation/docs/ExternalData/SubSelect-Sinogram.png){w=200px alt="SubSelect sinogram"}
+
+This script uses the SheppLogan phantom.
+
+```{literalinclude} SubSelect.sh
+```

applications/rtksubselect/SubSelect.sh

@@ -0,0 +1,10 @@
+# Create a simulated geometry
+rtksimulatedgeometry -n 180 -o geometry.xml
+
+# Create projections of the Shepp–Logan phantom
+rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --size 200
+
+# Subselect every 2nd projection and write matching geometry
+rtksubselect -p . -r projections.mha --geometry geometry.xml \
+    --out-geometry geometry_subset.xml --out-proj projections_subset.mha \
+    --first 30 --last 140 --step 2

applications/rtksubselect/rtksubselect.py

@@ -0,0 +1,147 @@
+#!/usr/bin/env python
+import argparse
+import itk
+from itk import RTK as rtk
+
+
+def build_parser():
+    # argument parsing
+    parser = rtk.RTKArgumentParser(
+        description="Subselect projections from a stack and write updated geometry."
+    )
+
+    parser.add_argument(
+        "--geometry", "-g", help="XML geometry file name", type=str, required=True
+    )
+    parser.add_argument(
+        "--out_geometry", help="Output geometry file name", type=str, required=True
+    )
+    parser.add_argument(
+        "--out_proj", help="Output projections stack file name", type=str, required=True
+    )
+
+    parser.add_argument(
+        "--first", "-f", help="First projection index", type=int, default=0
+    )
+    parser.add_argument("--last", "-l", help="Last projection index", type=int)
+    parser.add_argument(
+        "--step", "-s", help="Step between projections", type=int, default=1
+    )
+    parser.add_argument(
+        "--list",
+        help="List of projection indices to keep (0-based)",
+        type=int,
+        nargs="+",
+    )
+    # RTK specific groups
+    rtk.add_rtkinputprojections_group(parser)
+
+    # Parse the command line arguments
+    return parser
+
+
+def process(args_info: argparse.Namespace):
+    OutputPixelType = itk.F
+    Dimension = 3
+    OutputImageType = itk.Image[OutputPixelType, Dimension]
+
+    # Projections reader
+    reader = rtk.ProjectionsReader[OutputImageType].New()
+    rtk.SetProjectionsReaderFromArgParse(reader, args_info)
+
+    # Geometry
+    if args_info.verbose:
+        print(f"Reading geometry information from {args_info.geometry}...")
+    geometry = rtk.read_geometry(args_info.geometry)
+
+    # Compute the indices of the selected projections
+    indices = []
+    n = len(geometry.GetGantryAngles())
+    if args_info.last:
+        n = min(args_info.last, n)
+    if args_info.list:
+        for i in args_info.list:
+            indices.append(i)
+    else:
+        start = args_info.first
+        step = args_info.step
+        for noProj in range(start, n, step):
+            indices.append(noProj)
+
+    # Output RTK geometry object
+    outputGeometry = rtk.ThreeDCircularProjectionGeometry.New()
+
+    # Output projections object
+    source = rtk.ConstantImageSource[OutputImageType].New()
+    source.SetInformationFromImage(reader.GetOutput())
+    outputSize = itk.size(reader.GetOutput())
+    outputSize_list = [s for s in outputSize]
+    outputSize_list[Dimension - 1] = len(indices)
+    source.SetSize(outputSize_list)
+    source.SetConstant(0)
+    source.Update()
+
+    # Fill in the outputGeometry and the output projections
+    paste = itk.PasteImageFilter[OutputImageType].New()
+    paste.SetSourceImage(reader.GetOutput())
+    paste.SetDestinationImage(source.GetOutput())
+
+    for i, src_idx in enumerate(indices):
+        # If it is not the first projection, we need to use the output of the paste filter as input
+        if i:
+            pimg = paste.GetOutput()
+            pimg.DisconnectPipeline()
+            paste.SetDestinationImage(pimg)
+
+        sourceRegion = reader.GetOutput().GetLargestPossibleRegion()
+        sourceRegion.SetIndex(Dimension - 1, src_idx)
+        sourceRegion.SetSize(Dimension - 1, 1)
+        paste.SetSourceRegion(sourceRegion)
+
+        destinationIndex = reader.GetOutput().GetLargestPossibleRegion().GetIndex()
+        destinationIndex.SetElement(Dimension - 1, i)
+        paste.SetDestinationIndex(destinationIndex)
+
+        paste.Update()
+
+        # Fill in the output geometry object
+        outputGeometry.SetRadiusCylindricalDetector(
+            geometry.GetRadiusCylindricalDetector()
+        )
+        outputGeometry.AddProjectionInRadians(
+            geometry.GetSourceToIsocenterDistances()[src_idx],
+            geometry.GetSourceToDetectorDistances()[src_idx],
+            geometry.GetGantryAngles()[src_idx],
+            geometry.GetProjectionOffsetsX()[src_idx],
+            geometry.GetProjectionOffsetsY()[src_idx],
+            geometry.GetOutOfPlaneAngles()[src_idx],
+            geometry.GetInPlaneAngles()[src_idx],
+            geometry.GetSourceOffsetsX()[src_idx],
+            geometry.GetSourceOffsetsY()[src_idx],
+        )
+        outputGeometry.SetCollimationOfLastProjection(
+            geometry.GetCollimationUInf()[src_idx],
+            geometry.GetCollimationUSup()[src_idx],
+            geometry.GetCollimationVInf()[src_idx],
+            geometry.GetCollimationVSup()[src_idx],
+        )
+
+    # Geometry writer
+    if args_info.verbose:
+        print(f"Writing geometry information in {args_info.out_geometry}...")
+    rtk.write_geometry(outputGeometry, args_info.out_geometry)
+
+    # Write
+    if args_info.verbose:
+        print(f"Writing selected projections to {args_info.out_proj}...")
+    itk.imwrite(paste.GetOutput(), args_info.out_proj)
+
+
+def main(argv=None):
+    parser = build_parser()
+    args_info = parser.parse_args(argv)
+    process(args_info)
+
+
+if __name__ == "__main__":
+    main()

documentation/docs/ExternalData/SubSelect-Sinogram.png.sha512

@@ -0,0 +1 @@
+7d5c53c19ccbc227fde52d86fd2a58c381ee8825b8fef4e9da87ef7688d6304d11af2287081ed11f6f7508ef096ac50519ac8a7065d1b265ec65fa1ad0c26f45

pyproject.toml

@@ -70,6 +70,7 @@ rtkprojectshepploganphantom= "itk.rtkprojectshepploganphantom:main"
 rtkshowgeometry = "itk.rtkshowgeometry:main"
 rtksart = "itk.rtksart:main"
 rtksimulatedgeometry = "itk.rtksimulatedgeometry:main"
+rtksubselect = "itk.rtksubselect:main"
 rtktotalvariationdenoising = "itk.rtktotalvariationdenoising:main"
 rtkvarianobigeometry = "itk.rtkvarianobigeometry:main"
 rtkwaveletsdenoising = "itk.rtkwaveletsdenoising:main"

wrapping/__init_rtk__.py

@@ -52,6 +52,7 @@
     "rtkshowgeometry",
     "rtksart",
     "rtksimulatedgeometry",
+    "rtksubselect",
     "rtktotalvariationdenoising",
     "rtkvarianobigeometry",
     "rtkwaveletsdenoising",