Topic title says it all, really. Looking for a way to calculate the volume enclosed by a collection of polygons.

Now that is an interesting problem. I have not made a script for this yet unfortunately. Just to be clear to anyone who might attempt this, you mention convex hull in the title, but in the description you mention a collection of polygons. If the script writer assumes convexity would it work for your needs? (I have a feeling a concave shape would be much more difficult) I know in 2d the best way to get the area of a convex polygon is by breaking it down into triangles and adding the area of all the triangles. Perhaps there is a similar approach in 3d. In any case, I know that qhull (a convex hull code library) has the ability to tell you the volume. There is a python wrapper to qhull called pyhull. I would investigate this before writing your own.

Aha - thanks - I hadn't found qhull - I'll have to look into that.

My simple estimation approach was going to be this:

Find a point at the center of the volume. Then add up the volumes of the cones made by that point and each face. Volume of each cone is 1/3 × (height) × (area of face).

That would give me a rough estimate, and yes, it would have to be a convex shape for this to be even close.

But I'll go and take a look at qhull and see if that helps me out.

Your cone approach seems dead on to me! I knew there was a good 3d analogy to the 2d approach. :)

This is what I ended up with. It's a bit hacky but it actually comes out pretty damn close for most volumes. Build a convex hull and select all the polys, then run this script.

# ------------------------------------------------------------------------------------------------------------
# A script to guesstimate the enclosed volume of a convex hull.
# Process steps for each polygon:
# - find area of polygon
# - find centerpoint of polygon
# - calculate distance from the combined centerpoint of all polys in the select list, to centerpoint of polygon
# - calculate a rough cone/pyramid volume using this formula: ((height) × (area of face))/3
# Then add up all the individual volumes to give a rough estimate.
# For cubes, the volume is 100% correct.
# For spheres it comes out at 95% correct. (under-calculates by 5%)
# For complex convex hulls it comes out about 90% correct (under-calculated by about 10%)
# *** FOR CONCAVE OR OPEN SHAPES, CALCULATION IS INVALID ***
# Chris Longhurst / 2015
# ------------------------------------------------------------------------------------------------------------
# Set up hex masks for status box
MMBX_OK = 0x00000001
MMBX_OKCANCEL = 0x00000002
MMBX_YESNO = 0x00000004
MMBX_YESNOCANCEL = 0x00000008
MMBX_STATUS = 0x00000100
MMBX_WARNING = 0x00000200
MMBX_ERROR = 0x00000400
MMBX_QUESTION = 0x00000800

#Import math library
import math

# Determinant of matrix a
def det(a):
    return a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][2]*a[1][1]*a[2][0] - a[0][1]*a[1][0]*a[2][2] - a[0][0]*a[1][2]*a[2][1]

# Unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
    x = det([[1,a[1],a[2]], [1,b[1],b[2]], [1,c[1],c[2]]])
    y = det([[a[0],1,a[2]], [b[0],1,b[2]], [c[0],1,c[2]]])
    z = det([[a[0],a[1],1], [b[0],b[1],1], [c[0],c[1],1]])
    magnitude = (x**2 + y**2 + z**2)**.5
    return (x/magnitude, y/magnitude, z/magnitude)

# Dot product of vectors a and b
def dot(a, b):
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

# Cross product of vectors a and b
def cross(a, b):
    x = a[1] * b[2] - a[2] * b[1]
    y = a[2] * b[0] - a[0] * b[2]
    z = a[0] * b[1] - a[1] * b[0]
    return (x, y, z)

# Area of polygon poly
def calculate3DPolygonArea(rec):
    total = [0, 0, 0]
    numVtx=mgCountChild(rec)
    for i in range(1,numVtx+1):
        tempVtx = mgGetChildNth(rec,i)
        ok, x1, y1, z1 = mgGetVtxCoord(tempVtx)
        vtx1=[x1,y1,z1]
        if i is numVtx:
            tempVtx = mgGetChildNth(rec,1)
        else:
            tempVtx = mgGetChildNth(rec,i+1)
        ok, x2, y2, z2 = mgGetVtxCoord(tempVtx)
        vtx2=[x2,y2,z2]
        prod = cross(vtx1, vtx2)
        total[0] += prod[0]
        total[1] += prod[1]
        total[2] += prod[2]
    tempVtx1=mgGetChildNth(rec,1)
    tempVtx2=mgGetChildNth(rec,2)
    tempVtx3=mgGetChildNth(rec,3)
    ok, x1, y1, z1 = mgGetVtxCoord(tempVtx1)
    ok, x2, y2, z2 = mgGetVtxCoord(tempVtx2)
    ok, x3, y3, z3 = mgGetVtxCoord(tempVtx3)
    vtx1=[x1,y1,z1]
    vtx2=[x2,y2,z2]
    vtx3=[x3,y3,z3]
    result = dot(total, unit_normal(vtx1, vtx2, vtx3))
    return abs(result/2)

