AreaVolume script?

Lightwave has a plug-in called AreaVolume. It takes the selected polygons in the foreground layer and calculates their surface area and volume.
http://amber.rc.arizona.edu/lw/areavol.html

It is based on Ronald N. Goldman, “Area of Planar Polygons and Volume of Polyhedra,” in Graphics Gems II, which was suitably intimidating to me.

However, it did include the source code. Examining it, it seems like it would be easy enough to convert to Python. But I have been utterly mistaken before.


/*
======================================================================
areavol.c

Ernie Wright  29 Nov 96
SAS/C 6.51
MSVC 4.0

LightWave Modeler plug-in that calculates the surface area and volume
of the selected polygons in the foreground layer.  Based on Ronald N.
Goldman, "Area of Planar Polygons and Volume of Polyhedra," Graphics
Gems II, James Arvo ed., Academic Press, 1991 (ISBN 0-12-064480-0).
====================================================================== */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <splug.h>
#include <lwmod.h>


typedef struct { double x, y, z; } VEC3;

static EDError polycb( MeshEditOp *edit, const PolygonInfo *polygon );
static double vec_dot( VEC3 *a, VEC3 *b );
static VEC3 *vec_cross( VEC3 *a, VEC3 *b, VEC3 *c );

double volume, area;
int selection;


/*
======================================================================
AreaVolume()

Activation function.

INPUTS
   version     should be 1
   global      access to global functions
   local       a ModCommand with access to mesh edit functions
   serverData  not used

RESULTS
   Displays the area and volume of the selected polygons.

All of the calculation is done during the edit->polyScan() call.  The
number of selected polygons is obtained first.  If this is zero, the
callback operates on all polygons in the foreground layer, otherwise
it ignores polygons that weren't explicitly selected.
====================================================================== */

XCALL_( int )
AreaVolume( long version, GlobalFunc *global, ModCommand *local,
   void *serverData)
{
   MeshEditOp *edit;
   MessageFuncs *message;
   char buf[ 2 ][ 40 ];


   XCALL_INIT;
   if ( version != 1 ) return AFUNC_BADVERSION;

   volume = area = 0;

   edit = local->editBegin( 0, 0, OPSEL_USER );
   selection = edit->polyCount( edit->state, OPLYR_FG, EDCOUNT_SELECT );
   edit->polyScan( edit->state, polycb, edit, OPLYR_FG );
   edit->done( edit->state, EDERR_NONE, 0 );

   sprintf( buf[ 0 ], "Area:  %lg square meters", area );
   sprintf( buf[ 1 ], "Volume:  %lg cubic meters", fabs( volume ));

   if ( message = global( "Info Messages", GFUSE_TRANSIENT ))
      message->info( buf[ 0 ], buf[ 1 ] );

   return AFUNC_OK;
}


/*
======================================================================
polycb()

Polygon enumeration callback.

INPUTS
   edit        MeshEditOp with access to edit and query functions
   polygon     the current polygon

RESULTS
   Adds the contribution of the current polygon to the calculation of
   area and volume.

The methods presented in Goldman involve sums over all polygons, so
as each polygon is presented to polycb(), the relevant quantities are
simply added to the two global variables for area and volume.

This function ignores curves, detail polygons, and one- and two-point
polygons.  If polygons were explicitly selected by the user, it also
ignores unselected polygons.
====================================================================== */

static EDError polycb( MeshEditOp *edit, const PolygonInfo *polygon )
{
   PointInfo *pt;
   VEC3 vert[ 2 ], v0, norm, v = { 0 };
   double a;
   int i, j;


   if (( polygon->numPnts < 3 ) || ( polygon->flags & PPDF_CURVE ) ||
       ( polygon->flags & PPDF_DETAIL )) return EDERR_NONE;

   if ( selection && !( polygon->flags & PPDF_SELECT ))
      return EDERR_NONE;

   pt = edit->pointInfo( edit->state, polygon->points[ 0 ] );
   memcpy( &vert[ 0 ], pt->position, sizeof( VEC3 ));

   for ( i = 0; i < polygon->numPnts; i++ ) {
      j = ( i + 1 ) % polygon->numPnts;
      pt = edit->pointInfo( edit->state, polygon->points[ j ] );
      memcpy( &vert[ 1 ], pt->position, sizeof( VEC3 ));
      vec_cross( &vert[ 0 ], &vert[ 1 ], &v0 );
      v.x += v0.x;
      v.y += v0.y;
      v.z += v0.z;
      vert[ 0 ] = vert[ 1 ];
   }

   edit->polyNormal( edit->state, polygon->pol, ( double * ) &norm );
   area += a = .5 * fabs( vec_dot( &norm, &v ));
   volume += vec_dot( &vert[ 0 ], &norm ) * a / 3;

   return EDERR_NONE;
}


