Skip to content

Commit 19970ed

Browse files
Jerome Lesaintlesaintjerome
Jerome Lesaint
authored and
lesaintjerome
committed
WIP: Add Katsevich reconstruction from helical projections
First attempt for an implementation of Katsevich algorithm as described in the paper "Exact helical reconstruction using native cone-beam geometries", Noo et al., PMB, 2003. Both flat panel and curved detector geometries are implemented though the curved geometry is still under debugging. As for the flat panel implementation, the whole thing more or less works. Some strange "tiling" artifacts remain. The algorithm runs in 4 steps : Compute a derivative of the projections Cosine-weight Forward rebinning + Hilbert transform (from ACoussat PR) + backward rebinning Backprojection. First commit - 2 First commit - nth Add post-cosine weight for curved detector filtering Various debug stuff and computation of FOV radius Misc Update Hilbert transform with ACoussant PR + minor changes requiring future cleanup Clean code from debugging couts Add a Katsevich Reconstruction filter and corresponding application
1 parent 0e39558 commit 19970ed

53 files changed

Lines changed: 5664 additions & 2 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎.gitattributes‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
.git* export-ignore
22
# Custom attribute to mark sources as using our C++/C code style.
3-
[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat,uncrustify
3+
[attr]our-c-style whitespace=tab-in-indent,no-lf-at-eof hooks.style=KWStyle,clangformat
44

55
cmake/*.cxx our-c-style
66
include/*.h* our-c-style

‎.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,3 +33,6 @@ CMakeLists.txt.user*
3333

3434
# Mac System File
3535
.DS_Store
36+
37+
38+
.project

‎applications/CMakeLists.txt‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -102,6 +102,17 @@ add_subdirectory(rtkxradgeometry)
102102
add_subdirectory(rtkimagxgeometry)
103103
add_subdirectory(rtkorageometry)
104104
add_subdirectory(rtkbioscangeometry)
105+
106+
#All executables below are specific to the tomo reconstruction from projections aquired along a helical trajectory of the gantry
107+
add_subdirectory(rtkhelicalgeometry)
108+
add_subdirectory(rtkkatsevichderivative)
109+
add_subdirectory(rtkkatsevichforwardbinning)
110+
add_subdirectory(rtkhilbertonkappalines)
111+
add_subdirectory(rtkpilines)
112+
add_subdirectory(rtkkatsevich)
113+
add_subdirectory(rtkkatsevichtrash)
114+
add_subdirectory(rtkkatsevichrecons)
115+
105116
#=========================================================
106117

107118
#-----------------------------------------------------------------------------

‎applications/rtkhelicalgeometry/CMakeLists.txt‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
WRAP_GGO(rtkhelicalgeometry_GGO_C rtkhelicalgeometry.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo)
2+
add_executable(rtkhelicalgeometry rtkhelicalgeometry.cxx ${rtkhelicalgeometry_GGO_C})
3+
target_link_libraries(rtkhelicalgeometry RTK)
4+
5+
# Installation code
6+
if(NOT RTK_INSTALL_NO_EXECUTABLES)
7+
foreach(EXE_NAME rtkhelicalgeometry)
8+
install(TARGETS ${EXE_NAME}
9+
RUNTIME DESTINATION ${RTK_INSTALL_RUNTIME_DIR} COMPONENT Runtime
10+
LIBRARY DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
11+
ARCHIVE DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} COMPONENT Development)
12+
endforeach()
13+
14+
# Install Python application
15+
install(FILES rtkhelicalgeometry.py
16+
DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT PythonWheelRuntimeLibraries)
17+
endif()
18+

‎applications/rtkhelicalgeometry/rtkhelicalgeometry.cxx‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
/*=========================================================================
2+
*
3+
* Copyright RTK Consortium
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* http://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
#include "rtkhelicalgeometry_ggo.h"
20+
#include "rtkGgoFunctions.h"
21+
22+
#include "rtkThreeDHelicalProjectionGeometryXMLFileWriter.h"
23+
24+
int
25+
main(int argc, char * argv[])
26+
{
27+
GGO(rtkhelicalgeometry, args_info);
28+
29+
// RTK geometry object
30+
using GeometryType = rtk::ThreeDHelicalProjectionGeometry;
31+
GeometryType::Pointer geometry = GeometryType::New();
32+
33+
// Projection matrices
34+
for (int noProj = 0; noProj < args_info.nproj_arg; noProj++)
35+
{
36+
// Compute the angles
37+
double angular_gap = args_info.arc_arg / args_info.nproj_arg;
38+
double first_angle = 0.;
39+
if (!args_info.first_angle_given)
40+
{
41+
first_angle = -0.5 * angular_gap * (args_info.nproj_arg - 1);
42+
}
43+
else
44+
first_angle = args_info.first_angle_arg;
45+
46+
double angle = first_angle + noProj * angular_gap;
47+
48+
49+
// Compute the vertical displacement
50+
double vertical_coverage = args_info.arc_arg / 360 * args_info.pitch_arg;
51+
double vertical_gap = vertical_coverage / args_info.nproj_arg;
52+
double first_sy = 0.;
53+
if (!args_info.first_sy_given)
54+
{
55+
first_sy = -0.5 * vertical_gap * (args_info.nproj_arg - 1);
56+
}
57+
else
58+
{
59+
first_sy = args_info.first_sy_arg;
60+
}
61+
62+
double sy = first_sy + noProj * vertical_gap;
63+
64+
geometry->AddProjection(args_info.sid_arg, args_info.sdd_arg, angle, 0, sy, 0, 0, 0, sy);
65+
}
66+
67+
// Set cylindrical detector radius
68+
if (args_info.rad_cyl_given)
69+
geometry->SetRadiusCylindricalDetector(args_info.rad_cyl_arg);
70+
71+
72+
// Write
73+
rtk::ThreeDHelicalProjectionGeometryXMLFileWriter::Pointer xmlWriter =
74+
rtk::ThreeDHelicalProjectionGeometryXMLFileWriter::New();
75+
xmlWriter->SetFilename(args_info.output_arg);
76+
xmlWriter->SetObject(&(*geometry));
77+
TRY_AND_EXIT_ON_ITK_EXCEPTION(xmlWriter->WriteFile())
78+
79+
return EXIT_SUCCESS;
80+
}

‎applications/rtkhelicalgeometry/rtkhelicalgeometry.ggo‎

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

‎applications/rtkhelicalgeometry/rtkhelicalgeometry.py‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#!/usr/bin/env python
2+
import argparse
3+
import sys
4+
from itk import RTK as rtk
5+
6+
if __name__ == '__main__':
7+
# Argument parsing
8+
parser = argparse.ArgumentParser(description=
9+
"Creates an RTK helical geometry file from regular helical trajectory. See http://www.openrtk.org/Doxygen/DocGeo3D.html for more information.")
10+
11+
12+
parser.add_argument('--nproj', '-n', type=int, help='Number of projections')
13+
parser.add_argument('--output', '-o', help='Output file name')
14+
parser.add_argument('--verbose', '-v', type=bool, default=False, help='Verbose execution')
15+
parser.add_argument('--config', '-c', help='Config file')
16+
parser.add_argument('--first_angle', '-f', type=float, help='First angle in degrees')
17+
parser.add_argument('--first_sy', '-y', type=float, help='First vertical position (default = centered around 0)')
18+
parser.add_argument('--arc', '-a', type=float, default=360, help='Angular arc covevered by the acquisition in degrees')
19+
parser.add_argument('--pitch', '-p', type=float, default=200, help='Helix pitch (vertical displacement in one full (2pi) rotation')
20+
parser.add_argument('--sdd', type=float, default=1536, help='Source to detector distance (mm)')
21+
parser.add_argument('--sid', type=float, default=1000, help='Source to isocenter distance (mm)')
22+
parser.add_argument('--rad_cyl', type=float, default=0, help='Radius cylinder of cylindrical detector')
23+
24+
args = parser.parse_args()
25+
26+
if args.nproj is None or args.output is None :
27+
parser.print_help()
28+
sys.exit()
29+
30+
# Simulated Geometry
31+
GeometryType = rtk.ThreeDCircularProjectionGeometry
32+
geometry = GeometryType.New()
33+
34+
for noProj in range(0, args.nproj):
35+
36+
# Compute the angles
37+
angular_gap = args.arc/args.nproj
38+
if args.first_angle is None :
39+
first_angle = -0.5*angular_gap*(args.nproj-1)
40+
else :
41+
first_angle = args.first_angle
42+
43+
angle = first_angle + noProj * angular_gap
44+
45+
# Compute vertical positions
46+
vertical_coverage = args.arc/360.0*args.pitch
47+
vertical_gap = vertical_coverage/args.nproj
48+
if args.first_sy is None :
49+
first_sy = -0.5*vertical_gap*(args.nproj-1)
50+
else :
51+
first_sy = args.first_sy
52+
53+
sy = first_sy + noProj * vertical_gap
54+
55+
geometry.AddProjection(args.sid,
56+
args.sdd,
57+
angle,
58+
0.,
59+
sy,
60+
0.,
61+
0.,
62+
0.,
63+
sy)
64+
65+
geometry.SetRadiusCylindricalDetector(args.rad_cyl)
66+
67+
writer = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()
68+
writer.SetFilename(args.output)
69+
writer.SetObject(geometry)
70+
writer.WriteFile()

‎applications/rtkhilbertonkappalines/CMakeLists.txt‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
WRAP_GGO(rtkhilbertonkappalines_GGO_C rtkhilbertonkappalines.ggo ../rtkinputprojections_section.ggo ${RTK_BINARY_DIR}/rtkVersion.ggo)
2+
add_executable(rtkhilbertonkappalines rtkhilbertonkappalines.cxx ${rtkhilbertonkappalines_GGO_C})
3+
target_link_libraries(rtkhilbertonkappalines RTK)
4+
5+
# Installation code
6+
if(NOT RTK_INSTALL_NO_EXECUTABLES)
7+
foreach(EXE_NAME rtkhilbertonkappalines)
8+
install(TARGETS ${EXE_NAME}
9+
RUNTIME DESTINATION ${RTK_INSTALL_RUNTIME_DIR} COMPONENT Runtime
10+
LIBRARY DESTINATION ${RTK_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
11+
ARCHIVE DESTINATION ${RTK_INSTALL_ARCHIVE_DIR} COMPONENT Development)
12+
endforeach()
13+
endif()
14+

‎applications/rtkhilbertonkappalines/rtkhilbertonkappalines.cxx‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
/*=========================================================================
2+
*
3+
* Copyright RTK Consortium
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* http://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
#include "rtkhilbertonkappalines_ggo.h"
20+
#include "rtkGgoFunctions.h"
21+
#include "rtkConfiguration.h"
22+
23+
#include "rtkThreeDHelicalProjectionGeometryXMLFileReader.h"
24+
#include "rtkHilbertTransformOnKappaLinesImageFilter.h"
25+
#include "rtkProgressCommands.h"
26+
27+
#include <itkStreamingImageFilter.h>
28+
#include <itkImageRegionSplitterDirection.h>
29+
#include <itkImageFileWriter.h>
30+
31+
int
32+
main(int argc, char * argv[])
33+
{
34+
GGO(rtkhilbertonkappalines, args_info);
35+
36+
using OutputPixelType = float;
37+
constexpr unsigned int Dimension = 3;
38+
39+
using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
40+
#ifdef RTK_USE_CUDA
41+
using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
42+
#else
43+
using OutputImageType = CPUOutputImageType;
44+
#endif
45+
46+
// Projections reader
47+
using ReaderType = rtk::ProjectionsReader<OutputImageType>;
48+
ReaderType::Pointer reader = ReaderType::New();
49+
rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkhilbertonkappalines>(reader, args_info);
50+
51+
if (args_info.verbose_flag)
52+
std::cout << "Reading... " << std::endl;
53+
TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update())
54+
55+
// Geometry
56+
if (args_info.verbose_flag)
57+
std::cout << "Reading geometry information from " << args_info.geometry_arg << "..." << std::endl;
58+
rtk::ThreeDHelicalProjectionGeometryXMLFileReader::Pointer geometryReader;
59+
geometryReader = rtk::ThreeDHelicalProjectionGeometryXMLFileReader::New();
60+
geometryReader->SetFilename(args_info.geometry_arg);
61+
TRY_AND_EXIT_ON_ITK_EXCEPTION(geometryReader->GenerateOutputInformation())
62+
rtk::ThreeDHelicalProjectionGeometry::Pointer geometry;
63+
geometry = geometryReader->GetOutputObject();
64+
geometry->VerifyHelixParameters();
65+
66+
// Check on hardware parameter
67+
#ifndef RTK_USE_CUDA
68+
if (!strcmp(args_info.hardware_arg, "cuda"))
69+
{
70+
std::cerr << "The program has not been compiled with cuda option" << std::endl;
71+
return EXIT_FAILURE;
72+
}
73+
#endif
74+
75+
using HilbertTransformOnKappaLinesType =
76+
rtk::HilbertTransformOnKappaLinesImageFilter<OutputImageType, OutputImageType>;
77+
HilbertTransformOnKappaLinesType::Pointer fwd = HilbertTransformOnKappaLinesType::New();
78+
fwd->SetGeometry(geometry);
79+
fwd->SetInput(reader->GetOutput());
80+
81+
// Write
82+
using WriterType = itk::ImageFileWriter<CPUOutputImageType>;
83+
WriterType::Pointer writer = WriterType::New();
84+
writer->SetFileName(args_info.output_arg);
85+
writer->SetInput(fwd->GetOutput());
86+
87+
if (args_info.verbose_flag)
88+
std::cout << "Reconstructing and writing... " << std::endl;
89+
90+
TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
91+
92+
return EXIT_SUCCESS;
93+
}

‎applications/rtkhilbertonkappalines/rtkhilbertonkappalines.ggo‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
package "rtkhilbertonkappalines"
2+
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."
3+
4+
option "verbose" v "Verbose execution" flag off
5+
option "config" - "Config file" string no
6+
option "geometry" g "XML geometry file name" string yes
7+
option "output" o "Output file name" string yes
8+
option "hardware" - "Hardware used for computation" values="cpu","cuda" no default="cpu"
9+

0 commit comments

Comments
 (0)