diff --git a/data b/data
index 3bc228c83..40ebe0ab4 160000
--- a/data
+++ b/data
@@ -1 +1 @@
-Subproject commit 3bc228c83d14b7bb5c3c64df7f94112c836ca8b7
+Subproject commit 40ebe0ab4e893620339f928b33b8744dd9c111a7
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;
+}
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()
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..ff94c24f3
--- /dev/null
+++ b/docs/source/pyfigures/xray_3d_vec.py
@@ -0,0 +1,81 @@
+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.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
new file mode 100644
index 000000000..78ee46db7
--- /dev/null
+++ b/docs/source/pyfigures/xray_3d_vol.py
@@ -0,0 +1,94 @@
+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.subplots_adjust(-0.1, -0.1, 1, 1.07)
+fig.show()
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
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 3a5dc7d4a..b909f5855 100644
--- a/scico/linop/xray/__init__.py
+++ b/scico/linop/xray/__init__.py
@@ -22,12 +22,12 @@
 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
    :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
diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index aa530a6ef..294b73244 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -227,6 +227,86 @@ 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#projection-geometries>`__
+    "parallel3d" and "parallel3d_vec" are supported by this interface.
+
+    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
+       :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 `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.
+
+    In the "parallel3d" case, the source and detector rotate clockwise
+    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: 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:
+
+    .. 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) (direction of
+           increasing detector column index)
+       * - :math:`\mb{v}`
+         - 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
+       :show-source-link: False
+
+    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.)
+    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.
+
+    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, and multiple such
+    row vectors are stacked to specify the geometry for a set
+    of views.
     """
 
     def __init__(
@@ -238,11 +318,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.
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)