Вычислить подписанную площадь многих треугольников

Question

Мне нужно вычислить область со знаком из множества треугольников в 2D, заданную массивом фигур (2, 3, n) (координата x / y, узел в треугольнике, количество треугольников).Я ищу способ сделать это быстро, и лучшее, что я мог придумать, это

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

Другие попытки медленнее:

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.show(
    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)],
    logx=True,
    logy=True,
    )

Любые подсказки о том, каксделать его еще более эффективным?

Answer 1

Я несколько согласен с ответом @ 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)

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

Я думаю, что проблема не в вычислениях, а в том, как данные организованы в памяти.Если вам не нужно сохранять исходные данные, вы можете попытаться минимизировать количество временных объектов в памяти, используя функцию 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]