def findCenterPoint (db, parent, rec, i):
    global midX
    global midY
    global midZ
    global vertexCount
    type = mgGetCode (rec)
    if (type==fltVertex):
        ok, x1, y1, z1 = mgGetVtxCoord(rec)
        midX += x1
        midY += y1
        midZ += z1
        vertexCount+=1
    return MG_TRUE

def calculate3DPolygonMidPoint(rec):
    ok,box = mgGetBounds (rec)
    coord = mgBoxGetCenter (box)
    vtx = mgNewConstructVertex (editorContext, coord)
    ok, x, y, z = mgGetVtxCoord(vtx)
    return x,y,z

editorContext = mgNewEditorContext ("Calculate convex hull volume")
db = mgEditorGetDbRec (editorContext)

# Only run this tool if something is selected - specifically polygons
selectList = mgGetSelectList (db)
num = mgGetRecListCount (selectList)
if (num == 0):
    mgMessageDialog(None, "Convex hull volume guesstimator", "Nothing selected", MMBX_OK+MMBX_ERROR)
else:
    temprec,tempm = mgGetNextRecInList(selectList)
    tempitem = mgGetCode(temprec)
    if (tempitem != fltPolygon):
        mgMessageDialog(None, "Convex hull volume guesstimator", "No polys selected", MMBX_OK+MMBX_ERROR)
    else:
        # Reset the select list
        selectList = mgGetSelectList (db)
        totalVolume=0

        # Find a rough centerpoint and dump a construction vertex there
        # Hacky method - add all the vertex values and average them
        vertexCount=0
        midX=0
        midY=0
        midZ=0
        mgWalkRecList(selectList, findCenterPoint, None, 0, MWALK_VERTEX)
        averageX=midX/vertexCount
        averageY=midY/vertexCount
        averageZ=midZ/vertexCount
        norm = mgCoord3dZero()
        norm.x= averageX
        norm.y= averageY
        norm.z= averageZ
        tempVtx = mgNewConstructVertex (editorContext, norm)

        selectList = mgGetSelectList (db)
        num = mgGetRecListCount (selectList)
        rec,matrix = mgGetNextRecInList(selectList)
        while(rec):
            polygonArea=calculate3DPolygonArea(rec)
            polygonMidPointX,polygonMidPointY,polygonMidPointZ=calculate3DPolygonMidPoint(rec)
            tempHeight=math.sqrt( (polygonMidPointX-averageX)**2 + (polygonMidPointY-averageY)**2 + (polygonMidPointZ-averageZ)**2 )
            tempVolume = (tempHeight*polygonArea)/3
            totalVolume += tempVolume
            rec,matrix = mgGetNextRecInList(selectList)
        print "Total estimated volume: ",totalVolume
        message="Guessed volume: "+str(totalVolume)
        mgMessageDialog(None, "Convex hull volume guesstimator", message, MMBX_OK+MMBX_STATUS)

Hi Chris,

Completely unrelated to the convex hull code but I wonder why your script defines MMBX_OK, etc. Those symbols should be defined by the OpenFlight API. You should not have to re-define them in your script. What version of OpenFlight API / Script are you using?

I think I see where some of your error in volume is coming from... You are calculating the center of the polygon, then getting the distance to your center of hull from this point. This really isn't always your height unless the center of the polygon is the closest point to the center of hull. You need to find the distance to the polygon's plane from the center of hull.

An easy way to do this is use mgCoord3dProjectOnPlane to get your second point rather than finding a point in the polygon.

Re: message boxes - just force of habit. Did that ages ago and didn't realise I didn't need it any more.

Re: point on plane. Good idea - except is there an easy function to return the mgplaned of a given polygon record?

