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/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/README.md b/examples/ConjugateGradient/README.md new file mode 100644 index 000000000..48ae52683 --- /dev/null +++ b/examples/ConjugateGradient/README.md @@ -0,0 +1,39 @@ +# Conjugate gradient + +--- + +This example uses the Shepp–Logan phantom. + +## 3D + +```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/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/README.md b/examples/ForwardProjection/README.md new file mode 100644 index 000000000..a578bc07e --- /dev/null +++ b/examples/ForwardProjection/README.md @@ -0,0 +1,10 @@ +# Forward Projection + +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/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..6e5c0dfc5 --- /dev/null +++ b/examples/RayBoxIntersection/README.md @@ -0,0 +1,6 @@ +# Ray box Intersection + +This script uses the file [Box.txt](http://wiki.openrtk.org/images/f/fa/Box.txt) as input. + +```{literalinclude} Code.sh +``` \ 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..ef36a829b 100644 --- a/examples/index.md +++ b/examples/index.md @@ -5,4 +5,8 @@ Examples :maxdepth: 1 ./FDK/README +./ConjugateGradient/README +./ForwardProjection/README +./RayBoxIntersection/README +./WaterPreCorrection/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