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

Question

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

我试过像这样用数字计算体积：

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)

...基本上我对整个立方体进行积分，每个点根据它是否在半个 space 内获得值 1 或 0。这变得非常慢（每次调用超过几秒钟）并不断给我警告，如

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

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

Answer 1

如果我们假设 half-space 的边界由 $\{(x, y, z) \mid ax + by + cz + d = 0 \}$ 和 $c \not= 0$，并且感兴趣的 half-space 是平面下方（在 $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$中至少有一个必须是non-zero，$c = 0$的情况可以通过改变坐标来处理。

Answer 2

整合

不连续函数的积分是有问题的，尤其是在多维度上。需要进行一些初步工作，将问题简化为连续函数的积分。在这里，我计算出高度 (top-bottom) 作为 x 和 y 的函数，并为此使用 dblquad：它 returns in 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)

Monte Carlo 方法

随机选取样本点（Monte Carlo 方法）避免了不连续性问题：不连续函数的精度与连续函数的精度大致相同，误差以 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 class for these but so far (1.0.0) does not implement finding the volume of such an object. If you could find the vertices of the polytope, then the ConvexHull class 可用于计算体积。但实际上，SciPy 空间模块似乎没有帮助。也许在 SciPy 的未来版本中...

Answer 3

特征函数的积分在数学上是正确的，但不实用。这是因为大多数积分方案都设计为在某种程度上准确地积分 多项式 ，因此所有“相对平滑”的函数都相当好。然而，特征函数一点也不平滑。 Polynomial-style 整合会让你一事无成。

一种非常better-suited的方法是先构建域的离散化版本，然后简单地求和小四面体的体积。

可以在 3D 中进行离散化，例如，使用 pygalmesh (a project of mine interfacing CGAL)。下面的代码将 cut-off 立方体离散化为

您可以通过降低 max_cell_circumradius and/or 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

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

Volume of 3d shape using numerical integration with scipy

geometry

scipy

numerical-integration

整合

Monte Carlo 方法

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