FanBeam: Difference between revisions

From Openrtk
Jump to navigation Jump to search
No edit summary
No edit summary
Line 47: Line 47:
recon = fdk.Execute(recon,proj2D)
recon = fdk.Execute(recon,proj2D)
writer.SetFileName('fdk.mha')
writer.SetFileName('fdk.mha')
writer.Execute(recon)
</source>
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.
<source lang="python">
#!/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)
# The projection can be 1D
projarray = srtk.GetArrayFromImage(proj)
projarray = projarray.reshape([projarray.shape[0],1,projarray.shape[1]])
proj2D = srtk.GetImageFromArray(projarray)
proj2D.SetSpacing([411./projarray.shape[2]]*3)
proj2D.SetOrigin((np.array(proj2D.GetSize())-1)*np.array(proj2D.GetSpacing())*-0.5)
writer = srtk.ImageFileWriter()
writer.SetFileName('conebeam.mha')
writer.Execute(proj2D)
# 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 = proj2D*0.+1.
cg = srtk.ConjugateGradientConeBeamReconstructionFilter()
cg.SetGeometry(g)
cg.SetForwardProjectionFilter(0)
cg.SetBackProjectionFilter(1)
cg.SetNumberOfIterations(50)
recon = cg.Execute(recon,proj2D,mask)
writer.SetFileName('cg.mha')
writer.Execute(recon)
writer.Execute(recon)
</source>
</source>

Revision as of 04:59, 18 May 2017

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()
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 projection can be 1D
projarray = srtk.GetArrayFromImage(proj)
projarray = projarray.reshape([projarray.shape[0],1,projarray.shape[1]])
proj2D = srtk.GetImageFromArray(projarray)
proj2D.SetSpacing([411./projarray.shape[2]]*3)
proj2D.SetOrigin((np.array(proj2D.GetSize())-1)*np.array(proj2D.GetSpacing())*-0.5)
writer = srtk.ImageFileWriter()
writer.SetFileName('conebeam.mha')
writer.Execute(proj2D)

# 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 = proj2D*0.+1.
cg = srtk.ConjugateGradientConeBeamReconstructionFilter()
cg.SetGeometry(g)
cg.SetForwardProjectionFilter(0)
cg.SetBackProjectionFilter(1)
cg.SetNumberOfIterations(50)
recon = cg.Execute(recon,proj2D,mask)
writer.SetFileName('cg.mha')
writer.Execute(recon)