diff --git a/include/lsst/sphgeom/Region.h b/include/lsst/sphgeom/Region.h index d58725e..f142b90 100644 --- a/include/lsst/sphgeom/Region.h +++ b/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 clone() const = 0; @@ -148,6 +148,10 @@ class Region { static std::unique_ptr decode(std::uint8_t const * buffer, size_t n); ///@} + + /// `getRegions` returns a vector of Region. + static std::vector> getRegions(Region const ®ion); + ///@} }; }} // namespace lsst::sphgeom diff --git a/python/lsst/sphgeom/_region.cc b/python/lsst/sphgeom/_region.cc index 4f133d2..11cb15b 100644 --- a/python/lsst/sphgeom/_region.cc +++ b/python/lsst/sphgeom/_region.cc @@ -27,6 +27,7 @@ * along with this program. If not, see . */ #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_> &cls) { "region"_a); cls.def("encode", &python::encode); cls.def_static("decode", &python::decode, "bytes"_a); + cls.def_static("getRegions", Region::getRegions, "region"_a); } } // sphgeom diff --git a/src/PixelFinder.h b/src/PixelFinder.h index 942652d..718d184 100644 --- a/src/PixelFinder.h +++ b/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 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(&r))) { - Finder find(s, *c, level, maxRanges); + if (auto circle = dynamic_cast(&r)) { + Finder find(s, *circle, level, maxRanges); find(); - } else if ((e = dynamic_cast(&r))) { + } else if (auto ellipse = dynamic_cast(&r)) { Finder find( - s, e->getBoundingCircle(), level, maxRanges); + s, ellipse->getBoundingCircle(), level, maxRanges); find(); - } else if ((b = dynamic_cast(&r))) { - Finder find(s, *b, level, maxRanges); + } else if (auto box = dynamic_cast(&r)) { + Finder find(s, *box, level, maxRanges); find(); - } else { + } else if (auto polygon = dynamic_cast(&r)) { Finder find( - s, dynamic_cast(r), level, maxRanges); + s, *polygon, level, maxRanges); find(); + } else if (auto union_region = dynamic_cast(&r)) { + Region const ®ion1 = union_region->getOperand(0); + Region const ®ion2 = union_region->getOperand(1); + auto rs1 = findPixels(region1, maxRanges, level); + auto rs2 = findPixels(region2, maxRanges, level); + s = rs1.join(rs2); + } else if (auto intersection_region = dynamic_cast(&r)) { + Region const ®ion1 = intersection_region->getOperand(0); + Region const ®ion2 = intersection_region->getOperand(1); + auto rs1 = findPixels(region1, maxRanges, level); + auto rs2 = findPixels(region2, maxRanges, level); + s = rs1.intersection(rs2); + } else { + throw std::runtime_error(std::string("PixelFinder: Unsupported type ") + typeid(r).name()); } return s; } diff --git a/src/Region.cc b/src/Region.cc index 37d6cfb..7361149 100644 --- a/src/Region.cc +++ b/src/Region.cc @@ -73,4 +73,21 @@ std::unique_ptr Region::decode(std::uint8_t const * buffer, size_t n) { throw std::runtime_error("Byte-string is not an encoded Region"); } +std::vector> Region::getRegions(Region const ®ion) { + std::vector> result; + if (auto union_region = dynamic_cast(®ion)) { + for(int i = 0; i < 2; ++i) { + result.emplace_back(union_region->getOperand(i).clone()); + } + } else if(auto intersection_region = dynamic_cast(®ion)) { + 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 diff --git a/tests/test_CompoundRegion.py b/tests/test_CompoundRegion.py index da5ad5d..eb62dca 100644 --- a/tests/test_CompoundRegion.py +++ b/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() diff --git a/tests/test_HtmPixelization.py b/tests/test_HtmPixelization.py index 8393d51..f7d17ae 100644 --- a/tests/test_HtmPixelization.py +++ b/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):