From 71c708442ed5d7ce3dfe20b725665e978bae8885 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Fri, 12 Jul 2024 10:50:34 -0600
Subject: [PATCH 01/15] Initial work on improved astra 3d docs

---
 scico/linop/xray/astra.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index aa530a6ef..4453a71af 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -227,6 +227,18 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
     Perform tomographic projection (also called X-ray projection) of a
     volume at specified angles, using the
     `ASTRA toolbox <https://github.com/astra-toolbox/astra-toolbox>`_.
+    The `3D geometries <https://astra-toolbox.com/docs/geom3d.html>`__
+    "parallel3d" and "parallel3d_vec" are supported by this interface.
+
+    The reconstruction volume is fixed with respect to the coordinate
+    system, with the volume centred at the origin, and the unit length
+    (in arbitrary units) of the sides of the voxels define the scale
+    for all other dimensions in the source-volume-detector configuration.
+    Geometry axes x, y, and z correspond to volume array axes 0, 1, and 2
+    respectively
+
+    In the "parallel3d" case, the source and detector rotate clockwise
+    about the z axis in the x-y plane.
     """
 
     def __init__(

From 0d6278065bfd7084032e0e68001e6171aaf71a7f Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 12:27:18 -0600
Subject: [PATCH 02/15] Typo fix

---
 scico/linop/xray/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scico/linop/xray/__init__.py b/scico/linop/xray/__init__.py
index 3a5dc7d4a..cb47c9ec6 100644
--- a/scico/linop/xray/__init__.py
+++ b/scico/linop/xray/__init__.py
@@ -27,7 +27,7 @@
    :include-source: False
    :show-source-link: False
    :caption: Comparison of 2D X-ray projector geometries. The red arrows
-      are are directed towards the detector, which is oriented with pixel
+      are directed towards the detector, which is oriented with pixel
       indices ordered in the same direction as clockwise rotation (e.g.
       in the "scico" geometry, the :math:`\theta=0` projection
       corresponds to row sums ordered from the top to the bottom of the

From 1379278197cbb830d357af8d58e2b598de6235db Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 12:28:09 -0600
Subject: [PATCH 03/15] Improve 3d astra docs

---
 docs/source/pyfigures/xray_3d_ang.py | 72 +++++++++++++++++++++
 docs/source/pyfigures/xray_3d_vec.py | 80 ++++++++++++++++++++++++
 docs/source/pyfigures/xray_3d_vol.py | 93 ++++++++++++++++++++++++++++
 3 files changed, 245 insertions(+)
 create mode 100644 docs/source/pyfigures/xray_3d_ang.py
 create mode 100644 docs/source/pyfigures/xray_3d_vec.py
 create mode 100644 docs/source/pyfigures/xray_3d_vol.py

diff --git a/docs/source/pyfigures/xray_3d_ang.py b/docs/source/pyfigures/xray_3d_ang.py
new file mode 100644
index 000000000..5b0c677db
--- /dev/null
+++ b/docs/source/pyfigures/xray_3d_ang.py
@@ -0,0 +1,72 @@
+import numpy as np
+
+import matplotlib as mpl
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+
+mpl.rcParams["savefig.transparent"] = True
+
+
+c = 1.0 / np.sqrt(2.0)
+e = 1e-2
+style = "Simple, tail_width=0.5, head_width=4, head_length=8"
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))
+ax.set_aspect(1.0)
+ax.set_xlim(-1.1, 1.1)
+ax.set_ylim(-1.1, 1.1)
+ax.set_xticks(np.linspace(-1.0, 1.0, 5))
+ax.set_yticks(np.linspace(-1.0, 1.0, 5))
+ax.tick_params(axis="x", labelsize=12)
+ax.tick_params(axis="y", labelsize=12)
+ax.set_xlabel("$x$", fontsize=14)
+ax.set_ylabel("$y$", fontsize=14)
+
+plist = [
+    patches.FancyArrowPatch((0.0, -1.0), (0.0, -0.5), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((c, -c), (c / 2.0, -c / 2.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((1.0, 0.0), (0.5, 0.0), arrowstyle=style, color="r"),
+    patches.Arc((0.0, 0.0), 2.0, 2.0, theta1=-90, theta2=45.0, color="b", ls="dotted"),
+    patches.FancyArrowPatch((c + e, c - e), (c - e, c + e), arrowstyle=style, color="b"),
+]
+for p in plist:
+    ax.add_patch(p)
+ax.text(0.02, -0.75, r"$\theta=0$", color="r", fontsize=14)
+ax.text(
+    3 * c / 4 + 0.01,
+    -3 * c / 4 + 0.01,
+    r"$\theta=\frac{\pi}{4}$",
+    color="r",
+    fontsize=14,
+)
+ax.text(0.65, 0.05, r"$\theta=\frac{\pi}{2}$", color="r", fontsize=14)
+
+ax.plot((-0.375, 0.375), (1.0, 1.0), color="orange", lw=2)
+ax.arrow(
+    -0.375,
+    0.94,
+    0.75,
+    0.0,
+    color="orange",
+    lw=0.5,
+    ls="--",
+    head_width=0.03,
+    length_includes_head=True,
+)
+ax.text(0.0, 0.82, r"$\theta=0$", color="orange", ha="center", fontsize=14)
+
+ax.plot((-1.0, -1.0), (-0.375, 0.375), color="orange", lw=2)
+ax.arrow(
+    -0.94,
+    -0.375,
+    0.0,
+    0.75,
+    color="orange",
+    lw=0.5,
+    ls="--",
+    head_width=0.03,
+    length_includes_head=True,
+)
+ax.text(-0.9, 0.0, r"$\theta=\frac{\pi}{2}$", color="orange", ha="left", fontsize=14)
+
+fig.tight_layout()
+fig.show()
diff --git a/docs/source/pyfigures/xray_3d_vec.py b/docs/source/pyfigures/xray_3d_vec.py
new file mode 100644
index 000000000..af26c3645
--- /dev/null
+++ b/docs/source/pyfigures/xray_3d_vec.py
@@ -0,0 +1,80 @@
+import numpy as np
+
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+from matplotlib.patches import FancyArrowPatch
+from mpl_toolkits.mplot3d import proj3d
+
+mpl.rcParams["savefig.transparent"] = True
+
+
+# See https://github.com/matplotlib/matplotlib/issues/21688
+class Arrow3D(FancyArrowPatch):
+    def __init__(self, xs, ys, zs, *args, **kwargs):
+        FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
+        self._verts3d = xs, ys, zs
+
+    def do_3d_projection(self, renderer=None):
+        xs3d, ys3d, zs3d = self._verts3d
+        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M)
+        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
+
+        return np.min(zs)
+
+
+# Define vector components
+𝜃 = 10 * np.pi / 180.0  # angle in x-y plane (azimuth angle)
+𝛼 = 70 * np.pi / 180.0  # angle with z axis (zenith angle)
+𝛥p, 𝛥d = 0.3, 1.0
+d = (-𝛥d * np.sin(𝛼) * np.sin(𝜃), 𝛥d * np.sin(𝛼) * np.cos(𝜃), 𝛥d * np.cos(𝛼))
+u = (𝛥p * np.cos(𝜃), 𝛥p * np.sin(𝜃), 0.0)
+v = (𝛥p * np.cos(𝛼) * np.sin(𝜃), -𝛥p * np.cos(𝛼) * np.cos(𝜃), 𝛥p * np.sin(𝛼))
+
+# Location of text labels
+d_txtpos = np.array(d) + np.array([0, 0, -0.12])
+u_txtpos = np.array(d) + np.array(u) + np.array([0, 0, -0.1])
+v_txtpos = np.array(d) + np.array(v) + np.array([0, 0, 0.03])
+
+
+arrowstyle = "-|>,head_width=2.5,head_length=9"
+
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+# Set view
+ax.set_aspect("equal")
+ax.elev = 15
+ax.azim = -50
+ax.set_box_aspect(None, zoom=2)
+ax.set_xlim((-1.1, 1.1))
+ax.set_ylim((-1.1, 1.1))
+ax.set_zlim((-1.1, 1.1))
+
+# Disable shaded 3d axis grids
+ax.set_axis_off()
+
+# Draw central x,y,z axes and labels
+axis_crds = np.array([[-1, 1], [0, 0], [0, 0]])
+axis_lbls = ("$x$", "$y$", "$z$")
+for k in range(3):
+    crd = np.roll(axis_crds, k, axis=0)
+    ax.add_artist(
+        Arrow3D(
+            *crd.tolist(),
+            lw=1.5,
+            ls="--",
+            arrowstyle=arrowstyle,
+            color="black",
+        )
+    )
+    ax.text(*(1.05 * crd[:, 1]).tolist(), axis_lbls[k], fontsize=12)
+
+# Draw d, u, v and labels
+ax.quiver(0, 0, 0, *d, arrow_length_ratio=0.08, lw=2, color="blue")
+ax.quiver(*d, *u, arrow_length_ratio=0.08 / 𝛥p, lw=2, color="blue")
+ax.quiver(*d, *v, arrow_length_ratio=0.08 / 𝛥p, lw=2, color="blue")
+ax.text(*d_txtpos, r"$\mathbf{d}$", fontsize=12)
+ax.text(*u_txtpos, r"$\mathbf{u}$", fontsize=12)
+ax.text(*v_txtpos, r"$\mathbf{v}$", fontsize=12)
+
+fig.tight_layout()
+fig.show()
diff --git a/docs/source/pyfigures/xray_3d_vol.py b/docs/source/pyfigures/xray_3d_vol.py
new file mode 100644
index 000000000..2811933a6
--- /dev/null
+++ b/docs/source/pyfigures/xray_3d_vol.py
@@ -0,0 +1,93 @@
+import numpy as np
+
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+from matplotlib.patches import FancyArrowPatch
+from mpl_toolkits.mplot3d import proj3d
+
+mpl.rcParams["savefig.transparent"] = True
+
+
+# See https://github.com/matplotlib/matplotlib/issues/21688
+class Arrow3D(FancyArrowPatch):
+    def __init__(self, xs, ys, zs, *args, **kwargs):
+        FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
+        self._verts3d = xs, ys, zs
+
+    def do_3d_projection(self, renderer=None):
+        xs3d, ys3d, zs3d = self._verts3d
+        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M)
+        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
+
+        return np.min(zs)
+
+
+# Define vector components
+𝜃 = 10 * np.pi / 180.0  # angle in x-y plane (azimuth angle)
+𝛼 = 70 * np.pi / 180.0  # angle with z axis (zenith angle)
+𝛥p, 𝛥d = 0.3, 1.0
+d = (-𝛥d * np.sin(𝛼) * np.sin(𝜃), 𝛥d * np.sin(𝛼) * np.cos(𝜃), 𝛥d * np.cos(𝛼))
+u = (𝛥p * np.cos(𝜃), 𝛥p * np.sin(𝜃), 0.0)
+v = (𝛥p * np.cos(𝛼) * np.sin(𝜃), -𝛥p * np.cos(𝛼) * np.cos(𝜃), 𝛥p * np.sin(𝛼))
+
+# Location of text labels
+d_txtpos = np.array(d) + np.array([0, 0, -0.12])
+u_txtpos = np.array(d) + np.array(u) + np.array([0, 0, -0.1])
+v_txtpos = np.array(d) + np.array(v) + np.array([0, 0, 0.03])
+
+
+arrowstyle = "-|>,head_width=2.5,head_length=9"
+
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+# Set view
+ax.set_aspect("equal")
+ax.elev = 40
+ax.azim = -60
+ax.set_box_aspect(None, zoom=1.8)
+ax.set_xlim((-10.5, 10.5))
+ax.set_ylim((-10.5, 10.5))
+ax.set_zlim((-10.5, 10.5))
+
+# Disable shaded 3d axis grids
+ax.set_axis_off()
+
+# Draw central x,y,z axes and labels
+axis_crds = np.array([[-10, 10], [0, 0], [0, 0]])
+axis_lbls = ("$x$", "$y$", "$z$")
+for k in range(3):
+    crd = np.roll(axis_crds, k, axis=0)
+    ax.add_artist(
+        Arrow3D(
+            *crd.tolist(),
+            lw=1.5,
+            ls="--",
+            arrowstyle=arrowstyle,
+            color="black",
+        )
+    )
+    ax.text(*(1.05 * crd[:, 1]).tolist(), axis_lbls[k], fontsize=12)
+
+wx = 4
+wy = 3
+wz = 2
+bx = np.array([-wx, wx, wx, wx, -wx, -wx, -wx])
+by = np.array([-wy, -wy, wy, wy, wy, -wy, -wy])
+bz = np.array([-wz, -wz, -wz, wz, wz, wz, -wz])
+ax.plot(bx, by, bz, lw=2, color="blue")
+ax.plot(bx[0:3], by[0:3], -bz[0:3], lw=2, color="blue")
+bx = np.array([wx, wx])
+by = np.array([-wy, -wy])
+bz = np.array([-wz, wz])
+ax.plot(bx, by, bz, lw=2, color="blue")
+bx = np.array([-wx, -wx, wx])
+by = np.array([-wy, wy, wy])
+bz = np.array([-wz, -wz, -wz])
+ax.plot(bx, by, bz, lw=2, ls="--", color="blue")
+bx = np.array([-wx, -wx])
+by = np.array([wy, wy])
+bz = np.array([-wz, wz])
+ax.plot(bx, by, bz, lw=2, ls="--", color="blue")
+
+fig.tight_layout()
+fig.show()

From 2e507d7fd3998503dffd6ccbe523ceb5fb4d8ff5 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 18:08:51 -0600
Subject: [PATCH 04/15] Improve 3d astra docs

---
 scico/linop/xray/astra.py | 61 ++++++++++++++++++++++++++++++++-------
 1 file changed, 50 insertions(+), 11 deletions(-)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index 4453a71af..920b6a303 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -227,18 +227,58 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
     Perform tomographic projection (also called X-ray projection) of a
     volume at specified angles, using the
     `ASTRA toolbox <https://github.com/astra-toolbox/astra-toolbox>`_.
-    The `3D geometries <https://astra-toolbox.com/docs/geom3d.html>`__
+    The `3D geometries <https://astra-toolbox.com/docs/geom3d.html#projection-geometries>`__
     "parallel3d" and "parallel3d_vec" are supported by this interface.
 
     The reconstruction volume is fixed with respect to the coordinate
-    system, with the volume centred at the origin, and the unit length
-    (in arbitrary units) of the sides of the voxels define the scale
-    for all other dimensions in the source-volume-detector configuration.
-    Geometry axes x, y, and z correspond to volume array axes 0, 1, and 2
-    respectively
+    system, with the volume centred at the origin, as illustrated below:
+
+    .. plot:: pyfigures/xray_3d_vol.py
+       :align: center
+       :include-source: False
+       :show-source-link: False
+
+    The voxels sides have unit length (in arbitrary units), which defines
+    the scale for all other dimensions in the source-volume-detector
+    configuration. Geometry axes `x`, `y`, and `z` correspond to volume
+    array axes 0, 1, and 2 respectively.
 
     In the "parallel3d" case, the source and detector rotate clockwise
-    about the z axis in the x-y plane.
+    about the `z` axis in the `x`-`y` plane, as illustrated below:
+
+    .. plot:: pyfigures/xray_3d_ang.py
+       :align: center
+       :include-source: False
+       :show-source-link: False
+       :caption: The red arrows indicate the direction of the beam towards
+          the detector (orange) and the arrows parallel to the detector
+          indicate the direction of increasing pixel indices.
+
+    In the "parallel3d_vec" case, each view is determined by the following
+    vectors
+
+    .. list-table:: View definition vectors
+       :widths: 10 90
+
+       * - :math:`\mb{r}`
+         - Direction of the parallel beam
+       * - :math:`\mb{d}`
+         - Center of the detector
+       * - :math:`\mb{u}`
+         - Vector from detector pixel (0,0) to (0,1)
+       * - :math:`\mb{v}`
+         - Vector from detector pixel (0,0) to (1,0)
+
+    .. plot:: pyfigures/xray_3d_vec.py
+       :align: center
+       :include-source: False
+       :show-source-link: False
+
+    which are concatenated into a single row vector
+    :math:`(\mb{r}, \mb{d}, \mb{u}, \mb{v})` to form the full
+    geometry specification for a single view. Multiple such
+    row vectors are stacked to specify the geometry for a set
+    of views.
     """
 
     def __init__(
@@ -250,11 +290,10 @@ def __init__(
         vectors: Optional[np.ndarray] = None,
     ):
         """
-        This class supports both "parallel3d" and "parallel3d_vec" astra
-        `projection geometries <https://www.astra-toolbox.com/docs/geom3d.html#projection-geometries>`__.
         Keyword arguments `det_spacing` and `angles` should be specified
-        to use the former, and keyword argument `vectors` should be
-        specified to use the latter. These options are mutually exclusive.
+        to use the "parallel3d" geometry, and keyword argument `vectors`
+        should be specified to use the "parallel3d_vec" geometry. These
+        options are mutually exclusive.
 
         Args:
             input_shape: Shape of the input array.

From 05624f069124fce8f1163207d69d6a0c89ef1c5f Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 19:39:04 -0600
Subject: [PATCH 05/15] Improve 3d astra docs

---
 scico/linop/xray/astra.py | 41 ++++++++++++++++++++++++++++++---------
 1 file changed, 32 insertions(+), 9 deletions(-)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index 920b6a303..571928906 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -230,8 +230,8 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
     The `3D geometries <https://astra-toolbox.com/docs/geom3d.html#projection-geometries>`__
     "parallel3d" and "parallel3d_vec" are supported by this interface.
 
-    The reconstruction volume is fixed with respect to the coordinate
-    system, with the volume centred at the origin, as illustrated below:
+    The volume is fixed with respect to the coordinate system, centered
+    at the origin, as illustrated below:
 
     .. plot:: pyfigures/xray_3d_vol.py
        :align: center
@@ -241,7 +241,9 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
     The voxels sides have unit length (in arbitrary units), which defines
     the scale for all other dimensions in the source-volume-detector
     configuration. Geometry axes `x`, `y`, and `z` correspond to volume
-    array axes 0, 1, and 2 respectively.
+    array axes 0, 1, and 2 respectively. The projected array axes 0, 1,
+    and 2 correspond respectively to detector rows, views, and detector
+    columns.
 
     In the "parallel3d" case, the source and detector rotate clockwise
     about the `z` axis in the `x`-`y` plane, as illustrated below:
@@ -250,12 +252,21 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
        :align: center
        :include-source: False
        :show-source-link: False
-       :caption: The red arrows indicate the direction of the beam towards
+       :caption: Red arrows indicate the direction of the beam towards
           the detector (orange) and the arrows parallel to the detector
           indicate the direction of increasing pixel indices.
 
+    In this case the `z` axis is in the same direction as the
+    vertical/row axis of the detector and its projection corresponds to
+    a vertical line in the center of the horizontal/column detector axis.
+    Note that the view images must be displayed with the origin at the
+    bottom left (i.e. vertically inverted from the top left origin image
+    indexing convention) in order for the projections to correspond to
+    the positive up/negative down orientation of the `z` axis in the
+    figures here.
+
     In the "parallel3d_vec" case, each view is determined by the following
-    vectors
+    vectors:
 
     .. list-table:: View definition vectors
        :widths: 10 90
@@ -265,18 +276,30 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
        * - :math:`\mb{d}`
          - Center of the detector
        * - :math:`\mb{u}`
-         - Vector from detector pixel (0,0) to (0,1)
+         - Vector from detector pixel (0,0) to (0,1) (direction of
+           increasing detector column index)
        * - :math:`\mb{v}`
-         - Vector from detector pixel (0,0) to (1,0)
+         - Vector from detector pixel (0,0) to (1,0) (direction of
+           increasing detector row index)
 
     .. plot:: pyfigures/xray_3d_vec.py
        :align: center
        :include-source: False
        :show-source-link: False
 
-    which are concatenated into a single row vector
+    Vector :math:`\mb{r}` is not illustrated to avoid cluttering the
+    figure, but will typically be directed toward the center of the
+    detector. In practice, since the volume-detector distance does not
+    have a geometric effect for a parallel-beam configuration,
+    :math:`\mb{d}` may be set to the zero vector. Note that the view
+    images must be displayed with the origin at the bottom left (i.e.
+    vertically inverted from the top left origin image indexing
+    convention) in order for the row indexing of the projections to
+    correspond to the direction of :math:`\mb{v}` in the figure.
+
+    These vectors are concatenated into a single row vector
     :math:`(\mb{r}, \mb{d}, \mb{u}, \mb{v})` to form the full
-    geometry specification for a single view. Multiple such
+    geometry specification for a single view, and multiple such
     row vectors are stacked to specify the geometry for a set
     of views.
     """

From 7e5fb0f12e74b56be172065138df30e559c1a74f Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 19:41:09 -0600
Subject: [PATCH 06/15] Improve 3d astra docs

---
 scico/linop/xray/astra.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index 571928906..e4d13ca60 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -289,13 +289,14 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
 
     Vector :math:`\mb{r}` is not illustrated to avoid cluttering the
     figure, but will typically be directed toward the center of the
-    detector. In practice, since the volume-detector distance does not
-    have a geometric effect for a parallel-beam configuration,
-    :math:`\mb{d}` may be set to the zero vector. Note that the view
-    images must be displayed with the origin at the bottom left (i.e.
-    vertically inverted from the top left origin image indexing
-    convention) in order for the row indexing of the projections to
-    correspond to the direction of :math:`\mb{v}` in the figure.
+    detector (i.e. in the direction of :math:`\mb{d}` in the figure.) In
+    practice, since the volume-detector distance does not have a
+    geometric effect for a parallel-beam configuration, :math:`\mb{d}`
+    may be set to the zero vector. Note that the view images must be
+    displayed with the origin at the bottom left (i.e. vertically
+    inverted from the top left origin image indexing convention) in order
+    for the row indexing of the projections to correspond to the
+    direction of :math:`\mb{v}` in the figure.
 
     These vectors are concatenated into a single row vector
     :math:`(\mb{r}, \mb{d}, \mb{u}, \mb{v})` to form the full

From 85e721969d3d603d8ab044a98c824657ae3b93a4 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 15 Jul 2024 19:51:54 -0600
Subject: [PATCH 07/15] Reduce space around figures

---
 docs/source/pyfigures/xray_3d_vec.py | 1 +
 docs/source/pyfigures/xray_3d_vol.py | 1 +
 2 files changed, 2 insertions(+)

diff --git a/docs/source/pyfigures/xray_3d_vec.py b/docs/source/pyfigures/xray_3d_vec.py
index af26c3645..ff94c24f3 100644
--- a/docs/source/pyfigures/xray_3d_vec.py
+++ b/docs/source/pyfigures/xray_3d_vec.py
@@ -77,4 +77,5 @@ def do_3d_projection(self, renderer=None):
 ax.text(*v_txtpos, r"$\mathbf{v}$", fontsize=12)
 
 fig.tight_layout()
+fig.subplots_adjust(-0.1, -0.06, 1, 1)
 fig.show()
diff --git a/docs/source/pyfigures/xray_3d_vol.py b/docs/source/pyfigures/xray_3d_vol.py
index 2811933a6..78ee46db7 100644
--- a/docs/source/pyfigures/xray_3d_vol.py
+++ b/docs/source/pyfigures/xray_3d_vol.py
@@ -90,4 +90,5 @@ def do_3d_projection(self, renderer=None):
 ax.plot(bx, by, bz, lw=2, ls="--", color="blue")
 
 fig.tight_layout()
+fig.subplots_adjust(-0.1, -0.1, 1, 1.07)
 fig.show()

From 96acc70fd516959ebfc59392c3952f5a3771e1cf Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Tue, 16 Jul 2024 09:08:25 -0600
Subject: [PATCH 08/15] Improve figure caption style

---
 docs/source/_static/scico.css | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/docs/source/_static/scico.css b/docs/source/_static/scico.css
index eed587dc6..ea0618c0d 100644
--- a/docs/source/_static/scico.css
+++ b/docs/source/_static/scico.css
@@ -103,3 +103,15 @@ div.doctest.highlight-default {
 [data-theme=light] dl.py.property > dt {
   border-radius: 4px;
 }
+
+
+/* Style for figure captions */
+
+div.figure p.caption span.caption-text,
+figcaption span.caption-text {
+  font-size: var(--font-size--small);
+  margin-left: 5%;
+  margin-right: 5%;
+  display: inline-block;
+  text-align: justify;
+}

From f2150f5206af78f79b7cb7fea252def36442faed Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Tue, 16 Jul 2024 09:47:05 -0600
Subject: [PATCH 09/15] Improve 3d astra docs

---
 scico/linop/xray/astra.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index e4d13ca60..f994900d8 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -282,6 +282,9 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
          - Vector from detector pixel (0,0) to (1,0) (direction of
            increasing detector row index)
 
+    Note that the components of these vectors are in `x`, `y`, `z` order,
+    not the `z`, `y`, `x` order of the volume axes.
+
     .. plot:: pyfigures/xray_3d_vec.py
        :align: center
        :include-source: False

From bfa36daa3d9923ff0f87737888e289c96bbb98c4 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Fri, 19 Jul 2024 07:09:22 -0600
Subject: [PATCH 10/15] Typo fix

---
 scico/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scico/_version.py b/scico/_version.py
index cdfd89aa6..86c7b18fc 100644
--- a/scico/_version.py
+++ b/scico/_version.py
@@ -95,7 +95,7 @@ def package_version(split: bool = False) -> Union[str, Tuple[str, str]]:  # prag
     Returns:
         Package version string or tuple of strings.
     """
-    version = init_variable_assign_value("__version_")
+    version = init_variable_assign_value("__version__")
     # don't extend purely numeric version numbers, possibly ending with post<n>
     if re.match(r"^[0-9\.]+(post[0-9]+)?$", version):
         git_hash = None

From ca7251461a70d3a8411cba6e131993e66cf0264e Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Fri, 19 Jul 2024 11:39:35 -0600
Subject: [PATCH 11/15] Move matplotlib figures from submodule

---
 docs/source/pyfigures/cylindgrad.py   | 46 +++++++++++++++
 docs/source/pyfigures/polargrad.py    | 35 +++++++++++
 docs/source/pyfigures/spheregrad.py   | 47 +++++++++++++++
 docs/source/pyfigures/xray_2d_geom.py | 83 +++++++++++++++++++++++++++
 4 files changed, 211 insertions(+)
 create mode 100644 docs/source/pyfigures/cylindgrad.py
 create mode 100644 docs/source/pyfigures/polargrad.py
 create mode 100644 docs/source/pyfigures/spheregrad.py
 create mode 100644 docs/source/pyfigures/xray_2d_geom.py

diff --git a/docs/source/pyfigures/cylindgrad.py b/docs/source/pyfigures/cylindgrad.py
new file mode 100644
index 000000000..9b26220ac
--- /dev/null
+++ b/docs/source/pyfigures/cylindgrad.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+import scico.linop as scl
+from scico import plot
+
+input_shape = (7, 7, 7)
+centre = (np.array(input_shape) - 1) / 2
+end = np.array(input_shape) - centre
+g0, g1, g2 = np.mgrid[-centre[0] : end[0], -centre[1] : end[1], -centre[2] : end[2]]
+
+cg = scl.CylindricalGradient(input_shape=input_shape)
+
+ang = cg.coord[0]
+rad = cg.coord[1]
+axi = cg.coord[2]
+
+theta = np.arctan2(g0, g1)
+clr = theta
+# See https://stackoverflow.com/a/49888126
+clr = (clr.ravel() - clr.min()) / np.ptp(clr)
+clr = np.concatenate((clr, np.repeat(clr, 2)))
+clr = plot.plt.cm.plasma(clr)
+
+plot.plt.rcParams["savefig.transparent"] = True
+
+fig = plot.plt.figure(figsize=(20, 6))
+ax = fig.add_subplot(1, 3, 1, projection="3d")
+ax.quiver(g0, g1, g2, ang[0], ang[1], ang[2], colors=clr, length=0.9)
+ax.set_title("Angular local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+ax = fig.add_subplot(1, 3, 2, projection="3d")
+ax.quiver(g0, g1, g2, rad[0], rad[1], rad[2], colors=clr, length=0.9)
+ax.set_title("Radial local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+ax = fig.add_subplot(1, 3, 3, projection="3d")
+ax.quiver(g0, g1, g2, axi[0], axi[1], axi[2], colors=clr[0], length=0.9)
+ax.set_title("Axial local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+fig.tight_layout()
+fig.show()
diff --git a/docs/source/pyfigures/polargrad.py b/docs/source/pyfigures/polargrad.py
new file mode 100644
index 000000000..23da202f2
--- /dev/null
+++ b/docs/source/pyfigures/polargrad.py
@@ -0,0 +1,35 @@
+import numpy as np
+
+import scico.linop as scl
+from scico import plot
+
+input_shape = (21, 21)
+centre = (np.array(input_shape) - 1) / 2
+end = np.array(input_shape) - centre
+g0, g1 = np.mgrid[-centre[0] : end[0], -centre[1] : end[1]]
+
+pg = scl.PolarGradient(input_shape=input_shape)
+
+ang = pg.coord[0]
+rad = pg.coord[1]
+
+clr = (np.arctan2(ang[1], ang[0]) + np.pi) / (2 * np.pi)
+
+plot.plt.rcParams["image.cmap"] = "plasma"
+plot.plt.rcParams["savefig.transparent"] = True
+
+fig, ax = plot.plt.subplots(nrows=1, ncols=2, figsize=(13, 6))
+ax[0].quiver(g0, g1, ang[0], ang[1], clr)
+ax[0].set_title("Angular local coordinate axis")
+ax[0].set_xlabel("$x$")
+ax[0].set_ylabel("$y$")
+ax[0].xaxis.set_ticks((-10, -5, 0, 5, 10))
+ax[0].yaxis.set_ticks((-10, -5, 0, 5, 10))
+ax[1].quiver(g0, g1, rad[0], rad[1], clr)
+ax[1].set_title("Radial local coordinate axis")
+ax[1].set_xlabel("$x$")
+ax[1].set_ylabel("$y$")
+ax[1].xaxis.set_ticks((-10, -5, 0, 5, 10))
+ax[1].yaxis.set_ticks((-10, -5, 0, 5, 10))
+fig.tight_layout()
+fig.show()
diff --git a/docs/source/pyfigures/spheregrad.py b/docs/source/pyfigures/spheregrad.py
new file mode 100644
index 000000000..ea2149d92
--- /dev/null
+++ b/docs/source/pyfigures/spheregrad.py
@@ -0,0 +1,47 @@
+import numpy as np
+
+import scico.linop as scl
+from scico import plot
+
+input_shape = (7, 7, 7)
+centre = (np.array(input_shape) - 1) / 2
+end = np.array(input_shape) - centre
+g0, g1, g2 = np.mgrid[-centre[0] : end[0], -centre[1] : end[1], -centre[2] : end[2]]
+
+sg = scl.SphericalGradient(input_shape=input_shape)
+
+azi = sg.coord[0]
+pol = sg.coord[1]
+rad = sg.coord[2]
+
+theta = np.arctan2(g0, g1)
+phi = np.arctan2(np.sqrt(g0**2 + g1**2), g2)
+clr = theta * phi
+# See https://stackoverflow.com/a/49888126
+clr = (clr.ravel() - clr.min()) / np.ptp(clr)
+clr = np.concatenate((clr, np.repeat(clr, 2)))
+clr = plot.plt.cm.plasma(clr)
+
+plot.plt.rcParams["savefig.transparent"] = True
+
+fig = plot.plt.figure(figsize=(20, 6))
+ax = fig.add_subplot(1, 3, 1, projection="3d")
+ax.quiver(g0, g1, g2, azi[0], azi[1], azi[2], colors=clr, length=0.9)
+ax.set_title("Azimuthal local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+ax = fig.add_subplot(1, 3, 2, projection="3d")
+ax.quiver(g0, g1, g2, pol[0], pol[1], pol[2], colors=clr, length=0.9)
+ax.set_title("Polar local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+ax = fig.add_subplot(1, 3, 3, projection="3d")
+ax.quiver(g0, g1, g2, rad[0], rad[1], rad[2], colors=clr, length=0.9)
+ax.set_title("Radial local coordinate axis")
+ax.set_xlabel("$x$")
+ax.set_ylabel("$y$")
+ax.set_zlabel("$z$")
+fig.tight_layout()
+fig.show()
diff --git a/docs/source/pyfigures/xray_2d_geom.py b/docs/source/pyfigures/xray_2d_geom.py
new file mode 100644
index 000000000..c1714851f
--- /dev/null
+++ b/docs/source/pyfigures/xray_2d_geom.py
@@ -0,0 +1,83 @@
+import numpy as np
+
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+
+c = 1.0 / np.sqrt(2.0)
+e = 1e-2
+style = "Simple, tail_width=0.5, head_width=4, head_length=8"
+fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(21, 7))
+for n in range(3):
+    ax[n].set_aspect(1.0)
+    ax[n].set_xlim(-1.1, 1.1)
+    ax[n].set_ylim(-1.1, 1.1)
+    ax[n].set_xticks(np.linspace(-1.0, 1.0, 5))
+    ax[n].set_yticks(np.linspace(-1.0, 1.0, 5))
+    ax[n].tick_params(axis="x", labelsize=12)
+    ax[n].tick_params(axis="y", labelsize=12)
+    ax[n].set_xlabel("axis 1", fontsize=14)
+    ax[n].set_ylabel("axis 0", fontsize=14)
+
+# scico
+plist = [
+    patches.FancyArrowPatch((-1.0, 0.0), (-0.5, 0.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((-c, -c), (-c / 2.0, -c / 2.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch(
+        (
+            0.0,
+            -1.0,
+        ),
+        (0.0, -0.5),
+        arrowstyle=style,
+        color="r",
+    ),
+    patches.Arc((0.0, 0.0), 2.0, 2.0, theta1=180, theta2=-45.0, color="b", ls="dotted"),
+    patches.FancyArrowPatch((c - e, -c - e), (c + e, -c + e), arrowstyle=style, color="b"),
+]
+for p in plist:
+    ax[0].add_patch(p)
+ax[0].text(-0.88, 0.02, r"$\theta=0$", color="r", fontsize=14)
+ax[0].text(-3 * c / 4 - 0.01, -3 * c / 4 - 0.1, r"$\theta=\frac{\pi}{4}$", color="r", fontsize=14)
+ax[0].text(0.03, -0.8, r"$\theta=\frac{\pi}{2}$", color="r", fontsize=14)
+ax[0].set_title("scico", fontsize=14)
+
+# astra
+plist = [
+    patches.FancyArrowPatch((0.0, -1.0), (0.0, -0.5), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((c, -c), (c / 2.0, -c / 2.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((1.0, 0.0), (0.5, 0.0), arrowstyle=style, color="r"),
+    patches.Arc((0.0, 0.0), 2.0, 2.0, theta1=-90, theta2=45.0, color="b", ls="dotted"),
+    patches.FancyArrowPatch((c + e, c - e), (c - e, c + e), arrowstyle=style, color="b"),
+]
+for p in plist:
+    ax[1].add_patch(p)
+ax[1].text(0.02, -0.75, r"$\theta=0$", color="r", fontsize=14)
+ax[1].text(3 * c / 4 + 0.01, -3 * c / 4 + 0.01, r"$\theta=\frac{\pi}{4}$", color="r", fontsize=14)
+ax[1].text(0.65, 0.05, r"$\theta=\frac{\pi}{2}$", color="r", fontsize=14)
+ax[1].set_title("astra", fontsize=14)
+
+# svmbir
+plist = [
+    patches.FancyArrowPatch((-1.0, 0.0), (-0.5, 0.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch((-c, c), (-c / 2.0, c / 2.0), arrowstyle=style, color="r"),
+    patches.FancyArrowPatch(
+        (
+            0.0,
+            1.0,
+        ),
+        (0.0, 0.5),
+        arrowstyle=style,
+        color="r",
+    ),
+    patches.Arc((0.0, 0.0), 2.0, 2.0, theta1=45, theta2=180, color="b", ls="dotted"),
+    patches.FancyArrowPatch((c - e, c + e), (c + e, c - e), arrowstyle=style, color="b"),
+]
+for p in plist:
+    ax[2].add_patch(p)
+ax[2].text(-0.88, 0.02, r"$\theta=0$", color="r", fontsize=14)
+ax[2].text(-3 * c / 4 + 0.01, 3 * c / 4 + 0.01, r"$\theta=\frac{\pi}{4}$", color="r", fontsize=14)
+ax[2].text(0.03, 0.75, r"$\theta=\frac{\pi}{2}$", color="r", fontsize=14)
+ax[2].set_title("svmbir", fontsize=14)
+
+fig.tight_layout()
+fig.show()

From 0e97c5874ddc3a9f493664b91f9fe4c45eeb8fd9 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Fri, 19 Jul 2024 11:39:54 -0600
Subject: [PATCH 12/15] Move matplotlib figures from submodule

---
 data                         | 2 +-
 scico/linop/_grad.py         | 6 +++---
 scico/linop/xray/__init__.py | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/data b/data
index 3bc228c83..ece865c38 160000
--- a/data
+++ b/data
@@ -1 +1 @@
-Subproject commit 3bc228c83d14b7bb5c3c64df7f94112c836ca8b7
+Subproject commit ece865c38bef6cc414c25390994e759ccca7e12e
diff --git a/scico/linop/_grad.py b/scico/linop/_grad.py
index 416cc125d..67f9e3f7c 100644
--- a/scico/linop/_grad.py
+++ b/scico/linop/_grad.py
@@ -185,7 +185,7 @@ class PolarGradient(ProjectedGradient):
     directions, as described in :cite:`hossein-2024-total`. Local
     coordinate axes are illustrated in the figure below.
 
-    .. plot:: figures/polargrad.py
+    .. plot:: pyfigures/polargrad.py
        :align: center
        :include-source: False
        :show-source-link: False
@@ -283,7 +283,7 @@ class CylindricalGradient(ProjectedGradient):
     described in :cite:`hossein-2024-total`. The local coordinate axes
     are illustrated in the figure below.
 
-    .. plot:: figures/cylindgrad.py
+    .. plot:: pyfigures/cylindgrad.py
        :align: center
        :include-source: False
        :show-source-link: False
@@ -410,7 +410,7 @@ class SphericalGradient(ProjectedGradient):
     the approach described in :cite:`hossein-2024-total`. The local
     coordinate axes are illustrated in the figure below.
 
-    .. plot:: figures/spheregrad.py
+    .. plot:: pyfigures/spheregrad.py
        :align: center
        :include-source: False
        :show-source-link: False
diff --git a/scico/linop/xray/__init__.py b/scico/linop/xray/__init__.py
index cb47c9ec6..b909f5855 100644
--- a/scico/linop/xray/__init__.py
+++ b/scico/linop/xray/__init__.py
@@ -22,7 +22,7 @@
 these transforms uses a different convention for view angle directions,
 as illustrated in the figure below.
 
-.. plot:: figures/xray_2d_geom.py
+.. plot:: pyfigures/xray_2d_geom.py
    :align: center
    :include-source: False
    :show-source-link: False

From d7811dc4bc38123072baaa5cd80fd3af74534675 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Fri, 19 Jul 2024 19:04:59 -0600
Subject: [PATCH 13/15] Add test

---
 scico/test/linop/xray/test_astra.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/scico/test/linop/xray/test_astra.py b/scico/test/linop/xray/test_astra.py
index de8302331..5ba29123f 100644
--- a/scico/test/linop/xray/test_astra.py
+++ b/scico/test/linop/xray/test_astra.py
@@ -167,3 +167,8 @@ def test_angle_to_vector():
     det_spacing = [0.9, 1.5]
     vectors = angle_to_vector(det_spacing, angles)
     assert vectors.shape == (angles.size, 12)
+
+
+def test_ensure_writeable():
+    assert isinstance(_ensure_writeable(np.ones((2, 1))), np.ndarray)
+    assert isinstance(_ensure_writeable(snp.ones((2, 1))), np.ndarray)

From 286e36bf15f80d3bbbf36915ce09f3e4b716c6db Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 22 Jul 2024 13:44:59 -0600
Subject: [PATCH 14/15] Address PR review comments

---
 scico/linop/xray/astra.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index f994900d8..294b73244 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -240,7 +240,7 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
 
     The voxels sides have unit length (in arbitrary units), which defines
     the scale for all other dimensions in the source-volume-detector
-    configuration. Geometry axes `x`, `y`, and `z` correspond to volume
+    configuration. Geometry axes `z`, `y`, and `x` correspond to volume
     array axes 0, 1, and 2 respectively. The projected array axes 0, 1,
     and 2 correspond respectively to detector rows, views, and detector
     columns.
@@ -292,11 +292,12 @@ class XRayTransform3D(LinearOperator):  # pragma: no cover
 
     Vector :math:`\mb{r}` is not illustrated to avoid cluttering the
     figure, but will typically be directed toward the center of the
-    detector (i.e. in the direction of :math:`\mb{d}` in the figure.) In
-    practice, since the volume-detector distance does not have a
-    geometric effect for a parallel-beam configuration, :math:`\mb{d}`
-    may be set to the zero vector. Note that the view images must be
-    displayed with the origin at the bottom left (i.e. vertically
+    detector (i.e. in the direction of :math:`\mb{d}` in the figure.)
+    Since the volume-detector distance does not have a geometric effect
+    for a parallel-beam configuration, :math:`\mb{d}` may be set to the
+    zero vector when the detector and beam centers coincide (e.g., as in
+    the case of the "parallel3d" geometry). Note that the view images
+    must be displayed with the origin at the bottom left (i.e. vertically
     inverted from the top left origin image indexing convention) in order
     for the row indexing of the projections to correspond to the
     direction of :math:`\mb{v}` in the figure.

From 87c8d39afb6c977a3c9781218b9b9e333473bac0 Mon Sep 17 00:00:00 2001
From: Brendt Wohlberg <brendt@ieee.org>
Date: Mon, 22 Jul 2024 14:19:08 -0600
Subject: [PATCH 15/15] Update submodule

---
 data | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/data b/data
index ece865c38..40ebe0ab4 160000
--- a/data
+++ b/data
@@ -1 +1 @@
-Subproject commit ece865c38bef6cc414c25390994e759ccca7e12e
+Subproject commit 40ebe0ab4e893620339f928b33b8744dd9c111a7