I'm trying to write this function but the plane record is always coming back as 0,0,0,0

def calculate3DPolygonMidPoint(hullCenterVtx,rec):
    name=mgGetName(rec)
    print "Processing polygon ",name
    # Create an mgcoord3d from the centerpoint vertex
    hullCenterCoord=mgCoord3dZero()
    ok, hullCenterCoord.x, hullCenterCoord.y, hullCenterCoord.z = mgGetVtxCoord(hullCenterVtx)
    
    # Create an mgvectord to hold the polygon unit normal
    polygonUnitNormal=mgVectordZero()
    
    # Get the first three verts and feed them to the unit_normal function.
    # Stuff the results in the new mgvectord
    tempVtx1=mgGetChildNth(rec,1)
    tempVtx2=mgGetChildNth(rec,2)
    tempVtx3=mgGetChildNth(rec,3)
    ok, x1, y1, z1 = mgGetVtxCoord(tempVtx1)
    ok, x2, y2, z2 = mgGetVtxCoord(tempVtx2)
    ok, x3, y3, z3 = mgGetVtxCoord(tempVtx3)
    vtx1=[x1,y1,z1]
    vtx2=[x2,y2,z2]
    vtx3=[x3,y3,z3]
    polygonUnitNormal.x,polygonUnitNormal.y,polygonUnitNormal.z=unit_normal(vtx1, vtx2, vtx3)
    
    # Make a plane given the first vertex and the unit normal
    firstVertexCoord=mgCoord3dZero()
    ok, firstVertexCoord.x, firstVertexCoord.y, firstVertexCoord.z = mgGetVtxCoord(tempVtx1)
    tempPlane = mgMakePlaned(firstVertexCoord,polygonUnitNormal)
    print tempPlane.a,tempPlane.b,tempPlane.c,tempPlane.d
    # ^^^^^^ This line is always returning zero ^^^^^^
    pointOnPlane = mgCoord3dProjectOnPlane(hullCenterCoord, tempPlane)
    vtx = mgNewConstructVertex (editorContext, pointOnPlane)
    return pointOnPlane.x,pointOnPlane.y,pointOnPlane.z

I’m trying to write this function but the plane record is always coming back as 0,0,0,0

In addition to looking at the end result, try printing out the intermediate results to make sure the normal you are calculating is correct.

Oh, and there is a function to return the polygon normal so you don't have to calculate it. I put all this together to make a function to return the plane of a polygon:

def GetPolyPlane(poly):
    normal = mgvectord()
    b,normal.i,normal.j,normal.k = mgGetPolyNormal(poly)
    vtx = mgGetChild(poly)
    coord = mgcoord3d() 
    b,coord.x,coord.y,coord.z = mgGetVtxCoord(vtx)
    plane = mgMakePlaned(coord, normal)
    return plane

Saviour.

Is it just me or is the new API documentation (the chm file) not as easy to use as the old HTML one? I would likely have found mgGetPolyNormal if I'd been able to scroll through all the functions at a glance. Now with the CHM files you kinda have to know what to look for. :(

Anyway .. modified code, tested on some weird shapes with much better results .. I'm sure I can clean this up more but it does seem to do the trick.

# ------------------------------------------------------------------------------------------------------------
# A script to guesstimate the enclosed volume of a convex hull.
# Process steps for each polygon:
# - find area of polygon
# - find centerpoint of polygon
# - calculate distance from the combined centerpoint of all polys in the select list, to centerpoint of polygon
# - calculate a rough cone/pyramid volume using this formula: ((height) × (area of face))/3
# Then add up all the individual volumes to give a rough estimate.
# For cubes, the volume is 100% correct.
# For spheres it comes out at 95% correct. (under-calculates by 5%)
# For complex convex hulls it comes out about 90% correct (under-calculated by about 10%)
# *** FOR CONCAVE OR OPEN SHAPES, CALCULATION IS INVALID ***
# Chris Longhurst / 2015
# ------------------------------------------------------------------------------------------------------------
# Status box possible values
# MMBX_OK
# MMBX_OKCANCEL
# MMBX_YESNO
# MMBX_YESNOCANCEL
# MMBX_STATUS
# MMBX_WARNING
# MMBX_ERROR
# MMBX_QUESTION

#Import math library
import math

