Overlapping polygons in Python PIL

Question

Instead of overwriting the overlapping regions of multiple polygons by the value of the last polygon drawn, I would like to draw the mean value of these polygons. Is this possible in Python PIL?

The overlapping pixels in the example should have the value of 1.5.

In the full working program I have to draw about 100000 polygons (that may or may not intersect) on a very large grid which is the reason I use PIL instead of Numpy.

from PIL import Image, ImageDraw
import numpy as np
import matplotlib.pyplot as plt

img = Image.new('F', (50, 50), 0)

ImageDraw.Draw(img).polygon([(20, 20), (20, 40), (40, 30), (30, 20)],
                            fill=1., outline=None)
ImageDraw.Draw(img).polygon([(10, 5), (10, 25), (25, 25), (25, 10)],
                            fill=2., outline=None)

myimg = np.ma.masked_equal(np.array(img), 0.)
plt.imshow(myimg, interpolation="None")
plt.colorbar()
plt.show()

Answer 1

I suggest scikit-image, skimage.draw.polygon() returns coordinates in the polygon. Here is an example. Create some random polygon data first:

import pylab as pl
from random import randint
import numpy as np
from skimage import draw

W, H = 800, 600

def make_poly(x0, y0, r, n):
    a = np.linspace(0, np.pi*2, n, endpoint=False)
    x = x0 + r * np.cos(a)
    y = y0 + r * np.sin(a)
    return y, x

count_buf = np.zeros((H, W))
sum_buf = np.zeros((H, W))

N = 2000

polys = []
for i in range(N):
    x0, y0, r, n = randint(10, W-10), randint(10, H-10), randint(10, 50), randint(3, 10)
    polys.append((make_poly(x0, y0, r, n), randint(1, 10)))

Then draw the polygons:

for (y, x), v in polys:
    rr, cc = draw.polygon(y, x, (H, W))
    count_buf[rr, cc] += 1
    sum_buf[rr, cc] += v

mean_buf = np.zeros_like(sum_buf)
mask = count_buf > 0
mean_buf[mask] = sum_buf[mask] / count_buf[mask]

the time is about 1.5s on my pc to draw 2000 polygons with average radius 30 px.

Here is the result:

If you need better speed, you can copy the following code in scikit-image:

https://github.com/scikit-image/scikit-image/blob/master/skimage/draw/_draw.pyx#L189

https://github.com/scikit-image/scikit-image/blob/master/skimage/_shared/geometry.pyx#L7

and change the count_buf and sum_buf in the for loop if point_in_polygon() returns True.

Edit

Here is the Cython code:

%%cython
#cython: cdivision=True
#cython: boundscheck=False
#cython: wraparound=False

import numpy as np
cimport numpy as np
from libc.math cimport ceil

cdef unsigned char point_in_polygon(double[::1] xp, double[::1] yp,
                                           double x, double y):
    cdef Py_ssize_t i
    cdef unsigned char c = 0
    cdef Py_ssize_t j = xp.shape[0] - 1
    for i in range(xp.shape[0]):
        if (
            (((yp[i] <= y) and (y < yp[j])) or
            ((yp[j] <= y) and (y < yp[i])))
            and (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])
        ):
            c = not c
        j = i
    return c


cdef class PolygonAccumulator:

    cdef int width, height
    cdef int[:, ::1] count_buf
    cdef double[:, ::1] sum_buf

    def __cinit__(self, width, height):
        self.width = width
        self.height = height
        shape = (height, width)
        self.count_buf = np.zeros(shape, dtype=int)
        self.sum_buf = np.zeros(shape, dtype=float)

    def reset(self):
        self.count_buf[:, :] = 0
        self.sum_buf[:, :] = 0

    def add_polygon(self, ya, xa, double value):
        cdef Py_ssize_t minr = int(max(0, np.min(ya)))
        cdef Py_ssize_t maxr = int(ceil(np.max(ya)))
        cdef Py_ssize_t minc = int(max(0, np.min(xa)))
        cdef Py_ssize_t maxc = int(ceil(np.max(xa)))

        cdef double[::1] x = xa
        cdef double[::1] y = ya

        cdef Py_ssize_t r, c

        maxr = min(self.height - 1, maxr)
        maxc = min(self.width  - 1, maxc)

        for r in range(minr, maxr+1):
            for c in range(minc, maxc+1):
                if point_in_polygon(x, y, c, r):
                    self.count_buf[r, c] += 1
                    self.sum_buf[r, c] += value

    def mean(self):
        count_buf = self.count_buf.base
        sum_buf = self.sum_buf.base
        mean_buf = np.zeros_like(sum_buf)
        mask = count_buf > 0
        mean_buf[mask] = sum_buf[mask] / count_buf[mask]
        return mean_buf

To draw the polygons:

pa = PolygonAccumulator(800, 600)
for (y, x), value in polys:
    pa.add_polygon(y, x, value)
pl.imshow(pa.mean(), cmap="gray")

It's about 4.5x faster than skimage.draw.polygon()