Python：3D直线插值及体素网格欧氏距离计算技术需求

Got it, let's tackle this 3D distance transformation problem step by step. Here's a practical, efficient approach using Python with scipy and skimage—perfect for your (21,41,91) voxel grid use case:

Step 1: Mark Voxels Covered by the 3D Line #

First, we need to accurately mark all voxels that the line passes through. Since your input points are on a straight line, we can use a 3D line-drawing algorithm (like Bresenham's) to get precise voxel indices instead of relying on dense interpolation.

import numpy as np
from skimage.draw import line_nd

# Example: Your input line points (replace with your actual data)
line_points = np.array([[5, 10, 15], [15, 30, 75]])  # Works with multiple sequential points too

# Define your voxel grid shape
grid_shape = (21, 41, 91)

# Initialize mask: True = non-line voxels, False = line voxels
mask = np.ones(grid_shape, dtype=bool)

# Iterate over each segment of the line
for i in range(len(line_points) - 1):
    start = line_points[i].astype(int)
    end = line_points[i+1].astype(int)
    # Get all voxels along the 3D line segment
    rr, cc, zz = line_nd(start, end)
    
    # Filter out indices that are outside the grid bounds
    valid_mask = (
        (rr >= 0) & (rr < grid_shape[0]) &
        (cc >= 0) & (cc < grid_shape[1]) &
        (zz >= 0) & (zz < grid_shape[2])
    )
    
    # Mark these voxels as part of the line (set to False)
    mask[rr[valid_mask], cc[valid_mask], zz[valid_mask]] = False

Step 2: Compute 3D Euclidean Distance Transformation #

Now use skimage's distance_transform_edt to calculate the Euclidean distance from every voxel to the nearest line voxel. This function is optimized for exactly this scenario—it's fast even for 3D grids.

from skimage.morphology import distance_transform_edt

# Calculate distance grid: line voxels will have distance 0
distance_grid = distance_transform_edt(mask)

# Now distance_grid is a (21,41,91) array where each value is the Euclidean distance to the line

Notes for Edge Cases & Optimization #

• If your line is curved (not straight): If your input points are on a curved path instead of a straight line, replace the line-drawing step with spline interpolation to generate dense points, then map those to voxels:

  from scipy.interpolate import splprep, splev
  
  # Interpolate to get dense points along the curve
  tck, u = splprep(line_points.T, s=0)  # s=0 ensures we pass through all input points
  dense_u = np.linspace(u.min(), u.max(), 2000)  # Adjust density based on grid size
  dense_points = np.array(splev(dense_u, tck)).T
  
  # Convert dense points to voxel indices (round to nearest voxel)
  voxel_indices = np.round(dense_points).astype(int)
  valid_indices = (
      (voxel_indices[:,0] >=0) & (voxel_indices[:,0] < grid_shape[0]) &
      (voxel_indices[:,1] >=0) & (voxel_indices[:,1] < grid_shape[1]) &
      (voxel_indices[:,2] >=0) & (voxel_indices[:,2] < grid_shape[2])
  )
  mask[voxel_indices[valid_indices,0], voxel_indices[valid_indices,1], voxel_indices[valid_indices,2]] = False
  

• Voxel coordinate mapping: If your voxel grid doesn't align with integer coordinates (e.g., voxels have a non-unit size or offset), adjust the point-to-voxel conversion by scaling/translating your input points first:

  # Example: If voxels are 2x2x2 units, with origin at (1,1,1)
  voxel_size = np.array([2,2,2])
  origin = np.array([1,1,1])
  voxel_indices = np.round((line_points - origin) / voxel_size).astype(int)
  

This approach is efficient, leverages optimized library functions, and gives you exactly the distance grid you need.

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