Skip to content

Instantly share code, notes, and snippets.

@adler-j
Last active February 22, 2017 09:07
Show Gist options
  • Save adler-j/e9eddcf4ef919042a72178e523429142 to your computer and use it in GitHub Desktop.
Save adler-j/e9eddcf4ef919042a72178e523429142 to your computer and use it in GitHub Desktop.
Using odl projector, CPU vs CUDA
import odl
import numpy as np
# Discrete reconstruction space: discretized functions on the rectangle
# [-20, 20]^2 with 300 samples per dimension.
reco_space = odl.uniform_discr(
min_pt=[-1, -1], max_pt=[1, 1], shape=[256, 256])
# Make a parallel beam geometry with flat detector
# Angles: uniformly spaced, n = 180, min = 0, max = pi
angle_partition = odl.uniform_partition(0, np.pi, 180)
# Detector: uniformly sampled, n = 256, min = -1.152, max = 1.152
detector_partition = odl.uniform_partition(-1.152, 1.152, 256)
geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
# Create a discrete Shepp-Logan phantom (modified version)
phantom = odl.phantom.shepp_logan(reco_space, modified=True)
# ray transform aka forward projection. We use ASTRA CUDA/CPU backend.
ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
ray_trafo_cpu = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cpu')
# Create projection data by calling the ray transform on the phantom
proj_data = ray_trafo(phantom)
proj_data_cpu = ray_trafo_cpu(phantom)
# Back-projection can be done by simply calling the adjoint operator on the
# projection data (or any element in the projection space).
backproj = ray_trafo.adjoint(proj_data)
backproj_cpu = ray_trafo_cpu.adjoint(proj_data_cpu)
# Shows a slice of the phantom, projections, and reconstruction
phantom.show(title='Phantom')
proj_data.show(title='cuda sinogram')
backproj.show(title='cuda backprojection')
proj_data_cpu.show(title='cpu sinogram')
backproj_cpu.show(title='cpu backprojection')
(proj_data_cpu - proj_data).show(title='diff sinogram')
(backproj_cpu - backproj).show(title='diff backprojection')
@rshpeley
Copy link

I saw your code example here and decided to try it, but I get an error from reco_space = odl.uniform_discr()

TypeError: uniform_discr() missing 3 required positional arguments: 'min_pt', 'max_pt', and 'shape'

@adler-j
Copy link
Author

adler-j commented Feb 22, 2017

Thank you for the update, we did some renames of the min_corner etc parameters in ODL that broke this example. I have now updated it.

Please ask if you have further questions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment