RTK/Meetings/TrainingNov17
Geometry
Real
Edit and analyze the geometry file obtained from a real Elekta Synergy acquisition. Compare to the projection files using the command
vv --sequence *his
(VV help: F1).
Simulated
Use the command line application rtksimulatedgeometry to generate simulated geometries and analyze the outputs. (Help: rtksimulatedgeometry -h).
Basic operators
Here is a list of the RTK command lines you will have to run during the training, along with the data you will need. To start, download the following CT data and a relevant Cone-Beam CT geometry file.
Forward projection
To forward project the CT volume, following the geometry described in geometry.xml, run
rtkforwardprojections -i muCT.mha -o projections.mha -g geometry.xml --verbose --dimension 128 --spacing 2
Look at the result using vv
vv projections.mha
There are as many projections as described in the geometry.xml file. You can edit the geometry.xml file, remove all projections except one, and try again. Such a modified geometry file is available here. Run the following command line
rtkforwardprojections -i muCT.mha -o singleProjection.mha -g singleProjectionGeometry.xml --verbose --dimension 256 --spacing 1
You can check with vv that the file contains only one projection
vv singleProjection.mha
Back projection
Start by back projecting a single projection. Run
rtkbackprojections -v -p . -r singleProjection.mha -o singleBackProjection.mha -g singleProjectionGeometry.xml --dimension 128 --spacing 2
Observe the result with vv
vv singleBackProjection.mha
As you can see, the projection data seems to be smeared into the volume. Now back project all projections, by running
rtkbackprojections -v -p . -r projections.mha -o backProjection.mha -g geometry.xml --dimension 128 --spacing 2
Look at the result with vv. In order to compare it with the original image, run
vv backProjection.mha muCT.mha --linkall
Click in one of the display windows, hit 'Ctrl + G' to center the view on the tumor. Hit 'w' to automatically compute reasonable window and level values. You can navigate between "muCT.mha" and "backProjection.mha" by hitting the 'Tab' key.
Using different forward and back projection algorithms
Both rtkforwardprojections and rtkbackprojections can use several algorithms to perform the forward or back projections. For forward projections, you can try the following command lines
rtkforwardprojections -i muCT.mha -o josephProjections.mha -g geometry.xml --fp Joseph --dimension 128 --spacing 2
rtkforwardprojections -i muCT.mha -o rayCastInterpolatorProjections.mha -g geometry.xml --fp RayCastInterpolator --dimension 128 --spacing 2
And for back projection
rtkbackprojections -p . -r projections.mha -o voxelBasedBackProjection.mha -g geometry.xml --bp VoxelBasedBackProjection --dimension 128 --spacing 2
rtkbackprojections -p . -r projections.mha -o josephBackProjection.mha -g geometry.xml --bp Joseph --dimension 128 --spacing 2
rtkbackprojections -p . -r projections.mha -o normalizedJosephBackProjection.mha -g geometry.xml --bp NormalizedJoseph --dimension 128 --spacing 2
If you have a CUDA-capable graphics card, and have compiled RTK with RTK_USE_CUDA=ON, you can also try the following
rtkforwardprojections -i muCT.mha -o cudaRayCastProjections.mha -g geometry.xml --fp CudaRayCast
rtkbackprojections -p . -r projections.mha -o cudaRayCastBackProjection.mha -g geometry.xml --bp CudaRayCast
rtkbackprojections -p . -r projections.mha -o cudaBackProjection.mha -g geometry.xml --bp CudaBackProjection
Pre-processing
Pre-process the sequence of real .his projections
- Crop
- Bin
- Median
3D FBP reconstruction
Real
Reconstruct the sequence of (pre-processed) real projections with and without width-truncated correction.
rtkfdk -p projections -r '.*.his' -o fdk.mha -g geometry.xml --dimension 128 --spacing 2 --binning 4,4
rtkfdk -p projections -r '.*.his' -o fdk_pad.mha -g geometry.xml --dimension 128 --spacing 2 --binning 4,4 --pad 0.5
Simulated
Simulate Shepp Logan projections (rtkprojectshepplogan) for a short scan and reconstruct them:
rtksimulatedgeometry -a 220 -n 110 -o parker.xml
rtkprojectshepploganphantom -g parker.xml --dimension 128 --spacing 2 -o parker.mha --phantomscale 64
rtkfdk -p . -r parker.mha -o fdk_park.mha -g parker.xml --dimension 128
Idem for width-truncated projections:
rtksimulatedgeometry -n 110 -o widthtrunc.xml --proj_iso_x 80
rtkprojectshepploganphantom -g widthtrunc.xml --dimension 128 --spacing 2 -o widthtrunc.mha --phantomscale 64
rtkfdk -p . -r widthtrunc.mha -o fdk_wt.mha -g widthtrunc.xml --dimension 128
3D iterative reconstruction
Creating a very small dataset
Various reconstruction algorithms exist, which employ forward and back projection extensively. We will now start performing real reconstructions using some of the algorithms implemented in RTK. Since these algorithms are computationally much more demanding than the standard FDK, we propose to run them with low resolution data. You will be able to re-compute them at full resolution later, since the present page will remain available to you after the training. Note that, for the regularized reconstruction techniques, you will have to re-tune the regularization parameters yourself, since they depend (among many things) on resolution.
Convert and preprocess the real projections by running
rtkprojections -p projections/ -r .*.his -o downsampledProjections.mha --binning "4, 4"
The '--binning "4, 4"' option downsamples the projections by 4 along both axes. Since this is not enough to reach acceptable reconstruction times for the training, we will build a sub-dataset by selecting only one projection out of 5. You can do this by running
rtksubselect -p . -r downsampledProjections.mha -g geometry.xml --out_geometry subGeometry.xml --out_proj subProjections.mha -s 5
This creates a new geometry file and a new projections file, from which even a laptop should be able to perform iterative reconstruction in reasonable time.
Unregularized 3D iterative reconstruction algorithms
The most popular iterative reconstruction algorithm is the Simultaneous Algebraic Reconstruction Technique (SART). Try it by running
rtksart -p . -r subProjections.mha -g subGeometry.xml -o SART.mha --dimension 128 --spacing 3
Another popular technique uses the conjugate gradient algorithm
rtkconjugategradient -p . -r subProjections.mha -g subGeometry.xml -o CG.mha --dimension 128 --spacing 3 -n 10
You can compare the results by running
vv SART.mha CG.mha --linkall
The SART reconstruction should be sharper than the conjugate gradient one. If you are interested in understanding the reasons for this difference, try running
rtksart -p . -r subProjections.mha -g subGeometry.xml -o OSSART_10projsPerSubset.mha --dimension 128 --spacing 3 --nprojpersubset 10
Visualization is improved when one masks the field-of-view of the reconstruction:
rtkfieldofview -p . -r subProjections.mha -g subGeometry.xml -m -o FOVmask.mha --reconstruction SART.mha
rtkfieldofview -p . -r subProjections.mha -g subGeometry.xml -o SART_masked.mha --reconstruction SART.mha
The differences in image quality (sharpness, contrast, ...) and convergence speed, stability, ... between the various algorithms may be discussed briefly during the training.
Regularized 3D iterative reconstruction algorithms
RTK also includes algorithms that perform regularization during the reconstruction. You can try for example
rtkadmmtotalvariation -p . -r subProjections.mha -g subGeometry.xml -o admmTV.mha --dimension 128 --spacing 3 -n 10 --beta 1000 --alpha 1
and
rtkadmmwavelets -p . -r subProjections.mha -g subGeometry.xml -o admmWavelets.mha --dimension 128 --spacing 3 -n 10 --beta 1000 --alpha 1
Note that these algorithms have little interest when working with such low resolution data, since most of the noise and artifacts they are meant to remove have already been eliminated by undersampling.
Dual energy CT simulation and decomposition
The following script generates noiseless line integrals through a phantom composed of a large cylinder of water and a smaller cylinder of iodine, with a concentration of 0.001 g/ml. It requires CLITK, so if you do not have a working CLITK, you can just use the files provided in the data folder.
# Generate geometry for a single projection
rtksimulatedgeometry -n 1 -o singleProjGeo.xml
# Compute analytical projections of the PRH phantom
for i in water iodine
do
rtkdrawgeometricphantom --phantom $i.txt -o vol_$i.mha --dimension 256,64,256 --spacing 1,1,1
rtkforwardprojections -o proj_$i.mha -i vol_$i.mha -g singleProjGeo.xml --dimension 256,32,1 --spacing 1,1,1 --fp Joseph
clitkSum -i proj_$i.mha -o proj_$i.mha -d 2
done
# Merge the projections into a single file
clitkMergeSequence proj_water.mha proj_iodine.mha -o lineIntegrals.mha
clitk4DImageToNVectorImage -i lineIntegrals.mha -o lineIntegrals.mha
# Create a starting point for decompotion
clitkImageCreate -o zero_proj.mha --like proj_water.mha
clitkMergeSequence zero_proj.mha zero_proj.mha -o zero_proj.mha
clitk4DImageToNVectorImage -i zero_proj.mha -o initial.mha
This has generated the files lineIntegrals.mha and initial.mha (the starting point for decomposition, set to 0 everywhere).
Forward model
This will compute the expected energy in each of the two bins of the detector, for each pixel of the projection. It is the simplest possible simulation tool to compute dual energy projections.
# Compute dual energy projections
rtkdualenergyforwardmodel -i lineIntegrals.mha -o dualEnergyProjections.mha --high highEnergyEffectiveSpectrum.mha --low lowEnergyEffectiveSpectrum.mha -a water_iodine_attenuation.mha
Decomposition
And the following command computes the line integrals from the dual energy projections. It is this step that is performed on actual dual energy CT projection data extracted from a real scanner.
# Recover line integrals
rtkdualenergysimplexdecomposition -i initial.mha -o decomposed.mha --high highEnergyEffectiveSpectrum.mha --low lowEnergyEffectiveSpectrum.mha -a water_iodine_attenuation.mha --dual dualEnergyProjections.mha -n 10000 -r
All this is done with noiseless data, but one can simulate any kind of noise on the output of the forward model, and try to perform the decomposition.