# Determinant of matrix a
def det(a):
    return a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][2]*a[1][1]*a[2][0] - a[0][1]*a[1][0]*a[2][2] - a[0][0]*a[1][2]*a[2][1]

# Unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
    x = det([[1,a[1],a[2]], [1,b[1],b[2]], [1,c[1],c[2]]])
    y = det([[a[0],1,a[2]], [b[0],1,b[2]], [c[0],1,c[2]]])
    z = det([[a[0],a[1],1], [b[0],b[1],1], [c[0],c[1],1]])
    magnitude = (x**2 + y**2 + z**2)**.5
    return (x/magnitude, y/magnitude, z/magnitude)

# Dot product of vectors a and b
def dot(a, b):
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

# Cross product of vectors a and b
def cross(a, b):
    x = a[1] * b[2] - a[2] * b[1]
    y = a[2] * b[0] - a[0] * b[2]
    z = a[0] * b[1] - a[1] * b[0]
    return (x, y, z)

# Area of polygon poly
def calculate3DPolygonArea(rec):
    total = [0, 0, 0]
    numVtx=mgCountChild(rec)
    for i in range(1,numVtx+1):
        tempVtx = mgGetChildNth(rec,i)
        ok, x1, y1, z1 = mgGetVtxCoord(tempVtx)
        vtx1=[x1,y1,z1]
        if i is numVtx:
            tempVtx = mgGetChildNth(rec,1)
        else:
            tempVtx = mgGetChildNth(rec,i+1)
        ok, x2, y2, z2 = mgGetVtxCoord(tempVtx)
        vtx2=[x2,y2,z2]
        prod = cross(vtx1, vtx2)
        total[0] += prod[0]
        total[1] += prod[1]
        total[2] += prod[2]
    tempVtx1=mgGetChildNth(rec,1)
    tempVtx2=mgGetChildNth(rec,2)
    tempVtx3=mgGetChildNth(rec,3)
    ok, x1, y1, z1 = mgGetVtxCoord(tempVtx1)
    ok, x2, y2, z2 = mgGetVtxCoord(tempVtx2)
    ok, x3, y3, z3 = mgGetVtxCoord(tempVtx3)
    vtx1=[x1,y1,z1]
    vtx2=[x2,y2,z2]
    vtx3=[x3,y3,z3]
    result = dot(total, unit_normal(vtx1, vtx2, vtx3))
    return abs(result/2)

def findCenterPoint (db, parent, rec, i):
    global midX
    global midY
    global midZ
    global vertexCount
    type = mgGetCode (rec)
    if (type==fltVertex):
        ok, x1, y1, z1 = mgGetVtxCoord(rec)
        midX += x1
        midY += y1
        midZ += z1
        vertexCount+=1
    return MG_TRUE

def GetPolyPlane(poly):
    normal = mgvectord()
    b,normal.i,normal.j,normal.k = mgGetPolyNormal(poly)
    vtx = mgGetChild(poly)
    coord = mgcoord3d() 
    b,coord.x,coord.y,coord.z = mgGetVtxCoord(vtx)
    plane = mgMakePlaned(coord, normal)
    return plane

def calculate3DPolygonMidPoint(hullCenterVtx,rec):
    # Create an mgcoord3d from the centerpoint vertex
    hullCenterCoord=mgCoord3dZero()
    ok, hullCenterCoord.x, hullCenterCoord.y, hullCenterCoord.z = mgGetVtxCoord(hullCenterVtx)
    tempPlane = GetPolyPlane(rec)
    pointOnPlane = mgCoord3dProjectOnPlane(hullCenterCoord, tempPlane)
    vtx = mgNewConstructVertex (editorContext, pointOnPlane)
    return pointOnPlane.x,pointOnPlane.y,pointOnPlane.z

editorContext = mgNewEditorContext ("Calculate convex hull volume")
db = mgEditorGetDbRec (editorContext)

# Only run this tool if something is selected - specifically polygons
selectList = mgGetSelectList (db)
num = mgGetRecListCount (selectList)
if (num == 0):
    mgMessageDialog(None, "Convex hull volume guesstimator", "Nothing selected", MMBX_OK+MMBX_ERROR)
