diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index d8041ef4149..749dc98a09f 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -116,7 +116,7 @@ int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
-  int max_mats, int32_t* materials, double* volumes);
+  int max_mats, int32_t* materials, double* volumes, double* bboxes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 7ceb623dae1..852d301367e 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -87,8 +87,12 @@ namespace detail {
 
 class MaterialVolumes {
 public:
+  MaterialVolumes(int32_t* mats, double* vols, double* bboxes, int table_size)
+    : materials_(mats), volumes_(vols), bboxes_(bboxes), table_size_(table_size)
+  {}
+
   MaterialVolumes(int32_t* mats, double* vols, int table_size)
-    : materials_(mats), volumes_(vols), table_size_(table_size)
+    : MaterialVolumes(mats, vols, nullptr, table_size)
   {}
 
   //! Add volume for a given material in a mesh element
@@ -96,8 +100,11 @@ class MaterialVolumes {
   //! \param[in] index_elem Index of the mesh element
   //! \param[in] index_material Index of the material within the model
   //! \param[in] volume Volume to add
-  void add_volume(int index_elem, int index_material, double volume);
-  void add_volume_unsafe(int index_elem, int index_material, double volume);
+  //! \param[in] bbox Bounding box to union into the result (optional)
+  void add_volume(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
+  void add_volume_unsafe(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
 
   // Accessors
   int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; }
@@ -112,11 +119,23 @@ class MaterialVolumes {
     return volumes_[i * table_size_ + j];
   }
 
+  double& bboxes(int i, int j, int k)
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+  const double& bboxes(int i, int j, int k) const
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+
+  bool has_bboxes() const { return bboxes_ != nullptr; }
+
   bool table_full() const { return table_full_; }
 
 private:
   int32_t* materials_;      //!< material index (bins, table_size)
   double* volumes_;         //!< volume in [cm^3] (bins, table_size)
+  double* bboxes_;          //!< bounding boxes (bins, table_size, 6)
   int table_size_;          //!< Size of hash table for each mesh element
   bool table_full_ {false}; //!< Whether the hash table is full
 };
@@ -239,6 +258,19 @@ class Mesh {
   void material_volumes(int nx, int ny, int nz, int max_materials,
     int32_t* materials, double* volumes) const;
 
+  //! Determine volume and bounding boxes of materials within each mesh element
+  //
+  //! \param[in] nx Number of samples in x direction
+  //! \param[in] ny Number of samples in y direction
+  //! \param[in] nz Number of samples in z direction
+  //! \param[in] max_materials Maximum number of materials in a single mesh
+  //!                          element
+  //! \param[inout] materials Array storing material indices
+  //! \param[inout] volumes Array storing volumes
+  //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6)
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes, double* bboxes) const;
+
   //! Determine bounding box of mesh
   //
   //! \return Bounding box of mesh
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index 97277cbda84..57bbe437ffb 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -269,6 +269,7 @@ def step1_neutron_transport(
             # Compute material volume fractions on the mesh
             if mat_vol_kwargs is None:
                 mat_vol_kwargs = {}
+            mat_vol_kwargs.setdefault('bounding_boxes', True)
             self.results['mesh_material_volumes'] = mmv = comm.bcast(
                 self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
 
@@ -539,14 +540,15 @@ def step3_photon_transport(
     def get_decay_photon_source_mesh(
         self,
         time_index: int = -1
-    ) -> list[openmc.MeshSource]:
+    ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
-        This function creates N :class:`MeshSource` objects where N is the
-        maximum number of unique materials that appears in a single mesh
-        element. For each mesh element-material combination, and
-        IndependentSource instance is created with a spatial constraint limited
-        the sampled decay photons to the correct region.
+        For each mesh element-material combination, an
+        :class:`~openmc.IndependentSource` is created with a
+        :class:`~openmc.stats.Box` spatial distribution based on the bounding
+        box of the material within the mesh element. A material constraint is
+        also applied so that sampled source sites are limited to the correct
+        region.
 
         When the photon transport model is different from the neutron model, the
         photon MeshMaterialVolumes is used to determine whether an (element,
@@ -559,19 +561,15 @@ def get_decay_photon_source_mesh(
 
         Returns
         -------
-        list of openmc.MeshSource
-            A list of MeshSource objects, each containing IndependentSource
-            instances for the decay photons in the corresponding mesh element.
+        list of openmc.IndependentSource
+            A list of IndependentSource objects for the decay photons, one for
+            each mesh element-material combination with non-zero source strength.
 
         """
         mat_dict = self.neutron_model._get_all_materials()
 
-        # Some MeshSource objects will have empty positions; create a "null source"
-        # that is used for this case
-        null_source = openmc.IndependentSource(particle='photon', strength=0.0)
-
-        # List to hold sources for each MeshSource (length = N)
-        source_lists = []
+        # List to hold all sources
+        sources = []
 
         # Index in the overall list of activated materials
         index_mat = 0
@@ -594,7 +592,7 @@ def get_decay_photon_source_mesh(
                     if mat_id is not None
                 }
 
-            for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)):
+            for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
                 # Skip void volume
                 if mat_id is None:
                     continue
@@ -604,30 +602,27 @@ def get_decay_photon_source_mesh(
                     index_mat += 1
                     continue
 
-                # Check whether a new MeshSource object is needed
-                if j >= len(source_lists):
-                    source_lists.append([null_source]*n_elements)
-
                 # Get activated material composition
                 original_mat = materials[index_mat]
                 activated_mat = results[time_index].get_material(str(original_mat.id))
 
-                # Create decay photon source source
+                # Create decay photon source
                 energy = activated_mat.get_decay_photon_energy()
                 if energy is not None:
                     strength = energy.integral()
-                    source_lists[j][index_elem] = openmc.IndependentSource(
+                    space = openmc.stats.Box(*bbox)
+                    sources.append(openmc.IndependentSource(
+                        space=space,
                         energy=energy,
                         particle='photon',
                         strength=strength,
                         constraints={'domains': [mat_dict[mat_id]]}
-                    )
+                    ))
 
                 # Increment index of activated material
                 index_mat += 1
 
-        # Return list of mesh sources
-        return [openmc.MeshSource(self.domains, sources) for sources in source_lists]
+        return sources
 
     def load_results(self, path: PathLike):
         """Load results from a previous R2S calculation.
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 3a720f98a31..19e6f74d7ad 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -1,5 +1,5 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
+from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p,
                     create_string_buffer, c_size_t)
 from math import sqrt
 import sys
@@ -47,7 +47,8 @@
 _dll.openmc_mesh_bounding_box.restype = c_int
 _dll.openmc_mesh_bounding_box.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
-    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
+    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double,
+    c_void_p]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
@@ -188,6 +189,7 @@ def material_volumes(
             n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
             output: bool = True,
+            bounding_boxes: bool = False,
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
@@ -213,6 +215,11 @@ def material_volumes(
             Estimated maximum number of materials in any given mesh element.
         output : bool, optional
             Whether or not to show output.
+        bounding_boxes : bool, optional
+            Whether or not to compute an axis-aligned bounding box for each
+            (mesh element, material) combination. When enabled, the bounding
+            box encloses the ray-estimator prisms used for the volume
+            estimation.
 
         Returns
         -------
@@ -243,23 +250,36 @@ def material_volumes(
         table_size = slot_factor*max_materials
         materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
         volumes = np.zeros((n, table_size), dtype=np.float64)
+        bboxes = None
+        if bounding_boxes:
+            bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+            bboxes[..., 0:3] = np.inf
+            bboxes[..., 3:6] = -np.inf
 
         # Run material volume calculation
         while True:
             try:
+                bboxes_ptr = None
+                if bboxes is not None:
+                    bboxes_ptr = bboxes.ctypes.data_as(POINTER(c_double))
                 with quiet_dll(output):
                     _dll.openmc_mesh_material_volumes(
-                        self._index, nx, ny, nz, table_size, materials, volumes)
+                        self._index, nx, ny, nz, table_size, materials,
+                        volumes, bboxes_ptr)
             except AllocationError:
                 # Increase size of result array and try again
                 table_size *= 2
                 materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
                 volumes = np.zeros((n, table_size), dtype=np.float64)
+                if bounding_boxes:
+                    bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+                    bboxes[..., 0:3] = np.inf
+                    bboxes[..., 3:6] = -np.inf
             else:
                 # If no error, break out of loop
                 break
 
-        return MeshMaterialVolumes(materials, volumes)
+        return MeshMaterialVolumes(materials, volumes, bboxes)
 
     def get_plot_bins(
             self,
diff --git a/openmc/mesh.py b/openmc/mesh.py
index ce5218b5e97..296a9c1f53b 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -17,6 +17,7 @@
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from openmc.utility_funcs import change_directory
+from .bounding_box import BoundingBox
 from ._xml import get_elem_list, get_text
 from .mixin import IDManagerMixin
 from .surface import _BOUNDARY_TYPES
@@ -40,6 +41,11 @@ class MeshMaterialVolumes(Mapping):
         Array of shape (elements, max_materials) storing material IDs
     volumes : numpy.ndarray
         Array of shape (elements, max_materials) storing material volumes
+    bboxes : numpy.ndarray, optional
+        Array of shape (elements, max_materials, 6) storing axis-aligned
+        bounding boxes for each (element, material) combination with ordering
+        (xmin, ymin, zmin, xmax, ymax, zmax). Bounding boxes enclose the
+        ray-estimator prisms used to compute volumes.
 
     See Also
     --------
@@ -64,9 +70,30 @@ class MeshMaterialVolumes(Mapping):
     [(2, 31.87963824195591), (1, 6.129949130817542)]
 
     """
-    def __init__(self, materials: np.ndarray, volumes: np.ndarray):
+    def __init__(
+        self,
+        materials: np.ndarray,
+        volumes: np.ndarray,
+        bboxes: np.ndarray | None = None
+    ):
         self._materials = materials
         self._volumes = volumes
+        self._bboxes = bboxes
+
+        if self._bboxes is not None:
+            if self._bboxes.shape[:2] != self._materials.shape:
+                raise ValueError(
+                    'bboxes must have shape (elements, max_materials, 6) '
+                    'matching materials/volumes.'
+                )
+            if self._bboxes.shape[2] != 6:
+                raise ValueError(
+                    'bboxes must have shape (elements, max_materials, 6).'
+                )
+
+    @property
+    def has_bounding_boxes(self) -> bool:
+        return self._bboxes is not None
 
     @property
     def num_elements(self) -> int:
@@ -92,7 +119,11 @@ def __getitem__(self, material_id: int) -> np.ndarray:
             volumes[indices] = self._volumes[indices, i]
         return volumes
 
-    def by_element(self, index_elem: int) -> list[tuple[int | None, float]]:
+    def by_element(
+        self,
+        index_elem: int,
+        include_bboxes: bool = False
+    ) -> list[tuple[int | None, float] | tuple[int | None, float, BoundingBox | None]]:
         """Get a list of volumes for each material within a specific element.
 
         Parameters
@@ -102,15 +133,32 @@ def by_element(self, index_elem: int) -> list[tuple[int | None, float]]:
 
         Returns
         -------
-        list of tuple of (material ID, volume)
+        list of tuple
+            If ``include_bboxes`` is False (default), returns tuples of
+            (material ID, volume). If ``include_bboxes`` is True, returns
+            tuples of (material ID, volume, bounding box).
 
         """
         table_size = self._volumes.shape[1]
-        return [
-            (m if m > -1 else None, self._volumes[index_elem, i])
-            for i in range(table_size)
-            if (m := self._materials[index_elem, i]) != -2
-        ]
+        if include_bboxes and self._bboxes is None:
+            raise ValueError('Bounding boxes were not computed for this object.')
+
+        results = []
+        for i in range(table_size):
+            m = self._materials[index_elem, i]
+            if m == -2:
+                continue
+            mat_id = m if m > -1 else None
+            vol = self._volumes[index_elem, i]
+
+            if include_bboxes:
+                vals = self._bboxes[index_elem, i]
+                bbox = BoundingBox(vals[0:3], vals[3:6])
+                results.append((mat_id, vol, bbox))
+            else:
+                results.append((mat_id, vol))
+
+        return results
 
     def save(self, filename: PathLike):
         """Save material volumes to a .npz file.
@@ -120,8 +168,10 @@ def save(self, filename: PathLike):
         filename : path-like
             Filename where data will be saved
         """
-        np.savez_compressed(
-            filename, materials=self._materials, volumes=self._volumes)
+        kwargs = {'materials': self._materials, 'volumes': self._volumes}
+        if self._bboxes is not None:
+            kwargs['bboxes'] = self._bboxes
+        np.savez_compressed(filename, **kwargs)
 
     @classmethod
     def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes:
@@ -134,7 +184,8 @@ def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes:
 
         """
         filedata = np.load(filename)
-        return cls(filedata['materials'], filedata['volumes'])
+        bboxes = filedata['bboxes'] if 'bboxes' in filedata.files else None
+        return cls(filedata['materials'], filedata['volumes'], bboxes)
 
 
 class MeshBase(IDManagerMixin, ABC):
@@ -366,6 +417,7 @@ def material_volumes(
             model: openmc.Model,
             n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
+            bounding_boxes: bool = False,
             **kwargs
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
@@ -388,6 +440,11 @@ def material_volumes(
             the x, y, and z dimensions.
         max_materials : int, optional
             Estimated maximum number of materials in any given mesh element.
+        bounding_boxes : bool, optional
+            Whether to compute an axis-aligned bounding box for each
+            (mesh element, material) combination. When enabled, the bounding
+            box encloses the ray-estimator prisms used for the volume
+            estimation.
         **kwargs : dict
             Keyword arguments passed to :func:`openmc.lib.init`
 
@@ -419,7 +476,8 @@ def material_volumes(
 
             # Compute material volumes
             volumes = mesh.material_volumes(
-                n_samples, max_materials, output=kwargs['output'])
+                n_samples, max_materials, output=kwargs['output'],
+                bounding_boxes=bounding_boxes)
 
         # Restore original tallies
         model.tallies = original_tallies
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 9a7a8e7517b..b84a45a6a87 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1,6 +1,8 @@
 #include "openmc/mesh.h"
 #include <algorithm> // for copy, equal, min, min_element
 #include <cassert>
+#include <cstdint>        // for uint64_t
+#include <cstring>        // for memcpy
 #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
 #include <cmath>          // for ceil
 #include <cstddef>        // for size_t
@@ -140,6 +142,63 @@ inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired)
 #endif
 }
 
+// Helper function equivalent to std::bit_cast in C++20
+template<typename To, typename From>
+inline To bit_cast_value(const From& value)
+{
+  To out;
+  std::memcpy(&out, &value, sizeof(To));
+  return out;
+}
+
+inline void atomic_update_double(double* ptr, double value, bool is_min)
+{
+#if defined(__GNUC__) || defined(__clang__)
+  using may_alias_uint64_t [[gnu::may_alias]] = uint64_t;
+  auto* bits_ptr = reinterpret_cast<may_alias_uint64_t*>(ptr);
+  uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST);
+  double current = bit_cast_value<double>(current_bits);
+  while (is_min ? (value < current) : (value > current)) {
+    uint64_t desired_bits = bit_cast_value<uint64_t>(value);
+    uint64_t expected_bits = current_bits;
+    if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits,
+          false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) {
+      return;
+    }
+    current_bits = expected_bits;
+    current = bit_cast_value<double>(current_bits);
+  }
+
+#elif defined(_MSC_VER)
+  auto* bits_ptr = reinterpret_cast<volatile long long*>(ptr);
+  long long current_bits = *bits_ptr;
+  double current = bit_cast_value<double>(current_bits);
+  while (is_min ? (value < current) : (value > current)) {
+    long long desired_bits = bit_cast_value<long long>(value);
+    long long old_bits =
+      _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits);
+    if (old_bits == current_bits) {
+      return;
+    }
+    current_bits = old_bits;
+    current = bit_cast_value<double>(current_bits);
+  }
+
+#else
+#error "No compare-and-swap implementation available for this compiler."
+#endif
+}
+
+inline void atomic_max_double(double* ptr, double value)
+{
+  atomic_update_double(ptr, value, false);
+}
+
+inline void atomic_min_double(double* ptr, double value)
+{
+  atomic_update_double(ptr, value, true);
+}
+
 namespace detail {
 
 //==============================================================================
@@ -147,7 +206,7 @@ namespace detail {
 //==============================================================================
 
 void MaterialVolumes::add_volume(
-  int index_elem, int index_material, double volume)
+  int index_elem, int index_material, double volume, const BoundingBox* bbox)
 {
   // This method handles adding elements to the materials hash table,
   // implementing open addressing with linear probing. Consistency across
@@ -166,10 +225,18 @@ void MaterialVolumes::add_volume(
     // Non-atomic read of current material
     int32_t current_val = *slot_ptr;
 
-    // Found the desired material; accumulate volume
+    // Found the desired material; accumulate volume and bbox
     if (current_val == index_material) {
 #pragma omp atomic
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);
+        atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);
+        atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);
+        atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);
+        atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);
+        atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
 
@@ -185,6 +252,14 @@ void MaterialVolumes::add_volume(
       if (claimed_slot || (expected_val == index_material)) {
 #pragma omp atomic
         this->volumes(index_elem, slot) += volume;
+        if (bbox) {
+          atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);
+          atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);
+          atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);
+          atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);
+          atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);
+          atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);
+        }
         return;
       }
     }
@@ -195,7 +270,7 @@ void MaterialVolumes::add_volume(
 }
 
 void MaterialVolumes::add_volume_unsafe(
-  int index_elem, int index_material, double volume)
+  int index_elem, int index_material, double volume, const BoundingBox* bbox)
 {
   // Linear probe
   for (int attempt = 0; attempt < table_size_; ++attempt) {
@@ -207,9 +282,23 @@ void MaterialVolumes::add_volume_unsafe(
     // Read current material
     int32_t current_val = this->materials(index_elem, slot);
 
-    // Found the desired material; accumulate volume
+    // Found the desired material; accumulate volume and bbox
     if (current_val == index_material) {
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        this->bboxes(index_elem, slot, 0) =
+          std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);
+        this->bboxes(index_elem, slot, 1) =
+          std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);
+        this->bboxes(index_elem, slot, 2) =
+          std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);
+        this->bboxes(index_elem, slot, 3) =
+          std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);
+        this->bboxes(index_elem, slot, 4) =
+          std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);
+        this->bboxes(index_elem, slot, 5) =
+          std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
 
@@ -217,6 +306,20 @@ void MaterialVolumes::add_volume_unsafe(
     if (current_val == EMPTY) {
       this->materials(index_elem, slot) = index_material;
       this->volumes(index_elem, slot) += volume;
+      if (bbox) {
+        this->bboxes(index_elem, slot, 0) =
+          std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);
+        this->bboxes(index_elem, slot, 1) =
+          std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);
+        this->bboxes(index_elem, slot, 2) =
+          std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);
+        this->bboxes(index_elem, slot, 3) =
+          std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);
+        this->bboxes(index_elem, slot, 4) =
+          std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);
+        this->bboxes(index_elem, slot, 5) =
+          std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);
+      }
       return;
     }
   }
