Merge pull request #205 from mcara/fix-crash-malformed-poly

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 @@ def overlap_matrix(images):
  """
  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 @@ def _max_overlap_pair(images, enforce_user_order):
  # 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 @@ def _max_overlap_image(refimage, images, enforce_user_order):
  # 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 @@ def _guarded_intersection_area(self, wcsim):
  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 @@ def _guarded_intersection_area(self, wcsim):
  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 @@ def _guarded_intersection_area(self, wcsim):
  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):
  """