如何在内存中的3D Numpy数组高效填充指定三角形（Python）

Hey Luis,

Great question—filling a 3D triangle in such a large array efficiently is all about minimizing unnecessary computations and leveraging NumPy's vectorized operations instead of slow loops. Here's a step-by-step optimized approach that avoids iterating over your entire 1000x1200x1500 array:

Step 1: Define Core Geometry #

First, we'll compute the plane containing your triangle and handle edge cases like colinear points.

import numpy as np

# Your initial empty array
arrayXYZ = np.full((1000, 1200, 1500), np.uint8(0))

# Replace these with your actual triangle vertices
p1 = np.array([x1, y1, z1], dtype=int)
p2 = np.array([x2, y2, z2], dtype=int)
p3 = np.array([x3, y3, z3], dtype=int)

# Calculate the plane's normal vector (defines the triangle's orientation)
v1 = p2 - p1
v2 = p3 - p1
normal = np.cross(v1, v2)
a, b, c = normal

# Handle degenerate case: three points are colinear (no valid triangle)
if np.all(normal == 0):
    print("Points are colinear—no triangle to fill!")
    exit()

Step 2: Limit to Bounding Box #

Instead of checking every voxel in the huge array, we only work within the smallest box that encloses the triangle. This drastically cuts down the number of points we need to evaluate.

# Get min/max coordinates for each axis to define the bounding box
x_min, x_max = min(p1[0], p2[0], p3[0]), max(p1[0], p2[0], p3[0])
y_min, y_max = min(p1[1], p2[1], p3[1]), max(p1[1], p2[1], p3[1])
z_min, z_max = min(p1[2], p2[2], p3[2]), max(p1[2], p2[2], p3[2])

# Clamp bounds to your array's dimensions (avoids index errors if vertices are outside the array)
x_min = max(0, x_min)
x_max = min(arrayXYZ.shape[0]-1, x_max)
y_min = max(0, y_min)
y_max = min(arrayXYZ.shape[1]-1, y_max)
z_min = max(0, z_min)
z_max = min(arrayXYZ.shape[2]-1, z_max)

Step 3: Vectorized Plane & Triangle Membership Check #

We generate a grid of points within the bounding box, then filter down to only those that lie on the triangle's plane and inside the triangle itself. All operations use NumPy's optimized C-backed functions to avoid slow Python loops.

# Generate a grid of all points in the bounding box
x = np.arange(x_min, x_max + 1)
y = np.arange(y_min, y_max + 1)
z = np.arange(z_min, z_max + 1)
xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')  # Matches your array's (x,y,z) indexing

# Check which points lie on the triangle's plane (account for floating-point precision)
plane_eq = a * (xx - p1[0]) + b * (yy - p1[1]) + c * (zz - p1[2])
on_plane = np.isclose(plane_eq, 0, atol=1e-6)

# For points on the plane, check if they're inside the triangle using cross products
# Calculate vectors from triangle vertices to each grid point
ap = np.stack([xx - p1[0], yy - p1[1], zz - p1[2]], axis=-1)
bp = np.stack([xx - p2[0], yy - p2[1], zz - p2[2]], axis=-1)
cp = np.stack([xx - p3[0], yy - p3[1], zz - p3[2]], axis=-1)

# Compute cross products and dot with the normal to check orientation consistency
cross_ab_ap = np.cross(v1, ap)
dot1 = np.dot(cross_ab_ap, normal)
cross_bc_bp = np.cross(p3 - p2, bp)
dot2 = np.dot(cross_bc_bp, normal)
cross_ca_cp = np.cross(p1 - p3, cp)
dot3 = np.dot(cross_ca_cp, normal)

# Points are inside if all dot products have the same sign (or zero, for edges)
inside = ((dot1 >= 0) & (dot2 >= 0) & (dot3 >= 0)) | ((dot1 <= 0) & (dot2 <= 0) & (dot3 <= 0))

# Combine masks to get only valid voxels that are on the plane and inside the triangle
valid_voxels = on_plane & inside

# Extract indices of valid voxels
x_idx = xx[valid_voxels]
y_idx = yy[valid_voxels]
z_idx = zz[valid_voxels]

Step 4: Assign Values to Valid Voxels #

Finally, we update the array in one vectorized operation—no loops required.

# Set valid voxels to 255
arrayXYZ[x_idx, y_idx, z_idx] = np.uint8(255)

Key Optimizations #

• Bounding Box Reduction: We only process voxels near the triangle, skipping the vast majority of your array.
• Vectorized Operations: All checks use NumPy's optimized functions, which are orders of magnitude faster than Python loops.
• Precision Handling: np.isclose avoids missing points due to tiny floating-point errors in plane calculations.

If you need even more memory efficiency (e.g., if the bounding box is still large), we can solve the plane equation for one axis (e.g., compute z from x and y if the normal's z-component isn't zero) to work with a 2D grid instead of 3D. Just let me know!

内容的提问来源于stack exchange，提问作者Luis Gonçalves