else:
    temprec,tempm = mgGetNextRecInList(selectList)
    tempitem = mgGetCode(temprec)
    if (tempitem != fltPolygon):
        mgMessageDialog(None, "Convex hull volume guesstimator", "No polys selected", MMBX_OK+MMBX_ERROR)
    else:
        # Reset the select list
        selectList = mgGetSelectList (db)
        totalVolume=0

        # Find a rough centerpoint and dump a construction vertex there
        # Hacky method - add all the vertex values and average them
        vertexCount=0
        midX=0
        midY=0
        midZ=0
        mgWalkRecList(selectList, findCenterPoint, None, 0, MWALK_VERTEX)
        averageX=midX/vertexCount
        averageY=midY/vertexCount
        averageZ=midZ/vertexCount
        norm = mgCoord3dZero()
        norm.x= averageX
        norm.y= averageY
        norm.z= averageZ
        hullCenterVtx = mgNewConstructVertex (editorContext, norm)

        selectList = mgGetSelectList (db)
        num = mgGetRecListCount (selectList)
        rec,matrix = mgGetNextRecInList(selectList)
        while(rec):
            polygonArea=calculate3DPolygonArea(rec)
            polygonMidPointX,polygonMidPointY,polygonMidPointZ=calculate3DPolygonMidPoint(hullCenterVtx,rec)
            #print polygonMidPointX,polygonMidPointY,polygonMidPointZ
            tempHeight=math.sqrt( (polygonMidPointX-averageX)**2 + (polygonMidPointY-averageY)**2 + (polygonMidPointZ-averageZ)**2 )
            tempVolume = (tempHeight*polygonArea)/3
            totalVolume += tempVolume
            rec,matrix = mgGetNextRecInList(selectList)
        print "Total estimated volume: ",totalVolume
        message="Guessed volume: "+str(totalVolume)
        mgMessageDialog(None, "Convex hull volume guesstimator", message, MMBX_OK+MMBX_STATUS)

Is it just me or is the new API documentation (the chm file) not as easy to use as the old HTML one? I would likely have found mgGetPolyNormal if I’d been able to scroll through all the functions at a glance. Now with the CHM files you kinda have to know what to look for. :(

The HTML help files are still included. No need to use the CHM if you find them less usable.

The HTML help files are still included

I usually just use the script editor's "OpenFlight API Reference Set" menu item from the help menu to get directly to the html versions. Or the "Help On..." to get help about the current function. Keeps my head from exploding trying to memorize too much.

Is it just me or is the new API documentation (the chm file) not as easy to use as the old HTML one? I would likely have found mgGetPolyNormal if I’d been able to scroll through all the functions at a glance. Now with the CHM files you kinda have to know what to look for. :(

The HTML help files are still included. No need to use the CHM if you find them less usable.

Hmm. When I did the Openflight API 13 installation I got CHM and PDF files but no html. Was installed here : C:\Presagis\Openflight_API_13

Maybe I need to re-install - I might have missed something or checked a box somewhere along the line.

Anyway - thanks for the help on the plane equation thing - this gives me a good solution for the problem I was trying to solve now :D

The HTML help files are still included

I usually just use the script editor's "OpenFlight API Reference Set" menu item from the help menu to get directly to the html versions. Or the "Help On..." to get help about the current function. Keeps my head from exploding trying to memorize too much.

That option is greyed out in my script editor help menu (Openflight API Reference Set)

That option is greyed out in my script editor help menu (Openflight API Reference Set)

Hmmm...

Obviously you have the OpenFlight API installed. As you would expect, context help would not work if that was not the case. What version of Creator are you using now? Different versions use different strategies for "locating" the API docs:

Creator 14 and earlier use the env variable PRESAGIS_OPENFLIGHT_API to locate:

$PRESAGIS_OPENFLIGHT_API/docs/reference/OpenFlight_API_Reference_Set.htm.

This can be problematic if for some reason PRESAGIS_OPENFLIGHT_API is not set or set to the wrong version of the API

Creator 15 got a bit smarter and looks adjacent to its own installation folder for the OpenFlight API installation and then dives into the docs folder from there. This is less likely to get confused by missing environment variables and won't end up linking to a version of the API that might not match "this" version of Creator.

In any event, check the Creator status log after startup. There might be some diagnostic information there about why the API help could not be linked to properly.

Whatever version you are using the "links" in the Creator OpenFlight Script Editor are worth trying to fix. They will make your life much easier. You'll be able to, for example, move the text caret over the name of an API function, do Help>On and go straight to the help for that function.