From ea35072a1f199db05d6a62884862001202eefe58 Mon Sep 17 00:00:00 2001 From: Tim Jenness Date: Tue, 10 Sep 2024 07:37:35 -0700 Subject: [PATCH 1/6] Add a failing tests for HtmPixelization envelope with UnionRegion --- tests/test_HtmPixelization.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/tests/test_HtmPixelization.py b/tests/test_HtmPixelization.py index 8393d51c..8d8fb879 100644 --- a/tests/test_HtmPixelization.py +++ b/tests/test_HtmPixelization.py @@ -34,7 +34,7 @@ import unittest -from lsst.sphgeom import Angle, Circle, ConvexPolygon, HtmPixelization, RangeSet, UnitVector3d +from lsst.sphgeom import Angle, Circle, ConvexPolygon, HtmPixelization, RangeSet, UnionRegion, UnitVector3d class HtmPixelizationTestCase(unittest.TestCase): @@ -70,14 +70,25 @@ 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) + def test_index_to_string(self): strings = ["S0", "S1", "S2", "S3", "N0", "N1", "N2", "N3"] for i in range(8, 16): From fb504128cde3d68da0ea8cd03769993175aea782 Mon Sep 17 00:00:00 2001 From: Matthias Wittgen Date: Tue, 1 Oct 2024 16:06:16 -0700 Subject: [PATCH 2/6] Add UnionRegion handling in PixelFinder --- src/PixelFinder.h | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/PixelFinder.h b/src/PixelFinder.h index 942652d6..b04dec59 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,28 @@ 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 { + throw std::runtime_error(std::string("Unsupported type ") + typeid(r).name()); } return s; } From 75a7240222c3578b9223053e6f3052e4af49483d Mon Sep 17 00:00:00 2001 From: Tim Jenness Date: Wed, 2 Oct 2024 08:23:17 -0700 Subject: [PATCH 3/6] Test pixelization envelope with nested unions --- tests/test_HtmPixelization.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_HtmPixelization.py b/tests/test_HtmPixelization.py index 8d8fb879..a0449348 100644 --- a/tests/test_HtmPixelization.py +++ b/tests/test_HtmPixelization.py @@ -89,6 +89,12 @@ def test_envelope_and_interior(self): 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) + def test_index_to_string(self): strings = ["S0", "S1", "S2", "S3", "N0", "N1", "N2", "N3"] for i in range(8, 16): From 18ca8276a8b9ae15f64522ce805c148c24f12f95 Mon Sep 17 00:00:00 2001 From: Matthias Wittgen Date: Thu, 3 Oct 2024 15:07:38 -0700 Subject: [PATCH 4/6] Add support for IntersectionRegion --- src/PixelFinder.h | 8 +++++++- tests/test_HtmPixelization.py | 26 +++++++++++++++++++++++++- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/PixelFinder.h b/src/PixelFinder.h index b04dec59..92139902 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(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("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 a0449348..6f282eba 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, UnionRegion, UnitVector3d +from lsst.sphgeom import ( + Angle, + Circle, + ConvexPolygon, + HtmPixelization, + IntersectionRegion, + RangeSet, + UnionRegion, + UnitVector3d, +) class HtmPixelizationTestCase(unittest.TestCase): @@ -95,6 +104,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): From d98c7808629637d132fd9fc847430ff072526784 Mon Sep 17 00:00:00 2001 From: Matthias Wittgen Date: Sat, 12 Oct 2024 13:56:33 -0700 Subject: [PATCH 5/6] Add unit test for flatten --- include/lsst/sphgeom/CompoundRegion.h | 1 + include/lsst/sphgeom/Region.h | 6 +++++- python/lsst/sphgeom/_region.cc | 6 ++++++ src/Region.cc | 15 +++++++++++++++ tests/test_CompoundRegion.py | 11 +++++++++++ 5 files changed, 38 insertions(+), 1 deletion(-) diff --git a/include/lsst/sphgeom/CompoundRegion.h b/include/lsst/sphgeom/CompoundRegion.h index c4b6067a..0dbf47a3 100644 --- a/include/lsst/sphgeom/CompoundRegion.h +++ b/include/lsst/sphgeom/CompoundRegion.h @@ -88,6 +88,7 @@ class CompoundRegion : public Region { } static std::unique_ptr decode(std::uint8_t const *buffer, size_t n); ///@} + static void flatten(Region const ®ion, std::vector> &result); protected: diff --git a/include/lsst/sphgeom/Region.h b/include/lsst/sphgeom/Region.h index d58725ef..1283490f 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); ///@} + + /// `flatten` flattens a nested Region to a list of regions. + static void flatten(Region const ®ion, std::vector> &result); + ///@} }; }} // namespace lsst::sphgeom diff --git a/python/lsst/sphgeom/_region.cc b/python/lsst/sphgeom/_region.cc index 4f133d26..783330cd 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,11 @@ void defineClass(py::class_> &cls) { "region"_a); cls.def("encode", &python::encode); cls.def_static("decode", &python::decode, "bytes"_a); + cls.def_static("flatten", [](Region const ®ion) { + std::vector> result; + Region::flatten(region, result); + return result; + }, "region"_a); } } // sphgeom diff --git a/src/Region.cc b/src/Region.cc index 37d6cfb8..b8e32107 100644 --- a/src/Region.cc +++ b/src/Region.cc @@ -72,5 +72,20 @@ std::unique_ptr Region::decode(std::uint8_t const * buffer, size_t n) { } throw std::runtime_error("Byte-string is not an encoded Region"); } + void Region::flatten(Region const ®ion, std::vector> &result) { + if (auto union_region = dynamic_cast(®ion)) { + for(int i = 0; i < 2; ++i) { + union_region->getOperand(i); + flatten(union_region->getOperand(i), result); + } + } else if(auto intersection_region = dynamic_cast(®ion)) { + for(int i = 0; i < 2; ++i) { + intersection_region->getOperand(i); + flatten(intersection_region->getOperand(i), result); + } + } else { + result.emplace_back(region.clone()); + } + } }} // namespace lsst:sphgeom diff --git a/tests/test_CompoundRegion.py b/tests/test_CompoundRegion.py index da5ad5d7..3e3418b1 100644 --- a/tests/test_CompoundRegion.py +++ b/tests/test_CompoundRegion.py @@ -201,6 +201,17 @@ def testRelate(self): self.assertEqual(self.box.relate(self.instance), CONTAINS) self.assertEqual(self.faraway.relate(self.instance), DISJOINT) + def testFlatten(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) + c = UnionRegion(u1, u2) + d = Region.flatten(c) + self.assertEqual(d, [c1, b1, c2, b2]) + if __name__ == "__main__": unittest.main() From d008b965b9b1fd30f3207d7befdde763c135ef2b Mon Sep 17 00:00:00 2001 From: Matthias Wittgen Date: Mon, 14 Oct 2024 15:09:13 -0700 Subject: [PATCH 6/6] Add Region.getRegions() --- include/lsst/sphgeom/CompoundRegion.h | 1 - include/lsst/sphgeom/Region.h | 4 ++-- python/lsst/sphgeom/_region.cc | 6 +----- src/PixelFinder.h | 2 +- src/Region.cc | 28 ++++++++++++++------------- tests/test_CompoundRegion.py | 14 ++++++++++---- tests/test_HtmPixelization.py | 2 +- 7 files changed, 30 insertions(+), 27 deletions(-) diff --git a/include/lsst/sphgeom/CompoundRegion.h b/include/lsst/sphgeom/CompoundRegion.h index 0dbf47a3..c4b6067a 100644 --- a/include/lsst/sphgeom/CompoundRegion.h +++ b/include/lsst/sphgeom/CompoundRegion.h @@ -88,7 +88,6 @@ class CompoundRegion : public Region { } static std::unique_ptr decode(std::uint8_t const *buffer, size_t n); ///@} - static void flatten(Region const ®ion, std::vector> &result); protected: diff --git a/include/lsst/sphgeom/Region.h b/include/lsst/sphgeom/Region.h index 1283490f..f142b90e 100644 --- a/include/lsst/sphgeom/Region.h +++ b/include/lsst/sphgeom/Region.h @@ -149,8 +149,8 @@ class Region { static std::unique_ptr decode(std::uint8_t const * buffer, size_t n); ///@} - /// `flatten` flattens a nested Region to a list of regions. - static void flatten(Region const ®ion, std::vector> &result); + /// `getRegions` returns a vector of Region. + static std::vector> getRegions(Region const ®ion); ///@} }; diff --git a/python/lsst/sphgeom/_region.cc b/python/lsst/sphgeom/_region.cc index 783330cd..11cb15b1 100644 --- a/python/lsst/sphgeom/_region.cc +++ b/python/lsst/sphgeom/_region.cc @@ -70,11 +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("flatten", [](Region const ®ion) { - std::vector> result; - Region::flatten(region, result); - return result; - }, "region"_a); + cls.def_static("getRegions", Region::getRegions, "region"_a); } } // sphgeom diff --git a/src/PixelFinder.h b/src/PixelFinder.h index 92139902..718d184e 100644 --- a/src/PixelFinder.h +++ b/src/PixelFinder.h @@ -177,7 +177,7 @@ RangeSet findPixels(Region const & r, size_t maxRanges, int level) { auto rs1 = findPixels(region1, maxRanges, level); auto rs2 = findPixels(region2, maxRanges, level); s = rs1.join(rs2); - }else if (auto intersection_region = dynamic_cast(&r)) { + } 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); diff --git a/src/Region.cc b/src/Region.cc index b8e32107..73611492 100644 --- a/src/Region.cc +++ b/src/Region.cc @@ -72,20 +72,22 @@ std::unique_ptr Region::decode(std::uint8_t const * buffer, size_t n) { } throw std::runtime_error("Byte-string is not an encoded Region"); } - void Region::flatten(Region const ®ion, std::vector> &result) { - if (auto union_region = dynamic_cast(®ion)) { - for(int i = 0; i < 2; ++i) { - union_region->getOperand(i); - flatten(union_region->getOperand(i), result); - } - } else if(auto intersection_region = dynamic_cast(®ion)) { - for(int i = 0; i < 2; ++i) { - intersection_region->getOperand(i); - flatten(intersection_region->getOperand(i), result); - } - } else { - result.emplace_back(region.clone()); + +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 3e3418b1..eb62dca6 100644 --- a/tests/test_CompoundRegion.py +++ b/tests/test_CompoundRegion.py @@ -201,16 +201,22 @@ def testRelate(self): self.assertEqual(self.box.relate(self.instance), CONTAINS) self.assertEqual(self.faraway.relate(self.instance), DISJOINT) - def testFlatten(self): + 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) - c = UnionRegion(u1, u2) - d = Region.flatten(c) - self.assertEqual(d, [c1, b1, 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__": diff --git a/tests/test_HtmPixelization.py b/tests/test_HtmPixelization.py index 6f282eba..f7d17ae6 100644 --- a/tests/test_HtmPixelization.py +++ b/tests/test_HtmPixelization.py @@ -104,7 +104,7 @@ def test_envelope_and_interior(self): rsu2 = pixelization.envelope(union2) self.assertEqual(rsu2, rsu2 | rsu) - # CHeck with intersection + # 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)