Fix xy coordinates of SegmentationImage polygon vertices

CHANGES.rst

@@ -15,7 +15,7 @@ New Features
     much faster implementation of binary dilation. [#1638]
 
   - Added a ``scale`` keyword to the ``SegmentationImage.to_patches()``
-    method to scale the sizes of the polygon patches. [#1641]
+    method to scale the sizes of the polygon patches. [#1641, #1646]
 
 Bug Fixes
 ^^^^^^^^^
@@ -31,6 +31,9 @@ Bug Fixes
     raise an error if the ``contrast`` keyword was set to 1 (meaning no
     deblending). [#1636]
 
+  - Fixed an issue where the vertices of the ``SegmentationImage``
+    ``polygons`` were shifted by 0.5 pixels in both x and y. [#1646]
+
 API Changes
 ^^^^^^^^^^^
 

photutils/segmentation/core.py

@@ -1176,11 +1176,15 @@ def _geo_polygons(self):
         Each item in the list is tuple of (polygon, value) where the
         polygon is a GeoJSON-like dict and the value is the label from
         the segmentation image.
+
+        Note that the coordinates of these polygon vertices are in a
+        reference frame with the (0, 0) origin at the *lower-left*
+        corner of the lower-left pixel.
         """
         from rasterio.features import shapes
 
         polygons = list(shapes(self.data.astype('int32'), connectivity=8))
-        polygons.sort(key=lambda x: x[1])
+        polygons.sort(key=lambda x: x[1])  # sort in label order
 
         # do not include polygons for background (label = 0)
         return polygons[1:]
@@ -1191,11 +1195,16 @@ def polygons(self):
         A list of `Shapely <https://shapely.readthedocs.io/en/stable/>`_
         polygons representing each source segment.
         """
+        from shapely import transform
         from shapely.geometry import shape
 
         polygons = []
         for geo_poly in self._geo_polygons:
             polygons.append(shape(geo_poly[0]))
+        # shift the vertices so that the (0, 0) origin is at the
+        # center of the lower-left pixel
+        polygons = transform(polygons, lambda x: x - [0.5, 0.5])
+
         return polygons
 
     def to_patches(self, *, origin=(0, 0), scale=1.0, **kwargs):
@@ -1226,20 +1235,11 @@ def to_patches(self, *, origin=(0, 0), scale=1.0, **kwargs):
         patch_kwargs = {'edgecolor': 'white', 'facecolor': 'none'}
         patch_kwargs.update(kwargs)
 
-        # This is the shapely equivalent for patches instead of using
-        # self._geo_polygons below.
-        # patches = []
-        # for poly in self.polygons:
-        #     x = np.array(poly.exterior.coords.xy[0])
-        #     y = np.array(poly.exterior.coords.xy[1])
-        #     xy = np.column_stack((x, y)) - origin - np.array((0.5, 0.5))
-        #     patches.append(Polygon(xy, **patch_kwargs))
-
         patches = []
-        for geo_poly in self._geo_polygons:
-            xy = (np.array(geo_poly[0]['coordinates'][0]) - origin
-                  - np.array((0.5, 0.5)))
-            xy *= scale
+        for poly in self.polygons:
+            xy = np.array(poly.exterior.coords)
+            xy = scale * (xy + 0.5) - 0.5
+            xy -= origin
             patches.append(Polygon(xy, **patch_kwargs))
 
         return patches

photutils/segmentation/tests/test_core.py

@@ -404,6 +404,16 @@ def test_polygons(self):
         assert len(polygons) == self.segm.nlabels
         assert isinstance(polygons[0], Polygon)
 
+        data = np.zeros((5, 5), dtype=int)
+        data[2, 2] = 10
+        segm = SegmentationImage(data)
+        polygons = segm.polygons
+        assert len(polygons) == 1
+        verts = np.array(polygons[0].exterior.coords)
+        expected_verts = np.array([[1.5, 1.5], [1.5, 2.5], [2.5, 2.5],
+                                   [2.5, 1.5], [1.5, 1.5]])
+        assert_equal(verts, expected_verts)
+
     @pytest.mark.skipif(not HAS_RASTERIO, reason='rasterio is required')
     @pytest.mark.skipif(not HAS_MATPLOTLIB, reason='matplotlib is required')
     def test_patches(self):
@@ -417,7 +427,8 @@ def test_patches(self):
         patches2 = self.segm.to_patches(scale=scale)
         v1 = patches[0].get_verts()
         v2 = patches2[0].get_verts()
-        assert_allclose(v2, v1 * scale)
+        v3 = scale * (v1 + 0.5) - 0.5
+        assert_allclose(v2, v3)
 
         patches = self.segm.plot_patches(edgecolor='red')
         assert isinstance(patches[0], Polygon)