numpy - 如何从球坐标中的点云创建体积？

Question

我在球坐标中有两组离散点，每组代表对象的顶面和底面。

我正在尝试从这些点创建体积，以分离位于对象内部和外部的点。有什么建议在哪里看或使用哪个库？

蓝色和红色点代表顶部和底部表面。红点是通过以一定的恒定半径径向向下移动顶面产生的。

Answer 1

如果我是对的，蓝色和红色表面是网格化的（并且是防水的）。因此，对于每个点，您都可以从球体中心画线并寻找与网格的交点。这是通过找到两个三角形使线穿过它们来完成的（这可以通过仅查看角坐标来完成，使用三角形中的点公式），然后找到交点。然后将点分类为红色表面之前、蓝色表面之后或两者之间是一件容易的事。

对三角形的详尽搜索可能代价高昂。例如，您可以使用边界框或类似设备的层次结构来加速它。

Answer 2

这是一种自定义的修补方法，它可以在原始表面中点之间的平均距离远小于体积的厚度和表面轮廓上的不规则度的情况下工作。换句话说，有很多点描述了蓝色表面。

import matplotlib.pylab as plt

import numpy as np
from scipy.spatial import KDTree

# Generate a test surface:
theta = np.linspace(3, 1, 38)
phi = np.zeros_like(theta)
r = 1 + 0.1*np.sin(8*theta)

surface_points = np.stack((r, theta, phi), axis=1)  #  n x 3 array

# Generate test points:
x_span, y_span = np.linspace(-1, 0.7, 26), np.linspace(0.1, 1.2, 22)
x_grid, y_grid = np.meshgrid(x_span, y_span)

r_test = np.sqrt(x_grid**2 + y_grid**2).ravel()
theta_test = np.arctan2(y_grid, x_grid).ravel()
phi_test = np.zeros_like(theta_test)

test_points = np.stack((r_test, theta_test, phi_test), axis=1)  #  n x 3 array

# Determine if the test points are in the volume:

volume_thickness = 0.2  # Distance between the two surfaces
angle_threshold = 0.05  # Angular threshold to determine for a point 
                        # if the line from the origin to the point 
                        # go through the surface

# Get the nearest point: (replace the interpolation)
get_nearest_points = KDTree(surface_points[:, 1:]) # keep only the angles
# This is based on the cartesian distance,
# and therefore not enterily valid for the angle between points on a sphere
# It could be better to project the points on a unit shpere, and convert 
# all coordinates in cartesian frame in order to do the nearest point seach...

distance, idx = get_nearest_points.query(test_points[:, 1:])

go_through = distance < angle_threshold

nearest_surface_radius = surface_points[idx, 0]

is_in_volume = (go_through) & (nearest_surface_radius > test_points[:, 0]) \
                & (nearest_surface_radius - volume_thickness < test_points[:, 0])

not_in_volume = np.logical_not(is_in_volume)

# Graph;
plt.figure(figsize=(10, 7))
plt.polar(test_points[is_in_volume, 1], test_points[is_in_volume, 0], '.r',
          label='in volume');
plt.polar(test_points[not_in_volume, 1], test_points[not_in_volume, 0], '.k',
          label='not in volume', alpha=0.2);
plt.polar(test_points[go_through, 1], test_points[go_through, 0], '.g',
          label='go through', alpha=0.2);
plt.polar(surface_points[:, 1], surface_points[:, 0], '.b',
          label='surface');
plt.xlim([0, np.pi]); plt.grid(False);plt.legend();

二维情况下的结果图是：

这个想法是通过仅考虑方向而不是半径来寻找表面中最近的点的每个测试点。一旦找到这个“相同方向”的点，就可以测试该点是否沿径向 ( volume_thickness) 位于体积内部，以及使用参数 是否足够靠近表面angle_threshold。

我认为最好对蓝色表面进行网格化（非凸）并执行适当的插值，但我不知道 Scipy 方法。