diff --git a/src/scenic/core/regions.py b/src/scenic/core/regions.py
index c6b5b33f4..33fc987f8 100644
--- a/src/scenic/core/regions.py
+++ b/src/scenic/core/regions.py
@@ -1013,6 +1013,10 @@ def AABB(self):
 
     @cached_property
     def _transform(self):
+        """Transform from input mesh to final mesh.
+
+        :meta private:
+        """
         if self.dimensions is not None:
             scale = numpy.array(self.dimensions) / self._mesh.extents
         else:
@@ -1024,9 +1028,33 @@ def _transform(self):
         transform = compose_matrix(scale=scale, angles=angles, translate=self.position)
         return transform
 
+    @cached_property
+    def _shapeTransform(self):
+        """Transform from Shape mesh (scaled to unit dimensions) to final mesh.
+
+        :meta private:
+        """
+        if self.dimensions is not None:
+            scale = numpy.array(self.dimensions)
+        else:
+            scale = self._mesh.extents
+        if self.rotation is not None:
+            angles = self.rotation._trimeshEulerAngles()
+        else:
+            angles = None
+        transform = compose_matrix(scale=scale, angles=angles, translate=self.position)
+        return transform
+
     @cached_property
     def _boundingPolygonHull(self):
         assert not isLazy(self)
+        if self._shape:
+            raw = self._shape._multipoint
+            tr = self._shapeTransform
+            matrix = numpy.concatenate((tr[0:3, 0:3].flatten(), tr[0:3, 3]))
+            transformed = shapely.affinity.affine_transform(raw, matrix)
+            return transformed.convex_hull
+
         return shapely.multipoints(self.mesh.vertices).convex_hull
 
     @cached_property
@@ -1848,7 +1876,8 @@ def _circumradius(self):
         if self._scaledShape:
             return self._scaledShape._circumradius
         if self._shape:
-            scale = max(self.dimensions) if self.dimensions else 1
+            dims = self.dimensions or self._mesh.extents
+            scale = max(dims)
             return scale * self._shape._circumradius
 
         return numpy.max(numpy.linalg.norm(self.mesh.vertices, axis=1))
@@ -1863,7 +1892,7 @@ def _interiorPoint(self):
         if self._shape:
             raw = self._shape._interiorPoint
             homog = numpy.append(raw, [1])
-            return numpy.dot(self._transform, homog)[:3]
+            return numpy.dot(self._shapeTransform, homog)[:3]
 
         return findMeshInteriorPoint(self.mesh, num_samples=self.num_samples)
 
diff --git a/src/scenic/core/shapes.py b/src/scenic/core/shapes.py
index 6dc224388..db1c04deb 100644
--- a/src/scenic/core/shapes.py
+++ b/src/scenic/core/shapes.py
@@ -3,6 +3,7 @@
 from abc import ABC, abstractmethod
 
 import numpy
