Extract image segment from polygon, with skimage

I would like to get the sub-image that results from cutting out a polygon within an image.

I have an image in skimage, and I have a polygon in matplotlib.patches.

How to do it?

Below is what I have tried. I am not necessarily looking for a way similar to below, I am looking for the cleanest, most efficient implementation.

With this code, the polygon correctly overlays the part of the image I want to extract (but doesn't extract the segment of interest):

import numpy as np
import skimage.io as io
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon

I = io.imread(fp) # fp is path to image
plt.imshow(I)
ax = plt.gca()

polygons, color = [], []   
c = np.random.random((1, 3)).tolist()[0]
for seg in ann['segmentation']:
    poly = np.array(seg).reshape((len(seg)/2, 2))
    polygons.append(Polygon(poly, True,alpha=0.4))
    color.append(c)
p = PatchCollection(polygons, facecolors=color, edgecolors=(0,0,0,1), linewidths=3, alpha=0.4)

ax.add_collection(p)

But when I try to obtain the segmented image with this code, the overlay displays wrongly:

fig, ax = plt.subplots()
im = ax.imshow(I)
im.set_clip_path(polygon)
plt.axis('off')
plt.show()

It looks like the Y coordinates of the polygon just need to flipped (especially since image above shows Y axis ordered other way around), but that's not the case:

a = polygons[0].xy.copy()
a[:,1] = im._A.shape[0] - a[:,1]
newPoly = Polygon(a, True,alpha=0.4)
fig, ax = plt.subplots()
im = ax.imshow(I)
im.set_clip_path(newPoly)
plt.axis('off')
plt.show()

(In fact, there is not only an offset issue in X coordinates, there is even a scale issue in Y coordinates. I have no idea why)

I can't explain the strange behavior also. In any case recently in another question I've suggested a recipe that might help here (although I wouldn't call it the cleanest solution). With this (not very pretty) code:

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from matplotlib import path

class LineBuilder:
    def __init__(self, line,ax,color):
        self.line = line
        self.ax = ax
        self.color = color
        self.xs = []
        self.ys = []
        self.cid = line.figure.canvas.mpl_connect('button_press_event', self)
        self.counter = 0
        self.shape_counter = 0
        self.shape = {}
        self.precision = 10

    def __call__(self, event):
        if event.inaxes!=self.line.axes: return
        if self.counter == 0:
            self.xs.append(event.xdata)
            self.ys.append(event.ydata)
        if np.abs(event.xdata-self.xs[0])<=self.precision and np.abs(event.ydata-self.ys[0])<=self.precision and self.counter != 0:
            self.xs.append(self.xs[0])
            self.ys.append(self.ys[0])
            self.ax.scatter(self.xs,self.ys,s=120,color=self.color)
            self.ax.scatter(self.xs[0],self.ys[0],s=80,color='blue')
            self.ax.plot(self.xs,self.ys,color=self.color)
            self.line.figure.canvas.draw()
            self.shape[self.shape_counter] = [self.xs,self.ys]
            self.shape_counter = self.shape_counter + 1
            self.xs = []
            self.ys = []
            self.counter = 0
        else:
            if self.counter != 0:
                self.xs.append(event.xdata)
                self.ys.append(event.ydata)
            self.ax.scatter(self.xs,self.ys,s=120,color=self.color)
            self.ax.plot(self.xs,self.ys,color=self.color)
            self.line.figure.canvas.draw()
            self.counter = self.counter + 1

def create_shape_on_image(data,cmap='jet'):
    def change_shapes(shapes):
        new_shapes = {}
        for i in range(len(shapes)):
            l = len(shapes[i][1])
            new_shapes[i] = np.zeros((l,2),dtype='int')
            for j in range(l):
                new_shapes[i][j,0] = shapes[i][0][j]
                new_shapes[i][j,1] = shapes[i][1][j]
        return new_shapes
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_title('click to include shape markers (10 pixel precision to close the shape)')
    line = ax.imshow(data) 
    ax.set_xlim(0,data[:,:,0].shape[1])
    ax.set_ylim(0,data[:,:,0].shape[0])
    linebuilder = LineBuilder(line,ax,'red')
    plt.gca().invert_yaxis()
    plt.show()
    new_shapes = change_shapes(linebuilder.shape)
    return new_shapes

img = mpimg.imread('wm4HA.png')

shapes = create_shape_on_image(img)[0]

xx,yy = np.meshgrid(range(img.shape[0]),range(img.shape[1]))
shapes = np.hstack((shapes[:,1][:,np.newaxis],shapes[:,0][:,np.newaxis]))
p = path.Path(shapes)
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        if not p.contains_point((i,j)):
            img[i,j,:] = np.array([0,0,0,0])

plt.imshow(img)

plt.show()

I can build your desired result:

The code that is most important to you is this:

img = mpimg.imread('wm4HA.png')

shapes = create_shape_on_image(img)[0] # Here I'm calling a function to build a polygon.

xx,yy = np.meshgrid(range(img.shape[0]),range(img.shape[1]))
shapes = np.hstack((shapes[:,1][:,np.newaxis],shapes[:,0][:,np.newaxis]))
p = path.Path(shapes)
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        if not p.contains_point((i,j)):
            img[i,j,:] = np.array([0,0,0,0])

plt.imshow(img)

plt.show()

In this case I used a recipe to build a polygon by point and click: