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: 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
    

