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..bec72a3f5 --- /dev/null +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.h @@ -0,0 +1,120 @@ +/*========================================================================= + * + * 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. + * + * * \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 + * + * \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..bf538ddd8 --- /dev/null +++ b/include/rtkHilbertTransformOnKappaLinesImageFilter.hxx @@ -0,0 +1,71 @@ +/*========================================================================= + * + * 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() +{ + // Call the superclass' implementation of this method + Superclass::GenerateOutputInformation(); +} + +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..cabe7ca3a --- /dev/null +++ b/include/rtkKatsevichBackProjectionImageFilter.h @@ -0,0 +1,168 @@ +/*========================================================================= + * + * 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 +#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 rtk::BackProjectionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(KatsevichBackProjectionImageFilter); + + /** Standard class type alias. */ + using Self = KatsevichBackProjectionImageFilter; + using Superclass = rtk::BackProjectionImageFilter; + 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 = typename Superclass::ProjectionImageType; + 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 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 + {} + + /** 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..67a9dbf12 --- /dev/null +++ b/include/rtkKatsevichBackProjectionImageFilter.hxx @@ -0,0 +1,704 @@ +/*========================================================================= + * + * 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::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 = this->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 = this->template 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 = this->GetVolumeIndexToProjectionPhysicalPointMatrix(iProj); + itk::Matrix projPPToProjIndex = + this->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 +} + + +} // end namespace rtk + +#endif diff --git a/include/rtkKatsevichBackwardBinningImageFilter.h b/include/rtkKatsevichBackwardBinningImageFilter.h new file mode 100644 index 000000000..87f682f9e --- /dev/null +++ b/include/rtkKatsevichBackwardBinningImageFilter.h @@ -0,0 +1,136 @@ +/*========================================================================= + * + * 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 + * + * Resample the Hilbert transformed data back to the original v-coordinate, according + * to a specified geometry. + * + * \test + * + * \author Jerome Lesaint + * + * \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; + + 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..64b725e3c --- /dev/null +++ b/include/rtkKatsevichBackwardBinningImageFilter.hxx @@ -0,0 +1,233 @@ +/*========================================================================= + * + * 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::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 / D); // Consistent with Forward 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..b9e008a56 --- /dev/null +++ b/include/rtkKatsevichDerivativeImageFilter.h @@ -0,0 +1,130 @@ +/*========================================================================= + * + * 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 formula Eq. 46 (curved det) and Eq. 87 (flat panel) + * from 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..aee37704f --- /dev/null +++ b/include/rtkKatsevichForwardBinningImageFilter.h @@ -0,0 +1,131 @@ +/*========================================================================= + * + * 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 Rebinning along kappa-lines + * + * Rebin the projection data along kappa-lines using linear + * interpolation according to a specified helical geometry. + * + * \author Jerome Lesaint + * + * \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