Helical4cylindrical

applications/CMakeLists.txt

@@ -102,6 +102,17 @@ add_subdirectory(rtkxradgeometry)
 add_subdirectory(rtkimagxgeometry)
 add_subdirectory(rtkorageometry)
 add_subdirectory(rtkbioscangeometry)
+
+#All executables below are specific to the tomo reconstruction from projections aquired along a helical trajectory of the gantry
+add_subdirectory(rtkhelicalgeometry)
+add_subdirectory(rtkkatsevichderivative)
+add_subdirectory(rtkkatsevichforwardbinning)
+add_subdirectory(rtkhilbertonkappalines)
+add_subdirectory(rtkpilines)
+add_subdirectory(rtkkatsevich)
+add_subdirectory(rtkkatsevichtrash)
+add_subdirectory(rtkkatsevichrecons)
+
 #=========================================================
 
 #-----------------------------------------------------------------------------

applications/rtkhelicalgeometry/CMakeLists.txt

@@ -0,0 +1,18 @@
+WRAP_GGO(rtkhelicalgeometry_GGO_C rtkhelicalgeometry.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo)
+add_executable(rtkhelicalgeometry rtkhelicalgeometry.cxx ${rtkhelicalgeometry_GGO_C})
+target_link_libraries(rtkhelicalgeometry RTK)
+
+# Installation code
+if(NOT RTK_INSTALL_NO_EXECUTABLES)
+  foreach(EXE_NAME rtkhelicalgeometry)
+    install(TARGETS ${EXE_NAME}
+      RUNTIME DESTINATION ${RTK_INSTALL_RUNTIME_DIR} COMPONENT Runtime
+      LIBRARY DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
+      ARCHIVE DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} COMPONENT Development)
+  endforeach()
+
+  # Install Python application
+  install(FILES rtkhelicalgeometry.py
+    DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT PythonWheelRuntimeLibraries)
+endif()
+

applications/rtkhelicalgeometry/rtkhelicalgeometry.cxx

@@ -0,0 +1,80 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "rtkhelicalgeometry_ggo.h"
+#include "rtkGgoFunctions.h"
+
+#include "rtkThreeDHelicalProjectionGeometryXMLFileWriter.h"
+
+int
+main(int argc, char * argv[])
+{
+  GGO(rtkhelicalgeometry, args_info);
+
+  // RTK geometry object
+  using GeometryType = rtk::ThreeDHelicalProjectionGeometry;
+  GeometryType::Pointer geometry = GeometryType::New();
+
+  // Projection matrices
+  for (int noProj = 0; noProj < args_info.nproj_arg; noProj++)
+  {
+    // Compute the angles
+    double angular_gap = args_info.arc_arg / args_info.nproj_arg;
+    double first_angle = 0.;
+    if (!args_info.first_angle_given)
+    {
+      first_angle = -0.5 * angular_gap * (args_info.nproj_arg - 1);
+    }
+    else
+      first_angle = args_info.first_angle_arg;
+
+    double angle = first_angle + noProj * angular_gap;
+
+
+    // Compute the vertical displacement
+    double vertical_coverage = args_info.arc_arg / 360 * args_info.pitch_arg;
+    double vertical_gap = vertical_coverage / args_info.nproj_arg;
+    double first_sy = 0.;
+    if (!args_info.first_sy_given)
+    {
+      first_sy = -0.5 * vertical_gap * (args_info.nproj_arg - 1);
+    }
+    else
+    {
+      first_sy = args_info.first_sy_arg;
+    }
+
+    double sy = first_sy + noProj * vertical_gap;
+
+    geometry->AddProjection(args_info.sid_arg, args_info.sdd_arg, angle, 0, sy, 0, 0, 0, sy);
+  }
+
+  // Set cylindrical detector radius
+  if (args_info.rad_cyl_given)
+    geometry->SetRadiusCylindricalDetector(args_info.rad_cyl_arg);
+
+
+  // Write
+  rtk::ThreeDHelicalProjectionGeometryXMLFileWriter::Pointer xmlWriter =
+    rtk::ThreeDHelicalProjectionGeometryXMLFileWriter::New();
+  xmlWriter->SetFilename(args_info.output_arg);
+  xmlWriter->SetObject(&(*geometry));
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(xmlWriter->WriteFile())
+
+  return EXIT_SUCCESS;
+}

applications/rtkhelicalgeometry/rtkhelicalgeometry.ggo

