diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ea8ebdf..f57e9d4 100644 --- a/CHANGELOG.rst +++ b/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) diff --git a/tweakwcs/imalign.py b/tweakwcs/imalign.py index 51994de..6394168 100644 --- a/tweakwcs/imalign.py +++ b/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] diff --git a/tweakwcs/tests/test_wcsimage.py b/tweakwcs/tests/test_wcsimage.py index bc1af8e..420fd62 100644 --- a/tweakwcs/tests/test_wcsimage.py +++ b/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, ) diff --git a/tweakwcs/wcsimage.py b/tweakwcs/wcsimage.py index 691f1bf..cde5a63 100644 --- a/tweakwcs/wcsimage.py +++ b/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): """