numpy阵列索引的scanline算法三角网格化

Question

我需要访问位于一个小三角形区域的大型numpy数组的元素。使用matplotlib.path.contains_points或ImageDraw.Draw.Polygon的强力解决方案太慢了，因为它们返回一个布尔数组，而不是访问相关元素的索引列表。

因此，我决定实现一个scanline算法，该算法对三角形进行栅格化，并返回三角形内元素的索引数组。我在阳光“三角星光化”页面上使用了“标准算法”的出色描述，并进行了三次更改：

• 三角形的顶点不是返回数组的一部分(因为我不想在后续的处理步骤中使用这些顶点)
• 三角形所有边上的点是返回数组的一部分。这意味着，如果我处理两个相邻的三角形，边缘点将是两个三角形的一部分。我这么做是因为在我的例子中，重要的是不要错过任何点，所以我宁愿处理两次。
• 我没有对平底三角形和平顶三角形使用两个单独的函数，而是使用了一个向上或向下移动的函数。

代码可以工作，但我不确定这是否是一个特别优雅的实现。我有两个主要问题：

• 对于三角形顶点这样的小结构，我应该使用numpy数组和相关的函数(如argsort )，还是使用元组和列表来实现更有效(更好的实践)？
• 我使用嵌套的for循环作为算法的“心脏”。我知道，只有循环被矢量化时numpy才会发光，但在这种情况下，我不知道该如何做。对这个循环进行矢量化是明智的，如果是，我该如何做？
• 有比我目前使用的rasterize_horizontal_triangle表达式更优雅的方法来排序tri[1, 0], tri[2, 0] = tri[2, 0], tri[1, 0]中的两个水平顶点吗？

我不是一个非常有经验的程序员，更不用说专业的程序员了，所以我欢迎对代码的任何反馈，不管是关于风格、设计还是效率。

triangle_rasterization.py：

# Use real division everywhere
from __future__ import division

import numpy as np


def rasterize_triangle(tri):
    """
    Given a 3x2 numpy array TRI describing the integer vertices of a general
    triangle, return an array containing all the points that lie within this
    triangle or on the triangle's edge, but not the triangle vertices
    themselves.
    This code is based on the description given in
    http://www.sunshine2k.de/coding/java/TriangleRasterization/TriangleRasterization.html
    """
    # Sort by increasing y coordinate
    tri = tri[tri[:, 1].argsort()]

    # Check for triangles with horizontal edge
    if tri[1, 1] == tri[2, 1]:
        # Bottom is horizontal
        points = rasterize_flat_triangle(tri)
    elif tri[0, 1] == tri[1, 1]:
        # Top is horizontal
        points = rasterize_flat_triangle(tri[(2, 0, 1), :])
    else:
        # General triangle.
        # We'll split this into two triangles with horizontal edges and process
        # them separately.
        # Find the additional vertex that splits the triangle.
        helper_point = np.array([tri[0, 0] + (tri[1, 1] - tri[0, 1]) /
                                         (tri[2, 1] - tri[0, 1]) *
                                         (tri[2, 0] - tri[0, 0]),
                                         tri[1, 1]]).round()
        # Top triangle
        points = rasterize_flat_triangle(tri[(0, 1), :],
            helper_point=helper_point)
        # Bottom triangle
        points = np.vstack([points, rasterize_flat_triangle(tri[(2, 1), :],
            helper_point=helper_point)])
    return points


def rasterize_flat_triangle(tri, helper_point=None):
    '''
    Given a 3x2 numpy array TRI describing the vertices of a triangle where the
    second and third vertex have the same y coordinate, return an array
    containing all the points that lie within this triangle or
    on the triangle's edge, but not the triangle vertices themselves.
    Or, given a 2x2 numpy array TRI containing two vertices and HELPER_POINT
    containing the third vertex, again return the same points as before, but
    additionally return the helper_point as well (used when treating a
    general triangle that's split into two triangles with horizontal edges)
    '''
    # Is the triangle we're treating part of a split triangle?
    if helper_point is not None:
        tri = np.vstack([tri, helper_point])

    # Is the bottom or the top edge horizontal?
    ydir = np.sign(tri[1, 1] - tri[0, 1])

    # Make sure that the horizontal edge is left-right oriented
    if tri[1, 0] > tri[2, 0]:
        tri[1, 0], tri[2, 0] = tri[2, 0], tri[1, 0]

    # Find the inverse slope (dx/dy) for the two non-horizontal edges
    invslope1 = ydir * (tri[1, 0] - tri[0, 0]) / (tri[1, 1] - tri[0, 1])
    invslope2 = ydir * (tri[2, 0] - tri[0, 0]) / (tri[2, 1] - tri[0, 1])

    # Initialize the first scan line, which is one y-step below or above the
    # first vertex
    curx1 = tri[0, 0] + invslope1
    curx2 = tri[0, 0] + invslope2
    points = []

    # Step vertically. Don't include the first row, because that row only
    # contains the first vertex and we don't want to return the vertices
    for y in np.arange(tri[0, 1] + ydir, tri[1, 1], ydir):
        for x in np.arange(curx1.round(), curx2.round() + 1):
            points.extend([(x, y)])
        curx1 += invslope1
        curx2 += invslope2

    # If we're dealing with the first half of a split triangle, add the
    # helper point (because that's not a "real" vertex of the triangle)
    if helper_point is not None and ydir == 1:
        points.extend([tuple(helper_point)])

    # If we're not dealing with a split triangle, or if we're dealing with the
    # first half of a split triangle, add the last line (but without the end
    # points, because they're the vertices of the triangle
    if helper_point is None or ydir == 1:
        for x in np.arange(tri[1, 0] + 1, tri[2, 0]):
            points.extend([(x, tri[1, 1])])

    return np.array(points, dtype='int')