@@ -0,0 +1,16 @@
+package "rtkhelicalgeometry"
+purpose "Creates an RTK helical geometry file from regular helical trajectory. See http://www.openrtk.org/Doxygen/DocGeo3D.html for more information."
+
+option "verbose"     v "Verbose execution"                                                    flag   off
+option "config"      - "Config file"                                                          string no
+
+option "output"      o "Output file name"                                                     string yes
+option "first_angle" f "First angle in degrees (default = centered around 0)"                 double no
+option "first_sy"    y "First vertical position (default = centered around 0)"                double no 
+option "nproj"       n "Number of projections"                                                int    yes
+option "pitch"       p "Helix pitch (vertical displacement in one full (2pi) rotation"        double yes default="200"
+option "arc"         a "Angular arc covevered by the acquisition in degrees"                  double yes default="360"
+option "sdd"         - "Source to detector distance (mm)"                                     double no  default="1536"
+option "sid"         - "Source to isocenter distance (mm)"                                    double no  default="1000"
+option "rad_cyl"     - "Radius cylinder of cylindrical detector"                              double no  default="0"
+

applications/rtkhelicalgeometry/rtkhelicalgeometry.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+import argparse
+import sys
+from itk import RTK as rtk
+
+if __name__ == '__main__':
+  # Argument parsing
+  parser = argparse.ArgumentParser(description=
+    "Creates an RTK helical geometry file from regular helical trajectory. See http://www.openrtk.org/Doxygen/DocGeo3D.html for more information.")
+
+
+  parser.add_argument('--nproj', '-n', type=int, help='Number of projections')
+  parser.add_argument('--output', '-o', help='Output file name')
+  parser.add_argument('--verbose', '-v', type=bool, default=False, help='Verbose execution')
+  parser.add_argument('--config', '-c', help='Config file')
+  parser.add_argument('--first_angle', '-f', type=float, help='First angle in degrees')
+  parser.add_argument('--first_sy', '-y', type=float, help='First vertical position (default = centered around 0)')
+  parser.add_argument('--arc', '-a', type=float, default=360, help='Angular arc covevered by the acquisition in degrees')
+  parser.add_argument('--pitch', '-p', type=float, default=200, help='Helix pitch (vertical displacement in one full (2pi) rotation')
+  parser.add_argument('--sdd', type=float, default=1536, help='Source to detector distance (mm)')
+  parser.add_argument('--sid', type=float, default=1000, help='Source to isocenter distance (mm)')
+  parser.add_argument('--rad_cyl', type=float, default=0, help='Radius cylinder of cylindrical detector')
+
+  args = parser.parse_args()
+
+  if args.nproj is None or args.output is None :
+    parser.print_help()
+    sys.exit()
+
+  # Simulated Geometry
+  GeometryType = rtk.ThreeDCircularProjectionGeometry
+  geometry = GeometryType.New()
+
+  for noProj in range(0, args.nproj):
+
+    # Compute the angles
+    angular_gap = args.arc/args.nproj
+    if args.first_angle is None :
+      first_angle = -0.5*angular_gap*(args.nproj-1)
+    else :
+      first_angle = args.first_angle
+
+    angle = first_angle + noProj * angular_gap
+
+    # Compute vertical positions
+    vertical_coverage = args.arc/360.0*args.pitch
+    vertical_gap = vertical_coverage/args.nproj
+    if args.first_sy is None :
+      first_sy = -0.5*vertical_gap*(args.nproj-1)
+    else :
+      first_sy = args.first_sy
+
+    sy = first_sy + noProj * vertical_gap
+
+    geometry.AddProjection(args.sid,
+                           args.sdd,
+                           angle,
+                           0.,
+                           sy,
+                           0.,
+                           0.,
+                           0.,
+                           sy)
+
+  geometry.SetRadiusCylindricalDetector(args.rad_cyl)
+
+  writer = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()
+  writer.SetFilename(args.output)
+  writer.SetObject(geometry)
+  writer.WriteFile()

applications/rtkhilbertonkappalines/CMakeLists.txt

@@ -0,0 +1,14 @@
+WRAP_GGO(rtkhilbertonkappalines_GGO_C rtkhilbertonkappalines.ggo ../rtkinputprojections_section.ggo  ${RTK_BINARY_DIR}/rtkVersion.ggo)
+add_executable(rtkhilbertonkappalines rtkhilbertonkappalines.cxx ${rtkhilbertonkappalines_GGO_C})
+target_link_libraries(rtkhilbertonkappalines RTK)
+
+# Installation code
+if(NOT RTK_INSTALL_NO_EXECUTABLES)
+  foreach(EXE_NAME rtkhilbertonkappalines)
+    install(TARGETS ${EXE_NAME}
+      RUNTIME DESTINATION ${RTK_INSTALL_RUNTIME_DIR} COMPONENT Runtime
+      LIBRARY DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
+      ARCHIVE DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} COMPONENT Development)
+  endforeach()
+endif()
+

