使用 scipy 数值积分的 3d 形状体积

Question

我已经编写了一个计算立方体和半空间相交体积的函数，现在我正在为它编写测试。

我试过像这样用数字计算音量：

integral = scipy.integrate.tplquad(lambda z, y, x: int(Vector(x, y, z).dot(normal) < distance),
                                   -0.5, 0.5,
                                   lambda x: -0.5, lambda x: 0.5,
                                   lambda x, y: -0.5, lambda x, y: 0.5,
                                   epsabs=1e-5,
                                   epsrel=1e-5)

...基本上我整合了整个立方体，每个点根据它是否在半空间内得到值 1 或 0。 这变得非常慢（每次调用超过几秒钟）并不断给我警告，例如

scipy.integrate.quadpack.IntegrationWarning: The integral is probably divergent, or slowly convergent

有没有更好的方法来计算这个体积？

Answer 1

特征函数的积分在数学上是正确的，但不实用。这是因为大多数积分方案都设计为在某种程度上精确地积分 多项式，因此所有“相对平滑”的函数都相当好。然而，特征函数并不平滑。多项式式集成将无济于事。

更适合的方法是首先构建域的离散版本，然后简单地总结小四面体的体积。

例如，可以使用pygalmesh（我的一个与CGAL 接口的项目）进行3D 离散化。下面的代码将截止立方体离散化为

您可以通过减小max_cell_circumradius 和/或max_edge_size_at_feature_edges 来提高精度，但网格化需要更长的时间。此外，您可以指定“特征边缘”来精确解决相交边缘。即使像元大小最粗，这也会为您提供完全正确的结果。

import pygalmesh
import numpy


c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1])
h = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 4.0, 10.0)
u = pygalmesh.Intersection([c, h])
mesh = pygalmesh.generate_mesh(
    u, max_cell_circumradius=3.0e-2, max_edge_size_at_feature_edges=1.0e-2
)


def compute_tet_volumes(vertices, tets):
    cell_coords = vertices[tets]
    a = cell_coords[:, 1, :] - cell_coords[:, 0, :]
    b = cell_coords[:, 2, :] - cell_coords[:, 0, :]
    c = cell_coords[:, 3, :] - cell_coords[:, 0, :]
    # omega = <a, b x c>
    omega = numpy.einsum("ij,ij->i", a, numpy.cross(b, c))
    # https://en.wikipedia.org/wiki/Tetrahedron#Volume
    return abs(omega) / 6.0

vol = numpy.sum(compute_tet_volumes(mesh.points, mesh.get_cells_type("tetra")))
print(f"{vol:.8e}")

8.04956436e-01

Answer 2

整合

不连续函数的积分是有问题的，尤其是在多维中。需要一些初步工作，将问题简化为连续函数的积分。在这里，我将高度（上下）作为 x 和 y 的函数，并使用dblquad：它返回36.2 ms。

我将平面方程表示为a*x + b*y + c*z = distance。需要注意 c 的符号，因为平面可能是顶部或底部的一部分。

from scipy.integrate import dblquad
distance = 0.1
a, b, c = 3, -4, 2  # normal
zmin, zmax = -0.5, 0.5  # cube bounds

# preprocessing: make sure that c > 0
# by rearranging coordinates, and flipping the signs of all if needed

height = lambda y, x: min(zmax, max(zmin, (distance-a*x-b*y)/c)) - zmin
integral = dblquad(height, -0.5, 0.5, 
                   lambda x: -0.5, lambda x: 0.5,
                   epsabs=1e-5, epsrel=1e-5)

蒙特卡洛方法

随机选取样本点（蒙特卡洛方法）避免了不连续性问题：不连续函数的精度与连续函数大致相同，误差以1/sqrt(N) 的速率减小，其中 N 是样本点的数量。

polytope package 在内部使用它。有了它，计算可以像

import numpy as np
import polytope as pc
a, b, c = 3, 4, -5  # normal vector
distance = 0.1 
A = np.concatenate((np.eye(3), -np.eye(3), [[a, b, c]]), axis=0)
b = np.array(6*[0.5] + [distance])
p = pc.Polytope(A, b)
print(p.volume)

这里 A 和 b 将半空间编码为Ax<=b：第一个固定行用于立方体的面，最后一个用于平面。

为了更好地控制精度，要么自己实现 Monte-Carlo 方法（简单），要么使用mcint 包（差不多一样简单）。

多面体体积：线性代数的任务，而不是积分器

你想计算一个多面体的体积，一个由相交的半空间形成的凸体。这应该有一个代数解。 SciPy 有 HalfspaceIntersection 类，但到目前为止（1.0.0）还没有实现查找此类对象的体积。如果你能找到多面体的顶点，那么ConvexHull 类可以用来计算体积。但事实上，SciPy 空间模块似乎无济于事。也许在 SciPy 的未来版本中......

Answer 3

如果我们假设半空间的边界由 $\{(x, y, z) \mid ax + by + cz + d = 0 \}$ 给出，$c \not= 0$，并且感兴趣的半空间是在平面下方（在 $z$ 方向），那么你的积分由下式给出

scipy.integrate.tplquad(lambda z, y, x: 1,
                        -0.5, 0.5,
                        lambda x: -0.5, lambda x: 0.5,
                        lambda x, y: -0.5, lambda x, y: max(-0.5, min(0.5, -(b*y+a*x+d)/c)))

由于 $a$、$b$ 和 $c$ 中至少有一个必须为非零，因此 $c = 0$ 的情况可以通过改变坐标来处理。