@@ -333,6 +436,12 @@ vector<double> Mesh::volumes() const
 
 void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   int32_t* materials, double* volumes) const
+{
+  this->material_volumes(nx, ny, nz, table_size, materials, volumes, nullptr);
+}
+
+void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
+  int32_t* materials, double* volumes, double* bboxes) const
 {
   if (mpi::master) {
     header("MESH MATERIAL VOLUMES CALCULATION", 7);
@@ -351,7 +460,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   timer.start();
 
   // Create object for keeping track of materials/volumes
-  detail::MaterialVolumes result(materials, volumes, table_size);
+  detail::MaterialVolumes result(materials, volumes, bboxes, table_size);
+  bool compute_bboxes = bboxes != nullptr;
 
   // Determine bounding box
   auto bbox = this->bounding_box();
@@ -370,8 +480,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 #pragma omp parallel
   {
     // Preallocate vector for mesh indices and length fractions and particle
-    std::vector<int> bins;
-    std::vector<double> length_fractions;
+    vector<int> bins;
+    vector<double> length_fractions;
     Particle p;
 
     SourceSite site;
@@ -453,12 +563,36 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
             if (i_material != C_NONE) {
               i_material = model::materials[i_material]->id();
             }
+            double cumulative_frac = 0.0;
             for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
               int mesh_index = bins[i_bin];
               double length = distance * length_fractions[i_bin];
-
-              // Add volume to result
-              result.add_volume(mesh_index, i_material, length * d1 * d2);
+              double volume = length * d1 * d2;
+
+              if (compute_bboxes) {
+                double axis_start = r0[axis] + distance * cumulative_frac;
+                double axis_end = axis_start + length;
+                cumulative_frac += length_fractions[i_bin];
+
+                Position contrib_min = site.r;
+                Position contrib_max = site.r;
+
+                contrib_min[ax1] = site.r[ax1] - 0.5 * d1;
+                contrib_max[ax1] = site.r[ax1] + 0.5 * d1;
+                contrib_min[ax2] = site.r[ax2] - 0.5 * d2;
+                contrib_max[ax2] = site.r[ax2] + 0.5 * d2;
+                contrib_min[axis] = std::min(axis_start, axis_end);
+                contrib_max[axis] = std::max(axis_start, axis_end);
+
+                BoundingBox contrib_bbox {contrib_min, contrib_max};
+                contrib_bbox &= bbox;
+
+                result.add_volume(
+                  mesh_index, i_material, volume, &contrib_bbox);
+              } else {
+                // Add volume to result
+                result.add_volume(mesh_index, i_material, volume);
+              }
             }
 
             if (distance == max_distance)
@@ -505,10 +639,15 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   // Combine results from multiple MPI processes
   if (mpi::n_procs > 1) {
     int total = this->n_bins() * table_size;
+    int total_bbox = total * 6;
     if (mpi::master) {
       // Allocate temporary buffer for receiving data
-      std::vector<int32_t> mats(total);
-      std::vector<double> vols(total);
+      vector<int32_t> mats(total);
+      vector<double> vols(total);
+      vector<double> recv_bboxes;
+      if (compute_bboxes) {
+        recv_bboxes.resize(total_bbox);
+      }
 
       for (int i = 1; i < mpi::n_procs; ++i) {
         // Receive material indices and volumes from process i
@@ -516,6 +655,10 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           MPI_STATUS_IGNORE);
         MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,
           MPI_STATUS_IGNORE);
+        if (compute_bboxes) {
+          MPI_Recv(recv_bboxes.data(), total_bbox, MPI_DOUBLE, i, i,
+            mpi::intracomm, MPI_STATUS_IGNORE);
+        }
 
         // Combine with existing results; we can call thread unsafe version of
         // add_volume because each thread is operating on a different element
@@ -524,7 +667,18 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           for (int k = 0; k < table_size; ++k) {
             int index = index_elem * table_size + k;
             if (mats[index] != EMPTY) {
-              result.add_volume_unsafe(index_elem, mats[index], vols[index]);
+              if (compute_bboxes) {
+                int bbox_index = index * 6;
+                BoundingBox slot_bbox {
+                  {recv_bboxes[bbox_index + 0], recv_bboxes[bbox_index + 1],
+                    recv_bboxes[bbox_index + 2]},
+                  {recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4],
+                    recv_bboxes[bbox_index + 5]}};
+                result.add_volume_unsafe(
+                  index_elem, mats[index], vols[index], &slot_bbox);
+              } else {
+                result.add_volume_unsafe(index_elem, mats[index], vols[index]);
+              }
             }
           }
         }
@@ -533,6 +687,9 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
       // Send material indices and volumes to process 0
       MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm);
       MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
+      if (compute_bboxes) {
+        MPI_Send(bboxes, total_bbox, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
+      }
     }
   }
 
@@ -2428,14 +2585,14 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
 }
 
 extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny,
-  int nz, int table_size, int32_t* materials, double* volumes)
+  int nz, int table_size, int32_t* materials, double* volumes, double* bboxes)
 {
   if (int err = check_mesh(index))
     return err;
 
   try {
     model::meshes[index]->material_volumes(
-      nx, ny, nz, table_size, materials, volumes);
+      nx, ny, nz, table_size, materials, volumes, bboxes);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
     if (starts_with(e.what(), "Mesh")) {
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 9aca8b59656..c5855a7b051 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1,6 +1,8 @@
 from math import pi
 from tempfile import TemporaryDirectory
 from pathlib import Path
+import itertools
+import random
 
 import h5py
 import numpy as np
@@ -691,6 +693,49 @@ def test_mesh_material_volumes_serialize():
     assert new_volumes.by_element(3) == [(2, 1.0)]
 
 
+def test_mesh_material_volumes_serialize_with_bboxes():
+    materials = np.array([
+        [1, -1, -2],
+        [-1, -2, -2],
+        [2, 1, -2],
+        [2, -2, -2]
+    ])
+    volumes = np.array([
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0]
+    ])
+
+    # (xmin, ymin, zmin, xmax, ymax, zmax)
+    bboxes = np.empty((4, 3, 6))
+    bboxes[..., 0:3] = np.inf
+    bboxes[..., 3:6] = -np.inf
+    bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0]    # material 1
+    bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0]    # void
+    bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0]      # void
+    bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0]    # material 2
+    bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]       # material 1
+    bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0]    # material 2
+
+    mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes)
+    with TemporaryDirectory() as tmpdir:
+        path = f'{tmpdir}/volumes_bboxes.npz'
+        mmv.save(path)
+        loaded = openmc.MeshMaterialVolumes.from_npz(path)
+
+    assert loaded.has_bounding_boxes
+    first = loaded.by_element(0, include_bboxes=True)[0][2]
+    assert isinstance(first, openmc.BoundingBox)
+    np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0))
+    np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0))
+
+    second = loaded.by_element(0, include_bboxes=True)[1][2]
+    assert isinstance(second, openmc.BoundingBox)
+    np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0))
+    np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
+
+
 def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
     that overlaps with a vacuum boundary condition."""
@@ -718,6 +763,53 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model):
         assert evaluated[1] == pytest.approx(expected[1], rel=1e-2)
 
 
+def test_mesh_material_volumes_bounding_boxes():
+    # Create a model with 8 spherical cells at known locations with random radii
+    box = openmc.model.RectangularParallelepiped(
+        -10, 10, -10, 10, -10, 10, boundary_type='vacuum')
+
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+
+    sph_cells = []
+    for x, y, z in itertools.product((-5., 5.), repeat=3):
+        mat_i = mat.clone()
+        sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5))
+        sph_cells.append(openmc.Cell(region=-sph, fill=mat_i))
+    background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells]))
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry(sph_cells + [background])
+    model.settings.particles = 1000
+    model.settings.batches = 10
+
+    # Create a one-element mesh that encompasses the entire geometry
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-10., -10., -10.)
+    mesh.upper_right = (10., 10., 10.)
+    mesh.dimension = (1, 1, 1)
+
+    # Run material volume calculation with bounding boxes
+    n_samples = (400, 400, 400)
+    mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True)
+    assert mmv.has_bounding_boxes
+
+    # Create a mapping of material ID to bounding box
+    bbox_by_mat = {
+        mat_id: bbox
+        for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True)
+        if mat_id is not None and vol > 0.0
+    }
+
+    # Match the mesh ray spacing used for the bounding box estimator.
+    tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0]
+    for cell in sph_cells:
+        bbox = bbox_by_mat[cell.fill.id]
+        cell_bbox = cell.bounding_box
+        np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
+        np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
+
+
 def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell
     sphere = openmc.Sphere(r=100, boundary_type='vacuum')