使用scipy进行数值积分的3d形状体积 [英] Volume of 3d shape using numerical integration with scipy

问题描述

我已经编写了一个用于计算立方体和半个空间的交点体积的函数，现在我正在为其编写测试.

I have written a function for computing volume of intersection of a cube and a half-space and now I'm writing tests for it.

我试图像这样用数字计算体积:

I've tried computing the volume numerically like this:

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.这变得非常慢(每次调用要超过几秒钟)，并且不断向我发出类似

... basically I integrate over the whole cube and each point gets value 1 or 0 based on if it is inside the half space. This gets very slow (more than several seconds per invocation) and keeps giving me warnings like

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

是否有更好的方法来计算此体积?

Is there a better way to calculate this volume?

推荐答案

集成

不连续功能的集成是有问题的，尤其是在多维方面.需要一些初步的工作，将问题减少到一个连续函数的整数中.在这里，我计算出高度(上下)与x和y的关系，并为此使用 dblquad :它以 36.2 ms 返回.

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

I express the plane equations as a*x + b*y + c*z = distance. Some care is needed with the sign of c, as the plane could be a part of the top or of the bottom.

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软件包在内部使用它.有了它，计算就可以像

The polytope package uses it internally. With it, a computation could go as

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 :前一个固定行用于立方体的面，最后一个固定行用于平面.

Here A and b encode the halfspaces as Ax<=b: the first fix rows are for faces of the cube, the last is for the plane.

要对精度进行更多控制，请自己(轻松)实现Monte-Carlo方法，或使用>code> mcint 包(大致一样简单).

To have more control over precision, either implement Monte-Carlo method yourself (easy) or use mcint package (about as easy).

您要计算多面体的体积，该多面体由相交的半空间形成.这应该有一个代数解.SciPy具有 HalfspaceIntersection 类，但到目前为止(1.0.0)尚未实现查找此类对象的体积.如果您可以找到多面体的顶点，则 ConvexHull 类可用于计算体积.但是就目前而言，SciPy空间模块似乎无济于事.也许在SciPy的未来版本中...

You want to compute the volume of a polytope, a convex body formed by intersecting halfspaces. This ought to have an algebraic solution. SciPy has 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 could be used to compute the volume. But as is, it seems that SciPy spatial module is no help. Maybe in a future version of SciPy...

这篇关于使用scipy进行数值积分的3d形状体积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！