Python Numerical Integration for Region Volume

Question

For the program, I need an algorithm that allows you to quickly calculate the volume of a solid. This form is defined by a function that, for a given point P (x, y, z), returns 1 if P is a body point and 0 if P is not a body point.

I 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 this does not give me the following errors:

Warning (from the warning module): File "C: \ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py", line 321 warnings.warn (msg, IntegrationWarning) IntegrationWarning: The maximum number of units (50) was reached. If increasing the limit does not improve, it is recommended to analyze the integrand to determine the difficulties. If the position can be determined by local complexity (singularity, discontinuity) it is likely to benefit from dividing the interval and calling the integrator on the subbands. Perhaps you should use a special integrator.

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

Warning (from the warning module): File "C: \ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py", line 321 warnings.warn (msg, IntegrationWarning) IntegrationWarning: a rounding error was detected that prevents the required tolerance from achievements. The error may be underestimated.

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

So, of course, I was looking for "special purpose integrators", but could not find anything that could do what I needed.

Then I tried to write my own integration using the Monte Carlo method and tested it with the same form:

 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 it works too slow. For this program, I will use this digital integration about 100 times or so, and I will also do it on larger pieces, which will take minutes, if not an hour or two, at the speed that is being discussed now, not to mention what i want is better than 2 decimal places.

I tried implementing the MISER Monte Carlo method, but had some difficulties, and I'm still not sure how much faster it will be.

So, I ask if there are libraries that can do what I ask, or if there are any better algorithms that work several times faster (with the same accuracy). Any suggestions are welcome as I have been working on this for quite some time.

EDIT:

If I cannot get this to work in Python, I am open to switching to any other language that is compiled and has relatively lightweight GUI functionality. Any suggestions are welcome.

Answer 1

As others have already noted, it is difficult to find the volume of the region specified by the Boolean function. You can use pygalmesh (my small project), which is on top of CGAL and returns you a tetrahedral mesh. it

 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 now on, you can use various packages to retrieve the volume. One possibility: meshplex (another one from 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 you want more accuracy, reduce cell_size in the mesh generation above.

Answer 2

Your function is not a continuous function, it is difficult for me to integrate.

What 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 

exit:

 (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 

exit:

 4.186006622559839, 12.621690438955586 

Answer 3

This seems to be very complicated without giving a bit of a hint to the border:

 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