python - 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 倍左右