diff --git a/.gitignore b/.gitignore
index 21d5538cc..0396a82b8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -123,6 +123,9 @@ venv.bak/
 # Rope project settings
 .ropeproject
 
+# VS Code settings
+.vscode/
+
 # mkdocs documentation
 /site
 
diff --git a/CHANGES.rst b/CHANGES.rst
index bdcfd1295..8a7864843 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,6 +7,7 @@ Version 0.0.6   (unreleased)
 ----------------------------
 
 • Significant changes to ``linop.xray.astra`` API.
+• New integrated 3D X-ray transform via ``linop.xray.XRayTransform3D``.
 • New functional ``functional.IsotropicTVNorm`` and faster implementation
   of ``functional.AnisotropicTVNorm``.
 • New linear operators ``linop.ProjectedGradient``, ``linop.PolarGradient``,
diff --git a/data b/data
index 40ebe0ab4..790d0316d 160000
--- a/data
+++ b/data
@@ -1 +1 @@
-Subproject commit 40ebe0ab4e893620339f928b33b8744dd9c111a7
+Subproject commit 790d0316dd1fc77cc0ae0258d3ef20d158e3f258
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index 13572f3d7..58ba847f0 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -37,7 +37,9 @@ Computed Tomography
    examples/ct_astra_modl_train_foam2
    examples/ct_astra_odp_train_foam2
    examples/ct_astra_unet_train_foam2
-   examples/ct_projector_comparison
+   examples/ct_projector_comparison_2d
+   examples/ct_projector_comparison_3d
+   examples/ct_multi_cs_tv_admm
    examples/ct_multi_tv_admm
 
 Deconvolution
