From 4009b8605b9aa8db187390cf89789c000a9f0e0c Mon Sep 17 00:00:00 2001 From: axel garcia Date: Tue, 17 Dec 2024 11:30:20 +0100 Subject: [PATCH] DOC : Update new ReadTheDocs with content from the old documentation --- CodeContribution.md | 25 ++ INSTALLATION.md | 20 ++ .../rtkThreeDCircularProjectionGeometry.dox | 289 ------------------ documentation/docs/CMakeLists.txt | 1 + documentation/docs/Geometry.md | 217 +++++++++++++ documentation/docs/Tutorial.md | 22 ++ examples/4DRooster/Blurred.jpg.sha512 | 1 + examples/4DRooster/MotionMask.jpg.sha512 | 1 + examples/4DRooster/README.md | 225 ++++++++++++++ examples/4DRooster/Signal.jpg.sha512 | 1 + examples/AmsterdamShroud/Amsterdam.png.sha512 | 1 + .../Moving-Phantom-Sinogram.png.sha512 | 1 + examples/AmsterdamShroud/README.md | 13 + examples/ConjugateGradient/Code3D.sh | 21 ++ .../ConjugateGradient-Sinogram.png.sha512 | 1 + .../ConjugateGradient.png.sha512 | 1 + examples/ConjugateGradient/README.md | 42 +++ .../GammexPhantom.png.sha512 | 1 + examples/CreateGammexPhantom/README.md | 10 + .../Code.sh | 12 + .../Overlay.png.sha512 | 1 + .../README.md | 9 + .../Sinogram.png.sha512 | 1 + .../ElektaReconstruction/Elekta.png.sha512 | 1 + examples/ElektaReconstruction/README.md | 72 +++++ examples/ForwardProjection/Code.sh | 8 + .../POPI-Reconstruction.png.sha512 | 1 + .../POPI-Sinogram.png.sha512 | 1 + examples/ForwardProjection/README.md | 13 + .../Blurred.jpg.sha512 | 1 + .../Blurred_vs_mc.gif.sha512 | 1 + .../MotionMask.jpg.sha512 | 1 + .../README.md | 169 ++++++++++ .../Signal.jpg.sha512 | 1 + .../VectorField.gif.sha512 | 1 + examples/RayBoxIntersection/Code.sh | 11 + examples/RayBoxIntersection/README.md | 9 + .../RayBox-Sinogram.png.sha512 | 1 + examples/RayBoxIntersection/RayBox.png.sha512 | 1 + .../Code.sh | 12 + .../Overlay.png.sha512 | 1 + .../README.md | 9 + .../Sinogram.png.sha512 | 1 + examples/VarianReconstruction/README.md | 116 +++++++ .../VarianReconstruction/Varian.png.sha512 | 1 + .../VarianProBeam.png.sha512 | 1 + examples/WaterPreCorrection/README.md | 11 + examples/index.md | 12 + index.md | 11 +- .../SheppLogan_forbild.txt.md5 | 1 + 50 files changed, 1091 insertions(+), 293 deletions(-) create mode 100644 CodeContribution.md delete mode 100644 documentation/Doxygen/rtkThreeDCircularProjectionGeometry.dox create mode 100644 documentation/docs/Geometry.md create mode 100644 documentation/docs/Tutorial.md create mode 100644 examples/4DRooster/Blurred.jpg.sha512 create mode 100644 examples/4DRooster/MotionMask.jpg.sha512 create mode 100644 examples/4DRooster/README.md create mode 100644 examples/4DRooster/Signal.jpg.sha512 create mode 100644 examples/AmsterdamShroud/Amsterdam.png.sha512 create mode 100644 examples/AmsterdamShroud/Moving-Phantom-Sinogram.png.sha512 create mode 100644 examples/AmsterdamShroud/README.md create mode 100755 examples/ConjugateGradient/Code3D.sh create mode 100644 examples/ConjugateGradient/ConjugateGradient-Sinogram.png.sha512 create mode 100644 examples/ConjugateGradient/ConjugateGradient.png.sha512 create mode 100644 examples/ConjugateGradient/README.md create mode 100644 examples/CreateGammexPhantom/GammexPhantom.png.sha512 create mode 100644 examples/CreateGammexPhantom/README.md create mode 100644 examples/DaubechiesWaveletsRegularizedReconstruction/Code.sh create mode 100644 examples/DaubechiesWaveletsRegularizedReconstruction/Overlay.png.sha512 create mode 100644 examples/DaubechiesWaveletsRegularizedReconstruction/README.md create mode 100644 examples/DaubechiesWaveletsRegularizedReconstruction/Sinogram.png.sha512 create mode 100644 examples/ElektaReconstruction/Elekta.png.sha512 create mode 100644 examples/ElektaReconstruction/README.md create mode 100644 examples/ForwardProjection/Code.sh create mode 100644 examples/ForwardProjection/POPI-Reconstruction.png.sha512 create mode 100644 examples/ForwardProjection/POPI-Sinogram.png.sha512 create mode 100644 examples/ForwardProjection/README.md create mode 100644 examples/MotionCompensationReconstruction/Blurred.jpg.sha512 create mode 100644 examples/MotionCompensationReconstruction/Blurred_vs_mc.gif.sha512 create mode 100644 examples/MotionCompensationReconstruction/MotionMask.jpg.sha512 create mode 100644 examples/MotionCompensationReconstruction/README.md create mode 100644 examples/MotionCompensationReconstruction/Signal.jpg.sha512 create mode 100644 examples/MotionCompensationReconstruction/VectorField.gif.sha512 create mode 100644 examples/RayBoxIntersection/Code.sh create mode 100644 examples/RayBoxIntersection/README.md create mode 100644 examples/RayBoxIntersection/RayBox-Sinogram.png.sha512 create mode 100644 examples/RayBoxIntersection/RayBox.png.sha512 create mode 100644 examples/TotalVariationRegularizedReconstruction/Code.sh create mode 100644 examples/TotalVariationRegularizedReconstruction/Overlay.png.sha512 create mode 100644 examples/TotalVariationRegularizedReconstruction/README.md create mode 100644 examples/TotalVariationRegularizedReconstruction/Sinogram.png.sha512 create mode 100644 examples/VarianReconstruction/README.md create mode 100644 examples/VarianReconstruction/Varian.png.sha512 create mode 100644 examples/VarianReconstruction/VarianProBeam.png.sha512 create mode 100644 examples/WaterPreCorrection/README.md create mode 100644 test/Input/GeometricPhantom/SheppLogan_forbild.txt.md5 diff --git a/CodeContribution.md b/CodeContribution.md new file mode 100644 index 000000000..8d079510c --- /dev/null +++ b/CodeContribution.md @@ -0,0 +1,25 @@ +# Code contribution + +## Coding style + +RTK is based on ITK and aims at following its coding conventions. Any developer should follow these conventions when submitting new code or contributions to the existing one. We strongly recommend you to read thoroughly [ITK's style guide](http://www.itk.org/Wiki/ITK/Coding_Style_Guide). + +## Testing + +This section describes how to add/edit datasets for testing purposes for RTK. Datasets are not stored in the GIT repository for efficiency and also to avoid having large history due to binary files. Instead the files are stored on a [Girder](http://data.kitware.com) instance. Here's the recipe to add new datasets: + +1. Register/Login to Girder hosted at Kitware: [http://data.kitware.com](http://data.kitware.com) +2. Locate the RTK collection: [https://data.kitware.com/#collection/5a7706878d777f0649e04776](https://data.kitware.com/#collection/5a7706878d777f0649e04776) +3. Upload the new datasets in the appropriate folder. If you don't have the necessary privileges please email the mailing list +4. In the GIT repository, add in testing/Data a file with the exact filename of the original file **but with the .md5 extension**. Inside that file put the md5sum of the file on Girder. +5. When adding a test use the new macro 'RTK_ADD_TEST' instead of 'ADD_TEST' and specify the datasets you want CTest to download by appending the data to 'DATA{}'. For example: + +``` +RTK_ADD_TEST(NAME rtkimagxtest + COMMAND ${EXECUTABLE\_OUTPUT\_PATH}/rtkimagxtest + DATA{Data/Input/ImagX/raw.xml,raw.raw} + DATA{Data/Baseline/ImagX/attenuation.mha}) +``` +## Dashboard + +* The RTK dashboard is available at [RTK Dashboard](http://my.cdash.org/index.php?project=RTK) \ No newline at end of file diff --git a/INSTALLATION.md b/INSTALLATION.md index 2d3278819..cdd22608f 100644 --- a/INSTALLATION.md +++ b/INSTALLATION.md @@ -38,6 +38,26 @@ We also provide pre-compiled [CUDA](https://developer.nvidia.com/cuda-toolkit) p python -m pip install itk-rtk-cuda116 ``` +Docker +--------------------- + +Another installation solution is to use the Docker solution provided by Thomas Baudier: + +``` +docker pull tbaudier/rtk:v1.3.0 +docker run -ti --rm -e DISPLAY=$DISPLAY -v \[Documents\]:/home tbaudier/rtk:v1.3.0 bash +``` +Information on what is installed can be reached using the commands: +``` +docker images +docker ps -a +``` +To clean it after use, you can do: +``` +docker rm -f \[container id\] +docker rmi -f \[image id\] +``` + Getting started --------------- See [GettingStarted.md](GettingStarted.md). Your `CMakeLists.txt` can now use RTK when importing ITK as shown in the [FirstReconstruction's CMakeLists.txt](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction/CMakeLists.txt#L7). diff --git a/documentation/Doxygen/rtkThreeDCircularProjectionGeometry.dox b/documentation/Doxygen/rtkThreeDCircularProjectionGeometry.dox deleted file mode 100644 index 2ce39005c..000000000 --- a/documentation/Doxygen/rtkThreeDCircularProjectionGeometry.dox +++ /dev/null @@ -1,289 +0,0 @@ -/*! -\page DocGeo3D RTK 3D circular projection geometry - -\tableofcontents - -\section Purpose - -The purpose of this page is to describe the geometry format used in RTK to -relate a tomography to projection images. There is currently only one geometry -format, ThreeDCircularProjectionGeometry. - -\section Units - - - Degrees are used to store angles in the geometry objects. Angles are - wrapped between 0 and 360 degrees. - - No unit is enforced for distances but it is the responsibility of the - user to have a consistent unit for all distances (pixel and voxel spacings, - geometry parameters...). Millimeters are typically used in ITK and DICOM. - -\section ICS Image coordinate system - -An \p itk::Image<> contains information to convert voxel indices to physical -coordinates using its members \p m_Origin, \p m_Spacing and \p m_Direction. -Voxel coordinates are not used in RTK except for internal computation. The -conversion from voxel index coordinates to physical coordinates and the -dimensions of the images are out of the scope of this document. In the -following, origin refers to point with coordinates \f$\vec{0}\f$ mm (and not to -\p m_Origin in ITK). - -\section FCS Fixed coordinate system - -The fixed coordinate system \f$(x,y,z)\f$ in RTK is the coordinate system of the -tomography with the isocenter at the origin \f$(0,0,0)\f$. - -\section PG ProjectionGeometry - -This is the mother class for relating a TDimension-D tomography to a -(TDimension-1)-D projection image. It holds a vector of -(TDimension)x(TDimension+1) projection matrices accessible via -\p GetMatrices. The construction of those matrices is geometry dependent. - -\section TDCPG ThreeDCircularProjectionGeometry - -This class is meant to define a set of 2D projection images, acquired with a -flat panel along a circular trajectory, around a 3D tomography. The trajectory -does not have to be strictly circular but it is assumed in some reconstruction -algorithms that the rotation axis is y. The description of the geometry is -based on the international standard IEC 61217 which has been designed for -cone-beam imagers on isocentric radiotherapy systems but it can be used for any -3D circular trajectory. The fixed coordinate system of RTK and the fixed -coordinate system of IEC 61217 are the same. - -9 parameters are used per projection to define the position of the source and -the detector relatively to the fixed coordinate system. The 9 parameters are -set with the method \p AddProjection. Default values are provided for the -parameters which are not mandatory. Note that explicit names have been used but -this does not necessarily correspond to the value returned by the scanner which -can use its own parameterization. - -\subsection DO Detector orientation - -\subsubsection IDO Initial detector orientation - -With all parameters set to 0, the detector is normal to the z direction of the -fixed coordinate system, similarly to the x-ray image receptor in the IEC -61217. - -\subsubsection RO Rotation order - -Three rotation angles are used to define the -orientation of the detector. The ZXY convention of Euler angles is used for -detector orientation where GantryAngle is the rotation around y, -OutOfPlaneAngle the rotation around x and InPlaneAngle the rotation around z. -These three angles are detailed in the following. - -\subsubsection GantryAngle - -Gantry angle of the scanner. It corresponds to \f$\phi g\f$ in Section 2.3 of IEC -61217: - -> The rotation of the "g" system is defined by the rotation of coordinate axes -> Xg, Zg by an angle \f$\phi g\f$ about axis Yg (therefore about Yf of the "f" -> system). -> -> An increase in the value of \f$\phi g\f$ corresponds to a clockwise rotation of the -> GANTRY as viewed along the horizontal axis Yf from the ISOCENTRE towards the -> GANTRY. - -\subsubsection OutOfPlaneAngle - -Out of plane rotation of the flat panel complementary to the GantryAngle -rotation, i.e. with a rotation axis perpendicular to the gantry rotation axis -and parallel to the flat panel. It is optional with a default value equals to -0. There is no corresponding rotation in IEC 61217. After gantry rotation, the -rotation is defined by the rotation of the coordinate axes y and z about x. An -increase in the value of OutOfPlaneAngle corresponds to a counter-clockwise -rotation of the flat panel as viewed from a positive value along the x axis -towards the isocenter. - -\subsubsection InPlaneAngle - -In plane rotation of the 2D projection. It is optional with 0 as default value. -If OutOfPlaneAngle equals 0, it corresponds to \f$\theta r\f$ in Section 2.6 of IEC -61217: - - -> The rotation of the "r" system is defined by the rotation of the coordinate -> axes Xr, Yr about Zr (parallel to axis Zg) by an angle \f$\theta r\f$. -> -> An increase in the value of angle \f$\theta r\f$ corresponds to a counter-clockwise -> rotation of the X- RAY IMAGE RECEPTOR as viewed from the RADIATION SOURCE. - -\subsubsection RM Rotation matrix - -The rotation matrix in homogeneous coordinate is then (constructed with -\p itk::Euler3DTransform::ComputeMatrix() with opposite angles -because we rotate the volume coordinates instead of the scanner): - -\f[ - \begin{split} - M_R = & - \begin{pmatrix} - \cos(-InPlaneAngle) & -\sin(-InPlaneAngle) & 0 & 0\\ - \sin(-InPlaneAngle) & \cos(-InPlaneAngle) & 0 & 0\\ - 0 & 0 & 1 & 0\\ - 0 & 0 & 0 & 1 - \end{pmatrix}\\ - &\times - \begin{pmatrix} - 1 & 0 & 0 & 0\\ - 0 & \cos(-OutOfPlaneAngle) & -\sin(-OutOfPlaneAngle) & 0\\ - 0 & \sin(-OutOfPlaneAngle) & \cos(-OutOfPlaneAngle) & 0\\ - 0 & 0 & 0 & 1 - \end{pmatrix}\\ - &\times - \begin{pmatrix} - \cos(-GantryAngle) & 0 & \sin(-GantryAngle) & 0 \\ - 0 & 1 & 0 & 0 \\ - -\sin(-GantryAngle) & 0 & \cos(-GantryAngle) & 0 \\ - 0 & 0 & 0 & 1 - \end{pmatrix} - \end{split} -\f] - -\subsection Drawings - -The following drawing describes the parameters of the source and the detector -positions in the rotated coordinate system \f$(Rx,Ry,Rz)\f$ (i.e., oriented -according to the detector orientation), with its origin at the isocenter, when -all values are positive (but all distances can be negative in this geometry): - -\image html https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry.svg - -These 6 parameters are used to describe any source and detector positions. It -is simpler to understand the circular geometry when all Offset values equal 0 : - -\image html https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_aligned.svg - -\subsection SP Source position - -The source position is defined with respect to the isocenter with three -parameters, SourceOffsetX, SourceOffsetY and SourceToIsocenterDistance. -(SourceOffsetX,SourceOffsetY,SourceToIsocenterDistance) are the coordinates of -the source in the rotated coordinated system. In IEC 61217, -SourceToIsocenterDistance is the RADIATION SOURCE axis distance, SAD. -SourceOffsetX and SourceOffsetY are optional and zero by default. - -\subsection DP Detector position - -The detector position is defined with respect to the source with three -parameters: ProjectionOffsetX, ProjectionOffsetY and SourceToDetectorDistance. -(ProjectionOffsetX,ProjectionOffsetY,SourceToIsocenterDistance-SourceToDetectorDistance) -are the coordinates of the detector origin \f$(0,0)\f$ in the rotated -coordinated system. In IEC 61217, SourceToDetectorDistance is the RADIATION -SOURCE to IMAGE RECEPTION AREA distance, SID. ProjectionOffsetX and -ProjectionOffsetY are optional and zero by default. - -\subsection FM Final matrix - -Each matrix, accessible via \p GetMatrices, is constructed with: - -\f[ - \begin{split} - M_P = - &\begin{pmatrix} - 1 & 0 & SourceOffsetX-ProjectionOffsetX \\ - 0 & 1 & SourceOffsetY-ProjectionOffsetY \\ - 0 & 0 & 1 - \end{pmatrix}\\ - &\times - \begin{pmatrix} - -SourceToDetectorDistance & 0 & 0 & 0 \\ - 0 & -SourceToDetectorDistance & 0 & 0 \\ - 0 & 0 & 1 & -SourceToIsocenterDistance - \end{pmatrix}\\ - &\times - \begin{pmatrix} - 1 & 0 & 0 & -SourceOffsetX \\ - 0 & 1 & 0 & -SourceOffsetY \\ - 0 & 0 & 1 & 0 \\ - 0 & 0 & 0 & 1 - \end{pmatrix}\\ - &\times - M_R - \end{split} -\f] - -\subsection DR Detector radius - -In addition to flat panel detectors, some of the forward and back projectors in -RTK can handle cylindrical detectors. The radius of the cylindrical detector is -stored only once, as the variable RadiusCylindricalDetector. The default value -for RadiusCylindricalDetector is 0, and indicates that the detector is a flat -panel (i.e. infinite radius, but 0 is easier to deal with). When the value is -non-zero, then the flat detector is curved according to the radius and remain -tangent to the corresponding flat detector along the line defined by the -detector origin \f$(0,0)\f$ and second axis of the detector without accouting -for the parameters ProjectionOffsetX and ProjectionOffsetY. The latter two -allow to modify the Origin of each projection as is the case for a flat panel. -The cylindrical detector geometry is illustrated in the following scheme: - -\image html https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_cylindrical.svg - -This scheme is based on the previous one with all offsets equal 0 but this is not required. - -\subsection ParG Parallel geometry - -When SourceToDetectorDistance is set to 0, the geometry is assumed to be -parallel (i.e. infinite distance, but 0 is easier to deal with). The detector -is then flat. The rays are perpendicular to the detector plane which is -oriented similarly to the divergent geometry. The (plane) source is actually -placed at a distance SourceToIsocenterDistance from the isocenter and the -detector is placed symetrically around the origin \f$(0,0,0)\f$ at the same -SourceToIsocenterDistance. This is summarized in the following scheme: - -\image html https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_parallel.svg - -In this case, the projection matrix becomes: - -\f[ - M_P = - \begin{pmatrix} - 1 & 0 & 0 & -ProjectionOffsetX \\ - 0 & 1 & 0 & -ProjectionOffsetY \\ - 0 & 0 & 0 & 1 - \end{pmatrix} - \times M_R. -\f] - -\subsection XF XML file - -ThreeDCircularProjectionGeometry can be saved and loaded from an XML file. If -the parameter is equal to the default value for all projections, it is not -stored in the file. If it is equal for all projections but different from the -default value, it is stored once. Otherwise, it is stored for each projection. -The matrix is given for information. It is read and checked to be consistent -with the parameters but a manual modification of the file must consistently -modify both the parameters and the matrix. An example is given hereafter: - -\verbatim - - - - 1000 - 1536 - 1536 - - 271.847274780273 - -117.056503295898 - -1.01195001602173 - - -166.5093078829 0 -1531.42837748039 -117056.503295898 - -1.01142410874151 -1536 0.0326206557691505 -1011.95001602173 - -0.999480303105996 0 0.0322354417240802 -1000 - - - - 271.852905273438 - -117.056831359863 - -1.01187002658844 - - -166.660129424325 0 -1531.41199650136 -117056.831359863 - -1.01134095059569 -1536 0.0327174625589984 -1011.87002658844 - -0.999477130482326 0 0.0323336611415466 -1000 - - - -\endverbatim -*/ diff --git a/documentation/docs/CMakeLists.txt b/documentation/docs/CMakeLists.txt index 7d5c2e8d5..0dd31b063 100644 --- a/documentation/docs/CMakeLists.txt +++ b/documentation/docs/CMakeLists.txt @@ -12,6 +12,7 @@ if(RTK_BUILD_SPHINX) COMMAND ${CMAKE_COMMAND} -E copy "${RTK_SOURCE_DIR}/index.md" "${RTK_DOC_OUTPUT_DIR}/index.md" COMMAND ${CMAKE_COMMAND} -E copy "${RTK_SOURCE_DIR}/GettingStarted.md" "${RTK_DOC_OUTPUT_DIR}/GettingStarted.md" COMMAND ${CMAKE_COMMAND} -E copy "${RTK_SOURCE_DIR}/INSTALLATION.md" "${RTK_DOC_OUTPUT_DIR}/INSTALLATION.md" + COMMAND ${CMAKE_COMMAND} -E copy "${RTK_SOURCE_DIR}/CodeContribution.md" "${RTK_DOC_OUTPUT_DIR}/CodeContribution.md" COMMENT "Copying documentation sources" ) diff --git a/documentation/docs/Geometry.md b/documentation/docs/Geometry.md new file mode 100644 index 000000000..5e38873d7 --- /dev/null +++ b/documentation/docs/Geometry.md @@ -0,0 +1,217 @@ +# 3D Circular Projection Geometry + +## Purpose + +The purpose of this page is to describe the geometry format used in RTK to relate a tomography to projection images. There is currently only one geometry format, `ThreeDCircularProjectionGeometry`. + +## Units + +- Degrees are used to store angles in the geometry objects. Angles are wrapped between 0 and 360 degrees. +- No unit is enforced for distances but it is the responsibility of the user to have a consistent unit for all distances (pixel and voxel spacings, geometry parameters...). Millimeters are typically used in ITK and DICOM. + +## Image Coordinate System + +An `itk::Image<>` contains information to convert voxel indices to physical coordinates using its members `m_Origin`, `m_Spacing`, and `m_Direction`. Voxel coordinates are not used in RTK except for internal computation. The conversion from voxel index coordinates to physical coordinates and the dimensions of the images are out of the scope of this document. In the following, origin refers to point with coordinates $\vec{0}$ mm (and not to `m_Origin` in ITK). + +## Fixed Coordinate System + +The fixed coordinate system $(x,y,z)$ in RTK is the coordinate system of the tomography with the isocenter at the origin $(0,0,0)$. + +## ProjectionGeometry + +This is the mother class for relating a TDimension-D tomography to a (TDimension-1)-D projection image. It holds a vector of (TDimension)x(TDimension+1) projection matrices accessible via `GetMatrices`. The construction of those matrices is geometry dependent. + +## ThreeDCircularProjectionGeometry + +This class is meant to define a set of 2D projection images, acquired with a flat panel along a circular trajectory, around a 3D tomography. The trajectory does not have to be strictly circular but it is assumed in some reconstruction algorithms that the rotation axis is y. The description of the geometry is based on the international standard IEC 61217 which has been designed for cone-beam imagers on isocentric radiotherapy systems but it can be used for any 3D circular trajectory. The fixed coordinate system of RTK and the fixed coordinate system of IEC 61217 are the same. + +9 parameters are used per projection to define the position of the source and the detector relatively to the fixed coordinate system. The 9 parameters are set with the method `AddProjection`. Default values are provided for the parameters which are not mandatory. Note that explicit names have been used but this does not necessarily correspond to the value returned by the scanner which can use its own parameterization. + +### Detector Orientation + +#### Initial Detector Orientation + +With all parameters set to 0, the detector is normal to the z direction of the fixed coordinate system, similarly to the x-ray image receptor in the IEC 61217. + +#### Rotation Order + +Three rotation angles are used to define the orientation of the detector. The ZXY convention of Euler angles is used for detector orientation where `GantryAngle` is the rotation around y, `OutOfPlaneAngle` the rotation around x, and `InPlaneAngle` the rotation around z. These three angles are detailed in the following. + +#### GantryAngle + +Gantry angle of the scanner. It corresponds to $\phi_g$ in Section 2.3 of IEC 61217: + +> The rotation of the "g" system is defined by the rotation of coordinate axes Xg, Zg by an angle $\phi_g$ about axis Yg (therefore about Yf of the "f" system). +> +> An increase in the value of $\phi_g$ corresponds to a clockwise rotation of the GANTRY as viewed along the horizontal axis Yf from the ISOCENTER towards the GANTRY. + +#### OutOfPlaneAngle + +Out of plane rotation of the flat panel complementary to the GantryAngle rotation, i.e., with a rotation axis perpendicular to the gantry rotation axis and parallel to the flat panel. It is optional with a default value equal to 0. There is no corresponding rotation in IEC 61217. After gantry rotation, the rotation is defined by the rotation of the coordinate axes y and z about x. An increase in the value of `OutOfPlaneAngle` corresponds to a counter-clockwise rotation of the flat panel as viewed from a positive value along the x axis towards the isocenter. + +#### InPlaneAngle + +In plane rotation of the 2D projection. It is optional with 0 as default value. If `OutOfPlaneAngle` equals 0, it corresponds to $\theta_r$ in Section 2.6 of IEC 61217: + +> The rotation of the "r" system is defined by the rotation of the coordinate axes Xr, Yr about Zr (parallel to axis Zg) by an angle $\theta_r$. +> +> An increase in the value of angle $\theta_r$ corresponds to a counter-clockwise rotation of the X-RAY IMAGE RECEPTOR as viewed from the RADIATION SOURCE. + +#### Rotation Matrix + +The rotation matrix in homogeneous coordinates is then (constructed with `itk::Euler3DTransform::ComputeMatrix()` with opposite angles because we rotate the volume coordinates instead of the scanner): + +$$ +\begin{split} +M_R = & +\begin{pmatrix} +\cos(-InPlaneAngle) & -\sin(-InPlaneAngle) & 0 & 0\\ +\sin(-InPlaneAngle) & \cos(-InPlaneAngle) & 0 & 0\\ +0 & 0 & 1 & 0\\ +0 & 0 & 0 & 1 +\end{pmatrix}\\ +&\times +\begin{pmatrix} +1 & 0 & 0 & 0\\ +0 & \cos(-OutOfPlaneAngle) & -\sin(-OutOfPlaneAngle) & 0\\ +0 & \sin(-OutOfPlaneAngle) & \cos(-OutOfPlaneAngle) & 0\\ +0 & 0 & 0 & 1 +\end{pmatrix}\\ +&\times +\begin{pmatrix} +\cos(-GantryAngle) & 0 & \sin(-GantryAngle) & 0 \\ +0 & 1 & 0 & 0 \\ +-\sin(-GantryAngle) & 0 & \cos(-GantryAngle) & 0 \\ +0 & 0 & 0 & 1 +\end{pmatrix} +\end{split} +$$ + +### Drawings + +The following drawing describes the parameters of the source and the detector positions in the rotated coordinate system $(Rx,Ry,Rz)$ (i.e., oriented according to the detector orientation), with its origin at the isocenter, when all values are positive (but all distances can be negative in this geometry): + +![Drawing](https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry.svg) + +These 6 parameters are used to describe any source and detector positions. It is simpler to understand the circular geometry when all Offset values equal 0: + +![Drawing](https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_aligned.svg) + +### Source Position + +The source position is defined with respect to the isocenter with three parameters, `SourceOffsetX`, `SourceOffsetY`, and `SourceToIsocenterDistance`. `(SourceOffsetX, SourceOffsetY, SourceToIsocenterDistance)` are the coordinates of the source in the rotated coordinated system. In IEC 61217, `SourceToIsocenterDistance` is the RADIATION SOURCE axis distance, SAD. `SourceOffsetX` and `SourceOffsetY` are optional and zero by default. + +### Detector Position + +The detector position is defined with respect to the source with three parameters: `ProjectionOffsetX`, `ProjectionOffsetY`, and `SourceToDetectorDistance`. `(ProjectionOffsetX, ProjectionOffsetY, SourceToIsocenterDistance - SourceToDetectorDistance)` are the coordinates of the detector origin $(0,0)$ in the rotated coordinated system. In IEC 61217, `SourceToDetectorDistance` is the RADIATION SOURCE to IMAGE RECEPTION AREA distance, SID. `ProjectionOffsetX` and `ProjectionOffsetY` are optional and zero by default. + +### Final Matrix + +Each matrix, accessible via `GetMatrices`, is constructed with: + +$$ + \begin{split} + M_P = + &\begin{pmatrix} + 1 & 0 & SourceOffsetX-ProjectionOffsetX \\ + 0 & 1 & SourceOffsetY-ProjectionOffsetY \\ + 0 & 0 & 1 + \end{pmatrix}\\ + &\times + \begin{pmatrix} + -SourceToDetectorDistance & 0 & 0 & 0 \\ + 0 & -SourceToDetectorDistance & 0 & 0 \\ + 0 & 0 & 1 & -SourceToIsocenterDistance + \end{pmatrix}\\ + &\times + \begin{pmatrix} + 1 & 0 & 0 & -SourceOffsetX \\ + 0 & 1 & 0 & -SourceOffsetY \\ + 0 & 0 & 1 & 0 \\ + 0 & 0 & 0 & 1 + \end{pmatrix}\\ + &\times + M_R + \end{split} +$$ + +### Detector Radius + +In addition to flat panel detectors, some of the forward and back projectors in +RTK can handle cylindrical detectors. The radius of the cylindrical detector is +stored only once, as the variable `RadiusCylindricalDetector`. The default value +for `RadiusCylindricalDetector` is 0, and indicates that the detector is a flat +panel (i.e. infinite radius, but 0 is easier to deal with). When the value is +non-zero, then the flat detector is curved according to the radius and remain +tangent to the corresponding flat detector along the line defined by the +detector origin `(0,0)` and second axis of the detector without accouting +for the parameters `ProjectionOffsetX` and `ProjectionOffsetY`. The latter two +allow to modify the Origin of each projection as is the case for a flat panel. +The cylindrical detector geometry is illustrated in the following scheme: + +![Drawing](https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_cylindrical.svg) + +This scheme is based on the previous one with all offsets equal 0 but this is not required. + +### Parallel Geometry + +When `SourceToDetectorDistance` is set to 0, the geometry is assumed to be +parallel (i.e. infinite distance, but 0 is easier to deal with). The detector +is then flat. The rays are perpendicular to the detector plane which is +oriented similarly to the divergent geometry. The (plane) source is actually +placed at a distance `SourceToIsocenterDistance` from the isocenter and the +detector is placed symetrically around the origin `(0,0,0)` at the same +`SourceToIsocenterDistance`. This is summarized in the following scheme: + +![Drawing](https://www.openrtk.org/RTK/img/ThreeDCircularProjectionGeometry_parallel.svg) + +In this case, the projection matrix becomes: + +$$ + M_P = + \begin{pmatrix} + 1 & 0 & 0 & -ProjectionOffsetX \\ + 0 & 1 & 0 & -ProjectionOffsetY \\ + 0 & 0 & 0 & 1 + \end{pmatrix} + \times M_R. +$$ + +### XML File + +`ThreeDCircularProjectionGeometry` can be saved and loaded from an XML file. If +the parameter is equal to the default value for all projections, it is not +stored in the file. If it is equal for all projections but different from the +default value, it is stored once. Otherwise, it is stored for each projection. +The matrix is given for information. It is read and checked to be consistent +with the parameters but a manual modification of the file must consistently +modify both the parameters and the matrix. An example is given hereafter: + +```xml + + + + 1000 + 1536 + 1536 + + 271.847274780273 + -117.056503295898 + -1.01195001602173 + + -166.5093078829 0 -1531.42837748039 -117056.503295898 + -1.01142410874151 -1536 0.0326206557691505 -1011.95001602173 + -0.999480303105996 0 0.0322354417240802 -1000 + + + + 271.852905273438 + -117.056831359863 + -1.01187002658844 + + -166.660129424325 0 -1531.41199650136 -117056.831359863 + -1.01134095059569 -1536 0.0327174625589984 -1011.87002658844 + -0.999477130482326 0 0.0323336611415466 -1000 + + + diff --git a/documentation/docs/Tutorial.md b/documentation/docs/Tutorial.md new file mode 100644 index 000000000..c7af1f023 --- /dev/null +++ b/documentation/docs/Tutorial.md @@ -0,0 +1,22 @@ + +# Tutorials +------------------------------------------------------- +## Building a HelloWorld application + +RTK is a library, therefore it's meant to be integrated into application. This tutorial shows how to create a simple FirstReconstruction project that links with RTK. The source code for this tutorial is located in [RTK/examples/FirstReconstruction](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction). + +* First you need to create a [CMakeLists.txt](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction/CMakeLists.txt) + +```{literalinclude} ../../examples/FirstReconstruction/CMakeLists.txt +:language: cmake +``` +* Create a [FirstReconstruction.cxx](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction/FirstReconstruction.cxx) file +```{literalinclude} ../../examples/FirstReconstruction/FirstReconstruction.cxx +``` +* Run CMake on the FirstReconstruction directory and create a HelloWorld-bin, +* Configure and build the project using your favorite compiler, +* Run `FirstReconstruction image.mha geometry.xml`. If everything runs correctly, you should see a few messages ending with `Done!` and two new files in the current directory, image.mha and geometry.xml. image.mha is the reconstruction of a sphere from 360 projections and geometry.xml is the geometry file in the [RTK format](./Geometry.md). + +## Modifying a basic RTK application + +In [applications/rtktutorialapplication/](https://github.com/RTKConsortium/RTK/blob/master/applications/rtktutorialapplication), you will find a very basic RTK application that can be used as a starting point for building more complex applications. \ No newline at end of file diff --git a/examples/4DRooster/Blurred.jpg.sha512 b/examples/4DRooster/Blurred.jpg.sha512 new file mode 100644 index 000000000..12ddb174f --- /dev/null +++ b/examples/4DRooster/Blurred.jpg.sha512 @@ -0,0 +1 @@ +52a2f4d2559e0a60c5af7c50c5969ea4f9649573492df3f25e7e963b408bc25022ea8aec3a661ae1cf6a82119a249ef4c383d97c483371b0a5fa578314457c2d \ No newline at end of file diff --git a/examples/4DRooster/MotionMask.jpg.sha512 b/examples/4DRooster/MotionMask.jpg.sha512 new file mode 100644 index 000000000..c20004df8 --- /dev/null +++ b/examples/4DRooster/MotionMask.jpg.sha512 @@ -0,0 +1 @@ +13b801ac5e9833613e42abb9d9e7628fcde6f4dcf2d7e6c2a501e61a6b335d7b6766c512db547887a367f280943a9715934f6f21918117942e6e3ea99e39865d \ No newline at end of file diff --git a/examples/4DRooster/README.md b/examples/4DRooster/README.md new file mode 100644 index 000000000..ca7aea524 --- /dev/null +++ b/examples/4DRooster/README.md @@ -0,0 +1,225 @@ +# 4DROOSTER: Total variation-regularized 3D + time reconstruction + +RTK provides a tool to reconstruct a 3D + time image which does not require explicit motion estimation. Instead, it uses total variation regularization along both space and time. The implementation is based on a paper that we have published ([article](http://www.creatis.insa-lyon.fr/site/fr/publications/MORY-14)). You should read the article to understand the basics of the algorithm before trying to use the software. + +The algorithm requires a set of projection images with the associated RTK geometry, the respiratory phase of each projection image and a motion mask in which the region of the space where movement is expected is set to 1 and the rest is set to 0. Each piece of data is described in more details below and can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776). It is assumed that the breathing motion is periodic, which implies that the mechanical state of the chest depends only on the respiratory phase. + +## Contents +-------- + +* [1 Projection images](#projection-images) +* [2 Motion mask](#motion-mask) +* [3 Respiratory signal](#respiratory-signal) +* [4 4D ROOSTER for cone-beam CT reconstruction](#rooster-for-conebeam-ct-reconstruction) +* [5 Motion-Aware 4D ROOSTER (MA-ROOSTER)](#motion-aware-4d-rooster) + +## Projection images +================= + +This example is illustrated with a set of projection images of the [POPI patient](http://www.creatis.insa-lyon.fr/rio/popi-model_original_page). You can [download the projections](https://data.kitware.com/api/v1/item/5be99af88d777f2179a2e144/download) and the required tables of the Elekta database, [FRAME.DBF](https://data.kitware.com/api/v1/item/5be99a068d777f2179a2cf4f/download) and [IMAGE.DBF](https://data.kitware.com/api/v1/item/5be99a078d777f2179a2cf65/download). The dataset is first used to reconstruct a blurry image: + +``` +# Convert Elekta database to RTK geometry +rtkelektasynergygeometry + -o geometry.rtk + -f FRAME.DBF + -i IMAGE.DBF + -u 1.3.46.423632.141000.1169042526.68 + +# Reconstruct a 3D volume from all projection images. +# We implicitly assume that the patient's chest is static, which is wrong. +# Therefore the reconstruction will be blurry. This blurry reconstruction is +# not required for 4D ROOSTER, but there is usually no other way to generate +# a motion mask, which is required. Additionally, the blurry reconstruction +# can be used to initialize the 4D ROOSTER reconstruction. +rtkfdk + -p . + -r .*.his + -o fdk.mha + -g geometry.rtk + --hann 0.5 + --pad 1.0 + --dimension 160 + --spacing 2 +``` + +You should obtain something like that with [VV](http://vv.creatis.insa-lyon.fr/): + +![Blurred](Blurred.jpg){w=600px alt="Blurred image"} + +## Motion mask +=========== + +The next piece of data is a 3D motion mask: a volume that contains only zeros, where no movement is expected to occur, and ones, where movement is expected. Typically, during breathing, the whole chest moves, so it is safe to use the [patient mask](https://data.kitware.com/api/v1/item/5be99a418d777f2179a2ddf4/download) (red+green). However, restricting the movement to a smaller region, e.g. the rib cage, can help reconstruct this region more accurately and mitigate artifacts, so we might want to use the [rib cage mask](https://data.kitware.com/api/v1/item/5be99a2c8d777f2179a2dc4f/download) (red): + +![Mm](MotionMask.jpg){w=400px alt="Motion mask"} + +## Respiratory signal +================== + +The 4D ROOSTER algorithm requires that we associate each projection image with the instant of the respiratory cycle at which it has been acquired. We used the Amsterdam shroud solution of Lambert Zijp (described [here](http://www.creatis.insa-lyon.fr/site/fr/publications/RIT-12a)) which is implemented in RTK + +``` +rtkamsterdamshroud --path . + --regexp '.*.his' + --output shroud.mha + --unsharp 650 +rtkextractshroudsignal --input shroud.mha + --output signal.txt + --phase sphase.txt +``` + +to get the phase signal. Note that the phase must go from 0 to 1 where 0.3 corresponds to 30% in the respiratory cycle, i.e., frame 3 if you have a 10-frames 4D reconstruction or frame 6 if you have a 20-frames 4D reconstruction. The [resulting phase](https://data.kitware.com/api/v1/item/5be99af98d777f2179a2e160/download) is in green on top of the blue respiratory signal and the detected end-exhale peaks: + +![Signal](Signal.jpg){w=800px alt="Phase signal"} + +## ROOSTER for conebeam CT reconstruction +========================================== + +We now have all the pieces to perform a 3D + time reconstruction. The algorithm will perform "niter" iterations of the main loop, and run three inner loops at each iteration: + +* conjugate gradient reconstruction with "cgiter" iterations +* total-variation regularization in space, with parameter "gamma_space" (the higher, the more regularized) and "tviter" iterations +* total-variation regularization in time, with parameter "gamma_time" (the higher, the more regularized) and "tviter" iterations + +The number of iterations suggested here should work in most situations. The parameters gamma_space and gamma_time, on the other hand, must be adjusted carefully for each type of datasets (and, unfortunately, for each resolution). Unlike in analytical reconstruction methods like FDK, with 4D ROOSTER it is also very important that the reconstruction volume contains the whole patient's chest. You should therefore adapt the "dimension", "spacing" and "origin" parameters carefully, based on what you have observed on the blurry FDK reconstruction. + +``` +# Reconstruct from all projection images with 4D ROOSTER +rtkfourdrooster + -p . + -r .*.his + -o rooster.mha + -g geometry.rtk + --signal sphase.txt + --motionmask MotionMask.mha + --gamma_time 0.0001 + --gamma_space 0.0001 + --niter 30 + --cgiter 4 + --tviter 10 + --spacing 2 + --dimension 160 + --frames 5 +``` + +Depending on the resolution you choose, and even on powerful computers, the reconstruction time can range from a few minutes (very low resolution, typically 32³ * 5, only for tests) to many hours (standard resolution, typically 256³ * 10). To speed things up, it is recommended to use the CUDA version of the forward and back projectors by running this command instead: + +``` +# Reconstruct from all projection images with 4D ROOSTER, using CUDA forward and back projectors +rtkfourdrooster + -p . + -r .*.his + -o rooster.mha + -g geometry.rtk + --signal sphase.txt + --motionmask MotionMask.mha + --fp CudaRayCast + --bp CudaVoxelBased + --gamma_time 0.0001 + --gamma_space 0.0001 + --niter 30 + --cgiter 4 + --tviter 10 + --spacing 2 + --dimension 160 + --frames 5 +``` + +With a recent GPU, this should allow you to perform a standard resolution reconstruction in less than one hour. + +Note that the reconstructed volume in this example does not fully contain the attenuating object, causing hyper-attenuation artifacts on the borders of the result. To avoid these artifacts, reconstruct a larger volume (--dimension 256) should be fine. Note that you will have to resize your motion mask as well, as 3D the motion mask is expected to have the same size, spacing and origin as the first 3 dimensions of the 4D output. + +## Motion-Aware 4D Rooster +========================== + +4D ROOSTER doesn't require explicit motion information, but can take advantage of it to guide TV-regularization if motion information is available. Please refer to the tutorial on Motion-Compensated FDK to learn how to obtain a valid 4D Displacement Vector Field (DVF). Once you have it, simply adding the option --dvf "4D_DVF_Filename" to the list of rtkfourdrooster arguments will run MA-ROOSTER instead of 4D ROOSTER. + +``` +# Reconstruct from all projection images with MA ROOSTER +rtkfourdrooster + -p . + -r .*.his + -o rooster.mha + -g geometry.rtk + --signal sphase.txt + --motionmask MotionMask.mha + --gamma_time 0.0001 + --gamma_space 0.0001 + --niter 30 + --cgiter 4 + --tviter 10 + --spacing 2 + --dimension 160 + --frames 5 + --dvf deformationField_4D.mhd +``` + +Making use of the motion information adds to computation time. If you have compiled RTK with RTK_USE_CUDA = ON and have a working and CUDA-enabled nVidia GPU, it will automatically be used to speed up that part of the process. + +The article in which the theoretical foundations of MA-ROOSTER are presented (link to be added) contains results obtained on two patients. If you wish to reproduce these results, you can download all the necessary data here: + +* Original projections, log-transformed projections with the table removed, motion mask, respiratory signal, and transform matrices to change from CT to CBCT coordinates and back + * [Patient 1](https://data.kitware.com/api/v1/item/5be97d688d777f2179a28e39/download) + * [Patient 2](https://data.kitware.com/api/v1/item/5be99df68d777f2179a2e904/download) +* Inverse-consistent 4D Displacement Vector Fields, to the end-exhale phase and from the end-exhale phase + * [Patient 1](https://data.kitware.com/api/v1/item/5be989e68d777f2179a29e95/download) + * [Patient 2](https://data.kitware.com/api/v1/item/5be9a1388d777f2179a2f44d/download) + +Extract the data of each patient in a separate folder. From the folder containing the data of patient 1, run the following command line: + +``` +# Reconstruct patient 1 with MA ROOSTER +rtkfourdrooster + -p . + -r correctedProjs.mha + -o marooster.mha + -g geom.xml + --signal cphase.txt + --motionmask dilated_resampled_mm.mhd + --gamma_time 0.0002 + --gamma_space 0.00005 + --niter 10 + --cgiter 4 + --tviter 10 + --spacing "1, 1, 1, 1" + --dimension "220, 280, 370, 10" + --origin "-140, -140, -75, 0" + --frames 10 + --dvf toPhase50_4D.mhd + --idvf fromPhase50_4D.mhd +``` + +From the folder containing the data of patient 2, run the following command line (only the dimension and origin parameters are different): + +``` +# Reconstruct patient 2 with MA ROOSTER +rtkfourdrooster + -p . + -r correctedProjs.mha + -o marooster.mha + -g geom.xml + --signal cphase.txt + --motionmask dilated_resampled_mm.mhd + --gamma_time 0.0002 + --gamma_space 0.00005 + --niter 10 + --cgiter 4 + --tviter 10 + --spacing "1, 1, 1, 1" + --dimension "285, 270, 307, 10" + --origin "-167.5, -135, -205, 0" + --frames 10 + --dvf toPhase50_4D.mhd + --idvf fromPhase50_4D.mhd +``` + +Note the option "--idvf", which allows to provide the inverse DVF. It is used to inverse warp the 4D reconstruction after the temporal regularization. MA-ROOSTER will work with and without the inverse DVF, and yield almost the same results in both cases. Not using the inverse DVF is approximately two times slower, as it requires MA-ROOSTER to perform the inverse warping by an iterative method. + +Again, if you have a CUDA-enabled GPU (in this case with at least 3 GB of VRAM), and have compiled RTK with RTK_USE_CUDA = ON, you can add the "--bp CudaVoxelBased" and "--fp CudaRayCast" to speed up the computation by performing the forward and back projections on the GPU. + +You do not need the 4D planning CT data to perform the MA-ROOSTER reconstructions. It is only required to compute the DVFs, which can be downloaded above. We do provide it anyway, in case you want to use your own method, or the one described in Motion-Compensated FDK, to extract a DVF from it: + +* 4D planning CT + * [Patient 1](https://data.kitware.com/api/v1/item/5be98bd28d777f2179a2a279/download) + * [Patient 2](https://data.kitware.com/api/v1/item/5be9a1918d777f2179a2f568/download) diff --git a/examples/4DRooster/Signal.jpg.sha512 b/examples/4DRooster/Signal.jpg.sha512 new file mode 100644 index 000000000..c46a3d097 --- /dev/null +++ b/examples/4DRooster/Signal.jpg.sha512 @@ -0,0 +1 @@ +7c265ec19e50310585f65ae90147ca1972225f93f35b090e233acfa081b084e329ed5b2d31fad23ec43e64d281937fad213d1229f5b717a98592518f6d56a7c2 \ No newline at end of file diff --git a/examples/AmsterdamShroud/Amsterdam.png.sha512 b/examples/AmsterdamShroud/Amsterdam.png.sha512 new file mode 100644 index 000000000..0feab162e --- /dev/null +++ b/examples/AmsterdamShroud/Amsterdam.png.sha512 @@ -0,0 +1 @@ +bf4f5bc8e887d7faf7dc9e835814b123a84dd617ef83343590411651f35a03083cc2d71cff6f7a3d8fd5b4b0fdd517871335650183238814de769f49ca07210e \ No newline at end of file diff --git a/examples/AmsterdamShroud/Moving-Phantom-Sinogram.png.sha512 b/examples/AmsterdamShroud/Moving-Phantom-Sinogram.png.sha512 new file mode 100644 index 000000000..0993939b3 --- /dev/null +++ b/examples/AmsterdamShroud/Moving-Phantom-Sinogram.png.sha512 @@ -0,0 +1 @@ +b31fe75d2f3ef6289e6dbca4bf7dc2675eebd80bfd25d61b167727587b796cf9b371f05db51949904699ceb3bf5cc681265ac7b4aa9b89ab632a2213e7bcca1d \ No newline at end of file diff --git a/examples/AmsterdamShroud/README.md b/examples/AmsterdamShroud/README.md new file mode 100644 index 000000000..ee9e5e551 --- /dev/null +++ b/examples/AmsterdamShroud/README.md @@ -0,0 +1,13 @@ +# Amsterdam Shroud + +![sin](Moving-Phantom-Sinogram.png){w=400px alt="Moving phantom sinogram"} +![img](Amsterdam.png){w=400px alt="Amsterdam image"} + +Picture 1 shows the sinogram of the input and picture 2 the shroud image that was created using the command line below. + +The script uses the file [movingPhantom.mha](https://data.kitware.com/api/v1/file/5be99c428d777f2179a2e537/download) as input: + +``` +# Creating an Amsterdam Shroud image from a set of projections +rtkamsterdamshroud -p . -r movingPhantom.mha -o Amsterdam.mha + ``` \ No newline at end of file diff --git a/examples/ConjugateGradient/Code3D.sh b/examples/ConjugateGradient/Code3D.sh new file mode 100755 index 000000000..9a6b8c7b8 --- /dev/null +++ b/examples/ConjugateGradient/Code3D.sh @@ -0,0 +1,21 @@ + # Create a simulated geometry + rtksimulatedgeometry -n 180 -o geometry.xml + # You may add "--arc 200" to make the scan short or "--proj_iso_x 200" to offset the detector + + # Create projections of the phantom file + rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 + + # Reconstruct + rtkconjugategradient -p . -r projections.mha -o 3dcg.mha -g geometry.xml --spacing 2 --dimension 256 -n 20 + + # Create a reference volume for comparison + rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha + + # Perform least squares reconstruction + rtkconjugategradient -p . -r noisyLineIntegrals.mha -o LeastSquares.mha -g geom.xml -n 20 + +# # Perform weighted least squares reconstruction +# rtkconjugategradient -p . -r noisyLineIntegrals.mha -o WeightedLeastSquares.mha -g geom.xml -w weightsmap.mha -n 20 + +# # Perform preconditioned conjugate gradient reconstruction with weighted least squares cost function +# rtkconjugategradient -p . -r noisyLineIntegrals.mha -o WeightedLeastSquares.mha -g geom.xml -w weightsmap.mha -n 20 --preconditioned \ No newline at end of file diff --git a/examples/ConjugateGradient/ConjugateGradient-Sinogram.png.sha512 b/examples/ConjugateGradient/ConjugateGradient-Sinogram.png.sha512 new file mode 100644 index 000000000..1f4321204 --- /dev/null +++ b/examples/ConjugateGradient/ConjugateGradient-Sinogram.png.sha512 @@ -0,0 +1 @@ +6a1b6bd8a92f84cbd9a7e7385216f15bef58c814e182c3598a0938f2b4d839b23e7237cbd536d5d6902c0720c7694b6b0b47742dae7ed350306b7cc6314b910a \ No newline at end of file diff --git a/examples/ConjugateGradient/ConjugateGradient.png.sha512 b/examples/ConjugateGradient/ConjugateGradient.png.sha512 new file mode 100644 index 000000000..15f983f14 --- /dev/null +++ b/examples/ConjugateGradient/ConjugateGradient.png.sha512 @@ -0,0 +1 @@ +8fa29a5ddbe658203faea0b1b3d48146904f1e2362d44385681cb7a96e2eba5b9d8a99f26a9c73c1a8fb24fe6d19891d1db3c57489c68fecb540f522703cf6fa \ No newline at end of file diff --git a/examples/ConjugateGradient/README.md b/examples/ConjugateGradient/README.md new file mode 100644 index 000000000..d6ab0c2fe --- /dev/null +++ b/examples/ConjugateGradient/README.md @@ -0,0 +1,42 @@ +# Conjugate gradient + +--- + +This example uses the Shepp–Logan phantom. + +## 3D + +![sin](ConjugateGradient-Sinogram.png){w=200px alt="Conjugate sinogram"} +![img](ConjugateGradient.png){w=200px alt="ConjugateGradient image"} + +```shell + # Create a simulated geometry + rtksimulatedgeometry -n 180 -o geometry.xml + # You may add "--arc 200" to make the scan short or "--proj_iso_x 200" to offset the detector + + # Create projections of the phantom file + rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 + + # Reconstruct + rtkconjugategradient -p . -r projections.mha -o 3dcg.mha -g geometry.xml --spacing 2 --dimension 256 -n 20 + + # Create a reference volume for comparison + rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha +``` + +In the presence of noise, all projection data may not be equally reliable. The conjugate gradient algorithm can be modified to take this into account, and each pixel of the projections can be associated with a weight. The higher the weight, the more reliable the pixel data. Download [noisy projections](https://data.kitware.com/api/v1/item/5be99cdf8d777f2179a2e63d/download) and [the associated weights](https://data.kitware.com/api/v1/item/5be99d268d777f2179a2e6f8/download), as well as [the geometry](https://data.kitware.com/api/v1/item/5be99d268d777f2179a2e700/download), and run the following to compare the regular least squares reconstruction (without weights) and the weighted least squares reconstruction. + +```shell + # Perform least squares reconstruction + rtkconjugategradient -p . -r noisyLineIntegrals.mha -o LeastSquares.mha -g geom.xml -n 20 + + # Perform weighted least squares reconstruction + rtkconjugategradient -p . -r noisyLineIntegrals.mha -o WeightedLeastSquares.mha -g geom.xml -w weightsmap.mha -n 20 +``` + + +Taking the weights into account slows down convergence. This can be corrected by using a preconditioner in the conjugate gradient algorithm. The preconditioner is computed automatically from the weights map, you just need to activate the flag : +```shell + # Perform preconditioned conjugate gradient reconstruction with weighted least squares cost function + rtkconjugategradient -p . -r noisyLineIntegrals.mha -o WeightedLeastSquares.mha -g geom.xml -w weightsmap.mha -n 20 --preconditioned +``` \ No newline at end of file diff --git a/examples/CreateGammexPhantom/GammexPhantom.png.sha512 b/examples/CreateGammexPhantom/GammexPhantom.png.sha512 new file mode 100644 index 000000000..76e165def --- /dev/null +++ b/examples/CreateGammexPhantom/GammexPhantom.png.sha512 @@ -0,0 +1 @@ +51510a6083a3c6dfae0ffa229c0a98ac11f2bea2b14eee814658b5ef47a98ec3545d54951afd57ba103bbaa1db04ff7e7dcddfa862e78111904da7d079ee40bf \ No newline at end of file diff --git a/examples/CreateGammexPhantom/README.md b/examples/CreateGammexPhantom/README.md new file mode 100644 index 000000000..a362b6f56 --- /dev/null +++ b/examples/CreateGammexPhantom/README.md @@ -0,0 +1,10 @@ +# Create gammex phantom + +![img](GammexPhantom.png){w=400px alt="Gammex"} + +This script uses the file [Gammex.txt](https://data.kitware.com/api/v1/file/6762da8a290777363f95c293/download) as configuration file which creates a Gammex phantom. + +``` + # Create a 3D Gammex phantom + rtkdrawgeometricphantom --phantomfile Gammex.txt -o gammex.mha + ``` \ No newline at end of file diff --git a/examples/DaubechiesWaveletsRegularizedReconstruction/Code.sh b/examples/DaubechiesWaveletsRegularizedReconstruction/Code.sh new file mode 100644 index 000000000..5b4a2361a --- /dev/null +++ b/examples/DaubechiesWaveletsRegularizedReconstruction/Code.sh @@ -0,0 +1,12 @@ + # Create a simulated geometry + rtksimulatedgeometry -n 180 -o geometry.xml + # You may add "--arc 200" to make the scan short or "--proj_iso_x 200" to offset the detector + + # Create projections of the phantom file + rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 + + # Reconstruct + rtkadmmwavelets -p . -r projections.mha -o admmwavelets.mha -g geometry.xml --spacing 2 --dimension 256 --alpha 1 --beta 1000 -n 3 + + # Create a reference volume for comparison + rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha \ No newline at end of file diff --git a/examples/DaubechiesWaveletsRegularizedReconstruction/Overlay.png.sha512 b/examples/DaubechiesWaveletsRegularizedReconstruction/Overlay.png.sha512 new file mode 100644 index 000000000..15f983f14 --- /dev/null +++ b/examples/DaubechiesWaveletsRegularizedReconstruction/Overlay.png.sha512 @@ -0,0 +1 @@ +8fa29a5ddbe658203faea0b1b3d48146904f1e2362d44385681cb7a96e2eba5b9d8a99f26a9c73c1a8fb24fe6d19891d1db3c57489c68fecb540f522703cf6fa \ No newline at end of file diff --git a/examples/DaubechiesWaveletsRegularizedReconstruction/README.md b/examples/DaubechiesWaveletsRegularizedReconstruction/README.md new file mode 100644 index 000000000..c8d14e5aa --- /dev/null +++ b/examples/DaubechiesWaveletsRegularizedReconstruction/README.md @@ -0,0 +1,9 @@ +# Daubechies Wavelets Regularized Reconstruction + +![sin](Sinogram.png){w=300px alt="sinogram"} +![img](Overlay.png){w=300px alt="image"} + +This script uses the SheppLogan phantom + +```{literalinclude} Code.sh +``` \ No newline at end of file diff --git a/examples/DaubechiesWaveletsRegularizedReconstruction/Sinogram.png.sha512 b/examples/DaubechiesWaveletsRegularizedReconstruction/Sinogram.png.sha512 new file mode 100644 index 000000000..1f4321204 --- /dev/null +++ b/examples/DaubechiesWaveletsRegularizedReconstruction/Sinogram.png.sha512 @@ -0,0 +1 @@ +6a1b6bd8a92f84cbd9a7e7385216f15bef58c814e182c3598a0938f2b4d839b23e7237cbd536d5d6902c0720c7694b6b0b47742dae7ed350306b7cc6314b910a \ No newline at end of file diff --git a/examples/ElektaReconstruction/Elekta.png.sha512 b/examples/ElektaReconstruction/Elekta.png.sha512 new file mode 100644 index 000000000..dd929f0b7 --- /dev/null +++ b/examples/ElektaReconstruction/Elekta.png.sha512 @@ -0,0 +1 @@ +8cb55621686009fb44c03656aa978114ee689ef474afea0bfa20837c8b321eaedc3aebd3c6ba7592d555c215460b3c3a16034568d7c6a5f99352ca5edf65c953 \ No newline at end of file diff --git a/examples/ElektaReconstruction/README.md b/examples/ElektaReconstruction/README.md new file mode 100644 index 000000000..a759d10af --- /dev/null +++ b/examples/ElektaReconstruction/README.md @@ -0,0 +1,72 @@ +# Elekta Reconstruction + +Elekta provides easy access to raw data. The data and projection images are stored in a single directory which is user-configurable. The default location is `D:\db`. In this folder, there is a database in DBase format. Each table is contained in a `.DBF` file. RTK needs the `IMAGE.DBF` and `FRAME.DBF` tables. + +Patient data are stored in individual folders. By default, the name of each patient folder is `patient_ID` where `ID` is the patient ID. In these folders, one can access the planning CT in the `CT_SET` subfolder and the cone-beam projections in `IMAGES/img_DICOM_UID` subfolders, where `DICOM_UID` is the DICOM UID of the acquisition. The projection images are `.his` files. The reconstructed images are the `IMAGES/img_DICOM_UID/Reconstruction/*SCAN` files. + +## Reconstruction Steps + +The first step before proceeding with the reconstruction is to convert Elekta's database information into an RTK geometry file using a command-line tool. Follow these steps: + +### 1. Download Elekta Dataset + +Download the dataset from [Elekta-data](https://data.kitware.com/api/v1/item/5be973478d777f2179a26e1c/download). + +### 2. Convert Geometry + +Run the application to convert Elekta's geometry into RTK's format (DICOM_UID is contained in the subfolder name of the `.his` files): + +```bash +rtkelektasynergygeometry \ + --image_db IMAGE.DBF \ + --frame_db FRAME.DBF \ + --dicom_uid 1.3.46.423632.135428.1351013645.166 \ + -o elektaGeometry +``` + +Since XVI v5, the geometry is contained in a separate `_Frames.xml` file, which can be used as follows: + +```bash +rtkelektasynergygeometry \ + --xml _Frames.xml \ + -o elektaGeometry +``` + +An example of such a file is available in the test data [here](https://data.kitware.com/api/v1/item/5b179c898d777f15ebe201fd/download). + +### 3. Reconstruct Using RTK Applications + +Use the `rtkfdk` algorithm to reconstruct a single axial slice (e.g., slice 29.5) of the volume: + +```bash +rtkfdk \ + --lowmem \ + --geometry elektaGeometry \ + --path img_1.3.46.423632.135428.1351013645.166/ \ + --regexp '.*.his' \ + --output slice29.5.mha \ + --verbose \ + --spacing 0.25,0.25,0.25 \ + --dimension 1024,1,1024 \ + --origin -127.875,29.5,-127.875 +``` + +### 4. Apply the FOV Filter + +Apply the field-of-view (FOV) filter to mask out everything outside the FOV: + +```bash +rtkfieldofview \ + --geometry elektaGeometry \ + --path img_1.3.46.423632.135428.1351013645.166/ \ + --regexp '.*.his' \ + --reconstruction slice29.5.mha \ + --output slice29.5.mha \ + --verbose +``` + +### 5. Visualize the Result + +You can visualize the result using a viewer (e.g., VV). The resulting image should look like the following: + +![Elekta.jpg](Elekta.png){w=400px alt="Elekta snapshot"} diff --git a/examples/ForwardProjection/Code.sh b/examples/ForwardProjection/Code.sh new file mode 100644 index 000000000..f8b6a584e --- /dev/null +++ b/examples/ForwardProjection/Code.sh @@ -0,0 +1,8 @@ + # Create a simulated geometry + rtksimulatedgeometry -n 180 -o geometry.xml + + # Forward project + rtkforwardprojections -g geometry.xml -o projections.mha -i 00.mhd --spacing 2 --dimension 512 + + # Reconstruct in the same resolution as the original + rtkfdk -p . -r projections.mha -o fdk.mha -g geometry.xml --spacing\=0.976562,0.976562,2 --origin\=\-250,-250,-164.5 --dimension\=512,512,141 \ No newline at end of file diff --git a/examples/ForwardProjection/POPI-Reconstruction.png.sha512 b/examples/ForwardProjection/POPI-Reconstruction.png.sha512 new file mode 100644 index 000000000..8eab64487 --- /dev/null +++ b/examples/ForwardProjection/POPI-Reconstruction.png.sha512 @@ -0,0 +1 @@ +5261406f0891c4141006a4a4f75aaabe81e372dc380c27230b8b7355d13392ab43abc351de028f880c8df2bfdf352a5754bc8f575d6365fb4c9b66b48c547f37 \ No newline at end of file diff --git a/examples/ForwardProjection/POPI-Sinogram.png.sha512 b/examples/ForwardProjection/POPI-Sinogram.png.sha512 new file mode 100644 index 000000000..408224629 --- /dev/null +++ b/examples/ForwardProjection/POPI-Sinogram.png.sha512 @@ -0,0 +1 @@ +72cd6bd4bed25bb626ee5eee5aa142c3514b0ac8346def3887ba8eece4f9fc0379cd6d9a1e085848fb5e7ef2ed3ef9d1c322c9e9cfa14f49986a329956a24289 \ No newline at end of file diff --git a/examples/ForwardProjection/README.md b/examples/ForwardProjection/README.md new file mode 100644 index 000000000..838087e30 --- /dev/null +++ b/examples/ForwardProjection/README.md @@ -0,0 +1,13 @@ +# Forward Projection + +![sin](POPI-Sinogram.png){w=400px alt="POPI sinogram"} +![img](POPI-Reconstruction.png){w=400px alt="POPI reconstruction"} + +This script uses the files [00.mhd](http://www.creatis.insa-lyon.fr/~srit/POPI/MedPhys11/bl/mhd/00.mhd) and [00.raw](http://www.creatis.insa-lyon.fr/~srit/POPI/MedPhys11/bl/mhd/00.raw) of the [POPI](http://www.creatis.insa-lyon.fr/rio/popi-model/) as input. + +```{literalinclude} Code.sh +``` + +Note that the original file is in Hounsfield units which explains the negative values in the projection images since, e.g., the attenuation of air is -1000 HU. + +It is also worth of note that the file is oriented in the DICOM coordinate system although RTK uses the IEC 61217 which results in a rotation around the antero-posterior axis of the patient. This can be easily changed by modifying the TransformMatrix in the 00.mhd file. diff --git a/examples/MotionCompensationReconstruction/Blurred.jpg.sha512 b/examples/MotionCompensationReconstruction/Blurred.jpg.sha512 new file mode 100644 index 000000000..12ddb174f --- /dev/null +++ b/examples/MotionCompensationReconstruction/Blurred.jpg.sha512 @@ -0,0 +1 @@ +52a2f4d2559e0a60c5af7c50c5969ea4f9649573492df3f25e7e963b408bc25022ea8aec3a661ae1cf6a82119a249ef4c383d97c483371b0a5fa578314457c2d \ No newline at end of file diff --git a/examples/MotionCompensationReconstruction/Blurred_vs_mc.gif.sha512 b/examples/MotionCompensationReconstruction/Blurred_vs_mc.gif.sha512 new file mode 100644 index 000000000..6ac501ad4 --- /dev/null +++ b/examples/MotionCompensationReconstruction/Blurred_vs_mc.gif.sha512 @@ -0,0 +1 @@ +2c05f3a1f4c14db8b0280e59640aec0c6a856dbecba1e2d3db3b5db9a202cdfe2f40982531f97174652c210bea120b9867b83ca7e4110a48fe28ecf2f3ce776e \ No newline at end of file diff --git a/examples/MotionCompensationReconstruction/MotionMask.jpg.sha512 b/examples/MotionCompensationReconstruction/MotionMask.jpg.sha512 new file mode 100644 index 000000000..c20004df8 --- /dev/null +++ b/examples/MotionCompensationReconstruction/MotionMask.jpg.sha512 @@ -0,0 +1 @@ +13b801ac5e9833613e42abb9d9e7628fcde6f4dcf2d7e6c2a501e61a6b335d7b6766c512db547887a367f280943a9715934f6f21918117942e6e3ea99e39865d \ No newline at end of file diff --git a/examples/MotionCompensationReconstruction/README.md b/examples/MotionCompensationReconstruction/README.md new file mode 100644 index 000000000..684dd4b32 --- /dev/null +++ b/examples/MotionCompensationReconstruction/README.md @@ -0,0 +1,169 @@ +# Motion Compensation Reconstruction + +RTK provides the necessary tools to reconstruct an image with motion compensation. The implementation is based on two articles that we have published ([article 1](https://hal.archives-ouvertes.fr/hal-00443440) and [article 2](https://hal.archives-ouvertes.fr/hal-01967313)) but only the FDK-based motion-compensated CBCT reconstruction (analytic algorithm in article 1) and without optimization (very slow reconstruction compared to article 2). You should read the articles to understand the basics of the algorithm before trying to use the software. + +The algorithm requires a set of projection images with the associated RTK geometry, the respiratory phase of each projection image and the 4D motion vector field over a respiratory cycle in the cone-beam coordinate system. Each piece of data is described in more details below and can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776). It is assumed that we have a breathing motion that is cyclic and similar to that described by the vector field. Note that you could modify the code and create your own motion model if you want to, in which case you should probably [contact us](http://www.openrtk.org/RTK/project/contactus.html). + +## Contents + +- [1 Projection images](#projection-images) +- [2 Deformation vector field](#deformation-vector-field) +- [3 Respiratory signal](#respiratory-signal) +- [4 Motion-compensated cone-beam CT reconstruction](#motion-compensated-cone-beam-ct-reconstruction) + +### Projection images + +This example is illustrated with a set of projection images of the [POPI patient](http://www.creatis.insa-lyon.fr/rio/popi-model_original_page). This dataset has been used in the first previously-mentioned article. You can [download the projections](https://data.kitware.com/api/v1/item/5be99af88d777f2179a2e144/download) and the required tables of the Elekta database, [FRAME.DBF](https://data.kitware.com/api/v1/item/5be99a068d777f2179a2cf4f/download) and [IMAGE.DBF](https://data.kitware.com/api/v1/item/5be99a078d777f2179a2cf65/download). The dataset is first used to reconstruct a blurry image: + +```bash +# Convert Elekta database to RTK geometry +rtkelektasynergygeometry \ + -o geometry.rtk \ + -f FRAME.DBF \ + -i IMAGE.DBF \ + -u 1.3.46.423632.141000.1169042526.68 + +# Reconstruct from all projection images without any motion compensation +rtkfdk \ + -p . \ + -r .*.his \ + -o fdk.mha \ + -g geometry.rtk \ + --hann 0.5 \ + --pad 1.0 + +# Keep only the field-of-view of the image +rtkfieldofview \ + --reconstruction fdk.mha \ + --output fdk.mha \ + --geometry geometry.rtk \ + --path . \ + --regexp '.*.his' +``` + +You should obtain something like that with [VV](http://vv.creatis.insa-lyon.fr/): + +![Blurred](Blurred.jpg){w=600px alt="Blurred image"} + +### Deformation vector field + +The next piece of data is a 4D deformation vector field that describes a respiratory cycle. Typically, it can be obtained from the 4D planning CT with deformable image registration. Here, I have used [Elastix](http://elastix.lumc.nl/) with the [sliding module](http://elastix.lumc.nl/modelzoo/par0016) developed by Vivien Delmon. The registration uses a [patient mask](https://data.kitware.com/api/v1/item/5be99a408d777f2179a2dde8/download) (red+green) and a [motion mask](https://data.kitware.com/api/v1/item/5be99a088d777f2179a2cf6f/download) (red) as described in [Jef's publication](http://www.creatis.insa-lyon.fr/site/fr/publications/VAND-12): + +![Mm](MotionMask.jpg){w=400px alt="Motion mask"} + +The registration can easily be scripted, here with bash, where each phase image of the POPI 4D CT has been stored in files 00.mhd to 50.mhd: + +```bash +for i in $(seq -w 0 10 90) +do + mkdir $i + elastix -f 50.mhd \ + -m $i.mhd \ + -out $i \ + -labels mm_50.mha \ + -fMask patient_50.mha \ + -p Par0016.multibsplines.lung.sliding.txt +done +``` + +Deformable Image Registration is a complex and long process so you will have to be patient here. Note that the reference frame is phase 50% and it is registered to each phase from 0% to 90%. One subtle step is that the vector field is a displacement vector field, i.e., each vector is the local displacement of the point at its location. Since I ran the registration on the 4D planning CT, the coordinate system is not that of the cone-beam CT. In order to produce the vector field in the cone-beam coordinate system, I have used the following bash script that combines transformix and several "clitk" tools that are provided along with VV: + +```bash +# Create 4x4 matrix that describes the CT to CBCT change of coordinate system. +# This matrix is a combination of the knowledge of the isocenter position / axes orientation +# and a rigid alignment that has been performed with Elastix +echo "-0.0220916855767852 0.9996655273534405 -0.0134458487848415 -83.6625731437426197" >CT_CBCT.mat +echo " 0.0150924269790251 -0.0131141301144939 -0.9998000991394341 -4.0763571826687057" >>CT_CBCT.mat +echo " 0.9996420239647088 0.0222901999207823 0.0147976657359281 77.8903364738220034" >>CT_CBCT.mat +echo " 0.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000" >>CT_CBCT.mat + +# Transform 4x4 matrix that describes the transformation +# from planning CT to CBCT to a vector field +clitkMatrixTransformToVF --like 50.mhd \ + --matrix CT_CBCT.mat \ + --output CT_CBCT.mha + +# Inverse transformation. Also remove upper slices that are outside the +# planning CT CBCT_CT.mat is the inverse of CT_CBCT.mha +clitkMatrixInverse -i CT_CBCT.mat \ + -o CBCT_CT.mat +clitkMatrixTransformToVF --origin -127.5,-107.5,-127.5 \ + --spacing 1,1,1 \ + --size 256,236,256 \ + --matrix CBCT_CT.mat \ + --output CBCT_CT.mha + +# Go over each elastix output file, generate the vector field with +# transformix and compose with the two rigid vector fields +for i in $(seq -w 0 10 90) +do + transformix -in 50.mhd \ + -out $i \ + -tp $i/TransformParameters.0.txt \ + -def all -threads 16 + clitkComposeVF --input1 CBCT_CT.mha \ + --input2 $i/deformationField.mhd \ + --output $i/deformationField.mhd + clitkComposeVF --input1 $i/deformationField.mhd \ + --input2 CT_CBCT.mha \ + --output $i/deformationField.mhd +done +``` + +This is a bit complicated and there are probably other ways of doing this. For example, Vivien has resampled the planning CT frames on the CBCT coordinate system before doing the registrations, in which case you do not need to do all this. Just pick one of your choice but motion-compensated CBCT reconstruction requires a 4D vector field that is nicely displayed on top of a CBCT image, for example the fdk.mha that has been produced in the first step (the vector field is downsampled and displayed with VV): + +![Vf](VectorField.gif){w=400px alt="Vector field"} + +The elastix output files and the transformed 4D DVF are available [here](https://data.kitware.com/api/v1/item/5be99a058d777f2179a2cf42/download). + +### Respiratory signal + +The motion model requires that we associate each projection image with one frame of the 4D vector field. We used the Amsterdam shroud solution of Lambert Zijp (described [here](http://www.creatis.insa-lyon.fr/site/fr/publications/RIT-12a)) which is implemented in RTK + +```bash +rtkamsterdamshroud --path . \ + --regexp '.*.his' \ + --output shroud.mha \ + --unsharp 650 +rtkextractshroudsignal --input shroud.mha \ + --output signal.txt +``` + +Post-process with Matlab to obtain the phase signal, ensuring the phase ranges from 0 to 1 (e.g., 0.3 corresponds to 30% of the respiratory cycle). The resulting phase is visualized [here](https://data.kitware.com/api/v1/item/5be99af98d777f2179a2e160/download): + +![Signal](Signal.jpg){w=800px alt="Phase signal"} + +--- + +## Motion-Compensated Cone-Beam CT Reconstruction + +Gather all the pieces to perform motion-compensated reconstruction. Use the following commands: + +### Reconstruct with Motion Compensation +```bash +rtkfdk \ + -p . \ + -r .*.his \ + -o fdk.mha \ + -g geometry.rtk \ + --hann 0.5 \ + --pad 1.0 \ + --signal sphase.txt \ + --dvf deformationField_4D.mhd +``` + +### Apply the Field-of-View Filter +```bash +rtkfieldofview \ + --reconstruction fdk.mha \ + --output fdk.mha \ + --geometry geometry.rtk \ + --path . \ + --regexp '.*.his' +``` + +Toggle between uncorrected and motion-compensated reconstruction to appreciate the improvement: + +![Blurred vs mc.gif](Blurred_vs_mc.gif){w=400 alt="blurred vs motion compensation image"} + +The 4D vector field is constructed with phase 50% as a reference. Modify the reference image to reconstruct other phases, such as the time-average position. \ No newline at end of file diff --git a/examples/MotionCompensationReconstruction/Signal.jpg.sha512 b/examples/MotionCompensationReconstruction/Signal.jpg.sha512 new file mode 100644 index 000000000..c46a3d097 --- /dev/null +++ b/examples/MotionCompensationReconstruction/Signal.jpg.sha512 @@ -0,0 +1 @@ +7c265ec19e50310585f65ae90147ca1972225f93f35b090e233acfa081b084e329ed5b2d31fad23ec43e64d281937fad213d1229f5b717a98592518f6d56a7c2 \ No newline at end of file diff --git a/examples/MotionCompensationReconstruction/VectorField.gif.sha512 b/examples/MotionCompensationReconstruction/VectorField.gif.sha512 new file mode 100644 index 000000000..847eb63f1 --- /dev/null +++ b/examples/MotionCompensationReconstruction/VectorField.gif.sha512 @@ -0,0 +1 @@ +def878b887fa050f59978b9441050ea177736d1e5a227e8a6f3693ef4f8c9c7175e06cdaf71305958ed8b131686fea184fe2fb6810228fa11d69fa7b99759cb2 \ No newline at end of file diff --git a/examples/RayBoxIntersection/Code.sh b/examples/RayBoxIntersection/Code.sh new file mode 100644 index 000000000..5107895e6 --- /dev/null +++ b/examples/RayBoxIntersection/Code.sh @@ -0,0 +1,11 @@ + # Create a simulated geometry + rtksimulatedgeometry -n 360 -o geometry.xml + + # Create a box volume + rtkdrawgeometricphantom -o box.mha --spacing 1 --dimension 90 --phantomfile box.txt + + # Calculate radiological path of the box (ray intersection length) + rtkrayboxintersection -g geometry.xml -i box.mha -o rayboxintersection.mha --spacing 1 --dimension 256 + + # Reconstruct/Backproject the set of projections + rtkfdk -g geometry.xml -p . -r rayboxintersection.mha -o fdk.mha --spacing 1 --dimension 256 \ No newline at end of file diff --git a/examples/RayBoxIntersection/README.md b/examples/RayBoxIntersection/README.md new file mode 100644 index 000000000..2097404d8 --- /dev/null +++ b/examples/RayBoxIntersection/README.md @@ -0,0 +1,9 @@ +# Ray Box Intersection + +![sin](RayBox-Sinogram.png){w=400px alt="RayBox sinogram"} +![img](RayBox.png){w=400px alt="RayBox image"} + +This script uses the file [Box.txt](https://data.kitware.com/api/v1/file/67629a61290777363f95c273/download) as input. + +```{literalinclude} Code.sh +``` \ No newline at end of file diff --git a/examples/RayBoxIntersection/RayBox-Sinogram.png.sha512 b/examples/RayBoxIntersection/RayBox-Sinogram.png.sha512 new file mode 100644 index 000000000..1a4af28a2 --- /dev/null +++ b/examples/RayBoxIntersection/RayBox-Sinogram.png.sha512 @@ -0,0 +1 @@ +388cef2ce2c16b878a2100b87d27b85723019ba6e7f3137d675fe19c7625a61b194dd6d41b07d644c2e78099b33906af96bd45b592bd3a8c70e8935826226adf \ No newline at end of file diff --git a/examples/RayBoxIntersection/RayBox.png.sha512 b/examples/RayBoxIntersection/RayBox.png.sha512 new file mode 100644 index 000000000..dae9fada8 --- /dev/null +++ b/examples/RayBoxIntersection/RayBox.png.sha512 @@ -0,0 +1 @@ +fd0352b723b61dd7555c53aed396e2520fcc145628c71db7143c34e4be06858daa8091aa0492b5fc6a65b0ad3b784bb880c881c375da2b58a17546f89368edf6 \ No newline at end of file diff --git a/examples/TotalVariationRegularizedReconstruction/Code.sh b/examples/TotalVariationRegularizedReconstruction/Code.sh new file mode 100644 index 000000000..dc7d15c83 --- /dev/null +++ b/examples/TotalVariationRegularizedReconstruction/Code.sh @@ -0,0 +1,12 @@ + # Create a simulated geometry + rtksimulatedgeometry -n 180 -o geometry.xml + # You may add "--arc 200" to make the scan short or "--proj_iso_x 200" to offset the detector + + # Create projections of the phantom file + rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 + + # Reconstruct + rtkadmmtotalvariation -p . -r projections.mha -o admmtv.mha -g geometry.xml --spacing 2 --dimension 256 --alpha 1 --beta 1000 -n 3 + + # Create a reference volume for comparison + rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha \ No newline at end of file diff --git a/examples/TotalVariationRegularizedReconstruction/Overlay.png.sha512 b/examples/TotalVariationRegularizedReconstruction/Overlay.png.sha512 new file mode 100644 index 000000000..15f983f14 --- /dev/null +++ b/examples/TotalVariationRegularizedReconstruction/Overlay.png.sha512 @@ -0,0 +1 @@ +8fa29a5ddbe658203faea0b1b3d48146904f1e2362d44385681cb7a96e2eba5b9d8a99f26a9c73c1a8fb24fe6d19891d1db3c57489c68fecb540f522703cf6fa \ No newline at end of file diff --git a/examples/TotalVariationRegularizedReconstruction/README.md b/examples/TotalVariationRegularizedReconstruction/README.md new file mode 100644 index 000000000..f459fab72 --- /dev/null +++ b/examples/TotalVariationRegularizedReconstruction/README.md @@ -0,0 +1,9 @@ +# Total Variation Regularized Reconstruction + +![sin](Sinogram.png){w=300px alt="sinogram"} +![img](Overlay.png){w=300px alt="image"} + +This script uses the SheppLogan phantom + +```{literalinclude} Code.sh +``` \ No newline at end of file diff --git a/examples/TotalVariationRegularizedReconstruction/Sinogram.png.sha512 b/examples/TotalVariationRegularizedReconstruction/Sinogram.png.sha512 new file mode 100644 index 000000000..1f4321204 --- /dev/null +++ b/examples/TotalVariationRegularizedReconstruction/Sinogram.png.sha512 @@ -0,0 +1 @@ +6a1b6bd8a92f84cbd9a7e7385216f15bef58c814e182c3598a0938f2b4d839b23e7237cbd536d5d6902c0720c7694b6b0b47742dae7ed350306b7cc6314b910a \ No newline at end of file diff --git a/examples/VarianReconstruction/README.md b/examples/VarianReconstruction/README.md new file mode 100644 index 000000000..f4ada416d --- /dev/null +++ b/examples/VarianReconstruction/README.md @@ -0,0 +1,116 @@ +# Varian Reconstruction +----------- + +## Varian OBI Reconstruction + +The first step before proceeding with reconstruction is to convert Varian's geometry into RTK's format using a command-line tool. Follow these simple steps: + +### 1. Download Varian Dataset + +Download the dataset from [Varian-data](https://data.kitware.com/api/v1/item/5be94de88d777f2179a24de0/download). + +### 2. Convert Geometry + +Run the application to convert Varian's geometry into RTK's format: + +```bash +rtkvarianobigeometry \ + --xml_file ProjectionInfo.xml \ + --path Scan0/ \ + --regexp Proj_.*.hnd \ + -o geometry.xml +``` + +### 3. Reconstruct Using RTK Applications + +Reconstruct a slice (e.g., slice 30) of the volume using the `rtkfdk` algorithm: + +```bash +rtkfdk \ + --geometry geometry.xml \ + --regexp .*\.hnd \ + --path Scan0 \ + --output slice30.mha \ + --verbose \ + --spacing 0.25,0.25,0.25 \ + --dimension 1024,1,1024 \ + --origin -127.875,30,-127.875 +``` + +### 4. Apply the FOV Filter + +Apply the field-of-view (FOV) filter to discard everything outside the FOV: + +```bash +rtkfieldofview \ + --geometry geometry.xml \ + --regexp .*\.hnd \ + --path Scan0 \ + --reconstruction slice30.mha \ + --output slice30.mha \ + --verbose +``` + +### 5. Visualize the Result + +You can visualize the result using a viewer (e.g., VV). The resulting image should look like this: + +![Varian](Varian.png){w=400px alt="Varian snapshot"} + +--- + +## Varian ProBeam Reconstruction + +Follow these steps for the Varian ProBeam format: + +### 1. Download Dataset + +Download the dataset from [Varian-ProBeam-data](https://data.kitware.com/api/v1/item/5be94bef8d777f2179a24ae1/download). + +### 2. Convert Geometry + +Run the application to convert Varian ProBeam's geometry into RTK's format: + +```bash +rtkvarianprobeamgeometry \ + --xml_file Scan.xml \ + --path Acquisitions/733061622 \ + --regexp Proj_.*.xim \ + -o geometry.xml +``` + +### 3. Reconstruct Using RTK Applications + +Reconstruct a slice (e.g., slice 58) of the volume using the `rtkfdk` algorithm: + +```bash +rtkfdk \ + --geometry geometry.xml \ + --regexp .*\.xim \ + --path Acquisitions/733061622 \ + --output slice58.mha \ + --verbose \ + --spacing 0.25,0.25,0.25 \ + --dimension 1024,1,1024 \ + --origin -127.875,-58,-127.875 +``` + +### 4. Apply the FOV Filter + +Apply the field-of-view (FOV) filter to discard everything outside the FOV: + +```bash +rtkfieldofview \ + --geometry geometry.xml \ + --regexp .*\.xim \ + --path Acquisitions/733061622 \ + --reconstruction slice58.mha \ + --output slice58.mha \ + --verbose +``` + +### 5. Visualize the Result + +You can visualize the result using a viewer (e.g., VV). The resulting image should look like this: + +![VarianProBeam](VarianProBeam.png){w=400px alt="VarianProBeam snapshot"} \ No newline at end of file diff --git a/examples/VarianReconstruction/Varian.png.sha512 b/examples/VarianReconstruction/Varian.png.sha512 new file mode 100644 index 000000000..7247aa401 --- /dev/null +++ b/examples/VarianReconstruction/Varian.png.sha512 @@ -0,0 +1 @@ +32682b605de454017def2fd481c13f8f8ad7bf989025f0534538655e39c5e82dda7c41c00e558cfd8399ed15cc899f0dff8f8720b5bffb6a426e5c199ca82d0d \ No newline at end of file diff --git a/examples/VarianReconstruction/VarianProBeam.png.sha512 b/examples/VarianReconstruction/VarianProBeam.png.sha512 new file mode 100644 index 000000000..a7b7d433e --- /dev/null +++ b/examples/VarianReconstruction/VarianProBeam.png.sha512 @@ -0,0 +1 @@ +8e0e9cabae8aa32fb54f23e26791133c40766fda48b94c9145a50a4110ae26291f7aca8f1d29699595cdb1199aacb990ddf12d52a674cbbb82e8ef4e043fd0c3 \ No newline at end of file diff --git a/examples/WaterPreCorrection/README.md b/examples/WaterPreCorrection/README.md new file mode 100644 index 000000000..40ca7f6bc --- /dev/null +++ b/examples/WaterPreCorrection/README.md @@ -0,0 +1,11 @@ + +# WaterPreCorrection + +This example illustrates how to apply empirical cupping correction using the [algorithm of Kachelriess et al.](http://onlinelibrary.wiley.com/doi/10.1118/1.2188076/abstract) named [WaterPrecorrection](http://www.openrtk.org/Doxygen/classrtk_1_1WaterPrecorrectionImageFilter.html) in RTK. The example uses a Gate simulation using the [fixed forced detection actor](http://wiki.opengatecollaboration.org/index.php/Users_Guide:Tools_to_Interact_with_the_Simulation_:_Actors#Fixed_Forced_Detection_CT). + +The simulation implements a 120 kV beam, a detector with 512x3 pixels and an energy response curve. Only the primary beam is simulated. + +This version uses the [Python wrapping](https://pypi.org/project/itk-rtk/). The simulation files, the output projections and the processing script are available [here](https://data.kitware.com/api/v1/file/5d394cea877dfcc9022c922b/download). + +```{literalinclude} waterPreCorrection.py +``` diff --git a/examples/index.md b/examples/index.md index b09db27a0..f8734241a 100644 --- a/examples/index.md +++ b/examples/index.md @@ -5,4 +5,16 @@ Examples :maxdepth: 1 ./FDK/README +./ConjugateGradient/README +./ForwardProjection/README +./RayBoxIntersection/README +./WaterPreCorrection/README +./CreateGammexPhantom/README +./AmsterdamShroud/README +./ElektaReconstruction/README +./VarianReconstruction/README +./MotionCompensationReconstruction/README +./TotalVariationRegularizedReconstruction/README +./DaubechiesWaveletsRegularizedReconstruction/README +./4DRooster/README ``` \ No newline at end of file diff --git a/index.md b/index.md index fdeedccf4..38306c413 100644 --- a/index.md +++ b/index.md @@ -6,30 +6,33 @@ The Reconstruction Toolkit (RTK) is an open-source and cross-platform software f [![PyPI](https://img.shields.io/pypi/v/itk-rtk.svg)](https://pypi.python.org/pypi/itk-rtk) [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://github.com/RTKConsortium/RTK/blob/master/LICENSE.TXT) +```{toctree} +GettingStarted +``` ```{toctree} :maxdepth: 1 :caption: 💾 Download - INSTALLATION ``` ```{toctree} :maxdepth: 1 :caption: 📖 Learn - -GettingStarted documentation/docs/Phantom.md +documentation/docs/Geometry.md +documentation/docs/Tutorial.md examples/index ``` ```{toctree} :maxdepth: 1 :caption: 🔨 Develop - API Discussion Issue tracker +CodeContribution documentation/docs/README documentation/docs/Release + ``` \ No newline at end of file diff --git a/test/Input/GeometricPhantom/SheppLogan_forbild.txt.md5 b/test/Input/GeometricPhantom/SheppLogan_forbild.txt.md5 new file mode 100644 index 000000000..8267e2433 --- /dev/null +++ b/test/Input/GeometricPhantom/SheppLogan_forbild.txt.md5 @@ -0,0 +1 @@ +674da1cb2da68b3050ac6a01 \ No newline at end of file