applications/rtkhilbertonkappalines/rtkhilbertonkappalines.cxx

@@ -0,0 +1,93 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "rtkhilbertonkappalines_ggo.h"
+#include "rtkGgoFunctions.h"
+#include "rtkConfiguration.h"
+
+#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h"
+#include "rtkHilbertTransformOnKappaLinesImageFilter.h"
+#include "rtkProgressCommands.h"
+
+#include <itkStreamingImageFilter.h>
+#include <itkImageRegionSplitterDirection.h>
+#include <itkImageFileWriter.h>
+
+int
+main(int argc, char * argv[])
+{
+  GGO(rtkhilbertonkappalines, args_info);
+
+  using OutputPixelType = float;
+  constexpr unsigned int Dimension = 3;
+
+  using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
+#ifdef RTK_USE_CUDA
+  using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
+#else
+  using OutputImageType = CPUOutputImageType;
+#endif
+
+  // Projections reader
+  using ReaderType = rtk::ProjectionsReader<OutputImageType>;
+  ReaderType::Pointer reader = ReaderType::New();
+  rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkhilbertonkappalines>(reader, args_info);
+
+  if (args_info.verbose_flag)
+    std::cout << "Reading... " << std::endl;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update())
+
+  // Geometry
+  if (args_info.verbose_flag)
+    std::cout << "Reading geometry information from " << args_info.geometry_arg << "..." << std::endl;
+  rtk::ThreeDHelicalProjectionGeometryXMLFileReader::Pointer geometryReader;
+  geometryReader = rtk::ThreeDHelicalProjectionGeometryXMLFileReader::New();
+  geometryReader->SetFilename(args_info.geometry_arg);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(geometryReader->GenerateOutputInformation())
+  rtk::ThreeDHelicalProjectionGeometry::Pointer geometry;
+  geometry = geometryReader->GetOutputObject();
+  geometry->VerifyHelixParameters();
+
+  // Check on hardware parameter
+#ifndef RTK_USE_CUDA
+  if (!strcmp(args_info.hardware_arg, "cuda"))
+  {
+    std::cerr << "The program has not been compiled with cuda option" << std::endl;
+    return EXIT_FAILURE;
+  }
+#endif
+
+  using HilbertTransformOnKappaLinesType =
+    rtk::HilbertTransformOnKappaLinesImageFilter<OutputImageType, OutputImageType>;
+  HilbertTransformOnKappaLinesType::Pointer fwd = HilbertTransformOnKappaLinesType::New();
+  fwd->SetGeometry(geometry);
+  fwd->SetInput(reader->GetOutput());
+
+  // Write
+  using WriterType = itk::ImageFileWriter<CPUOutputImageType>;
+  WriterType::Pointer writer = WriterType::New();
+  writer->SetFileName(args_info.output_arg);
+  writer->SetInput(fwd->GetOutput());
+
+  if (args_info.verbose_flag)
+    std::cout << "Reconstructing and writing... " << std::endl;
+
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
+
+  return EXIT_SUCCESS;
+}

applications/rtkhilbertonkappalines/rtkhilbertonkappalines.ggo

@@ -0,0 +1,9 @@
+package "rtkhilbertonkappalines"
+purpose "Compute the Hilbert transform of a stack of projections over kappa- lines (see Noo et al., PMB, 2003). Data is rebinned to apply a linear 1D Hilbert transform."
+
+option "verbose"    v "Verbose execution"                                           flag                         off
+option "config"     - "Config file"                                                 string                       no
+option "geometry"   g "XML geometry file name"                                      string                       yes
+option "output"     o "Output file name"                                            string                       yes
+option "hardware"   - "Hardware used for computation"                               values="cpu","cuda"          no   default="cpu"
+

applications/rtkkatsevich/CMakeLists.txt

@@ -0,0 +1,14 @@
+WRAP_GGO(rtkkatsevich_GGO_C rtkkatsevich.ggo ../rtkinputprojections_section.ggo ../rtk3Doutputimage_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo)
+add_executable(rtkkatsevich rtkkatsevich.cxx ${rtkkatsevich_GGO_C})
+target_link_libraries(rtkkatsevich RTK)
+
+# Installation code
+if(NOT RTK_INSTALL_NO_EXECUTABLES)
+  foreach(EXE_NAME rtkkatsevich)
+    install(TARGETS ${EXE_NAME}
+      RUNTIME DESTINATION ${RTK_INSTALL_RUNTIME_DIR} COMPONENT Runtime
+      LIBRARY DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
+      ARCHIVE DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} COMPONENT Development)
+  endforeach()
+endif()
+