Implement flattenting of CompoundRegion.

Costructor of CompoundRegions will now flatten nested CompoundRegions of the
same type.

include/lsst/sphgeom/CompoundRegion.h

@@ -106,6 +106,10 @@ class CompoundRegion : public Region {
     // Check if we have any operands.
     bool empty() const { return _operands.empty(); }
 
+    // Flatten vector of regions in-place.
+    template <typename Compound>
+    void flatten_operands();
+
 private:
     std::vector<std::unique_ptr<Region>> _operands;
 };
@@ -118,7 +122,8 @@ class UnionRegion : public CompoundRegion {
 public:
     static constexpr std::uint8_t TYPE_CODE = 'u';
 
-    using CompoundRegion::CompoundRegion;
+    /// Construct by taking ownership of operands.
+    explicit UnionRegion(std::vector<std::unique_ptr<Region>> operands);
 
     // Region interface.
     std::unique_ptr<Region> clone() const override { return std::make_unique<UnionRegion>(*this); }
@@ -157,7 +162,8 @@ class IntersectionRegion : public CompoundRegion {
 public:
     static constexpr std::uint8_t TYPE_CODE = 'i';
 
-    using CompoundRegion::CompoundRegion;
+    /// Construct by taking ownership of operands.
+    explicit IntersectionRegion(std::vector<std::unique_ptr<Region>> operands);
 
     // Region interface.
     std::unique_ptr<Region> clone() const override { return std::make_unique<IntersectionRegion>(*this); }

src/CompoundRegion.cc

@@ -32,6 +32,7 @@
 
 #include "lsst/sphgeom/CompoundRegion.h"
 
+#include <algorithm>
 #include <iostream>
 #include <stdexcept>
 
@@ -91,6 +92,24 @@ CompoundRegion::CompoundRegion(CompoundRegion const &other)
     }
 }
 
+// Flatten vector of regions in-place.
+template <typename Compound>
+void CompoundRegion::flatten_operands() {
+    for (size_t i = 0; i != _operands.size(); ) {
+        if (auto compound = dynamic_cast<Compound*>(_operands[i].get())) {
+            // Move all regions from this operand, then remove it.
+            std::move(
+                compound->_operands.begin(),
+                compound->_operands.end(),
+                std::inserter(_operands, _operands.begin() + i + 1)
+            );
+            _operands.erase(_operands.begin() + i);
+        } else {
+            ++ i;
+        }
+    }
+}
+
 Relationship CompoundRegion::relate(Box const &b) const { return relate(static_cast<Region const &>(b)); }
 Relationship CompoundRegion::relate(Circle const &c) const { return relate(static_cast<Region const &>(c)); }
 Relationship CompoundRegion::relate(ConvexPolygon const &p) const { return relate(static_cast<Region const &>(p)); }
@@ -143,6 +162,12 @@ std::unique_ptr<CompoundRegion> CompoundRegion::decode(std::uint8_t const *buffe
     }
 }
 
+UnionRegion::UnionRegion(std::vector<std::unique_ptr<Region>> operands)
+    : CompoundRegion(std::move(operands))
+{
+    flatten_operands<UnionRegion>();
+}
+
 bool UnionRegion::isEmpty() const {
     // It can be empty when there are no operands or all operands are empty.
     for (auto&& operand: operands()) {
@@ -241,6 +266,12 @@ TriState UnionRegion::overlaps(Ellipse const &e) const {
     return overlaps(static_cast<Region const&>(e));
 }
 
+IntersectionRegion::IntersectionRegion(std::vector<std::unique_ptr<Region>> operands)
+    : CompoundRegion(std::move(operands))
+{
+    flatten_operands<IntersectionRegion>();
+}
+
 bool IntersectionRegion::isEmpty() const {
     // Intersection is harder to decide - the only clear case is when there are
     // no operands, which we declare to be equivalent to full sphere. Other

tests/test_CompoundRegion.py

@@ -261,10 +261,16 @@ def testGetRegion(self):
         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])
+        self.assertEqual(Region.getRegions(i1), [c1, b1])
+        self.assertEqual(Region.getRegions(u1), [c1, b1])
+        # Compounds of compounds will be flattened, order preserved.
+        self.assertEqual(Region.getRegions(ir), [c1, b1, c2, b2])
+        self.assertEqual(Region.getRegions(ur), [c1, b1, c2, b2])
+
+        # TODO: This test fails because CompoundRegion does not define
+        # equality operator, and it is non-trivial to add one.
+        # ur2 = UnionRegion(u1, i1, u2)
+        # self.assertEqual(Region.getRegions(ur2), [c1, b1, i1, c2, b2])
 
 
 if __name__ == "__main__":