How to convert NumPy arrays obtained from cv2.findContours to Shapely polygons?

Question

I am using CV2 to find contours from an image and then converting them into polygons using Shapely. I am currently stuck because when I try putting one of the contour arrays into Polygon() from Shapely it throws an unspecified error.

I have double-checked that I imported everything I needed, and that creating a Shapely polygon works when I manually enter the array coordinate points.

Here is the problematic section of the code:

contours, hierarchy = cv2.findContours(thresh, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
testcontour = contours[1]

ply = Polygon(testcontour)

Where the list of contours looks like this:

contours = [np.array([[[700, 700]],
                      [[700, 899]],
                      [[899, 899]],
                      [[899, 700]]]), 
            np.array([[[774, 775]], 
                      [[775, 774]],
                      [[824, 774]],
                      [[825, 775]],
                      [[825, 824]],
                      [[824, 825]],
                      [[775, 825]],
                      [[774, 824]]]), 
            np.array([[[200, 200]],
                      [[200, 399]],
                      [[399, 399]],
                      [[399, 200]]]), 
            np.array([[[274, 275]],
                      [[275, 274]],
                      [[324, 274]],
                      [[325, 275]],
                      [[325, 324]],
                      [[324, 325]],
                      [[275, 325]],
                      [[274, 324]]])]

The error I get is:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-65-4124f49b42e1> in <module>
----> 1 ply = Polygon(testcontour)

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\geometry\polygon.py in __init__(self, shell, holes)
    238 
    239         if shell is not None:
--> 240             ret = geos_polygon_from_py(shell, holes)
    241             if ret is not None:
    242                 self._geom, self._ndim = ret

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\geometry\polygon.py in geos_polygon_from_py(shell, holes)
    492 
    493     if shell is not None:
--> 494         ret = geos_linearring_from_py(shell)
    495         if ret is None:
    496             return None

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

AssertionError:

score 21 · Accepted Answer · answered Sep 17 '19 at 10:41

The problem is that for some reason cv2.findContours returns each contour as a 3D NumPy array with one redundant dimension:

>>> contours[1]
array([[[774, 775]],
       [[775, 774]],
       [[824, 774]],
       [[825, 775]],
       [[825, 824]],
       [[824, 825]],
       [[775, 825]],
       [[774, 824]]])

but Shapely expects a 2D array in this form (see the docs):

array([[774, 775],
       [775, 774],
       [824, 774],
       [825, 775],
       [825, 824],
       [824, 825],
       [775, 825],
       [774, 824]])

So, what we can do is to use np.squeeze to remove that redundant dimension, and use the result to obtain our polygon:

import numpy as np
from shapely.geometry import Polygon

contour = np.squeeze(contours[1])
polygon = Polygon(contour)
print(polygon.wkt)
# POLYGON ((774 775, 775 774, 824 774, 825 775, 825 824, 824 825, 775 825, 774 824, 774 775))

In case if you want to convert all of the contours at once, I would do it like this:

contours = map(np.squeeze, contours)  # removing redundant dimensions
polygons = map(Polygon, contours)  # converting to Polygons
multipolygon = MultiPolygon(polygons)  # putting it all together in a MultiPolygon

The resulting multipolygon will look like this:

And to get the second polygon from here you would just write:

my_polygon = multipolygon[1]
print(my_polygon.wkt)
# POLYGON ((774 775, 775 774, 824 774, 825 775, 825 824, 824 825, 775 825, 774 824, 774 775))

score 4 · Answer 2 · answered Feb 20 '23 at 14:21

Below is code which takes into account also hierarchy of contours. It recursively walks through the hierarchy structure and adds/subtracts child contours in alternating manner.

The code also handles special cases when OpenCV returns trivial contours (single point) or polygons which are not shapely-valid.

import numpy as np
import cv2
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt

# Read and threshold image
img = cv2.imread('/path/to/image')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
_,img_th = cv2.threshold(gray,185,255,cv2.THRESH_BINARY)

# Find contours and hierarchy
contours, hierarchy = cv2.findContours(img_th, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
hierarchy = np.squeeze(hierarchy)
print("Number of contours found = " + str(len(contours)))

# Plot image with contours
cv2.drawContours(img, contours, -1, (255, 0, 0), 3)
plt.imshow(img,origin='lower')
plt.show()

# Define function to merge individual contours into one MultiPolygon
def merge_polygons(polygon:MultiPolygon,idx:int,add:bool) -> MultiPolygon:
    """
    polygon: Main polygon to which a new polygon is added
    idx: Index of contour
    add: If this contour should be added (True) or subtracted (False)
    """

    # Get contour from global list of contours
    contour = np.squeeze(contours[idx])

    # cv2.findContours() sometimes returns a single point -> skip this case
    if len(contour) > 2:
        # Convert contour to shapely polygon
        new_poly = Polygon(contour)

        # Not all polygons are shapely-valid (self intersection, etc.)
        if not new_poly.is_valid:
            # Convert invalid polygon to valid
            new_poly = new_poly.buffer(0)

        # Merge new polygon with the main one
        if add: polygon = polygon.union(new_poly)
        else:   polygon = polygon.difference(new_poly)

    # Check if current polygon has a child
    child_idx = hierarchy[idx][2]
    if child_idx >= 0:
        # Call this function recursively, negate `add` parameter
        polygon = merge_polygons(polygon,child_idx,not add)

    # Check if there is some next polygon at the same hierarchy level
    next_idx = hierarchy[idx][0]
    if next_idx >= 0:
        # Call this function recursively
        polygon = merge_polygons(polygon,next_idx,add)

    return polygon

# Call the function with an initial empty polygon and start from contour 0
polygon = merge_polygons(MultiPolygon(),0,True)

Results: Sample image with contours and final shapely multipolygon.

How to convert NumPy arrays obtained from cv2.findContours to Shapely polygons?

2 Answers2

Linked