ENH: Add Python examples and associated documentation

Co-authored-by: Cyril Mory <cyril.mory@gmail.com>

Author: Simon Rit

applications/README.md

@@ -20,6 +20,10 @@ The following are examples using RTK applications:
 ./rtkfourdconjugategradient/README.md
 ./rtkfourdsart/README.md
 ./rtkfourdrooster/README.md
+./rtkspectralforwardmodel/README.md
+./rtkspectralsimplexdecomposition/README.md
+./rtkspectralrooster/README.md
+./rtkspectralonestep/README.md
 ./rtkshowgeometry/README.md
 ./rtkosem/README.md
 ./rtksart/README.md

applications/rtkspectralforwardmodel/README.md

@@ -0,0 +1,31 @@
+# Two-step spectral CT reconstruction
+
+This page describes an application for two-step spectral CT reconstruction. Please start by reading the [general documentation on spectral CT](../../documentation/docs/multi_energy.md).
+
+# Spectral forward model
+
+RTK provides a command line application to simulate the measured projections in spectral CT from the decomposed projections, the incident spectrum, the material attenuations and the detector response.
+It is the "forward operator" in the [inverse problem](https://en.wikipedia.org/wiki/Inverse_problem) that RTK solves to obtain the decomposed projections from the measured projections.
+
+## Dataset
+
+The data used in this documentation page can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69a005bdf94d0dc14645a412).
+It is a simulation through a cylinder of water with two inserts, one empty, the other filled with cortical bone.
+
+## Command line application
+
+```
+# Simulate the measured projections (photon counts on the spectral detector)
+# from the decomposed projections, the incident spectrum,
+# the material attenuations and the detector response.
+
+rtkspectralforwardmodel \
+  -o cl_counts.mha \
+  -i decomposed.mha \
+  -d drm.mha \
+  --incident no_vector_spectrum.mha \
+  -a mat_basis.mha \
+  -t 20,40,60,80,100
+```
+
+`cl_counts.mha` should be identical to `counts.mha`. You can check it with [VV](http://vv.creatis.insa-lyon.fr/) using the overlay feature.

applications/rtkspectralonestep/README.md

@@ -0,0 +1,52 @@
+# One-step spectral CT reconstruction
+
+Please start by reading the [general documentation on spectral CT](../../documentation/docs/multi_energy.md).
+
+RTK provides a command line application to perform one-step reconstruction. It uses the measured projections, the incident spectrum, the material attenuations and the detector response for the spectral part, and the geometry and an initial reconstructed volume for the tomographic part.
+Providing a mask of the region of the reconstructed volume where the image object is expected to be found, and outside of which there should be nothing, helps the algorithm converge.
+The method is described in [this paper](https://ieeexplore.ieee.org/document/7979598).
+
+## Dataset
+
+The data used in this documentation page can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69a005bdf94d0dc14645a412).
+It is a simulation through a cylinder of water with two inserts, one empty, the other filled with cortical bone.
+
+## Command line application
+
+```
+# One-step spectral reconstruction
+
+rtkspectralonestep \
+  -s counts.mha \
+  --incident spectrum.mha \
+  -a mat_basis.mha \
+  -d drm.mha \
+  -g geometry.xml \
+  -i decomposed_init.mha \
+  -o cl_onestep.mha \
+  -t 20,40,60,80,100 \
+  --mask mask.mha \
+  -n 20
+```
+
+`cl_onestep.mha` should look similar to `recon_ref.mha`. You can check it with [VV](http://vv.creatis.insa-lyon.fr/) using the overlay feature. Note that `recon_ref.mha` is a 4D image, while `cl_onestep.mha` is a VectorImage, so their data is not organized the same way on disk, but VV will open both a 4D images and overlay them correctly.
+
+To speed things up, it is recommended to use the CUDA version of the forward and back projectors by running this command instead:
+
+```
+# One-step spectral reconstruction using CUDA accelerations
+
+rtkspectralonestep \
+  -s counts.mha \
+  --incident spectrum.mha \
+  -a mat_basis.mha \
+  -d drm.mha \
+  -g geometry.xml \
+  -i decomposed_init.mha \
+  -o cl_onestep.mha \
+  -t 20,40,60,80,100 \
+  --mask mask.mha \
+  -n 20 \
+  -f CudaRayCast \
+  -b CudaVoxelBased
+```

applications/rtkspectralrooster/README.md

@@ -0,0 +1,53 @@
+# Two-step spectral CT reconstruction
+
+This page describes an application for two-step spectral CT reconstruction. Please start by reading the [general documentation on spectral CT](../../documentation/docs/multi_energy.md).
+
+Once measured projections have been decomposed into material projections, any reconstruction algorithm can be used to reconstruct each material. In RTK, simply splitting the material projections file (saved as an `itkVectorImage`) into one set of projections per material, and then using `rtkfdk`, `rtksart` or `rtkconjugategradient` will work.
+
+# Spectral ROOSTER
+
+RTK provides a command line application to jointly reconstruct all materials at once, allowing to apply spatial regularization in each material, but also regularization along the material dimension, as reconstructed material volumes sometimes share the same spatial structure.
+Since `rtkspectralrooster` uses the [ROOSTER filter](https://www.openrtk.org/Doxygen/classrtk_1_1FourDROOSTERConeBeamReconstructionFilter.html), originally designed to reconstruct 3D+time volume sequences, the material dimension is referred to as time in the command-line help and parameter names.
+
+## Dataset
+
+The data used in this documentation page can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69a005bdf94d0dc14645a412).
+It is a simulation through a cylinder of water with two inserts, one empty, the other filled with cortical bone.
+
+## Command line application
+
+```
+# Reconstruct a set of material volumes from all material projections,
+# applying Total Nuclear Variation regularization
+
+rtkspectralrooster \
+  -o cl_spectralrooster.mha \
+  -i mat_recon.mha \
+  -p decomposed.mha \
+  -g geometry.xml \
+  --niter 4 \
+  --cgiter 2 \
+  --tviter 2 \
+  --gamma_tnv 1
+```
+
+`cl_spectralrooster.mha` should look similar to `recon_ref.mha`. You can check it with [VV](http://vv.creatis.insa-lyon.fr/) using the overlay feature. Note that `recon_ref.mha` is a 4D image, while `cl_spectralrooster.mha` is a VectorImage, so their data is not organized the same way on disk, but VV will open both a 4D images and overlay them correctly.
+To speed things up, it is recommended to use the CUDA version of the forward and back projectors, and by using CUDA in the conjugate gradient descent algorithm, by running this command instead:
+
+```
+# Reconstruct a set of material volumes from all material projections,
+# applying Total Nuclear Variation regularization, using CUDA accelerations
+
+rtkspectralrooster \
+  -o cl_spectralrooster.mha \
+  -i mat_recon.mha \
+  -p decomposed.mha \
+  -g geometry.xml \
+  --niter 4 \
+  --cgiter 2 \
+  --tviter 2 \
+  --gamma_tnv 1
+  --cudacg \
+  --fp CudaRayCast \
+  --bp CudaVoxelBased
+```

applications/rtkspectralsimplexdecomposition/README.md

@@ -0,0 +1,33 @@
+# Two-step spectral CT reconstruction
+
+This page describes an application for two-step spectral CT reconstruction. Please start by reading the [general documentation on spectral CT](../../documentation/docs/multi_energy.md).
+
+# Spectral simplex decomposition
+
+RTK provides a command line application to decompose the measured projections into one set of projections per material assumed compose the imaged object. It uses the incident spectrum, the material attenuations and the detector response.
+The method is described in [this paper](https://iopscience.iop.org/article/10.1088/0031-9155/53/15/002).
+
+## Dataset
+
+The data used in this documentation page can be downloaded using [Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69a005bdf94d0dc14645a412).
+It is a simulation through a cylinder of water with two inserts, one empty, the other filled with cortical bone.
+
+## Command line application
+
+```
+# Decompose the measured projections (photon counts on the spectral detector)
+# into material projections, using the incident spectrum,
+# the material attenuations and the detector response.
+
+rtkspectralsimplexdecomposition \
+  -o cl_decomposed.mha \
+  -i zero_decomposed.mha \
+  -s counts.mha \
+  -d drm.mha \
+  --incident no_vector_spectrum.mha \
+  -a mat_basis.mha \
+  -t 20,40,60,80,100
+
+```
+
+`cl_decomposed.mha` should be identical to `decomposed.mha`. You can check it with [VV](http://vv.creatis.insa-lyon.fr/) using the overlay feature.

examples/README.md

@@ -22,4 +22,6 @@ providing efficient implementations for high-performance applications.
 ./InlineReconstruction/README.md
 ./AddNoise/README.md
 ./GeometricPhantom/README.md
+./SpectralOneStep/README.md
+./SpectralTwoStep/README.md
 ```

examples/SpectralOneStep/README.md

@@ -0,0 +1,13 @@
+# One-step spectral reconstruction
+
+This example shows how to generate simulated spectral x-ray projections of an in silico phantom, and then perform a one-step reconstruction using [Mechlem *et al*'s reconstruction](https://doi.org/10.1109/TMI.2017.2726687) (ordered subset separable quadratic surrogate using [Weidinger *et al*'s surrogate](https://doi.org/10.1155/2016/5871604)).
+
+`````{tab-set}
+
+````{tab-item} Python
+
+```{literalinclude} ./SpectralOneStep.py
+:language: python
+```
+````
+`````