下面是一个小测试脚本：

import triangle_rasterization as tr
import matplotlib.pyplot as plt

triangleA = np.array([[5,15],[5,1],[14,8]])
pointsA = tr.rasterize_triangle(triangleA)

triangleB = np.array([[14,8],[5,1],[18, 1]])
pointsB = tr.rasterize_triangle(triangleB)

triangleC = np.array([[5,15],[14,15],[14,8]])
pointsC = tr.rasterize_triangle(triangleC)

array = np.zeros([20,20])
array[pointsA[:,1], pointsA[:,0]] = 1
array[pointsB[:,1], pointsB[:,0]] = 2
array[pointsC[:,1], pointsC[:,0]] = 3

plt.imshow(array, interpolation='none')
plt.scatter(*triangleA.T, c='white')
plt.scatter(*triangleB.T, c='white')
plt.scatter(*triangleC.T, c='white')

Answer 1

1. 在rasterize_triangle中，docstring说:返回一个包含所有…点的数组，但是我认为写得更清楚:返回一个包含点…坐标的数组
2. 在Numpy中，返回坐标的函数通常更方便地返回数组元组，而不是多维数组。例如，参见numpy.diag_indices或numpy.triu_indices。这是因为您可以使用一个坐标数组的元组来索引一个数组：>>> a= np.arange(16).reshape(4,4) >>> np.diag_indices(4，4) ( array (0，1，2，3)，array(0，1，2，3)) >>> a_数组(0、5、10、15)，返回n×2数组不太方便:必须对结果调用numpy.unravel_index，然后才能使用它来索引数组。
3. 在rasterize_triangle中，docstring说:这段代码基于…文档中给出的描述，最好是从用户的角度编写(我如何使用这个函数？我要通过什么论据？)，但这句话是从实现者的角度来看的，最好把它写成评论。
4. 函数rasterize_flat_triangle有一个长而复杂的docstring，这表明它的接口不太正确。调用此函数的方法有两种:传递2×2数组和辅助点(在这种情况下，数组和辅助点被堆叠成三角形，如果底部是水平的，则包含辅助点的扫描线包括在结果中)，或者传递3×2数组而没有帮助点。如果您让调用方担心输入的堆叠以及包含或排除包含辅助点的scanline，它将简化界面。这只发生在一个地方(在rasterize_triangle中，三角形是分裂的)，所以这将是一个整体简化。
5. 没有示例或测试用例。这是一种自我包含的函数，在这种功能中，doctest可以很好地工作。
6. 辅助点的x坐标是使用以下表达式计算的: tri0，0 + (tri1，1 - tri0，1) / (tri2，1 - tri0，1) * (tri2，0 - tri0，0)，但我认为numpy.interp会使意图变得更清楚(这一点是沿着点0和2之间的直线插值的)：np.interp(tri1，1，tri(0，2)，1，tri(0，2)，0)。
7. 三角形的底部是水平的(如果有必要的话)如下: tri(2，0，1)，但是我认为numpy.roll会使这个操作更清晰:np.roll(tri1，0)
8. 要使rasterize_flat_triangle工作，必须是第二顶点和第三顶点具有相同的y坐标，因此我将断言如下: assert 1，1 == tri2，1
9. 在rasterize_flat_triangle中，结果作为np.array(points, dtype='int')返回。这使得将此函数与其他数据类型(如'int64' )一起使用变得很尴尬。考虑将结果的数据类型作为关键字参数(默认为'int')。