Inside/Outside Surface Determination

[One poster writes:]

The user picks a "Seed polygon" and I recursively "grow" a surface starting at that polygon, a neighboring (child) polygon is considered to be flipped if the child traverses the two shared vertices in the same direction as the parent.

[Another poster replies:]

If you are not able to interactively pick such a seed polygon, there is an automatic way, but it takes time to write. The idea is an extension of a point in polygon algorithm. Basically, to find if a point is within a polygon, you can extend a ray from that point to infinity, if the ray crosses an even number of edges of the polygon, the point is outside (where 0 is defined to be even). If it crosses an odd number, it is inside.

A better method:

Consistently orient all of the polygons (you need to know which edges are shared to do this, and this isn't possible for some objects -- e.g. mobius strips, Klein bottles, etc).

Now, choose the point 'P' in the middle of the bounding box of the object. For each triangle in the object, compute a signed volume for the tetrahedron formed by the triangle and P. Arrange your calculation so that the area will be positive if P is left-hand-side of the triangle and negative if P is on the right-hand-side of the triangle. (if your triangle is ABC, then doing (AB.cross.BC).dot.P will have this property).

Add up all of the volumes. If the result is positive, then the normals are oriented 'outside'. If the result is negative, then the normals are oriented 'inside'. If the result is zero or very close to zero, then the object is flat or has just as many concave parts as convex parts.

This will always work for completely enclosed objects, and does the right thing for surfaces-- it chooses the orientation that marks the surface 'most convex'. It works for self-intersecting objects.

Here's the code I use: (If you have Inventor 1, this is in the 'ivnorm' code:)

    int i, j;

    int total_v = 0;
    SbVec3f average(0.0, 0.0, 0.0);

    for (j = 0; j < length(); j++)
    {
        Face *f = (*this)[j];
        if (f->degenerate) continue;

        for (i = 0; i < f->nv; i++)
        {
            average += verts[f->v[i]];
            ++total_v;
        }
    }
    average /= (float) total_v;

    float result = 0.0;

    for (j = 0; j < length(); j++)
    {
        Face *f = (*this)[j];
        if (f->degenerate) continue;

        for (i = 1; i < f->nv-1; i++)
        {
            SbVec3f v1 = verts[f->v[0]] - average;
            SbVec3f v2 = verts[f->v[i]] - average;
            SbVec3f v3 = verts[f->v[i+1]] - average;
            
            float t = (v1.cross(v2)).dot(v3);
            if (f->orientation == Face::CCW)
            {
                result += t;
            }
            else if (f->orientation == Face::CW)
            {
                result -= t;
            }
            else
            {
                assert(0);
            }
        }
    }
    return result > 0;