Conjugate gradient¶
3D
This script uses the SheppLogan phantom.
# 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 10
# Create a reference volume for comparison
rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha
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 (https://data.kitware.com/#item/67d1ff45c6dec2fc9c534d0c/download).
# Generate geometry for fan-beam setup
rtksimulatedgeometry -n 720 -o geometry.xml --arc 360
# 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
# Perform Conjugate Gradient reconstruction
rtkconjugategradient -p . -r projections.mha -o cg.mha -g geometry.xml --spacing 2 --dimension 256 1 256 -n 10
# Create a reference volume for comparison
rtkdrawgeometricphantom --spacing 2 --dimension=256,1,256 --phantomfile SheppLogan-2d.txt -o ref.mha --phantomscale=256,1,256
Noisy Reconstruction
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 and the associated weights, as well as the geometry, and run the following to compare the regular least squares reconstruction (without weights) and the weighted least squares reconstruction.
# 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
# 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 10
# Create a reference volume for comparison
rtkdrawshepploganphantom --spacing 2 --dimension 256 -o ref.mha