Start a new topic

## Has anyone written a script that can calculate the enclosed volume of a convex hull?

Original Post by: chrisell Wed Dec 2 20:38:03 2015

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

Original Post by: ChrisRogers Thu Dec 3 15:44:19 2015

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.

Original Post by: chrisell Thu Dec 3 15:49:43 2015

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.

Original Post by: ChrisRogers Thu Dec 3 16:12:05 2015

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

Original Post by: chrisell Thu Dec 3 18:01:57 2015

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 boxMMBX_OK = 0x00000001MMBX_OKCANCEL = 0x00000002MMBX_YESNO = 0x00000004MMBX_YESNOCANCEL = 0x00000008MMBX_STATUS = 0x00000100MMBX_WARNING = 0x00000200MMBX_ERROR = 0x00000400MMBX_QUESTION = 0x00000800#Import math libraryimport math# Determinant of matrix adef det(a):    return a*a*a + a*a*a + a*a*a - a*a*a - a*a*a - a*a*a# Unit normal vector of plane defined by points a, b, and cdef unit_normal(a, b, c):    x = det([[1,a,a], [1,b,b], [1,c,c]])    y = det([[a,1,a], [b,1,b], [c,1,c]])    z = det([[a,a,1], [b,b,1], [c,c,1]])    magnitude = (x**2 + y**2 + z**2)**.5    return (x/magnitude, y/magnitude, z/magnitude)# Dot product of vectors a and bdef dot(a, b):    return a*b + a*b + a*b# Cross product of vectors a and bdef cross(a, b):    x = a * b - a * b    y = a * b - a * b    z = a * b - a * b    return (x, y, z)# Area of polygon polydef 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 += prod        total += prod        total += prod    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_TRUEdef calculate3DPolygonMidPoint(rec):    ok,box = mgGetBounds (rec)    coord = mgBoxGetCenter (box)    vtx = mgNewConstructVertex (editorContext, coord)    ok, x, y, z = mgGetVtxCoord(vtx)    return x,y,zeditorContext = mgNewEditorContext ("Calculate convex hull volume")db = mgEditorGetDbRec (editorContext)# Only run this tool if something is selected - specifically polygonsselectList = 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)`
Original Post by: SteveThompson Thu Dec 3 19:39:54 2015

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?

Original Post by: ChrisRogers Thu Dec 3 20:29:56 2015

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.

Original Post by: chrisell Thu Dec 3 21:54:14 2015

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`
Original Post by: SteveThompson Thu Dec 3 22:35:27 2015

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`
Original Post by: chrisell Thu Dec 3 22:57:04 2015

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 libraryimport math# Determinant of matrix adef det(a):    return a*a*a + a*a*a + a*a*a - a*a*a - a*a*a - a*a*a# Unit normal vector of plane defined by points a, b, and cdef unit_normal(a, b, c):    x = det([[1,a,a], [1,b,b], [1,c,c]])    y = det([[a,1,a], [b,1,b], [c,1,c]])    z = det([[a,a,1], [b,b,1], [c,c,1]])    magnitude = (x**2 + y**2 + z**2)**.5    return (x/magnitude, y/magnitude, z/magnitude)# Dot product of vectors a and bdef dot(a, b):    return a*b + a*b + a*b# Cross product of vectors a and bdef cross(a, b):    x = a * b - a * b    y = a * b - a * b    z = a * b - a * b    return (x, y, z)# Area of polygon polydef 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 += prod        total += prod        total += prod    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_TRUEdef 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 planedef 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.zeditorContext = mgNewEditorContext ("Calculate convex hull volume")db = mgEditorGetDbRec (editorContext)# Only run this tool if something is selected - specifically polygonsselectList = 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)`
Original Post by: SteveThompson Thu Dec 3 23:01:50 2015

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.

Original Post by: ChrisRogers Thu Dec 3 23:05:45 2015

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.

Original Post by: chrisell Thu Dec 3 23:07:50 2015

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

Original Post by: chrisell Thu Dec 3 23:09:08 2015

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)

Original Post by: SteveThompson Thu Dec 3 23:51:50 2015

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. 