+import shapely
 import trimesh
 from trimesh.transformations import (
     concatenate_matrices,
@@ -72,6 +73,10 @@ def _circumradius(self):
     def _interiorPoint(self):
         return findMeshInteriorPoint(self.mesh)
 
+    @cached_property
+    def _multipoint(self):
+        return shapely.multipoints(self.mesh.vertices)
+
 
 ###################################################################################################
 # 3D Shape Classes
diff --git a/tests/core/test_regions.py b/tests/core/test_regions.py
index 580224508..03b832e6e 100644
--- a/tests/core/test_regions.py
+++ b/tests/core/test_regions.py
@@ -14,6 +14,11 @@
 from tests.utils import deprecationTest, sampleSceneFrom
 
 
+def assertPolygonsEqual(p1, p2, prec=1e-6):
+    assert p1.difference(p2).area == pytest.approx(0, abs=prec)
+    assert p2.difference(p1).area == pytest.approx(0, abs=prec)
+
+
 def sample_ignoring_rejections(region, num_samples):
     samples = []
     for _ in range(num_samples):
@@ -340,7 +345,23 @@ def test_mesh_intersects():
     assert not r1.getSurfaceRegion().intersects(r2.getSurfaceRegion())
 
 
-def test_mesh_circumradius(getAssetPath):
+def test_mesh_boundingPolygon(getAssetPath, pytestconfig):
+    r = BoxRegion(dimensions=(8, 6, 2)).difference(BoxRegion(dimensions=(2, 2, 3)))
+    bp = r.boundingPolygon
+    poly = shapely.geometry.Polygon(
+        [(-4, 3), (4, 3), (4, -3), (-4, -3)], [[(-1, 1), (1, 1), (1, -1), (-1, -1)]]
+    )
+    assertPolygonsEqual(bp.polygons, poly)
+
+    shape = MeshShape(BoxRegion(dimensions=(1, 2, 3)).mesh)
+    o = Orientation.fromEuler(0, 0, math.pi / 4)
+    r = MeshVolumeRegion(shape.mesh, dimensions=(2, 4, 2), rotation=o, _shape=shape)
+    bp = r.boundingPolygon
+    sr2 = math.sqrt(2)
+    poly = shapely.geometry.Polygon([(-sr2, 2), (sr2, 2), (sr2, -2), (-sr2, -2)])
+    assertPolygonsEqual(bp.polygons, poly)
+
+    samples = 50 if pytestconfig.getoption("--fast") else 200
     r1 = BoxRegion(dimensions=(1, 2, 3), position=(4, 5, 6))
     bo = Orientation.fromEuler(math.pi / 4, math.pi / 4, math.pi / 4)
     r2 = MeshVolumeRegion(r1.mesh, position=(15, 20, 5), rotation=bo, _scaledShape=r1)
@@ -349,6 +370,29 @@ def test_mesh_circumradius(getAssetPath):
     shape = MeshShape.fromFile(planePath)
     r4 = MeshVolumeRegion(shape.mesh, dimensions=(0.5, 2, 1.5), rotation=bo, _shape=shape)
     for reg in (r1, r2, r3, r4):
+        bp = reg.boundingPolygon
+        pts = trimesh.sample.volume_mesh(reg.mesh, samples)
+        assert all(bp.containsPoint(pt) for pt in pts)
+
+
+def test_mesh_circumradius(getAssetPath):
+    r1 = BoxRegion(dimensions=(1, 2, 3), position=(4, 5, 6))
+
+    bo = Orientation.fromEuler(math.pi / 4, math.pi / 4, math.pi / 4)
+    r2 = MeshVolumeRegion(r1.mesh, position=(15, 20, 5), rotation=bo, _scaledShape=r1)
+
+    planePath = getAssetPath("meshes/classic_plane.obj.bz2")
+    r3 = MeshVolumeRegion.fromFile(planePath, dimensions=(20, 20, 10))
+
+    shape = MeshShape.fromFile(planePath)
+    r4 = MeshVolumeRegion(shape.mesh, dimensions=(0.5, 2, 1.5), rotation=bo, _shape=shape)
+
+    r = BoxRegion(dimensions=(1, 2, 3)).difference(BoxRegion(dimensions=(0.5, 1, 1)))
+    shape = MeshShape(r.mesh)
+    scaled = MeshVolumeRegion(shape.mesh, dimensions=(6, 5, 4)).mesh
+    r5 = MeshVolumeRegion(scaled, position=(-10, -5, 30), rotation=bo, _shape=shape)
+
+    for reg in (r1, r2, r3, r4, r5):
         pos = reg.position
         d = 2.01 * reg._circumradius
         assert SpheroidRegion(dimensions=(d, d, d), position=pos).containsRegion(reg)
@@ -374,6 +418,12 @@ def test_mesh_interiorPoint():
     r3 = MeshVolumeRegion(shape.mesh, position=(-10, -5, 30), rotation=bo, _shape=shape)
     regions.append(r3)
 
+    r = BoxRegion(dimensions=(1, 2, 3)).difference(BoxRegion(dimensions=(0.5, 1, 1)))
+    shape = MeshShape(r.mesh)
+    scaled = MeshVolumeRegion(shape.mesh, dimensions=(0.1, 0.1, 0.1)).mesh
+    r4 = MeshVolumeRegion(scaled, position=(-10, -5, 30), rotation=bo, _shape=shape)
+    regions.append(r4)
+
     for reg in regions:
         cp = reg._interiorPoint
         # N.B. _containsPointExact can fail with embreex installed!