Python Numerical Integration for Volume of Region

Question

For a program, I need an algorithm to very quickly compute the volume of a solid. This shape is specified by a function that, given a point P(x,y,z), returns 1 if P is a point of the solid and 0 if P is not a point of the solid.

I have tried using numpy using the following test:

import numpy
from scipy.integrate import *
def integrand(x,y,z):
    if x**2. + y**2. + z**2. <=1.:
        return 1.
    else:
        return 0.
g=lambda x: -2.
f=lambda x: 2.
q=lambda x,y: -2.
r=lambda x,y: 2.
I=tplquad(integrand,-2.,2.,g,f,q,r)
print I

but it fails giving me the following errors:

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is assumed that the requested tolerance cannot be achieved, and that the returned result (if full_output = 1) is the best which can be obtained.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The integral is probably divergent, or slowly convergent.

So, naturally, I looked for "special-purpose integrators", but could not find any that would do what I needed.

Then, I tried writing my own integration using the Monte Carlo method and tested it with the same shape:

import random

# Monte Carlo Method
def get_volume(f,(x0,x1),(y0,y1),(z0,z1),prec=0.001,init_sample=5000):
    xr=(x0,x1)
    yr=(y0,y1)
    zr=(z0,z1)
    vdomain=(x1-x0)*(y1-y0)*(z1-z0)
    def rand((p0,p1)):
        return p0+random.random()*(p1-p0)
    vol=0.
    points=0.
    s=0.  # sum part of variance of f
    err=0.
    percent=0
    while err>prec or points<init_sample:
        p=(rand(xr),rand(yr),rand(zr))
        rpoint=f(p)
        vol+=rpoint
        points+=1
        s+=(rpoint-vol/points)**2
        if points>1:
            err=vdomain*(((1./(points-1.))*s)**0.5)/(points**0.5)
        if err>0:
            if int(100.*prec/err)>=percent+1:
                percent=int(100.*prec/err)
                print percent,'% complete\n  error:',err
    print int(points),'points used.'
    return vdomain*vol/points
f=lambda (x,y,z): ((x**2)+(y**2)<=4.) and ((z**2)<=9.) and ((x**2)+(y**2)>=0.25)
print get_volume(f,(-2.,2.),(-2.,2.),(-2.,2.))

but this works too slowly. For this program I will be using this numerical integration about 100 times or so, and I will also be doing it on larger shapes, which will take minutes if not an hour or two at the rate it goes now, not to mention that I want a better precision than 2 decimal places.

I have tried implementing a MISER Monte Carlo method, but was having some difficulties and I'm still unsure how much faster it would be.

So, I am asking if there are any libraries that can do what I am asking, or if there are any better algorithms which work several times faster (for the same accuracy). Any suggestions are welcome, as I've been working on this for quite a while now.

EDIT:

If I cannot get this working in Python, I am open to switching to any other language that is both compilable and has relatively easy GUI functionality. Any suggestions are welcome.

Answer 1

As the others have already noted, finding the volume of domain that's given by a Boolean function is hard. You could used pygalmesh (a small project of mine) which sits on top of CGAL and gives you back a tetrahedral mesh. This

import numpy
import pygalmesh
import meshplex


class Custom(pygalmesh.DomainBase):
    def __init__(self):
        super(Custom, self).__init__()
        return

    def eval(self, x):
        return (x[0]**2 + x[1]**2 + x[2]**2) - 1.0

    def get_bounding_sphere_squared_radius(self):
        return 2.0


mesh = pygalmesh.generate_mesh(Custom(), cell_size=1.0e-1)

gives you

From there on out, you can use a variety of packages for extracting the volume. One possibility: meshplex (another one out of my zoo):

import meshplex
mp = meshplex.MeshTetra(mesh.points, mesh.cells["tetra"])
print(numpy.sum(mp.cell_volumes))

gives

4.161777

which is close enough to the true value 4/3 pi = 4.18879020478.... If want more precision, decrease the cell_size in the mesh generation above.

Answer 2

Your function is not a continuous function, I think it's difficult to do the integration.

How about:

import numpy as np

def sphere(x,y,z):
    return x**2 + y**2 + z**2 <= 1

x, y, z = np.random.uniform(-2, 2, (3, 2000000))

sphere(x, y, z).mean() * (4**3), 4/3.0*np.pi

output:

(4.1930560000000003, 4.1887902047863905)

Or VTK:

from tvtk.api import tvtk

n = 151
r = 2.0
x0, x1 = -r, r
y0, y1 = -r, r
z0, z1 = -r, r
X,Y,Z = np.mgrid[x0:x1:n*1j, y0:y1:n*1j, z0:z1:n*1j]
s = sphere(X, Y, Z)

img = tvtk.ImageData(spacing=((x1-x0)/(n-1), (y1-y0)/(n-1), (z1-z0)/(n-1)), 
                     origin=(x0, y0, z0), dimensions=(n, n, n))

img.point_data.scalars = s.astype(float).ravel()

blur = tvtk.ImageGaussianSmooth(input=img) 
blur.set_standard_deviation(1)
contours = tvtk.ContourFilter(input = blur.output) 

contours.set_value(0, 0.5)
mp = tvtk.MassProperties(input = contours.output)
mp.volume, mp.surface_area

output:

4.186006622559839, 12.621690438955586

Answer 3

It seems to be very difficult without giving a little hint on the boundary:

import numpy as np
from scipy.integrate import *

def integrand(z,y,x):
    return 1. if x**2 + y**2 + z**2 <= 1. else 0.

g=lambda x: -2
h=lambda x: 2
q=lambda x,y: -np.sqrt(max(0, 1-x**2-y**2))
r=lambda x,y: np.sqrt(max(0, 1-x**2-y**2))
I=tplquad(integrand,-2.,2.,g,h,q,r)

print I