2

Given a set of data (x,y,z) with z>0 (around 10^4 points), I want to find the volume under the surface it describes using python. All these data points are randomly generated and unsorted, but from the 3D plot the surface is smooth and vanishing. See here plotted data

So far, I have checked the following related topics
Using Convex Hull : _Among other errors this one: QhullError: QH6114 qhull precision error: initial simplex is not convex.
Trapezoidal: Although it runs, I don't trust in the result, since even generating half sphere shell the volume was wrong.
Other Convex Hull: Same errors as before.

I also though I could make a grid and interpolate and then compute the double integral

# Given independents numpy.ndarray x2, y2, and z2
# Set up a regular grid of interpolation points
x1i = np.linspace(x2.min(), x2.max(), 1000)
y1i = np.linspace(y2.min(), y2.max(), 1000)
x1i, y1i = np.meshgrid(x1i, y1i)

# Interpolate; 
z1i = scipy.interpolate.griddata((x2, y2), z2, (x1i, y1i), method='linear')

But then I don't know how to implement dblquad given the griddata

integrate.dblquad(z1i, min(x2),max(x2),lambda x: min(y2), lambda x: max(y2), epsabs=1.49e-08, epsrel=1.49e-08)

doesn't work.

How can I do it?

Daniel Santos
  • 14,328
  • 21
  • 91
  • 174
ktitimbo
  • 73
  • 6
  • Is the surface radially symmetric (as it appears in the image?). The whole problem would be much simpler if you could reformulate in a 2-D z-r space and the integrate over 2pi radians. – Daniel F Nov 16 '17 at 14:03
  • For this realization its projection on xy is a circumference, however for other experiments it turns to be an ellipsis. – ktitimbo Nov 16 '17 at 14:37
  • Maybe I'm missing something. If you already made the 2D interpolation you can calculate the volume at any given cell by multiplying the height by the square of the cell size (if all cells are above 0). The result from griddata is not a function (which would be the input for dblquad). – armatita Nov 16 '17 at 15:48

2 Answers2

0

I suggest you to use Monte Carlo method for this kind of application. Since you know the data points , it is easy to use Monte Carlo method compared to double integrals. One caveat is that, this only gives approximate volume under the surface. With increasing the number of random points you can increase the accuracy. I would suggest you to implement this in google colab as it requires high computation power. For more info:https://www.youtube.com/watch?v=6dprNJVaGPA&t=2s

Bharat Prem
  • 1
0

This is an old post by maybe someone is still interested.

The volume can be computed by summing up the volumes of all tetrahedra of the Delaunay triangulation. This is efficiently computed using the determinant of three corner vectors which is the signed (could be neg!) volume of the resulting parallelepiped. Turns out the desired tetrahedral volume is just 1/6 of that.

(see https://www.youtube.com/watch?v=xgGdrTH6WGw)

This gives the volume of the convex hull. I had a more complex problem where I had a cubic volume of points which I transformed into another space in which the outer surface was not a complex hull. However, since I had both sets of points, I compute the complex hull in the original cube, then use those tetrahedra to compute the volume. (Note this assumes the transformation between GT and GP doesn't distort the points too much so as to change their order.) If you just want the convex hull volume, GT and GP are the same in the code below.

def getGamutVol(GT, GP): # Calculates volume of a set of points. ASSUMES SAME SPATIAL ORDER IN SPACE IN GT and GP # GT [N x 3] Points from which to get triangulation # GP [N x 3] Points on which to calc volume

tri             = Delaunay(GT).simplices
p0              = GP[tri[:,0],:]
p1              = GP[tri[:,1],:]
p2              = GP[tri[:,2],:]
p3              = GP[tri[:,3],:]
                                        # Form vectors between triangle corners
v1              = p1-p0
v2              = p2-p0
v3              = p3-p0
                                        # Calc triple products (signed volume of parallelpiped)
dets            = np.linalg.det(np.dstack([v1, v2, v3]))
vol             = np.sum(np.abs(dets))/6

return vol
Kevin Fine
  • 1