diff --git a/source/blender/blenkernel/BKE_pbvh_api.hh b/source/blender/blenkernel/BKE_pbvh_api.hh
index ba9306c8e7d..66e9e4ebbe1 100644
--- a/source/blender/blenkernel/BKE_pbvh_api.hh
+++ b/source/blender/blenkernel/BKE_pbvh_api.hh
@@ -13,6 +13,7 @@
 #include <string>
 
 #include "BLI_bit_group_vector.hh"
+#include "BLI_bounds_types.hh"
 #include "BLI_compiler_compat.h"
 #include "BLI_function_ref.hh"
 #include "BLI_index_mask.hh"
@@ -310,7 +311,7 @@ bool BKE_pbvh_has_faces(const PBVH *pbvh);
 /**
  * Get the PBVH root's bounding box.
  */
-void BKE_pbvh_bounding_box(const PBVH *pbvh, float min[3], float max[3]);
+blender::Bounds<blender::float3> BKE_pbvh_bounding_box(const PBVH *pbvh);
 
 void BKE_pbvh_sync_visibility_from_verts(PBVH *pbvh, Mesh *me);
 
@@ -391,8 +392,8 @@ blender::Vector<int> BKE_pbvh_node_calc_face_indices(const PBVH &pbvh, const PBV
  */
 int BKE_pbvh_num_faces(const PBVH *pbvh);
 
-void BKE_pbvh_node_get_BB(const PBVHNode *node, float bb_min[3], float bb_max[3]);
-void BKE_pbvh_node_get_original_BB(const PBVHNode *node, float bb_min[3], float bb_max[3]);
+blender::Bounds<blender::float3> BKE_pbvh_node_get_BB(const PBVHNode *node);
+blender::Bounds<blender::float3> BKE_pbvh_node_get_original_BB(const PBVHNode *node);
 
 float BKE_pbvh_node_get_tmin(const PBVHNode *node);
 
@@ -425,7 +426,7 @@ void BKE_pbvh_update_mask(PBVH *pbvh);
 void BKE_pbvh_update_vertex_data(PBVH *pbvh, int flags);
 void BKE_pbvh_update_visibility(PBVH *pbvh);
 void BKE_pbvh_update_normals(PBVH *pbvh, SubdivCCG *subdiv_ccg);
-void BKE_pbvh_redraw_BB(PBVH *pbvh, float bb_min[3], float bb_max[3]);
+blender::Bounds<blender::float3> BKE_pbvh_redraw_BB(PBVH *pbvh);
 blender::IndexMask BKE_pbvh_get_grid_updates(const PBVH *pbvh,
                                              blender::Span<const PBVHNode *> nodes,
                                              blender::IndexMaskMemory &memory);
diff --git a/source/blender/blenkernel/intern/pbvh.cc b/source/blender/blenkernel/intern/pbvh.cc
index 01a59c92bad..e8f344966a7 100644
--- a/source/blender/blenkernel/intern/pbvh.cc
+++ b/source/blender/blenkernel/intern/pbvh.cc
@@ -13,6 +13,7 @@
 #include "BLI_array_utils.hh"
 #include "BLI_bit_span_ops.hh"
 #include "BLI_bitmap.h"
+#include "BLI_bounds.hh"
 #include "BLI_math_geom.h"
 #include "BLI_math_matrix.h"
 #include "BLI_math_vector.h"
@@ -47,6 +48,7 @@
 #include "pbvh_intern.hh"
 
 using blender::BitGroupVector;
+using blender::Bounds;
 using blender::float3;
 using blender::MutableSpan;
 using blender::Span;
@@ -81,93 +83,44 @@ struct PBVHIter {
   int stackspace;
 };
 
-void BB_reset(BB *bb)
+/** Create invalid bounds for use with #math::min_max. */
+static Bounds<float3> negative_bounds()
 {
-  bb->bmin[0] = bb->bmin[1] = bb->bmin[2] = FLT_MAX;
-  bb->bmax[0] = bb->bmax[1] = bb->bmax[2] = -FLT_MAX;
-}
-
-void BB_expand(BB *bb, const float co[3])
-{
-  for (int i = 0; i < 3; i++) {
-    bb->bmin[i] = min_ff(bb->bmin[i], co[i]);
-    bb->bmax[i] = max_ff(bb->bmax[i], co[i]);
-  }
-}
-
-void BB_expand_with_bb(BB *bb, const BB *bb2)
-{
-  for (int i = 0; i < 3; i++) {
-    bb->bmin[i] = min_ff(bb->bmin[i], bb2->bmin[i]);
-    bb->bmax[i] = max_ff(bb->bmax[i], bb2->bmax[i]);
-  }
-}
-
-int BB_widest_axis(const BB *bb)
-{
-  float dim[3];
-
-  for (int i = 0; i < 3; i++) {
-    dim[i] = bb->bmax[i] - bb->bmin[i];
-  }
-
-  if (dim[0] > dim[1]) {
-    if (dim[0] > dim[2]) {
-      return 0;
-    }
-
-    return 2;
-  }
-
-  if (dim[1] > dim[2]) {
-    return 1;
-  }
-
-  return 2;
-}
-
-void BBC_update_centroid(BBC *bbc)
-{
-  for (int i = 0; i < 3; i++) {
-    bbc->bcentroid[i] = (bbc->bmin[i] + bbc->bmax[i]) * 0.5f;
-  }
+  return {float3(std::numeric_limits<float>::max()), float3(std::numeric_limits<float>::lowest())};
 }
 
 namespace blender::bke::pbvh {
 
 void update_node_bounds_mesh(const Span<float3> positions, PBVHNode &node)
 {
-  BB vb;
-  BB_reset(&vb);
+  Bounds<float3> bounds = negative_bounds();
   for (const int vert : node.vert_indices) {
-    BB_expand(&vb, positions[vert]);
+    math::min_max(positions[vert], bounds.min, bounds.max);
   }
-  node.vb = vb;
+  node.vb = bounds;
 }
 
 void update_node_bounds_grids(const CCGKey &key, const Span<CCGElem *> grids, PBVHNode &node)
 {
-  BB vb;
-  BB_reset(&vb);
+  Bounds<float3> bounds = negative_bounds();
   for (const int grid : node.prim_indices) {
     for (const int i : IndexRange(key.grid_area)) {
-      BB_expand(&vb, CCG_elem_offset_co(&key, grids[grid], i));
+      math::min_max(float3(CCG_elem_offset_co(&key, grids[grid], i)), bounds.min, bounds.max);
     }
   }
-  node.vb = vb;
+  node.vb = bounds;
 }
 
 void update_node_bounds_bmesh(PBVHNode &node)
 {
-  BB vb;
-  BB_reset(&vb);
+  Bounds<float3> bounds = negative_bounds();
   for (const BMVert *vert : node.bm_unique_verts) {
-    BB_expand(&vb, vert->co);
+    math::min_max(float3(vert->co), bounds.min, bounds.max);
   }
   for (const BMVert *vert : node.bm_other_verts) {
-    BB_expand(&vb, vert->co);
+    math::min_max(float3(vert->co), bounds.min, bounds.max);
   }
-  node.vb = vb;
+  node.vb = bounds;
 }
 
 }  // namespace blender::bke::pbvh
@@ -175,6 +128,7 @@ void update_node_bounds_bmesh(PBVHNode &node)
 /* Not recursive */
 static void update_node_vb(PBVH *pbvh, PBVHNode *node)
 {
+  using namespace blender;
   using namespace blender::bke::pbvh;
   if (node->flag & PBVH_Leaf) {
     switch (pbvh->header.type) {
@@ -190,24 +144,11 @@ static void update_node_vb(PBVH *pbvh, PBVHNode *node)
     }
   }
   else {
-    BB vb;
-    BB_reset(&vb);
-    BB_expand_with_bb(&vb, &pbvh->nodes[node->children_offset].vb);
-    BB_expand_with_bb(&vb, &pbvh->nodes[node->children_offset + 1].vb);
-    node->vb = vb;
+    node->vb = bounds::merge(pbvh->nodes[node->children_offset].vb,
+                             pbvh->nodes[node->children_offset + 1].vb);
   }
 }
 
-// void BKE_pbvh_node_BB_reset(PBVHNode *node)
-//{
-//  BB_reset(&node->vb);
-//}
-//
-// void BKE_pbvh_node_BB_expand(PBVHNode *node, float co[3])
-//{
-//  BB_expand(&node->vb, co);
-//}
-
 static bool face_materials_match(const int *material_indices,
                                  const bool *sharp_faces,
                                  const int a,
@@ -234,9 +175,10 @@ static int partition_prim_indices(blender::MutableSpan<int> prim_indices,
                                   int hi,
                                   int axis,
                                   float mid,
-                                  const Span<BBC> prim_bbc,
+                                  const Span<Bounds<float3>> prim_bounds,
                                   const Span<int> prim_to_face_map)
 {
+  using namespace blender;
   for (int i = lo; i < hi; i++) {
     prim_scratch[i - lo] = prim_indices[i];
   }
@@ -246,7 +188,8 @@ static int partition_prim_indices(blender::MutableSpan<int> prim_indices,
 
   while (i1 < hi) {
     const int face_i = prim_to_face_map[prim_scratch[i2]];
-    bool side = prim_bbc[prim_scratch[i2]].bcentroid[axis] >= mid;
+    const Bounds<float3> &bounds = prim_bounds[prim_scratch[i2]];
+    const bool side = math::midpoint(bounds.min[axis], bounds.max[axis]) >= mid;
 
     while (i1 < hi && prim_to_face_map[prim_scratch[i2]] == face_i) {
       prim_indices[side ? hi2-- : lo2++] = prim_scratch[i2];
@@ -371,11 +314,16 @@ static void build_mesh_leaf_node(const Span<int> corner_verts,
   BKE_pbvh_node_mark_rebuild_draw(node);
 }
 
-static void update_vb(PBVH *pbvh, PBVHNode *node, const Span<BBC> prim_bbc, int offset, int count)
+static void update_vb(const Span<int> prim_indices,
+                      PBVHNode *node,
+                      const Span<Bounds<float3>> prim_bounds,
+                      int offset,
+                      int count)
 {
-  BB_reset(&node->vb);
-  for (int i = offset + count - 1; i >= offset; i--) {
-    BB_expand_with_bb(&node->vb, (BB *)(&prim_bbc[pbvh->prim_indices[i]]));
+  using namespace blender;
+  node->vb = prim_bounds[prim_indices[offset]];
+  for (const int i : IndexRange(offset, count).drop_front(1)) {
+    node->vb = bounds::merge(node->vb, prim_bounds[prim_indices[i]]);
   }
   node->orig_vb = node->vb;
 }
@@ -432,7 +380,7 @@ static void build_leaf(PBVH *pbvh,
                        const Span<int> looptri_faces,
                        const bool *hide_poly,
                        int node_index,
-                       const Span<BBC> prim_bbc,
+                       const Span<Bounds<float3>> prim_bounds,
                        int offset,
                        int count)
 {
@@ -442,7 +390,7 @@ static void build_leaf(PBVH *pbvh,
   node.prim_indices = pbvh->prim_indices.as_span().slice(offset, count);
 
   /* Still need vb for searches */
-  update_vb(pbvh, &node, prim_bbc, offset, count);
+  update_vb(pbvh->prim_indices, &node, prim_bounds, offset, count);
 
   if (!pbvh->looptri.is_empty()) {
     build_mesh_leaf_node(
@@ -544,18 +492,18 @@ static void build_sub(PBVH *pbvh,
                       const int *material_indices,
                       const bool *sharp_faces,
                       int node_index,
-                      BB *cb,
-                      const Span<BBC> prim_bbc,
+                      const Bounds<float3> *cb,
+                      const Span<Bounds<float3>> prim_bounds,
                       int offset,
                       int count,
                       int *prim_scratch,
                       int depth)
 {
+  using namespace blender;
   const Span<int> prim_to_face_map = pbvh->header.type == PBVH_FACES ?
                                          looptri_faces :
                                          pbvh->subdiv_ccg->grid_to_face_map;
   int end;
-  BB cb_backing;
 
   if (!prim_scratch) {
     prim_scratch = static_cast<int *>(MEM_malloc_arrayN(pbvh->totprim, sizeof(int), __func__));
@@ -573,7 +521,7 @@ static void build_sub(PBVH *pbvh,
                  looptri_faces,
                  hide_poly,
                  node_index,
-                 prim_bbc,
+                 prim_bounds,
                  offset,
                  count);
 
@@ -590,18 +538,21 @@ static void build_sub(PBVH *pbvh,
   pbvh_grow_nodes(pbvh, pbvh->nodes.size() + 2);
 
   /* Update parent node bounding box */
-  update_vb(pbvh, &pbvh->nodes[node_index], prim_bbc, offset, count);
+  update_vb(pbvh->prim_indices, &pbvh->nodes[node_index], prim_bounds, offset, count);
 
+  Bounds<float3> cb_backing;
   if (!below_leaf_limit) {
     /* Find axis with widest range of primitive centroids */
     if (!cb) {
-      cb = &cb_backing;
-      BB_reset(cb);
+      cb_backing = negative_bounds();
       for (int i = offset + count - 1; i >= offset; i--) {
-        BB_expand(cb, prim_bbc[pbvh->prim_indices[i]].bcentroid);
+        const int prim = pbvh->prim_indices[i];
+        const float3 center = math::midpoint(prim_bounds[prim].min, prim_bounds[prim].max);
+        math::min_max(center, cb_backing.min, cb_backing.max);
       }
+      cb = &cb_backing;
     }
-    const int axis = BB_widest_axis(cb);
+    const int axis = math::dominant_axis(cb->max - cb->min);
 
     /* Partition primitives along that axis */
     end = partition_prim_indices(pbvh->prim_indices,
@@ -609,8 +560,8 @@ static void build_sub(PBVH *pbvh,
                                  offset,
                                  offset + count,
                                  axis,
-                                 (cb->bmax[axis] + cb->bmin[axis]) * 0.5f,
-                                 prim_bbc,
+                                 math::midpoint(cb->min[axis], cb->max[axis]),
+                                 prim_bounds,
                                  prim_to_face_map);
   }
   else {
@@ -633,7 +584,7 @@ static void build_sub(PBVH *pbvh,
             sharp_faces,
             pbvh->nodes[node_index].children_offset,
             nullptr,
-            prim_bbc,
+            prim_bounds,
             offset,
             end - offset,
             prim_scratch,
@@ -647,7 +598,7 @@ static void build_sub(PBVH *pbvh,
             sharp_faces,
             pbvh->nodes[node_index].children_offset + 1,
             nullptr,
-            prim_bbc,
+            prim_bounds,
             end,
             offset + count - end,
             prim_scratch,
@@ -665,8 +616,8 @@ static void pbvh_build(PBVH *pbvh,
                        const bool *hide_poly,
                        const int *material_indices,
                        const bool *sharp_faces,
-                       BB *cb,
-                       const Span<BBC> prim_bbc,
+                       const Bounds<float3> *cb,
+                       const Span<Bounds<float3>> prim_bounds,
                        int totprim)
 {
   if (totprim != pbvh->totprim) {
@@ -688,7 +639,7 @@ static void pbvh_build(PBVH *pbvh,
             sharp_faces,
             0,
             cb,
-            prim_bbc,
+            prim_bounds,
             0,
             totprim,
             nullptr,
@@ -779,6 +730,7 @@ void BKE_pbvh_update_mesh_pointers(PBVH *pbvh, Mesh *mesh)
 
 void BKE_pbvh_build_mesh(PBVH *pbvh, Mesh *mesh)
 {
+  using namespace blender;
   const int totvert = mesh->totvert;
   const int looptri_num = poly_to_tri_count(mesh->faces_num, mesh->totloop);
   MutableSpan<float3> vert_positions = mesh->vert_positions_for_write();
@@ -811,32 +763,25 @@ void BKE_pbvh_build_mesh(PBVH *pbvh, Mesh *mesh)
   pbvh->faces_num = mesh->faces_num;
 
   /* For each face, store the AABB and the AABB centroid */
-  blender::Array<BBC> prim_bbc(looptri_num);
-  BB cb;
-  BB_reset(&cb);
-  cb = blender::threading::parallel_reduce(
+  Array<Bounds<float3>> prim_bounds(looptri_num);
+  const Bounds<float3> cb = threading::parallel_reduce(
       looptris.index_range(),
       1024,
-      cb,
-      [&](const blender::IndexRange range, const BB &init) {
-        BB current = init;
+      negative_bounds(),
+      [&](const IndexRange range, const Bounds<float3> &init) {
+        Bounds<float3> current = init;
         for (const int i : range) {
           const MLoopTri &lt = looptris[i];
-          BBC *bbc = &prim_bbc[i];
-          BB_reset((BB *)bbc);
-          for (int j = 0; j < 3; j++) {
-            BB_expand((BB *)bbc, vert_positions[corner_verts[lt.tri[j]]]);
-          }
-          BBC_update_centroid(bbc);
-          BB_expand(&current, bbc->bcentroid);
+          Bounds<float3> &bounds = prim_bounds[i];
+          bounds = {vert_positions[corner_verts[lt.tri[0]]]};
+          math::min_max(vert_positions[corner_verts[lt.tri[1]]], bounds.min, bounds.max);
+          math::min_max(vert_positions[corner_verts[lt.tri[2]]], bounds.min, bounds.max);
+          const float3 center = math::midpoint(prim_bounds[i].min, prim_bounds[i].max);
+          math::min_max(center, current.min, current.max);
         }
         return current;
       },
-      [](const BB &a, const BB &b) {
-        BB current = a;
-        BB_expand_with_bb(&current, &b);
-        return current;
-      });
+      [](const Bounds<float3> &a, const Bounds<float3> &b) { return bounds::merge(a, b); });
 
   if (looptri_num) {
     const bool *hide_poly = static_cast<const bool *>(
@@ -853,7 +798,7 @@ void BKE_pbvh_build_mesh(PBVH *pbvh, Mesh *mesh)
                material_indices,
                sharp_faces,
                &cb,
-               prim_bbc,
+               prim_bounds,
                looptri_num);
 
 #ifdef TEST_PBVH_FACE_SPLIT
@@ -873,6 +818,7 @@ void BKE_pbvh_build_mesh(PBVH *pbvh, Mesh *mesh)
 
 void BKE_pbvh_build_grids(PBVH *pbvh, const CCGKey *key, Mesh *me, SubdivCCG *subdiv_ccg)
 {
+  using namespace blender;
   const int gridsize = key->grid_size;
   const Span<CCGElem *> grids = subdiv_ccg->grids;
 
@@ -898,32 +844,26 @@ void BKE_pbvh_build_grids(PBVH *pbvh, const CCGKey *key, Mesh *me, SubdivCCG *su
   pbvh->mesh = me;
 
   /* For each grid, store the AABB and the AABB centroid */
-  blender::Array<BBC> prim_bbc(grids.size());
-  BB cb;
-  BB_reset(&cb);
-  cb = blender::threading::parallel_reduce(
+  Array<Bounds<float3>> prim_bounds(grids.size());
+  const Bounds<float3> cb = threading::parallel_reduce(
       grids.index_range(),
       1024,
-      cb,
-      [&](const blender::IndexRange range, const BB &init) {
-        BB current = init;
+      negative_bounds(),
+      [&](const IndexRange range, const Bounds<float3> &init) {
+        Bounds<float3> current = init;
         for (const int i : range) {
           CCGElem *grid = grids[i];
-          BBC *bbc = &prim_bbc[i];
-          BB_reset((BB *)bbc);
-          for (int j = 0; j < gridsize * gridsize; j++) {
-            BB_expand((BB *)bbc, CCG_elem_offset_co(key, grid, j));
+          prim_bounds[i] = negative_bounds();
+          for (const int j : IndexRange(key->grid_area)) {
+            const float3 position = float3(CCG_elem_offset_co(key, grid, j));
+            math::min_max(position, prim_bounds[i].min, prim_bounds[i].max);
           }
-          BBC_update_centroid(bbc);
-          BB_expand(&current, bbc->bcentroid);
+          const float3 center = math::midpoint(prim_bounds[i].min, prim_bounds[i].max);
+          math::min_max(center, current.min, current.max);
         }
         return current;
       },
-      [](const BB &a, const BB &b) {
-        BB current = a;
-        BB_expand_with_bb(&current, &b);
-        return current;
-      });
+      [](const Bounds<float3> &a, const Bounds<float3> &b) { return bounds::merge(a, b); });
 
   if (!grids.is_empty()) {
     const int *material_indices = static_cast<const int *>(
@@ -931,7 +871,7 @@ void BKE_pbvh_build_grids(PBVH *pbvh, const CCGKey *key, Mesh *me, SubdivCCG *su
     const bool *sharp_faces = (const bool *)CustomData_get_layer_named(
         &me->face_data, CD_PROP_BOOL, "sharp_face");
     pbvh_build(
-        pbvh, {}, {}, {}, nullptr, material_indices, sharp_faces, &cb, prim_bbc, grids.size());
+        pbvh, {}, {}, {}, nullptr, material_indices, sharp_faces, &cb, prim_bounds, grids.size());
 
 #ifdef TEST_PBVH_FACE_SPLIT
     test_face_boundaries(pbvh);
@@ -1784,29 +1724,25 @@ void BKE_pbvh_update_visibility(PBVH *pbvh)
   }
 }
 
-void BKE_pbvh_redraw_BB(PBVH *pbvh, float bb_min[3], float bb_max[3])
+Bounds<float3> BKE_pbvh_redraw_BB(PBVH *pbvh)
 {
+  using namespace blender;
   if (pbvh->nodes.is_empty()) {
-    return;
+    return {};
   }
+  Bounds<float3> bounds = negative_bounds();
+
   PBVHIter iter;
-  PBVHNode *node;
-  BB bb;
-
-  BB_reset(&bb);
-
   pbvh_iter_begin(&iter, pbvh, {});
-
+  PBVHNode *node;
   while ((node = pbvh_iter_next(&iter, PBVH_Leaf))) {
     if (node->flag & PBVH_UpdateRedraw) {
-      BB_expand_with_bb(&bb, &node->vb);
+      bounds = bounds::merge(bounds, node->vb);
     }
   }
-
   pbvh_iter_end(&iter);
 
-  copy_v3_v3(bb_min, bb.bmin);
-  copy_v3_v3(bb_max, bb.bmax);
+  return bounds;
 }
 
 blender::IndexMask BKE_pbvh_get_grid_updates(const PBVH *pbvh,
@@ -1840,16 +1776,12 @@ bool BKE_pbvh_has_faces(const PBVH *pbvh)
   return (pbvh->totprim != 0);
 }
 
-void BKE_pbvh_bounding_box(const PBVH *pbvh, float min[3], float max[3])
+Bounds<float3> BKE_pbvh_bounding_box(const PBVH *pbvh)
 {
   if (pbvh->nodes.is_empty()) {
-    zero_v3(min);
-    zero_v3(max);
-    return;
+    return float3(0);
   }
-  const BB *bb = &pbvh->nodes[0].vb;
-  copy_v3_v3(min, bb->bmin);
-  copy_v3_v3(max, bb->bmax);
+  return pbvh->nodes.first().vb;
 }
 
 const CCGKey *BKE_pbvh_get_grid_key(const PBVH *pbvh)
@@ -2104,16 +2036,14 @@ Span<int> BKE_pbvh_node_get_grid_indices(const PBVHNode &node)
   return node.prim_indices;
 }
 
-void BKE_pbvh_node_get_BB(const PBVHNode *node, float bb_min[3], float bb_max[3])
+Bounds<float3> BKE_pbvh_node_get_BB(const PBVHNode *node)
 {
-  copy_v3_v3(bb_min, node->vb.bmin);
-  copy_v3_v3(bb_max, node->vb.bmax);
+  return node->vb;
 }
 
-void BKE_pbvh_node_get_original_BB(const PBVHNode *node, float bb_min[3], float bb_max[3])
+Bounds<float3> BKE_pbvh_node_get_original_BB(const PBVHNode *node)
 {
-  copy_v3_v3(bb_min, node->orig_vb.bmin);
-  copy_v3_v3(bb_max, node->orig_vb.bmax);
+  return node->orig_vb;
 }
 
 blender::MutableSpan<PBVHProxyNode> BKE_pbvh_node_get_proxies(PBVHNode *node)
@@ -2156,20 +2086,10 @@ struct RaycastData {
 
 static bool ray_aabb_intersect(PBVHNode *node, RaycastData *rcd)
 {
-  const float *bb_min, *bb_max;
-
   if (rcd->original) {
-    /* BKE_pbvh_node_get_original_BB */
-    bb_min = node->orig_vb.bmin;
-    bb_max = node->orig_vb.bmax;
+    return isect_ray_aabb_v3(&rcd->ray, node->orig_vb.min, node->orig_vb.max, &node->tmin);
   }
-  else {
-    /* BKE_pbvh_node_get_BB */
-    bb_min = node->vb.bmin;
-    bb_max = node->vb.bmax;
-  }
-
-  return isect_ray_aabb_v3(&rcd->ray, bb_min, bb_max, &node->tmin);
+  return isect_ray_aabb_v3(&rcd->ray, node->vb.min, node->vb.max, &node->tmin);
 }
 
 void BKE_pbvh_raycast(PBVH *pbvh,
@@ -2534,17 +2454,18 @@ void BKE_pbvh_clip_ray_ortho(
     return;
   }
   float rootmin_start, rootmin_end;
-  float bb_min_root[3], bb_max_root[3], bb_center[3], bb_diff[3];
+  Bounds<float3> bb_root;
+  float bb_center[3], bb_diff[3];
   IsectRayAABB_Precalc ray;
   float ray_normal_inv[3];
   float offset = 1.0f + 1e-3f;
   const float offset_vec[3] = {1e-3f, 1e-3f, 1e-3f};
 
   if (original) {
-    BKE_pbvh_node_get_original_BB(&pbvh->nodes.first(), bb_min_root, bb_max_root);
+    bb_root = BKE_pbvh_node_get_original_BB(&pbvh->nodes.first());
   }
   else {
-    BKE_pbvh_node_get_BB(&pbvh->nodes.first(), bb_min_root, bb_max_root);
+    bb_root = BKE_pbvh_node_get_BB(&pbvh->nodes.first());
   }
 
   /* Calc rough clipping to avoid overflow later. See #109555. */
@@ -2553,8 +2474,8 @@ void BKE_pbvh_clip_ray_ortho(
   float a[3], b[3], min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {FLT_MIN, FLT_MIN, FLT_MIN};
 
   /* Compute AABB bounds rotated along ray_normal.*/
-  copy_v3_v3(a, bb_min_root);
-  copy_v3_v3(b, bb_max_root);
+  copy_v3_v3(a, bb_root.min);
+  copy_v3_v3(b, bb_root.max);
   mul_m3_v3(mat, a);
   mul_m3_v3(mat, b);
   minmax_v3v3_v3(min, max, a);
@@ -2563,7 +2484,7 @@ void BKE_pbvh_clip_ray_ortho(
   float cent[3];
 
   /* Find midpoint of aabb on ray. */
-  mid_v3_v3v3(cent, bb_min_root, bb_max_root);
+  mid_v3_v3v3(cent, bb_root.min, bb_root.max);
   float t = line_point_factor_v3(cent, ray_start, ray_end);
   interp_v3_v3v3(cent, ray_start, ray_end, t);
 
@@ -2574,17 +2495,17 @@ void BKE_pbvh_clip_ray_ortho(
 
   /* Slightly offset min and max in case we have a zero width node
    * (due to a plane mesh for instance), or faces very close to the bounding box boundary. */
-  mid_v3_v3v3(bb_center, bb_max_root, bb_min_root);
+  mid_v3_v3v3(bb_center, bb_root.max, bb_root.min);
   /* Diff should be same for both min/max since it's calculated from center. */
-  sub_v3_v3v3(bb_diff, bb_max_root, bb_center);
+  sub_v3_v3v3(bb_diff, bb_root.max, bb_center);
   /* Handles case of zero width bb. */
   add_v3_v3(bb_diff, offset_vec);
-  madd_v3_v3v3fl(bb_max_root, bb_center, bb_diff, offset);
-  madd_v3_v3v3fl(bb_min_root, bb_center, bb_diff, -offset);
+  madd_v3_v3v3fl(bb_root.max, bb_center, bb_diff, offset);
+  madd_v3_v3v3fl(bb_root.min, bb_center, bb_diff, -offset);
 
   /* Final projection of start ray. */
   isect_ray_aabb_v3_precalc(&ray, ray_start, ray_normal);
-  if (!isect_ray_aabb_v3(&ray, bb_min_root, bb_max_root, &rootmin_start)) {
+  if (!isect_ray_aabb_v3(&ray, bb_root.min, bb_root.max, &rootmin_start)) {
     return;
   }
 
@@ -2592,7 +2513,7 @@ void BKE_pbvh_clip_ray_ortho(
   mul_v3_v3fl(ray_normal_inv, ray_normal, -1.0);
   isect_ray_aabb_v3_precalc(&ray, ray_end, ray_normal_inv);
   /* Unlikely to fail exiting if entering succeeded, still keep this here. */
-  if (!isect_ray_aabb_v3(&ray, bb_min_root, bb_max_root, &rootmin_end)) {
+  if (!isect_ray_aabb_v3(&ray, bb_root.min, bb_root.max, &rootmin_end)) {
     return;
   }
 
@@ -2626,13 +2547,13 @@ static bool nearest_to_ray_aabb_dist_sq(PBVHNode *node, FindNearestRayData *rcd)
 
   if (rcd->original) {
     /* BKE_pbvh_node_get_original_BB */
-    bb_min = node->orig_vb.bmin;
-    bb_max = node->orig_vb.bmax;
+    bb_min = node->orig_vb.min;
+    bb_max = node->orig_vb.max;
   }
   else {
     /* BKE_pbvh_node_get_BB */
-    bb_min = node->vb.bmin;
-    bb_max = node->vb.bmax;
+    bb_min = node->vb.min;
+    bb_max = node->vb.max;
   }
 
   float co_dummy[3], depth;
@@ -2811,8 +2732,7 @@ enum PlaneAABBIsect {
  * Returns true if the AABB is at least partially within the frustum
  * (ok, not a real frustum), false otherwise.
  */
-static PlaneAABBIsect test_frustum_aabb(const float bb_min[3],
-                                        const float bb_max[3],
+static PlaneAABBIsect test_frustum_aabb(const Bounds<float3> &bounds,
                                         const PBVHFrustumPlanes *frustum)
 {
   PlaneAABBIsect ret = ISECT_INSIDE;
@@ -2823,12 +2743,12 @@ static PlaneAABBIsect test_frustum_aabb(const float bb_min[3],
 
     for (int axis = 0; axis < 3; axis++) {
       if (planes[i][axis] < 0) {
-        vmin[axis] = bb_min[axis];
-        vmax[axis] = bb_max[axis];
+        vmin[axis] = bounds.min[axis];
+        vmax[axis] = bounds.max[axis];
       }
       else {
-        vmin[axis] = bb_max[axis];
-        vmax[axis] = bb_min[axis];
+        vmin[axis] = bounds.max[axis];
+        vmax[axis] = bounds.min[axis];
       }
     }
 
@@ -2845,22 +2765,12 @@ static PlaneAABBIsect test_frustum_aabb(const float bb_min[3],
 
 bool BKE_pbvh_node_frustum_contain_AABB(const PBVHNode *node, const PBVHFrustumPlanes *data)
 {
-  const float *bb_min, *bb_max;
-  /* BKE_pbvh_node_get_BB */
-  bb_min = node->vb.bmin;
-  bb_max = node->vb.bmax;
-
-  return test_frustum_aabb(bb_min, bb_max, data) != ISECT_OUTSIDE;
+  return test_frustum_aabb(node->vb, data) != ISECT_OUTSIDE;
 }
 
 bool BKE_pbvh_node_frustum_exclude_AABB(const PBVHNode *node, const PBVHFrustumPlanes *data)
 {
-  const float *bb_min, *bb_max;
-  /* BKE_pbvh_node_get_BB */
-  bb_min = node->vb.bmin;
-  bb_max = node->vb.bmax;
-
-  return test_frustum_aabb(bb_min, bb_max, data) != ISECT_INSIDE;
+  return test_frustum_aabb(node->vb, data) != ISECT_INSIDE;
 }
 
 void BKE_pbvh_update_normals(PBVH *pbvh, SubdivCCG *subdiv_ccg)
@@ -2979,7 +2889,7 @@ void BKE_pbvh_draw_debug_cb(PBVH *pbvh,
       continue;
     }
 
-    draw_fn(&node, user_data, node.vb.bmin, node.vb.bmax, node.flag);
+    draw_fn(&node, user_data, node.vb.min, node.vb.max, node.flag);
   }
 }
 
diff --git a/source/blender/blenkernel/intern/pbvh_bmesh.cc b/source/blender/blenkernel/intern/pbvh_bmesh.cc
index bb007d5a749..ec4515c1bd5 100644
--- a/source/blender/blenkernel/intern/pbvh_bmesh.cc
+++ b/source/blender/blenkernel/intern/pbvh_bmesh.cc
@@ -8,10 +8,12 @@
 
 #include "MEM_guardedalloc.h"
 
+#include "BLI_bounds.hh"
 #include "BLI_ghash.h"
 #include "BLI_heap_simple.h"
 #include "BLI_math_geom.h"
 #include "BLI_math_vector.h"
+#include "BLI_math_vector.hh"
 #include "BLI_memarena.h"
 #include "BLI_span.hh"
 #include "BLI_utildefines.h"
@@ -32,6 +34,8 @@
 static CLG_LogRef LOG = {"pbvh.bmesh"};
 
 using blender::Array;
+using blender::Bounds;
+using blender::float3;
 using blender::IndexRange;
 using blender::Span;
 using blender::Vector;
@@ -118,6 +122,11 @@ static void pbvh_bmesh_verify(PBVH *pbvh);
   } \
   ((void)0)
 
+static Bounds<float3> negative_bounds()
+{
+  return {float3(std::numeric_limits<float>::max()), float3(std::numeric_limits<float>::lowest())};
+}
+
 static std::array<BMEdge *, 3> bm_edges_from_tri(BMesh *bm, const blender::Span<BMVert *> v_tri)
 {
   return {
@@ -205,10 +214,11 @@ static void pbvh_bmesh_node_finalize(PBVH *pbvh,
                                      const int cd_vert_node_offset,
                                      const int cd_face_node_offset)
 {
+  using namespace blender;
   PBVHNode *n = &pbvh->nodes[node_index];
   bool has_visible = false;
 
-  BB_reset(&n->vb);
+  n->vb = negative_bounds();
 
   for (BMFace *f : n->bm_faces) {
     /* Update ownership of faces. */
@@ -230,7 +240,7 @@ static void pbvh_bmesh_node_finalize(PBVH *pbvh,
         }
       }
       /* Update node bounding box. */
-      BB_expand(&n->vb, v->co);
+      math::min_max(float3(v->co), n->vb.min, n->vb.max);
     } while ((l_iter = l_iter->next) != l_first);
 
     if (!BM_elem_flag_test(f, BM_ELEM_HIDDEN)) {
@@ -251,8 +261,11 @@ static void pbvh_bmesh_node_finalize(PBVH *pbvh,
 }
 
 /** Recursively split the node if it exceeds the leaf_limit. */
-static void pbvh_bmesh_node_split(PBVH *pbvh, const Span<BBC> bbc_array, int node_index)
+static void pbvh_bmesh_node_split(PBVH *pbvh,
+                                  const Span<Bounds<float3>> face_bounds,
+                                  int node_index)
 {
+  using namespace blender;
   const int cd_vert_node_offset = pbvh->cd_vert_node_offset;
   const int cd_face_node_offset = pbvh->cd_face_node_offset;
   PBVHNode *n = &pbvh->nodes[node_index];
@@ -264,17 +277,16 @@ static void pbvh_bmesh_node_split(PBVH *pbvh, const Span<BBC> bbc_array, int nod
   }
 
   /* Calculate bounding box around primitive centroids. */
-  BB cb;
-  BB_reset(&cb);
+  Bounds<float3> cb = negative_bounds();
   for (BMFace *f : n->bm_faces) {
-    const BBC *bbc = &bbc_array[BM_elem_index_get(f)];
-
-    BB_expand(&cb, bbc->bcentroid);
+    const int i = BM_elem_index_get(f);
+    const float3 center = math::midpoint(face_bounds[i].min, face_bounds[i].max);
+    math::min_max(center, cb.min, cb.max);
   }
 
   /* Find widest axis and its midpoint. */
-  const int axis = BB_widest_axis(&cb);
-  const float mid = (cb.bmax[axis] + cb.bmin[axis]) * 0.5f;
+  const int axis = math::dominant_axis(cb.max - cb.min);
+  const float mid = math::midpoint(cb.max[axis], cb.min[axis]);
 
   /* Add two new child nodes. */
   const int children = pbvh->nodes.size();
@@ -293,9 +305,8 @@ static void pbvh_bmesh_node_split(PBVH *pbvh, const Span<BBC> bbc_array, int nod
 
   /* Partition the parent node's faces between the two children. */
   for (BMFace *f : n->bm_faces) {
-    const BBC *bbc = &bbc_array[BM_elem_index_get(f)];
-
-    if (bbc->bcentroid[axis] < mid) {
+    const int i = BM_elem_index_get(f);
+    if (math::midpoint(face_bounds[i].min[axis], face_bounds[i].max[axis]) < mid) {
       c1->bm_faces.add(f);
     }
     else {
@@ -347,22 +358,22 @@ static void pbvh_bmesh_node_split(PBVH *pbvh, const Span<BBC> bbc_array, int nod
   n->flag &= ~PBVH_Leaf;
 
   /* Recurse. */
-  pbvh_bmesh_node_split(pbvh, bbc_array, children);
-  pbvh_bmesh_node_split(pbvh, bbc_array, children + 1);
+  pbvh_bmesh_node_split(pbvh, face_bounds, children);
+  pbvh_bmesh_node_split(pbvh, face_bounds, children + 1);
 
   /* Array maybe reallocated, update current node pointer */
   n = &pbvh->nodes[node_index];
 
   /* Update bounding box. */
-  BB_reset(&n->vb);
-  BB_expand_with_bb(&n->vb, &pbvh->nodes[n->children_offset].vb);
-  BB_expand_with_bb(&n->vb, &pbvh->nodes[n->children_offset + 1].vb);
+  n->vb = bounds::merge(pbvh->nodes[n->children_offset].vb,
+                        pbvh->nodes[n->children_offset + 1].vb);
   n->orig_vb = n->vb;
 }
 
 /** Recursively split the node if it exceeds the leaf_limit. */
 static bool pbvh_bmesh_node_limit_ensure(PBVH *pbvh, int node_index)
 {
+  using namespace blender;
   PBVHNode &node = pbvh->nodes[node_index];
   const int faces_num = node.bm_faces.size();
   if (faces_num <= pbvh->leaf_limit) {
@@ -374,21 +385,19 @@ static bool pbvh_bmesh_node_limit_ensure(PBVH *pbvh, int node_index)
   pbvh->draw_cache_invalid = true;
 
   /* For each BMFace, store the AABB and AABB centroid. */
-  Array<BBC> bbc_array(faces_num);
+  Array<Bounds<float3>> face_bounds(faces_num);
 
   int i = 0;
   for (BMFace *f : node.bm_faces) {
-    BBC *bbc = &bbc_array[i];
+    face_bounds[i] = negative_bounds();
 
-    BB_reset((BB *)bbc);
     BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
     BMLoop *l_iter = l_first;
     do {
-      BB_expand((BB *)bbc, l_iter->v->co);
+      math::min_max(float3(l_iter->v->co), face_bounds[i].min, face_bounds[i].max);
     } while ((l_iter = l_iter->next) != l_first);
-    BBC_update_centroid(bbc);
 
-    /* So we can do direct lookups on 'bbc_array'. */
+    /* So we can do direct lookups on 'face_bounds'. */
     BM_elem_index_set(f, i); /* set_dirty! */
     i++;
   }
@@ -396,7 +405,7 @@ static bool pbvh_bmesh_node_limit_ensure(PBVH *pbvh, int node_index)
   /* Likely this is already dirty. */
   pbvh->header.bm->elem_index_dirty |= BM_FACE;
 
-  pbvh_bmesh_node_split(pbvh, bbc_array, node_index);
+  pbvh_bmesh_node_split(pbvh, face_bounds, node_index);
 
   return true;
 }
@@ -1955,9 +1964,13 @@ struct FastNodeBuildInfo {
  * This function is multi-thread-able since each invocation applies
  * to a sub part of the arrays.
  */
-static void pbvh_bmesh_node_limit_ensure_fast(
-    PBVH *pbvh, BMFace **nodeinfo, BBC *bbc_array, FastNodeBuildInfo *node, MemArena *arena)
+static void pbvh_bmesh_node_limit_ensure_fast(PBVH *pbvh,
+                                              BMFace **nodeinfo,
+                                              Bounds<float3> *face_bounds,
+                                              FastNodeBuildInfo *node,
+                                              MemArena *arena)
 {
+  using namespace blender;
   FastNodeBuildInfo *child1, *child2;
 
   if (node->totface <= pbvh->leaf_limit) {
@@ -1965,20 +1978,19 @@ static void pbvh_bmesh_node_limit_ensure_fast(
   }
 
   /* Calculate bounding box around primitive centroids. */
-  BB cb;
-  BB_reset(&cb);
+  Bounds<float3> cb = negative_bounds();
   for (int i = 0; i < node->totface; i++) {
     BMFace *f = nodeinfo[i + node->start];
-    BBC *bbc = &bbc_array[BM_elem_index_get(f)];
-
-    BB_expand(&cb, bbc->bcentroid);
+    const int face_index = BM_elem_index_get(f);
+    const float3 center = math::midpoint(face_bounds[face_index].min, face_bounds[face_index].max);
+    math::min_max(center, cb.min, cb.max);
   }
 
   /* Initialize the children. */
 
   /* Find widest axis and its midpoint. */
-  const int axis = BB_widest_axis(&cb);
-  const float mid = (cb.bmax[axis] + cb.bmin[axis]) * 0.5f;
+  const int axis = math::dominant_axis(cb.max - cb.min);
+  const float mid = math::midpoint(cb.max[axis], cb.min[axis]);
 
   int num_child1 = 0, num_child2 = 0;
 
@@ -1986,17 +1998,18 @@ static void pbvh_bmesh_node_limit_ensure_fast(
   const int end = node->start + node->totface;
   for (int i = node->start; i < end - num_child2; i++) {
     BMFace *f = nodeinfo[i];
-    BBC *bbc = &bbc_array[BM_elem_index_get(f)];
+    const int face_i = BM_elem_index_get(f);
 
-    if (bbc->bcentroid[axis] > mid) {
+    if (math::midpoint(face_bounds[face_i].min[axis], face_bounds[face_i].max[axis]) > mid) {
       int i_iter = end - num_child2 - 1;
       int candidate = -1;
       /* Found a face that should be part of another node, look for a face to substitute with. */
 
       for (; i_iter > i; i_iter--) {
         BMFace *f_iter = nodeinfo[i_iter];
-        const BBC *bbc_iter = &bbc_array[BM_elem_index_get(f_iter)];
-        if (bbc_iter->bcentroid[axis] <= mid) {
+        const int face_iter_i = BM_elem_index_get(f_iter);
+        if (math::midpoint(face_bounds[face_iter_i].min[axis],
+                           face_bounds[face_iter_i].max[axis]) <= mid) {
           candidate = i_iter;
           break;
         }
@@ -2048,13 +2061,17 @@ static void pbvh_bmesh_node_limit_ensure_fast(
   child2->start = node->start + num_child1;
   child1->child1 = child1->child2 = child2->child1 = child2->child2 = nullptr;
 
-  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, bbc_array, child1, arena);
-  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, bbc_array, child2, arena);
+  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, face_bounds, child1, arena);
+  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, face_bounds, child2, arena);
 }
 
-static void pbvh_bmesh_create_nodes_fast_recursive(
-    PBVH *pbvh, BMFace **nodeinfo, BBC *bbc_array, FastNodeBuildInfo *node, int node_index)
+static void pbvh_bmesh_create_nodes_fast_recursive(PBVH *pbvh,
+                                                   BMFace **nodeinfo,
+                                                   Bounds<float3> *face_bounds,
+                                                   FastNodeBuildInfo *node,
+                                                   int node_index)
 {
+  using namespace blender;
   PBVHNode *n = &pbvh->nodes[node_index];
   /* Two cases, node does not have children or does have children. */
   if (node->child1) {
@@ -2063,16 +2080,15 @@ static void pbvh_bmesh_create_nodes_fast_recursive(
     n->children_offset = children_offset;
     pbvh_grow_nodes(pbvh, pbvh->nodes.size() + 2);
     pbvh_bmesh_create_nodes_fast_recursive(
-        pbvh, nodeinfo, bbc_array, node->child1, children_offset);
+        pbvh, nodeinfo, face_bounds, node->child1, children_offset);
     pbvh_bmesh_create_nodes_fast_recursive(
-        pbvh, nodeinfo, bbc_array, node->child2, children_offset + 1);
+        pbvh, nodeinfo, face_bounds, node->child2, children_offset + 1);
 
     n = &pbvh->nodes[node_index];
 
     /* Update bounding box. */
-    BB_reset(&n->vb);
-    BB_expand_with_bb(&n->vb, &pbvh->nodes[n->children_offset].vb);
-    BB_expand_with_bb(&n->vb, &pbvh->nodes[n->children_offset + 1].vb);
+    n->vb = bounds::merge(pbvh->nodes[n->children_offset].vb,
+                          pbvh->nodes[n->children_offset + 1].vb);
     n->orig_vb = n->vb;
   }
   else {
@@ -2086,13 +2102,12 @@ static void pbvh_bmesh_create_nodes_fast_recursive(
     n->flag = PBVH_Leaf;
     n->bm_faces.reserve(node->totface);
 
-    BB_reset(&n->vb);
+    n->vb = face_bounds[node->start];
 
     const int end = node->start + node->totface;
 
     for (int i = node->start; i < end; i++) {
       BMFace *f = nodeinfo[i];
-      BBC *bbc = &bbc_array[BM_elem_index_get(f)];
 
       /* Update ownership of faces. */
       n->bm_faces.add_new(f);
@@ -2119,7 +2134,7 @@ static void pbvh_bmesh_create_nodes_fast_recursive(
         has_visible = true;
       }
 
-      BB_expand_with_bb(&n->vb, (BB *)bbc);
+      n->vb = bounds::merge(n->vb, face_bounds[BM_elem_index_get(f)]);
     }
 
     BLI_assert(n->vb.bmin[0] <= n->vb.bmax[0] && n->vb.bmin[1] <= n->vb.bmax[1] &&
@@ -2149,6 +2164,7 @@ void BKE_pbvh_build_bmesh(PBVH *pbvh,
                           const int cd_vert_node_offset,
                           const int cd_face_node_offset)
 {
+  using namespace blender;
   pbvh->header.bm = bm;
 
   BKE_pbvh_bmesh_detail_size_set(pbvh, 0.75);
@@ -2162,7 +2178,8 @@ void BKE_pbvh_build_bmesh(PBVH *pbvh,
   BKE_pbvh_update_bmesh_offsets(pbvh, cd_vert_node_offset, cd_face_node_offset);
 
   /* bounding box array of all faces, no need to recalculate every time. */
-  BBC *bbc_array = static_cast<BBC *>(MEM_mallocN(sizeof(BBC) * bm->totface, "BBC"));
+  Bounds<float3> *face_bounds = static_cast<Bounds<float3> *>(
+      MEM_mallocN(sizeof(Bounds<float3>) * bm->totface, "Bounds<float3>"));
   BMFace **nodeinfo = static_cast<BMFace **>(
       MEM_mallocN(sizeof(*nodeinfo) * bm->totface, "nodeinfo"));
   MemArena *arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "fast PBVH node storage");
@@ -2171,17 +2188,15 @@ void BKE_pbvh_build_bmesh(PBVH *pbvh,
   BMFace *f;
   int i;
   BM_ITER_MESH_INDEX (f, &iter, bm, BM_FACES_OF_MESH, i) {
-    BBC *bbc = &bbc_array[i];
+    face_bounds[i] = negative_bounds();
+
     BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
     BMLoop *l_iter = l_first;
-
-    BB_reset((BB *)bbc);
     do {
-      BB_expand((BB *)bbc, l_iter->v->co);
+      math::min_max(float3(l_iter->v->co), face_bounds[i].min, face_bounds[i].max);
     } while ((l_iter = l_iter->next) != l_first);
-    BBC_update_centroid(bbc);
 
-    /* so we can do direct lookups on 'bbc_array' */
+    /* so we can do direct lookups on 'face_bounds' */
     BM_elem_index_set(f, i); /* set_dirty! */
     nodeinfo[i] = f;
     BM_ELEM_CD_SET_INT(f, cd_face_node_offset, DYNTOPO_NODE_NONE);
@@ -2199,7 +2214,7 @@ void BKE_pbvh_build_bmesh(PBVH *pbvh,
   rootnode.totface = bm->totface;
 
   /* Start recursion, assign faces to nodes accordingly. */
-  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, bbc_array, &rootnode, arena);
+  pbvh_bmesh_node_limit_ensure_fast(pbvh, nodeinfo, face_bounds, &rootnode, arena);
 
   /* We now have all faces assigned to a node,
    * next we need to assign those to the gsets of the nodes. */
@@ -2208,10 +2223,10 @@ void BKE_pbvh_build_bmesh(PBVH *pbvh,
   pbvh->nodes.append({});
 
   /* Take root node and visit and populate children recursively. */
-  pbvh_bmesh_create_nodes_fast_recursive(pbvh, nodeinfo, bbc_array, &rootnode, 0);
+  pbvh_bmesh_create_nodes_fast_recursive(pbvh, nodeinfo, face_bounds, &rootnode, 0);
 
   BLI_memarena_free(arena);
-  MEM_freeN(bbc_array);
+  MEM_freeN(face_bounds);
   MEM_freeN(nodeinfo);
 }
 
diff --git a/source/blender/blenkernel/intern/pbvh_intern.hh b/source/blender/blenkernel/intern/pbvh_intern.hh
index a7fe9552e73..d1937f9b8a9 100644
--- a/source/blender/blenkernel/intern/pbvh_intern.hh
+++ b/source/blender/blenkernel/intern/pbvh_intern.hh
@@ -5,6 +5,7 @@
 #pragma once
 
 #include "BLI_array.hh"
+#include "BLI_bounds_types.hh"
 #include "BLI_math_vector_types.hh"
 #include "BLI_set.hh"
 #include "BLI_span.hh"
@@ -25,16 +26,6 @@ struct MLoopTri;
 struct BMVert;
 struct BMFace;
 
-/* Axis-aligned bounding box */
-struct BB {
-  float bmin[3], bmax[3];
-};
-
-/* Axis-aligned bounding box with centroid */
-struct BBC {
-  float bmin[3], bmax[3], bcentroid[3];
-};
-
 /* NOTE: this structure is getting large, might want to split it into
  * union'd structs */
 struct PBVHNode {
@@ -42,8 +33,8 @@ struct PBVHNode {
   blender::draw::pbvh::PBVHBatches *draw_batches = nullptr;
 
   /* Voxel bounds */
-  BB vb = {};
-  BB orig_vb = {};
+  blender::Bounds<blender::float3> vb = {};
+  blender::Bounds<blender::float3> orig_vb = {};
 
   /* For internal nodes, the offset of the children in the PBVH
    * 'nodes' array. */
@@ -218,20 +209,6 @@ struct PBVH {
 
 /* pbvh.cc */
 
-void BB_reset(BB *bb);
-/**
- * Expand the bounding box to include a new coordinate.
- */
-void BB_expand(BB *bb, const float co[3]);
-/**
- * Expand the bounding box to include another bounding box.
- */
-void BB_expand_with_bb(BB *bb, const BB *bb2);
-void BBC_update_centroid(BBC *bbc);
-/**
- * Return 0, 1, or 2 to indicate the widest axis of the bounding box.
- */
-int BB_widest_axis(const BB *bb);
 void pbvh_grow_nodes(PBVH *bvh, int totnode);
 bool ray_face_intersection_quad(const float ray_start[3],
                                 IsectRayPrecalc *isect_precalc,
diff --git a/source/blender/blenkernel/intern/pbvh_pixels.cc b/source/blender/blenkernel/intern/pbvh_pixels.cc
index e1f088967a7..94b79009e57 100644
--- a/source/blender/blenkernel/intern/pbvh_pixels.cc
+++ b/source/blender/blenkernel/intern/pbvh_pixels.cc
@@ -105,10 +105,9 @@ static void split_thread_job(TaskPool *__restrict pool, void *taskdata);
 static void split_pixel_node(
     PBVH *pbvh, SplitNodePair *split, Image *image, ImageUser *image_user, SplitQueueData *tdata)
 {
-  BB cb;
   PBVHNode *node = &split->node;
 
-  cb = node->vb;
+  const Bounds<float3> cb = node->vb;
 
   if (count_node_pixels(*node) <= pbvh->pixel_leaf_limit || split->depth >= pbvh->depth_limit) {
     BKE_pbvh_pixels_node_data_get(split->node).rebuild_undo_regions();
@@ -116,8 +115,8 @@ static void split_pixel_node(
   }
 
   /* Find widest axis and its midpoint */
-  const int axis = BB_widest_axis(&cb);
-  const float mid = (cb.bmax[axis] + cb.bmin[axis]) * 0.5f;
+  const int axis = math::dominant_axis(cb.max - cb.min);
+  const float mid = (cb.max[axis] + cb.min[axis]) * 0.5f;
 
   node->flag = (PBVHNodeFlags)(int(node->flag) & int(~PBVH_TexLeaf));
 
@@ -134,10 +133,10 @@ static void split_pixel_node(
   child2->flag = PBVH_TexLeaf;
 
   child1->vb = cb;
-  child1->vb.bmax[axis] = mid;
+  child1->vb.max[axis] = mid;
 
   child2->vb = cb;
-  child2->vb.bmin[axis] = mid;
+  child2->vb.min[axis] = mid;
 
   NodeData &data = BKE_pbvh_pixels_node_data_get(split->node);
 
diff --git a/source/blender/blenlib/BLI_bounds_types.hh b/source/blender/blenlib/BLI_bounds_types.hh
index fb627a70453..5613c5433b4 100644
--- a/source/blender/blenlib/BLI_bounds_types.hh
+++ b/source/blender/blenlib/BLI_bounds_types.hh
@@ -13,6 +13,9 @@ namespace blender {
 template<typename T> struct Bounds {
   T min;
   T max;
+  Bounds() = default;
+  Bounds(const T &value) : min(value), max(value) {}
+  Bounds(const T &min, const T &max) : min(min), max(max) {}
 };
 
 }  // namespace blender
diff --git a/source/blender/draw/engines/eevee_next/eevee_sync.cc b/source/blender/draw/engines/eevee_next/eevee_sync.cc
index 2b1c2e40d26..abd820d2b84 100644
--- a/source/blender/draw/engines/eevee_next/eevee_sync.cc
+++ b/source/blender/draw/engines/eevee_next/eevee_sync.cc
@@ -230,11 +230,9 @@ bool SyncModule::sync_sculpt(Object *ob,
 
   /* Use a valid bounding box. The PBVH module already does its own culling, but a valid */
   /* bounding box is still needed for directional shadow tile-map bounds computation. */
-  float3 min, max;
-  BKE_pbvh_bounding_box(ob_ref.object->sculpt->pbvh, min, max);
-  float3 center = (min + max) * 0.5;
-  float3 half_extent = max - center;
-  half_extent += inflate_bounds;
+  const Bounds<float3> bounds = BKE_pbvh_bounding_box(ob_ref.object->sculpt->pbvh);
+  const float3 center = math::midpoint(bounds.min, bounds.max);
+  const float3 half_extent = bounds.max - center + inflate_bounds;
   inst_.manager->update_handle_bounds(res_handle, center, half_extent);
 
   inst_.manager->extract_object_attributes(res_handle, ob_ref, material_array.gpu_materials);
diff --git a/source/blender/editors/sculpt_paint/sculpt.cc b/source/blender/editors/sculpt_paint/sculpt.cc
index 343c010549b..8131e6af53f 100644
--- a/source/blender/editors/sculpt_paint/sculpt.cc
+++ b/source/blender/editors/sculpt_paint/sculpt.cc
@@ -1652,17 +1652,16 @@ static void sculpt_extend_redraw_rect_previous(Object *ob, rcti *rect)
 
 bool SCULPT_get_redraw_rect(ARegion *region, RegionView3D *rv3d, Object *ob, rcti *rect)
 {
+  using namespace blender;
   PBVH *pbvh = ob->sculpt->pbvh;
-  float bb_min[3], bb_max[3];
-
   if (!pbvh) {
     return false;
   }
 
-  BKE_pbvh_redraw_BB(pbvh, bb_min, bb_max);
+  const Bounds<float3> bounds = BKE_pbvh_redraw_BB(pbvh);
 
   /* Convert 3D bounding box to screen space. */
-  if (!paint_convert_bb_to_rect(rect, bb_min, bb_max, region, rv3d, ob)) {
+  if (!paint_convert_bb_to_rect(rect, bounds.min, bounds.max, region, rv3d, ob)) {
     return false;
   }
 
@@ -2718,6 +2717,7 @@ void SCULPT_calc_vertex_displacement(SculptSession *ss,
 
 bool SCULPT_search_sphere(PBVHNode *node, SculptSearchSphereData *data)
 {
+  using namespace blender;
   const float *center;
   float nearest[3];
   if (data->center) {
@@ -2726,7 +2726,7 @@ bool SCULPT_search_sphere(PBVHNode *node, SculptSearchSphereData *data)
   else {
     center = data->ss->cache ? data->ss->cache->location : data->ss->cursor_location;
   }
-  float t[3], bb_min[3], bb_max[3];
+  float t[3];
 
   if (data->ignore_fully_ineffective) {
     if (BKE_pbvh_node_fully_hidden_get(node)) {
@@ -2737,19 +2737,15 @@ bool SCULPT_search_sphere(PBVHNode *node, SculptSearchSphereData *data)
     }
   }
 
-  if (data->original) {
-    BKE_pbvh_node_get_original_BB(node, bb_min, bb_max);
-  }
-  else {
-    BKE_pbvh_node_get_BB(node, bb_min, bb_max);
-  }
+  const Bounds<float3> bounds = (data->original) ? BKE_pbvh_node_get_original_BB(node) :
+                                                   BKE_pbvh_node_get_BB(node);
 
   for (int i = 0; i < 3; i++) {
-    if (bb_min[i] > center[i]) {
-      nearest[i] = bb_min[i];
+    if (bounds.min[i] > center[i]) {
+      nearest[i] = bounds.min[i];
     }
-    else if (bb_max[i] < center[i]) {
-      nearest[i] = bb_max[i];
+    else if (bounds.max[i] < center[i]) {
+      nearest[i] = bounds.max[i];
     }
     else {
       nearest[i] = center[i];
@@ -2763,24 +2759,19 @@ bool SCULPT_search_sphere(PBVHNode *node, SculptSearchSphereData *data)
 
 bool SCULPT_search_circle(PBVHNode *node, SculptSearchCircleData *data)
 {
-  float bb_min[3], bb_max[3];
-
+  using namespace blender;
   if (data->ignore_fully_ineffective) {
     if (BKE_pbvh_node_fully_masked_get(node)) {
       return false;
     }
   }
 
-  if (data->original) {
-    BKE_pbvh_node_get_original_BB(node, bb_min, bb_max);
-  }
-  else {
-    BKE_pbvh_node_get_BB(node, bb_min, bb_min);
-  }
+  const Bounds<float3> bounds = (data->original) ? BKE_pbvh_node_get_original_BB(node) :
+                                                   BKE_pbvh_node_get_BB(node);
 
   float dummy_co[3], dummy_depth;
   const float dist_sq = dist_squared_ray_to_aabb_v3(
-      data->dist_ray_to_aabb_precalc, bb_min, bb_max, dummy_co, &dummy_depth);
+      data->dist_ray_to_aabb_precalc, bounds.min, bounds.max, dummy_co, &dummy_depth);
 
   /* Seems like debug code.
    * Maybe this function can just return true if the node is not fully masked. */
@@ -5459,10 +5450,7 @@ void SCULPT_flush_update_step(bContext *C, SculptUpdateType update_flags)
        * the object's evaluated geometry bounding box is necessary because sculpt strokes don't
        * cause an object reevaluation. */
       BKE_mesh_tag_positions_changed_no_normals(mesh);
-
-      Bounds<float3> bounds;
-      BKE_pbvh_bounding_box(ob->sculpt->pbvh, bounds.min, bounds.max);
-      mesh->bounds_set_eager(bounds);
+      mesh->bounds_set_eager(BKE_pbvh_bounding_box(ob->sculpt->pbvh));
       if (ob->runtime->bounds_eval) {
         ob->runtime->bounds_eval = mesh->bounds_min_max();
       }
diff --git a/source/blender/editors/sculpt_paint/sculpt_detail.cc b/source/blender/editors/sculpt_paint/sculpt_detail.cc
index 4e9789e410c..1b8611f32de 100644
--- a/source/blender/editors/sculpt_paint/sculpt_detail.cc
+++ b/source/blender/editors/sculpt_paint/sculpt_detail.cc
@@ -11,6 +11,7 @@
 #include "BLI_blenlib.h"
 #include "BLI_math_matrix.h"
 #include "BLI_math_rotation.h"
+#include "BLI_math_vector.hh"
 
 #include "BLT_translation.h"
 
@@ -85,11 +86,10 @@ static bool sculpt_and_dynamic_topology_poll(bContext *C)
 
 static int sculpt_detail_flood_fill_exec(bContext *C, wmOperator *op)
 {
+  using namespace blender;
   Sculpt *sd = CTX_data_tool_settings(C)->sculpt;
   Object *ob = CTX_data_active_object(C);
   SculptSession *ss = ob->sculpt;
-  float size;
-  float bb_min[3], bb_max[3], center[3], dim[3];
 
   blender::Vector<PBVHNode *> nodes = blender::bke::pbvh::search_gather(ss->pbvh, {});
 
@@ -101,11 +101,10 @@ static int sculpt_detail_flood_fill_exec(bContext *C, wmOperator *op)
     BKE_pbvh_node_mark_topology_update(node);
   }
   /* Get the bounding box, its center and size. */
-  BKE_pbvh_bounding_box(ob->sculpt->pbvh, bb_min, bb_max);
-  add_v3_v3v3(center, bb_min, bb_max);
-  mul_v3_fl(center, 0.5f);
-  sub_v3_v3v3(dim, bb_max, bb_min);
-  size = max_fff(dim[0], dim[1], dim[2]);
+  const Bounds<float3> bounds = BKE_pbvh_bounding_box(ob->sculpt->pbvh);
+  const float3 center = math::midpoint(bounds.min, bounds.max);
+  const float3 dim = bounds.max - bounds.min;
+  const float size = math::reduce_max(dim);
 
   /* Update topology size. */
   float object_space_constant_detail = 1.0f /
diff --git a/source/blender/editors/sculpt_paint/sculpt_ops.cc b/source/blender/editors/sculpt_paint/sculpt_ops.cc
index 61d597ccb0f..203c16153e4 100644
--- a/source/blender/editors/sculpt_paint/sculpt_ops.cc
+++ b/source/blender/editors/sculpt_paint/sculpt_ops.cc
@@ -11,6 +11,8 @@
 
 #include "BLI_ghash.h"
 #include "BLI_gsqueue.h"
+#include "BLI_math_matrix.hh"
+#include "BLI_math_vector.hh"
 #include "BLI_task.h"
 #include "BLI_utildefines.h"
 
@@ -310,6 +312,7 @@ static void sculpt_init_session(Main *bmain, Depsgraph *depsgraph, Scene *scene,
 
 void SCULPT_ensure_valid_pivot(const Object *ob, Scene *scene)
 {
+  using namespace blender;
   UnifiedPaintSettings *ups = &scene->toolsettings->unified_paint_settings;
   const SculptSession *ss = ob->sculpt;
 
@@ -320,11 +323,9 @@ void SCULPT_ensure_valid_pivot(const Object *ob, Scene *scene)
 
   /* No valid pivot? Use bounding box center. */
   if (ups->average_stroke_counter == 0 || !ups->last_stroke_valid) {
-    float location[3], max[3];
-    BKE_pbvh_bounding_box(ss->pbvh, location, max);
-
-    interp_v3_v3v3(location, location, max, 0.5f);
-    mul_m4_v3(ob->object_to_world, location);
+    const Bounds<float3> bounds = BKE_pbvh_bounding_box(ob->sculpt->pbvh);
+    const float3 center = math::midpoint(bounds.min, bounds.max);
+    const float3 location = math::transform_point(float4x4(ob->object_to_world), center);
 
     copy_v3_v3(ups->average_stroke_accum, location);
     ups->average_stroke_counter = 1;