FanBeam

This example illustrates how to do a FBP fan-beam reconstruction using RTK. Since RTK only has 3D (back-projector), we create a 2D projection from the fan-beam projection. The pixels values of the first and the third rows are actually useless but required to let the backprojection work.

#!/usr/bin/env python

import SimpleRTK as srtk
import numpy as np

# Create a geometry
g = srtk.ThreeDCircularProjectionGeometry()
for i in range(720):
    g.AddProjection(1000,1500,i*.5)

# Simulate fan-beam acquisition
proj  = srtk.ConstantImageSource()
proj.SetSize([1024,1,720])
proj.SetSpacing([0.5]*3)
proj.SetOrigin((np.array(proj.GetSize())-1)*np.array(proj.GetSpacing())*-0.5)
proj = proj.Execute()
sl = srtk.SheppLoganPhantomFilter()
sl.SetGeometry(g)
proj = sl.Execute(proj)

writer = srtk.ImageFileWriter()
writer.SetFileName('fanbeam.mha')
writer.Execute(proj)

# RTK currently has 3D (back-)projector. To overcome this, we add two rows
projarray = srtk.GetArrayFromImage(proj)
projarray = np.concatenate([projarray,projarray,projarray], 1)
projarray[:,0,:] = 1e10 # Just to illustrate that this row values are not used
projarray[:,2,:] = 1e10 # Just to illustrate that this row values are not used
proj2D = srtk.GetImageFromArray(projarray)
proj2D.SetSpacing(proj.GetSpacing())
proj2D.SetOrigin((np.array(proj2D.GetSize())-1)*np.array(proj2D.GetSpacing())*-0.5)
writer.SetFileName('conebeam.mha')
writer.Execute(proj2D)

# Reconstruct using FDK
recon  = srtk.ConstantImageSource()
recon.SetSize([512,1,512])
recon.SetSpacing([0.5]*3)
recon.SetOrigin((np.array(recon.GetSize())-1)*np.array(recon.GetSpacing())*-0.5)
recon = recon.Execute()
fdk = srtk.FDKConeBeamReconstructionFilter()
fdk.SetGeometry(g)
recon = fdk.Execute(recon,proj2D)
writer.SetFileName('fdk.mha')
writer.Execute(recon)

This example now illustrates how to do a 2D conjugate gradient fan-beam reconstruction with RTK. Since the projector is ray-based, it's the opposite: the projection can be 1D but the volume must be 3D and cannot be 2D as in the previous example. The backprojector must then be the adjoint of the projector.

#!/usr/bin/env python
 
import SimpleRTK as srtk
import numpy as np
 
# Create a geometry
g = srtk.ThreeDCircularProjectionGeometry()
#g.SetRadiusCylindricalDetector(1500)
for i in range(720):
    g.AddProjection(1000,1500,i*.5)

# Simulate fan-beam acquisition
proj  = srtk.ConstantImageSource()
proj.SetSize([1024,1,720])
proj.SetSpacing([0.5]*3)
proj.SetOrigin((np.array(proj.GetSize())-1)*np.array(proj.GetSpacing())*-0.5)
proj = proj.Execute()
sl = srtk.SheppLoganPhantomFilter()
sl.SetGeometry(g)
proj = sl.Execute(proj)
 
writer = srtk.ImageFileWriter()
writer.SetFileName('fanbeam.mha')
writer.Execute(proj)

# The ray-based (back-)projector needs a 3D image even if only the central slice (in front of the 1D projection) is used
recon  = srtk.ConstantImageSource()
recon.SetSize([512,3,512])
recon.SetSpacing([0.5]*3)
recon.SetOrigin((np.array(recon.GetSize())-1)*np.array(recon.GetSpacing())*-0.5)
recon.SetConstant(0)

recon = recon.Execute()
mask = proj*0.+1.
cg = srtk.ConjugateGradientConeBeamReconstructionFilter()
cg.SetGeometry(g)
cg.SetForwardProjectionFilter(0)
cg.SetBackProjectionFilter(1)
cg.SetNumberOfIterations(50)
cg.SetDisableDisplacedDetectorFilter(True)
recon = cg.Execute(recon,proj,mask)
writer.SetFileName('cg.mha')
writer.Execute(recon)