3D volume calculation

Using the convex hull makes sense but you should do so via the operation Triangulate3D. Here is an example:

make/n=(33,3) ddd=enoise(1)  // create sample data

triangulate3d/out=2 ddd  // perform the triangulation and get tetrahedra output

To plot this data use:

Window gizmo0() : GizmoPlot

    PauseUpdate; Silent 1       // building window...

    // Building Gizmo 7 window...

    NewGizmo/W=(35,45,550,505)

    ModifyGizmo startRecMacro=700

    ModifyGizmo scalingOption=63

    AppendToGizmo Scatter=root:ddd,name=scatter0

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ scatterColorType,0}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ markerType,0}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ sizeType,0}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ rotationType,0}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ Shape,2}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ size,0.25}

    ModifyGizmo ModifyObject=scatter0,objectType=scatter,property={ color,1,0,0,1}

    AppendToGizmo Axes=boxAxes,name=axes0

    ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}

    ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisColor,0,0,0,1}

    ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3}

    ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3}

    ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3}

    ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}

    AppendToGizmo Path=root:M_TetraPath,name=path0

    ModifyGizmo ModifyObject=path0,objectType=path,property={ pathColorType,1}

    ModifyGizmo ModifyObject=path0,objectType=path,property={ lineWidthType,0}

    ModifyGizmo ModifyObject=path0,objectType=path,property={ pathColor,0,0,1,1}

    ModifyGizmo setDisplayList=0, object=scatter0

    ModifyGizmo setDisplayList=1, object=axes0

    ModifyGizmo setDisplayList=2, object=path0

    ModifyGizmo autoscaling=1

    ModifyGizmo currentGroupObject=""

    ModifyGizmo showInfo

    ModifyGizmo infoWindow={551,23,1368,320}

    ModifyGizmo zoomFactor=1.180000

    ModifyGizmo endRecMacro

    ModifyGizmo SETQUATERNION={-0.497650,0.659761,0.323963,-0.460555}

EndMacro

If you inspect the wave M_TetraPath you will see that it contains triangle vertices separated by NaN rows. In each triangle the last vertex repeats the first in order to close the path for drawing. If you collect every 4 consecutive triangles you have a full tetrahedron.

Given the list of tetrahedra, you can compute the volume by summing the individual volumes of all the tetrahedra because by definition, the triangulation decomposes the volume into non-overlapping tetrahedra. One way to compute the volume of an arbitrary tetrahedron is using Igor's function MatrixDet(waveT)/6 where waveT is a 2D wave constructed from the vertices of the tetrahedron as follows:

Make/O/N=(4,4) waveT={{x1,x2,x3,x4},{y1,y2,y3,y4},{z1,z2,z3,z4},{1,1,1,1}}

for a matrix of the form:
x1,y1,z1,1
x2,y2,z2,1
x3,y3,z3,1
x4,y4,z4,1

I hope this helps,

A.G.
WaveMetrics, Inc.

Note: in IP7 (after the nightly build) you can get the result directly using the commands:

 Triangulate3D/VOL ddd

Print V_value

February 15, 2016 at 10:04 am - Permalink

sjr51

Thanks A.G.

The new flag in the nightly build works great.

Edit: deleted my comment about timing since this referred to my own procedure and isn't of interest to anyone else.

February 16, 2016 at 11:15 am - Permalink