# Volume of Region Contained by 3D Data Points

I have a dataset consisting of (x,y,z) data that describe a surface, and I need to find the volume of the region bounded by that surface. The data do not follow any mathematical equation, but the volume is roughly of the shape of a parabola that has been rotated about the vertical axis. I have created dummy data to here to represent my scenario:

xydata = Table[{x, x^2 + 8*Cos[x]}, {x, 0, 10, 0.2}];
xyzdata = Flatten[Table[{xydata[[n, 1]]*Cos[\[Theta]],xydata[[n, 1]]*Sin[\[Theta]], xydata[[n,2]]}, {n,1,Length[xydata]}, {\[Theta], 0, 2 Pi, Pi/30}], 1];
xyzdata = DeleteDuplicates[xyzdata];


These example data looks like this:

ListPlot3D[xyzdata]


How do I find the volume enclosed by this dataset? I am suspecting that interpolation might be the easiest path forward,

int = Interpolation[xyzdata,InterpolationOrder->1]


and then numerically integrating the interpolated function, but is there a clever way to find the volume using the Volume function or something similar?

• Do you mean the volume of the object with the hole closed, or the volume of the region under the surface? Commented Sep 13, 2023 at 20:02
• You might try creating an ImplicitRegion where z is between this surface and the maximal z and then taking the Volume or RegionMeasure of that region. Commented Sep 13, 2023 at 20:41

# Volume under the surface

{zbot, ztop} = CoordinateBounds[xyzdata][[3]];

int = Interpolation[xyzdata, InterpolationOrder -> 1];

vbot = NIntegrate[
Boole[x^2 + y^2 < 10^2] (int[x, y] - zbot),
{x, -10, 10}, {y, -10, 10}, Method -> "LocalAdaptive"]

15414.7


# Volume of the enclosing surface

### Method 1

Volume[Cylinder[{{0, 0, zbot}, {0, 0, ztop}}, 10]] - vbot

13936.1


### Method 2

mr = RepairMesh[DiscretizeGraphics[ListPlot3D[xyzdata, Mesh -> None]], "HoleEdges"];
bmr = BoundaryMeshRegion[MeshCoordinates[mr], MeshCells[mr, 2]];
Volume[bmr]

13936.2

• Thank you for this reply, it looks super helpful. For method 1, you set the inequality and integration limits around 10 because that's what the example data were generated to, correct? If possible, could you edit your answer to use data-defined limits rather than hardcoded numbers? Thank you again!
– John
Commented Sep 15, 2023 at 12:41
• I also assume then that the membership condition in the Boole needs to be arbitrary too? Commented Sep 15, 2023 at 15:09
• If I understand how that code works then yes I think it would need to be arbitrary.
– John
Commented Sep 15, 2023 at 16:54
• Any way we could access your actual data points? Commented Sep 15, 2023 at 17:49
• Unfortunately I am not allowed to do that. The Boole inequality doesn't work for me because the XY radius is not quite as uniform as the example dataset shows. However, the data do have a very clear well-defined limit at the upper z bound. If I just shift the entire dataset down so that the top of the dataset is at z=0 and then integrate the interpolation without the Boolean, that would find the volume of the shape correctly, wouldn't it?
– John
Commented Sep 18, 2023 at 14:20

Just to illustrate closeness of other answers to analytic answer for this volume of revolution:

g[x_] := x^2 + 8 Cos [x];
p[u_, v_] := {u Cos[v], u Sin[v], g[u]};
below = Integrate[
Cross[D[p[u, v], u], D[p[u, v], v]] . {0, 0, g[u]}, {u, 0, 10}, {v,
0, 2 Pi}]
min = x /. Minimize[{g[x], 0 <= x <= 10}, x][[2]];
cylinder = Pi 10^2 (g[10] - g[min]) // N
below // N
cylinder - below

Show[ParametricPlot3D[p[u, v], {u, 0, 10}, {v, 0, 2 Pi}],
ParametricPlot3D[{10 Cos[u], 10 Sin[u], v}, {u, 0, 2 Pi}, {v, g[min],
g[10]}, PlotStyle -> {LightBlue, Opacity[0.5]}, Mesh -> None],
BoxRatios -> 1]


First result is volume below bounding surface, second result of volume of bounding cylinder, third is numerical value of volume below bounding surface and fourth result approximation of volume of solid of revolution 14015.9

• Use the xyzdata and int in the question. We define a ImplicitRegion bounded by that surface.
{{x1, x2}, {y1, y2}, {z1, z2}} = CoordinateBounds[xyzdata];
reg1 = ImplicitRegion[
int[x, y] <= z <= z2, {{x, x1, x2}, {y, y1, y2}, {z, z1, z2}}];
dreg1 =DiscretizeRegion[reg1, BoxRatios -> 1, MaxCellMeasure -> .1]
dreg1 // Volume


13926.6.

The Mathematica documentation actually has a relevant example that works very well and is similar to what others have been saying. This is the full code for my solution using the example data above:

xydata = Table[{x, x^2 + 8*Cos[x]}, {x, 0, 10, 0.2}];
xyzdata = Flatten[Table[{xydata[[n, 1]]*Cos[\[Theta]],xydata[[n, 1]]*Sin[\
[Theta]], xydata[[n,2]]}, {n,1,Length[xydata]}, {\[Theta], 0, 2 Pi, Pi/30}], 1];
xyzdata = DeleteDuplicates[xyzdata];

r = ConvexHullMesh[xyzdata];
volume = NIntegrate[1,{x,y,z}\[Element] r]


I tested this method out on some known volumes (paraboloids for example) and got the correct answer for volume. A few notes:

• The NIntegrate solution asymptotically approached the true value as I turned up the sample rate in the generated xyzdata.
• ConvexHullMesh was far less sensitive to duplicate points than the BoundaryRegionMesh and Interpolate functions other responses had suggested and tended to be fairly computationally fast.