/*
======================================================================
vec_dot()
vec_cross()

Dot and cross products.
====================================================================== */

static double vec_dot( VEC3 *a, VEC3 *b )
{
   return a->x * b->x + a->y * b->y + a->z * b->z;
}


static VEC3 *vec_cross( VEC3 *a, VEC3 *b, VEC3 *c )
{
   c->x = a->y * b->z - a->z * b->y;
   c->y = a->z * b->x - a->x * b->z;
   c->z = a->x * b->y - a->y * b->x;
   return c;
}


ServerRecord ServerDesc[] = {
   { "CommandSequence", "AreaVolume", AreaVolume },
   { NULL }
};

I know there is a branch of Blender that includes such a calculation, does it use a similar algorithm?

area and angle drawing is actually in bf-blender and will be in the next release

calculating the volume of closed manifold [solid] meshes is pretty trivial, but I forget the exact method now [it was summing the volume of the shape bounded by 3 verts [a triangular face] and the origin, negative if the front side of the face faces the origin]

That shouldn’t be hard to code, but it will be slow for large models. I know because I wrote a wings plugin to calculate surface area and it is slow.

scorpius,

surface area calculations in the flatten script by cambo is pretty fast as far as i can tell, you might have a look at his calculation method.

LetterRip

Are we talking tri - Area here?

I’m looking at it now and it is using the same ‘Cross Product’ method. The question is if the Blender.Mathutils.CrossVecs() function is faster than the Python equivalent.

Erlang and Python are similar in terms of speed. Both can call external libs, so if ‘Mathutils’ is written in C, then there might be a speed difference.

scorpius,

mathutils is written in C.

LetterRip

Then it is faster. In fact, 2x faster. The following calculates the surface area of triangulated mesh objects:

import Blender, vector, time

object = Blender.Object.GetSelected()
objname = object[0].name
meshname = object[0].data.name
mesh = Blender.NMesh.GetRaw(meshname)

def tri_area1(face):
    v1, v2, v3 = face.v
    v1, v2, v3 = v1.co, v2.co, v3.co
    e1 = v1-v2
    e2 = v3-v2
    cp = Blender.Mathutils.CrossVecs(e1,e2)
    area = cp.length/2
    return area

def tri_area2(face):
    v1, v2, v3 = face.v
    v1, v2, v3 = v1.co, v2.co, v3.co
    e1 = vector.vsub(v1,v2)
    e2 = vector.vsub(v3,v2)
    cp = vector.vcross(e1, e2)
    area = vector.vlen(cp)/2
    #volume = vector.vdot(v1, cp) * area / 3
    return area

print
start = time.clock()
print "Area1:", sum([tri_area1(face) for face in mesh.faces]),
end = time.clock()
print "in %.2f %s" % (end-start, "seconds")

start = time.clock()
print "Area2:", sum([tri_area2(face) for face in mesh.faces]),
end = time.clock()
print "in %.2f %s" % (end-start, "seconds")

Here’s a script that calculates the area and volume for closed & triangulated polyhedrons only. Try it on a cube and see that the area and volume should be exacly 24 and 8, but they are not! That is, until you select and snap all verts to the grid. Maybe blender’s cube object has some floating point error.

I also tried it with my geodome script and, as the subdivision level is increased, the area and volume approach 4pi and 4/3pi, respectively.

The only thing that remains is to calculate correctly for a torus type object, where some faces point to the origin.

import Blender, time  

object = Blender.Object.GetSelected()
objname = object[0].name
meshname = object[0].data.name
mesh = Blender.NMesh.GetRaw(meshname)

def area_vol(face):
	v1, v2, v3 = face.v
	v1, v2, v3 = v1.co, v2.co, v3.co
	e1 = v1-v2
	e2 = v3-v2
	cp = Blender.Mathutils.CrossVecs(e1,e2)
	area = cp.length/2
	bc = Blender.Mathutils.CrossVecs(v2,v3)
	volume = Blender.Mathutils.DotVecs(v1, bc) / 6	  # Vol = (1/6)|A.(B x C)|
	return area, volume

print
start = time.clock()
both = [area_vol(face) for face in mesh.faces]
print "Area:",   sum([a for a,v in both]), "square Blender units"
print "Volume:", sum([v for a,v in both]),  "cubic Blender units"
end = time.clock()
print "in %.2f %s" % (end-start, "seconds"), meshname