Python Numerical Integration for Volume of Region

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 numpyimport pygalmeshimport meshplexclass 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.0mesh = 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 meshplexmp = 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.

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

How about:

import numpy as npdef sphere(x,y,z):    return x**2 + y**2 + z**2 <= 1x, 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 tvtkn = 151r = 2.0x0, x1 = -r, ry0, y1 = -r, rz0, z1 = -r, rX,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

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

import numpy as npfrom scipy.integrate import *def integrand(z,y,x):    return 1. if x**2 + y**2 + z**2 <= 1. else 0.g=lambda x: -2h=lambda x: 2q=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