Python PIL中的重叠多边形

Question

我不想用绘制的最后一个多边形的值覆盖多个多边形的重叠区域，而是绘制这些多边形的平均值。 这在 Python PIL 中可行吗？

示例中的重叠像素的值应为 1.5。

在完整的工作程序中，我必须在一个非常大的网格上绘制大约 100000 个多边形（可能相交或不相交），这就是我使用 PIL 而不是 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

我建议scikit-image、skimage.draw.polygon() 返回多边形中的坐标。这是一个例子。先创建一些随机多边形数据：

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)))

然后绘制多边形：

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]

在我的电脑上绘制 2000 个平均半径为 30 像素的多边形大约需要 1.5 秒。

结果如下：

如果你需要更好的速度，可以在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

如果 point_in_polygon() 返回 True，则在 for 循环中更改 count_buf 和 sum_buf。

编辑

这是 Cython 代码：

%%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

绘制多边形：

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

比 skimage.draw.polygon() 快大约 4.5 倍