Merge pull request #532 from fedor-goncharov/tutorial-ctscan-correction

remark on classical sinograms added + jacobian correction; added two …

README.md

@@ -148,3 +148,4 @@ A list of video tutorials to learn more about PyLops:
 * Aniket Singh Rawat, dikwickley
 * Rohan Babbar, rohanbabbar04
 * Wei Zhang, ZhangWeiGeo
+* Fedor Goncharov, fedor-goncharov

docs/source/credits.rst

@@ -19,4 +19,5 @@ Contributors
 *  `Juan Daniel Romero <https://github.com/jdromerom>`_, jdromerom
 *  `Aniket Singh Rawat <https://github.com/dikwickley>`_, dikwickley
 *  `Rohan Babbar <https://github.com/rohanbabbar04>`_, rohanbabbar04
-*  `Wei Zhang <https://github.com/ZhangWeiGeo>`_, ZhangWeiGeo
+*  `Wei Zhang <https://github.com/ZhangWeiGeo>`_, ZhangWeiGeo
+*  `Fedor Goncharov <https://github.com/fedor-goncharov>`_, fedor-goncharov

tutorials/ctscan.py

@@ -28,11 +28,17 @@
 # such a type:
 #
 # .. math::
-#    t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)}
+#    t\sin(\theta) + x\cos(\theta) = r,
 #
 # where :math:`\theta` is the angle between the x-axis (:math:`x`) and
 # the perpendicular to the summation line and :math:`r` is the distance
-# from the origin of the summation line.
+# from the origin of the summation line. Radon transform in CT
+# corresponds to the integral of the input image along the straight line above.
+# To implement the integration in PyLops we simply need to express
+# :math:`t(r,\theta;x)` which is given by:
+#
+# .. math::
+#    t(r,\theta; x) = \tan\left(\frac{\pi}{2}-\theta\right)x + \frac{r}{\sin(\theta)}.
 
 
 @jit(nopython=True)
@@ -44,6 +50,13 @@ def radoncurve(x, r, theta):
     )
 
 
+###############################################################################
+# Note that in the above implementation we added centering :math:`t \mapsto t - n_y/2` and
+# :math:`r \mapsto r - n_y/2` so that the origin of the integration lines is exactly in the
+# center of the image (centering for :math:`x` is not needed because we will use
+# ``centeredh=True`` in the constructor of ``Radon2D``).
+
+
 x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T
 x = x / x.max()
 nx, ny = x.shape
@@ -81,14 +94,36 @@ def radoncurve(x, r, theta):
 axs[1].set_title("Data")
 axs[1].axis("tight")
 axs[2].imshow(xrec.T, cmap="gray")
-axs[2].set_title("Adjoint model")
+axs[2].set_title("Adjoint model in PyLops")
 axs[2].axis("tight")
 fig.tight_layout()
 
+###############################################################################
+# Note that our raw data ``y`` *does not represent exactly* classical sinograms
+# in medical imaging. Integration along curves in the adjoint form of
+# :func:`pylops.signalprocessing.Radon2D` is performed with respect to
+# :math:`dx`, whereas canonically it is assumed to be with respect to the natural
+# parametrization :math:`dl = \sqrt{(dx)^2 + (dt)^2}`. To retrieve back the
+# classical sinogram we have to divide data by the jacobian
+# :math:`j(x,l) = \left\vert dx/dl \right\vert = |\sin(\theta)|`.
+
+sinogram = np.divide(
+    y.T, np.abs(np.sin(theta) + 1e-15)
+)  # small shift to avoid zero-division
+fig, axs = plt.subplots(1, 2, figsize=(10, 4))
+axs[0].imshow(y.T, cmap="gray")
+axs[0].set_title("Data")
+axs[0].axis("tight")
+axs[1].imshow(sinogram, cmap="gray")
+axs[1].set_title("Sinogram in medical imaging")
+axs[1].axis("tight")
+fig.tight_layout()
 
 ###############################################################################
-# Finally we take advantage of our different solvers and try to invert the
-# modelling operator both in a least-squares sense and using TV-reg.
+# From now on, we will not pursue further working with the "true sinogram", instead
+# we will reconstruct the original phantom directly from ``y``. For this we take advantage
+# of our different solvers and try to invert the modelling operator both in a
+# least-squares sense and using TV-reg.
 Dop = [
     pylops.FirstDerivative(
         (nx, ny), axis=0, edge=True, kind="backward", dtype=np.float64