布尔联合两个网格仅适用于椭圆体

Question

我想在街道（平面）中创建一个坑洞（椭圆体）。经过一些测试，我意识到布尔联合可能是一个可以使用的（因为我想结合两个网格但不希望街道掩盖坑洼）。
当我使用 Sphere 尝试这个时，它工作得很好。
但是，我使用的是 ParameticEllipsoid。在那里我经常收到错误“vtkMath::Jacobi：提取特征函数时出错”，并且想要的效果仅适用于布尔差值而不是并集。我想它必须与 vtk/py 及其几何对象创建有关。

虽然采用布尔差值会带来想要的结果，但它会花费很多时间（我不会介意太多）并且很少起作用（并且仅当椭圆体未与其他对象啮合时）。

有没有办法避免这种情况？我的想法是提取网格的点（NumPy 数组）并计算它们的并集，但是我做不到。

def add_pothole():
    street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=100, j_resolution=100)

    points = street_mesh.points
    for i in range(4):
        pothole_mesh = pv.ParametricEllipsoid(10,5,3)
        # alternatively to get the half ellipsoid: (but boolean difference seemed to only work with the closed one
        # pothole_mesh = pv.ParametricEllipsoid(10,5,3, max_v=pi /2)
        pothole_mesh.translate(choice(points))
    
        # "add" the pothole to the street
        street_mesh = pothole_mesh.boolean_difference(street_mesh.triangulate())
    
    street_mesh.plot()

或者，我考虑直接定义椭圆形凹陷，使用高度作为平面内位置的函数。是否可以在不使用布尔运算的情况下执行此操作？

Answer 1

为了解决您当前的问题，您当前的代码存在两个问题：

与Sphere不同，ParametricEllipsoid由于VTK构造它的方式有很多杂散点，所以最终结果在技术上不是一个封闭的曲面。这是我计划在 pyvista 中修复的已知问题。现在你可以通过在你的椭圆体上调用 clean(tolerance=1e-5) 来解决这个问题。这应该会阻止那些 vtkMath 错误的弹出。
在每次调用 boolean_difference 后，从坑洞中减去街道将反转平面的法线，这将使您在平面的交替两侧出现颠簸，这可能不是您想要的。您可以通过在循环结束时调用 street_mesh.flip_normals() 来解决此问题。

通过这两项更改，您的代码将运行，尽管速度很慢（我真的不知道为什么，但我也不知道布尔运算在幕后是如何工作的）：

from random import choice
import numpy as np
import pyvista as pv

def add_pothole():
    street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=100, j_resolution=100)
    pothole_template = pv.ParametricEllipsoid(10, 5, 3).clean(tolerance=1e-5).triangulate(

    points = street_mesh.points
    for i in range(4):
        pothole_mesh = pothole_template.copy()
        pothole_mesh.translate(choice(points))

        # "add" the pothole to the street
        street_mesh = pothole_mesh.boolean_difference(street_mesh.triangulate())
        street_mesh.flip_normals()

    street_mesh.plot()

add_pothole()

但是要解决您的根本问题，如果您的街道一开始就是 z(x, y) 表面，则您不必一开始就这样做。考虑到 Perlin 噪声，您可能从一个带有从 Perlin 噪声采样的标量的平面开始，然后调用 warp_by_scalar() 来实现噪声调制。好吧，你可以只在相同的标量上添加椭圆形的凸起，然后在最后一步进行扭曲（与完全相同的方式）。

为此，您必须计算椭圆体的 z(x, y) 参数化。这不是微不足道的，但也不是很难。更广为人知的是，围绕原点半径为 R 的球体的方程是

x^2 + y^2 + z^2= R^2.

除以 R^2:

更好（并且规范）

(x/R)^2 + (y/R)^2 + (z/R)^2 = 1.

获得椭圆体的方法是线性变换每个轴，改变它们各自的半径（因此它们不再都是 R）。这是椭圆体的隐式方程：

(x/a)^2 + (y/b)^2 + (z/c)^2 = 1.

如果你想要椭圆体有一个非平凡的中心，你必须平移每个坐标：

(x - x0)^2/a^2 + (y - y0)^2/b^2 + (z - z0)^2/c^2 = 1.

这是 3 维椭球体的一般形式，假设其轴与笛卡尔坐标系的轴定向。

嗯，这很好，因为我们可以重新排列它以获得

(z - z0)^2/c^2 = 1 - (x - x0)^2/a^2 - (y - y0)^2/b^2
(z - z0)^2 = c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2
z = z0 +- sqrt(c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2)

正负部分对应于椭圆体不是 z(x, y) 函数的事实，因为对于每个 (x, y) 点，您有两个可能值。但是你可以用两个函数构造一个椭圆体：一个顶面和一个底面。

最后回到你的问题，你可以随机选择一个 (x0, y0, z0) 点，然后选择上面椭圆体的底面（z = z0 - sqrt(...) 的那个）。您唯一需要担心的是，对于没有椭圆体的 (x, y) 点，您会得到一个虚数平方根。所以你必须先过滤掉椭圆体内部的点。执行此操作的最简单方法是计算平方根，并丢弃 NaN：

import numpy as np
import pyvista as pv

# denser plane for sharper ellipsoid
street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=1000, j_resolution=1000)
street_mesh['Normals'] *= -1  # make normals point up for warping later
# add Perlin noise for example
freq = [0.03, 0.03, 0.03]
noise = pv.perlin_noise(5, freq, [0, 0, 0])
street_mesh['height'] = [noise.EvaluateFunction(point) for point in street_mesh.points]

# the relevant part: ellipsoid function
x0, y0, z0 = street_mesh.points[street_mesh.points.size//6, :]  # or whatever random point you want
x, y, _ = street_mesh.points.T  # two 1d arrays of coordinates
a, b, c = 10, 5, 3  # semi-axes
ellipsoid_fun = z0 - np.sqrt(c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2)  # RuntimeWarning
keep_inds = ~np.isnan(ellipsoid_fun)
street_mesh['height'][keep_inds] += ellipsoid_fun[keep_inds]
street_mesh.set_active_scalars('height')

# warp by 'height' and do a quick plot
street_mesh.warp_by_scalar().plot()

以上将发出有关虚数平方根的警告（导致数据中出现 NaN）。如果这让您感到困扰，您可以明确地将其静音。或者，我们可以按照 LBYL 并在计算平方根之前自己检查值：

arg = c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2
keep_inds = arg >= 0
ellipsoid = z0 - np.sqrt(arg[keep_inds])
street_mesh['height'][keep_inds] = ellipsoid

这是结果，使用增加的平面密度获得更尖锐的椭圆体：

（我最初建议您可以使用 numpy 广播计算一次坐下的每个坑洞，但这实际上并不容易，甚至不可能，因为 NaN 过滤打破了这一点。）

布尔联合两个网格仅适用于椭圆体

Boolean union two meshes only works for ellipsoid

python

3d

open3d

pyvista