Examples¶
Comparison of parallel-beam reconstruction methods¶
This example illustrates the performance of the
different reconstruction techniques for a parallel-beam
geometry. The left column shows the reconstruction of
the original image and the right column shows the reconstruction
of the corresponding binary images. Note that the
SART process could be sped-up by computing an
initial guess with a non-iterative method and
setting it with the initial
keyword argument.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | from matplotlib import pylab as plt
import numpy as np
import radontea
from radontea.logo import get_original
N = 55 # image size
A = 13 # number of sinogram angles
ITA = 10 # number of iterations a
ITB = 100 # number of iterations b
angles = np.linspace(0, np.pi, A)
im = get_original(N)
sino = radontea.radon_parallel(im, angles)
fbp = radontea.backproject(sino, angles)
fintp = radontea.fourier_map(sino, angles).real
sarta = radontea.sart(sino, angles, iterations=ITA)
sartb = radontea.sart(sino, angles, iterations=ITB)
im2 = (im >= (im.max() / 5)) * 255
sino2 = radontea.radon_parallel(im2, angles)
fbp2 = radontea.backproject(sino2, angles)
fintp2 = radontea.fourier_map(sino2, angles).real
sarta2 = radontea.sart(sino2, angles, iterations=ITA)
sartb2 = radontea.sart(sino2, angles, iterations=ITB)
plt.figure(figsize=(8, 22))
pltkw = {"vmin": -20,
"vmax": 280}
plt.subplot(6, 2, 1, title="original image")
plt.imshow(im, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 2, title="binary image")
plt.imshow(im2, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 3, title="sinogram (from original)")
plt.imshow(sino)
plt.axis('off')
plt.subplot(6, 2, 4, title="sinogram (from binary)")
plt.imshow(sino2)
plt.axis('off')
plt.subplot(6, 2, 5, title="filtered backprojection")
plt.imshow(fbp, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 6, title="filtered backprojection")
plt.imshow(fbp2, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 7, title="Fourier interpolation")
plt.imshow(fintp, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 8, title="Fourier interpolation")
plt.imshow(fintp2, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 9, title="SART ({} iterations)".format(ITA))
plt.imshow(sarta, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 10, title="SART ({} iterations)".format(ITA))
plt.imshow(sarta2, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 11, title="SART ({} iterations)".format(ITB))
plt.imshow(sartb, **pltkw)
plt.axis('off')
plt.subplot(6, 2, 12, title="SART ({} iterations)".format(ITB))
plt.imshow(sartb2, **pltkw)
plt.axis('off')
plt.tight_layout()
plt.show()
|
Volumetric data reconstruction benchmark¶
This simple example can be used to quantify the speed-up due to
multiprocessing on the hardware used. Note that
multiprocessing.cpu_count
does not return the number of
physical cores.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | from multiprocessing import cpu_count
import time
import numpy as np
import radontea as rt
A = 70 # number of angles
N = 128 # detector size x
M = 24 # detector size y (number of slices)
# generate random data
sino0 = np.random.random((A, N)) # for 2d example
sino = np.random.random((A, M, N)) # for 3d example
sino[:, 0, :] = sino0
angles = np.linspace(0, np.pi, A) # for both
a = time.time()
data1 = rt.backproject_3d(sino, angles, ncpus=1)
print("time on 1 core: {:.2f} s".format(time.time() - a))
a = time.time()
data2 = rt.backproject_3d(sino, angles, ncpus=cpu_count())
print("time on {} cores: {:.2f} s".format(cpu_count(), time.time() - a))
assert np.all(data1 == data2), "2D and 3D results don't match"
|
Progress monitoring with progression¶
The progress of the reconstruction algorithms in radontea can be tracked with other packages such as progression.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | from multiprocessing import cpu_count
import numpy as np
import progression as pr
import radontea as rt
A = 55 # number of angles
N = 128 # detector size x
M = 24 # detector size y (number of slices)
# generate random data
sino0 = np.random.random((A, N))
sino = np.random.random((A, M, N))
sino[:, 0, :] = sino0
angles = np.linspace(0, np.pi, A)
count = pr.UnsignedIntValue()
max_count = pr.UnsignedIntValue()
with pr.ProgressBar(count=count,
max_count=max_count,
interval=0.3
) as pb:
pb.start()
rt.volume_recon(func2d=rt.sart, sinogram=sino, angles=angles,
ncpus=cpu_count(), count=count, max_count=max_count)
|