计算许多三角形的有符号面积

Question

我需要计算 2D 中许多三角形的 signed area，由形状为 (2, 3, n) 的 numpy 数组给出（x/y 坐标，三角形中的节点，三角形的数量）。我正在寻找一种快速完成的方法，到目前为止我能想到的最好的方法是

import numpy
import perfplot


def six(p):
    return (
        +p[0][2] * p[1][0]
        + p[0][0] * p[1][1]
        + p[0][1] * p[1][2]
        - p[0][2] * p[1][1]
        - p[0][0] * p[1][2]
        - p[0][1] * p[1][0]
    ) / 2


def mix(p):
    return (
        +p[0][2] * (p[1][0] - p[1][1])
        + p[0][0] * (p[1][1] - p[1][2])
        + p[0][1] * (p[1][2] - p[1][0])
    ) / 2


def mix2(p):
    p1 = p[1] - p[1][[1, 2, 0]]
    return (+p[0][2] * p1[0] + p[0][0] * p1[1] + p[0][1] * p1[2]) / 2


def cross(p):
    e1 = p[:, 1] - p[:, 0]
    e2 = p[:, 2] - p[:, 0]
    return (e1[0] * e2[1] - e1[1] * e2[0]) / 2


def einsum(p):
    return (
        +numpy.einsum("ij,ij->j", p[0][[2, 0, 1]], p[1][[0, 1, 2]])
        - numpy.einsum("ij,ij->j", p[0][[2, 0, 1]], p[1][[1, 2, 0]])
    ) / 2


def einsum2(p):
    return numpy.einsum("ij,ij->j", p[0][[2, 0, 1]], p[1] - p[1][[1, 2, 0]]) / 2


def einsum3(p):
    return (
        numpy.einsum(
            "ij,ij->j", numpy.roll(p[0], 1, axis=0), p[1] - numpy.roll(p[1], 2, axis=0)
        )
        / 2
    )


perfplot.save(
    "out.png",
    setup=lambda n: numpy.random.rand(2, 3, n),
    kernels=[six, mix, mix2, cross, einsum, einsum2, einsum3],
    n_range=[2 ** k for k in range(19)],
)

关于如何提高效率的任何提示？

Answer 1

我认为问题不在于计算，而在于数据在内存中的组织方式。如果您不需要保留原始数据，您可以尝试使用inplace 函数尽量减少内存中临时对象的数量：

def mymix(p):

    p[0][1] -= p[0][0] # x1 = x1 - x0
    p[0][2] -= p[0][0] # x2 = x2 - x0
    p[1][1] -= p[1][0] # y1 = y1 - y0
    p[1][2] -= p[1][0] # y2 = y2 - y0

    p[0][1] *= p[1][2] # x1 = (x1-x0) * (y2-y0)
    p[0][2] *= p[1][1] # x2 = (x2-x0) * (y1-y0)

    p[0][1] -= p[0][2] # r = (x1-x0) * (y2-y0) - (x2-x0) * (y1-y0)
    p[0][1] /= 2

    return p[0][1]

Answer 2

您可以使用 Numba 或 Cython 来优化此类问题。最好的方法可能是在从这些三角形中获取角点时直接计算三角形的有符号面积。即使使用单线程版本，您的内存带宽也可能成为那些少数加法和乘法的限制因素。

代码

import numba as nb
import numpy 

@nb.njit(fastmath=True)
def mix_nb(p):
    assert p.shape[0]==2
    assert p.shape[1]==3

    res=numpy.empty(p.shape[2],dtype=p.dtype)

    for i in range(p.shape[2]):
      res[i]=(+p[0,2,i] * (p[1,0,i] - p[1,1,i])+ p[0,0,i] * (p[1,1,i] - p[1,2,i])+ p[0,1,i] * (p[1,2,i] - p[1,0,i])) /2.

    return res

#Compile
p= np.random.rand(2, 3, 10000)
A=mix_nb(p)

import perfplot

#Benchmark
perfplot.show(
    setup=lambda n: np.random.rand(2, 3, n),
    kernels=[six, mix, mix2, cross, einsum, einsum2, einsum3,mix_nb],
    n_range=[2**k for k in range(19)],
    logx=True,
    logy=True,
    )

结果

Answer 3

这是 Cython 中的一个简单版本，为了完整起见，我将其包含在此处：

import numpy as np
from libc.stdlib cimport malloc, calloc, realloc, free

def calculate_signed_areas(double[:,:,::1] triangles):
    cdef:
        int i, n
        double area
        double[::1] areas
        double x0, x1, x2
        double y0, y1, y2

    n = triangles.shape[2]
    areas = <double[:n]>malloc(n * sizeof(double))
    for i in range(n):
        x0 = triangles[0,0,i]
        y0 = triangles[1,0,i]
        x1 = triangles[0,1,i]
        y1 = triangles[1,1,i]
        x2 = triangles[0,2,i]
        y2 = triangles[1,2,i]
        area = (
            + x2 * (y0 - y1)
            + x0 * (y1 - y2)
            + x1 * (y2 - y0)
        )/2.0
        areas[i] = area
    return np.asarray(areas)

Answer 4

我有点同意@TrapII 的回答：这是内存使用的问题。但是，我认为数据如何存储在内存中更为重要。让我们做以下实验：

In [1]: import numpy as np

In [2]: p = np.random.random((1000,3,2))

In [3]: p2 = np.array(np.moveaxis(p, (1, 0), (0, 2)).reshape((6, 1000)), order='C')

In [4]: def mix(p):
   ...:     return (
   ...:           p[:,2,0] * (p[:,0,1] - p[:,1,1])
   ...:         + p[:,0,0] * (p[:,1,1] - p[:,2,1])
   ...:         + p[:,1,0] * (p[:,2,1] - p[:,0,1])
   ...:         ) / 2.0
   ...:     

In [5]: def cross2(p2):
   ...:     e0x, e0y, e1x, e1y, e2x, e2y = p2
   ...:     return 0.5 * ((e1x - e0x)*(e2y - e0y) - (e1y - e0y)*(e2x - e0x))
   ...: 

In [6]: %timeit mix(p)
18.5 µs ± 351 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [7]: %timeit cross2(p2)
9.03 µs ± 90.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)