From 19970edbe14e0ab6a80294ffc99f3d40274f80cf Mon Sep 17 00:00:00 2001 From: Jerome Lesaint Date: Fri, 26 Mar 2021 16:08:50 +0100 Subject: [PATCH 1/3] WIP: Add Katsevich reconstruction from helical projections First attempt for an implementation of Katsevich algorithm as described in the paper "Exact helical reconstruction using native cone-beam geometries", Noo et al., PMB, 2003. Both flat panel and curved detector geometries are implemented though the curved geometry is still under debugging. As for the flat panel implementation, the whole thing more or less works. Some strange "tiling" artifacts remain. The algorithm runs in 4 steps : Compute a derivative of the projections Cosine-weight Forward rebinning + Hilbert transform (from ACoussat PR) + backward rebinning Backprojection. First commit - 2 First commit - nth Add post-cosine weight for curved detector filtering Various debug stuff and computation of FOV radius Misc Update Hilbert transform with ACoussant PR + minor changes requiring future cleanup Clean code from debugging couts Add a Katsevich Reconstruction filter and corresponding application --- .gitattributes | 2 +- .gitignore | 3 + applications/CMakeLists.txt | 11 + .../rtkhelicalgeometry/CMakeLists.txt | 18 + .../rtkhelicalgeometry/rtkhelicalgeometry.cxx | 80 ++ .../rtkhelicalgeometry/rtkhelicalgeometry.ggo | 16 + .../rtkhelicalgeometry/rtkhelicalgeometry.py | 70 ++ .../rtkhilbertonkappalines/CMakeLists.txt | 14 + .../rtkhilbertonkappalines.cxx | 93 ++ .../rtkhilbertonkappalines.ggo | 9 + applications/rtkkatsevich/CMakeLists.txt | 14 + applications/rtkkatsevich/rtkkatsevich.cxx | 198 +++++ applications/rtkkatsevich/rtkkatsevich.ggo | 12 + .../rtkkatsevichderivative/CMakeLists.txt | 14 + .../rtkkatsevichderivative.cxx | 92 ++ .../rtkkatsevichderivative.ggo | 9 + .../rtkkatsevichforwardbinning/CMakeLists.txt | 14 + .../rtkkatsevichforwardbinning.cxx | 92 ++ .../rtkkatsevichforwardbinning.ggo | 9 + .../rtkkatsevichrecons/CMakeLists.txt | 14 + .../rtkkatsevichrecons/rtkkatsevichrecons.cxx | 115 +++ .../rtkkatsevichrecons/rtkkatsevichrecons.ggo | 12 + applications/rtkkatsevichtrash/CMakeLists.txt | 14 + .../rtkkatsevichtrash/rtkkatsevichtrash.cxx | 166 ++++ .../rtkkatsevichtrash/rtkkatsevichtrash.ggo | 10 + applications/rtkpilines/CMakeLists.txt | 14 + applications/rtkpilines/rtkpilines.cxx | 153 ++++ applications/rtkpilines/rtkpilines.ggo | 12 + .../rtksimulatedgeometry.cxx | 1 + include/rtkFFTHilbertImageFilter.h | 112 +++ include/rtkFFTHilbertImageFilter.hxx | 101 +++ ...kHilbertTransformOnKappaLinesImageFilter.h | 108 +++ ...ilbertTransformOnKappaLinesImageFilter.hxx | 76 ++ .../rtkKatsevichBackProjectionImageFilter.h | 190 ++++ .../rtkKatsevichBackProjectionImageFilter.hxx | 840 ++++++++++++++++++ .../rtkKatsevichBackwardBinningImageFilter.h | 142 +++ ...rtkKatsevichBackwardBinningImageFilter.hxx | 255 ++++++ ...rtkKatsevichConeBeamReconstructionFilter.h | 215 +++++ ...kKatsevichConeBeamReconstructionFilter.hxx | 192 ++++ include/rtkKatsevichDerivativeImageFilter.h | 128 +++ include/rtkKatsevichDerivativeImageFilter.hxx | 253 ++++++ .../rtkKatsevichForwardBinningImageFilter.h | 135 +++ .../rtkKatsevichForwardBinningImageFilter.hxx | 232 +++++ include/rtkPILineImageFilter.h | 134 +++ include/rtkPILineImageFilter.hxx | 215 +++++ include/rtkThreeDCircularProjectionGeometry.h | 2 +- include/rtkThreeDHelicalProjectionGeometry.h | 162 ++++ ...eDHelicalProjectionGeometryXMLFileReader.h | 127 +++ ...eDHelicalProjectionGeometryXMLFileWriter.h | 92 ++ src/CMakeLists.txt | 3 + src/rtkThreeDHelicalProjectionGeometry.cxx | 249 ++++++ ...HelicalProjectionGeometryXMLFileReader.cxx | 186 ++++ ...HelicalProjectionGeometryXMLFileWriter.cxx | 236 +++++ 53 files changed, 5664 insertions(+), 2 deletions(-) create mode 100644 applications/rtkhelicalgeometry/CMakeLists.txt create mode 100644 applications/rtkhelicalgeometry/rtkhelicalgeometry.cxx create mode 100644 applications/rtkhelicalgeometry/rtkhelicalgeometry.ggo create mode 100644 applications/rtkhelicalgeometry/rtkhelicalgeometry.py create mode 100644 applications/rtkhilbertonkappalines/CMakeLists.txt create mode 100644 applications/rtkhilbertonkappalines/rtkhilbertonkappalines.cxx create mode 100644 applications/rtkhilbertonkappalines/rtkhilbertonkappalines.ggo create mode 100644 applications/rtkkatsevich/CMakeLists.txt create mode 100644 applications/rtkkatsevich/rtkkatsevich.cxx create mode 100644 applications/rtkkatsevich/rtkkatsevich.ggo create mode 100644 applications/rtkkatsevichderivative/CMakeLists.txt create mode 100644 applications/rtkkatsevichderivative/rtkkatsevichderivative.cxx create mode 100644 applications/rtkkatsevichderivative/rtkkatsevichderivative.ggo create mode 100644 applications/rtkkatsevichforwardbinning/CMakeLists.txt create mode 100644 applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.cxx create mode 100644 applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.ggo create mode 100644 applications/rtkkatsevichrecons/CMakeLists.txt create mode 100644 applications/rtkkatsevichrecons/rtkkatsevichrecons.cxx create mode 100644 applications/rtkkatsevichrecons/rtkkatsevichrecons.ggo create mode 100644 applications/rtkkatsevichtrash/CMakeLists.txt create mode 100644 applications/rtkkatsevichtrash/rtkkatsevichtrash.cxx create mode 100644 applications/rtkkatsevichtrash/rtkkatsevichtrash.ggo create mode 100644 applications/rtkpilines/CMakeLists.txt create mode 100644 applications/rtkpilines/rtkpilines.cxx create mode 100644 applications/rtkpilines/rtkpilines.ggo create mode 100644 include/rtkFFTHilbertImageFilter.h create mode 100644 include/rtkFFTHilbertImageFilter.hxx create mode 100644 include/rtkHilbertTransformOnKappaLinesImageFilter.h create mode 100644 include/rtkHilbertTransformOnKappaLinesImageFilter.hxx create mode 100644 include/rtkKatsevichBackProjectionImageFilter.h create mode 100644 include/rtkKatsevichBackProjectionImageFilter.hxx create mode 100644 include/rtkKatsevichBackwardBinningImageFilter.h create mode 100644 include/rtkKatsevichBackwardBinningImageFilter.hxx create mode 100644 include/rtkKatsevichConeBeamReconstructionFilter.h create mode 100644 include/rtkKatsevichConeBeamReconstructionFilter.hxx create mode 100644 include/rtkKatsevichDerivativeImageFilter.h create mode 100644 include/rtkKatsevichDerivativeImageFilter.hxx create mode 100644 include/rtkKatsevichForwardBinningImageFilter.h create mode 100644 include/rtkKatsevichForwardBinningImageFilter.hxx create mode 100644 include/rtkPILineImageFilter.h create mode 100644 include/rtkPILineImageFilter.hxx create mode 100644 include/rtkThreeDHelicalProjectionGeometry.h create mode 100644 include/rtkThreeDHelicalProjectionGeometryXMLFileReader.h create mode 100644 include/rtkThreeDHelicalProjectionGeometryXMLFileWriter.h create mode 100644 src/rtkThreeDHelicalProjectionGeometry.cxx create mode 100644 src/rtkThreeDHelicalProjectionGeometryXMLFileReader.cxx create mode 100644 src/rtkThreeDHelicalProjectionGeometryXMLFileWriter.cxx diff --git a/.gitattributes b/.gitattributes index d702b5308..1a279742b 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,6 +1,6 @@ .git* export-ignore # Custom attribute to mark sources as using our C++/C code style. -[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat,uncrustify +[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat cmake/*.cxx our-c-style include/*.h* our-c-style diff --git a/.gitignore b/.gitignore index 479088e28..bde0b8fdd 100644 --- a/.gitignore +++ b/.gitignore @@ -33,3 +33,6 @@ CMakeLists.txt.user* # Mac System File .DS_Store + + +.project diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index dd1e0932d..730fddf3b 100644 --- a/applications/CMakeLists.txt +++ b/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) + #========================================================= #----------------------------------------------------------------------------- diff --git a/applications/rtkhelicalgeometry/CMakeLists.txt b/applications/rtkhelicalgeometry/CMakeLists.txt new file mode 100644 index 000000000..dca3fc510 --- /dev/null +++ b/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() + diff --git a/applications/rtkhelicalgeometry/rtkhelicalgeometry.cxx b/applications/rtkhelicalgeometry/rtkhelicalgeometry.cxx new file mode 100644 index 000000000..a7f3000d6 --- /dev/null +++ b/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; +} diff --git a/applications/rtkhelicalgeometry/rtkhelicalgeometry.ggo b/applications/rtkhelicalgeometry/rtkhelicalgeometry.ggo new file mode 100644 index 000000000..78e96a7e6 --- /dev/null +++ b/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" + diff --git a/applications/rtkhelicalgeometry/rtkhelicalgeometry.py b/applications/rtkhelicalgeometry/rtkhelicalgeometry.py new file mode 100644 index 000000000..bc6a62f46 --- /dev/null +++ b/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() diff --git a/applications/rtkhilbertonkappalines/CMakeLists.txt b/applications/rtkhilbertonkappalines/CMakeLists.txt new file mode 100644 index 000000000..1a5bbfa20 --- /dev/null +++ b/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() + diff --git a/applications/rtkhilbertonkappalines/rtkhilbertonkappalines.cxx b/applications/rtkhilbertonkappalines/rtkhilbertonkappalines.cxx new file mode 100644 index 000000000..8fb63715d --- /dev/null +++ b/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 +#include +#include + +int +main(int argc, char * argv[]) +{ + GGO(rtkhilbertonkappalines, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(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; + HilbertTransformOnKappaLinesType::Pointer fwd = HilbertTransformOnKappaLinesType::New(); + fwd->SetGeometry(geometry); + fwd->SetInput(reader->GetOutput()); + + // Write + using WriterType = itk::ImageFileWriter; + 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; +} diff --git a/applications/rtkhilbertonkappalines/rtkhilbertonkappalines.ggo b/applications/rtkhilbertonkappalines/rtkhilbertonkappalines.ggo new file mode 100644 index 000000000..92e1a2f86 --- /dev/null +++ b/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" + diff --git a/applications/rtkkatsevich/CMakeLists.txt b/applications/rtkkatsevich/CMakeLists.txt new file mode 100644 index 000000000..b778e7951 --- /dev/null +++ b/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() + diff --git a/applications/rtkkatsevich/rtkkatsevich.cxx b/applications/rtkkatsevich/rtkkatsevich.cxx new file mode 100644 index 000000000..d46a063f1 --- /dev/null +++ b/applications/rtkkatsevich/rtkkatsevich.cxx @@ -0,0 +1,198 @@ +/*========================================================================= + * + * 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 "rtkkatsevich_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkCyclicDeformationImageFilter.h" +#include "rtkProgressCommands.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + + +int +main(int argc, char * argv[]) +{ + GGO(rtkkatsevich, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(reader, args_info); + + if (!args_info.lowmem_flag) + { + if (args_info.verbose_flag) + std::cout << "Reading... " << std::endl; + TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update()) + } + + using DerivativeFilterType = rtk::KatsevichDerivativeImageFilter; + using ForwardBinningType = rtk::KatsevichForwardBinningImageFilter; + using HilbertFilterType = rtk::FFTHilbertImageFilter; + using BackwardBinningType = rtk::KatsevichBackwardBinningImageFilter; + using BackProjectionType = rtk::KatsevichBackProjectionImageFilter; + + + // 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(); + + // Katsevich filtering steps + std::cout << "Starting derivative..." << std::endl; + DerivativeFilterType::Pointer deriv = DerivativeFilterType::New(); + deriv->SetGeometry(geometry); + deriv->SetInput(reader->GetOutput()); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + + if (args_info.writefiles_flag) + { + writer->SetFileName("deriv.mha"); + writer->SetInput(deriv->GetOutput()); + writer->Update(); + } + + + std::cout << "Starting forward binning..." << std::endl; + ForwardBinningType::Pointer fwdbin = ForwardBinningType::New(); + fwdbin->SetGeometry(geometry); + fwdbin->SetInput(deriv->GetOutput()); + + if (args_info.writefiles_flag) + { + writer->SetFileName("forward.mha"); + fwdbin->Update(); + OutputImageType::Pointer im = fwdbin->GetOutput(); + // OutputImageType::SpacingType sp = im->GetSpacing(); + // sp[1] = 1.; + // im->SetSpacing(sp); + // std::cout << "spacing " << im->GetSpacing() << std::endl; + writer->SetInput(im); + writer->Update(); + } + + std::cout << "Starting Hilbert..." << std::endl; + HilbertFilterType::Pointer hilbert = HilbertFilterType::New(); + // hilbert->SetGeometry(geometryReader->GetOutputObject()); + hilbert->SetInput(fwdbin->GetOutput()); + hilbert->SetPixelShift(0.5); + + + if (args_info.writefiles_flag) + { + writer->SetFileName("hilbert.mha"); + hilbert->Update(); + OutputImageType::Pointer im = hilbert->GetOutput(); + // OutputImageType::SpacingType sp = im->GetSpacing(); + // sp[1] = 1.; + // im->SetSpacing(sp); + writer->SetInput(im); + writer->Update(); + } + + std::cout << "Starting backward binning..." << std::endl; + BackwardBinningType::Pointer bwdbin = BackwardBinningType::New(); + bwdbin->SetGeometry(geometryReader->GetOutputObject()); + bwdbin->SetInput(hilbert->GetOutput()); + + if (args_info.writefiles_flag) + { + writer->SetFileName("backward.mha"); + writer->SetInput(bwdbin->GetOutput()); + writer->Update(); + } + // Create reconstructed image + using ConstantImageSourceType = rtk::ConstantImageSource; + ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New(); + rtk::SetConstantImageSourceFromGgo(constantImageSource, args_info); + + // Katsevich back-projection + std::cout << "Starting bp..." << std::endl; + BackProjectionType::Pointer bp = BackProjectionType::New(); + bp->SetGeometry(geometryReader->GetOutputObject()); + bp->SetInput(0, constantImageSource->GetOutput()); + bp->SetInput(1, bwdbin->GetOutput()); + + if (args_info.writefiles_flag) + { + writer->SetFileName("bp.mha"); + writer->SetInput(bp->GetOutput()); + writer->Update(); + } + + +// 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 + + //// Streaming depending on streaming capability of writer + // using StreamerType = itk::StreamingImageFilter; + // StreamerType::Pointer streamerBP = StreamerType::New(); + // streamerBP->SetInput(bp->GetOutput()); + // streamerBP->SetNumberOfStreamDivisions(args_info.divisions_arg); + // itk::ImageRegionSplitterDirection::Pointer splitter = itk::ImageRegionSplitterDirection::New(); + // splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis + // streamerBP->SetRegionSplitter(splitter); + + // Write + // using WriterType = itk::ImageFileWriter; + // WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output_arg); + writer->SetInput(bp->GetOutput()); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + return EXIT_SUCCESS; +} diff --git a/applications/rtkkatsevich/rtkkatsevich.ggo b/applications/rtkkatsevich/rtkkatsevich.ggo new file mode 100644 index 000000000..39ca65a64 --- /dev/null +++ b/applications/rtkkatsevich/rtkkatsevich.ggo @@ -0,0 +1,12 @@ +package "rtkfdk" +purpose "Reconstructs a 3D volume from a sequence of projections [Feldkamp, David, Kress, 1984]." + +option "verbose" v "Verbose execution" flag off +option "writefiles" w "Write intermediate files to disk" 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" +option "lowmem" l "Load only one projection per thread in memory" flag off +option "divisions" d "Streaming option: number of stream divisions of the CT" int no default="1" +option "subsetsize" - "Streaming option: number of projections processed at a time" int no default="16" diff --git a/applications/rtkkatsevichderivative/CMakeLists.txt b/applications/rtkkatsevichderivative/CMakeLists.txt new file mode 100644 index 000000000..80f496f19 --- /dev/null +++ b/applications/rtkkatsevichderivative/CMakeLists.txt @@ -0,0 +1,14 @@ +WRAP_GGO(rtkkatsevichderivative_GGO_C rtkkatsevichderivative.ggo ../rtkinputprojections_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo) +add_executable(rtkkatsevichderivative rtkkatsevichderivative.cxx ${rtkkatsevichderivative_GGO_C}) +target_link_libraries(rtkkatsevichderivative RTK) + +# Installation code +if(NOT RTK_INSTALL_NO_EXECUTABLES) + foreach(EXE_NAME rtkkatsevichderivative) + 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() + diff --git a/applications/rtkkatsevichderivative/rtkkatsevichderivative.cxx b/applications/rtkkatsevichderivative/rtkkatsevichderivative.cxx new file mode 100644 index 000000000..0aa032250 --- /dev/null +++ b/applications/rtkkatsevichderivative/rtkkatsevichderivative.cxx @@ -0,0 +1,92 @@ +/*========================================================================= + * + * 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 "rtkkatsevichderivative_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkKatsevichDerivativeImageFilter.h" +#include "rtkProgressCommands.h" + +#include +#include +#include + +int +main(int argc, char * argv[]) +{ + GGO(rtkkatsevichderivative, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(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 KatsevichDerivativeType = rtk::KatsevichDerivativeImageFilter; + KatsevichDerivativeType::Pointer deriv = KatsevichDerivativeType::New(); + deriv->SetGeometry(geometry); + deriv->SetInput(reader->GetOutput()); + + // Write + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output_arg); + writer->SetInput(deriv->GetOutput()); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + return EXIT_SUCCESS; +} diff --git a/applications/rtkkatsevichderivative/rtkkatsevichderivative.ggo b/applications/rtkkatsevichderivative/rtkkatsevichderivative.ggo new file mode 100644 index 000000000..7a5cad707 --- /dev/null +++ b/applications/rtkkatsevichderivative/rtkkatsevichderivative.ggo @@ -0,0 +1,9 @@ +package "rtkkatsevichderivative" +purpose "Computes the derivative of a stack of projections wrt to the projection index (rotation angle). See Noo et al., PMB, 2003" + +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" + diff --git a/applications/rtkkatsevichforwardbinning/CMakeLists.txt b/applications/rtkkatsevichforwardbinning/CMakeLists.txt new file mode 100644 index 000000000..aa8347cd9 --- /dev/null +++ b/applications/rtkkatsevichforwardbinning/CMakeLists.txt @@ -0,0 +1,14 @@ +WRAP_GGO(rtkkatsevichforwardbinning_GGO_C rtkkatsevichforwardbinning.ggo ../rtkinputprojections_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo) +add_executable(rtkkatsevichforwardbinning rtkkatsevichforwardbinning.cxx ${rtkkatsevichforwardbinning_GGO_C}) +target_link_libraries(rtkkatsevichforwardbinning RTK) + +# Installation code +if(NOT RTK_INSTALL_NO_EXECUTABLES) + foreach(EXE_NAME rtkkatsevichforwardbinning) + 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() + diff --git a/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.cxx b/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.cxx new file mode 100644 index 000000000..26a59b196 --- /dev/null +++ b/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.cxx @@ -0,0 +1,92 @@ +/*========================================================================= + * + * 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 "rtkkatsevichforwardbinning_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkKatsevichForwardBinningImageFilter.h" +#include "rtkProgressCommands.h" + +#include +#include +#include + +int +main(int argc, char * argv[]) +{ + GGO(rtkkatsevichforwardbinning, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(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 KatsevichForwardBinningType = rtk::KatsevichForwardBinningImageFilter; + KatsevichForwardBinningType::Pointer fwd = KatsevichForwardBinningType::New(); + fwd->SetGeometry(geometry); + fwd->SetInput(reader->GetOutput()); + + // Write + using WriterType = itk::ImageFileWriter; + 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; +} diff --git a/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.ggo b/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.ggo new file mode 100644 index 000000000..a7c704082 --- /dev/null +++ b/applications/rtkkatsevichforwardbinning/rtkkatsevichforwardbinning.ggo @@ -0,0 +1,9 @@ +package "rtkkatsevichforwardbinning" +purpose "Re-bin a stack of projections to the psi-variable. See Noo et al., PMB, 2003" + +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" + diff --git a/applications/rtkkatsevichrecons/CMakeLists.txt b/applications/rtkkatsevichrecons/CMakeLists.txt new file mode 100644 index 000000000..9f6902112 --- /dev/null +++ b/applications/rtkkatsevichrecons/CMakeLists.txt @@ -0,0 +1,14 @@ +WRAP_GGO(rtkkatsevichrecons_GGO_C rtkkatsevichrecons.ggo ../rtkinputprojections_section.ggo ../rtk3Doutputimage_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo) +add_executable(rtkkatsevichrecons rtkkatsevichrecons.cxx ${rtkkatsevichrecons_GGO_C}) +target_link_libraries(rtkkatsevichrecons RTK) + +# Installation code +if(NOT RTK_INSTALL_NO_EXECUTABLES) + foreach(EXE_NAME rtkkatsevichrecons) + 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() + diff --git a/applications/rtkkatsevichrecons/rtkkatsevichrecons.cxx b/applications/rtkkatsevichrecons/rtkkatsevichrecons.cxx new file mode 100644 index 000000000..fe7b181ac --- /dev/null +++ b/applications/rtkkatsevichrecons/rtkkatsevichrecons.cxx @@ -0,0 +1,115 @@ +/*========================================================================= + * + * 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 "rtkkatsevichrecons_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkCyclicDeformationImageFilter.h" +#include "rtkProgressCommands.h" + +#include +#include +#include + +#include + + +int +main(int argc, char * argv[]) +{ + GGO(rtkkatsevichrecons, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(reader, args_info); + + if (!args_info.lowmem_flag) + { + 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(); + + // Create reconstructed image + using ConstantImageSourceType = rtk::ConstantImageSource; + ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New(); + rtk::SetConstantImageSourceFromGgo(constantImageSource, + args_info); + + using KatsevichReconsFilterType = rtk::KatsevichConeBeamReconstructionFilter; + KatsevichReconsFilterType::Pointer kats = KatsevichReconsFilterType::New(); + kats->SetGeometry(geometry); + kats->SetInput(0, constantImageSource->GetOutput()); + kats->SetInput(1, reader->GetOutput()); + + +// 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 + + //// Streaming depending on streaming capability of writer + // using StreamerType = itk::StreamingImageFilter; + // StreamerType::Pointer streamerBP = StreamerType::New(); + // streamerBP->SetInput(bp->GetOutput()); + // streamerBP->SetNumberOfStreamDivisions(args_info.divisions_arg); + // itk::ImageRegionSplitterDirection::Pointer splitter = itk::ImageRegionSplitterDirection::New(); + // splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis + // streamerBP->SetRegionSplitter(splitter); + + // Write + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output_arg); + writer->SetInput(kats->GetOutput()); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + return EXIT_SUCCESS; +} diff --git a/applications/rtkkatsevichrecons/rtkkatsevichrecons.ggo b/applications/rtkkatsevichrecons/rtkkatsevichrecons.ggo new file mode 100644 index 000000000..bb30328a3 --- /dev/null +++ b/applications/rtkkatsevichrecons/rtkkatsevichrecons.ggo @@ -0,0 +1,12 @@ +package "rtkkatsevichrecons" +purpose "Reconstructs a 3D volume from a sequence of projections acquired along a helical trajectory [Noo et al., PMB, 2003]." + +option "verbose" v "Verbose execution" flag off +option "writefiles" w "Write intermediate files to disk" 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" +option "lowmem" l "Load only one projection per thread in memory" flag off +option "divisions" d "Streaming option: number of stream divisions of the CT" int no default="1" +option "subsetsize" - "Streaming option: number of projections processed at a time" int no default="16" diff --git a/applications/rtkkatsevichtrash/CMakeLists.txt b/applications/rtkkatsevichtrash/CMakeLists.txt new file mode 100644 index 000000000..69163fafd --- /dev/null +++ b/applications/rtkkatsevichtrash/CMakeLists.txt @@ -0,0 +1,14 @@ +WRAP_GGO(rtkkatsevichtrash_GGO_C rtkkatsevichtrash.ggo ../rtkinputprojections_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo) +add_executable(rtkkatsevichtrash rtkkatsevichtrash.cxx ${rtkkatsevichtrash_GGO_C}) +target_link_libraries(rtkkatsevichtrash RTK) + +# Installation code +if(NOT RTK_INSTALL_NO_EXECUTABLES) + foreach(EXE_NAME rtkkatsevichtrash) + 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() + diff --git a/applications/rtkkatsevichtrash/rtkkatsevichtrash.cxx b/applications/rtkkatsevichtrash/rtkkatsevichtrash.cxx new file mode 100644 index 000000000..d4d892801 --- /dev/null +++ b/applications/rtkkatsevichtrash/rtkkatsevichtrash.cxx @@ -0,0 +1,166 @@ +/*========================================================================= + * + * 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 "rtkkatsevichtrash_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkProgressCommands.h" +#include "rtkFFTHilbertImageFilter.h" +#include + + +#include +#include +#include + +int +main(int argc, char * argv[]) +{ + GGO(rtkkatsevichtrash, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + using PILineRealType = float; + using PILineImageType = itk::Image, Dimension>; + using PILineImagePointer = typename PILineImageType::Pointer; + using PILineImageFilterType = rtk::PILineImageFilter; + using PILinePointer = typename PILineImageFilterType::Pointer; + + + //// Projections reader + using ReaderType = rtk::ProjectionsReader; + ReaderType::Pointer reader = ReaderType::New(); + rtk::SetProjectionsReaderFromGgo(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(); + + OutputImageType::Pointer cst = OutputImageType::New(); + OutputImageType::RegionType region; + OutputImageType::IndexType index; + OutputImageType::SizeType size; + OutputImageType::SpacingType spacing; + OutputImageType::PointType origin; + spacing.Fill(1.); + index.Fill(0); + size[0] = 256; + size[1] = 1; + size[2] = 256; + origin[0] = -0.5 * 255; + origin[1] = 0.; + origin[2] = -0.5 * 255; + region.SetSize(size); + region.SetIndex(index); + cst->SetRegions(region); + cst->SetSpacing(spacing); + cst->SetOrigin(origin); + cst->Allocate(); + cst->FillBuffer(0.); + + PILineImageFilterType::Pointer pil = PILineImageFilterType::New(); + pil->SetInput(cst); + pil->SetGeometry(geometry); + pil->Update(); + + + OutputImageType::Pointer pil1 = OutputImageType::New(); + OutputImageType::Pointer pil2 = OutputImageType::New(); + pil1->SetRegions(region); + pil1->SetSpacing(spacing); + pil1->SetOrigin(origin); + pil1->Allocate(); + pil2->SetRegions(region); + pil2->SetSpacing(spacing); + pil2->SetOrigin(origin); + pil2->Allocate(); + + using PILineIterator = itk::ImageRegionConstIterator; + using OutputRegionIterator = itk::ImageRegionIteratorWithIndex; + + + PILineIterator itPIL(pil->GetOutput(),region); + OutputRegionIterator itPIL1(pil1,region); + OutputRegionIterator itPIL2(pil2,region); + + itPIL.GoToBegin(); + itPIL1.GoToBegin(); + itPIL2.GoToBegin(); + itk::Vector pil_bounds; + + while (!itPIL.IsAtEnd()) + { + pil_bounds = itPIL.Get(); + itPIL1.Set(pil_bounds[0]); + itPIL2.Set(pil_bounds[1]); + ++itPIL; + ++itPIL1; + ++itPIL2; + } + + + // 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 + + + // Write + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output_arg); + writer->SetInput(pil1); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + writer->SetFileName(args_info.input_arg); + writer->SetInput(pil2); + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + + return EXIT_SUCCESS; +} diff --git a/applications/rtkkatsevichtrash/rtkkatsevichtrash.ggo b/applications/rtkkatsevichtrash/rtkkatsevichtrash.ggo new file mode 100644 index 000000000..4ff30fd88 --- /dev/null +++ b/applications/rtkkatsevichtrash/rtkkatsevichtrash.ggo @@ -0,0 +1,10 @@ +package "rtkkatsevichderivative" +purpose "Computes the derivative of a stack of projections wrt to the projection index (rotation angle). See Noo et al., PMB, 2003" + +option "verbose" v "Verbose execution" flag off +option "config" - "Config file" string no +option "geometry" g "XML geometry file name" string yes +option "input" i "Input file name" string yes +option "output" o "Output file name" string yes +option "hardware" - "Hardware used for computation" values="cpu","cuda" no default="cpu" + diff --git a/applications/rtkpilines/CMakeLists.txt b/applications/rtkpilines/CMakeLists.txt new file mode 100644 index 000000000..0e40f3af2 --- /dev/null +++ b/applications/rtkpilines/CMakeLists.txt @@ -0,0 +1,14 @@ +WRAP_GGO(rtkpilines_GGO_C rtkpilines.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo) +add_executable(rtkpilines rtkpilines.cxx ${rtkpilines_GGO_C}) +target_link_libraries(rtkpilines RTK) + +# Installation code +if(NOT RTK_INSTALL_NO_EXECUTABLES) + foreach(EXE_NAME rtkpilines) + 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() + diff --git a/applications/rtkpilines/rtkpilines.cxx b/applications/rtkpilines/rtkpilines.cxx new file mode 100644 index 000000000..35e5933a6 --- /dev/null +++ b/applications/rtkpilines/rtkpilines.cxx @@ -0,0 +1,153 @@ +/*========================================================================= + * + * 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 "rtkpilines_ggo.h" +#include "rtkGgoFunctions.h" +#include "rtkConfiguration.h" + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkPILineImageFilter.h" +#include "rtkProgressCommands.h" + +#include +#include +#include +#include +#include +#include + +int +main(int argc, char * argv[]) +{ + GGO(rtkpilines, args_info); + + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using InputImageType = itk::Image; + // using OutputImageType = itk::VectorImage; + using OutputImageType = itk::Image, Dimension>; + + // Projections reader + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(args_info.input_arg); + + 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 PILineFilterType = rtk::PILineImageFilter; + PILineFilterType::Pointer pilines = PILineFilterType::New(); + pilines->SetGeometry(geometry); + pilines->SetInput(reader->GetOutput()); + pilines->Update(); + + if (!strcmp(args_info.filetype_arg, "2")) + { + + InputImageType::RegionType region = reader->GetOutput()->GetLargestPossibleRegion(); + + OutputImageType::Pointer outImage = pilines->GetOutput(); + + InputImageType::Pointer im1 = InputImageType::New(); + im1->CopyInformation(reader->GetOutput()); + im1->SetRegions(region); + im1->Allocate(); + + InputImageType::Pointer im2 = InputImageType::New(); + im2->CopyInformation(reader->GetOutput()); + im2->SetRegions(region); + im2->Allocate(); + + using RegionIterator = itk::ImageRegionIterator; + using OutRegionIterator = itk::ImageRegionIterator; + + RegionIterator it1(im1, im1->GetLargestPossibleRegion()); + RegionIterator it2(im2, im2->GetLargestPossibleRegion()); + OutRegionIterator out(pilines->GetOutput(), pilines->GetOutput()->GetLargestPossibleRegion()); + + for (it1.GoToBegin(), it2.GoToBegin(), out.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++out) + { + using VectorType = itk::Vector; + VectorType vector = out.Get(); + + it1.Set(vector[0]); + it2.Set(vector[1]); + } + + // Write + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output1_arg); + writer->SetInput(im1); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + + writer->SetFileName(args_info.output2_arg); + writer->SetInput(im2); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + } + else if (!strcmp(args_info.filetype_arg, "1")) + { + // Write + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(args_info.output1_arg); + writer->SetInput(pilines->GetOutput()); + + if (args_info.verbose_flag) + std::cout << "Reconstructing and writing... " << std::endl; + + TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update()) + } + + else + { + std::cerr << "The filetype option was not properly set." << std::endl; + std::cerr << "1 : one vector image file. 2 : 2 single valued images." << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/applications/rtkpilines/rtkpilines.ggo b/applications/rtkpilines/rtkpilines.ggo new file mode 100644 index 000000000..cd666cf67 --- /dev/null +++ b/applications/rtkpilines/rtkpilines.ggo @@ -0,0 +1,12 @@ +package "rtkpilines" +purpose "Compute the PI-lines intervalle of a 3D volume. (see Noo et al., PMB, 2003). The output is a 3D 2D-vector image" + +option "verbose" v "Verbose execution" flag off +option "config" - "Config file" string no +option "geometry" g "XML geometry file name" string yes +option "input" i "Input file name" string yes +option "output1" o "Output file name" string yes +option "output2" p "Output file name" string yes +option "filetype" - "Type of output file. 1 : one single vector image. 2 : 2 sep files" string no default="2" +option "hardware" - "Hardware used for computation" values="cpu","cuda" no default="cpu" + diff --git a/applications/rtksimulatedgeometry/rtksimulatedgeometry.cxx b/applications/rtksimulatedgeometry/rtksimulatedgeometry.cxx index b6e5004f1..ae65dbf2f 100644 --- a/applications/rtksimulatedgeometry/rtksimulatedgeometry.cxx +++ b/applications/rtksimulatedgeometry/rtksimulatedgeometry.cxx @@ -49,6 +49,7 @@ main(int argc, char * argv[]) if (args_info.rad_cyl_given) geometry->SetRadiusCylindricalDetector(args_info.rad_cyl_arg); + // Write rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New(); diff --git a/include/rtkFFTHilbertImageFilter.h b/include/rtkFFTHilbertImageFilter.h new file mode 100644 index 000000000..5522e50c4 --- /dev/null +++ b/include/rtkFFTHilbertImageFilter.h @@ -0,0 +1,112 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkFFTHilbertImageFilter_h +#define rtkFFTHilbertImageFilter_h + +#include +#include "rtkConfiguration.h" +#include "rtkFFTProjectionsConvolutionImageFilter.h" +#include "rtkMacro.h" + +namespace rtk +{ + +/** \class FFTHilbertImageFilter + * \brief Implements the Hilbert transform. + * + * \test rtkhilbertfiltertest.cxx + * + * \author Aurélien Coussat + * + * \ingroup RTK ImageToImageFilter + */ + +template +class ITK_EXPORT FFTHilbertImageFilter + : public rtk::FFTProjectionsConvolutionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(FFTHilbertImageFilter); + + /** Standard class type alias. */ + using Self = FFTHilbertImageFilter; + using Superclass = rtk::FFTProjectionsConvolutionImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Some convenient type alias. */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using FFTPrecisionType = TFFTPrecision; + using IndexType = typename InputImageType::IndexType; + using SizeType = typename InputImageType::SizeType; + + using FFTInputImageType = typename Superclass::FFTInputImageType; + using FFTInputImagePointer = typename FFTInputImageType::Pointer; + using FFTOutputImageType = typename Superclass::FFTOutputImageType; + using FFTOutputImagePointer = typename FFTOutputImageType::Pointer; + + /** Standard New method. */ + itkNewMacro(Self); + + itkGetMacro(PixelShift, double); + // The Set macro is redefined to clear the current FFT kernel when a parameter + // is modified. + virtual void + SetPixelShift(const double _arg) + { + itkDebugMacro("setting PixelShift to " << _arg); + CLANG_PRAGMA_PUSH + CLANG_SUPPRESS_Wfloat_equal if (this->m_PixelShift != _arg) + { + this->m_PixelShift = _arg; + this->Modified(); + this->m_KernelFFT = nullptr; + } + CLANG_PRAGMA_POP + } + + /** Runtime information support. */ + itkTypeMacro(FFTHilbertImageFilter, FFTProjectionsConvolutionImageFilter); + +protected: + FFTHilbertImageFilter() = default; + ~FFTHilbertImageFilter() override = default; + + void + GenerateOutputInformation() override; + + /** Create and return a pointer to one line of the Hilbert kernel in Fourier space. + * Used in generate data functions. */ + void + UpdateFFTProjectionsConvolutionKernel(const SizeType s) override; + +private: + SizeType m_PreviousKernelUpdateSize; + double m_PixelShift; + +}; // end of class + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkFFTHilbertImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkFFTHilbertImageFilter.hxx b/include/rtkFFTHilbertImageFilter.hxx new file mode 100644 index 000000000..5477b2895 --- /dev/null +++ b/include/rtkFFTHilbertImageFilter.hxx @@ -0,0 +1,101 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkFFTHilbertImageFilter_hxx +#define rtkFFTHilbertImageFilter_hxx + +#include "rtkFFTHilbertImageFilter.h" +#include "itkForwardFFTImageFilter.h" + +namespace rtk +{ + +template +void +FFTHilbertImageFilter::GenerateOutputInformation() +{ + FFTProjectionsConvolutionImageFilter::GenerateOutputInformation(); + + auto origin = this->GetOutput()->GetOrigin(); + double spacing = this->GetOutput()->GetSpacing()[0]; + origin[0] += m_PixelShift * spacing; + this->GetOutput()->SetOrigin(origin); +} + +template +void +FFTHilbertImageFilter::UpdateFFTProjectionsConvolutionKernel(const SizeType s) +{ + if (this->m_KernelFFT.GetPointer() != nullptr && s == this->m_PreviousKernelUpdateSize) + { + return; + } + m_PreviousKernelUpdateSize = s; + + const int width = s[0]; + const int height = s[1]; + + // Allocate kernel + SizeType size; + size.Fill(1); + size[0] = width; + FFTInputImagePointer kernel = FFTInputImageType::New(); + kernel->SetRegions(size); + kernel->Allocate(); + kernel->FillBuffer(0.); + + itk::ImageRegionIterator it(kernel, kernel->GetLargestPossibleRegion()); + + // Compute band-limited kernel in space domain + double spacing = this->GetInput()->GetSpacing()[0]; + IndexType ix; + + while (!it.IsAtEnd()) + { + ix = it.GetIndex(); + + double x = (ix[0] + m_PixelShift) * spacing; + if (x > (size[0] / 2) * spacing) + x -= size[0] * spacing; + + if (x == 0.) + { + it.Set(0.); + } + else + { + double v = spacing * (1. - std::cos(itk::Math::pi * x / spacing)) / (itk::Math::pi * x); + it.Set(v); + } + + ++it; + } + + // FFT kernel + using FFTType = itk::RealToHalfHermitianForwardFFTImageFilter; + typename FFTType::Pointer fftK = FFTType::New(); + fftK->SetInput(kernel); + fftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK->Update(); + + this->m_KernelFFT = fftK->GetOutput(); + this->m_KernelFFT->DisconnectPipeline(); +} + +} // end namespace rtk +#endif diff --git a/include/rtkHilbertTransformOnKappaLinesImageFilter.h b/include/rtkHilbertTransformOnKappaLinesImageFilter.h new file mode 100644 index 000000000..484eda2bd --- /dev/null +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.h @@ -0,0 +1,108 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkHilbertTransformOnKappaLinesImageFilter_h +#define rtkHilbertTransformOnKappaLinesImageFilter_h + +#include "rtkKatsevichForwardBinningImageFilter.h" +#include "rtkFFTHilbertImageFilter.h" +#include "rtkKatsevichBackwardBinningImageFilter.h" + +#include + +#include + + +namespace rtk +{ + +/** \class katsHilbertTransformOnKappaLinesImageFilter + * \brief Implements the Hilbert transform as described in Noo et al., PMB, 2003. + * It is 3-step : + * - Resample the data in v on the psi angle values + * - Apply 1D Hilbert transform on the resampled psi-direction. + * - Resample the data back to the usual vertical coordinates. + * + * \author Jerome Lesaint and Alexandre Esa + * + * \test katsHilbertTransformOnKappaLines.cxx, katsKatsevitchHelicalReconstructionImageFilterTets.cxx. + * + * \ingroup + */ +template +class ITK_EXPORT HilbertTransformOnKappaLinesImageFilter : public itk::ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(HilbertTransformOnKappaLinesImageFilter); + + using Self = HilbertTransformOnKappaLinesImageFilter; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + + /** Method for creation through object factory */ + itkNewMacro(Self); + + /** Run-time type information */ + itkTypeMacro(HilbertTransformOnKappaLinesImageFilter, ImageToImageFilter); + + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + + using GeometryType = rtk::ThreeDHelicalProjectionGeometry; + using GeometryPointer = GeometryType::Pointer; + + itkGetMacro(Geometry, GeometryPointer); + itkSetMacro(Geometry, GeometryPointer); + +protected: + HilbertTransformOnKappaLinesImageFilter(); + ~HilbertTransformOnKappaLinesImageFilter() override = default; + +protected: + using ForwardBinningType = KatsevichForwardBinningImageFilter; + using FFTHilbertType = FFTHilbertImageFilter; + using BackwardBinningType = KatsevichBackwardBinningImageFilter; + + void + GenerateData() override; + + void + GenerateOutputInformation() override; + + /** Display */ + void + PrintSelf(std::ostream & os, itk::Indent indent) const override; + +private: + typename ForwardBinningType::Pointer m_ForwardFilter; + typename FFTHilbertType::Pointer m_HilbertFilter; + typename BackwardBinningType::Pointer m_BackwardFilter; + + GeometryPointer m_Geometry; + +}; // end of class + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkHilbertTransformOnKappaLinesImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx b/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx new file mode 100644 index 000000000..efd43669a --- /dev/null +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx @@ -0,0 +1,76 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkHilbertTransformOnKappaLinesImageFilter_hxx +#define rtkHilbertTransformOnKappaLinesImageFilter_hxx + +#include "rtkHilbertTransformOnKappaLinesImageFilter.h" + + +namespace rtk +{ + +template +void +HilbertTransformOnKappaLinesImageFilter::GenerateOutputInformation() +{ + // Mettre a jour la taille de l'image de sortie. Ainsi que l'origine. + // Sur le modele de itk::ShrinkImageFilter + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); + + // std::cout << "(Hilbert transf) Input Size : " << this->GetInput()->GetLargestPossibleRegion().GetSize(); + // std::cout << " Inpu torigin : " << this->GetInput()->GetOrigin() << std::endl; +} + +template +HilbertTransformOnKappaLinesImageFilter::HilbertTransformOnKappaLinesImageFilter() +{ + m_ForwardFilter = ForwardBinningType::New(); + m_HilbertFilter = FFTHilbertType::New(); + m_HilbertFilter->SetInput(m_ForwardFilter->GetOutput()); + m_BackwardFilter = BackwardBinningType::New(); + m_BackwardFilter->SetInput(m_HilbertFilter->GetOutput()); +} + +template +void +HilbertTransformOnKappaLinesImageFilter::GenerateData() +{ + typename InputImageType::Pointer input = InputImageType::New(); + input->Graft(const_cast(this->GetInput())); + m_ForwardFilter->SetInput(input); + m_ForwardFilter->SetGeometry(this->m_Geometry); + + m_BackwardFilter->SetGeometry(this->m_Geometry); + m_BackwardFilter->GraftOutput(this->GetOutput()); + m_BackwardFilter->Update(); + this->GraftOutput(m_BackwardFilter->GetOutput()); +} + +template +void +HilbertTransformOnKappaLinesImageFilter::PrintSelf(std::ostream & os, + itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "This is a Hilbert Image Filter" << std::endl; +} + +} // namespace rtk +#endif diff --git a/include/rtkKatsevichBackProjectionImageFilter.h b/include/rtkKatsevichBackProjectionImageFilter.h new file mode 100644 index 000000000..6e14b5b91 --- /dev/null +++ b/include/rtkKatsevichBackProjectionImageFilter.h @@ -0,0 +1,190 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichBackProjectionImageFilter_h +#define rtkKatsevichBackProjectionImageFilter_h + +#include + +#include +#include + +#include +#include + +#include +#include + +namespace rtk +{ + +/** \class KatsevichBackProjectionImageFilter + * \brief 3D backprojection + * + * Backprojects a stack of projection images (input 1) in a 3D volume (input 0) + * using linear interpolation according to a specified geometry. The operation + * is voxel-based, meaning that the center of each voxel is projected in the + * projection images to determine the interpolation location. + * + * \test + * + * \author Simon Rit, Jerome Lesaint + * + * \ingroup RTK Projector + */ +template +class KatsevichBackProjectionImageFilter : public itk::InPlaceImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichBackProjectionImageFilter); + + /** Standard class type alias. */ + using Self = KatsevichBackProjectionImageFilter; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Convenient type alias. */ + using InputPixelType = typename TInputImage::PixelType; + using InternalInputPixelType = typename TInputImage::InternalPixelType; + using OutputImageRegionType = typename TOutputImage::RegionType; + using GeometryType = rtk::ThreeDHelicalProjectionGeometry; + using GeometryPointer = typename GeometryType::Pointer; + using ProjectionMatrixType = typename GeometryType::MatrixType; + using ProjectionImageType = itk::Image; + using ProjectionImagePointer = typename ProjectionImageType::Pointer; + using PILineRealType = float; + using PILineImageType = itk::Image, TInputImage::ImageDimension>; + using PILineImagePointer = typename PILineImageType::Pointer; + using PILineImageFilterType = PILineImageFilter; + using PILinePointer = typename PILineImageFilterType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(KatsevichBackProjectionImageFilter, itk::ImageToImageFilter); + + // virtual void + // SetPILines(PILineImagePointer image) ; + // + // virtual const PILineImagePointer + // GetPILines() const; + + + /** Get / Set the object pointer to projection geometry */ + itkGetModifiableObjectMacro(Geometry, GeometryType); + itkSetObjectMacro(Geometry, GeometryType); + + ///** Get / Set the object pointer to PILines */ + // itkGetMacro(PILines, PILinePointer); + // itkSetMacro(PILines, PILinePointer); + + /** Get / Set the transpose flag for 2D projections (optimization trick) */ + itkGetMacro(Transpose, bool); + itkSetMacro(Transpose, bool); + +protected: + KatsevichBackProjectionImageFilter(); + ~KatsevichBackProjectionImageFilter() override = default; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + /** Apply changes to the input image requested region. */ + void + GenerateInputRequestedRegion() override; + + void + BeforeThreadedGenerateData() override; + + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + + /** Special case when the detector is cylindrical and centered on source */ + virtual void + CylindricalDetectorCenteredOnSourceBackprojection( + const OutputImageRegionType & region, + const ProjectionMatrixType & volIndexToProjPP, + const itk::Matrix & projPPToProjIndex, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda); + + /** Optimized version when the rotation is parallel to X, i.e. matrix[1][0] + and matrix[2][0] are zeros. */ + virtual void + OptimizedBackprojectionX(const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda); + + /** Optimized version when the rotation is parallel to Y, i.e. matrix[1][1] + and matrix[2][1] are zeros. */ + virtual void + OptimizedBackprojectionY(const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda); + + /** The two inputs should not be in the same space so there is nothing + * to verify. */ + void + VerifyInputInformation() const override + {} + + /** The input is a stack of projections, we need to interpolate in one projection + for efficiency during interpolation. Use of itk::ExtractImageFilter is + not threadsafe in ThreadedGenerateData, this one is. The output can be multiplied by a constant. + The function is templated to allow getting an itk::CudaImage. */ + template + typename TProjectionImage::Pointer + GetProjection(const unsigned int iProj); + + /** Creates iProj index to index projection matrices with current inputs + instead of the physical point to physical point projection matrix provided by Geometry */ + ProjectionMatrixType + GetIndexToIndexProjectionMatrix(const unsigned int iProj); + + ProjectionMatrixType + GetVolumeIndexToProjectionPhysicalPointMatrix(const unsigned int iProj); + + itk::Matrix + GetProjectionPhysicalPointToProjectionIndexMatrix(const unsigned int iProj); + + /** RTK geometry object */ + GeometryPointer m_Geometry; + PILineImagePointer m_PILines; + + +private: + /** Flip projection flag: infludences GetProjection and + GetIndexToIndexProjectionMatrix for optimization */ + bool m_Transpose{ false }; +}; + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkKatsevichBackProjectionImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkKatsevichBackProjectionImageFilter.hxx b/include/rtkKatsevichBackProjectionImageFilter.hxx new file mode 100644 index 000000000..361b68147 --- /dev/null +++ b/include/rtkKatsevichBackProjectionImageFilter.hxx @@ -0,0 +1,840 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichBackProjectionImageFilter_hxx +#define rtkKatsevichBackProjectionImageFilter_hxx + +#include "math.h" + +#include "rtkKatsevichBackProjectionImageFilter.h" + +#include + +#include +#include +#include +#include + +namespace rtk +{ + + +template +KatsevichBackProjectionImageFilter::KatsevichBackProjectionImageFilter() +{ + m_Geometry = nullptr; + m_PILines = nullptr; + this->SetNumberOfRequiredInputs(2); + this->SetInPlace(true); +} + + +// template +// void +// KatsevichBackProjectionImageFilter::SetPILines(PILineImageType * image) +//{ +// // Process object is not const-correct so the const_cast is required here +// this->ProcessObject::SetNthInput(1, const_cast(image)); +//} +// +// template +// const typename KatsevichBackProjectionImageFilter::PILineImageType * +// KatsevichBackProjectionImageFilter::GetPILines() const +//{ +// return itkDynamicCastInDebugMode(m_PILines); +//} + + +template +void +KatsevichBackProjectionImageFilter::GenerateInputRequestedRegion() +{ + // Input 0 is the volume in which we backproject + typename Superclass::InputImagePointer inputPtr0 = const_cast(this->GetInput(0)); + if (!inputPtr0) + return; + inputPtr0->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); + + // Input 1 is the stack of projections to backproject + typename Superclass::InputImagePointer inputPtr1 = const_cast(this->GetInput(1)); + if (!inputPtr1) + return; + + // Geometry size check + const unsigned int Dimension = TInputImage::ImageDimension; + const int lastProjIndex = this->GetInput(1)->GetLargestPossibleRegion().GetIndex(Dimension - 1) + + this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension - 1); + if ((int)this->m_Geometry->GetMatrices().size() < lastProjIndex) + { + itkExceptionMacro(<< "Mismatch between the number of projections and the geometry entries. " + << "Geometry has " << this->m_Geometry->GetMatrices().size() + << " entries, which is less than the " + << "last index of the projections stack, i.e., " << lastProjIndex << "."); + } + + typename TInputImage::RegionType reqRegion = inputPtr1->GetLargestPossibleRegion(); + if (m_Geometry.GetPointer() == nullptr || m_Geometry->GetRadiusCylindricalDetector() != 0) + { + inputPtr1->SetRequestedRegion(inputPtr1->GetLargestPossibleRegion()); + return; + } + + itk::ContinuousIndex cornerInf; + itk::ContinuousIndex cornerSup; + cornerInf[0] = itk::NumericTraits::max(); + cornerSup[0] = itk::NumericTraits::NonpositiveMin(); + cornerInf[1] = itk::NumericTraits::max(); + cornerSup[1] = itk::NumericTraits::NonpositiveMin(); + cornerInf[2] = reqRegion.GetIndex(2); + cornerSup[2] = reqRegion.GetIndex(2) + reqRegion.GetSize(2); + + // Go over each projection + const unsigned int nProj = this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension - 1); + const unsigned int iFirstProj = this->GetInput(1)->GetLargestPossibleRegion().GetIndex(Dimension - 1); + this->SetTranspose(false); + for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) + { + // Extract the current slice + ProjectionMatrixType matrix = GetIndexToIndexProjectionMatrix(iProj); + + // Check which part of the projection image will be backprojected in the + // volume. + double firstPerspFactor = 0.; + int c[Dimension] = { 0 }; + for (c[2] = 0; c[2] < 2; c[2]++) + for (c[1] = 0; c[1] < 2; c[1]++) + for (c[0] = 0; c[0] < 2; c[0]++) + { + // Compute corner index + const double eps = 1e-4; + itk::ContinuousIndex index; + for (unsigned int i = 0; i < Dimension; i++) + { + index[i] = this->GetOutput()->GetRequestedRegion().GetIndex()[i] - eps; + index[i] += c[i] * (this->GetOutput()->GetRequestedRegion().GetSize(i) - 1 + 2 * eps); + } + + // Compute projection index + itk::ContinuousIndex point; + for (unsigned int i = 0; i < Dimension - 1; i++) + { + point[i] = matrix[i][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + point[i] += matrix[i][j] * index[j]; + } + + // Apply perspective + double perspFactor = matrix[Dimension - 1][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + perspFactor += matrix[Dimension - 1][j] * index[j]; + perspFactor = 1 / perspFactor; + for (unsigned int i = 0; i < Dimension - 1; i++) + point[i] = point[i] * perspFactor; + + // Check if corners all have the same perspective factor sign. + // If not, source is too close for easily computing a smaller requested + // region than the largest possible one. + if (c[0] + c[1] + c[2] == 0) + firstPerspFactor = perspFactor; + if (perspFactor * firstPerspFactor < 0.) // Change of sign + { + inputPtr1->SetRequestedRegion(inputPtr1->GetLargestPossibleRegion()); + return; + } + // Look for extremas on projection to calculate requested region + for (int i = 0; i < 2; i++) + { + cornerInf[i] = std::min(cornerInf[i], point[i]); + cornerSup[i] = std::max(cornerSup[i], point[i]); + } + } + } + + + // std::cout << "Corner inf : " << cornerInf << "Corner sup " << cornerSup << std::endl; + reqRegion.SetIndex(0, itk::Math::floor(cornerInf[0])); + reqRegion.SetIndex(1, itk::Math::floor(cornerInf[1])); + reqRegion.SetSize(0, itk::Math::ceil(cornerSup[0] - reqRegion.GetIndex(0)) + 1); + reqRegion.SetSize(1, itk::Math::ceil(cornerSup[1] - reqRegion.GetIndex(1)) + 1); + + if (reqRegion.Crop(inputPtr1->GetLargestPossibleRegion())) + { + inputPtr1->SetRequestedRegion(reqRegion); + } + else + { + inputPtr1->SetRequestedRegion(inputPtr1->GetLargestPossibleRegion()); + } +} + +template +void +KatsevichBackProjectionImageFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull() || !this->m_Geometry->GetTheGeometryIsVerified()) + itkExceptionMacro(<< "Geometry has not been set or not been checked"); +} + + +template +void +KatsevichBackProjectionImageFilter::BeforeThreadedGenerateData() +{ + this->SetTranspose(true); + + // Check if detector is cylindrical + double radius = m_Geometry->GetRadiusCylindricalDetector(); + if ((radius != 0) && (radius != this->m_Geometry->GetSourceToDetectorDistances()[0])) + { + itkGenericExceptionMacro(<< "Voxel-based back projector can currently handle a cylindrical detector only when it " + "is centered on the source. " + << "Detector radius is " << radius << ", should be " + << this->m_Geometry->GetSourceToDetectorDistances()[0]); + } + + typename PILineImageFilterType::Pointer pil = PILineImageFilterType::New(); + pil->SetInput(this->GetOutput()); + pil->SetGeometry(this->m_Geometry); + pil->Update(); + this->m_PILines = pil->GetOutput(); +} + +/** + * GenerateData performs the accumulation + */ +template +void +KatsevichBackProjectionImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + const unsigned int nProj = this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension - 1); + const unsigned int iFirstProj = this->GetInput(1)->GetLargestPossibleRegion().GetIndex(Dimension - 1); + + // Create interpolator, could be any interpolation + using InterpolatorType = itk::LinearInterpolateImageFunction; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + + // Iterators on volume input and output + using InputRegionIterator = itk::ImageRegionConstIterator; + using PILineIterator = itk::ImageRegionConstIterator; + using OutputRegionIterator = itk::ImageRegionIteratorWithIndex; + + InputRegionIterator itIn(this->GetInput(0), outputRegionForThread); + PILineIterator itPIL(this->m_PILines, outputRegionForThread); + OutputRegionIterator itOut(this->GetOutput(), outputRegionForThread); + + // Initialize output region with input region in case the filter is not in + // place + if (this->GetInput() != this->GetOutput()) + { + itIn.GoToBegin(); + while (!itIn.IsAtEnd()) + { + itOut.Set(itIn.Get()); + ++itIn; + ++itOut; + } + } + + // Continuous index at which we interpolate + itk::ContinuousIndex pointProj; + + // Get First Angle Values to map idx_proj -> angles + double firstAngle = m_Geometry->GetHelicalAngles()[0]; + double delta_lambda = m_Geometry->GetHelixAngularGap(); + double lambda = 0.; + + // Go over each projections + for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) + { + // Extract the current slice + ProjectionImagePointer projection = GetProjection(iProj); + + ProjectionMatrixType matrix = this->GetIndexToIndexProjectionMatrix(iProj); + interpolator->SetInputImage(projection); + + lambda = firstAngle + iProj * delta_lambda; + + // Cylindrical detector centered on source case + if (m_Geometry->GetRadiusCylindricalDetector() != 0) + { + ProjectionMatrixType volIndexToProjPP = GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); + itk::Matrix projPPToProjIndex = + GetProjectionPhysicalPointToProjectionIndexMatrix(iProj); + CylindricalDetectorCenteredOnSourceBackprojection( + outputRegionForThread, volIndexToProjPP, projPPToProjIndex, projection, lambda, delta_lambda); + continue; + } + + // Optimized version + if (fabs(matrix[1][0]) < 1e-10 && fabs(matrix[2][0]) < 1e-10) + { + OptimizedBackprojectionX(outputRegionForThread, matrix, projection, lambda, delta_lambda); + continue; + } + if (fabs(matrix[1][1]) < 1e-10 && fabs(matrix[2][1]) < 1e-10) + { + OptimizedBackprojectionY(outputRegionForThread, matrix, projection, lambda, delta_lambda); + continue; + } + + // Go over each voxel + itOut.GoToBegin(); + itPIL.GoToBegin(); + while (!itOut.IsAtEnd()) + { + // Compute projection index + for (unsigned int i = 0; i < Dimension - 1; i++) + { + pointProj[i] = matrix[i][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + pointProj[i] += matrix[i][j] * itOut.GetIndex()[j]; + } + + // Apply perspective + double perspFactor = matrix[Dimension - 1][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + perspFactor += matrix[Dimension - 1][j] * itOut.GetIndex()[j]; + perspFactor = 1 / perspFactor; + for (unsigned int i = 0; i < Dimension - 1; i++) + pointProj[i] = pointProj[i] * perspFactor; + + // Interpolate if in projection + if (interpolator->IsInsideBuffer(pointProj)) + { + + itk::Vector piline = itPIL.Get(); + typename TOutputImage::InternalPixelType sb = piline[0]; + typename TOutputImage::InternalPixelType st = piline[1]; + + double rho = 0; + double d_in = (lambda - sb) / delta_lambda; + double d_out = (st - lambda) / delta_lambda; + // Apply PiLineConditions + + //if (lambda < sb - delta_lambda) + // rho = 0; + //else if (lambda < sb) + // rho = 0.5 * (d_in + 1) * (d_in + 1); + //else if (lambda < sb + delta_lambda) + // rho = 0.5 + d_in + 0.5 * d_in * d_in; + //else if (lambda < st - delta_lambda) + // rho = 1; + //else if (lambda < st) + // rho = 0.5 + d_out + 0.5 * d_out * d_out; + //else if (lambda < st + delta_lambda) + // rho = 0.5 * (1 + d_out * (1 + d_out)); + //else + // rho = 0; + + if (lambda < sb) + rho = 0.; + else if (lambda > sb && lambda < st) + rho = 1.; + else + rho = 0.; + + // Compute normlazation coefficient wstar + typename TOutputImage::PointType outPoint; + this->GetOutput()->TransformIndexToPhysicalPoint(itOut.GetIndex(), outPoint); + double sid = this->GetGeometry()->GetHelixRadius(); + double wstar = sid - outPoint[0] * sin(lambda) - outPoint[2] * cos(lambda); + + itOut.Set(itOut.Get() + + rho * delta_lambda * interpolator->EvaluateAtContinuousIndex(pointProj) / (wstar * 2 * M_PI)); + } + + ++itOut; + ++itPIL; + } + } +} + + +template +void +KatsevichBackProjectionImageFilter::CylindricalDetectorCenteredOnSourceBackprojection( + const OutputImageRegionType & region, + const ProjectionMatrixType & volIndexToProjPP, + const itk::Matrix & projPPToProjIndex, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda) +{ + using OutputRegionIterator = itk::ImageRegionIteratorWithIndex; + OutputRegionIterator itOut(this->GetOutput(), region); + + using PILineIterator = itk::ImageRegionConstIterator; + PILineIterator itPIL(this->m_PILines, region); + + + const unsigned int Dimension = TInputImage::ImageDimension; + + // Create interpolator, could be any interpolation + using InterpolatorType = itk::LinearInterpolateImageFunction; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + interpolator->SetInputImage(projection); + + // Get radius of the cylindrical detector + double radius = m_Geometry->GetRadiusCylindricalDetector(); + + // Continuous index at which we interpolate + itk::ContinuousIndex pointProj, pointProjIdx; + + + // Go over each voxel + itOut.GoToBegin(); + itPIL.GoToBegin(); + while (!itOut.IsAtEnd()) + { + // Compute projection index + for (unsigned int i = 0; i < Dimension - 1; i++) + { + pointProj[i] = volIndexToProjPP[i][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + pointProj[i] += volIndexToProjPP[i][j] * itOut.GetIndex()[j]; + } + + // Apply perspective + double perspFactor = volIndexToProjPP[Dimension - 1][Dimension]; + for (unsigned int j = 0; j < Dimension; j++) + perspFactor += volIndexToProjPP[Dimension - 1][j] * itOut.GetIndex()[j]; + perspFactor = 1 / perspFactor; + for (unsigned int i = 0; i < Dimension - 1; i++) + pointProj[i] = pointProj[i] * perspFactor; + + // Apply correction for cylindrical centered on source + const double u = pointProj[0]; + pointProj[0] = radius * atan2(u, radius); + pointProj[1] = pointProj[1] * radius / sqrt(radius * radius + u * u); + + // Convert to projection index + for (unsigned int i = 0; i < Dimension - 1; i++) + { + pointProjIdx[i] = projPPToProjIndex[i][Dimension - 1]; + for (unsigned int j = 0; j < Dimension - 1; j++) + pointProjIdx[i] += projPPToProjIndex[i][j] * pointProj[j]; + } + + // Interpolate if in projection + if (interpolator->IsInsideBuffer(pointProjIdx)) + { + typename TOutputImage::PixelType v = interpolator->EvaluateAtContinuousIndex(pointProjIdx); + itk::Vector piline = itPIL.Get(); + typename TOutputImage::InternalPixelType sb = piline[0]; + typename TOutputImage::InternalPixelType st = piline[1]; + + double rho = 0; + double d_in = (lambda - sb) / delta_lambda; + double d_out = (st - lambda) / delta_lambda; + // Apply PiLineConditions + + //if (lambda < sb - delta_lambda) + // rho = 0; + //else if (lambda < sb) + // rho = 0.5 * (d_in + 1) * (d_in + 1); + //else if (lambda < sb + delta_lambda) + // rho = 0.5 + d_in + 0.5 * d_in * d_in; + //else if (lambda < st - delta_lambda) + // rho = 1; + //else if (lambda < st) + // rho = 0.5 + d_out + 0.5 * d_out * d_out; + //else if (lambda < st + delta_lambda) + // rho = 0.5 * (1 + d_out * (1 + d_out)); + //else + // rho = 0; + if (lambda < sb) + rho = 0.; + else if (lambda > sb && lambda < st) + rho = 1.; + else + rho = 0.; + + // Compute normlazation coefficient wstar + typename TOutputImage::PointType outPoint; + this->GetOutput()->TransformIndexToPhysicalPoint(itOut.GetIndex(), outPoint); + double sid = this->GetGeometry()->GetHelixRadius(); + double wstar = sid - outPoint[0] * sin(lambda) - outPoint[2] * cos(lambda); + + itOut.Set(itOut.Get() + + rho * delta_lambda * interpolator->EvaluateAtContinuousIndex(pointProj) / (wstar * 2 * M_PI)); + } + + ++itOut; + ++itPIL; + } +} + + +template +void +KatsevichBackProjectionImageFilter::OptimizedBackprojectionX( + const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda) +{ + typename ProjectionImageType::SizeType pSize = projection->GetBufferedRegion().GetSize(); + typename ProjectionImageType::IndexType pIndex = projection->GetBufferedRegion().GetIndex(); + typename TOutputImage::SizeType vBufferSize = this->GetOutput()->GetBufferedRegion().GetSize(); + typename TOutputImage::IndexType vBufferIndex = this->GetOutput()->GetBufferedRegion().GetIndex(); + typename TInputImage::InternalPixelType * pProj = nullptr; + typename TOutputImage::InternalPixelType *pVol = nullptr, *pVolZeroPointer = nullptr; + typename TOutputImage::IndexType indexOutput; + typename TOutputImage::PointType outPoint; + double sid = this->GetGeometry()->GetHelixRadius(); + + + // Pointers in memory to index (0,0,0) which do not necessarily exist + pVolZeroPointer = this->GetOutput()->GetBufferPointer(); + pVolZeroPointer -= vBufferIndex[0] + vBufferSize[0] * (vBufferIndex[1] + vBufferSize[1] * vBufferIndex[2]); + + using PILineIterator = itk::ImageRegionConstIterator; + PILineIterator itPIL(this->m_PILines, region); + itPIL.GoToBegin(); + + + // Continuous index at which we interpolate + double u = NAN, v = NAN, w = NAN; + int ui = 0, vi = 0; + double du = NAN; + + for (int k = region.GetIndex(2); k < region.GetIndex(2) + (int)region.GetSize(2); k++) + { + for (int j = region.GetIndex(1); j < region.GetIndex(1) + (int)region.GetSize(1); j++) + { + int i = region.GetIndex(0); + u = matrix[0][0] * i + matrix[0][1] * j + matrix[0][2] * k + matrix[0][3]; + v = matrix[1][1] * j + matrix[1][2] * k + matrix[1][3]; + w = matrix[2][1] * j + matrix[2][2] * k + matrix[2][3]; + + // Apply perspective + w = 1 / w; + u = u * w - pIndex[0]; + v = v * w - pIndex[1]; + du = w * matrix[0][0]; + + using ComponentType = typename itk::PixelTraits::ValueType; + ComponentType u1, u2, v1, v2; + vi = itk::Math::floor(v); + if (vi >= 0 && vi < (int)pSize[1] - 1) + { + v1 = v - vi; + v2 = 1.0 - v1; + + pProj = projection->GetBufferPointer() + vi * pSize[0]; + pVol = pVolZeroPointer + i + vBufferSize[0] * (j + k * vBufferSize[1]); + + // Innermost loop + for (; i < (region.GetIndex(0) + (int)region.GetSize(0)); i++, u += du, pVol++) + { + ui = itk::Math::floor(u); + if (ui >= 0 && ui < (int)pSize[0] - 1) + { + itk::Vector piline = itPIL.Get(); + typename TOutputImage::InternalPixelType sb = piline[0]; + typename TOutputImage::InternalPixelType st = piline[1]; + + double rho = 0; + double d_in = (lambda - sb) / delta_lambda; + double d_out = (st - lambda) / delta_lambda; + // Apply PiLineConditions + //if (lambda < sb - delta_lambda) + // rho = 0; + //else if (lambda < sb) + // rho = 0.5 * (d_in + 1) * (d_in + 1); + //else if (lambda < sb + delta_lambda) + // rho = 0.5 + d_in + 0.5 * d_in * d_in; + //else if (lambda < st - delta_lambda) + // rho = 1; + //else if (lambda < st) + // rho = 0.5 + d_out + 0.5 * d_out * d_out; + //else if (lambda < st + delta_lambda) + // rho = 0.5 * (1 + d_out * (1 + d_out)); + //else + // rho = 0; + + if (lambda < sb) + rho = 0.; + else if (lambda > sb && lambda < st) + rho = 1.; + else + rho = 0.; + + + this->GetOutput()->TransformIndexToPhysicalPoint(indexOutput, outPoint); + double wstar = sid - outPoint[0] * sin(lambda) - outPoint[2] * cos(lambda); + u1 = u - ui; + u2 = 1.0 - u1; + *pVol += delta_lambda * rho * + (v2 * (u2 * *(pProj + ui) + u1 * *(pProj + ui + 1)) + + v1 * (u2 * *(pProj + ui + pSize[0]) + u1 * *(pProj + ui + pSize[0] + 1))) / + (wstar * 2 * M_PI); + } + ++itPIL; + } // i + } + } // j + } // k +} + +template +void +KatsevichBackProjectionImageFilter::OptimizedBackprojectionY( + const OutputImageRegionType & region, + const ProjectionMatrixType & matrix, + const ProjectionImagePointer projection, + const double lambda, + const double delta_lambda) + +{ + typename ProjectionImageType::SizeType pSize = projection->GetBufferedRegion().GetSize(); + typename ProjectionImageType::IndexType pIndex = projection->GetBufferedRegion().GetIndex(); + typename TOutputImage::SizeType vBufferSize = this->GetOutput()->GetBufferedRegion().GetSize(); + typename TOutputImage::IndexType vBufferIndex = this->GetOutput()->GetBufferedRegion().GetIndex(); + typename TInputImage::InternalPixelType * pProj = nullptr; + typename TOutputImage::InternalPixelType *pVol = nullptr, *pVolZeroPointer = nullptr; + typename TOutputImage::IndexType indexOutput; + typename TOutputImage::PointType outPoint; + double sid = this->GetGeometry()->GetHelixRadius(); + + + // Pointers in memory to index (0,0,0) which do not necessarily exist + pVolZeroPointer = this->GetOutput()->GetBufferPointer(); + pVolZeroPointer -= vBufferIndex[0] + vBufferSize[0] * (vBufferIndex[1] + vBufferSize[1] * vBufferIndex[2]); + + using PILineIterator = itk::ImageRegionConstIterator; + PILineIterator itPIL(this->m_PILines, region); + itPIL.GoToBegin(); + + // Continuous index at which we interpolate + double u = NAN, v = NAN, w = NAN; + int ui = 0, vi = 0; + double du = NAN; + + for (int k = region.GetIndex(2); k < region.GetIndex(2) + (int)region.GetSize(2); k++) + { + for (int i = region.GetIndex(0); i < region.GetIndex(0) + (int)region.GetSize(0); i++) + { + int j = region.GetIndex(1); + u = matrix[0][0] * i + matrix[0][1] * j + matrix[0][2] * k + matrix[0][3]; + v = matrix[1][0] * i + matrix[1][2] * k + matrix[1][3]; + w = matrix[2][0] * i + matrix[2][2] * k + matrix[2][3]; + + // Apply perspective + w = 1 / w; + u = u * w - pIndex[0]; + v = v * w - pIndex[1]; + du = w * matrix[0][1]; + + + vi = itk::Math::floor(v); + if (vi >= 0 && vi < (int)pSize[1] - 1) + { + pVol = pVolZeroPointer + i + vBufferSize[0] * (j + k * vBufferSize[1]); + for (; j < (region.GetIndex(1) + (int)region.GetSize(1)); j++, pVol += vBufferSize[0], u += du) + { + ui = itk::Math::floor(u); + if (ui >= 0 && ui < (int)pSize[0] - 1) + { + indexOutput[0] = i; + indexOutput[1] = j; + indexOutput[2] = k; + using ComponentType = typename itk::PixelTraits::ValueType; + ComponentType u1, u2, v1, v2; + pProj = projection->GetBufferPointer() + vi * pSize[0] + ui; + v1 = v - vi; + v2 = 1.0 - v1; + u1 = u - ui; + u2 = 1.0 - u1; + + itk::Vector piline = itPIL.Get(); + typename TOutputImage::InternalPixelType sb = piline[0]; + typename TOutputImage::InternalPixelType st = piline[1]; + + + double rho = 0; + double d_in = (lambda - sb) / delta_lambda; + double d_out = (st - lambda) / delta_lambda; + + // Apply PiLineConditions + //if (lambda < sb - delta_lambda) + // rho = 0; + //else if (lambda < sb) + // rho = 0.5 * (d_in + 1) * (d_in + 1); + //else if (lambda < sb + delta_lambda) + // rho = 0.5 + d_in + 0.5 * d_in * d_in; + //else if (lambda < st - delta_lambda) + // rho = 1; + //else if (lambda < st) + // rho = 0.5 + d_out + 0.5 * d_out * d_out; + //else if (lambda < st + delta_lambda) + // rho = 0.5 * (1 + d_out * (1 + d_out)); + //else + // rho = 0; + + if (lambda < sb) + rho = 0.; + else if (lambda > sb && lambda < st) + rho = 1.; + else + rho = 0.; + + this->GetOutput()->TransformIndexToPhysicalPoint(indexOutput, outPoint); + double wstar = sid - outPoint[0] * sin(lambda) - outPoint[2] * cos(lambda); + + + *pVol += rho * delta_lambda * + (v2 * (u2 * *(pProj) + u1 * *(pProj + 1)) + + v1 * (u2 * *(pProj + pSize[0]) + u1 * *(pProj + pSize[0] + 1))) / + (wstar * 2 * M_PI); + } + ++itPIL; + } // j + } + } // i + } // k +} + +template +template +typename TProjectionImage::Pointer +KatsevichBackProjectionImageFilter::GetProjection(const unsigned int iProj) +{ + + typename Superclass::InputImagePointer stack = const_cast(this->GetInput(1)); + + const int iProjBuff = stack->GetBufferedRegion().GetIndex(ProjectionImageType::ImageDimension); + + typename TProjectionImage::Pointer projection = TProjectionImage::New(); + typename TProjectionImage::RegionType region; + typename TProjectionImage::SpacingType spacing; + typename TProjectionImage::PointType origin; + + for (unsigned int i = 0; i < TProjectionImage::ImageDimension; i++) + { + origin[i] = stack->GetOrigin()[i]; + spacing[i] = stack->GetSpacing()[i]; + region.SetSize(i, stack->GetBufferedRegion().GetSize()[i]); + region.SetIndex(i, stack->GetBufferedRegion().GetIndex()[i]); + } + if (this->GetTranspose()) + { + typename TProjectionImage::SizeType size = region.GetSize(); + typename TProjectionImage::IndexType index = region.GetIndex(); + std::swap(size[0], size[1]); + std::swap(index[0], index[1]); + std::swap(origin[0], origin[1]); + std::swap(spacing[0], spacing[1]); + region.SetSize(size); + region.SetIndex(index); + } + projection->SetSpacing(spacing); + projection->SetOrigin(origin); + projection->SetRegions(region); + projection->Allocate(); + + const unsigned int npixels = projection->GetBufferedRegion().GetNumberOfPixels(); + const InternalInputPixelType * pi = stack->GetBufferPointer() + (iProj - iProjBuff) * npixels; + InternalInputPixelType * po = projection->GetBufferPointer(); + + // Transpose projection for optimization + if (this->GetTranspose()) + { + for (unsigned int j = 0; j < region.GetSize(0); j++, po -= npixels - 1) + for (unsigned int i = 0; i < region.GetSize(1); i++, po += region.GetSize(0)) + *po = *pi++; + } + else + for (unsigned int i = 0; i < npixels; i++) + *po++ = *pi++; + + return projection; +} + +template +typename KatsevichBackProjectionImageFilter::ProjectionMatrixType +KatsevichBackProjectionImageFilter::GetIndexToIndexProjectionMatrix(const unsigned int iProj) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + + ProjectionMatrixType VolumeIndexToProjectionPhysicalPointMatrix = + GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); + + itk::Matrix ProjectionPhysicalPointToProjectionIndexMatrix = + GetProjectionPhysicalPointToProjectionIndexMatrix(iProj); + + return ProjectionMatrixType(ProjectionPhysicalPointToProjectionIndexMatrix.GetVnlMatrix() * + VolumeIndexToProjectionPhysicalPointMatrix.GetVnlMatrix()); +} + +template +typename KatsevichBackProjectionImageFilter::ProjectionMatrixType +KatsevichBackProjectionImageFilter::GetVolumeIndexToProjectionPhysicalPointMatrix( + const unsigned int iProj) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + + itk::Matrix matrixVol = + GetIndexToPhysicalPointMatrix(this->GetOutput()); + + return ProjectionMatrixType(this->m_Geometry->GetMagnificationMatrices()[iProj].GetVnlMatrix() * + this->m_Geometry->GetSourceTranslationMatrices()[iProj].GetVnlMatrix() * + this->m_Geometry->GetRotationMatrices()[iProj].GetVnlMatrix() * matrixVol.GetVnlMatrix()); +} + +template +itk::Matrix +KatsevichBackProjectionImageFilter::GetProjectionPhysicalPointToProjectionIndexMatrix( + const unsigned int iProj) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + + itk::Matrix matrixStackProj = + GetPhysicalPointToIndexMatrix(this->GetInput(1)); + + itk::Matrix matrixProj; + matrixProj.SetIdentity(); + for (unsigned int i = 0; i < Dimension - 1; i++) + { + matrixProj[i][Dimension - 1] = matrixStackProj[i][Dimension]; + for (unsigned int j = 0; j < Dimension - 1; j++) + matrixProj[i][j] = matrixStackProj[i][j]; + } + + // Transpose projection for optimization + itk::Matrix matrixFlip; + matrixFlip.SetIdentity(); + if (this->GetTranspose()) + { + std::swap(matrixFlip[0][0], matrixFlip[0][1]); + std::swap(matrixFlip[1][0], matrixFlip[1][1]); + } + + return itk::Matrix( + matrixFlip.GetVnlMatrix() * matrixProj.GetVnlMatrix() * + this->m_Geometry->GetProjectionTranslationMatrices()[iProj].GetVnlMatrix()); +} + +} // end namespace rtk + +#endif diff --git a/include/rtkKatsevichBackwardBinningImageFilter.h b/include/rtkKatsevichBackwardBinningImageFilter.h new file mode 100644 index 000000000..2fc45904c --- /dev/null +++ b/include/rtkKatsevichBackwardBinningImageFilter.h @@ -0,0 +1,142 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichBackwardBinningImageFilter_h +#define rtkKatsevichBackwardBinningImageFilter_h + +#include + +#include +#include + +#include +#include +#include "itkImageSliceIteratorWithIndex.h" + +#include +#include + +namespace rtk +{ + +/** \class KatsevichBackwardBinningImageFilter + * \brief + * + * Backprojects a stack of projection images (input 1) in a 3D volume (input 0) + * using linear interpolation according to a specified geometry. The operation + * is voxel-based, meaning that the center of each voxel is projected in the + * projection images to determine the interpolation location. + * + * \test rtkfovtest.cxx + * + * \author Simon Rit + * + * \ingroup RTK Projector + */ +template +class KatsevichBackwardBinningImageFilter : public itk::ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichBackwardBinningImageFilter); + + /** Standard class type alias. */ + using Self = KatsevichBackwardBinningImageFilter; + using OutputImageType = TOutputImage; + using InputImageType = TInputImage; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + using InternalInputPixelType = typename TInputImage::InternalPixelType; + using OutputImageRegionType = typename TOutputImage::RegionType; + using OutputImageIndexType = typename TOutputImage::IndexType; + using OutputImageSizeType = typename TOutputImage::SizeType; + using InputImageRegionType = typename TInputImage::RegionType; + using InputImageIndexType = typename TInputImage::IndexType; + using InputImageSizeType = typename TInputImage::SizeType; + using InputSliceIteratorType = itk::ImageSliceIteratorWithIndex; + using OutputSliceIteratorType = itk::ImageSliceIteratorWithIndex; + using ConstantImageType = rtk::ConstantImageSource; + + + using GeometryType = rtk::ThreeDHelicalProjectionGeometry; + using GeometryPointer = GeometryType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(KatsevichBackwardBinningImageFilter, itk::ImageToImageFilter); + + /** Get / Set the object pointer to projection geometry */ + itkGetMacro(Geometry, GeometryPointer); + itkSetMacro(Geometry, GeometryPointer); + + /** Get / Set the transpose flag for 2D projections (optimization trick) */ + itkGetMacro(Transpose, bool); + itkSetMacro(Transpose, bool); + +protected: + KatsevichBackwardBinningImageFilter() = default; + ~KatsevichBackwardBinningImageFilter() override = default; + + /** Set up requested region for inputs*/ + void + GenerateInputRequestedRegion() override; + + void + GenerateOutputInformation() override; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + + /** The two inputs should not be in the same space so there is nothing + * to verify. */ + void + VerifyInputInformation() const override + {} + + /** The input is a stack of projections, we need to interpolate in one projection + for efficiency during interpolation. Use of itk::ExtractImageFilter is + not threadsafe in ThreadedGenerateData, this one is. The output can be multiplied by a constant. + The function is templated to allow getting an itk::CudaImage. */ + template + typename TProjectionImage::Pointer + GetProjection(const unsigned int iProj); + + /** RTK geometry object */ + GeometryPointer m_Geometry; + +private: + /** Flip projection flag: infludences GetProjection and + GetIndexToIndexProjectionMatrix for optimization */ + bool m_Transpose{ false }; +}; + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkKatsevichBackwardBinningImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkKatsevichBackwardBinningImageFilter.hxx b/include/rtkKatsevichBackwardBinningImageFilter.hxx new file mode 100644 index 000000000..b0a79e4f5 --- /dev/null +++ b/include/rtkKatsevichBackwardBinningImageFilter.hxx @@ -0,0 +1,255 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichBackwardBinningImageFilter_hxx +#define rtkKatsevichBackwardBinningImageFilter_hxx + +#include "math.h" + +#include + +#include + +#include +#include +#include +#include + +namespace rtk +{ + +template +void +KatsevichBackwardBinningImageFilter::GenerateInputRequestedRegion() +{ + // // Input 0 is Output ie result of forward binning + // typename Superclass::InputImagePointer inputPtr0 = const_cast(this->GetInput()); + // std::cout << "Input largest : " << this->GetInput()->GetLargestPossibleRegion() << std::endl; + // std::cout << "inputPtr0 : " << inputPtr0->GetRequestedRegion() << std::endl; + // if (!inputPtr0) + // return; + // inputPtr0->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); + // + // //TO DELETE : SPACING AND ORIGIN ALWAYS (1,1,1) AND (0,0,0) -> NOT OKAY !!!! + // //std::cout<<"Input size: "<GetInput()->GetLargestPossibleRegion().GetSize()<GetInput()->GetSpacing()<GetInput()->GetOrigin()<GetLargest = inputPtr0->GetRequested before assignment. + +template +void +KatsevichBackwardBinningImageFilter::GenerateOutputInformation() +{ + // Update output size and origin. + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); + + // Get pointers to the input and output + const InputImageType * inputPtr = this->GetInput(); + OutputImageType * outputPtr = this->GetOutput(); + + itkAssertInDebugAndIgnoreInReleaseMacro(inputPtr); + itkAssertInDebugAndIgnoreInReleaseMacro(outputPtr != nullptr); + + unsigned int i; + const typename TInputImage::SpacingType & inputSpacing = inputPtr->GetSpacing(); + const typename TInputImage::PointType & inputOrigin = inputPtr->GetOrigin(); + const typename TInputImage::SizeType & inputSize = inputPtr->GetLargestPossibleRegion().GetSize(); + const typename TInputImage::IndexType & inputStartIndex = inputPtr->GetLargestPossibleRegion().GetIndex(); + typename TOutputImage::IndexType outputStartIndex; + + typename TOutputImage::SpacingType outputSpacing(inputSpacing); + typename TOutputImage::SizeType outputSize(inputSize); + typename TOutputImage::PointType outputOrigin(inputOrigin); + + outputSize[1] = (int)(inputSize[1] - 1) / 2; + outputSpacing[1] = outputSpacing[0]; + outputOrigin[1] = -0.5 * outputSpacing[1] * (outputSize[1] - 1); + + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetIndex(inputStartIndex); + outputLargestPossibleRegion.SetSize(outputSize); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + outputPtr->SetOrigin(outputOrigin); + outputPtr->SetSpacing(outputSpacing); +} + +template +void +KatsevichBackwardBinningImageFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull() || !this->m_Geometry->GetTheGeometryIsVerified()) + itkExceptionMacro(<< "Geometry has not been set or not been verified"); +} + + +template +void +KatsevichBackwardBinningImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + + // Iterator over output projections (in v) + OutputSliceIteratorType itOut(this->GetOutput(), outputRegionForThread); + itOut.SetFirstDirection(1); + itOut.SetSecondDirection(0); + + // Detector details + unsigned int nbproj = outputRegionForThread.GetSize()[2]; + unsigned int det_nu = outputRegionForThread.GetSize()[0]; + unsigned int det_nv = outputRegionForThread.GetSize()[1]; + + unsigned int idx_proj = outputRegionForThread.GetIndex()[2]; + + InputImageRegionType requiredInputRegion; + InputImageIndexType inIndex = this->GetInput()->GetLargestPossibleRegion().GetIndex(); + InputImageSizeType inSize = this->GetInput()->GetLargestPossibleRegion().GetSize(); + inIndex[2] = idx_proj; + inSize[2] = nbproj; + requiredInputRegion.SetIndex(inIndex); + requiredInputRegion.SetSize(inSize); + + InputSliceIteratorType itIn(this->GetInput(), requiredInputRegion); + itIn.SetFirstDirection(1); + itIn.SetSecondDirection(0); + + double u = 0; + double v_current = 0; + double v_psi = 0; + double psi = 0; + + const double u_spacing = this->GetInput()->GetSpacing()[0]; + const double psi_spacing = this->GetInput()->GetSpacing()[1]; + const double v_spacing = this->GetOutput()->GetSpacing()[1]; + const double u_orig = this->GetInput()->GetOrigin()[0]; + const double psi_orig = this->GetInput()->GetOrigin()[1]; + const double v_orig = this->GetOutput()->GetOrigin()[1]; + + // Geometry details + const double D = m_Geometry->GetHelixSourceToDetectorDist(); + const double P = m_Geometry->GetHelixPitch(); + const double R = m_Geometry->GetHelixRadius(); + + const double alpha_m = + atan(0.5 * u_spacing * (det_nu - 1) / D); // Générique mais doit être cohérent avec le fwd binning. + const double r = R * sin(alpha_m); + const int L = this->GetInput()->GetLargestPossibleRegion().GetSize(1); + + unsigned int idx_v = 0; + double coeff_s = 0., coeff_i = 0.; + itIn.GoToBegin(); + itOut.GoToBegin(); + while (!itIn.IsAtEnd()) + { + while (!itIn.IsAtEndOfSlice()) + { + // Vectors of output values and coeffs accumulator + std::vector coeffs(det_nv, 0.); + std::vector outputVals(det_nv, 0.); + while (!itIn.IsAtEndOfLine()) + { + const typename TInputImage::IndexType idx = itIn.GetIndex(); + psi = psi_orig + idx[1] * psi_spacing; + u = u_orig + idx[0] * u_spacing; + + // Compute v corresponding to current psi value + if (psi != 0) + v_psi = (psi + (psi / tan(psi)) * (u / D)) * D * P / (R * 2 * M_PI); + else + v_psi = (u / D) * D * P / (R * 2 * M_PI); + + // Raise error if v outside detector range + if (v_psi < v_orig - 0.5 * v_spacing || v_psi >= v_orig + v_spacing * (det_nv - 1 + 0.5)) + { + itkExceptionMacro(<< "The v_k(psi) value is outside v-range. v_psi :" << v_psi << " v_orig, v_spacing " + << v_orig << " " << v_spacing << "psi : " << psi << " det_nv : " << det_nv); + } + else if (v_psi < v_orig) // v is before first pixel center but within the detector range + { + idx_v = itk::Math::floor((v_psi - v_orig) / v_spacing); + v_current = v_orig + idx_v * v_spacing; + coeff_s = (v_psi - v_current) / v_spacing; + // Accumulate + outputVals[0] += coeff_s * itIn.Get(); + coeffs[0] += coeff_s; + } + else if (v_psi >= + v_orig + v_spacing * (det_nv - 1)) // v is after last pixel center but within the detector range + { + idx_v = itk::Math::floor((v_psi - v_orig) / v_spacing); + v_current = v_orig + idx_v * v_spacing; + coeff_s = (v_psi - v_current) / v_spacing; + coeff_i = 1 - coeff_s; + outputVals[det_nv - 1] += coeff_i * itIn.Get(); + coeffs[det_nv - 1] += coeff_i; + } + else // standard case + { + idx_v = itk::Math::floor((v_psi - v_orig) / v_spacing); + v_current = v_orig + idx_v * v_spacing; + coeff_s = (v_psi - v_current) / v_spacing; + coeff_i = 1 - coeff_s; + + // Accumulate + outputVals[idx_v] += coeff_i * itIn.Get(); + outputVals[idx_v + 1] += coeff_s * itIn.Get(); + coeffs[idx_v] += coeff_i; + coeffs[idx_v + 1] += coeff_s; + } + // Move iterator forward + ++itIn; + } + + // Divide output by the sum of all coeffs contributions. + while (!itOut.IsAtEndOfLine()) + { + OutputImageIndexType index = itOut.GetIndex(); + int idx = index[1]; + if (coeffs[idx] == 0.) + { + itOut.Set(0.); + } + else + { + itOut.Set(outputVals[idx] / coeffs[idx]); + } + if (this->m_Geometry->GetRadiusCylindricalDetector() != 0.) + { + double cosalpha = cos((u_orig + index[0] * u_spacing) / this->m_Geometry->GetRadiusCylindricalDetector()); + itOut.Set(itOut.Get() * cosalpha); + } + + ++itOut; + } + itIn.NextLine(); + itOut.NextLine(); + } + itIn.NextSlice(); + itOut.NextSlice(); + } +} + +} // end namespace rtk + +#endif diff --git a/include/rtkKatsevichConeBeamReconstructionFilter.h b/include/rtkKatsevichConeBeamReconstructionFilter.h new file mode 100644 index 000000000..2f676c335 --- /dev/null +++ b/include/rtkKatsevichConeBeamReconstructionFilter.h @@ -0,0 +1,215 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichConeBeamReconstructionFilter_h +#define rtkKatsevichConeBeamReconstructionFilter_h + +#include +#include "rtkKatsevichDerivativeImageFilter.h" +#include "rtkKatsevichForwardBinningImageFilter.h" +#include "rtkFFTHilbertImageFilter.h" +#include "rtkKatsevichBackwardBinningImageFilter.h" +#include "rtkKatsevichBackProjectionImageFilter.h" + +#include "rtkConfiguration.h" + + +namespace rtk +{ + +/** \class KatsevichConeBeamReconstructionFilter + * \brief Implements Feldkamp, David and Kress cone-beam reconstruction + * + * KatsevichConeBeamReconstructionFilter is a mini-pipeline filter which combines + * the different steps of the Katsevich cone-beam reconstruction filter: + * - rtk::KatsevichDerivativeImageFilter for derivative and cosine-weighting of the projections, + * - rtk::KatsevichForwardBinningImageFilter, rtk::FFTHilbertImageFilter an d rtk::KatsevichBackwardBinningImageFilter + * for Hilbert filtering, + * - rtk::KatsevichBackProjectionImageFilter for backprojection. + * The input stack of projections is processed piece by piece (the size is + * controlled with ProjectionSubsetSize) via the use of itk::ExtractImageFilter + * to extract sub-stacks. + * Note though that this behaviour only mimicks the FDK filter. The ProjectionSubsetSize is set + * to the total number of projections (because the derivative filter does not process projections + * independently) + * + * \dot + * digraph KatsevichConeBeamReconstructionFilter { + * node [shape=box]; + * 1 [ label="rtk::KatsevichDerivativeImageFilter" URL="\ref rtk::KatsevichDerivativeImageFilter"]; + * 2 [ label="rtk::KatsevichForwardBinningImageFilter" URL="\ref rtk::KatsevichForwardBinningImageFilter"]; + * 3 [ label="rtk::FFTHilbertImageFilter" URL="\ref rtk::FFTHilbertImageFilter"]; + * 4 [ label="rtk::KatsevichBackwardBinningImageFilter" URL="\ref rtk::KatsevichBackwardBinningImageFilter"]; + * 5 [ label="rtk::KatsevichBackProjectionImageFilter" URL="\ref rtk::KatsevichBackProjectionImageFilter"]; + * 1 -> 2; + * 2 -> 3; + * 3 -> 4; + * 4 -> 5; + * } + * \enddot + * + * \test + * + * \author Jerome Lesaint + * + * \ingroup RTK ReconstructionAlgorithm + */ +template +class ITK_TEMPLATE_EXPORT KatsevichConeBeamReconstructionFilter + : public itk::InPlaceImageFilter +{ +public: +#if ITK_VERSION_MAJOR == 5 && ITK_VERSION_MINOR == 1 + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichConeBeamReconstructionFilter); +#else + ITK_DISALLOW_COPY_AND_MOVE(KatsevichConeBeamReconstructionFilter); +#endif + + /** Standard class type alias. */ + using Self = KatsevichConeBeamReconstructionFilter; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Some convenient type alias. */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + + /** Typedefs of each subfilter of this composite filter */ + using ExtractFilterType = itk::ExtractImageFilter; + using DerivativeFilterType = rtk::KatsevichDerivativeImageFilter; + using ForwardFilterType = rtk::KatsevichForwardBinningImageFilter; + using HilbertFilterType = rtk::FFTHilbertImageFilter; + using BackwardFilterType = rtk::KatsevichBackwardBinningImageFilter; + using BackProjectionFilterType = rtk::KatsevichBackProjectionImageFilter; + using BackProjectionFilterPointer = typename BackProjectionFilterType::Pointer; + + /** Standard New method. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(KatsevichConeBeamReconstructionFilter, itk::ImageToImageFilter); + + /** Get / Set the object pointer to projection geometry */ + itkGetModifiableObjectMacro(Geometry, ThreeDHelicalProjectionGeometry); + itkSetObjectMacro(Geometry, ThreeDHelicalProjectionGeometry); + + /** Get pointer to the extract filter used by the Katsevich reconstruction */ + typename DerivativeFilterType::Pointer + GetExtractFilter() + { + return m_ExtractFilter; + } + + /** Get pointer to the derivative filter used by the Katsevich reconstruction */ + typename DerivativeFilterType::Pointer + GetDerivativeFilter() + { + return m_DerivativeFilter; + } + + /** Get pointer to the forward binning filter used by the Katsevich reconstruction */ + typename ForwardFilterType::Pointer + GetForwardFilter() + { + return m_ForwardFilter; + } + + /** Get pointer to the Hilbert filter used by the Katsevich reconstruction */ + typename HilbertFilterType::Pointer + GetHilbertFilter() + { + return m_HilbertFilter; + } + + /** Get pointer to the forward binning filter used by the Katsevich reconstruction */ + typename BackwardFilterType::Pointer + GetBackwardFilter() + { + return m_BackwardFilter; + } + + /** Get pointer to the backprojection filter used by the Katsevich reconstruction */ + typename BackProjectionFilterType::Pointer + GetBackProjectionFilter() + { + return m_BackProjectionFilter; + } + + //////////////////////////// + // JL : the default here is set to the total number of projections. + // In another words, the extract filter does NOTHING ! + //////////////////////////// + /** Get / Set the number of cone-beam projection images processed + simultaneously. Default is 4. */ + itkGetMacro(ProjectionSubsetSize, unsigned int); + itkSetMacro(ProjectionSubsetSize, unsigned int); + + ///** Get / Set and init the backprojection filter. The set function takes care + // * of initializing the mini-pipeline and the ramp filter must therefore be + // * created before calling this set function. */ + // itkGetMacro(BackProjectionFilter, BackProjectionFilterPointer); + // virtual void + // SetBackProjectionFilter(const BackProjectionFilterPointer _arg); + +protected: + KatsevichConeBeamReconstructionFilter(); + ~KatsevichConeBeamReconstructionFilter() override = default; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + void + GenerateInputRequestedRegion() override; + + void + GenerateOutputInformation() override; + + void + GenerateData() override; + + /** The two inputs should not be in the same space so there is nothing + * to verify. */ + void + VerifyInputInformation() const override + {} + + /** Pointers to each subfilter of this composite filter */ + typename ExtractFilterType::Pointer m_ExtractFilter; + typename DerivativeFilterType::Pointer m_DerivativeFilter; + typename ForwardFilterType::Pointer m_ForwardFilter; + typename HilbertFilterType::Pointer m_HilbertFilter; + typename BackwardFilterType::Pointer m_BackwardFilter; + BackProjectionFilterPointer m_BackProjectionFilter; + +private: + /** Number of projections processed at a time. */ + unsigned int m_ProjectionSubsetSize{ 16 }; + + /** Geometry propagated to subfilters of the mini-pipeline. */ + ThreeDHelicalProjectionGeometry::Pointer m_Geometry; +}; // end of class + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkKatsevichConeBeamReconstructionFilter.hxx" +#endif + +#endif diff --git a/include/rtkKatsevichConeBeamReconstructionFilter.hxx b/include/rtkKatsevichConeBeamReconstructionFilter.hxx new file mode 100644 index 000000000..ea3bfc922 --- /dev/null +++ b/include/rtkKatsevichConeBeamReconstructionFilter.hxx @@ -0,0 +1,192 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichConeBeamReconstructionFilter_hxx +#define rtkKatsevichConeBeamReconstructionFilter_hxx + +#include "rtkKatsevichConeBeamReconstructionFilter.h" + +#include + +namespace rtk +{ + +template +KatsevichConeBeamReconstructionFilter::KatsevichConeBeamReconstructionFilter() +{ + this->SetNumberOfRequiredInputs(2); + + // Create each filter of the composite filter + m_ExtractFilter = ExtractFilterType::New(); + m_DerivativeFilter = DerivativeFilterType::New(); + m_ForwardFilter = ForwardFilterType::New(); + m_HilbertFilter = HilbertFilterType::New(); + m_BackwardFilter = BackwardFilterType::New(); + m_BackProjectionFilter = BackProjectionFilterType::New(); + + // Permanent internal connections + m_DerivativeFilter->SetInput(m_ExtractFilter->GetOutput()); + m_ForwardFilter->SetInput(m_DerivativeFilter->GetOutput()); + m_HilbertFilter->SetInput(m_ForwardFilter->GetOutput()); + m_HilbertFilter->SetPixelShift(0.5); // Hard-coded. Make it a parameter ? + m_BackwardFilter->SetInput(m_HilbertFilter->GetOutput()); + + // Bypass the SetBackProjectionFunction (since only the Katsevich BP is allowed) + // Need to set Input 1. + // this->SetBackProjectionFilter(BackProjectionFilterType::New()); + m_BackProjectionFilter->SetInput(1, m_BackwardFilter->GetOutput()); + + // Default parameters + m_ExtractFilter->SetDirectionCollapseToSubmatrix(); + // m_DerivativeFilter->InPlaceOn(); + + // Default to one projection per subset when FFTW is not available +#if !defined(USE_FFTWD) + if (typeid(TFFTPrecision).name() == typeid(double).name()) + m_ProjectionSubsetSize = 2; +#endif +#if !defined(USE_FFTWF) + if (typeid(TFFTPrecision).name() == typeid(float).name()) + m_ProjectionSubsetSize = 2; +#endif +} + +template +void +KatsevichConeBeamReconstructionFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull()) + itkExceptionMacro(<< "Geometry has not been set."); +} + +template +void +KatsevichConeBeamReconstructionFilter::GenerateInputRequestedRegion() +{ + typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); + if (!inputPtr) + return; + + // SR: is this useful? + m_ExtractFilter->SetInput(this->GetInput(1)); + m_BackProjectionFilter->SetInput(0, this->GetInput(0)); + m_BackProjectionFilter->SetInPlace(this->GetInPlace()); + m_BackProjectionFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); + m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); +} + +template +void +KatsevichConeBeamReconstructionFilter::GenerateOutputInformation() +{ + const unsigned int Dimension = this->InputImageDimension; + + m_DerivativeFilter->SetGeometry(m_Geometry); + m_ForwardFilter->SetGeometry(m_Geometry); + m_BackwardFilter->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; + projRegion = this->GetInput(1)->GetLargestPossibleRegion(); + unsigned int firstStackSize = (unsigned int)projRegion.GetSize(Dimension - 1); + projRegion.SetSize(Dimension - 1, firstStackSize); + m_ExtractFilter->SetExtractionRegion(projRegion); + + // Run composite filter update + m_BackProjectionFilter->SetInput(0, this->GetInput(0)); + m_BackProjectionFilter->SetInPlace(this->GetInPlace()); + m_ExtractFilter->SetInput(this->GetInput(1)); + m_BackProjectionFilter->UpdateOutputInformation(); + + // Update output information + this->GetOutput()->SetOrigin(m_BackProjectionFilter->GetOutput()->GetOrigin()); + this->GetOutput()->SetSpacing(m_BackProjectionFilter->GetOutput()->GetSpacing()); + this->GetOutput()->SetDirection(m_BackProjectionFilter->GetOutput()->GetDirection()); + this->GetOutput()->SetLargestPossibleRegion(m_BackProjectionFilter->GetOutput()->GetLargestPossibleRegion()); +} + +template +void +KatsevichConeBeamReconstructionFilter::GenerateData() +{ + const unsigned int Dimension = this->InputImageDimension; + + // The backprojection can work on a smaller stack of projections, not the full stack + // Here, it takes the whole stack. + typename ExtractFilterType::InputImageRegionType subsetRegion; + subsetRegion = this->GetInput(1)->GetLargestPossibleRegion(); + unsigned int nProj = subsetRegion.GetSize(Dimension - 1); + m_ProjectionSubsetSize = nProj; + + // The progress accumulator tracks the progress of the pipeline + // Each filter is equally weighted across all iterations of the stack + itk::ProgressAccumulator::Pointer progress = itk::ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + auto frac = (1.0f / 5) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize); + progress->RegisterInternalFilter(m_DerivativeFilter, frac); + progress->RegisterInternalFilter(m_ForwardFilter, frac); + progress->RegisterInternalFilter(m_HilbertFilter, frac); + progress->RegisterInternalFilter(m_BackwardFilter, frac); + progress->RegisterInternalFilter(m_BackProjectionFilter, frac); + + for (unsigned int i = 0; i < nProj; i += m_ProjectionSubsetSize) + { + // After the first bp update, we need to use its output as input. + if (i) + { + typename TInputImage::Pointer pimg = m_BackProjectionFilter->GetOutput(); + pimg->DisconnectPipeline(); + m_BackProjectionFilter->SetInput(pimg); + + // Change projection subset + subsetRegion.SetIndex(Dimension - 1, i); + subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i)); + m_ExtractFilter->SetExtractionRegion(subsetRegion); + + // This is required to reset the full pipeline + m_BackProjectionFilter->GetOutput()->UpdateOutputInformation(); + m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion(); + } + m_BackProjectionFilter->Update(); + } + + this->GraftOutput(m_BackProjectionFilter->GetOutput()); + this->GenerateOutputInformation(); +} + +// template +// void +// KatsevichConeBeamReconstructionFilter::SetBackProjectionFilter( +// const BackProjectionFilterPointer _arg) +//{ +// itkDebugMacro("setting BackProjectionFilter to " << _arg); +// if (this->m_BackProjectionFilter != _arg) +// { +// this->m_BackProjectionFilter = _arg; +// m_BackProjectionFilter->SetInput(1, m_RampFilter->GetOutput()); +// this->Modified(); +// } +//} + +} // end namespace rtk + +#endif // rtkKatsevichConeBeamReconstructionFilter_hxx diff --git a/include/rtkKatsevichDerivativeImageFilter.h b/include/rtkKatsevichDerivativeImageFilter.h new file mode 100644 index 000000000..113d6c3c2 --- /dev/null +++ b/include/rtkKatsevichDerivativeImageFilter.h @@ -0,0 +1,128 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichDerivativeImageFilter_h +#define rtkKatsevichDerivativeImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" +#include "itkIndent.h" + +#include "itkConstShapedNeighborhoodIterator.h" +#include "itkImageRegionIterator.h" + +#include "itkImageSliceIteratorWithIndex.h" + +#include "rtkThreeDHelicalProjectionGeometry.h" + + +namespace rtk +{ + +/** \class KatsevichDerivativeImageFilter + * \brief Computes the derivatives of the projections wrt the projection index (the rotation angle). + * + * This filter implements the derivative as described in Noo et al., PMB, 2003 + * + * \author Jerome Lesaint + * + * \ingroup + */ +template +class ITK_EXPORT KatsevichDerivativeImageFilter : public itk::ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichDerivativeImageFilter); + + /** Standard class type alias. */ + using Self = KatsevichDerivativeImageFilter; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(KatsevichDerivativeImageFilter, ImageToImageFilter); + + /** Typedef to images */ + using OutputImageType = TOutputImage; + using InputImageType = TInputImage; + using OutputImagePointer = typename OutputImageType::Pointer; + using InputImagePointer = typename InputImageType::Pointer; + using InputImageConstPointer = typename InputImageType::ConstPointer; + using OutputImageRegionType = typename TOutputImage::RegionType; + + using ConstShapedNeighborhoodIteratorType = itk::ConstShapedNeighborhoodIterator; + using IteratorType = itk::ImageRegionIterator; + using SliceIteratorType = itk::ImageSliceIteratorWithIndex; + using IndexType = typename InputImageType::IndexType; + + + using OutputPixelType = typename TOutputImage::PixelType; + using OutputInternalPixelType = typename TOutputImage::InternalPixelType; + using InputPixelType = typename TInputImage::PixelType; + using InputInternalPixelType = typename TInputImage::InternalPixelType; + static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension; + + /** Get/ Set geometry structure */ + itkGetMacro(Geometry, ::rtk::ThreeDHelicalProjectionGeometry::Pointer); + itkSetObjectMacro(Geometry, ::rtk::ThreeDHelicalProjectionGeometry); + + + void + GenerateOutputInformation() override; + + void + GenerateInputRequestedRegion() override; + +protected: + KatsevichDerivativeImageFilter(); + ~KatsevichDerivativeImageFilter() override = default; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + void + PrintSelf(std::ostream & os, itk::Indent indent) const override; + + /** Standard pipeline method. While this class does not implement a + * ThreadedGenerateData(), its GenerateData() delegates all + * calculations to an NeighborhoodOperatorImageFilter. Since the + * NeighborhoodOperatorImageFilter is multithreaded, this filter is + * multithreaded by default. */ + void + GenerateData() override; + + // void + // DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + +private: + ThreeDHelicalProjectionGeometry::Pointer m_Geometry; +}; + + +} // end namespace rtk + +#ifndef rtk_MANUAL_INSTANTIATION +# include "rtkKatsevichDerivativeImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkKatsevichDerivativeImageFilter.hxx b/include/rtkKatsevichDerivativeImageFilter.hxx new file mode 100644 index 000000000..e8d91d008 --- /dev/null +++ b/include/rtkKatsevichDerivativeImageFilter.hxx @@ -0,0 +1,253 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichDerivativeImageFilter_hxx +#define rtkKatsevichDerivativeImageFilter_hxx + +#include "rtkKatsevichDerivativeImageFilter.h" + +#include "itkNeighborhoodOperatorImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkProgressAccumulator.h" +#include "itkObjectFactory.h" + + +namespace rtk +{ + +/** + * Constructor + */ +template +KatsevichDerivativeImageFilter::KatsevichDerivativeImageFilter() +{ + this->SetNumberOfRequiredInputs(1); +} + +template +void +KatsevichDerivativeImageFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull() || !this->m_Geometry->GetTheGeometryIsVerified()) + itkExceptionMacro(<< "Geometry has not been set or not been checked"); +} + +template +void +KatsevichDerivativeImageFilter::GenerateOutputInformation() +{ + // Update output size and origin + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); + + // Get pointers to the input and output + const InputImageType * inputPtr = this->GetInput(); + OutputImageType * outputPtr = this->GetOutput(); + + itkAssertInDebugAndIgnoreInReleaseMacro(inputPtr); + itkAssertInDebugAndIgnoreInReleaseMacro(outputPtr != nullptr); + + const typename TInputImage::SpacingType & inputSpacing = inputPtr->GetSpacing(); + const typename TInputImage::PointType & inputOrigin = inputPtr->GetOrigin(); + const typename TInputImage::SizeType & inputSize = inputPtr->GetLargestPossibleRegion().GetSize(); + const typename TInputImage::IndexType & inputStartIndex = inputPtr->GetLargestPossibleRegion().GetIndex(); + typename TOutputImage::IndexType outputStartIndex; + + typename TOutputImage::SpacingType outputSpacing(inputSpacing); + typename TOutputImage::SizeType outputSize(inputSize); + typename TOutputImage::PointType outputOrigin(inputOrigin); + + for (int i = 0; i < 3; i++) + { + outputOrigin[i] += 0.5 * inputSpacing[i]; + // If cylindrical detector, sampling in v is unchanged. In case of a flat panel + // derivative is also computed at interlaced positions in v. + if (!(this->m_Geometry->GetRadiusCylindricalDetector() != 0. && i == 1)) + { + outputSize[i] -= 1; + } + if (outputSize[i] < 1) + { + itkExceptionMacro(<< "Output size is 0 on dimension " << i); + } + } + + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetIndex(inputStartIndex); + outputLargestPossibleRegion.SetSize(outputSize); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + outputPtr->SetOrigin(outputOrigin); + outputPtr->SetSpacing(outputSpacing); +} + +template +void +KatsevichDerivativeImageFilter::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); +} + +template +void +KatsevichDerivativeImageFilter::GenerateData() +{ + + // Get pointers to the input and output + const InputImageType * inputPtr = this->GetInput(); + OutputImageType * outputPtr = this->GetOutput(); + + if (this->m_Geometry->GetRadiusCylindricalDetector() != 0) // Curved detector + { + // ShapedIterator on InputImage + typename ConstShapedNeighborhoodIteratorType::OffsetType offset1 = { { 1, 0, 1 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset2 = { { 1, 0, 0 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset3 = { { 0, 0, 1 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset4 = { { 0, 0, 0 } }; + + typename ConstShapedNeighborhoodIteratorType::RadiusType radius; + radius.Fill(1); + + ConstShapedNeighborhoodIteratorType itIn(radius, inputPtr, inputPtr->GetRequestedRegion()); + itIn.ActivateOffset(offset1); + itIn.ActivateOffset(offset2); + itIn.ActivateOffset(offset3); + itIn.ActivateOffset(offset4); + + outputPtr->SetRegions(outputPtr->GetLargestPossibleRegion()); + outputPtr->Allocate(); + + IteratorType itOut(outputPtr, outputPtr->GetRequestedRegion()); + + typename InputImageType::SpacingType spacing = inputPtr->GetSpacing(); + typename InputImageType::PointType origin = inputPtr->GetOrigin(); + typename InputImageType::PixelType angular_gap = this->m_Geometry->GetHelixAngularGap(); + float D = this->m_Geometry->GetHelixSourceToDetectorDist(); + + for (itIn.GoToBegin(), itOut.GoToBegin(); !itIn.IsAtEnd(); ++itIn, ++itOut) + { + OutputPixelType t1 = 0.5 * (1 / spacing[0] + 1 / angular_gap) * itIn.GetPixel(offset1); + OutputPixelType t2 = 0.5 * (1 / spacing[0] - 1 / angular_gap) * itIn.GetPixel(offset2); + OutputPixelType t3 = 0.5 * (-1 / spacing[0] + 1 / angular_gap) * itIn.GetPixel(offset3); + OutputPixelType t4 = 0.5 * (-1 / spacing[0] - 1 / angular_gap) * itIn.GetPixel(offset4); + + // Length-correction weighting + typename InputImageType::IndexType index = itIn.GetIndex(); + OutputPixelType v = origin[1] + spacing[1] * index[1]; + + OutputPixelType length_correction = D / (sqrt(pow(D, 2) + pow(v, 2))); + itOut.Set(length_correction * (t1 + t2 + t3 + t4)); + } + } + else // Flat panel detector + { + + // ShapedIterator on InputImage + typename ConstShapedNeighborhoodIteratorType::OffsetType offset000 = { { 0, 0, 0 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset100 = { { 1, 0, 0 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset010 = { { 0, 1, 0 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset001 = { { 0, 0, 1 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset110 = { { 1, 1, 0 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset101 = { { 1, 0, 1 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset011 = { { 0, 1, 1 } }; + typename ConstShapedNeighborhoodIteratorType::OffsetType offset111 = { { 1, 1, 1 } }; + + typename ConstShapedNeighborhoodIteratorType::RadiusType radius; + radius.Fill(1); + + ConstShapedNeighborhoodIteratorType itIn(radius, inputPtr, inputPtr->GetRequestedRegion()); + itIn.ActivateOffset(offset000); + itIn.ActivateOffset(offset100); + itIn.ActivateOffset(offset010); + itIn.ActivateOffset(offset001); + itIn.ActivateOffset(offset110); + itIn.ActivateOffset(offset101); + itIn.ActivateOffset(offset011); + itIn.ActivateOffset(offset111); + + outputPtr->SetRegions(outputPtr->GetLargestPossibleRegion()); + outputPtr->Allocate(); + + IteratorType itOut(outputPtr, outputPtr->GetRequestedRegion()); + + typename InputImageType::SpacingType spacing = inputPtr->GetSpacing(); + typename InputImageType::PointType origin = inputPtr->GetOrigin(); + typename InputImageType::PixelType angular_gap = this->m_Geometry->GetHelixAngularGap(); + float D = this->m_Geometry->GetHelixSourceToDetectorDist(); + + for (itIn.GoToBegin(), itOut.GoToBegin(); !itIn.IsAtEnd(); ++itIn, ++itOut) + { + OutputPixelType p000 = itIn.GetPixel(offset000); + OutputPixelType p100 = itIn.GetPixel(offset100); + OutputPixelType p010 = itIn.GetPixel(offset010); + OutputPixelType p001 = itIn.GetPixel(offset001); + OutputPixelType p110 = itIn.GetPixel(offset110); + OutputPixelType p101 = itIn.GetPixel(offset101); + OutputPixelType p011 = itIn.GetPixel(offset011); + OutputPixelType p111 = itIn.GetPixel(offset111); + + typename InputImageType::IndexType index = itIn.GetIndex(); + OutputPixelType u = origin[0] + spacing[0] * (index[0] + 0.5); // Add half-pixel for interlaced positions. + OutputPixelType v = origin[1] + spacing[1] * (index[1] + 0.5); // + + + OutputPixelType dl = 0.25 * (1 / angular_gap) * (p001 - p000 + p011 - p010 + p101 - p100 + p111 - p110); + OutputPixelType dutmp = 0.25 * (1 / spacing[0]) * (p100 - p000 + p110 - p010 + p101 - p001 + p111 - p011); + OutputPixelType dvtmp = 0.25 * (1 / spacing[1]) * (p010 - p000 + p110 - p100 + p011 - p001 + p111 - p101); + + OutputPixelType length_correction = D / (sqrt(pow(D, 2) + pow(u, 2) + pow(v, 2))); + itOut.Set(length_correction * (dl + (pow(u, 2) + pow(D, 2)) / D * dutmp + u * v / D * dvtmp)); + } + } + + this->GraftOutput(outputPtr); + + // Update geometry to shift gantry angle by half angular step + rtk::ThreeDHelicalProjectionGeometry::Pointer new_geometry = rtk::ThreeDHelicalProjectionGeometry::New(); + const std::vector angles = this->GetGeometry()->GetHelicalAngles(); + const double angular_gap = this->GetGeometry()->GetHelixAngularGap(); + const double sid = this->GetGeometry()->GetHelixRadius(); + const double sdd = this->GetGeometry()->GetHelixSourceToDetectorDist(); + const std::vector dy = this->GetGeometry()->GetSourceOffsetsY(); + const double vertical_gap = this->GetGeometry()->GetHelixVerticalGap(); + + double new_angle(0); + double new_dy(0); + + for (int i = 0; i < angles.size(); i++) + { + new_angle = angles[i] + angular_gap * 0.5; + new_dy = dy[i] + vertical_gap * 0.5; + new_geometry->AddProjectionInRadians(sid, sdd, new_angle, 0, new_dy, 0, 0, 0, new_dy); + } + new_geometry->VerifyHelixParameters(); + this->SetGeometry(new_geometry); +} + + +template +void +KatsevichDerivativeImageFilter::PrintSelf(std::ostream & os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + +} // end namespace rtk + +#endif diff --git a/include/rtkKatsevichForwardBinningImageFilter.h b/include/rtkKatsevichForwardBinningImageFilter.h new file mode 100644 index 000000000..61008706e --- /dev/null +++ b/include/rtkKatsevichForwardBinningImageFilter.h @@ -0,0 +1,135 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichForwardBinningImageFilter_h +#define rtkKatsevichForwardBinningImageFilter_h + +#include + +#include +#include + +#include + +#include +#include + +namespace rtk +{ + +/** \class rtkKatsevichForwardBinningImageFilter + * \brief 3D backprojection + * + * Backprojects a stack of projection images (input 1) in a 3D volume (input 0) + * using linear interpolation according to a specified geometry. The operation + * is voxel-based, meaning that the center of each voxel is projected in the + * projection images to determine the interpolation location. + * + * \test rtkfovtest.cxx + * + * \author Simon Rit + * + * \ingroup RTK Projector + */ +template +class KatsevichForwardBinningImageFilter : public itk::ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichForwardBinningImageFilter); + + /** Standard class type alias. */ + using Self = KatsevichForwardBinningImageFilter; + using OutputImageType = TOutputImage; + using InputImageType = TInputImage; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + using InputPixelType = typename TInputImage::PixelType; + using InternalInputPixelType = typename TInputImage::InternalPixelType; + using OutputImageRegionType = typename TOutputImage::RegionType; + using OutputImageIndexType = typename TOutputImage::IndexType; + + + using GeometryType = rtk::ThreeDHelicalProjectionGeometry; + using GeometryPointer = typename GeometryType::Pointer; + using ProjectionMatrixType = typename GeometryType::MatrixType; + using ProjectionImageType = itk::Image; + using ProjectionImagePointer = typename ProjectionImageType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(KatsevichForwardBinningImageFilter, itk::ImageToImageFilter); + + /** Get / Set the object pointer to projection geometry */ + itkGetMacro(Geometry, GeometryPointer); + itkSetMacro(Geometry, GeometryPointer); + + /** Get / Set the transpose flag for 2D projections (optimization trick) */ + itkGetMacro(Transpose, bool); + itkSetMacro(Transpose, bool); + +protected: + KatsevichForwardBinningImageFilter() = default; + ~KatsevichForwardBinningImageFilter() override = default; + + /** Set up requested region for inputs*/ + void + GenerateInputRequestedRegion() override; + + void + GenerateOutputInformation() override; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + + /** The two inputs should not be in the same space so there is nothing + * to verify. */ + void + VerifyInputInformation() const override + {} + + /** The input is a stack of projections, we need to interpolate in one projection + for efficiency during interpolation. Use of itk::ExtractImageFilter is + not threadsafe in ThreadedGenerateData, this one is. The output can be multiplied by a constant. + The function is templated to allow getting an itk::CudaImage. */ + template + typename TProjectionImage::Pointer + GetProjection(const unsigned int iProj); + + /** RTK geometry object */ + GeometryPointer m_Geometry; + +private: + /** Flip projection flag: infludences GetProjection and + GetIndexToIndexProjectionMatrix for optimization */ + bool m_Transpose{ false }; +}; + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkKatsevichForwardBinningImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkKatsevichForwardBinningImageFilter.hxx b/include/rtkKatsevichForwardBinningImageFilter.hxx new file mode 100644 index 000000000..f8bad60e3 --- /dev/null +++ b/include/rtkKatsevichForwardBinningImageFilter.hxx @@ -0,0 +1,232 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkKatsevichForwardBinningImageFilter_hxx +#define rtkKatsevichForwardBinningImageFilter_hxx + +#include "math.h" + +#include + +#include +#include +#include +#include + +#include + +namespace rtk +{ + +template +void +KatsevichForwardBinningImageFilter::GenerateInputRequestedRegion() +{ + // // Input 0 is Output ie result of forward binning + typename Superclass::InputImagePointer inputPtr0 = const_cast(this->GetInput()); + if (!inputPtr0) + return; + inputPtr0->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); +} +// Not necessary as GetInput()->GetLargest = inputPtr0->GetRequested before assignment. + +template +void +KatsevichForwardBinningImageFilter::GenerateOutputInformation() +{ + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); + + // Get pointers to the input and output + const InputImageType * inputPtr = this->GetInput(); + OutputImageType * outputPtr = this->GetOutput(); + + itkAssertInDebugAndIgnoreInReleaseMacro(inputPtr); + itkAssertInDebugAndIgnoreInReleaseMacro(outputPtr != nullptr); + + unsigned int i; + const typename TInputImage::SpacingType & inputSpacing = inputPtr->GetSpacing(); + const typename TInputImage::PointType & inputOrigin = inputPtr->GetOrigin(); + const typename TInputImage::SizeType & inputSize = inputPtr->GetLargestPossibleRegion().GetSize(); + const typename TInputImage::IndexType & inputStartIndex = inputPtr->GetLargestPossibleRegion().GetIndex(); + + typename TOutputImage::SpacingType outputSpacing(inputSpacing); + typename TOutputImage::SizeType outputSize(inputSize); + typename TOutputImage::PointType outputOrigin(inputOrigin); + typename TOutputImage::IndexType outputStartIndex; + + outputSize[1] = 2 * inputSize[1] + 1; + + if (this->m_Geometry->GetRadiusCylindricalDetector() == 0.) + { // Flat panel detector + + + // Compute sampling in psi + const double D = m_Geometry->GetHelixSourceToDetectorDist(); + const double R = m_Geometry->GetHelixRadius(); + const double alpha_m = atan(0.5 * inputSpacing[0] * inputSize[0] / D); + const double r = R * sin(alpha_m); + const double psi_spacing = (M_PI + 2 * alpha_m) / (outputSize[1] - 1); + + outputSpacing[1] = psi_spacing; + outputOrigin[1] = -0.5 * psi_spacing * (outputSize[1] - 1); + + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetIndex(inputStartIndex); + outputLargestPossibleRegion.SetSize(outputSize); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + outputPtr->SetOrigin(outputOrigin); + outputPtr->SetSpacing(outputSpacing); + + } + else + { // Cylindrical detector + + double radius = this->m_Geometry->GetRadiusCylindricalDetector(); + // Compute sampling in psi + double alpha_m = 0.5 * inputSize[0] * inputSpacing[0] / radius; + double psi_spacing = (M_PI + 2 * alpha_m) / (outputSize[1] - 1); + outputSpacing[1] = psi_spacing; + outputOrigin[1] = -0.5 * psi_spacing * (outputSize[1] - 1); + + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetIndex(inputStartIndex); + outputLargestPossibleRegion.SetSize(outputSize); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + outputPtr->SetSpacing(outputSpacing); + outputPtr->SetOrigin(outputOrigin); + } +} + +template +void +KatsevichForwardBinningImageFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull() || !this->m_Geometry->GetTheGeometryIsVerified()) + itkExceptionMacro(<< "Geometry has not been set or not been checked"); +} + + +template +void +KatsevichForwardBinningImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + + const unsigned int Dimension = TInputImage::ImageDimension; + const unsigned int nProj = this->GetInput()->GetLargestPossibleRegion().GetSize(Dimension - 1); + + // Create interpolator, could be any interpolation + using InterpolatorType = itk::LinearInterpolateImageFunction; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + interpolator->SetInputImage(this->GetInput()); + + using OutputRegionIterator = itk::ImageRegionIteratorWithIndex; + OutputRegionIterator itOut(this->GetOutput(), outputRegionForThread); + + // Continuous index at which we interpolate + itk::ContinuousIndex pointProj; + + const typename TOutputImage::SpacingType & outputSpacing = this->GetOutput()->GetSpacing(); + const typename TOutputImage::PointType & outputOrigin = this->GetOutput()->GetOrigin(); + const typename TOutputImage::SizeType & outputSize = this->GetOutput()->GetLargestPossibleRegion().GetSize(); + const typename TOutputImage::IndexType & outputStartIndex = this->GetOutput()->GetLargestPossibleRegion().GetIndex(); + + const typename TInputImage::SpacingType & inputSpacing = this->GetInput()->GetSpacing(); + const typename TInputImage::PointType & inputOrigin = this->GetInput()->GetOrigin(); + const typename TInputImage::SizeType & inputSize = this->GetInput()->GetLargestPossibleRegion().GetSize(); + + + // Parameters to compute psi and v + const double R = m_Geometry->GetHelixRadius(); + const double D = m_Geometry->GetHelixSourceToDetectorDist(); + const double alpha_m = atan(0.5 * inputSpacing[0] * inputSize[0] / D); + const double r = R * sin(alpha_m); + + const int L = this->GetOutput()->GetLargestPossibleRegion().GetSize(1); + const double delta_psi = (M_PI + 2 * alpha_m) / (L - 1); + const double P = m_Geometry->GetHelixPitch(); + + + for (itOut.GoToBegin(); !itOut.IsAtEnd(); ++itOut) + { + if (this->m_Geometry->GetRadiusCylindricalDetector() == 0.) // Flat panel + { + typename TOutputImage::IndexType index = itOut.GetIndex(); + double psi = outputOrigin[1] + index[1] * outputSpacing[1]; + double u = outputOrigin[0] + index[0] * outputSpacing[0]; + double v = D * P / (2 * M_PI * R); + if (psi != 0.) + { + v *= (psi + (psi / tan(psi)) * u / D); + } + else + { + v *= u / D; + } + itk::ContinuousIndex point; + point[0] = index[0]; + point[1] = (v - inputOrigin[1]) / inputSpacing[1]; + point[2] = index[2]; + + // Interpolate + if (interpolator->IsInsideBuffer(point)) + { + itOut.Set(interpolator->EvaluateAtContinuousIndex(point)); + } + } + else // Curved detector + { + + double radius = this->m_Geometry->GetRadiusCylindricalDetector(); + if (radius != D) + { + itkExceptionMacro(<< "Implementation currently only handles cylindrical detector centered on source (i.e. " + "radius = source to det distance"); + } + typename TOutputImage::IndexType index = itOut.GetIndex(); + double alpha = (outputOrigin[0] + index[0] * outputSpacing[0]) / radius; + double psi = outputOrigin[1] + index[1] * outputSpacing[1]; + double v = D * P / (2 * M_PI * R); + if (psi != 0.) + { + v *= (psi * cos(alpha) + (psi / tan(psi)) * sin(alpha)); + } + else + { + v *= sin(alpha); + } + itk::ContinuousIndex point; + point[0] = index[0]; + point[1] = (v - inputOrigin[1]) / inputSpacing[1]; + point[2] = index[2]; + + // Interpolate + if (interpolator->IsInsideBuffer(point)) + { + itOut.Set(interpolator->EvaluateAtContinuousIndex(point)); + } + } + } +} + +} // end namespace rtk + +#endif diff --git a/include/rtkPILineImageFilter.h b/include/rtkPILineImageFilter.h new file mode 100644 index 000000000..95c59ca9e --- /dev/null +++ b/include/rtkPILineImageFilter.h @@ -0,0 +1,134 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkPILineImageFilter_h +#define rtkPILineImageFilter_h + +#include + +#include +#include + +#include +#include +#include "itkImageRegionIteratorWithIndex.h" +#include + +#include +#include + +namespace rtk +{ + +/** \class PILineImageFilter + * \brief + * + * Computes the extremities of the PI-Line for each voxel of the input image. + * The lower bound is computed with a Powell-based minimization. + * + * \test + * + * \author Jérome Lesaint + * + * \ingroup + */ +template +class PILineImageFilter : public itk::ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(PILineImageFilter); + + /** Standard class type alias. */ + using Self = PILineImageFilter; + using OutputImageType = TOutputImage; + using InputImageType = TInputImage; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + using InternalInputPixelType = typename TInputImage::InternalPixelType; + using OutputImageRegionType = typename TOutputImage::RegionType; + using OutputImageIndexType = typename TOutputImage::IndexType; + using OutputImageSizeType = typename TOutputImage::SizeType; + using InputImageRegionType = typename TInputImage::RegionType; + using InputImageIndexType = typename TInputImage::IndexType; + using InputImageSizeType = typename TInputImage::SizeType; + using InputIteratorType = itk::ImageRegionIteratorWithIndex; + using OutputIteratorType = itk::ImageRegionIteratorWithIndex; + using ConstantImageType = rtk::ConstantImageSource; + + + using GeometryType = rtk::ThreeDHelicalProjectionGeometry; + using GeometryPointer = GeometryType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(PILineImageFilter, itk::ImageToImageFilter); + + /** Get / Set the object pointer to projection geometry */ + itkGetModifiableObjectMacro(Geometry, GeometryType); + itkSetObjectMacro(Geometry, GeometryType); + + /** Get / Set the transpose flag for 2D projections (optimization trick) */ + itkGetMacro(Transpose, bool); + itkSetMacro(Transpose, bool); + +protected: + PILineImageFilter() = default; + ~PILineImageFilter() override = default; + + /** Checks that inputs are correctly set. */ + void + VerifyPreconditions() ITKv5_CONST override; + + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; + + /** The two inputs should not be in the same space so there is nothing + * to verify. */ + void + VerifyInputInformation() const override + {} + + /** The input is a stack of projections, we need to interpolate in one projection + for efficiency during interpolation. Use of itk::ExtractImageFilter is + not threadsafe in ThreadedGenerateData, this one is. The output can be multiplied by a constant. + The function is templated to allow getting an itk::CudaImage. */ + template + typename TProjectionImage::Pointer + GetProjection(const unsigned int iProj); + + /** RTK geometry object */ + GeometryPointer m_Geometry; + +private: + /** Flip projection flag: infludences GetProjection and + GetIndexToIndexProjectionMatrix for optimization */ + bool m_Transpose{ false }; +}; + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkPILineImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkPILineImageFilter.hxx b/include/rtkPILineImageFilter.hxx new file mode 100644 index 000000000..c5f31acfa --- /dev/null +++ b/include/rtkPILineImageFilter.hxx @@ -0,0 +1,215 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkPILineImageFilter_hxx +#define rtkPILineImageFilter_hxx + +#include "math.h" + +#include + +#include + +#include +#include +#include +#include + +#include + +namespace rtk +{ + +int POWELL_CALLS_TO_GET_VALUE = 0; + +class PowellPILineCostFunction : public itk::SingleValuedCostFunction +{ +public: + using Self = PowellPILineCostFunction; + using Superclass = itk::SingleValuedCostFunction; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + itkNewMacro(Self); + itkTypeMacro(PowellPILineCostFunction, SingleValuedCostFunction); + + /** Get and Set macro*/ + itkGetMacro(DistanceToOrigin, double); + itkSetMacro(DistanceToOrigin, double); + + itkGetMacro(Gamma, double); + itkSetMacro(Gamma, double); + + itkGetMacro(HelixPitch, double); + itkSetMacro(HelixPitch, double); + + itkGetMacro(HelixRadius, double); + itkSetMacro(HelixRadius, double); + + itkGetMacro(AxialPosition, double); + itkSetMacro(AxialPosition, double); + + enum + { + SpaceDimension = 1 + }; + + using ParametersType = Superclass::ParametersType; + using DerivativeType = Superclass::DerivativeType; + using MeasureType = Superclass::MeasureType; + +protected: + PowellPILineCostFunction(); + + + void + GetDerivative(const ParametersType &, DerivativeType &) const override + {} + + MeasureType + GetValue(const ParametersType & parameters) const override + { + ++POWELL_CALLS_TO_GET_VALUE; + + double s = parameters[0]; + + double tmp1 = m_HelixRadius - m_DistanceToOrigin * cos(m_Gamma - s); + double alpha = atan2(m_DistanceToOrigin * sin(m_Gamma - s), tmp1); + double tmp3 = (1 + (pow(m_DistanceToOrigin, 2) - pow(m_HelixRadius, 2)) / (2 * m_HelixRadius * tmp1)); + + MeasureType measure = pow(m_AxialPosition - m_HelixPitch * ((M_PI - 2 * alpha) * tmp3 + s), 2); + + return measure; + } + + unsigned int + GetNumberOfParameters() const override + { + return SpaceDimension; + } + +private: + double m_DistanceToOrigin; + double m_Gamma; + double m_HelixPitch; + double m_HelixRadius; + double m_AxialPosition; +}; + +PowellPILineCostFunction::PowellPILineCostFunction() {} + +template +void +PILineImageFilter::VerifyPreconditions() ITKv5_CONST +{ + this->Superclass::VerifyPreconditions(); + + if (this->m_Geometry.IsNull() || !this->m_Geometry->GetTheGeometryIsVerified()) + itkExceptionMacro(<< "Geometry has not been set or not been verified"); +} + + +template +void +PILineImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + const unsigned int Dimension = TInputImage::ImageDimension; + + // Iterators + OutputIteratorType itOut(this->GetOutput(), outputRegionForThread); + InputIteratorType itIn(this->GetInput(), outputRegionForThread); + + typename InputImageType::SpacingType spacing = this->GetInput()->GetSpacing(); + typename InputImageType::PointType origin = this->GetInput()->GetOrigin(); + + for (itOut.GoToBegin(), itIn.GoToBegin(); !itOut.IsAtEnd(); ++itOut, ++itIn) + { + OutputImageIndexType index = itOut.GetIndex(); + double x = origin[0] + spacing[0] * index[0]; + double y = origin[1] + spacing[1] * index[1]; + double z = origin[2] + spacing[2] * index[2]; + + double r = sqrt(pow(x, 2) + pow(z, 2)); + double gamma = atan2(x, z); + double h = this->GetGeometry()->GetHelixPitch() / (2 * M_PI); + double R = this->GetGeometry()->GetHelixRadius(); + + using OptimizerType = itk::PowellOptimizer; + + // Declaration of an itkOptimizer + OptimizerType::Pointer itkOptimizer = OptimizerType::New(); + + + // Declaration of the CostFunction + PowellPILineCostFunction::Pointer cost_fun = PowellPILineCostFunction::New(); + cost_fun->SetDistanceToOrigin(r); + cost_fun->SetGamma(gamma); + cost_fun->SetHelixPitch(h); + cost_fun->SetHelixRadius(R); + cost_fun->SetAxialPosition(y); + + + itkOptimizer->SetCostFunction(cost_fun); + + + using ParametersType = PowellPILineCostFunction::ParametersType; + + // We start at 0 + ParametersType initialPosition(1); // 1D minimization + + initialPosition[0] = 0.; + + itkOptimizer->SetMaximize(false); + itkOptimizer->SetStepLength(1); + itkOptimizer->SetStepTolerance(0.001); + itkOptimizer->SetValueTolerance(0.01); + itkOptimizer->SetMaximumIteration(1000); + + itkOptimizer->SetInitialPosition(initialPosition); + + try + { + itkOptimizer->StartOptimization(); + } + catch (const itk::ExceptionObject & e) + { + std::cout << "Exception thrown ! " << std::endl; + std::cout << "An error occurred during Optimization" << std::endl; + std::cout << "Location = " << e.GetLocation() << std::endl; + std::cout << "Description = " << e.GetDescription() << std::endl; + // return EXIT_FAILURE; + } + + double finalPosition = itkOptimizer->GetCurrentPosition()[0]; + + // We set the two bounds of the PI Line + using VectorType = itk::Vector; + VectorType vector; + // Lower bound + vector[0] = finalPosition; + // Upper bound + double tmp = atan2(r * sin(gamma - finalPosition), R - r * cos(gamma - finalPosition)); + vector[1] = finalPosition + M_PI - 2 * tmp; + + itOut.Set(vector); + } +} + +} // end namespace rtk + +#endif diff --git a/include/rtkThreeDCircularProjectionGeometry.h b/include/rtkThreeDCircularProjectionGeometry.h index fc289acbb..a2e505018 100644 --- a/include/rtkThreeDCircularProjectionGeometry.h +++ b/include/rtkThreeDCircularProjectionGeometry.h @@ -106,7 +106,7 @@ class RTK_EXPORT ThreeDCircularProjectionGeometry : public ProjectionGeometry<3> * @param detectorRowVector absolute direction vector indicating the * orientation of the detector's rows r (sometimes referred to as v1) * @param detectorColumnVector absolute direction vector indicating the - * orientation of the detector's columns c (sometimes referred to as v2) + * orientation of the detector's columns c (sometimes refered to as v2) * @return TRUE if the projection could be added to the RTK projections list */ bool diff --git a/include/rtkThreeDHelicalProjectionGeometry.h b/include/rtkThreeDHelicalProjectionGeometry.h new file mode 100644 index 000000000..a053fd79a --- /dev/null +++ b/include/rtkThreeDHelicalProjectionGeometry.h @@ -0,0 +1,162 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkThreeDHelicalProjectionGeometry_h +#define rtkThreeDHelicalProjectionGeometry_h + +#include "RTKExport.h" +#include "rtkThreeDCircularProjectionGeometry.h" + +namespace rtk +{ +/** \class ThreeDHelicalProjectionGeometry + * \brief Projection geometry for a source and a 2-D flat panel. + * + * The source and the detector rotate around a helix paremeterized + * with the SourceToDetectorDistance and the SourceToIsocenterDistance and a helical pitch. + * The position of each projection along the helix is parameterized + * by the HelicalAngle. + * The detector is assumed to have no horizontal offset (ProjectionOffsetX = 0). It can be also rotated with + * InPlaneAngles and OutOfPlaneAngles. All angles are in radians except for the function AddProjection that takes angles + * in degrees. The source is assumed to have no horizontal offset (SourceOffsetX = 0). + * + * If SDD equals 0., then one is dealing with a parallel geometry. + * If m_RadiusCylindricalDetector != 0 : the geometry assumes a cylindrical detector. + * + * More information is provided in \ref DocGeo3D. + * + * \author Jerome Lesaint + * + * \ingroup RTK ProjectionGeometry + */ + +class RTK_EXPORT ThreeDHelicalProjectionGeometry : public ThreeDCircularProjectionGeometry +{ +public: +#if ITK_VERSION_MAJOR == 5 && ITK_VERSION_MINOR == 1 + ITK_DISALLOW_COPY_AND_ASSIGN(ThreeDHelicalProjectionGeometry); +#else + ITK_DISALLOW_COPY_AND_MOVE(ThreeDHelicalProjectionGeometry); +#endif + + using Self = ThreeDHelicalProjectionGeometry; + using Superclass = ThreeDCircularProjectionGeometry; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + using VectorType = itk::Vector; + using HomogeneousVectorType = itk::Vector; + using TwoDHomogeneousMatrixType = itk::Matrix; + using ThreeDHomogeneousMatrixType = itk::Matrix; + using PointType = itk::Point; + using Matrix3x3Type = itk::Matrix; + using HomogeneousProjectionMatrixType = Superclass::MatrixType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Add projection to geometry. One projection is defined with the rotation + * angle in degrees and the in-plane translation of the detector in physical + * units (e.g. mm). The rotation axis is assumed to be (0,1,0). + */ + + /** Only the AddProjection in radians is modified to populate the helical angle vector. */ + virtual void + AddProjectionInRadians(const double sid, + const double sdd, + const double gantryAngle, + const double projOffsetX = 0., + const double projOffsetY = 0., + const double outOfPlaneAngle = 0., + const double inPlaneAngle = 0., + const double sourceOffsetX = 0., + const double sourceOffsetY = 0.) override; + + + /** Empty the geometry object. */ + void + Clear() override; + + /** Get the vector of geometry parameters (one per projection). Angles are + * in radians.*/ + /** Get the vector of geometry parameters (one per projection). Angles are + * in radians.*/ + const std::vector & + GetHelicalAngles() const + { + return this->m_HelicalAngles; + } + + const double & + GetHelixPitch() const + { + return this->m_HelixPitch; + } + + const double & + GetHelixVerticalGap() const + { + return this->m_HelixVerticalGap; + } + + const double & + GetHelixAngularGap() const + { + return this->m_HelixAngularGap; + } + + const double & + GetHelixRadius() const + { + return this->m_HelixRadius; + } + + const double & + GetHelixSourceToDetectorDist() const + { + return this->m_HelixSourceToDetectorDist; + } + + const bool & + GetTheGeometryIsVerified() const + { + return this->m_TheGeometryIsVerified; + } + + bool + VerifyHelixParameters(); + +protected: + ThreeDHelicalProjectionGeometry(); + ~ThreeDHelicalProjectionGeometry() override = default; + + /** Some additional helix-specific global parameteres. + * Helical angles are NOT converted between 0 and 2*pi. + */ + double m_HelixRadius; + double m_HelixSourceToDetectorDist; + double m_HelixVerticalGap; + double m_HelixAngularGap; + double m_HelixPitch; + std::vector m_HelicalAngles; + bool m_TheGeometryIsVerified; +}; +} // namespace rtk + + +#endif // __rtkThreeDHelicalProjectionGeometry_h diff --git a/include/rtkThreeDHelicalProjectionGeometryXMLFileReader.h b/include/rtkThreeDHelicalProjectionGeometryXMLFileReader.h new file mode 100644 index 000000000..fabbf5134 --- /dev/null +++ b/include/rtkThreeDHelicalProjectionGeometryXMLFileReader.h @@ -0,0 +1,127 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkThreeDHelicalProjectionGeometryXMLFileReader_h +#define rtkThreeDHelicalProjectionGeometryXMLFileReader_h + +#ifdef _MSC_VER +# pragma warning(disable : 4786) +#endif + +#include +#include "rtkThreeDHelicalProjectionGeometry.h" +#include "RTKExport.h" + +namespace rtk +{ + +/** \class ThreeDHelicalProjectionGeometryXMLFileReader + * + * Reads an XML-format file containing geometry for reconstruction + * + * \test rtkgeometryfiletest.cxx, rtkvariantest.cxx, rtkxradtest.cxx, + * rtkdigisenstest.cxx, rtkelektatest.cxx + * + * \author Simon Rit + * + * \ingroup RTK IOFilters + */ +class RTK_EXPORT ThreeDHelicalProjectionGeometryXMLFileReader : public itk::XMLReader +{ +public: +#if ITK_VERSION_MAJOR == 5 && ITK_VERSION_MINOR == 1 + ITK_DISALLOW_COPY_AND_ASSIGN(ThreeDHelicalProjectionGeometryXMLFileReader); +#else + ITK_DISALLOW_COPY_AND_MOVE(ThreeDHelicalProjectionGeometryXMLFileReader); +#endif + + /** Standard type alias */ + using Self = ThreeDHelicalProjectionGeometryXMLFileReader; + using Superclass = itk::XMLReader; + using Pointer = itk::SmartPointer; + + /** Convenient type alias */ + using GeometryType = ThreeDHelicalProjectionGeometry; + using GeometryPointer = GeometryType::Pointer; + + /** Latest version */ + static const unsigned int CurrentVersion = 3; + + /** Run-time type information (and related methods). */ + itkTypeMacro(ThreeDHelicalProjectionGeometryXMLFileReader, itk::XMLFileReader); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Determine if a file can be read */ + int + CanReadFile(const char * name) override; + + /** Get smart pointer to projection geometry. */ + itkGetModifiableObjectMacro(Geometry, GeometryType); + +protected: + ThreeDHelicalProjectionGeometryXMLFileReader(); + ~ThreeDHelicalProjectionGeometryXMLFileReader() override = default; + + /** Callback function -- called from XML parser with start-of-element + * information. + */ + void + StartElement(const char * name, const char ** atts) override; + + void + StartElement(const char * name); + + void + EndElement(const char * name) override; + + void + CharacterDataHandler(const char * inData, int inLength) override; + +private: + GeometryPointer m_Geometry{ GeometryType::New() }; + + std::string m_CurCharacterData{ "" }; + + /** Projection parameters */ + double m_InPlaneAngle{ 0. }; + double m_OutOfPlaneAngle{ 0. }; + double m_GantryAngle{ 0. }; + double m_HelicalAngle{ 0. }; + double m_SourceToIsocenterDistance{ 0. }; + double m_SourceOffsetX{ 0. }; + double m_SourceOffsetY{ 0. }; + double m_SourceToDetectorDistance{ 0. }; + double m_ProjectionOffsetX{ 0. }; + double m_ProjectionOffsetY{ 0. }; + double m_HelixPitch{ 0. }; + double m_HelixRadius{ 0. }; + double m_CollimationUInf{ std::numeric_limits::max() }; + double m_CollimationUSup{ std::numeric_limits::max() }; + double m_CollimationVInf{ std::numeric_limits::max() }; + double m_CollimationVSup{ std::numeric_limits::max() }; + + /** Projection matrix */ + ThreeDHelicalProjectionGeometry::MatrixType m_Matrix; + + /** File format version */ + unsigned int m_Version{ 0 }; +}; +} // namespace rtk +#endif diff --git a/include/rtkThreeDHelicalProjectionGeometryXMLFileWriter.h b/include/rtkThreeDHelicalProjectionGeometryXMLFileWriter.h new file mode 100644 index 000000000..214b13666 --- /dev/null +++ b/include/rtkThreeDHelicalProjectionGeometryXMLFileWriter.h @@ -0,0 +1,92 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef rtkThreeDHelicalProjectionGeometryXMLFileWriter_h +#define rtkThreeDHelicalProjectionGeometryXMLFileWriter_h + +#ifdef _MSC_VER +# pragma warning(disable : 4786) +#endif + +#include "RTKExport.h" +#include +#include "rtkThreeDHelicalProjectionGeometry.h" + +namespace rtk +{ + +/** \class ThreeDHelicalProjectionGeometryXMLFileWriter + * + * Writes an XML-format file containing geometry for reconstruction + * + * \author Simon Rit + * + * \ingroup RTK IOFilters + */ +class RTK_EXPORT ThreeDHelicalProjectionGeometryXMLFileWriter + : public itk::XMLWriterBase +{ +public: +#if ITK_VERSION_MAJOR == 5 && ITK_VERSION_MINOR == 1 + ITK_DISALLOW_COPY_AND_ASSIGN(ThreeDHelicalProjectionGeometryXMLFileWriter); +#else + ITK_DISALLOW_COPY_AND_MOVE(ThreeDHelicalProjectionGeometryXMLFileWriter); +#endif + + /** standard type alias */ + using Superclass = itk::XMLWriterBase; + using Self = ThreeDHelicalProjectionGeometryXMLFileWriter; + using Pointer = itk::SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ThreeDHelicalProjectionGeometryXMLFileWriter, itk::XMLWriterBase); + + /** Test whether a file is writable. */ + int + CanWriteFile(const char * name) override; + + /** Actually write out the file in question */ + int + WriteFile() override; + +protected: + ThreeDHelicalProjectionGeometryXMLFileWriter() = default; + ~ThreeDHelicalProjectionGeometryXMLFileWriter() override = default; + + /** If all values are equal in v, write first value (if not 0.) in + output file with parameter value s and return true. Return false + otherwise. + */ + bool + WriteGlobalParameter(std::ofstream & output, + const std::string & indent, + const std::vector & v, + const std::string & s, + bool convertToDegrees = false, + double defval = 0.); + + /** Write projection specific parameter with name s. */ + void + WriteLocalParameter(std::ofstream & output, const std::string & indent, const double & v, const std::string & s); +}; +} // namespace rtk + +#endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4fff512bc..031929e91 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,6 +40,9 @@ set(RTK_SRCS rtkThreeDCircularProjectionGeometry.cxx rtkThreeDCircularProjectionGeometryXMLFileReader.cxx rtkThreeDCircularProjectionGeometryXMLFileWriter.cxx + rtkThreeDHelicalProjectionGeometry.cxx + rtkThreeDHelicalProjectionGeometryXMLFileReader.cxx + rtkThreeDHelicalProjectionGeometryXMLFileWriter.cxx rtkResourceProbesCollector.cxx rtkVarianObiGeometryReader.cxx rtkVarianObiXMLFileReader.cxx diff --git a/src/rtkThreeDHelicalProjectionGeometry.cxx b/src/rtkThreeDHelicalProjectionGeometry.cxx new file mode 100644 index 000000000..d77ec4d30 --- /dev/null +++ b/src/rtkThreeDHelicalProjectionGeometry.cxx @@ -0,0 +1,249 @@ +/*========================================================================= + * + * 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 "math.h" + +#include "rtkThreeDHelicalProjectionGeometry.h" +#include "rtkMacro.h" + +#include +#include + +#include +#include + +rtk::ThreeDHelicalProjectionGeometry::ThreeDHelicalProjectionGeometry() +{ + m_HelixRadius = 0.; + m_HelixSourceToDetectorDist = 0.; + m_HelixVerticalGap = 0.; + m_HelixAngularGap = 0.; + m_HelixPitch = 0.; + m_TheGeometryIsVerified = false; +} + + +void +rtk::ThreeDHelicalProjectionGeometry::AddProjectionInRadians(const double sid, + const double sdd, + const double gantryAngle, + const double projOffsetX, + const double projOffsetY, + const double outOfPlaneAngle, + const double inPlaneAngle, + const double sourceOffsetX, + const double sourceOffsetY) +{ + // Check parallel / divergent projections consistency + if (!m_GantryAngles.empty()) + { + if (sdd == 0. && m_SourceToDetectorDistances[0] != 0.) + { + itkGenericExceptionMacro( + << "Cannot add a parallel projection in a 3D geometry object containing divergent projections"); + } + if (sdd != 0. && m_SourceToDetectorDistances[0] == 0.) + { + itkGenericExceptionMacro( + << "Cannot add a divergent projection in a 3D geometry object containing parallel projections"); + } + } + + // Detector orientation parameters + m_GantryAngles.push_back(ConvertAngleBetween0And2PIRadians(gantryAngle)); + m_HelicalAngles.push_back(gantryAngle); // No conversion mod 2*pi here for a helix. + m_OutOfPlaneAngles.push_back(ConvertAngleBetween0And2PIRadians(outOfPlaneAngle)); + m_InPlaneAngles.push_back(ConvertAngleBetween0And2PIRadians(inPlaneAngle)); + + // Source position parameters + m_SourceToIsocenterDistances.push_back(sid); + m_SourceOffsetsX.push_back(sourceOffsetX); + m_SourceOffsetsY.push_back(sourceOffsetY); + + // Detector position parameters + m_SourceToDetectorDistances.push_back(sdd); + m_ProjectionOffsetsX.push_back(projOffsetX); + m_ProjectionOffsetsY.push_back(projOffsetY); + + // Compute sub-matrices + AddProjectionTranslationMatrix( + ComputeTranslationHomogeneousMatrix(sourceOffsetX - projOffsetX, sourceOffsetY - projOffsetY)); + AddMagnificationMatrix(ComputeProjectionMagnificationMatrix(-sdd, -sid)); + AddRotationMatrix(ComputeRotationHomogeneousMatrix(-outOfPlaneAngle, -gantryAngle, -inPlaneAngle)); + AddSourceTranslationMatrix(ComputeTranslationHomogeneousMatrix(-sourceOffsetX, -sourceOffsetY, 0.)); + + Superclass::MatrixType matrix; + matrix = this->GetProjectionTranslationMatrices().back().GetVnlMatrix() * + this->GetMagnificationMatrices().back().GetVnlMatrix() * + this->GetSourceTranslationMatrices().back().GetVnlMatrix() * + this->GetRotationMatrices().back().GetVnlMatrix(); + this->AddMatrix(matrix); + + // Calculate source angle + VectorType z; + z.Fill(0.); + z[2] = 1.; + HomogeneousVectorType sph = GetSourcePosition(m_GantryAngles.size() - 1); + sph[1] = 0.; // Project position to central plane + VectorType sp(&(sph[0])); + sp.Normalize(); + double a = acos(sp * z); + if (sp[0] > 0.) + a = 2. * itk::Math::pi - a; + m_SourceAngles.push_back(ConvertAngleBetween0And2PIRadians(a)); + + // Default collimation (uncollimated) + m_CollimationUInf.push_back(std::numeric_limits::max()); + m_CollimationUSup.push_back(std::numeric_limits::max()); + m_CollimationVInf.push_back(std::numeric_limits::max()); + m_CollimationVSup.push_back(std::numeric_limits::max()); + + this->Modified(); +} + + +void +rtk::ThreeDHelicalProjectionGeometry::Clear() +{ + Superclass::Clear(); + + m_GantryAngles.clear(); + m_HelicalAngles.clear(); + m_OutOfPlaneAngles.clear(); + m_InPlaneAngles.clear(); + m_SourceAngles.clear(); + m_SourceToIsocenterDistances.clear(); + m_SourceOffsetsX.clear(); + m_SourceOffsetsY.clear(); + m_SourceToDetectorDistances.clear(); + m_ProjectionOffsetsX.clear(); + m_ProjectionOffsetsY.clear(); + m_CollimationUInf.clear(); + m_CollimationUSup.clear(); + m_CollimationVInf.clear(); + m_CollimationVSup.clear(); + + m_ProjectionTranslationMatrices.clear(); + m_MagnificationMatrices.clear(); + m_RotationMatrices.clear(); + m_SourceTranslationMatrices.clear(); + this->Modified(); +} + +bool +rtk::ThreeDHelicalProjectionGeometry::VerifyHelixParameters() +{ + // Check that sid is constant + std::vector v = this->GetSourceToIsocenterDistances(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "Helical Traj. : SourceToIsocenterDistance must be constant."); + return false; + } + m_HelixRadius = v[0]; + + // Check that sdd is constant + v = this->GetSourceToDetectorDistances(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "(Helical Traj. : SourceToDetectorDistance must be constant."); + return false; + } + m_HelixSourceToDetectorDist = v[0]; + + // Check that DetectorOffsetX is constant + v = this->GetProjectionOffsetsX(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "Helical Traj. : ProjectionOffsetX must be constant"); + return false; + } + // Check that SourceOffsetX is constant + v = this->GetSourceOffsetsX(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "Helical Traj. : SourceOffsetX must be constant"); + return false; + } + + // Check that SourceOffsetsY and DetectorOffsetsY are equal + v = this->GetProjectionOffsetsY(); + std::vector w = this->GetSourceOffsetsY(); + for (size_t i = 0; i < w.size(); i++) + if (!itk::Math::AlmostEquals(v[i], w[i])) + { + itkGenericExceptionMacro(<< "Helical Traj. : DetectorOffsetsY must all be equal to SourceOffsetsY"); + return false; + } + + // Check that vertical gaps are constant + for (size_t i = 1; i < v.size() - 1; i++) + if (!itk::Math::FloatAlmostEqual(v[i + 1] - v[i], v[i] - v[i - 1], 4, 1e-10)) + { + itkGenericExceptionMacro(<< "Helical Traj. : Vertical gaps must all be equal"); + return false; + } + m_HelixVerticalGap = v[1] - v[0]; + + // Check that InPlaneAnles are al ZERO + v = this->GetInPlaneAngles(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "Helical Traj. : InPlaneAngles must all be equal to ZERO"); + return false; + } + if (!itk::Math::AlmostEquals(v[0], 0.)) + { + itkGenericExceptionMacro(<< "Helical Traj. : InPlaneAngles must all be equal to ZERO"); + return false; + } + // Check that OutOfPlaneAngles are all ZERO + v = this->GetOutOfPlaneAngles(); + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + { + itkGenericExceptionMacro(<< "Helical Traj. : OutOfPlaneAngles must all be equal to ZERO"); + return false; + } + if (!itk::Math::AlmostEquals(v[0], 0.)) + { + itkGenericExceptionMacro(<< "Helical Traj. : OutOfPlaneAngles must all be equal to ZERO"); + return false; + } + + // Check that angluar gap is constant + v = this->GetHelicalAngles(); + for (size_t i = 1; i < v.size(); i++) + if (!itk::Math::FloatAlmostEqual(v[i] - v[i - 1], v[1] - v[0], 4, 1e-8)) + { + std::cout << "Gap 1 : " << v[i] - v[i - 1] << " Gap 2 : " << v[i + 1] - v[i] << std::endl; + itkGenericExceptionMacro(<< "Helical Traj. : Angular gaps must be constant"); + return false; + } + m_HelixAngularGap = v[1] - v[0]; + m_HelixPitch = m_HelixVerticalGap / m_HelixAngularGap * (2 * itk::Math::pi); + + m_TheGeometryIsVerified = true; + + return true; +} diff --git a/src/rtkThreeDHelicalProjectionGeometryXMLFileReader.cxx b/src/rtkThreeDHelicalProjectionGeometryXMLFileReader.cxx new file mode 100644 index 000000000..499dc5774 --- /dev/null +++ b/src/rtkThreeDHelicalProjectionGeometryXMLFileReader.cxx @@ -0,0 +1,186 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef _rtkThreeDHelicalProjectionGeometryXMLFileReader_cxx +#define _rtkThreeDHelicalProjectionGeometryXMLFileReader_cxx + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" + +#include +#include +#include + +#include + +namespace rtk +{ + +ThreeDHelicalProjectionGeometryXMLFileReader::ThreeDHelicalProjectionGeometryXMLFileReader() +{ + this->m_OutputObject = &(*m_Geometry); +} + +int +ThreeDHelicalProjectionGeometryXMLFileReader::CanReadFile(const char * name) +{ + if (!itksys::SystemTools::FileExists(name) || itksys::SystemTools::FileIsDirectory(name) || + itksys::SystemTools::FileLength(name) == 0) + return 0; + return 1; +} + +void +ThreeDHelicalProjectionGeometryXMLFileReader::StartElement(const char * name, const char ** atts) +{ + m_CurCharacterData = ""; + this->StartElement(name); + + // Check on last version of file format. Warning if not. + if (std::string(name) == "RTKThreeDHelicalGeometry") + { + while ((*atts) != nullptr) + { + if (std::string(atts[0]) == "version") + m_Version = atoi(atts[1]); + atts += 2; + } + // Version 3 is backward compatible with version 2 + if (m_Version != this->CurrentVersion && !(m_Version == 2 && this->CurrentVersion == 3)) + itkGenericExceptionMacro(<< "Incompatible version of input geometry (v" << m_Version + << ") with current geometry (v" << this->CurrentVersion + << "). You must re-generate your geometry file again."); + this->m_OutputObject->Clear(); + } +} + +void +ThreeDHelicalProjectionGeometryXMLFileReader::StartElement(const char * itkNotUsed(name)) +{} + +void +ThreeDHelicalProjectionGeometryXMLFileReader::EndElement(const char * name) +{ + if (itksys::SystemTools::Strucmp(name, "InPlaneAngle") == 0) + m_InPlaneAngle = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "GantryAngle") == 0 || + itksys::SystemTools::Strucmp(name, "Angle") == 0) // Second one for backward compatibility + m_GantryAngle = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "HelicalAngle") == 0) + m_HelicalAngle = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "HelixPitch") == 0) + m_HelixPitch = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "Helixradius") == 0) + m_HelixRadius = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "OutOfPlaneAngle") == 0) + m_OutOfPlaneAngle = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "SourceToIsocenterDistance") == 0) + m_SourceToIsocenterDistance = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "SourceOffsetX") == 0) + m_SourceOffsetX = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "SourceOffsetY") == 0) + m_SourceOffsetY = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "SourceToDetectorDistance") == 0) + m_SourceToDetectorDistance = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "ProjectionOffsetX") == 0) + m_ProjectionOffsetX = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "ProjectionOffsetY") == 0) + m_ProjectionOffsetY = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "RadiusCylindricalDetector") == 0) + { + double radiusCylindricalDetector = atof(this->m_CurCharacterData.c_str()); + this->m_OutputObject->SetRadiusCylindricalDetector(radiusCylindricalDetector); + } + + if (itksys::SystemTools::Strucmp(name, "CollimationUInf") == 0) + m_CollimationUInf = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "CollimationUSup") == 0) + m_CollimationUSup = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "CollimationVInf") == 0) + m_CollimationVInf = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "CollimationVSup") == 0) + m_CollimationVSup = atof(this->m_CurCharacterData.c_str()); + + if (itksys::SystemTools::Strucmp(name, "Matrix") == 0) + { + std::istringstream iss(this->m_CurCharacterData); + double value = 0.; + for (unsigned int i = 0; i < m_Matrix.RowDimensions; i++) + for (unsigned int j = 0; j < m_Matrix.ColumnDimensions; j++) + { + iss >> value; + m_Matrix[i][j] = value; + } + } + + if (itksys::SystemTools::Strucmp(name, "Projection") == 0) + { + this->m_OutputObject->AddProjection(m_SourceToIsocenterDistance, + m_SourceToDetectorDistance, + m_HelicalAngle, + m_ProjectionOffsetX, + m_ProjectionOffsetY, + m_OutOfPlaneAngle, + m_InPlaneAngle, + m_SourceOffsetX, + m_SourceOffsetY); + + this->m_OutputObject->SetCollimationOfLastProjection( + m_CollimationUInf, m_CollimationUSup, m_CollimationVInf, m_CollimationVSup); + + for (unsigned int i = 0; i < m_Matrix.RowDimensions; i++) + for (unsigned int j = 0; j < m_Matrix.ColumnDimensions; j++) + { + // Tolerance can not be vcl_numeric_limits::epsilon(), too strict + // 0.001 is a random choice to catch "large" inconsistencies + if (fabs(m_Matrix[i][j] - m_OutputObject->GetMatrices().back()[i][j]) > 0.001) + { + itkGenericExceptionMacro(<< "Matrix and parameters are not consistent." << std::endl + << "Read matrix from geometry file: " << std::endl + << m_Matrix << std::endl + << "Computed matrix from parameters:" << std::endl + << m_OutputObject->GetMatrices().back()); + } + } + } +} + +void +ThreeDHelicalProjectionGeometryXMLFileReader::CharacterDataHandler(const char * inData, int inLength) +{ + for (int i = 0; i < inLength; i++) + m_CurCharacterData = m_CurCharacterData + inData[i]; +} + +} // namespace rtk + +#endif diff --git a/src/rtkThreeDHelicalProjectionGeometryXMLFileWriter.cxx b/src/rtkThreeDHelicalProjectionGeometryXMLFileWriter.cxx new file mode 100644 index 000000000..cfcea70a4 --- /dev/null +++ b/src/rtkThreeDHelicalProjectionGeometryXMLFileWriter.cxx @@ -0,0 +1,236 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +#ifndef _rtkThreeDHelicalProjectionGeometryXMLFileWriter_cxx +#define _rtkThreeDHelicalProjectionGeometryXMLFileWriter_cxx + +#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h" +#include "rtkThreeDHelicalProjectionGeometryXMLFileWriter.h" + +#include +#include +#include + +#include + +namespace rtk +{ + +int +ThreeDHelicalProjectionGeometryXMLFileWriter::CanWriteFile(const char * name) +{ + std::ofstream output(name); + + if (output.fail()) + return false; + return true; +} + +int +ThreeDHelicalProjectionGeometryXMLFileWriter::WriteFile() +{ + if (this->m_InputObject->GetGantryAngles().empty()) + itkGenericExceptionMacro(<< "Geometry object is empty, cannot write it"); + + this->m_InputObject->VerifyHelixParameters(); + + std::ofstream output(this->m_Filename.c_str()); + constexpr int maxDigits = 15; + + output.precision(maxDigits); + std::string indent(" "); + + this->WriteStartElement("?xml version=\"1.0\"?", output); + output << std::endl; + this->WriteStartElement("!DOCTYPE RTKGEOMETRY", output); + output << std::endl; + std::ostringstream startWithVersion; + startWithVersion << "RTKThreeDHelicalGeometry version=\"" + << ThreeDHelicalProjectionGeometryXMLFileReader::CurrentVersion << '"'; + this->WriteStartElement(startWithVersion.str().c_str(), output); + output << std::endl; + + // First, we test for each of the 9 parameters per projection if it's constant + // over all projection images except GantryAngle which is supposed to be different + // for all projections. If 0. for OutOfPlaneAngle, InPlaneAngle, projection and source + // offsets X and Y, it is not written (default value). + bool bSIDGlobal = WriteGlobalParameter( + output, indent, this->m_InputObject->GetSourceToIsocenterDistances(), "SourceToIsocenterDistance"); + bool bSDDGlobal = WriteGlobalParameter( + output, indent, this->m_InputObject->GetSourceToDetectorDistances(), "SourceToDetectorDistance"); + + // Helix specific parameters + WriteLocalParameter(output, indent, this->m_InputObject->GetHelixPitch(), "HelixPitch"); + WriteLocalParameter(output, indent, this->m_InputObject->GetHelixRadius(), "HelixRadius"); + + bool bSourceXGlobal = WriteGlobalParameter(output, indent, this->m_InputObject->GetSourceOffsetsX(), "SourceOffsetX"); + bool bSourceYGlobal = WriteGlobalParameter(output, indent, this->m_InputObject->GetSourceOffsetsY(), "SourceOffsetY"); + bool bProjXGlobal = + WriteGlobalParameter(output, indent, this->m_InputObject->GetProjectionOffsetsX(), "ProjectionOffsetX"); + bool bProjYGlobal = + WriteGlobalParameter(output, indent, this->m_InputObject->GetProjectionOffsetsY(), "ProjectionOffsetY"); + bool bInPlaneGlobal = + WriteGlobalParameter(output, indent, this->m_InputObject->GetInPlaneAngles(), "InPlaneAngle", true); + bool bOutOfPlaneGlobal = + WriteGlobalParameter(output, indent, this->m_InputObject->GetOutOfPlaneAngles(), "OutOfPlaneAngle", true); + + bool bCollimationUInf = WriteGlobalParameter(output, + indent, + this->m_InputObject->GetCollimationUInf(), + "CollimationUInf", + false, + std::numeric_limits::max()); + + bool bCollimationUSup = WriteGlobalParameter(output, + indent, + this->m_InputObject->GetCollimationUSup(), + "CollimationUSup", + false, + std::numeric_limits::max()); + + bool bCollimationVInf = WriteGlobalParameter(output, + indent, + this->m_InputObject->GetCollimationVInf(), + "CollimationVInf", + false, + std::numeric_limits::max()); + + bool bCollimationVSup = WriteGlobalParameter(output, + indent, + this->m_InputObject->GetCollimationVSup(), + "CollimationVSup", + false, + std::numeric_limits::max()); + + const double radius = this->m_InputObject->GetRadiusCylindricalDetector(); + if (0. != radius) + WriteLocalParameter(output, indent, radius, "RadiusCylindricalDetector"); + + // Second, write per projection parameters (if corresponding parameter is not global) + const double radiansToDegrees = 45. / std::atan(1.); + for (unsigned int i = 0; i < this->m_InputObject->GetMatrices().size(); i++) + { + output << indent; + this->WriteStartElement("Projection", output); + output << std::endl; + + // Only the GantryAngle and the HelicalAngle are necessarily projection specific + WriteLocalParameter(output, indent, radiansToDegrees * this->m_InputObject->GetGantryAngles()[i], "GantryAngle"); + WriteLocalParameter(output, indent, radiansToDegrees * this->m_InputObject->GetHelicalAngles()[i], "HelicalAngle"); + if (!bSIDGlobal) + WriteLocalParameter( + output, indent, this->m_InputObject->GetSourceToIsocenterDistances()[i], "SourceToIsocenterDistance"); + if (!bSDDGlobal) + WriteLocalParameter( + output, indent, this->m_InputObject->GetSourceToDetectorDistances()[i], "SourceToDetectorDistance"); + if (!bSourceXGlobal) + WriteLocalParameter(output, indent, this->m_InputObject->GetSourceOffsetsX()[i], "SourceOffsetX"); + if (!bSourceYGlobal) + WriteLocalParameter(output, indent, this->m_InputObject->GetSourceOffsetsY()[i], "SourceOffsetY"); + if (!bProjXGlobal) + WriteLocalParameter(output, indent, this->m_InputObject->GetProjectionOffsetsX()[i], "ProjectionOffsetX"); + if (!bProjYGlobal) + WriteLocalParameter(output, indent, this->m_InputObject->GetProjectionOffsetsY()[i], "ProjectionOffsetY"); + if (!bInPlaneGlobal) + WriteLocalParameter( + output, indent, radiansToDegrees * this->m_InputObject->GetInPlaneAngles()[i], "InPlaneAngle"); + if (!bOutOfPlaneGlobal) + WriteLocalParameter( + output, indent, radiansToDegrees * this->m_InputObject->GetOutOfPlaneAngles()[i], "OutOfPlaneAngle"); + + if (!bCollimationUInf) + WriteLocalParameter(output, indent, this->m_InputObject->GetCollimationUInf()[i], "CollimationUInf"); + + if (!bCollimationUSup) + WriteLocalParameter(output, indent, this->m_InputObject->GetCollimationUSup()[i], "CollimationUSup"); + + if (!bCollimationVInf) + WriteLocalParameter(output, indent, this->m_InputObject->GetCollimationVInf()[i], "CollimationVInf"); + + if (!bCollimationVSup) + WriteLocalParameter(output, indent, this->m_InputObject->GetCollimationVSup()[i], "CollimationVSup"); + + // Matrix + output << indent << indent; + this->WriteStartElement("Matrix", output); + output << std::endl; + for (unsigned int j = 0; j < 3; j++) + { + output << indent << indent << indent; + for (unsigned int k = 0; k < 4; k++) + output << std::setw(maxDigits + 4) << this->m_InputObject->GetMatrices()[i][j][k] << ' '; + output.seekp(-1, std::ios_base::cur); + output << std::endl; + } + output << indent << indent; + this->WriteEndElement("Matrix", output); + output << std::endl; + + output << indent; + this->WriteEndElement("Projection", output); + output << std::endl; + } + + this->WriteEndElement("RTKThreeDHelicalGeometry", output); + output << std::endl; + + return 0; +} + +bool +ThreeDHelicalProjectionGeometryXMLFileWriter::WriteGlobalParameter(std::ofstream & output, + const std::string & indent, + const std::vector & v, + const std::string & s, + bool convertToDegrees, + double defval) +{ + // Test if all values in vector v are equal. Return false if not. + for (size_t i = 0; i < v.size(); i++) + if (v[i] != v[0]) + return false; + + // Write value in file if not 0. + if (defval != v[0]) + { + double val = v[0]; + if (convertToDegrees) + val *= 45. / std::atan(1.); + + WriteLocalParameter(output, indent, val, s); + } + return true; +} + +void +ThreeDHelicalProjectionGeometryXMLFileWriter::WriteLocalParameter(std::ofstream & output, + const std::string & indent, + const double & v, + const std::string & s) +{ + std::string ss(s); + output << indent << indent; + this->WriteStartElement(ss, output); + output << v; + this->WriteEndElement(ss, output); + output << std::endl; +} + +} // namespace rtk + +#endif From 5fe19c2236b08390002d63718aa7c1717e82bf30 Mon Sep 17 00:00:00 2001 From: Jerome Lesaint Date: Sat, 29 May 2021 13:16:04 +0200 Subject: [PATCH 2/3] WIP: Response to PR review Point by point response to PR review. --- .gitattributes | 2 +- .gitignore | 3 --- ...kHilbertTransformOnKappaLinesImageFilter.h | 14 +++++++++- ...ilbertTransformOnKappaLinesImageFilter.hxx | 5 ---- .../rtkKatsevichBackwardBinningImageFilter.h | 14 +++------- ...rtkKatsevichBackwardBinningImageFilter.hxx | 26 ++----------------- include/rtkKatsevichDerivativeImageFilter.h | 4 ++- .../rtkKatsevichForwardBinningImageFilter.h | 12 +++------ 8 files changed, 27 insertions(+), 53 deletions(-) diff --git a/.gitattributes b/.gitattributes index 1a279742b..d702b5308 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,6 +1,6 @@ .git* export-ignore # Custom attribute to mark sources as using our C++/C code style. -[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat +[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat,uncrustify cmake/*.cxx our-c-style include/*.h* our-c-style diff --git a/.gitignore b/.gitignore index bde0b8fdd..479088e28 100644 --- a/.gitignore +++ b/.gitignore @@ -33,6 +33,3 @@ CMakeLists.txt.user* # Mac System File .DS_Store - - -.project diff --git a/include/rtkHilbertTransformOnKappaLinesImageFilter.h b/include/rtkHilbertTransformOnKappaLinesImageFilter.h index 484eda2bd..bec72a3f5 100644 --- a/include/rtkHilbertTransformOnKappaLinesImageFilter.h +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.h @@ -38,9 +38,21 @@ namespace rtk * - Apply 1D Hilbert transform on the resampled psi-direction. * - Resample the data back to the usual vertical coordinates. * + * * \dot + * digraph HilbertTransformOnKappaLinesImageFilter { + * node [shape=box]; + * 1 [ label="rtk::KatsevichForwardBinningImageFilter" URL="\ref rtk::KatsevichForwardBinningImageFilter"]; + * 2 [ label="rtk::FFTHilbertImageFilter" URL="\ref rtk::FFTHilbertImageFilter"]; + * 3 [ label="rtk::KatsevichBackwardBinningImageFilter" URL="\ref rtk::KatsevichBackwardBinningImageFilter"]; + * 1 -> 2; + * 2 -> 3; + * } + * \enddot + * + * * \author Jerome Lesaint and Alexandre Esa * - * \test katsHilbertTransformOnKappaLines.cxx, katsKatsevitchHelicalReconstructionImageFilterTets.cxx. + * \test * * \ingroup */ diff --git a/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx b/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx index efd43669a..bf538ddd8 100644 --- a/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx @@ -29,13 +29,8 @@ template void HilbertTransformOnKappaLinesImageFilter::GenerateOutputInformation() { - // Mettre a jour la taille de l'image de sortie. Ainsi que l'origine. - // Sur le modele de itk::ShrinkImageFilter // Call the superclass' implementation of this method Superclass::GenerateOutputInformation(); - - // std::cout << "(Hilbert transf) Input Size : " << this->GetInput()->GetLargestPossibleRegion().GetSize(); - // std::cout << " Inpu torigin : " << this->GetInput()->GetOrigin() << std::endl; } template diff --git a/include/rtkKatsevichBackwardBinningImageFilter.h b/include/rtkKatsevichBackwardBinningImageFilter.h index 2fc45904c..87f682f9e 100644 --- a/include/rtkKatsevichBackwardBinningImageFilter.h +++ b/include/rtkKatsevichBackwardBinningImageFilter.h @@ -37,14 +37,12 @@ namespace rtk /** \class KatsevichBackwardBinningImageFilter * \brief * - * Backprojects a stack of projection images (input 1) in a 3D volume (input 0) - * using linear interpolation according to a specified geometry. The operation - * is voxel-based, meaning that the center of each voxel is projected in the - * projection images to determine the interpolation location. + * Resample the Hilbert transformed data back to the original v-coordinate, according + * to a specified geometry. * - * \test rtkfovtest.cxx + * \test * - * \author Simon Rit + * \author Jerome Lesaint * * \ingroup RTK Projector */ @@ -96,10 +94,6 @@ class KatsevichBackwardBinningImageFilter : public itk::ImageToImageFilter -void -KatsevichBackwardBinningImageFilter::GenerateInputRequestedRegion() -{ - // // Input 0 is Output ie result of forward binning - // typename Superclass::InputImagePointer inputPtr0 = const_cast(this->GetInput()); - // std::cout << "Input largest : " << this->GetInput()->GetLargestPossibleRegion() << std::endl; - // std::cout << "inputPtr0 : " << inputPtr0->GetRequestedRegion() << std::endl; - // if (!inputPtr0) - // return; - // inputPtr0->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); - // - // //TO DELETE : SPACING AND ORIGIN ALWAYS (1,1,1) AND (0,0,0) -> NOT OKAY !!!! - // //std::cout<<"Input size: "<GetInput()->GetLargestPossibleRegion().GetSize()<GetInput()->GetSpacing()<GetInput()->GetOrigin()<GetLargest = inputPtr0->GetRequested before assignment. - template void KatsevichBackwardBinningImageFilter::GenerateOutputInformation() @@ -109,7 +88,7 @@ KatsevichBackwardBinningImageFilter::DynamicThreadedG const OutputImageRegionType & outputRegionForThread) { const unsigned int Dimension = TInputImage::ImageDimension; - + // Iterator over output projections (in v) OutputSliceIteratorType itOut(this->GetOutput(), outputRegionForThread); itOut.SetFirstDirection(1); @@ -151,8 +130,7 @@ KatsevichBackwardBinningImageFilter::DynamicThreadedG const double P = m_Geometry->GetHelixPitch(); const double R = m_Geometry->GetHelixRadius(); - const double alpha_m = - atan(0.5 * u_spacing * (det_nu - 1) / D); // Générique mais doit être cohérent avec le fwd binning. + const double alpha_m = atan(0.5 * u_spacing * det_nu / D); // Consistent with Forward binning. const double r = R * sin(alpha_m); const int L = this->GetInput()->GetLargestPossibleRegion().GetSize(1); diff --git a/include/rtkKatsevichDerivativeImageFilter.h b/include/rtkKatsevichDerivativeImageFilter.h index 113d6c3c2..b9e008a56 100644 --- a/include/rtkKatsevichDerivativeImageFilter.h +++ b/include/rtkKatsevichDerivativeImageFilter.h @@ -37,7 +37,9 @@ namespace rtk /** \class KatsevichDerivativeImageFilter * \brief Computes the derivatives of the projections wrt the projection index (the rotation angle). * - * This filter implements the derivative as described in Noo et al., PMB, 2003 + * This filter implements the derivative formula Eq. 46 (curved det) and Eq. 87 (flat panel) + * from Noo et al., PMB, 2003 + * * * \author Jerome Lesaint * diff --git a/include/rtkKatsevichForwardBinningImageFilter.h b/include/rtkKatsevichForwardBinningImageFilter.h index 61008706e..aee37704f 100644 --- a/include/rtkKatsevichForwardBinningImageFilter.h +++ b/include/rtkKatsevichForwardBinningImageFilter.h @@ -33,16 +33,12 @@ namespace rtk { /** \class rtkKatsevichForwardBinningImageFilter - * \brief 3D backprojection + * \brief Rebinning along kappa-lines * - * Backprojects a stack of projection images (input 1) in a 3D volume (input 0) - * using linear interpolation according to a specified geometry. The operation - * is voxel-based, meaning that the center of each voxel is projected in the - * projection images to determine the interpolation location. + * Rebin the projection data along kappa-lines using linear + * interpolation according to a specified helical geometry. * - * \test rtkfovtest.cxx - * - * \author Simon Rit + * \author Jerome Lesaint * * \ingroup RTK Projector */ From 2236dadd21d5d934bfa8b27548e3861dbaf70fc1 Mon Sep 17 00:00:00 2001 From: Jerome Lesaint Date: Sat, 29 May 2021 14:33:44 +0200 Subject: [PATCH 3/3] WIP: Make KatsevichBackProj inherit from rtkBackProj To avoid redundancy in the code, the rtk::KatsevichBackProjectionImageFilter inherits from rtk::BackProjectionImageFilter. WARNING : Compilation threw a couple of warning due to Geometry getter/setter (implemented in macros) whcih overrides the mother function, without saying it. But don't know how to deal with that. --- .../rtkKatsevichBackProjectionImageFilter.h | 30 +-- .../rtkKatsevichBackProjectionImageFilter.hxx | 200 +++--------------- 2 files changed, 36 insertions(+), 194 deletions(-) diff --git a/include/rtkKatsevichBackProjectionImageFilter.h b/include/rtkKatsevichBackProjectionImageFilter.h index 6e14b5b91..cabe7ca3a 100644 --- a/include/rtkKatsevichBackProjectionImageFilter.h +++ b/include/rtkKatsevichBackProjectionImageFilter.h @@ -24,6 +24,7 @@ #include #include +#include #include #include @@ -48,14 +49,14 @@ namespace rtk * \ingroup RTK Projector */ template -class KatsevichBackProjectionImageFilter : public itk::InPlaceImageFilter +class KatsevichBackProjectionImageFilter : public rtk::BackProjectionImageFilter { public: ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichBackProjectionImageFilter); /** Standard class type alias. */ using Self = KatsevichBackProjectionImageFilter; - using Superclass = itk::ImageToImageFilter; + using Superclass = rtk::BackProjectionImageFilter; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; @@ -66,7 +67,7 @@ class KatsevichBackProjectionImageFilter : public itk::InPlaceImageFilter; + using ProjectionImageType = typename Superclass::ProjectionImageType; using ProjectionImagePointer = typename ProjectionImageType::Pointer; using PILineRealType = float; using PILineImageType = itk::Image, TInputImage::ImageDimension>; @@ -91,10 +92,6 @@ class KatsevichBackProjectionImageFilter : public itk::InPlaceImageFilter - typename TProjectionImage::Pointer - GetProjection(const unsigned int iProj); - - /** Creates iProj index to index projection matrices with current inputs - instead of the physical point to physical point projection matrix provided by Geometry */ - ProjectionMatrixType - GetIndexToIndexProjectionMatrix(const unsigned int iProj); - - ProjectionMatrixType - GetVolumeIndexToProjectionPhysicalPointMatrix(const unsigned int iProj); - - itk::Matrix - GetProjectionPhysicalPointToProjectionIndexMatrix(const unsigned int iProj); - /** RTK geometry object */ GeometryPointer m_Geometry; PILineImagePointer m_PILines; diff --git a/include/rtkKatsevichBackProjectionImageFilter.hxx b/include/rtkKatsevichBackProjectionImageFilter.hxx index 361b68147..67a9dbf12 100644 --- a/include/rtkKatsevichBackProjectionImageFilter.hxx +++ b/include/rtkKatsevichBackProjectionImageFilter.hxx @@ -43,23 +43,6 @@ KatsevichBackProjectionImageFilter::KatsevichBackProj this->SetInPlace(true); } - -// template -// void -// KatsevichBackProjectionImageFilter::SetPILines(PILineImageType * image) -//{ -// // Process object is not const-correct so the const_cast is required here -// this->ProcessObject::SetNthInput(1, const_cast(image)); -//} -// -// template -// const typename KatsevichBackProjectionImageFilter::PILineImageType * -// KatsevichBackProjectionImageFilter::GetPILines() const -//{ -// return itkDynamicCastInDebugMode(m_PILines); -//} - - template void KatsevichBackProjectionImageFilter::GenerateInputRequestedRegion() @@ -110,7 +93,7 @@ KatsevichBackProjectionImageFilter::GenerateInputRequ for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) { // Extract the current slice - ProjectionMatrixType matrix = GetIndexToIndexProjectionMatrix(iProj); + ProjectionMatrixType matrix = this->GetIndexToIndexProjectionMatrix(iProj); // Check which part of the projection image will be backprojected in the // volume. @@ -266,7 +249,7 @@ KatsevichBackProjectionImageFilter::DynamicThreadedGe for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) { // Extract the current slice - ProjectionImagePointer projection = GetProjection(iProj); + ProjectionImagePointer projection = this->template GetProjection(iProj); ProjectionMatrixType matrix = this->GetIndexToIndexProjectionMatrix(iProj); interpolator->SetInputImage(projection); @@ -276,9 +259,9 @@ KatsevichBackProjectionImageFilter::DynamicThreadedGe // Cylindrical detector centered on source case if (m_Geometry->GetRadiusCylindricalDetector() != 0) { - ProjectionMatrixType volIndexToProjPP = GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); + ProjectionMatrixType volIndexToProjPP = this->GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); itk::Matrix projPPToProjIndex = - GetProjectionPhysicalPointToProjectionIndexMatrix(iProj); + this->GetProjectionPhysicalPointToProjectionIndexMatrix(iProj); CylindricalDetectorCenteredOnSourceBackprojection( outputRegionForThread, volIndexToProjPP, projPPToProjIndex, projection, lambda, delta_lambda); continue; @@ -330,19 +313,19 @@ KatsevichBackProjectionImageFilter::DynamicThreadedGe double d_out = (st - lambda) / delta_lambda; // Apply PiLineConditions - //if (lambda < sb - delta_lambda) + // if (lambda < sb - delta_lambda) // rho = 0; - //else if (lambda < sb) + // else if (lambda < sb) // rho = 0.5 * (d_in + 1) * (d_in + 1); - //else if (lambda < sb + delta_lambda) + // else if (lambda < sb + delta_lambda) // rho = 0.5 + d_in + 0.5 * d_in * d_in; - //else if (lambda < st - delta_lambda) + // else if (lambda < st - delta_lambda) // rho = 1; - //else if (lambda < st) + // else if (lambda < st) // rho = 0.5 + d_out + 0.5 * d_out * d_out; - //else if (lambda < st + delta_lambda) + // else if (lambda < st + delta_lambda) // rho = 0.5 * (1 + d_out * (1 + d_out)); - //else + // else // rho = 0; if (lambda < sb) @@ -447,19 +430,19 @@ KatsevichBackProjectionImageFilter::CylindricalDetect double d_out = (st - lambda) / delta_lambda; // Apply PiLineConditions - //if (lambda < sb - delta_lambda) + // if (lambda < sb - delta_lambda) // rho = 0; - //else if (lambda < sb) + // else if (lambda < sb) // rho = 0.5 * (d_in + 1) * (d_in + 1); - //else if (lambda < sb + delta_lambda) + // else if (lambda < sb + delta_lambda) // rho = 0.5 + d_in + 0.5 * d_in * d_in; - //else if (lambda < st - delta_lambda) + // else if (lambda < st - delta_lambda) // rho = 1; - //else if (lambda < st) + // else if (lambda < st) // rho = 0.5 + d_out + 0.5 * d_out * d_out; - //else if (lambda < st + delta_lambda) + // else if (lambda < st + delta_lambda) // rho = 0.5 * (1 + d_out * (1 + d_out)); - //else + // else // rho = 0; if (lambda < sb) rho = 0.; @@ -558,19 +541,19 @@ KatsevichBackProjectionImageFilter::OptimizedBackproj double d_in = (lambda - sb) / delta_lambda; double d_out = (st - lambda) / delta_lambda; // Apply PiLineConditions - //if (lambda < sb - delta_lambda) + // if (lambda < sb - delta_lambda) // rho = 0; - //else if (lambda < sb) + // else if (lambda < sb) // rho = 0.5 * (d_in + 1) * (d_in + 1); - //else if (lambda < sb + delta_lambda) + // else if (lambda < sb + delta_lambda) // rho = 0.5 + d_in + 0.5 * d_in * d_in; - //else if (lambda < st - delta_lambda) + // else if (lambda < st - delta_lambda) // rho = 1; - //else if (lambda < st) + // else if (lambda < st) // rho = 0.5 + d_out + 0.5 * d_out * d_out; - //else if (lambda < st + delta_lambda) + // else if (lambda < st + delta_lambda) // rho = 0.5 * (1 + d_out * (1 + d_out)); - //else + // else // rho = 0; if (lambda < sb) @@ -677,19 +660,19 @@ KatsevichBackProjectionImageFilter::OptimizedBackproj double d_out = (st - lambda) / delta_lambda; // Apply PiLineConditions - //if (lambda < sb - delta_lambda) + // if (lambda < sb - delta_lambda) // rho = 0; - //else if (lambda < sb) + // else if (lambda < sb) // rho = 0.5 * (d_in + 1) * (d_in + 1); - //else if (lambda < sb + delta_lambda) + // else if (lambda < sb + delta_lambda) // rho = 0.5 + d_in + 0.5 * d_in * d_in; - //else if (lambda < st - delta_lambda) + // else if (lambda < st - delta_lambda) // rho = 1; - //else if (lambda < st) + // else if (lambda < st) // rho = 0.5 + d_out + 0.5 * d_out * d_out; - //else if (lambda < st + delta_lambda) + // else if (lambda < st + delta_lambda) // rho = 0.5 * (1 + d_out * (1 + d_out)); - //else + // else // rho = 0; if (lambda < sb) @@ -715,125 +698,6 @@ KatsevichBackProjectionImageFilter::OptimizedBackproj } // k } -template -template -typename TProjectionImage::Pointer -KatsevichBackProjectionImageFilter::GetProjection(const unsigned int iProj) -{ - - typename Superclass::InputImagePointer stack = const_cast(this->GetInput(1)); - - const int iProjBuff = stack->GetBufferedRegion().GetIndex(ProjectionImageType::ImageDimension); - - typename TProjectionImage::Pointer projection = TProjectionImage::New(); - typename TProjectionImage::RegionType region; - typename TProjectionImage::SpacingType spacing; - typename TProjectionImage::PointType origin; - - for (unsigned int i = 0; i < TProjectionImage::ImageDimension; i++) - { - origin[i] = stack->GetOrigin()[i]; - spacing[i] = stack->GetSpacing()[i]; - region.SetSize(i, stack->GetBufferedRegion().GetSize()[i]); - region.SetIndex(i, stack->GetBufferedRegion().GetIndex()[i]); - } - if (this->GetTranspose()) - { - typename TProjectionImage::SizeType size = region.GetSize(); - typename TProjectionImage::IndexType index = region.GetIndex(); - std::swap(size[0], size[1]); - std::swap(index[0], index[1]); - std::swap(origin[0], origin[1]); - std::swap(spacing[0], spacing[1]); - region.SetSize(size); - region.SetIndex(index); - } - projection->SetSpacing(spacing); - projection->SetOrigin(origin); - projection->SetRegions(region); - projection->Allocate(); - - const unsigned int npixels = projection->GetBufferedRegion().GetNumberOfPixels(); - const InternalInputPixelType * pi = stack->GetBufferPointer() + (iProj - iProjBuff) * npixels; - InternalInputPixelType * po = projection->GetBufferPointer(); - - // Transpose projection for optimization - if (this->GetTranspose()) - { - for (unsigned int j = 0; j < region.GetSize(0); j++, po -= npixels - 1) - for (unsigned int i = 0; i < region.GetSize(1); i++, po += region.GetSize(0)) - *po = *pi++; - } - else - for (unsigned int i = 0; i < npixels; i++) - *po++ = *pi++; - - return projection; -} - -template -typename KatsevichBackProjectionImageFilter::ProjectionMatrixType -KatsevichBackProjectionImageFilter::GetIndexToIndexProjectionMatrix(const unsigned int iProj) -{ - const unsigned int Dimension = TInputImage::ImageDimension; - - ProjectionMatrixType VolumeIndexToProjectionPhysicalPointMatrix = - GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); - - itk::Matrix ProjectionPhysicalPointToProjectionIndexMatrix = - GetProjectionPhysicalPointToProjectionIndexMatrix(iProj); - - return ProjectionMatrixType(ProjectionPhysicalPointToProjectionIndexMatrix.GetVnlMatrix() * - VolumeIndexToProjectionPhysicalPointMatrix.GetVnlMatrix()); -} - -template -typename KatsevichBackProjectionImageFilter::ProjectionMatrixType -KatsevichBackProjectionImageFilter::GetVolumeIndexToProjectionPhysicalPointMatrix( - const unsigned int iProj) -{ - const unsigned int Dimension = TInputImage::ImageDimension; - - itk::Matrix matrixVol = - GetIndexToPhysicalPointMatrix(this->GetOutput()); - - return ProjectionMatrixType(this->m_Geometry->GetMagnificationMatrices()[iProj].GetVnlMatrix() * - this->m_Geometry->GetSourceTranslationMatrices()[iProj].GetVnlMatrix() * - this->m_Geometry->GetRotationMatrices()[iProj].GetVnlMatrix() * matrixVol.GetVnlMatrix()); -} - -template -itk::Matrix -KatsevichBackProjectionImageFilter::GetProjectionPhysicalPointToProjectionIndexMatrix( - const unsigned int iProj) -{ - const unsigned int Dimension = TInputImage::ImageDimension; - - itk::Matrix matrixStackProj = - GetPhysicalPointToIndexMatrix(this->GetInput(1)); - - itk::Matrix matrixProj; - matrixProj.SetIdentity(); - for (unsigned int i = 0; i < Dimension - 1; i++) - { - matrixProj[i][Dimension - 1] = matrixStackProj[i][Dimension]; - for (unsigned int j = 0; j < Dimension - 1; j++) - matrixProj[i][j] = matrixStackProj[i][j]; - } - - // Transpose projection for optimization - itk::Matrix matrixFlip; - matrixFlip.SetIdentity(); - if (this->GetTranspose()) - { - std::swap(matrixFlip[0][0], matrixFlip[0][1]); - std::swap(matrixFlip[1][0], matrixFlip[1][1]); - } - - return itk::Matrix( - matrixFlip.GetVnlMatrix() * matrixProj.GetVnlMatrix() * - this->m_Geometry->GetProjectionTranslationMatrices()[iProj].GetVnlMatrix()); -} } // end namespace rtk