diff --git a/examples/scriptcheck.sh b/examples/scriptcheck.sh
index 498668585..38da4fd6e 100755
--- a/examples/scriptcheck.sh
+++ b/examples/scriptcheck.sh
@@ -90,6 +90,10 @@ for f in $SCRIPTPATH/scripts/*.py; do
         printf "%s\n" skipped
         continue
     fi
+    if [ $SKIP_GPU -eq 1 ] && grep -q 'ct_projector_comparison_3d' <<< $f; then
+        printf "%s\n" skipped
+        continue
+    fi
 
     # Create temporary copy of script with all algorithm maxiter values set
     # to small number and final input statements commented out.
diff --git a/examples/scripts/README.rst b/examples/scripts/README.rst
index 481e75ed0..446186a76 100644
--- a/examples/scripts/README.rst
+++ b/examples/scripts/README.rst
@@ -39,8 +39,10 @@ Computed Tomography
       CT Training and Reconstructions with ODP
    `ct_astra_unet_train_foam2.py <ct_astra_unet_train_foam2.py>`_
       CT Training and Reconstructions with UNet
-   `ct_projector_comparison.py <ct_projector_comparison.py>`_
-      X-ray Transform Comparison
+   `ct_projector_comparison_2d.py <ct_projector_comparison_2d.py>`_
+      2D X-ray Transform Comparison
+   `ct_projector_comparison_3d.py <ct_projector_comparison_3d.py>`_
+      3D X-ray Transform Comparison
    `ct_multi_tv_admm.py <ct_multi_tv_admm.py>`_
       TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors)
 
diff --git a/examples/scripts/ct_multi_tv_admm.py b/examples/scripts/ct_multi_tv_admm.py
index 87fc86865..957116799 100644
--- a/examples/scripts/ct_multi_tv_admm.py
+++ b/examples/scripts/ct_multi_tv_admm.py
@@ -27,7 +27,7 @@
 
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
-from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir
+from scico.linop.xray import XRayTransform2D, astra, svmbir
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -54,9 +54,7 @@
     "svmbir": svmbir.XRayTransform(
         x_gt.shape, 2 * np.pi - angles, det_count, delta_pixel=1.0, delta_channel=det_spacing
     ),  # svmbir
-    "scico": XRayTransform(
-        Parallel2dProjector((N, N), angles, det_count=det_count, dx=1 / det_spacing)
-    ),  # scico
+    "scico": XRayTransform2D((N, N), angles, det_count=det_count, dx=1 / det_spacing),  # scico
 }
 
 
diff --git a/examples/scripts/ct_projector_comparison.py b/examples/scripts/ct_projector_comparison_2d.py
similarity index 94%
rename from examples/scripts/ct_projector_comparison.py
rename to examples/scripts/ct_projector_comparison_2d.py
index b0e3b68b9..2e5d02d3f 100644
--- a/examples/scripts/ct_projector_comparison.py
+++ b/examples/scripts/ct_projector_comparison_2d.py
@@ -6,10 +6,10 @@
 
 
 r"""
-X-ray Transform Comparison
-==========================
+2D X-ray Transform Comparison
+=============================
 
-This example compares SCICO's native X-ray transform algorithm
+This example compares SCICO's native 2D X-ray transform algorithm
 to that of the ASTRA toolbox.
 """
 
@@ -22,7 +22,7 @@
 
 import scico.linop.xray.astra as astra
 from scico import plot
-from scico.linop import Parallel2dProjector, XRayTransform
+from scico.linop.xray import XRayTransform2D
 from scico.util import Timer
 
 """
@@ -46,7 +46,7 @@
 
 projectors = {}
 timer.start("scico_init")
-projectors["scico"] = XRayTransform(Parallel2dProjector((N, N), angles))
+projectors["scico"] = XRayTransform2D((N, N), angles)
 timer.stop("scico_init")
 
 timer.start("astra_init")
diff --git a/examples/scripts/ct_projector_comparison_3d.py b/examples/scripts/ct_projector_comparison_3d.py
new file mode 100644
index 000000000..222715d99
--- /dev/null
+++ b/examples/scripts/ct_projector_comparison_3d.py
@@ -0,0 +1,200 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# This file is part of the SCICO package. Details of the copyright
+# and user license can be found in the 'LICENSE.txt' file distributed
+# with the package.
+
+
+r"""
+3D X-ray Transform Comparison
+=============================
+
+This example shows how to define a SCICO native 3D X-ray transform using
+ASTRA toolbox conventions and vice versa.
+"""
+
+import numpy as np
+
+import jax
+import jax.numpy as jnp
+
+import scico.linop.xray.astra as astra
+from scico import plot
+from scico.examples import create_block_phantom
+from scico.linop.xray import XRayTransform3D
+from scico.util import ContextTimer, Timer
+
+"""
+Create a ground truth image and set detector dimensions.
+"""
+N = 64
+# use rectangular volume to check whether axes are handled correctly
+in_shape = (N + 1, N + 2, N + 3)
+x = create_block_phantom(in_shape)
+x = jnp.array(x)
+
+# use rectangular detector to check whether axes are handled correctly
+out_shape = (N, N + 1)
+
+
+"""
+Set up SCICO projection.
+"""
+num_angles = 3
+
+
+rot_X = 90.0 - 16.0
+rot_Y = np.linspace(0, 180, num_angles, endpoint=False)
+angles = np.stack(np.broadcast_arrays(rot_X, rot_Y), axis=-1)
+matrices = XRayTransform3D.matrices_from_euler_angles(
+    in_shape, out_shape, "XY", angles, degrees=True
+)
+
+"""
+Specify geometry using SCICO conventions and project.
+"""
+num_repeats = 3
+
+timer_scico = Timer()
+with ContextTimer(timer_scico, "init"):
+    H_scico = XRayTransform3D(in_shape, matrices, out_shape)
+
+with ContextTimer(timer_scico, "first_fwd"):
+    y_scico = H_scico @ x
+    jax.block_until_ready(y_scico)
+
+with ContextTimer(timer_scico, "avg_fwd"):
+    for _ in range(num_repeats):
+        y_scico = H_scico @ x
+        jax.block_until_ready(y_scico)
+timer_scico.td["avg_fwd"] /= num_repeats
+
+with ContextTimer(timer_scico, "first_back"):
+    HTy_scico = H_scico.T @ y_scico
+
+with ContextTimer(timer_scico, "avg_back"):
+    for _ in range(num_repeats):
+        HTy_scico = H_scico.T @ y_scico
+        jax.block_until_ready(HTy_scico)
+timer_scico.td["avg_back"] /= num_repeats
+
+
+"""
+Convert SCICO geometry to ASTRA and project.
+"""
+
+vectors_from_scico = astra.convert_from_scico_geometry(in_shape, matrices, out_shape)
+
+timer_astra = Timer()
+with ContextTimer(timer_astra, "init"):
+    H_astra_from_scico = astra.XRayTransform3D(
+        input_shape=in_shape, det_count=out_shape, vectors=vectors_from_scico
+    )
+
+with ContextTimer(timer_astra, "first_fwd"):
+    y_astra_from_scico = H_astra_from_scico @ x
+    jax.block_until_ready(y_astra_from_scico)
+
+with ContextTimer(timer_astra, "avg_fwd"):
+    for _ in range(num_repeats):
+        y_astra_from_scico = H_astra_from_scico @ x
+        jax.block_until_ready(y_astra_from_scico)
+timer_astra.td["avg_fwd"] /= num_repeats
+
+with ContextTimer(timer_astra, "first_back"):
+    HTy_astra_from_scico = H_astra_from_scico.T @ y_astra_from_scico
+
+with ContextTimer(timer_astra, "avg_back"):
+    for _ in range(num_repeats):
+        HTy_astra_from_scico = H_astra_from_scico.T @ y_astra_from_scico
+        jax.block_until_ready(HTy_astra_from_scico)
+timer_astra.td["avg_back"] /= num_repeats
+
+
+"""
+Specify geometry with ASTRA conventions and project.
+"""
+
+angles = np.random.rand(num_angles) * 180  # random projection angles
+det_spacing = [1.0, 1.0]
+vectors = astra.angle_to_vector(det_spacing, angles)
+
+H_astra = astra.XRayTransform3D(input_shape=in_shape, det_count=out_shape, vectors=vectors)
+
+y_astra = H_astra @ x
+HTy_astra = H_astra.T @ y_astra
+
+
+"""
+Convert ASTRA geometry to SCICO and project.
+"""
+
+P_from_astra = astra._astra_to_scico_geometry(H_astra.vol_geom, H_astra.proj_geom)
+H_scico_from_astra = XRayTransform3D(in_shape, P_from_astra, out_shape)
+
+y_scico_from_astra = H_scico_from_astra @ x
+HTy_scico_from_astra = H_scico_from_astra.T @ y_scico_from_astra
+
+
+"""
+Print timing results.
+"""
+print(f"init         astra    {timer_astra.td['init']:.2e} s")
+print(f"init         scico    {timer_scico.td['init']:.2e} s")
+print("")
+for tstr in ("first", "avg"):
+    for dstr in ("fwd", "back"):
+        for timer, pstr in zip((timer_astra, timer_scico), ("astra", "scico")):
+            print(f"{tstr:5s}  {dstr:4s}  {pstr}    {timer.td[tstr + '_' + dstr]:.2e} s")
+        print()
+
+
+"""
+Show projections.
+"""
+fig, ax = plot.subplots(nrows=3, ncols=2, figsize=(8, 10))
+plot.imview(y_scico[0], title="SCICO projections", cbar=None, fig=fig, ax=ax[0, 0])
+plot.imview(y_scico[1], cbar=None, fig=fig, ax=ax[1, 0])
+plot.imview(y_scico[2], cbar=None, fig=fig, ax=ax[2, 0])
+plot.imview(y_astra_from_scico[:, 0], title="ASTRA projections", cbar=None, fig=fig, ax=ax[0, 1])
+plot.imview(y_astra_from_scico[:, 1], cbar=None, fig=fig, ax=ax[1, 1])
+plot.imview(y_astra_from_scico[:, 2], cbar=None, fig=fig, ax=ax[2, 1])
+fig.suptitle("Using SCICO conventions")
+fig.tight_layout()
+fig.show()
+
+fig, ax = plot.subplots(nrows=3, ncols=2, figsize=(8, 10))
+plot.imview(y_scico_from_astra[0], title="SCICO projections", cbar=None, fig=fig, ax=ax[0, 0])
+plot.imview(y_scico_from_astra[1], cbar=None, fig=fig, ax=ax[1, 0])
+plot.imview(y_scico_from_astra[2], cbar=None, fig=fig, ax=ax[2, 0])
+plot.imview(y_astra[:, 0], title="ASTRA projections", cbar=None, fig=fig, ax=ax[0, 1])
+plot.imview(y_astra[:, 1], cbar=None, fig=fig, ax=ax[1, 1])
+plot.imview(y_astra[:, 2], cbar=None, fig=fig, ax=ax[2, 1])
+fig.suptitle("Using ASTRA conventions")
+fig.tight_layout()
+fig.show()
+
+
+"""
+Show back projections.
+"""
+fig, ax = plot.subplots(nrows=1, ncols=2, figsize=(8, 5))
+plot.imview(HTy_scico[N // 2], title="SCICO back projection", cbar=None, fig=fig, ax=ax[0])
+plot.imview(
+    HTy_astra_from_scico[N // 2], title="ASTRA back projection", cbar=None, fig=fig, ax=ax[1]
+)
+fig.suptitle("Using SCICO conventions")
+fig.tight_layout()
+fig.show()
+
+fig, ax = plot.subplots(nrows=1, ncols=2, figsize=(8, 5))
+plot.imview(
+    HTy_scico_from_astra[N // 2], title="SCICO back projection", cbar=None, fig=fig, ax=ax[0]
+)
+plot.imview(HTy_astra[N // 2], title="ASTRA back projection", cbar=None, fig=fig, ax=ax[1])
+fig.suptitle("Using ASTRA conventions")
+fig.tight_layout()
+fig.show()
+
+
+input("\nWaiting for input to close figures and exit")
diff --git a/examples/scripts/ct_tv_admm.py b/examples/scripts/ct_tv_admm.py
index 4b77eeb07..ec48d4eaa 100644
--- a/examples/scripts/ct_tv_admm.py
+++ b/examples/scripts/ct_tv_admm.py
@@ -29,7 +29,7 @@
 
 import scico.numpy as snp
 from scico import functional, linop, loss, metric, plot
-from scico.linop.xray import Parallel2dProjector, XRayTransform
+from scico.linop.xray import XRayTransform2D
 from scico.optimize.admm import ADMM, LinearSubproblemSolver
 from scico.util import device_info
 
@@ -46,7 +46,7 @@
 """
 n_projection = 45  # number of projections
 angles = np.linspace(0, np.pi, n_projection) + np.pi / 2.0  # evenly spaced projection angles
-A = XRayTransform(Parallel2dProjector((N, N), angles))  # CT projection operator
+A = XRayTransform2D((N, N), angles)  # CT projection operator
 y = A @ x_gt  # sinogram
 
 
diff --git a/examples/scripts/index.rst b/examples/scripts/index.rst
index ee4973097..4f05ba2fc 100644
--- a/examples/scripts/index.rst
+++ b/examples/scripts/index.rst
@@ -24,7 +24,8 @@ Computed Tomography
    - ct_astra_modl_train_foam2.py
    - ct_astra_odp_train_foam2.py
    - ct_astra_unet_train_foam2.py
-   - ct_projector_comparison.py
+   - ct_projector_comparison_2d.py
+   - ct_projector_comparison_3d.py
    - ct_multi_tv_admm.py
 
 Deconvolution
diff --git a/scico/examples.py b/scico/examples.py
index 076069562..955b7d5d3 100644
--- a/scico/examples.py
+++ b/scico/examples.py
@@ -12,10 +12,14 @@
 import os
 import tempfile
 import zipfile
+from functools import partial
 from typing import List, Optional, Tuple, Union
 
 import numpy as np
 
+import jax
+import jax.numpy as jnp
+
 import imageio.v3 as iio
 
 import scico.numpy as snp
@@ -518,7 +522,7 @@ def create_conv_sparse_phantom(Nx: int, Nnz: int) -> Tuple[np.ndarray, np.ndarra
 
 
 def create_tangle_phantom(nx: int, ny: int, nz: int) -> snp.Array:
-    """Construct a volume phantom.
+    """Construct a 3D phantom using the tangle function.
 
     Args:
         nx: x-size of output.
@@ -551,6 +555,44 @@ def create_tangle_phantom(nx: int, ny: int, nz: int) -> snp.Array:
     return (values < 2.0).astype(float)
 
 
+@partial(jax.jit, static_argnums=0)
+def create_block_phantom(out_shape: Shape) -> snp.Array:
+    """Construct a blocky 3D phantom.
+
+    Args:
+        out_shape: desired phantom shape.
+
+    Returns:
+        Phantom.
+
+    """
+    # make the phantom at a low resolution
+    low_res = jnp.array(
+        [
+            [
+                [0.0, 0.0, 0.0],
+                [0.0, 0.0, 0.0],
+                [0.0, 0.0, 0.0],
+            ],
+            [
+                [0.0, 1.0, 0.0],
+                [0.0, 1.0, 0.0],
+                [1.0, 1.0, 0.0],
+            ],
+            [
+                [0.0, 1.0, 0.0],
+                [0.0, 0.0, 0.0],
+                [0.0, 0.0, 0.0],
+            ],
+        ]
+    )
+    positions = jnp.stack(
+        jnp.meshgrid(*[jnp.linspace(-0.5, 2.5, s) for s in out_shape], indexing="ij")
+    )
+    indices = jnp.round(positions).astype(int)
+    return low_res[indices[0], indices[1], indices[2]]
+
+
 def spnoise(
     img: Union[np.ndarray, snp.Array], nfrac: float, nmin: float = 0.0, nmax: float = 1.0
 ) -> Union[np.ndarray, snp.Array]:
diff --git a/scico/linop/__init__.py b/scico/linop/__init__.py
index d04104e81..2d6f5cd5d 100644
--- a/scico/linop/__init__.py
+++ b/scico/linop/__init__.py
@@ -25,7 +25,6 @@
 from ._matrix import MatrixOperator
 from ._stack import DiagonalReplicated, DiagonalStack, VerticalStack, linop_over_axes
 from ._util import jacobian, operator_norm, power_iteration, valid_adjoint
-from .xray import Parallel2dProjector, XRayTransform
 
 __all__ = [
     "CircularConvolve",
@@ -51,8 +50,6 @@
     "Sum",
     "Transpose",
     "LinearOperator",
-    "XRayTransform",
-    "Parallel2dProjector",
     "ComposedLinearOperator",
     "linop_from_function",
     "linop_over_axes",
diff --git a/scico/linop/xray/__init__.py b/scico/linop/xray/__init__.py
index 906fb2c50..75c66368b 100644
--- a/scico/linop/xray/__init__.py
+++ b/scico/linop/xray/__init__.py
@@ -15,12 +15,16 @@
 transform that is the appropriate mathematical model for beam attenuation
 based imaging in three or more dimensions.
 
-SCICO includes its own integrated 2D X-ray transform, and also provides
-interfaces to those implemented in the
+SCICO includes its own integrated 2D and 3D X-ray transforms, and also
+provides interfaces to those implemented in the
 `ASTRA toolbox <https://github.com/astra-toolbox/astra-toolbox>`_
-and the `svmbir <https://github.com/cabouman/svmbir>`_ package. Each of
-these transforms uses a different convention for view angle directions,
-as illustrated in the figure below.
+and the `svmbir <https://github.com/cabouman/svmbir>`_ package.
+
+
+**2D Transforms**
+
+The SCICO, ASTRA, and svmbir transforms use different conventions for
+view angle directions, as illustrated in the figure below.
 
 .. plot:: pyfigures/xray_2d_geom.py
    :align: center
@@ -43,13 +47,33 @@
    \theta_{\text{svmbir}} &= 2 \pi - \theta_{\text{scico}} \;.
    \end{aligned}
 
+
+**3D Transforms**
+
+There are more significant differences in the interfaces for the 3D SCICO
+and ASTRA transforms. The SCICO 3D transform :class:`.xray.XRayTransform3D`
+defines the projection geometry in terms of a set of projection matrices,
+while the geometry for the ASTRA 3D transform
+:class:`.astra.XRayTransform3D` may either be specified in terms of a set
+of view angles, or via a more general set of vectors specifying projection
+direction and detector orientation. A number of support functions are
+provided for convering between these conventions.
+
+Note that the SCICO transform is implemented in JAX and can be run on
+both CPU and GPU devices, while the ASTRA transform is implemented in
+CUDA, and can only be run on GPU devices.
 """
 
 import sys
 
-from ._xray import Parallel2dProjector, XRayTransform
+from ._xray import XRayTransform2D, XRayTransform3D
 
 __all__ = [
-    "XRayTransform",
-    "Parallel2dProjector",
+    "XRayTransform2D",
+    "XRayTransform3D",
 ]
+
+
+# Imported items in __all__ appear to originate in top-level xray module
+for name in __all__:
+    getattr(sys.modules[__name__], name).__module__ = __name__
diff --git a/scico/linop/xray/_xray.py b/scico/linop/xray/_xray.py
index 645f5ca49..04bdabe88 100644
--- a/scico/linop/xray/_xray.py
+++ b/scico/linop/xray/_xray.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# Copyright (C) 2023 by SCICO Developers
+# Copyright (C) 2023-2024 by SCICO Developers
 # All rights reserved. BSD 3-clause License.
 # This file is part of the SCICO package. Details of the copyright and
 # user license can be found in the 'LICENSE' file distributed with the
@@ -18,35 +18,15 @@
 import jax.numpy as jnp
 from jax.typing import ArrayLike
 
+import scico.numpy as snp
 from scico.numpy.util import is_scalar_equiv
 from scico.typing import Shape
+from scipy.spatial.transform import Rotation
 
 from .._linop import LinearOperator
 
 
-class XRayTransform(LinearOperator):
-    """X-ray transform operator.
-
-    Wrap an X-ray projector object in a SCICO :class:`LinearOperator`.
-    """
-
-    def __init__(self, projector):
-        r"""
-        Args:
-            projector: instance of an X-ray projector object to wrap,
-                currently the only option is
-                :class:`Parallel2dProjector`
-        """
-        self.projector = projector
-        self._eval = projector.project
-
-        super().__init__(
-            input_shape=projector.im_shape,
-            output_shape=(len(projector.angles), projector.det_count),
-        )
-
-
-class Parallel2dProjector:
+class XRayTransform2D(LinearOperator):
     """Parallel ray, single axis, 2D X-ray projector.
 
     This implementation approximates the projection of each rectangular
@@ -68,7 +48,7 @@ class Parallel2dProjector:
 
     def __init__(
         self,
-        im_shape: Shape,
+        input_shape: Shape,
         angles: ArrayLike,
         x0: Optional[ArrayLike] = None,
         dx: Optional[ArrayLike] = None,
@@ -77,29 +57,29 @@ def __init__(
     ):
         r"""
         Args:
-            im_shape: Shape of input array.
+            input_shape: Shape of input array.
             angles: (num_angles,) array of angles in radians. Viewing an
                 (M, N) array as a matrix with M rows and N columns, an
                 angle of 0 corresponds to summing rows, an angle of pi/2
                 corresponds to summing columns, and an angle of pi/4
                 corresponds to summing along antidiagonals.
             x0: (x, y) position of the corner of the pixel `im[0,0]`. By
-                default, `(-im_shape / 2, -im_shape / 2)`.
+                default, `(-input_shape / 2, -input_shape / 2)`.
             dx: Image pixel side length in x- and y-direction. Should be
                 <= 1.0 in each dimension. By default, [1.0, 1.0].
             y0: Location of the edge of the first detector bin. By
                 default, `-det_count / 2`
             det_count: Number of elements in detector. If ``None``,
-                defaults to the size of the diagonal of `im_shape`.
+                defaults to the size of the diagonal of `input_shape`.
         """
-        self.im_shape = im_shape
+        self.input_shape = input_shape
         self.angles = angles
 
-        self.nx = np.array(im_shape)
+        self.nx = tuple(input_shape)
         if dx is None:
-            dx = np.full(2, np.sqrt(2) / 2)
+            dx = 2 * (np.sqrt(2) / 2,)
         if is_scalar_equiv(dx):
-            dx = np.full(2, dx)
+            dx = 2 * (dx,)
         self.dx = dx
 
         # check projected pixel width assumption
@@ -115,124 +95,369 @@ def __init__(
             )
 
         if x0 is None:
-            x0 = -(self.nx * self.dx) / 2
+            x0 = -(np.array(self.nx) * self.dx) / 2
         self.x0 = x0
 
         if det_count is None:
-            det_count = int(np.ceil(np.linalg.norm(im_shape)))
+            det_count = int(np.ceil(np.linalg.norm(input_shape)))
         self.det_count = det_count
         self.ny = det_count
+        self.output_shape = (len(angles), det_count)
 
         if y0 is None:
             y0 = -self.ny / 2
         self.y0 = y0
         self.dy = 1.0
 
-    def project(self, im):
+        super().__init__(
+            input_shape=self.input_shape,
+            output_shape=self.output_shape,
+            eval_fn=self.project,
+            adj_fn=self.back_project,
+        )
+
+    def project(self, im: ArrayLike) -> snp.Array:
         """Compute X-ray projection."""
-        return _project(im, self.x0, self.dx, self.y0, self.ny, self.angles)
+        return XRayTransform2D._project(im, self.x0, self.dx, self.y0, self.ny, self.angles)
 
-    def back_project(self, y):
+    def back_project(self, y: ArrayLike) -> snp.Array:
         """Compute X-ray back projection"""
-        return _back_project(y, self.x0, self.dx, self.nx, self.y0, self.angles)
-
-
-@partial(jax.jit, static_argnames=["ny"])
-def _project(im, x0, dx, y0, ny, angles):
-    r"""
-    Args:
-        im: Input array, (M, N).
-        x0: (x, y) position of the corner of the pixel im[0,0].
-        dx: Pixel side length in x- and y-direction. Units are such
-            that the detector bins have length 1.0.
-        y0: Location of the edge of the first detector bin.
-        ny: Number of detector bins.
-        angles: (num_angles,) array of angles in radians. Pixels are
-            projected onto units vectors pointing in these directions.
-    """
-    nx = im.shape
-    inds, weights = _calc_weights(x0, dx, nx, angles, y0)
-    # Handle out of bounds indices. In the .at call, inds >= y0 are
-    # ignored, while inds < 0 wrap around. So we set inds < 0 to ny.
-    inds = jnp.where(inds >= 0, inds, ny)
-
-    y = (
-        jnp.zeros((len(angles), ny))
-        .at[jnp.arange(len(angles)).reshape(-1, 1, 1), inds]
-        .add(im * weights)
-    )
-
-    y = y.at[jnp.arange(len(angles)).reshape(-1, 1, 1), inds + 1].add(im * (1 - weights))
-
-    return y
-
-
-@partial(jax.jit, static_argnames=["nx"])
-def _back_project(y, x0, dx, nx, y0, angles):
-    r"""
-    Args:
-        y: Input projection, (num_angles, N).
-        x0: (x, y) position of the corner of the pixel im[0,0].
-        dx: Pixel side length in x- and y-direction. Units are such
-            that the detector bins have length 1.0.
-        nx: Shape of back projection.
-        y0: Location of the edge of the first detector bin.
-        angles: (num_angles,) array of angles in radians. Pixels are
-            projected onto units vectors pointing in these directions.
-    """
-    ny = y.shape[1]
-    inds, weights = _calc_weights(x0, dx, nx, angles, y0)
-    # Handle out of bounds indices. In the .at call, inds >= y0 are
-    # ignored, while inds < 0 wrap around. So we set inds < 0 to ny.
-    inds = jnp.where(inds >= 0, inds, ny)
+        return XRayTransform2D._back_project(y, self.x0, self.dx, self.nx, self.y0, self.angles)
 
-    # the idea: [y[0, inds[0]], y[1, inds[1]], ...]
-    HTy = jnp.sum(y[jnp.arange(len(angles)).reshape(-1, 1, 1), inds] * weights, axis=0)
-    HTy = HTy + jnp.sum(
-        y[jnp.arange(len(angles)).reshape(-1, 1, 1), inds + 1] * (1 - weights), axis=0
-    )
+    @staticmethod
+    @partial(jax.jit, static_argnames=["ny"])
+    def _project(
+        im: ArrayLike, x0: ArrayLike, dx: ArrayLike, y0: float, ny: int, angles: ArrayLike
+    ) -> snp.Array:
+        r"""
+        Args:
+            im: Input array, (M, N).
+            x0: (x, y) position of the corner of the pixel im[0,0].
+            dx: Pixel side length in x- and y-direction. Units are such
+                that the detector bins have length 1.0.
+            y0: Location of the edge of the first detector bin.
+            ny: Number of detector bins.
+            angles: (num_angles,) array of angles in radians. Pixels are
+                projected onto unit vectors pointing in these directions.
+        """
+        nx = im.shape
+        inds, weights = XRayTransform2D._calc_weights(x0, dx, nx, angles, y0)
+        # Handle out of bounds indices. In the .at call, inds >= y0 are
+        # ignored, while inds < 0 wrap around. So we set inds < 0 to ny.
+        inds = jnp.where(inds >= 0, inds, ny)
+
+        y = (
+            jnp.zeros((len(angles), ny))
+            .at[jnp.arange(len(angles)).reshape(-1, 1, 1), inds]
+            .add(im * weights)
+        )
+
+        y = y.at[jnp.arange(len(angles)).reshape(-1, 1, 1), inds + 1].add(im * (1 - weights))
+
+        return y
+
+    @staticmethod
+    @partial(jax.jit, static_argnames=["nx"])
+    def _back_project(
+        y: ArrayLike, x0: ArrayLike, dx: ArrayLike, nx: Shape, y0: float, angles: ArrayLike
+    ) -> ArrayLike:
+        r"""
+        Args:
+            y: Input projection, (num_angles, N).
+            x0: (x, y) position of the corner of the pixel im[0,0].
+            dx: Pixel side length in x- and y-direction. Units are such
+                that the detector bins have length 1.0.
+            nx: Shape of back projection.
+            y0: Location of the edge of the first detector bin.
+            angles: (num_angles,) array of angles in radians. Pixels are
+                projected onto units vectors pointing in these directions.
+        """
+        ny = y.shape[1]
+        inds, weights = XRayTransform2D._calc_weights(x0, dx, nx, angles, y0)
+        # Handle out of bounds indices. In the .at call, inds >= y0 are
+        # ignored, while inds < 0 wrap around. So we set inds < 0 to ny.
+        inds = jnp.where(inds >= 0, inds, ny)
+
+        # the idea: [y[0, inds[0]], y[1, inds[1]], ...]
+        HTy = jnp.sum(y[jnp.arange(len(angles)).reshape(-1, 1, 1), inds] * weights, axis=0)
+        HTy = HTy + jnp.sum(
+            y[jnp.arange(len(angles)).reshape(-1, 1, 1), inds + 1] * (1 - weights), axis=0
+        )
+
+        return HTy
+
+    @staticmethod
+    @partial(jax.jit, static_argnames=["nx"])
+    @partial(jax.vmap, in_axes=(None, None, None, 0, None))
+    def _calc_weights(
+        x0: ArrayLike, dx: ArrayLike, nx: Shape, angle: float, y0: float
+    ) -> snp.Array:
+        """
+
+        Args:
+            x0: Location of the corner of the pixel im[0,0].
+            dx: Pixel side length in x- and y-direction. Units are such
+                that the detector bins have length 1.0.
+            nx: Input image shape.
+            angle: (num_angles,) array of angles in radians. Pixels are
+                projected onto units vectors pointing in these directions.
+                (This argument is `vmap`ed.)
+            y0: Location of the edge of the first detector bin.
+        """
+        u = [jnp.cos(angle), jnp.sin(angle)]
+        Px0 = x0[0] * u[0] + x0[1] * u[1] - y0
+        Pdx = [dx[0] * u[0], dx[1] * u[1]]
+        Pxmin = jnp.min(jnp.array([Px0, Px0 + Pdx[0], Px0 + Pdx[1], Px0 + Pdx[0] + Pdx[1]]))
+
+        Px = (
+            Pxmin
+            + Pdx[0] * jnp.arange(nx[0]).reshape(-1, 1)
+            + Pdx[1] * jnp.arange(nx[1]).reshape(1, -1)
+        )
+
+        # detector bin inds
+        inds = jnp.floor(Px).astype(int)
+
+        # weights
+        Pdx = jnp.array(u) * jnp.array(dx)
+        diag1 = jnp.abs(Pdx[0] + Pdx[1])
+        diag2 = jnp.abs(Pdx[0] - Pdx[1])
+        w = jnp.max(jnp.array([diag1, diag2]))
+        f = jnp.min(jnp.array([diag1, diag2]))
+
+        width = (w + f) / 2
+        distance_to_next = 1 - (Px - inds)  # always in (0, 1]
+        weights = jnp.minimum(distance_to_next, width) / width
+
+        return inds, weights
+
+
+class XRayTransform3D(LinearOperator):
+    r"""General-purpose, 3D, parallel ray X-ray projector.
+
+    For each view, the projection geometry is specified by an array
+    with shape (2, 4) that specifies a :math:`2 \times 3` projection
+    matrix and a :math:`2 \times 1` offset vector. Denoting the matrix
+    by :math:`\mathbf{M}` and the offset by :math:`\mathbf{t}`, a voxel at array
+    index `(i, j, k)` has its center projected to the detector coordinates
+
+    .. math::
+        \mathbf{M} \begin{bmatrix}
+        i + \frac{1}{2} \\ j + \frac{1}{2} \\ k + \frac{1}{2}
+        \end{bmatrix} + \mathbf{t} \,.
+
+    The detector pixel at index `(i, j)` covers detector coordinates
+    :math:`[i+1) \times [j+1)`.
+
+    :meth:`XRayTransform3D.matrices_from_euler_angles` can help to
+    make these geometry arrays.
 
-    return HTy
 
 
-@partial(jax.jit, static_argnames=["nx"])
-@partial(jax.vmap, in_axes=(None, None, None, 0, None))
-def _calc_weights(x0, dx, nx, angle, y0):
-    """
 
-    Args:
-        x0: Location of the corner of the pixel im[0,0].
-        dx: Pixel side length in x- and y-direction. Units are such
-            that the detector bins have length 1.0.
-        nx: Input image shape.
-        angle: (num_angles,) array of angles in radians. Pixels are
-            projected onto units vectors pointing in these directions.
-            (This argument is `vmap`ed.)
-        y0: Location of the edge of the first detector bin.
     """
-    u = [jnp.cos(angle), jnp.sin(angle)]
-    Px0 = x0[0] * u[0] + x0[1] * u[1] - y0
-    Pdx = [dx[0] * u[0], dx[1] * u[1]]
-    Pxmin = jnp.min(jnp.array([Px0, Px0 + Pdx[0], Px0 + Pdx[1], Px0 + Pdx[0] + Pdx[1]]))
-
-    Px = (
-        Pxmin
-        + Pdx[0] * jnp.arange(nx[0]).reshape(-1, 1)
-        + Pdx[1] * jnp.arange(nx[1]).reshape(1, -1)
-    )
-
-    # detector bin inds
-    inds = jnp.floor(Px).astype(int)
-
-    # weights
-    Pdx = jnp.array(u) * jnp.array(dx)
-    diag1 = jnp.abs(Pdx[0] + Pdx[1])
-    diag2 = jnp.abs(Pdx[0] - Pdx[1])
-    w = jnp.max(jnp.array([diag1, diag2]))
-    f = jnp.min(jnp.array([diag1, diag2]))
-
-    width = (w + f) / 2
-    distance_to_next = 1 - (Px - inds)  # always in (0, 1]
-    weights = jnp.minimum(distance_to_next, width) / width
-
-    return inds, weights
+
+    def __init__(
+        self,
+        input_shape: Shape,
+        matrices: ArrayLike,
+        det_shape: Shape,
+    ):
+        r"""
+        Args:
+            input_shape: Shape of input image.
+            matrices: (num_views, 2, 4) array of homogeneous projection matrices.
+            det_shape: Shape of detector.
+        """
+
+        self.input_shape: Shape = input_shape
+        self.matrices = matrices
+        self.det_shape = det_shape
+        self.output_shape = (len(matrices), *det_shape)
+        super().__init__(
+            input_shape=input_shape,
+            output_shape=self.output_shape,
+            eval_fn=self.project,
+            adj_fn=self.back_project,
+        )
+
+    def project(self, im: ArrayLike) -> snp.Array:
+        """Compute X-ray projection."""
+        return XRayTransform3D._project(im, self.matrices, self.det_shape)
+
+    def back_project(self, proj: ArrayLike) -> snp.Array:
+        """Compute X-ray back projection"""
+        return XRayTransform3D._back_project(proj, self.matrices, self.input_shape)
+
+    @staticmethod
+    def _project(im: ArrayLike, matrices: ArrayLike, det_shape: Shape) -> snp.Array:
+        r"""
+        Args:
+            im: Input image.
+            matrix: (num_views, 2, 4) array of homogeneous projection matrices.
+            det_shape: Shape of detector.
+        """
+        MAX_SLICE_LEN = 10
+        slice_offsets = list(range(0, im.shape[0], MAX_SLICE_LEN))
+
+        num_views = len(matrices)
+        proj = jnp.zeros((num_views,) + det_shape, dtype=im.dtype)
+        for view_ind, matrix in enumerate(matrices):
+            for slice_offset in slice_offsets:
+                proj = proj.at[view_ind].set(
+                    XRayTransform3D._project_single(
+                        im[slice_offset : slice_offset + MAX_SLICE_LEN],
+                        matrix,
+                        proj[view_ind],
+                        slice_offset=slice_offset,
+                    )
+                )
+        return proj
+
+    @staticmethod
+    @partial(jax.jit, donate_argnames="proj")
+    def _project_single(
+        im: ArrayLike, matrix: ArrayLike, proj: ArrayLike, slice_offset: int = 0
+    ) -> snp.Array:
+        r"""
+        Args:
+            im: Input image.
+            matrix: (2, 4) homogeneous projection matrix.
+            det_shape: Shape of detector.
+        """
+
+        ul_ind, ul_weight, ur_weight, ll_weight, lr_weight = XRayTransform3D._calc_weights(
+            im.shape, matrix, proj.shape, slice_offset
+        )
+        proj = proj.at[ul_ind[0], ul_ind[1]].add(ul_weight * im, mode="drop")
+        proj = proj.at[ul_ind[0] + 1, ul_ind[1]].add(ur_weight * im, mode="drop")
+        proj = proj.at[ul_ind[0], ul_ind[1] + 1].add(ll_weight * im, mode="drop")
+        proj = proj.at[ul_ind[0] + 1, ul_ind[1] + 1].add(lr_weight * im, mode="drop")
+        return proj
+
+    @staticmethod
+    def _back_project(proj: ArrayLike, matrices: ArrayLike, input_shape: Shape) -> ArrayLike:
+        r"""
+        Args:
+            proj: Input (set of) projection(s).
+            matrix: (num_views, 2, 4) array of homogeneous projection matrices.
+            input_shape: Shape of desired back projection.
+        """
+        MAX_SLICE_LEN = 10
+        slice_offsets = list(range(0, input_shape[0], MAX_SLICE_LEN))
+
+        HTy = jnp.zeros(input_shape, dtype=proj.dtype)
+        for view_ind, matrix in enumerate(matrices):
+            for slice_offset in slice_offsets:
+                HTy = HTy.at[slice_offset : slice_offset + MAX_SLICE_LEN].set(
+                    XRayTransform3D._back_project_single(
+                        proj[view_ind],
+                        matrix,
+                        HTy[slice_offset : slice_offset + MAX_SLICE_LEN],
+                        slice_offset=slice_offset,
+                    )
+                )
+                HTy.block_until_ready()  # prevent OOM
+
+        return HTy
+
+    @staticmethod
+    @partial(jax.jit, donate_argnames="HTy")
+    def _back_project_single(
+        y: ArrayLike, matrix: ArrayLike, HTy: ArrayLike, slice_offset: int = 0
+    ) -> snp.Array:
+        ul_ind, ul_weight, ur_weight, ll_weight, lr_weight = XRayTransform3D._calc_weights(
+            HTy.shape, matrix, y.shape, slice_offset
+        )
+        HTy = HTy + y[ul_ind[0], ul_ind[1]] * ul_weight
+        HTy = HTy + y[ul_ind[0] + 1, ul_ind[1]] * ur_weight
+        HTy = HTy + y[ul_ind[0], ul_ind[1] + 1] * ll_weight
+        HTy = HTy + y[ul_ind[0] + 1, ul_ind[1] + 1] * lr_weight
+        return HTy
+
+    @staticmethod
+    def _calc_weights(
+        input_shape: Shape, matrix: snp.Array, output_shape: Shape, slice_offset: int = 0
+    ) -> snp.Array:
+        # pixel (0, 0, 0) has its center at (0.5, 0.5, 0.5)
+        x = jnp.mgrid[: input_shape[0], : input_shape[1], : input_shape[2]] + 0.5  # (3, ...)
+        x = x.at[0].add(slice_offset)
+
+        Px = jnp.stack(
+            (
+                matrix[0, 0] * x[0] + matrix[0, 1] * x[1] + matrix[0, 2] * x[2] + matrix[0, 3],
+                matrix[1, 0] * x[0] + matrix[1, 1] * x[1] + matrix[1, 2] * x[2] + matrix[1, 3],
+            )
+        )  # (2, ...)
+
+        # calculate weight on 4 intersecting pixels
+        w = 0.5  # assumed <= 1.0
+        left_edge = Px - w / 2
+        to_next = jnp.minimum(jnp.ceil(left_edge) - left_edge, w)
+        ul_ind = jnp.floor(left_edge).astype("int32")
+        ul_ind = jnp.where(ul_ind < 0, max(output_shape), ul_ind)  # otherwise negative values wrap
+
+        ul_weight = to_next[0] * to_next[1] * (1 / w**2)
+        ur_weight = (w - to_next[0]) * to_next[1] * (1 / w**2)
+        ll_weight = to_next[0] * (w - to_next[1]) * (1 / w**2)
+        lr_weight = (w - to_next[0]) * (w - to_next[1]) * (1 / w**2)
+
+        return ul_ind, ul_weight, ur_weight, ll_weight, lr_weight
+
+    @staticmethod
+    def matrices_from_euler_angles(
+        input_shape: Shape,
+        output_shape: Shape,
+        seq: str,
+        angles: ArrayLike,
+        degrees: bool = False,
+        voxel_spacing: ArrayLike = None,
+        det_spacing: ArrayLike = None,
+    ) -> snp.Array:
+        """
+        Create a set of projection matrices from Euler angles. The
+        input voxels will undergo the specified rotation and then be
+        projected onto the global xy-plane.
+
+        Args:
+            input_shape: Shape of input image.
+            output_shape: Shape of output (detector).
+            str: Sequence of axes for rotation. Up to 3 characters belonging to the set {'X', 'Y', 'Z'}
+                for intrinsic rotations, or {'x', 'y', 'z'} for extrinsic rotations. Extrinsic and
+                intrinsic rotations cannot be mixed in one function call.
+            angles: (num_views, N), N = 1, 2, or 3 Euler angles.
+            degrees: If ``True``, angles are in degrees, otherwise radians. Default: ``True``, radians.
+            voxel_spacing: (3,) array giving the spacing of image
+                voxels.  Default: `[1.0, 1.0, 1.0]`. Experimental.
+            det_spacing: (2,) array giving the spacing of detector
+                pixels.  Default: `[1.0, 1.0]`. Experimental.
+
+
+        Returns:
+            (num_views, 2, 4) array of homogeneous projection matrices.
+        """
+
+        if voxel_spacing is None:
+            voxel_spacing = np.ones(3)
+
+        if det_spacing is None:
+            det_spacing = np.ones(2)
+
+        # make projection matrix: form a rotation matrix and chop off the last row
+        matrices = Rotation.from_euler(seq, angles, degrees=degrees).as_matrix()
+        matrices = matrices[:, :2, :]  # (num_views, 2, 3)
+
+        # handle scaling
+        M_voxel = np.diag(voxel_spacing)  # (3, 3)
+        M_det = np.diag(1 / np.array(det_spacing))  # (2, 2)
+
+        # idea: M_det * M * M_voxel, but with a leading batch dimension
+        matrices = np.einsum("vmn,nn->vmn", matrices, M_voxel)
+        matrices = np.einsum("mm,vmn->vmn", M_det, matrices)
+
+        # add translation to line up the centers
+        x0 = np.array(input_shape) / 2
+        t = -np.einsum("vmn,n->vm", matrices, x0) + np.array(output_shape) / 2
+        matrices = snp.concatenate((matrices, t[..., np.newaxis]), axis=2)
+
+        return matrices
diff --git a/scico/linop/xray/astra.py b/scico/linop/xray/astra.py
index b7eec9191..3f50df6f9 100644
--- a/scico/linop/xray/astra.py
+++ b/scico/linop/xray/astra.py
@@ -20,8 +20,10 @@
 from typing import List, Optional, Sequence, Tuple, Union
 
 import numpy as np
+import numpy.typing
 
 import jax
+from jax.typing import ArrayLike
 
 try:
     import astra
@@ -41,9 +43,11 @@
     # Monkey patching required because latest astra release uses old module path for Iterable
     collections.Iterable = collections.abc.Iterable  # type: ignore
 
-from scico.typing import Shape
+from scico.linop import LinearOperator
+from scico.typing import Shape, TypeAlias
 
-from .._linop import LinearOperator
+VolumeGeometry: TypeAlias = dict
+ProjectionGeometry: TypeAlias = dict
 
 
 def set_astra_gpu_index(idx: Union[int, Sequence[int]]):
@@ -55,6 +59,128 @@ def set_astra_gpu_index(idx: Union[int, Sequence[int]]):
     astra.set_gpu_index(idx)
 
 
+def _project_coords(
+    x_volume: np.ndarray, vol_geom: VolumeGeometry, proj_geom: ProjectionGeometry
+) -> np.ndarray:
+    """
+    Transform volume (logical) coordinates into world coordinates based
+    on ASTRA geometry objects.
+
+    Args:
+        x_volume: (..., 3) vector(s) of volume (AKA logical) coordinates
+        vol_geom: ASTRA volume geometry object.
+        proj_geom: ASTRA projection geometry object.
+    """
+    det_shape = (proj_geom["DetectorRowCount"], proj_geom["DetectorColCount"])
+    x_world = volume_coords_to_world_coords(x_volume, vol_geom=vol_geom)
+    x_dets = []
+    for vec in proj_geom["Vectors"]:
+        ray, d, u, v = vec[0:3], vec[3:6], vec[6:9], vec[9:12]
+        x_det = project_world_coordinates(x_world, ray, d, u, v, det_shape)
+        x_dets.append(x_det)
+
+    return np.stack(x_dets)
+
+
+def project_world_coordinates(
+    x: np.ndarray,
+    ray: np.typing.ArrayLike,
+    d: np.typing.ArrayLike,
+    u: np.typing.ArrayLike,
+    v: np.typing.ArrayLike,
+    det_shape: Sequence[int],
+) -> np.ndarray:
+    """Project world coordinates along ray into the specified basis.
+
+    Project world coordinates along `ray` into the basis described by `u`
+    and `v` with center `d`.
+
+    Args:
+        x: (..., 3) vector(s) of world coordinates.
+        ray: (3,) ray direction
+        d: (3,) center of the detector
+        u: (3,) vector from detector pixel (0,0) to (0,1), columns, x
+        v: (3,) vector from detector pixel (0,0) to (1,0), rows, y
+
+    Returns:
+        (..., 2) vector(s) in the detector coordinates
+
+    """
+    Phi = np.stack((ray, u, v), axis=1)
+    x = x - d  # express with respect to detector center
+    alpha = np.linalg.pinv(Phi) @ x[..., :, np.newaxis]  # (3,3) times <stack of> (3,1)
+    alpha = alpha[..., 0]  # squash from (..., 3, 1) to (..., 3)
+    Palpha = alpha[..., 1:]  # throw away ray coordinate
+    det_center_idx = (
+        np.array(det_shape)[::-1] / 2 - 0.5
+    )  # center of length-2 is index 0.5, length-3 -> index 1
+    ind_xy = Palpha + det_center_idx
+    ind_ij = ind_xy[..., ::-1]
+    return ind_ij
+
+
+def volume_coords_to_world_coords(idx: np.ndarray, vol_geom: VolumeGeometry) -> np.ndarray:
+    """Convert a volume coordinate into a world coordinate.
+
+    Convert a volume coordinate into a world coordinate using ASTRA
+    conventions.
+
+    Args:
+        idx: (..., 2) or (..., 3) vector(s) of index coordinates.
+        vol_geom: ASTRA volume geometry object.
+
+    Returns:
+        (..., 2) or (..., 3) vector(s) of world coordinates.
+
+    """
+    if "GridSliceCount" not in vol_geom:
+        return _volume_index_to_astra_world_2d(idx, vol_geom)
+
+    return _volume_index_to_astra_world_3d(idx, vol_geom)
+
+
+def _volume_index_to_astra_world_2d(idx: np.ndarray, vol_geom: VolumeGeometry) -> np.ndarray:
+    """Convert a 2D volume coordinate into a 2D world coordinate."""
+    coord = idx[..., [2, 1]]  # x:col, y:row,
+    nx = np.array(  # (x, y) order
+        (
+            vol_geom["GridColCount"],
+            vol_geom["GridRowCount"],
+        )
+    )
+    opt = vol_geom["option"]
+    dx = np.array(
+        (
+            (opt["WindowMaxX"] - opt["WindowMinX"]) / nx[0],
+            (opt["WindowMaxY"] - opt["WindowMinY"]) / nx[1],
+        )
+    )
+    center_coord = nx / 2 - 0.5  # center of length-2 is index 0.5, center of length-3 is index 1
+    return (coord - center_coord) * dx
+
+
+def _volume_index_to_astra_world_3d(idx: np.ndarray, vol_geom: VolumeGeometry) -> np.ndarray:
+    """Convert a 3D volume coordinate into a 3D world coordinate."""
+    coord = idx[..., [2, 1, 0]]  # x:col, y:row, z:slice
+    nx = np.array(  # (x, y, z) order
+        (
+            vol_geom["GridColCount"],
+            vol_geom["GridRowCount"],
+            vol_geom["GridSliceCount"],
+        )
+    )
+    opt = vol_geom["option"]
+    dx = np.array(
+        (
+            (opt["WindowMaxX"] - opt["WindowMinX"]) / nx[0],
+            (opt["WindowMaxY"] - opt["WindowMinY"]) / nx[1],
+            (opt["WindowMaxZ"] - opt["WindowMinZ"]) / nx[2],
+        )
+    )
+    center_coord = nx / 2 - 0.5  # center of length-2 is index 0.5, center of length-3 is index 1
+    return (coord - center_coord) * dx
+
+
 class XRayTransform2D(LinearOperator):
     r"""2D parallel beam X-ray transform based on the ASTRA toolbox.
 
@@ -221,6 +347,113 @@ def f(sino):
         return jax.pure_callback(f, jax.ShapeDtypeStruct(self.input_shape, self.input_dtype), sino)
 
 
+def convert_from_scico_geometry(
+    in_shape: Shape, matrices: ArrayLike, det_shape: Shape
+) -> np.ndarray:
+    """Convert SCICO projection matrices into ASTRA "parallel3d_vec" vectors.
+
+    For 3D arrays,
+    in ASTRA, the dimensions go (slices, rows, columns) and (z, y, x);
+    in SCICO, the dimensions go (x, y, z).
+
+    In ASTRA, the x-grid (recon) is centered on the origin and the y-grid (projection) can move.
+    In SCICO, the x-grid origin is the center of x[0, 0, 0], the y-grid origin is the center
+    of y[0, 0].
+
+    See section "parallel3d_vec" in the
+    `astra documentation <https://astra-toolbox.com/docs/geom3d.html#projection-geometries>`__.
+
+    Args:
+        in_shape: Shape of input image.
+        matrices: (num_angles, 2, 4) array of homogeneous projection matrices.
+        det_shape: Shape of detector.
+
+    Returns:
+        (num_angles, 12) vector array in the ASTRA "parallel3d_vec" convention.
+    """
+    # ray is perpendicular to projection axes
+    ray = np.cross(matrices[:, 0, :3], matrices[:, 1, :3])
+    # detector center comes from lifting the center index to 3D
+    y_center = (np.array(det_shape) - 1) / 2
+    x_center = (
+        np.einsum("...mn,n->...m", matrices[..., :3], (np.array(in_shape) - 1) / 2)
+        + matrices[..., 3]
+    )
+    d = np.einsum("...mn,...m->...n", matrices[..., :3], y_center - x_center)  # (V, 2, 3) x (V, 2)
+    u = matrices[:, 1, :3]
+    v = matrices[:, 0, :3]
+
+    # handle different axis conventions
+    ray = ray[:, [2, 1, 0]]
+    d = d[:, [2, 1, 0]]
+    u = u[:, [2, 1, 0]]
+    v = v[:, [2, 1, 0]]
+
+    vectors = np.concatenate((ray, d, u, v), axis=1)  # (v, 12)
+    return vectors
+
+
+def _astra_to_scico_geometry(vol_geom: VolumeGeometry, proj_geom: ProjectionGeometry) -> np.ndarray:
+    """Convert ASTRA geometry objects into a SCICO projection matrix.
+
+    Convert ASTRA volume and projection geometry into a SCICO X-ray
+    projection matrix, assuming "parallel3d_vec" format.
+
+    The approach is to locate 3 points in the volume domain,
+    deduce the corresponding projection locations, and, then, solve a
+    linear system to determine the affine relationship between them.
+
+    Args:
+        vol_geom: ASTRA volume geometry object.
+        proj_geom: ASTRA projection geometry object.
+
+    Returns:
+        (num_angles, 2, 4) array of homogeneous projection matrices.
+
+    """
+    x_volume = np.concatenate((np.zeros((1, 3)), np.eye(3)), axis=0)  # (4, 3)
+    x_dets = _project_coords(x_volume, vol_geom, proj_geom)  # (1, 4, 2)
+
+    x_volume_aug = np.concatenate((x_volume, np.ones((4, 1))), axis=1)  # (4, 4)
+    matrices = []
+    for x_det in x_dets:
+        M = np.linalg.solve(x_volume_aug, x_det).T
+        np.testing.assert_allclose(M @ x_volume_aug[0], x_det[0])
+        matrices.append(M)
+
+    return np.stack(matrices)
+
+
+def convert_to_scico_geometry(
+    input_shape: Shape,
+    det_count: Tuple[int, int],
+    det_spacing: Optional[Tuple[float, float]] = None,
+    angles: Optional[np.ndarray] = None,
+    vectors: Optional[np.ndarray] = None,
+) -> np.ndarray:
+    """Convert ASTRA geometry specificiation to a SCICO projection matrix.
+
+    Convert ASTRA volume and projection geometry into a SCICO X-ray
+    projection matrix, assuming "parallel3d_vec" format.
+
+    The approach is to locate 3 points in the volume domain,
+    deduce the corresponding projection locations, and, then, solve a
+    linear system to determine the affine relationship between them.
+
+    Args:
+        vol_geom: ASTRA volume geometry object.
+        proj_geom: ASTRA projection geometry object.
+
+    Returns:
+        (num_angles, 2, 4) array of homogeneous projection matrices.
+
+    """
+    vol_geom, proj_geom = XRayTransform3D.create_astra_geometry(
+        input_shape, det_count, det_spacing=det_spacing, angles=angles, vectors=vectors
+    )
+    return _astra_to_scico_geometry(vol_geom, proj_geom)
+
+
 class XRayTransform3D(LinearOperator):  # pragma: no cover
     r"""3D parallel beam X-ray transform based on the ASTRA toolbox.
 
@@ -354,32 +587,25 @@ def __init__(
             raise ValueError("Expected det_count to be a tuple with 2 elements.")
         if angles is not None:
             Nview = angles.size
-            self.angles: np.ndarray = np.array(angles)
+            self.angles: Optional[np.ndarray] = np.array(angles)
+            self.vectors: Optional[np.ndarray] = None
         else:
             assert vectors is not None
             Nview = vectors.shape[0]
-            self.vectors: np.ndarray = np.array(vectors)
+            self.vectors = np.array(vectors)
+            self.angles = None
         output_shape: Shape = (det_count[0], Nview, det_count[1])
 
         self.det_count = det_count
         assert isinstance(det_count, (list, tuple))
-        if angles is not None:
-            assert det_spacing is not None
-            self.proj_geom = astra.create_proj_geom(
-                "parallel3d",
-                det_spacing[0],
-                det_spacing[1],
-                det_count[0],
-                det_count[1],
-                self.angles,
-            )
-        else:
-            self.proj_geom = astra.create_proj_geom(
-                "parallel3d_vec", det_count[0], det_count[1], self.vectors
-            )
-
         self.input_shape: tuple = input_shape
-        self.vol_geom = astra.create_vol_geom(input_shape[1], input_shape[2], input_shape[0])
+        self.vol_geom, self.proj_geom = self.create_astra_geometry(
+            input_shape,
+            det_count,
+            det_spacing=det_spacing,
+            angles=self.angles,
+            vectors=self.vectors,
+        )
 
         # Wrap our non-jax function to indicate we will supply fwd/rev mode functions
         self._eval = jax.custom_vjp(self._proj)
@@ -396,6 +622,53 @@ def __init__(
             jit=False,
         )
 
+    @staticmethod
+    def create_astra_geometry(
+        input_shape: Shape,
+        det_count: Tuple[int, int],
+        det_spacing: Optional[Tuple[float, float]] = None,
+        angles: Optional[np.ndarray] = None,
+        vectors: Optional[np.ndarray] = None,
+    ) -> Tuple[VolumeGeometry, ProjectionGeometry]:
+        """Create ASTRA 3D geometry objects.
+
+        Keyword arguments `det_spacing` and `angles` should be specified
+        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.
+            det_count: Number of detector elements. See the
+               `astra documentation <https://www.astra-toolbox.com/docs/geom3d.html#projection-geometries>`__
+               for more information.
+            det_spacing: Spacing between detector elements. See the
+               `astra documentation <https://www.astra-toolbox.com/docs/geom3d.html#projection-geometries>`__
+               for more information.
+            angles: Array of projection angles in radians.
+            vectors: Array of geometry specification vectors.
+
+        Returns:
+            A tuple `(vol_geom, proj_geom)` of ASTRA volume geometry and
+            projection geometry objects.
+        """
+        vol_geom = astra.create_vol_geom(input_shape[1], input_shape[2], input_shape[0])
+        if angles is not None:
+            assert det_spacing is not None
+            proj_geom = astra.create_proj_geom(
+                "parallel3d",
+                det_spacing[0],
+                det_spacing[1],
+                det_count[0],
+                det_count[1],
+                angles,
+            )
+        else:
+            proj_geom = astra.create_proj_geom(
+                "parallel3d_vec", det_count[0], det_count[1], vectors
+            )
+        return vol_geom, proj_geom
+
     def _proj(self, x: jax.Array) -> jax.Array:
         # apply the forward projector and generate a sinogram
 
diff --git a/scico/test/linop/xray/test_astra.py b/scico/test/linop/xray/test_astra.py
index 55c6206ac..c94314fdd 100644
--- a/scico/test/linop/xray/test_astra.py
+++ b/scico/test/linop/xray/test_astra.py
@@ -104,7 +104,9 @@ def test_grad(testobj):
     A = testobj.A
     x = testobj.x
     g = lambda x: jax.numpy.linalg.norm(A(x)) ** 2
-    np.testing.assert_allclose(scico.grad(g)(x), 2 * A.adj(A(x)), rtol=get_tol())
+    np.testing.assert_allclose(
+        scico.grad(g)(x), 2 * A.adj(A(x)), atol=get_tol() * x.max(), rtol=np.inf
+    )
 
 
 def test_adjoint_grad(testobj):
@@ -174,6 +176,110 @@ def test_angle_to_vector():
     assert vectors.shape == (angles.size, 12)
 
 
+## conversion functions
+@pytest.fixture(scope="module")
+def test_geometry():
+    """
+    In this geometry, if vol[i, j, k]==1, we expect proj[j-2, k-1]==1.
+
+    Because:
+    - We project along z, i.e. `ray=(0,0,1)`, i.e., we remove axis=0.
+    - We set `v=(0, 1, 0)`, so detector rows go with y axis, axis=1.
+    - We set `u=(1, 0, 0)`, so detector columns go with x axis, axis=2.
+    - We shift the detector by (x=1, y=2, z=3) <-> i-3, j-2, k-1
+    """
+    in_shape = (30, 31, 32)
+    # in ASTRA terminology:
+    n_rows = in_shape[1]  # y
+    n_cols = in_shape[2]  # x
+    n_slices = in_shape[0]  # z
+    vol_geom = scico.linop.xray.astra.astra.create_vol_geom(n_rows, n_cols, n_slices)
+
+    assert vol_geom["option"]["WindowMinX"] == -n_cols / 2
+    assert vol_geom["option"]["WindowMinY"] == -n_rows / 2
+    assert vol_geom["option"]["WindowMinZ"] == -n_slices / 2
+
+    # project along z, axis=0
+    det_row_count = n_rows
+    det_col_count = n_cols
+    ray = (0, 0, 1)
+    d = (1, 2, 3)  # axis=2 offset by 1, axis=1 offset by 2, axis=0 offset by 3
+    u = (1, 0, 0)  # increments columns, goes with X
+    v = (0, 1, 0)  # increments rows, goes with Y
+    vectors = np.array(ray + d + u + v)[np.newaxis, :]
+    proj_geom = scico.linop.xray.astra.astra.create_proj_geom(
+        "parallel3d_vec", det_row_count, det_col_count, vectors
+    )
+
+    return vol_geom, proj_geom
+
+
+@pytest.mark.skipif(jax.devices()[0].platform != "gpu", reason="GPU required for test")
+def test_projection_convention(test_geometry):
+    """
+    If vol[i, j, k]==1, test that astra puts proj[j-2, k-1]==1.
+
+    See `test_geometry` for the setup.
+    """
+    vol_geom, proj_geom = test_geometry
+    in_shape = scico.linop.xray.astra.astra.functions.geom_size(vol_geom)
+    vol = np.zeros(in_shape)
+
+    i, j, k = [np.random.randint(0, s) for s in in_shape]
+    vol[i, j, k] = 1.0
+
+    proj_id, proj = scico.linop.xray.astra.astra.create_sino3d_gpu(vol, proj_geom, vol_geom)
+    scico.linop.xray.astra.astra.data3d.delete(proj_id)
+    proj = proj[:, 0, :]  # get first view
+    assert len(np.unique(proj) == 2)
+
+    idx_proj_i, idx_proj_j = np.nonzero(proj)
+    np.testing.assert_array_equal(idx_proj_i, j - 2)
+    np.testing.assert_array_equal(idx_proj_j, k - 1)
+
+
+def test_project_coords(test_geometry):
+    """
+    If vol[i, j, k]==1, test that we predict proj[j-2, k-1]==1.
+
+    See `test_geometry` for the setup and `test_projection_convention`
+    for proof ASTRA works this way.
+    """
+    vol_geom, proj_geom = test_geometry
+    in_shape = scico.linop.xray.astra.astra.functions.geom_size(vol_geom)
+    x_vol = np.array([np.random.randint(0, s) for s in in_shape])
+    x_proj_gt = np.array(
+        [[x_vol[1] - 2, x_vol[2] - 1]]
+    )  # projection along slices removes first index
+    x_proj = scico.linop.xray.astra._project_coords(x_vol, vol_geom, proj_geom)
+    np.testing.assert_array_equal(x_proj_gt, x_proj)
+
+
+def test_convert_to_scico_geometry(test_geometry):
+    """
+    Basic regression test, `test_project_coords` tests the logic.
+    """
+    vol_geom, proj_geom = test_geometry
+    matrices_truth = scico.linop.xray.astra._astra_to_scico_geometry(vol_geom, proj_geom)
+    truth = np.array([[[0.0, 1.0, 0.0, -2.0], [0.0, 0.0, 1.0, -1.0]]])
+    np.testing.assert_allclose(matrices_truth, truth)
+
+
+def test_convert_from_scico_geometry(test_geometry):
+    """
+    Basic regression test, `test_project_coords` tests the logic.
+    """
+    in_shape = (30, 31, 32)
+    matrices = np.array([[[0.0, 1.0, 0.0, -2.0], [0.0, 0.0, 1.0, -1.0]]])
+    det_shape = (31, 32)
+    vectors = scico.linop.xray.astra.convert_from_scico_geometry(in_shape, matrices, det_shape)
+
+    _, proj_geom_truth = test_geometry
+    # skip testing element 5, as it is detector center along the ray and doesn't matter
+    np.testing.assert_allclose(vectors[0, :5], proj_geom_truth["Vectors"][0, :5])
+    np.testing.assert_allclose(vectors[0, 6:], proj_geom_truth["Vectors"][0, 6:])
+
+
 def test_ensure_writeable():
     assert isinstance(_ensure_writeable(np.ones((2, 1))), np.ndarray)
     assert isinstance(_ensure_writeable(snp.ones((2, 1))), np.ndarray)
diff --git a/scico/test/linop/xray/test_xray.py b/scico/test/linop/xray/test_xray.py
index 6d9c2ba39..cd7c0dcdd 100644
--- a/scico/test/linop/xray/test_xray.py
+++ b/scico/test/linop/xray/test_xray.py
@@ -1,9 +1,11 @@
+import numpy as np
+
 import jax.numpy as jnp
 
 import pytest
 
 import scico
-from scico.linop import Parallel2dProjector, XRayTransform
+from scico.linop.xray import XRayTransform2D, XRayTransform3D
 
 
 @pytest.mark.filterwarnings("error")
@@ -11,22 +13,18 @@ def test_init():
     input_shape = (3, 3)
 
     # no warning with default settings, even at 45 degrees
-    H = XRayTransform(Parallel2dProjector(input_shape, jnp.array([jnp.pi / 4])))
+    H = XRayTransform2D(input_shape, jnp.array([jnp.pi / 4]))
 
     # no warning if we project orthogonally with oversized pixels
-    H = XRayTransform(Parallel2dProjector(input_shape, jnp.array([0]), dx=jnp.array([1, 1])))
+    H = XRayTransform2D(input_shape, jnp.array([0]), dx=jnp.array([1, 1]))
 
     # warning if the projection angle changes
     with pytest.warns(UserWarning):
-        H = XRayTransform(
-            Parallel2dProjector(input_shape, jnp.array([0.1]), dx=jnp.array([1.1, 1.1]))
-        )
+        H = XRayTransform2D(input_shape, jnp.array([0.1]), dx=jnp.array([1.1, 1.1]))
 
     # warning if the pixels get any larger
     with pytest.warns(UserWarning):
-        H = XRayTransform(
-            Parallel2dProjector(input_shape, jnp.array([0]), dx=jnp.array([1.1, 1.1]))
-        )
+        H = XRayTransform2D(input_shape, jnp.array([0]), dx=jnp.array([1.1, 1.1]))
 
 
 def test_apply():
@@ -37,13 +35,13 @@ def test_apply():
     angles = jnp.linspace(0, jnp.pi, num=num_angles, endpoint=False)
 
     # general projection
-    H = XRayTransform(Parallel2dProjector(x.shape, angles))
+    H = XRayTransform2D(x.shape, angles)
     y = H @ x
     assert y.shape[0] == (num_angles)
 
     # fixed det_count
     det_count = 14
-    H = XRayTransform(Parallel2dProjector(x.shape, angles, det_count=det_count))
+    H = XRayTransform2D(x.shape, angles, det_count=det_count)
     y = H @ x
     assert y.shape[1] == det_count
 
@@ -56,7 +54,7 @@ def test_apply_adjoint():
     angles = jnp.linspace(0, jnp.pi, num=num_angles, endpoint=False)
 
     # general projection
-    H = XRayTransform(Parallel2dProjector(x.shape, angles))
+    H = XRayTransform2D(x.shape, angles)
     y = H @ x
     assert y.shape[0] == (num_angles)
 
@@ -68,6 +66,49 @@ def test_apply_adjoint():
 
     # fixed det_length
     det_count = 14
-    H = XRayTransform(Parallel2dProjector(x.shape, angles, det_count=det_count))
+    H = XRayTransform2D(x.shape, angles, det_count=det_count)
     y = H @ x
     assert y.shape[1] == det_count
+
+
+def test_3d_scaling():
+    x = jnp.zeros((4, 4, 1))
+    x = x.at[1:3, 1:3, 0].set(1.0)
+
+    input_shape = x.shape
+    output_shape = x.shape[:2]
+
+    # default spacing
+    M = XRayTransform3D.matrices_from_euler_angles(input_shape, output_shape, "X", [0.0])
+    H = XRayTransform3D(input_shape, matrices=M, det_shape=output_shape)
+    # fmt: off
+    truth = jnp.array(
+        [[[0.0, 0.0, 0.0, 0.0],
+          [0.0, 1.0, 1.0, 0.0],
+          [0.0, 1.0, 1.0, 0.0],
+          [0.0, 0.0, 0.0, 0.0]]]
+    )  # fmt: on
+    np.testing.assert_allclose(H @ x, truth)
+
+    # bigger voxels in the x (first index) direction
+    M = XRayTransform3D.matrices_from_euler_angles(
+        input_shape, output_shape, "X", [0.0], voxel_spacing=[2.0, 1.0, 1.0]
+    )
+    H = XRayTransform3D(input_shape, matrices=M, det_shape=output_shape)
+    # fmt: off
+    truth = jnp.array(
+        [[[0. , 0.5, 0.5, 0. ],
+          [0. , 0.5, 0.5, 0. ],
+          [0. , 0.5, 0.5, 0. ],
+          [0. , 0.5, 0.5, 0. ]]]
+    )  # fmt: on
+    np.testing.assert_allclose(H @ x, truth)
+
+    # bigger detector pixels in the x (first index) direction
+    M = XRayTransform3D.matrices_from_euler_angles(
+        input_shape, output_shape, "X", [0.0], det_spacing=[2.0, 1.0]
+    )
+    H = XRayTransform3D(input_shape, matrices=M, det_shape=output_shape)
+    # fmt: off
+    truth = None  # fmt: on  # TODO: Check this case more closely.
+    # np.testing.assert_allclose(H @ x, truth)
diff --git a/scico/test/operator/test_operator.py b/scico/test/operator/test_operator.py
index f76343ff1..b4f5762f0 100644
--- a/scico/test/operator/test_operator.py
+++ b/scico/test/operator/test_operator.py
@@ -264,19 +264,19 @@ def setup_method(self):
     def test_jvp(self):
         Fu, JFuv = self.F.jvp(self.u, self.v)
         np.testing.assert_allclose(Fu, self.F(self.u))
-        np.testing.assert_allclose(JFuv, self.fmx @ self.v, rtol=1e-6)
+        np.testing.assert_allclose(JFuv, self.fmx @ self.v, atol=1e-6, rtol=0.0)
 
     def test_vjp_conj(self):
         Fu, G = self.F.vjp(self.u, conjugate=True)
         JFTw = G(self.w)
         np.testing.assert_allclose(Fu, self.F(self.u))
-        np.testing.assert_allclose(JFTw, self.fmx.T @ self.w, rtol=1e-6)
+        np.testing.assert_allclose(JFTw, self.fmx.T @ self.w, atol=1e-6, rtol=0.0)
 
     def test_vjp_noconj(self):
         Fu, G = self.F.vjp(self.u, conjugate=False)
         JFTw = G(self.w)
         np.testing.assert_allclose(Fu, self.F(self.u))
-        np.testing.assert_allclose(JFTw, self.fmx.T @ self.w, rtol=1e-6)
+        np.testing.assert_allclose(JFTw, self.fmx.T @ self.w, atol=1e-6, rtol=0.0)
 
 
 class TestJacobianProdComplex: