Update polygon test and add explanation.

lsst · Nov 5, 2021 · 018c791 · 018c791
1 parent 7e580e4
commit 018c791
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 12 deletions.
diff --git a/src/geom/polygon/Polygon.cc b/src/geom/polygon/Polygon.cc
@@ -121,7 +121,12 @@ void addSubSampledEdge(std::vector<LsstPoint>& vertices,  // Vector of points to
 /// @internal Calculate area of overlap between polygon and pixel
 double pixelOverlap(BoostPolygon const& poly, int const x, int const y) {
     std::vector<BoostPolygon> overlap;  // Overlap between pixel and polygon
-    LsstBox const pixel(lsst::geom::Point2D(x - 0.5, y - 0.5), lsst::geom::Point2D(x + 0.5, y + 0.5));
+    LsstBox const pixel(lsst::geom::Point2D(static_cast<double>(x) - 0.5, static_cast<double>(y) - 0.5),
+                        lsst::geom::Point2D(static_cast<double>(x) + 0.5, static_cast<double>(y) + 0.5));
+    // Note that the output of boost::geometry::intersection depends
+    // on values down to the precision limit, so minor variations
+    // in poly input can lead to surprisingly large variations in the
+    // output overlap regions and area computation.
     boost::geometry::intersection(poly, pixel, overlap);
     double area = 0.0;
     for (auto const &i : overlap) {

diff --git a/tests/test_polygon.py b/tests/test_polygon.py
@@ -255,18 +255,27 @@ def testImage(self):
         """Test Polygon.createImage"""
         if display:
             import lsst.afw.display as afwDisplay
+
         for i, num in enumerate(range(3, 30)):
-            poly = self.polygon(num, 25, 75, 75)
-            box = lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
-                                  lsst.geom.Extent2I(115, 115))
-            image = poly.createImage(box)
-            if display:
-                disp = afwDisplay.Display(frame=i + 1)
-                disp.mtv(image, title=f"Polygon nside={num}")
-                for p1, p2 in poly.getEdges():
-                    disp.line((p1, p2))
-            self.assertAlmostEqual(
-                image.getArray().sum()/poly.calculateArea(), 1.0, 6)
+            # We need to test at different centers, because the
+            # boost::intersection algorithm depends sensitively on
+            # input floating-point values, and you will get different
+            # aliasing depending on the central pixel value when
+            # generating the polygon from numpy values.
+            for cent in [-75, -50, -25, 0, 25, 50, 75]:
+                poly = self.polygon(num, 25, cent, cent)
+                box = lsst.geom.Box2I(lsst.geom.Point2I(cent - 60, cent - 60),
+                                      lsst.geom.Extent2I(115, 115))
+                image = poly.createImage(box)
+                if display:
+                    disp = afwDisplay.Display(frame=i + 1)
+                    disp.mtv(image, title=f"Polygon nside={num}")
+                    for p1, p2 in poly.getEdges():
+                        disp.line((p1, p2))
+                # Some computations of the image area have such large aliasing
+                # in boost::intersection that the precision required here is 0.025.
+                self.assertFloatsAlmostEqual(
+                    image.getArray().sum(), poly.calculateArea(), rtol=0.025)
 
     def testTransform(self):
         """Test constructor for Polygon involving transforms"""