如何从球坐标中的点云创建体积？答案

【问题标题】：How to create volume from point cloud in spherical coordinates?如何从球坐标中的点云创建体积？
【发布时间】：2018-08-27 10:20:50
【问题描述】：

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

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

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

【问题讨论】：

凸包不起作用吗？
@DrBwts 我使用凸包创建了体积，但它没有正确区分体积内部和外部的点。这里是凸包算法的 [结果] (imgur.com/Z90Cqff)。在图像中：3 作为位于两个表面之外的点。
@GTR 你能解释一下我们在图表上看到的内容吗？蓝色和红色表面？网格？看起来顶部（蓝色？）和底部（红色？）表面实际上是球形但半径不同？半径测试是否不足以测试该点是否介于“蓝色”和“红色”之间？
把图片和解释放在OP里会很有帮助
@DrBwts 我的意思是在convex set 的意义上，“对于区域内的每一对点，连接这对点的直线段上的每个点也在区域内”。例如，沿着看起来像海岸线的球面的边界不是凸的：在此上的凸包算法将缩短所有重新进入的部分？

标签： numpy computational-geometry spherical-coordinate scipy-spatial

【解决方案1】：

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

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

【讨论】：

【解决方案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 方法。

【讨论】：