Search code examples
pythonscikit-image

Correctly closing of polygons generated using skimage.measure.find_contours()

I'm currently using skimage.measure.find_contours() to find contours on a surface. Now that I've found the contours I need to able to find the area enclosed within them.

When all of the vertices are within the data set this is fine as a have a fully enclosed polygon. However, how do I ensure the polygon is fully enclosed if the contour breaches the edge of the surface, either at an edge or at a corner? When this happens I would like to use the edge of the surface as additional vertices to close off the polygon. For example in the following image, with contours shown, you can see that the contours end at the edge of the image, how do I close them up? Also in the example of the brown contour, which is just a single line, I don't think I want an area returned, how would I single out this case?

Contour example

I know I can check for enclosed contours/polygons by checking if the last vertices of the polygon is the same as the first.

I have code for calculating the area inside a polygon, taken from here

def find_area(array):
    a = 0
    ox,oy = array[0]
    for x,y in array[1:]:
        a += (x*oy-y*ox)
        ox,oy = x,y
    return -a/2

I just need help in closing off the polygons. And checking for the different cases that might occur.

Thanks

Update:

After applying the solution suggested by @soupault I have this code:

import numpy as np
import matplotlib.pyplot as plt

from skimage import measure


# Construct some test data
x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]
r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))

# Coordinates of point of interest
pt = [(49,75)]

# Apply thresholding to the surface
threshold = 0.8
blobs = r > threshold

# Make a labelled image based on the thresholding regions
blobs_labels = measure.label(blobs, background = 0)

# Show the thresholded regions
plt.figure()
plt.imshow(blobs_labels, cmap='spectral')

# Apply regionprops to charactersie each of the regions
props = measure.regionprops(blobs_labels, intensity_image = r)

# Loop through each region in regionprops, identify if the point of interest is
# in that region. If so, plot the region and print it's area.
plt.figure()
plt.imshow(r, cmap='Greys')
plt.plot(pt[0][0], pt[0][1],'rx')
for prop in props:
    coords = prop.coords
    if np.sum(np.all(coords[:,[1,0]] == pt[0], axis=1)):
        plt.plot(coords[:,1],coords[:,0],'r.')
        print(prop.area)

This solution assumes that each pixel is 1x1 in size. In my real data solution this isn't the case so I have also applied the following function to apply linear interpolation to the data. I believe you can also apply a similar function to make the area of each pixel smaller and increase the resolution of the data.

import numpy as np
from scipy import interpolate

def interpolate_patch(x,y,patch):
    x_interp = np.arange(np.ceil(x[0]), x[-1], 1)
    y_interp = np.arange(np.ceil(y[0]), y[-1], 1)

    f = interpolate.interp2d(x, y, patch, kind='linear')

    patch_interp = f(x_interp, y_interp)

    return x_interp, y_interp, patch_interp
Solution

  • If you need to measure the properties of different regions, it is natural to start with finding the regions (not contours).

    The algorithm will be the following, in this case:

    1. Prepare a labeled image:

      1.a Either fill the areas between different contour lines with the different colors;

      1.b Or apply some image thresholding function, and then run skimage.measure.label (http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.label);

    2. Execute regionprops using the very labeled image as an input (http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops);

    3. Iterate over regions in regionprops and calculate the desired parameters (area, perimeter, etc).

    Once you identified the regions in your image via regionprops, you can call .coords for each of them to get the enclosed contour.