DM-46239: Add pixelization envelopes for Intersection/UnionRegion

include/lsst/sphgeom/Region.h

@@ -86,7 +86,7 @@ class UnitVector3d;
 /// Region types outside of this library.
 class Region {
 public:
-    virtual ~Region() {}
+    virtual ~Region() = default;
 
     /// `clone` returns a deep copy of this region.
     virtual std::unique_ptr<Region> clone() const = 0;
@@ -148,6 +148,10 @@ class Region {
 
     static std::unique_ptr<Region> decode(std::uint8_t const * buffer, size_t n);
     ///@}
+
+    /// `getRegions` returns a vector of Region.
+    static std::vector<std::unique_ptr<Region>> getRegions(Region const &region);
+    ///@}
 };
 
 }} // namespace lsst::sphgeom

python/lsst/sphgeom/_region.cc

@@ -27,6 +27,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 #include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
 #include "pybind11/numpy.h"
 
 #include "lsst/sphgeom/python.h"
@@ -69,6 +70,7 @@ void defineClass(py::class_<Region, std::unique_ptr<Region>> &cls) {
             "region"_a);
     cls.def("encode", &python::encode);
     cls.def_static("decode", &python::decode<Region>, "bytes"_a);
+    cls.def_static("getRegions", Region::getRegions, "region"_a);
 }
 
 }  // sphgeom

src/PixelFinder.h

@@ -34,9 +34,10 @@
 /// \brief This file provides a base class for pixel finders.
 
 #include "lsst/sphgeom/RangeSet.h"
-
+#include "lsst/sphgeom/CompoundRegion.h"
 #include "ConvexPolygonImpl.h"
 
+#include <typeinfo>
 
 namespace lsst {
 namespace sphgeom {
@@ -156,23 +157,34 @@ template <
 >
 RangeSet findPixels(Region const & r, size_t maxRanges, int level) {
     RangeSet s;
-    Circle const * c = nullptr;
-    Ellipse const * e = nullptr;
-    Box const * b = nullptr;
-    if ((c = dynamic_cast<Circle const *>(&r))) {
-        Finder<Circle, InteriorOnly> find(s, *c, level, maxRanges);
+    if (auto circle = dynamic_cast<Circle const *>(&r)) {
+        Finder<Circle, InteriorOnly> find(s, *circle, level, maxRanges);
         find();
-    } else if ((e = dynamic_cast<Ellipse const *>(&r))) {
+    } else if (auto ellipse = dynamic_cast<Ellipse const *>(&r)) {
         Finder<Circle, InteriorOnly> find(
-            s, e->getBoundingCircle(), level, maxRanges);
+            s, ellipse->getBoundingCircle(), level, maxRanges);
         find();
-    } else if ((b = dynamic_cast<Box const *>(&r))) {
-        Finder<Box, InteriorOnly> find(s, *b, level, maxRanges);
+    } else if (auto box = dynamic_cast<Box const *>(&r)) {
+        Finder<Box, InteriorOnly> find(s, *box, level, maxRanges);
         find();
-    } else {
+    } else if (auto polygon = dynamic_cast<ConvexPolygon const *>(&r)) {
         Finder<ConvexPolygon, InteriorOnly> find(
-            s, dynamic_cast<ConvexPolygon const &>(r), level, maxRanges);
+                s, *polygon, level, maxRanges);
         find();
+    } else if (auto union_region = dynamic_cast<UnionRegion const *>(&r)) {
+        Region const &region1 = union_region->getOperand(0);
+        Region const &region2 = union_region->getOperand(1);
+        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("PixelFinder: Unsupported type ") + typeid(r).name());
     }
     return s;
 }

src/Region.cc

@@ -73,4 +73,21 @@ std::unique_ptr<Region> Region::decode(std::uint8_t const * buffer, size_t n) {
     throw std::runtime_error("Byte-string is not an encoded Region");
 }
 
