add example (close #1)

docs/examples.rst

@@ -0,0 +1,5 @@
+========
+Examples
+========
+
+.. fancy_include:: comparison_parallel.py

docs/index.rst

@@ -11,6 +11,7 @@ This is the documentaion of radontea version |release|.
 
    introduction
    code_reference
+   examples
 
 
 Indices and tables

examples/comparison_parallel.py

@@ -0,0 +1,93 @@
+"""Comparison of parallel reconstruction methods
+
+This example illustrates the performance of the
+different reconstruction techniques. 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.
+"""
+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(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(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()

examples/generate_example_images.py

@@ -0,0 +1,33 @@
+import os
+import os.path as op
+import sys
+
+import matplotlib.pylab as plt
+
+thisdir = op.dirname(op.abspath(__file__))
+sys.path.insert(0, op.dirname(thisdir))
+
+DPI = 80
+
+
+if __name__ == "__main__":
+    # Do not display example plots
+    plt.show = lambda: None
+    files = os.listdir(thisdir)
+    files = [f for f in files if f.endswith(".py")]
+    files = [f for f in files if not f == op.basename(__file__)]
+    files = sorted([op.join(thisdir, f) for f in files])
+
+    for f in files:
+        fname = f[:-3] + ".jpg"
+        if not op.exists(fname):
+            exec_str = open(f).read()
+            if exec_str.count("plt.show()"):
+                exec(exec_str)
+                plt.savefig(fname, dpi=DPI)
+                print("Image created: '{}'".format(fname))
+            else:
+                print("No image: '{}'".format(fname))
+        else:
+            print("Image skipped (already exists): '{}'".format(fname))
+        plt.close()