3D volume calculation

I have a 3D point cloud and would like to find its volume. My strategy has been to perform a ConvexHull on my tripletwave (x,y,z coords). I have M_Hull and can visualise it in the Gizmo, but I am stuck as to how to get out a volume for it. Any help appreciated. Maybe there is a much simpler way to do this, that I am missing?

Igor

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