Great question! When working with sparse masks, using vtkSampleFunction to calculate normals for every voxel in the full 3D grid is a huge waste of resources—you're spending cycles on points you don't need. Instead, you can directly compute normals only for the masked points using vtkImplicitPolyDataDistance's built-in method for evaluating gradients, which maps directly to surface normals.

Step-by-Step Solution #

1. Generate Coordinates for Masked Points #

First, you need to get the 3D coordinates of exactly the points your mask specifies. Since vtkSampleFunction samples points on a regular grid aligned with your bounds, we can replicate that grid's points using NumPy, then filter them with your mask:

import numpy as np

# Your existing bounds and sample dimensions
x_min, x_max, y_min, y_max, z_min, z_max = bounds
nx, ny, nz = (your_sample_dimensions)

# Calculate step sizes (matches vtkSampleFunction's sampling)
dx = (x_max - x_min) / (nx - 1)
dy = (y_max - y_min) / (ny - 1)
dz = (z_max - z_min) / (nz - 1)

# Generate the full regular grid of points
x = np.linspace(x_min, x_max, nx)
y = np.linspace(y_min, y_max, ny)
z = np.linspace(z_min, z_max, nz)
xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
all_points = np.stack([xx.ravel(), yy.ravel(), zz.ravel()], axis=1)

# Filter to keep only masked points
masked_points = all_points[mask]

2. Compute Normals for Masked Points #

Use vtkImplicitPolyDataDistance's EvaluateFunctionAndGradient method to get the gradient (which is the unnormalized normal) for each masked point. We'll normalize it to match the unit normals returned by vtkSampleFunction:

import vtk

# Assume `implicit` is your pre-initialized vtkImplicitPolyDataDistance instance
normals = []
for point in masked_points:
    # Prepare variables to capture output
    distance = vtk.mutable(0.0)
    gradient = [0.0, 0.0, 0.0]
    
    # Compute distance and gradient (gradient = unnormalized normal)
    implicit.EvaluateFunctionAndGradient(point, distance, gradient)
    
    # Convert to NumPy array, normalize, and handle near-zero gradients
    grad_np = np.array(gradient)
    norm = np.linalg.norm(grad_np)
    if norm > 1e-8:  # Avoid division by zero for points far from the surface
        grad_np /= norm
    
    normals.append(grad_np)

# Convert list to NumPy array (matches your original `normals` output)
normals = np.array(normals)

3. Speed Up with Parallel Processing (Optional) #

If your mask still includes a large number of points, you can parallelize the computation to leverage multiple CPU cores. Note that VTK objects can be finicky with multiprocessing, so we'll recreate the implicit distance object in each process:

from multiprocessing import Pool

def compute_normal(args):
    point, poly_data = args
    # Reinitialize the implicit distance object in the process
    implicit = vtk.vtkImplicitPolyDataDistance()
    implicit.SetInput(poly_data)
    
    distance = vtk.mutable(0.0)
    gradient = [0.0, 0.0, 0.0]
    implicit.EvaluateFunctionAndGradient(point, distance, gradient)
    
    grad_np = np.array(gradient)
    norm = np.linalg.norm(grad_np)
    if norm > 1e-8:
        grad_np /= norm
    return grad_np

# Assume `poly_data` is your original input PolyData grid
with Pool() as p:
    # Pass both the point and PolyData to each worker
    normals = p.map(compute_normal, [(point, poly_data) for point in masked_points])

normals = np.array(normals)

Key Notes #

• Direction Alignment: The gradient returned by EvaluateFunctionAndGradient points away from the PolyData surface. If you find the normals are reversed compared to vtkSampleFunction's output, simply negate the gradient array (grad_np *= -1).
• Precision: The normalization step ensures you get unit normals, just like vtkSampleFunction produces. The 1e-8 threshold skips normalization for points extremely far from the surface where the gradient is negligible.
• Efficiency: This approach only computes normals for your masked points, so runtime scales with the size of your sparse mask instead of the full 3D grid—this can lead to massive speedups when your mask covers a small fraction of the grid.

内容的提问来源于stack exchange，提问作者Daniel Bichou