.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/03_algorithms/03_run_fbp.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_03_algorithms_03_run_fbp.py: 2D non-TOF filtered back projection (FBP) of Poisson data ========================================================= This example demonstrates the run 2D filtered back projection (FBP) on pre-corrected Poisson emission data. .. GENERATED FROM PYTHON SOURCE LINES 10-18 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import parallelproj.projectors from parallelproj import to_numpy_array from utils import RadonDisk, RadonObjectSequence from scipy.ndimage import gaussian_filter .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: Python from array_utils import suggest_array_backend_and_device # To use a specific backend and/or device, replace the None arguments, e.g.: # xp, dev = suggest_array_backend_and_device(backend="numpy", dev="cpu") or by setting xp and dev manually xp, dev = suggest_array_backend_and_device(None, None) .. rst-class:: sphx-glr-script-out .. code-block:: none Using array API: array_api_compat.torch, device: cpu .. GENERATED FROM PYTHON SOURCE LINES 27-29 Input parameters ---------------- .. GENERATED FROM PYTHON SOURCE LINES 29-46 .. code-block:: Python # add Poisson noise to the sinogram add_noise = True # total "prompt" counts of the emission sinogram total_counts = 1e7 # sigma of the Gaussian filter to smooth sinogram in radial # direction to simulate limited resolution sino_res = 1.0 # linear attenuation coefficient of phantom in 1/cm mu = 0.1 # number of MLEM iterations to run num_iter = 50 # maximum angle of the sinogram in radians phi_max = xp.pi # undersampling factor of the sinogram in the angular direction phi_undersampling_factor = 1 .. GENERATED FROM PYTHON SOURCE LINES 47-49 Setup of arrays --------------- .. GENERATED FROM PYTHON SOURCE LINES 49-66 .. code-block:: Python num_rad = 201 num_phi = int(0.5 * num_rad * xp.pi * (phi_max / xp.pi)) + 1 num_phi = num_phi // phi_undersampling_factor r = xp.linspace(-30, 30, num_rad, device=dev, dtype=xp.float32) phi = xp.linspace(0, phi_max, num_phi, endpoint=False, device=dev, dtype=xp.float32) R, PHI = xp.meshgrid(r, phi, indexing="ij") X0, X1 = xp.meshgrid(r, r, indexing="ij") x = xp.linspace(float(xp.min(r)), float(xp.max(r)), 1001, device=dev, dtype=xp.float32) X0hr, X1hr = xp.meshgrid(x, x, indexing="ij") print(f"num rad: {num_rad}") print(f"phi max: {180*phi_max/xp.pi:.2f} deg") print(f"delta phi: {180*float(phi[1]-phi[0])/xp.pi:.2f} deg") .. rst-class:: sphx-glr-script-out .. code-block:: none num rad: 201 phi max: 180.00 deg delta phi: 0.57 deg .. GENERATED FROM PYTHON SOURCE LINES 67-71 Define an object with known Radon transform ------------------------------------------- We setup a combination of (scaled) disks as a simple phantom. .. GENERATED FROM PYTHON SOURCE LINES 71-105 .. code-block:: Python disk0 = RadonDisk(xp, dev, 8.0) disk0.amplitude = 1.0 disk0.s0 = 3.0 disk1 = RadonDisk(xp, dev, 2.0) disk1.amplitude = 0.5 disk1.x1_offset = 4.67 disk2 = RadonDisk(xp, dev, 1.4) disk2.amplitude = -0.5 disk2.x0_offset = -10.0 disk3 = RadonDisk(xp, dev, 0.93) disk3.amplitude = -0.5 disk3.x1_offset = -4.67 disk4 = RadonDisk(xp, dev, 0.67) disk4.amplitude = 1.0 disk4.x1_offset = -4.67 radon_object = RadonObjectSequence([disk0, disk1, disk2, disk3, disk4]) fig, ax = plt.subplots(tight_layout=True) ax.imshow( to_numpy_array(radon_object.values(X0hr, X1hr).T), cmap="Greys", origin="lower", ) ax.set_xlabel(r"$x_0$") ax.set_ylabel(r"$x_1$") ax.set_title("true object", fontsize="medium") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_001.png :alt: true object :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 106-108 Calculate the radon transform of the object ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 108-111 .. code-block:: Python rt_transform = radon_object.radon_transform(R, PHI) .. GENERATED FROM PYTHON SOURCE LINES 112-114 Simulate the effect of attenuation and Poisson noise ---------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 114-148 .. code-block:: Python sens_sino = xp.exp(-mu * disk0.radon_transform(R, PHI)) contam = xp.full( rt_transform.shape, 0.1 * xp.mean(sens_sino * rt_transform), device=dev ) emis_sino = sens_sino * rt_transform if sino_res > 0: for i in range(num_phi): emis_sino[:, i] = xp.asarray( gaussian_filter( to_numpy_array(emis_sino[:, i]), sino_res, ), device=dev, ) emis_sino = emis_sino + contam count_fac = total_counts / float(xp.sum(emis_sino)) emis_sino *= count_fac contam *= count_fac if add_noise: emis_sino = xp.asarray( np.random.poisson(to_numpy_array(emis_sino)).astype(np.float32), device=dev, ) # pre-correct sinogram - needed to run FBP pre_corrected_sino = (emis_sino - contam) / sens_sino .. GENERATED FROM PYTHON SOURCE LINES 149-184 .. code-block:: Python ext_sino = [float(xp.min(r)), float(xp.max(r)), float(xp.min(phi)), float(xp.max(phi))] fig, ax = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True) ax[0].imshow( to_numpy_array(rt_transform.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) ax[1].imshow( to_numpy_array(emis_sino.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) ax[2].imshow( to_numpy_array(pre_corrected_sino.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) for axx in ax.ravel(): axx.set_xlabel(r"$s$") axx.set_ylabel(r"$\phi$") ax[0].set_title("radon transform of object", fontsize="medium") ax[1].set_title("emission sinogram", fontsize="medium") ax[2].set_title("pre-corrected sinogram", fontsize="medium") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_002.png :alt: radon transform of object, emission sinogram, pre-corrected sinogram :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 185-187 Setup of the ramp filter ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 187-202 .. code-block:: Python n_filter = r.shape[0] r_shift = xp.arange(n_filter, device=dev, dtype=xp.float32) - n_filter // 2 f = xp.zeros(n_filter, device=dev, dtype=xp.float32) f[r_shift != 0] = -1 / (xp.pi**2 * r_shift[r_shift != 0] ** 2) f[(r_shift % 2) == 0] = 0 f[r_shift == 0] = 0.25 fig, ax = plt.subplots(tight_layout=True) ax.plot(to_numpy_array(r_shift), to_numpy_array(f), ".-") ax.set_xlabel(r"$s$") ax.set_ylabel(r"$f$") ax.set_title("ramp filter", fontsize="medium") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_003.png :alt: ramp filter :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-205 Ramp filtering of the pre-corrected emission sinogram ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 205-242 .. code-block:: Python filtered_pre_corrected_sino = 1.0 * pre_corrected_sino for i in range(num_phi): filtered_pre_corrected_sino[:, i] = xp.asarray( np.convolve( to_numpy_array(filtered_pre_corrected_sino[:, i]), to_numpy_array(f), mode="same", ), device=dev, ) fig, ax = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True) ax[0].imshow( to_numpy_array(pre_corrected_sino.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) ax[1].imshow( to_numpy_array(filtered_pre_corrected_sino.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) for axx in ax.ravel(): axx.set_xlabel(r"$s$") axx.set_ylabel(r"$\phi$") ax[0].set_title("pre-corrected sinogram", fontsize="medium") ax[1].set_title("filtered pre-corr. sinogram", fontsize="medium") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_004.png :alt: pre-corrected sinogram, filtered pre-corr. sinogram :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 243-245 Define a projector and run filtered back projection (FBP) --------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 245-264 .. code-block:: Python proj = parallelproj.projectors.ParallelViewProjector2D( (num_rad, num_rad), r, -phi, 2 * float(xp.max(r)), (float(xp.min(r)), float(xp.min(r))), (float(r[1] - r[0]), float(r[1] - r[0])), ) filtered_back_proj = proj.adjoint(filtered_pre_corrected_sino) fig, ax = plt.subplots(tight_layout=True) ax.imshow(to_numpy_array(filtered_back_proj).T, cmap="Greys", origin="lower") ax.set_xlabel(r"$x_0$") ax.set_ylabel(r"$x_1$") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_005.png :alt: 03 run fbp :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 265-267 Run iterative MLEM reconstruction for comparison ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 267-280 .. code-block:: Python x_mlem = xp.ones((num_rad, num_rad), device=dev, dtype=xp.float32) sens_img = proj.adjoint(sens_sino) for i in range(num_iter): exp = sens_sino * proj(x_mlem) + contam ratio = emis_sino / exp ratio_back = proj.adjoint(sens_sino * ratio) update_img = ratio_back / sens_img x_mlem *= update_img .. GENERATED FROM PYTHON SOURCE LINES 281-283 Visualize the results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 283-327 .. code-block:: Python ext_img = [float(xp.min(r)), float(xp.max(r)), float(xp.min(r)), float(xp.max(r))] fig, ax = plt.subplots(1, 4, figsize=(16, 4), tight_layout=True) ax[0].imshow( to_numpy_array(radon_object.values(X0hr, X1hr).T), cmap="Greys", extent=ext_img, origin="lower", ) ax[1].imshow( to_numpy_array(emis_sino.T), cmap="Greys", aspect=20, extent=ext_sino, origin="lower", ) ax[2].imshow( to_numpy_array(filtered_back_proj.T), cmap="Greys", extent=ext_img, origin="lower", ) ax[3].imshow( to_numpy_array(x_mlem.T), cmap="Greys", extent=ext_img, origin="lower", ) ax[0].set_xlabel(r"$x_0$") ax[0].set_ylabel(r"$x_1$") ax[1].set_xlabel(r"$s$") ax[1].set_ylabel(r"$\phi$") ax[2].set_xlabel(r"$x_0$") ax[2].set_ylabel(r"$x_1$") ax[3].set_xlabel(r"$x_0$") ax[3].set_ylabel(r"$x_1$") ax[0].set_title("true object", fontsize="medium") ax[1].set_title("emission sino", fontsize="medium") ax[2].set_title("filtered back projection", fontsize="medium") ax[3].set_title(f"MLEM {num_iter} it.", fontsize="medium") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_006.png :alt: true object, emission sino, filtered back projection, MLEM 50 it. :srcset: /auto_examples/03_algorithms/images/sphx_glr_03_run_fbp_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 18.278 seconds) .. _sphx_glr_download_auto_examples_03_algorithms_03_run_fbp.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03_run_fbp.ipynb <03_run_fbp.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03_run_fbp.py <03_run_fbp.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03_run_fbp.zip <03_run_fbp.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_