+std::vector<std::unique_ptr<Region>> Region::getRegions(Region const &region) {
+    std::vector<std::unique_ptr<Region>> result;
+    if (auto union_region = dynamic_cast<UnionRegion const *>(&region)) {
+        for(int i = 0; i < 2; ++i) {
+            result.emplace_back(union_region->getOperand(i).clone());
+        }
+    } else if(auto intersection_region = dynamic_cast<IntersectionRegion const *>(&region)) {
+        for(int i = 0; i < 2; ++i) {
+            intersection_region->getOperand(i);
+            result.emplace_back(intersection_region->getOperand(i).clone());
+        }
+    } else {
+        result.emplace_back(region.clone());
+    }
+    return result;
+}
+
 }} // namespace lsst:sphgeom

tests/test_CompoundRegion.py

@@ -201,6 +201,23 @@ def testRelate(self):
         self.assertEqual(self.box.relate(self.instance), CONTAINS)
         self.assertEqual(self.faraway.relate(self.instance), DISJOINT)
 
+    def testGetRegion(self):
+        c1 = Circle(UnitVector3d(0.0, 0.0, 1.0), 1.0)
+        c2 = Circle(UnitVector3d(1.0, 0.0, 1.0), 2.0)
+        b1 = Box.fromDegrees(90, 0, 180, 45)
+        b2 = Box.fromDegrees(135, 15, 135, 30)
+        u1 = UnionRegion(c1, b1)
+        u2 = UnionRegion(c2, b2)
+        i1 = IntersectionRegion(c1, b1)
+        i2 = IntersectionRegion(c2, b2)
+        ur = UnionRegion(u1, u2)
+        ir = IntersectionRegion(i1, i2)
+        self.assertEqual(Region.getRegions(c1), [c1])
+        self.assertEqual(Region.getRegions(Region.getRegions(ir)[0]), [c1, b1])
+        self.assertEqual(Region.getRegions(Region.getRegions(ir)[1]), [c2, b2])
+        self.assertEqual(Region.getRegions(Region.getRegions(ur)[0]), [c1, b1])
+        self.assertEqual(Region.getRegions(Region.getRegions(ur)[1]), [c2, b2])
+
 
 if __name__ == "__main__":
     unittest.main()

tests/test_HtmPixelization.py

@@ -34,7 +34,16 @@
 
 import unittest
 
-from lsst.sphgeom import Angle, Circle, ConvexPolygon, HtmPixelization, RangeSet, UnitVector3d
+from lsst.sphgeom import (
+    Angle,
+    Circle,
+    ConvexPolygon,
+    HtmPixelization,
+    IntersectionRegion,
+    RangeSet,
+    UnionRegion,
+    UnitVector3d,
+)
 
 
 class HtmPixelizationTestCase(unittest.TestCase):
@@ -70,14 +79,46 @@ def test_level(self):
     def test_envelope_and_interior(self):
         pixelization = HtmPixelization(3)
         c = Circle(UnitVector3d(1, 1, 1), Angle.fromDegrees(0.1))
-        rs = pixelization.envelope(c)
-        self.assertTrue(rs == RangeSet(0x3FF))
+        rs1 = pixelization.envelope(c)
+        self.assertEqual(rs1, RangeSet(0x3FF))
         rs = pixelization.envelope(c, 1)
-        self.assertTrue(rs == RangeSet(0x3FF))
+        self.assertEqual(rs, RangeSet(0x3FF))
         self.assertTrue(rs.isWithin(pixelization.universe()))
         rs = pixelization.interior(c)
         self.assertTrue(rs.empty())
 
+        # Create a second region for a union.
+        s3 = UnitVector3d(1.0, -1.0, 1.0)  # Center of S3
+        c2 = Circle(s3, 1e-8)
+        rs2 = pixelization.envelope(c2)
+        self.assertEqual(rs2.ranges(), [(831, 832)])
+
+        # Try again with a union.
+        union = UnionRegion(c, c2)
+        rsu = pixelization.envelope(union)
+        self.assertEqual(rsu, rs1 | rs2)
+
+        # Check that nested unions also work.
+        c3 = Circle(s3, 2e-8)
+        union2 = UnionRegion(union, c3)
+        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):