diff --git a/src/PixelFinder.h b/src/PixelFinder.h
index b04dec5..9213990 100644
--- a/src/PixelFinder.h
+++ b/src/PixelFinder.h
@@ -177,8 +177,14 @@ RangeSet findPixels(Region const & r, size_t maxRanges, int level) {
         auto rs1 = findPixels<Finder, InteriorOnly>(region1, maxRanges, level);
         auto rs2 = findPixels<Finder, InteriorOnly>(region2, maxRanges, level);
         s = rs1.join(rs2);
+    }else if (auto intersection_region = dynamic_cast<IntersectionRegion const *>(&r)) {
+        Region const &region1 = intersection_region->getOperand(0);
+        Region const &region2 = intersection_region->getOperand(1);
+        auto rs1 = findPixels<Finder, InteriorOnly>(region1, maxRanges, level);
+        auto rs2 = findPixels<Finder, InteriorOnly>(region2, maxRanges, level);
+        s = rs1.intersection(rs2);
     } else {
-        throw std::runtime_error(std::string("Unsupported type ") + typeid(r).name());
+        throw std::runtime_error(std::string("PixelFinder: Unsupported type ") + typeid(r).name());
     }
     return s;
 }
diff --git a/tests/test_HtmPixelization.py b/tests/test_HtmPixelization.py
index a044934..d34fe33 100644
--- a/tests/test_HtmPixelization.py
+++ b/tests/test_HtmPixelization.py
@@ -34,7 +34,8 @@
 
 import unittest
 
-from lsst.sphgeom import Angle, Circle, ConvexPolygon, HtmPixelization, RangeSet, UnionRegion, UnitVector3d
+from lsst.sphgeom import Angle, Circle, ConvexPolygon, HtmPixelization, IntersectionRegion
+from lsst.sphgeom import RangeSet, UnionRegion, UnitVector3d
 
 
 class HtmPixelizationTestCase(unittest.TestCase):
@@ -95,6 +96,21 @@ def test_envelope_and_interior(self):
         rsu2 = pixelization.envelope(union2)
         self.assertEqual(rsu2, rsu2 | rsu)
 
+        # CHeck with intersection
+        c4 = Circle(UnitVector3d(1.0, 1.0, 2.0), 1)
+        c5 = Circle(UnitVector3d(1.0, 1.0, 2.5), 0.5)
+        rs4 = pixelization.envelope(c4)
+        rs5 = pixelization.envelope(c5)
+        intersection = IntersectionRegion(c4, c5)
+        rsi = pixelization.envelope(intersection)
+        self.assertEqual(rsi, rs4 & rs5)
+
+        # Check that nested intersection also work.
+        c6 = Circle(UnitVector3d(1.0, 1.0, 2.0), 2)
+        intersection2 = IntersectionRegion(intersection, c6)
+        rsi2 = pixelization.envelope(intersection2)
+        self.assertEqual(rsi2, rsi2 & rsi)
+
     def test_index_to_string(self):
         strings = ["S0", "S1", "S2", "S3", "N0", "N1", "N2", "N3"]
         for i in range(8, 16):