Use a more robust intersection area computation

CHANGELOG.rst

@@ -4,17 +4,24 @@
 Release Notes
 =============
 
-.. 0.8.8 (unreleased)
+.. 0.8.9 (unreleased)
    ==================
 
-0.8.8 (unreleased)
-==================
+0.8.8 (17-July-2024)
+====================
+
+- Use a more robust algorithm for computing intersection polygon area that
+  ignores intersection polygons that raise ``MalformedPolygonError`` in the
+  ``spherical_geometry`` package. This may result in sub-optimal alignment
+  *order* but in practice, it should have minimal effect
+  on the end result. [#281]
 
 - ``align_wcs`` now will raise a custom exception of type ``NotEnoughCatalogs``
   when there are not enough input catalogs to perform alignment. [#203]
-
-- ``XYXYMatch`` now will raise a custom exception of type ``MatchSourceConfusionError``
-  when multipe reference sources match a single input source. [#204]
+
+- ``XYXYMatch`` now will raise a custom exception of type
+  ``MatchSourceConfusionError`` when multipe reference sources match a single
+  input source. [#204]
 
 
 0.8.7 (29-March-2024)

tweakwcs/imalign.py

@@ -777,11 +777,20 @@
     """
     nimg = len(images)
     m = np.zeros((nimg, nimg), dtype=np.double)
+    n_malformed = 0
     for i in range(nimg):
         for j in range(i + 1, nimg):
-            area = images[i]._guarded_intersection_area(images[j])
+            area, nf = images[i]._guarded_intersection_area(images[j])
+            n_malformed += nf
             m[j, i] = area
             m[i, j] = area
+
+    if n_malformed:
+        log.warning(
+            "MalformedPolygonError in spherical_geometry. Computed overlap "
+            "area is not accurate. Alignment order may be sub-optimal."
+        )
+
     return m
 
 
@@ -947,7 +956,7 @@
         # Also, when ref. catalog is static - revert to old tweakreg behavior
         im1 = images.pop(0)  # reference image
         im2 = images.pop(0)
-        overlap_area = im1._guarded_intersection_area(im2)
+        overlap_area = im1._guarded_intersection_area(im2)[0]
         return im1, im2, overlap_area
 
     m = overlap_matrix(images)
@@ -1018,7 +1027,19 @@
         # revert to old tweakreg behavior
         return images.pop(0), None
 
-    overlap_area = [refimage.intersection_area(im) for im in images]
+    n_malformed = 0
+    overlap_area = []
+    for im in images:
+        area, nf = refimage._guarded_intersection_area(im)
+        overlap_area.append(area)
+        n_malformed += nf
+
+    if n_malformed:
+        log.warning(
+            "MalformedPolygonError in spherical_geometry. Computed overlap "
+            "area is not accurate. Alignment order may be sub-optimal."
+        )
+
     idx = np.argmax(overlap_area)
 
     return images.pop(idx), overlap_area[idx]

tweakwcs/tests/test_wcsimage.py

@@ -120,7 +120,7 @@ def test_wcsimcat_intersections(mock_fits_wcs, rect_imcat):
 
 def test_wcsimcat_guarded_intersection_area(mock_fits_wcs, rect_imcat):
     assert np.allclose(
-        rect_imcat._guarded_intersection_area(rect_imcat),
+        rect_imcat._guarded_intersection_area(rect_imcat)[0],
         2.9904967391303217e-12, atol=0.0, rtol=5.0e-4
     )
 
@@ -241,7 +241,7 @@ def test_wcsgroupcat_intersections(mock_fits_wcs, rect_imcat):
 def test_wcsgroupcat_guarded_intersection_area(mock_fits_wcs, rect_imcat):
     g = WCSGroupCatalog(rect_imcat)
     assert np.allclose(
-        g._guarded_intersection_area(g), 2.9904967391303217e-12,
+        g._guarded_intersection_area(g)[0], 2.9904967391303217e-12,
         atol=0.0, rtol=5.0e-4
     )
 
@@ -458,7 +458,9 @@ def test_wcsrefcat_guarded_intersection_area(mock_fits_wcs, rect_imcat):
     ref.calc_tanp_xy(rect_imcat.corrector)
 
     assert np.allclose(
-        ref._guarded_intersection_area(ref), 2.9902125220360176e-10, atol=0.0,
+        ref._guarded_intersection_area(ref)[0],
+        2.9902125220360176e-10,
+        atol=0.0,
         rtol=0.005,
     )
 

tweakwcs/wcsimage.py

@@ -350,24 +350,32 @@
         intersections fail due to a bug/limitation of ``spherical_geometry``
         then the area of the valid intersections will be returned.
         If images do not intersect or intersection fails, 0 will be returned.
+
+        In addition to the area, this function returns the number of times
+        a failure in the ``spherical_geometry`` package occurred resulting in a
+        ``MalformedPolygonError`` error.
+
+        Return value: ``(area, nfailures)``
         """
         if isinstance(wcsim, (WCSImageCatalog, RefCatalog)):
             try:
-                return np.fabs(self.intersection(wcsim).area())
+                return np.fabs(self.intersection(wcsim).area()), 0
             except MalformedPolygonError:
-                return 0.0
+                return 0.0, 1
 
         else:
             # this is bug workaround:
             area = 0.0
+            nfailures = 0
             for wim in wcsim:
                 try:
                     area += np.fabs(
                         self.polygon.intersection(wim.polygon).area()
                     )
                 except MalformedPolygonError:
+                    nfailures += 1
                     continue
-            return area
+            return area, nfailures
 
     def _calc_chip_bounding_polygon(self, stepsize=None):
         """
@@ -656,8 +664,20 @@
         intersections fail due to a bug/limitation of ``spherical_geometry``
         then the area of the valid intersections will be returned.
         If images do not intersect or intersection fails, 0 will be returned.
+
+        In addition to the area, this function returns the number of times
+        a failure in the ``spherical_geometry`` package occurred resulting in a
+        ``MalformedPolygonError`` error.
+
+        Return value: ``(area, nfailures)``
         """
-        return sum(im._guarded_intersection_area(wcsim) for im in self._images)
+        area = 0
+        nfailures = 0
+        for im in self._images:
+            a, nf = im._guarded_intersection_area(wcsim)
+            area += a
+            nfailures += nf
+        return area, nfailures
 
     def _aproximate_bb(self):
         if not self._images:
@@ -1551,24 +1571,33 @@
         intersections fail due to a bug/limitation of ``spherical_geometry``
         then the area of the valid intersections will be returned.
         If images do not intersect or intersection fails, 0 will be returned.
+
+        In addition to the area, this function returns the number of times
+        a failure in the ``spherical_geometry`` package occurred resulting in a
+        ``MalformedPolygonError`` error.
+
+        Return value: ``(area, nfailures)``
         """
         if isinstance(wcsim, (WCSImageCatalog, RefCatalog)):
             try:
-                return np.fabs(self.intersection(wcsim).area())
+                return np.fabs(self.intersection(wcsim).area()), 0
             except MalformedPolygonError:
-                return 0.0
+                return 0.0, 1
 
         else:
             # this is bug workaround:
             area = 0.0
+            nfailures = 0
             for wim in wcsim:
                 try:
                     area += np.fabs(
                         self.polygon.intersection(wim.polygon).area()
                     )
                 except MalformedPolygonError:
+                    nfailures += 1
                     continue
-            return area
+
+            return area, nfailures
 
     def _calc_cat_convex_hull(self):
         """