From 274242cf51aa98bb0b96fe82731040d5c74779c8 Mon Sep 17 00:00:00 2001 From: axel garcia Date: Fri, 27 Dec 2024 18:13:24 +0100 Subject: [PATCH] DOC: Update new ReadTheDocs with content from the old documentation --- CodeContribution.md | 25 +++ GettingStarted.md | 2 +- INSTALLATION.md | 2 +- README.md | 2 +- conf.py | 2 +- documentation/docs/CMakeLists.txt | 1 + documentation/docs/Tutorial.md | 21 ++ examples/4DRooster/Blurred.jpg.sha512 | 1 + examples/4DRooster/MotionMask.jpg.sha512 | 1 + examples/4DRooster/README.md | 211 ++++++++++++++++++ examples/4DRooster/Signal.jpg.sha512 | 1 + examples/AmsterdamShroud/Amsterdam.png.sha512 | 1 + .../Moving-Phantom-Sinogram.png.sha512 | 1 + examples/AmsterdamShroud/README.md | 13 ++ .../ConjugateGradient-Sinogram.png.sha512 | 1 + .../ConjugateGradient.png.sha512 | 1 + .../ConjugateGradient/ConjugateGradient2D.sh | 11 + .../ConjugateGradient/ConjugateGradient3D.sh | 21 ++ examples/ConjugateGradient/README.md | 49 ++++ .../GammexPhantom.png.sha512 | 1 + examples/CreateGammexPhantom/README.md | 10 + .../DaubechiesWavelets.sh | 12 + .../Overlay.png.sha512 | 1 + .../README.md | 9 + .../Sinogram.png.sha512 | 1 + .../ElektaReconstruction/Elekta.png.sha512 | 1 + examples/ElektaReconstruction/README.md | 72 ++++++ examples/FDK/{Code2D.sh => FDK2D.sh} | 4 +- examples/FDK/{Code3D.sh => FDK3D.sh} | 4 +- examples/FDK/README.md | 4 +- .../ForwardProjection/ForwardProjection.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 | 162 ++++++++++++++ .../Signal.jpg.sha512 | 1 + .../VectorField.gif.sha512 | 1 + examples/RayBoxIntersection/README.md | 9 + .../RayBox-Sinogram.png.sha512 | 1 + examples/RayBoxIntersection/RayBox.png.sha512 | 1 + .../RayBoxIntersection/RayBoxInteraction.sh | 11 + .../Overlay.png.sha512 | 1 + .../README.md | 9 + .../Sinogram.png.sha512 | 1 + ...TotalVariationRegularizedReconstruction.sh | 12 + examples/VarianReconstruction/README.md | 115 ++++++++++ .../VarianReconstruction/Varian.png.sha512 | 1 + .../VarianProBeam.png.sha512 | 1 + .../WaterPreCorrection/Profile.png.sha512 | 1 + examples/WaterPreCorrection/README.md | 13 ++ .../WaterPreCorrection/WaterPreCorrection.py | 82 +++++++ examples/index.md | 18 +- index.md | 14 +- .../SheppLogan_forbild.txt.md5 | 1 + 57 files changed, 949 insertions(+), 18 deletions(-) create mode 100644 CodeContribution.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 100644 examples/ConjugateGradient/ConjugateGradient-Sinogram.png.sha512 create mode 100644 examples/ConjugateGradient/ConjugateGradient.png.sha512 create mode 100644 examples/ConjugateGradient/ConjugateGradient2D.sh create mode 100755 examples/ConjugateGradient/ConjugateGradient3D.sh 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/DaubechiesWavelets.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 rename examples/FDK/{Code2D.sh => FDK2D.sh} (57%) rename examples/FDK/{Code3D.sh => FDK3D.sh} (62%) create mode 100644 examples/ForwardProjection/ForwardProjection.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/README.md create mode 100644 examples/RayBoxIntersection/RayBox-Sinogram.png.sha512 create mode 100644 examples/RayBoxIntersection/RayBox.png.sha512 create mode 100644 examples/RayBoxIntersection/RayBoxInteraction.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/TotalVariationRegularizedReconstruction/TotalVariationRegularizedReconstruction.sh 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/Profile.png.sha512 create mode 100644 examples/WaterPreCorrection/README.md create mode 100644 examples/WaterPreCorrection/WaterPreCorrection.py 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/GettingStarted.md b/GettingStarted.md index d62b648ae..ab62e7b8a 100644 --- a/GettingStarted.md +++ b/GettingStarted.md @@ -7,6 +7,6 @@ Here are suggested steps for the RTK beginner. 2. If you are not already familiar with ITK, [get started with ITK](https://github.com/InsightSoftwareConsortium/ITK/blob/master/GettingStarted.md). - 3. Check out the [FirstReconstruction](https://github.com/SimonRit/RTK/blob/master/examples/FirstReconstruction) example (both the [CMake](https://github.com/SimonRit/RTK/blob/master/examples/FirstReconstruction/CMakeLists.txt) and the [C++](https://github.com/SimonRit/RTK/blob/master/examples/FirstReconstruction/FirstReconstruction.cxx) codes). Many other examples are on the wiki using the applications (corresponding source code in the [applications](https://github.com/SimonRit/RTK/blob/master/applications) subdirectory). + 3. Check out the [FirstReconstruction](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction) example (both the [CMake](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction/CMakeLists.txt) and the [C++](https://github.com/RTKConsortium/RTK/blob/master/examples/FirstReconstruction/FirstReconstruction.cxx) codes). Many other examples are in the [documentation](https://docs.openrtk.org/) using the applications (corresponding source code in the [applications subdirectory](https://github.com/RTKConsortium/RTK/blob/master/applications)). 4. Ask question to the [user mailing list](https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users). diff --git a/INSTALLATION.md b/INSTALLATION.md index 88083782d..dd4fabd81 100644 --- a/INSTALLATION.md +++ b/INSTALLATION.md @@ -7,7 +7,7 @@ RTK is a module of [ITK](https://www.itk.org), the Insight Toolkit. Follow the i * `Module_RTK`: Activates RTK download and compilation. Default is `OFF`. Turn it `ON` to activate RTK or compile RTK independently (see below). * `Module_RTK_GIT_TAG`: Git tag for the RTK download. By default, the RTK version which is downloaded and compiled is the one given in the [RTK.remote.cmake](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Remote/RTK.remote.cmake). Change this option to build another version. For example, you can change it to `master` to build the latest RTK version. RTK is only maintained to be backward compatible with the latest ITK release and ITK master branch. -* `RTK_BUILD_APPLICATIONS`: Activates the compilation of RTK's command line tools. Although RTK is mainly a toolkit, we also provide several command line tools for doing most of the available processing. These command line tools use [gengetopt](https://www.gnu.org/software/gengetopt/gengetopt.html). Several examples are available on the [Applications](https://wiki.openrtk.org/index.php/RTK_wiki_help#Applications) section of the [wiki](http://wikiopenrtk.org). +* `RTK_BUILD_APPLICATIONS`: Activates the compilation of RTK's command line tools. Although RTK is mainly a toolkit, we also provide several command line tools for doing most of the available processing. These command line tools use [gengetopt](https://www.gnu.org/software/gengetopt/gengetopt.html). Several examples are available in the [documentation](http://docs.openrtk.org). * `RTK_USE_CUDA`: Activates CUDA computation. Default is `ON` if CMake has automatically found the CUDA package and a CUDA-compatible GPU, and `OFF` otherwise. * `RTK_CUDA_VERSION`: Specifies an exact version of the CUDA toolkit which must be used. If unspecified, RTK only checks if the found version is recent enough. * `RTK_CUDA_PROJECTIONS_SLAB_SIZE`: Set the number of projections processed at once in CUDA processing. Default is 16. diff --git a/README.md b/README.md index 61d1353c2..69caef0ba 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ Links * [Download](https://www.openrtk.org/RTK/resources/software.html) * [Mailing List](https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users) * [Getting Started](GettingStarted.md) -* [Help](https://wiki.openrtk.org) +* [Help](https://docs.openrtk.org) * [Issue tracking](https://github.com/RTKConsortium/RTK/issues) diff --git a/conf.py b/conf.py index d202dfdd9..fe406c95c 100644 --- a/conf.py +++ b/conf.py @@ -15,7 +15,7 @@ def setup(app): # -- Project information ----------------------------------------------------- project = 'RTK' -copyright = f'{date.today().year}, Simon Rit' +copyright = f'{date.today().year}, RTK Consortium' author = 'RTK Consortium' # The full version, including alpha/beta/rc tags 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/Tutorial.md b/documentation/docs/Tutorial.md new file mode 100644 index 000000000..1baa62cad --- /dev/null +++ b/documentation/docs/Tutorial.md @@ -0,0 +1,21 @@ +# 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..a87203be8 --- /dev/null +++ b/examples/4DRooster/README.md @@ -0,0 +1,211 @@ +# 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. + +## 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/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/ConjugateGradient2D.sh b/examples/ConjugateGradient/ConjugateGradient2D.sh new file mode 100644 index 000000000..45d2b5387 --- /dev/null +++ b/examples/ConjugateGradient/ConjugateGradient2D.sh @@ -0,0 +1,11 @@ +# Generate geometry for fan-beam setup +rtksimulatedgeometry -n 720 -o geometry.xml --arc 360 + +# Simulate projections using the Shepp-Logan phantom +rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 0.5 --dimension 1024,3,720 + +# Create an initial 3D volume for reconstruction +rtkdrawshepploganphantom --spacing 0.5 --dimension 512,1,512 -o initial_volume.mha --phantomscale=512,1,512 + +# Perform Conjugate Gradient reconstruction +rtkconjugategradient -p . -r projections.mha -o cg.mha -g geometry.xml --spacing 0.5 --dimension 512 3 512 -n 50 diff --git a/examples/ConjugateGradient/ConjugateGradient3D.sh b/examples/ConjugateGradient/ConjugateGradient3D.sh new file mode 100755 index 000000000..dd600869e --- /dev/null +++ b/examples/ConjugateGradient/ConjugateGradient3D.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..bbe7f0526 --- /dev/null +++ b/examples/ConjugateGradient/README.md @@ -0,0 +1,49 @@ +# 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 +``` + +## 2D + +The same reconstruction can be performed using the original 2D Shepp-Logan phantom. +RTK can perform 2D reconstructions through images wide of 1 pixel in the y direction. +The following script performs the same reconstruction as above in a 2D environment and use the Shepp–Logan phantom. + +```{literalinclude} ConjugateGradient2D.sh +``` 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/DaubechiesWavelets.sh b/examples/DaubechiesWaveletsRegularizedReconstruction/DaubechiesWavelets.sh new file mode 100644 index 000000000..5b4a2361a --- /dev/null +++ b/examples/DaubechiesWaveletsRegularizedReconstruction/DaubechiesWavelets.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..fc19a0c41 --- /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} DaubechiesWavelets.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/FDK/Code2D.sh b/examples/FDK/FDK2D.sh similarity index 57% rename from examples/FDK/Code2D.sh rename to examples/FDK/FDK2D.sh index d02b8a858..c7cc0cac0 100644 --- a/examples/FDK/Code2D.sh +++ b/examples/FDK/FDK2D.sh @@ -3,10 +3,10 @@ rtksimulatedgeometry -n 180 -o geometry.xml # Create projections of the phantom file # Note the sinogram being 3 pixels wide in the y direction to allow back-projection interpolation in a 2D image -rtkprojectgeometricphantom -g geometry.xml -o projections.mha --spacing 2 --dimension=512,3,512 --phantomfile SheppLogan-2d.txt --phantomscale=256,1,256 +rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension=512,3,512 --phantomscale=256,1,256 # Reconstruct rtkfdk -p . -r projections.mha -o fdk.mha -g geometry.xml --spacing 2 --dimension 256,1,256 # Create a reference volume for comparison -rtkdrawgeometricphantom --spacing 2 --dimension=256,1,256 --phantomfile SheppLogan-2d.txt -o ref.mha --phantomscale=256,1,256 +rtkdrawshepploganphantom --spacing 2 --dimension=256,1,256 -o ref.mha --phantomscale=256,1,256 diff --git a/examples/FDK/Code3D.sh b/examples/FDK/FDK3D.sh similarity index 62% rename from examples/FDK/Code3D.sh rename to examples/FDK/FDK3D.sh index c0176aa3a..0cb49656c 100644 --- a/examples/FDK/Code3D.sh +++ b/examples/FDK/FDK3D.sh @@ -3,10 +3,10 @@ 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 -rtkprojectgeometricphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 --phantomfile SheppLogan.txt +rtkprojectshepploganphantom -g geometry.xml -o projections.mha --spacing 2 --dimension 256 # Reconstruct rtkfdk -p . -r projections.mha -o fdk.mha -g geometry.xml --spacing 2 --dimension 256 # Create a reference volume for comparison -rtkdrawgeometricphantom --spacing 2 --dimension 256 --phantomfile SheppLogan.txt -o ref.mha +rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha diff --git a/examples/FDK/README.md b/examples/FDK/README.md index d8f16f92a..a487d4284 100644 --- a/examples/FDK/README.md +++ b/examples/FDK/README.md @@ -9,7 +9,7 @@ Reconstruction of the Shepp–Logan phantom using Feldkamp, David and Kress cone This script uses the file [SheppLogan.txt](https://data.kitware.com/api/v1/item/5b179c938d777f15ebe2020b/download) as input. -```{literalinclude} Code3D.sh +```{literalinclude} FDK3D.sh ``` ## 2D @@ -21,5 +21,5 @@ The same reconstruction can be performed using the original 2D Shepp-Logan phant RTK can perform 2D reconstructions through images wide of 1 pixel in the y direction. The following script performs the same reconstruction as above in a 2D environment and uses the [2D Shepp-Logan](http://wiki.openrtk.org/images/7/73/SheppLogan-2d.txt) phantom as input. -```{literalinclude} Code2D.sh +```{literalinclude} FDK2D.sh ``` diff --git a/examples/ForwardProjection/ForwardProjection.sh b/examples/ForwardProjection/ForwardProjection.sh new file mode 100644 index 000000000..3fd607017 --- /dev/null +++ b/examples/ForwardProjection/ForwardProjection.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 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..6eab7f39e --- /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} ForwardProjection.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..c2e4ec4f2 --- /dev/null +++ b/examples/MotionCompensationReconstruction/README.md @@ -0,0 +1,162 @@ +# 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). + +## 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/README.md b/examples/RayBoxIntersection/README.md new file mode 100644 index 000000000..c9b7023bd --- /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} RayBoxInteraction.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/RayBoxIntersection/RayBoxInteraction.sh b/examples/RayBoxIntersection/RayBoxInteraction.sh new file mode 100644 index 000000000..5107895e6 --- /dev/null +++ b/examples/RayBoxIntersection/RayBoxInteraction.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/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..bc4e6660b --- /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} TotalVariationRegularizedReconstruction.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/TotalVariationRegularizedReconstruction/TotalVariationRegularizedReconstruction.sh b/examples/TotalVariationRegularizedReconstruction/TotalVariationRegularizedReconstruction.sh new file mode 100644 index 000000000..dc7d15c83 --- /dev/null +++ b/examples/TotalVariationRegularizedReconstruction/TotalVariationRegularizedReconstruction.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/VarianReconstruction/README.md b/examples/VarianReconstruction/README.md new file mode 100644 index 000000000..55bcbdf2f --- /dev/null +++ b/examples/VarianReconstruction/README.md @@ -0,0 +1,115 @@ +# 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/Profile.png.sha512 b/examples/WaterPreCorrection/Profile.png.sha512 new file mode 100644 index 000000000..c542ff6fa --- /dev/null +++ b/examples/WaterPreCorrection/Profile.png.sha512 @@ -0,0 +1 @@ +0722bc916860fca6a98ebaff9bb064d9153336f7fc18746329c557514f88450242361476075b1df296ac870ea966df7db825665fb8dc51f629650e0013d0b027 \ No newline at end of file diff --git a/examples/WaterPreCorrection/README.md b/examples/WaterPreCorrection/README.md new file mode 100644 index 000000000..9aeb92309 --- /dev/null +++ b/examples/WaterPreCorrection/README.md @@ -0,0 +1,13 @@ +# WaterPreCorrection + +This example illustrates how to apply empirical cupping correction using the [algorithm of Kachelriess et al.](https://www.doi.org/doi/10.1118/1.2188076/abstract) named [WaterPrecorrection](https://www.openrtk.org/Doxygen/classrtk_1_1WaterPrecorrectionImageFilter.html) in RTK. The example uses a Gate simulation using the [fixed forced detection actor](https://opengate.readthedocs.io/en/latest/tools_to_interact_with_the_simulation_actors.html#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. + +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 +``` +The resulting central profiles are + +![Profile](Profile.png){w=800px alt="Profile"} diff --git a/examples/WaterPreCorrection/WaterPreCorrection.py b/examples/WaterPreCorrection/WaterPreCorrection.py new file mode 100644 index 000000000..c49dfe256 --- /dev/null +++ b/examples/WaterPreCorrection/WaterPreCorrection.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python + +import itk +from itk import RTK as rtk +import os +import matplotlib.pyplot as plt +import numpy as np +import scipy.optimize +import glob + +# Script parameters +directory = "output/" +spacing = [ 0.5 ] * 3 +size = [ 512, 1, 512 ] +order = 6 +reference_attenuation = 0.02 +origin = [ -(size[0]/2+0.5)*spacing[0], 0., -(size[2]/2+0.5)*spacing[2]] + +# List of filenames +file_names = list() +for file in os.listdir(directory): + if file.startswith("attenuation") and file.endswith(".mha"): + file_names.append(directory + file) + +# Read in full geometry +geometry_reader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New() +geometry_reader.SetFilename ( "output/geometry.xml" ) +geometry_reader.GenerateOutputInformation() +geometry = geometry_reader.GetGeometry() + +# Crate template image +ImageType = itk.Image[itk.F,3] +projections_reader = rtk.ProjectionsReader[ImageType].New(file_names=file_names) +projections = projections_reader.GetOutput() +constant_image_filter = rtk.ConstantImageSource[ImageType].New(origin=origin, spacing=spacing, size=size) +constant_image = constant_image_filter.GetOutput() +template = rtk.draw_ellipsoid_image_filter(constant_image, density=reference_attenuation, axis=[100, 0, 100]) +itk.imwrite(template, 'template.mha') +template = itk.array_from_image(template).flatten() + +# Create weights (where the template should match) +weights = rtk.draw_ellipsoid_image_filter(constant_image, density=1., axis=[125, 0, 125]) +weights = rtk.draw_ellipsoid_image_filter(weights, density=-1., axis=[102, 0, 102]) +weights = rtk.draw_ellipsoid_image_filter(weights, density=1., axis=[98, 0, 98]) +itk.imwrite(weights, 'weights.mha') +weights = itk.array_from_image(weights).flatten() + +# Create reconstructed images +wpcoeffs = np.zeros((order+1)) +fdks = [None] * (order+1) +for o in range(0,order+1): + wpcoeffs[o-1] = 0. + wpcoeffs[o] = 1. + water_precorrection = rtk.water_precorrection_image_filter(projections, coefficients=wpcoeffs, in_place=False) + fdk = rtk.fdk_cone_beam_reconstruction_filter(constant_image, water_precorrection, geometry=geometry) + itk.imwrite(fdk, f'fdk{o}.mha') + fdks[o] = itk.array_from_image(fdk).flatten() + +# Create and solve the linear system of equation B.c= a to find the coeffs c +a = np.zeros((order+1)) +B = np.zeros((order+1,order+1)) +for i in range(0,order+1): + a[i] = np.sum(weights * fdks[i] * template) + for j in np.arange(i,order+1): + B[i,j] = np.sum(weights * fdks[i] * fdks[j]) + B[j,i] = B[i,j] +c = np.linalg.solve(B,a) + +water_precorrection = rtk.water_precorrection_image_filter(projections, coefficients=c) +fdk = rtk.fdk_cone_beam_reconstruction_filter(constant_image, water_precorrection, geometry=geometry) +itk.imwrite(fdk, 'fdk.mha' ) + +fdk = itk.imread('fdk.mha') +fdk1 = itk.imread('fdk1.mha') + +plt.plot(itk.array_from_image(fdk)[:,0,256], label='Corrected') +plt.plot(itk.array_from_image(fdk1)[:,0,256], label='Uncorrected') +plt.legend() +plt.xlabel('Pixel number') +plt.ylabel('Attenuation') +plt.xlim(0,512) +plt.savefig('profile.png') diff --git a/examples/index.md b/examples/index.md index b09db27a0..0e06d8da5 100644 --- a/examples/index.md +++ b/examples/index.md @@ -1,8 +1,22 @@ -Examples +Command-line examples ======== +RTK provides command line applications that can be built from the C++ code by turning `ON` the `RTK_BUILD_APPLICATIONS` CMake option. A few of these applications have also been translated to Python and integrated in the [Pypi package](https://pypi.org/project/itk-rtk/). The options of each command line application can be listed with the `--help option`. + ```{toctree} :maxdepth: 1 ./FDK/README -``` \ No newline at end of file +./ConjugateGradient/README +./ForwardProjection/README +./RayBoxIntersection/README +./CreateGammexPhantom/README +./AmsterdamShroud/README +./ElektaReconstruction/README +./VarianReconstruction/README +./MotionCompensationReconstruction/README +./TotalVariationRegularizedReconstruction/README +./DaubechiesWaveletsRegularizedReconstruction/README +./4DRooster/README +./WaterPreCorrection/README +``` diff --git a/index.md b/index.md index fdeedccf4..668bd9d5b 100644 --- a/index.md +++ b/index.md @@ -6,30 +6,32 @@ 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 -examples/index +documentation/docs/Tutorial.md +examples/index.md ``` ```{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