nurbsubdiv.c

[Game.Programming].Academic - Graphics Gems (6 books source code)

/*
 * NurbSubdiv.c - Perform adaptive subdivision on a NURB surface.
 *
 * John Peterson
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "nurbs.h"
#include "drawing.h"

#define EPSILON 0.0000001   /* Used to determine when things are too small. */
#define DIVPT( p, dn ) { ((p).x) /= (dn); ((p).y) /= (dn); ((p).z) /= (dn); }

double SubdivTolerance;	    /* Size, in pixels of facets produced */

#define maxV(surf) ((surf)->numV-1L)
#define maxU(surf) ((surf)->numU-1L)

/*
 * Split a knot vector at the center, by adding multiplicity k knots near
 * the middle of the parameter range.  Tries to start with an existing knot,
 * but will add a new knot value if there's nothing in "the middle" (e.g.,
 * a Bezier curve).
 */
static long
SplitKV( double * srckv,
	 double ** destkv,
	 long * splitPt,    /* Where the knot interval is split */
	 long m, long k )
{
    long i, last;
    long middex, extra, same;	/* "middex" ==> "middle index" */
    double midVal;

    extra = 0L;
    last = (m + k);

    middex = last / 2;
    midVal = srckv[middex];

    /* Search forward and backward to see if multiple knot is already there */

    i = middex+1L;
    same = 1L;
    while ((i < last) && (srckv[i] == midVal)) {
	i++;
	same++;
    }

    i = middex-1L;
    while ((i > 0L) && (srckv[i] == midVal)) {
	i--;
	middex--;	    /* middex is start of multiple knot */
	same++;
    }

    if (i <= 0L)	    /* No knot in middle, must create it */
    {
	midVal = (srckv[0L] + srckv[last]) / 2.0;
	middex = last / 2L;
	while (srckv[middex + 1L] < midVal)
	    middex++;
	same = 0L;
    }

    extra = k - same;
    CHECK( *destkv = (double *) malloc( (long) (sizeof( double )
					* (m+k+extra+1L) ) ) );

    if (same < k)	    /* Must add knots */
    {
	for (i = 0L; i <= middex; i++)
	    (*destkv)[i] = srckv[i];

	for (i = middex+1L; i <= middex+extra; i++)
	    (*destkv)[i] = midVal;

	for (i = middex + k - same + 1L; i <= m + k + extra; i++)
	    (*destkv)[i] = srckv[i - extra];
    }
    else
    {
	for (i = 0L; i <= m + k; i++)
	(*destkv)[i] = srckv[i];
    }

    *splitPt = (extra < k) ? middex - 1L : middex;
    return( extra );
}

/*
 * Given a line defined by firstPt and lastPt, project midPt onto
 * that line.  Used for fixing "cracks".
 */
static void
ProjectToLine( Point3 * firstPt, Point3 * lastPt, Point3 * midPt )
{
    Point3 base, v0, vm;
    double fraction, denom;

    base = *firstPt;

    (void) V3Sub( lastPt, &base, &v0 );
    (void) V3Sub( midPt, &base, &vm );

    denom = V3SquaredLength( &v0 );
    fraction = (denom == 0.0) ? 0.0 : (V3Dot( &v0, &vm ) / denom);

    midPt->x = base.x + fraction * v0.x;
    midPt->y = base.y + fraction * v0.y;
    midPt->z = base.z + fraction * v0.z;
}

/*
 * If a normal has collapsed to zero (normLen == 0.0) then try
 * and fix it by looking at its neighbors.  If all the neighbors
 * are sick, then re-compute them from the plane they form.
 * If that fails too, then we give up...
 */
static void
FixNormals( SurfSample * s0, SurfSample * s1, SurfSample * s2 )
{
    Boolean goodnorm;
    long i, j, ok;
    double dist;
    SurfSample * V[3];
    Point3 norm;

    V[0] = s0; V[1] = s1; V[2] = s2;

    /* Find a reasonable normal */
    for (ok = 0, goodnorm = FALSE;
	(ok < 3L) && !(goodnorm = (V[ok]->normLen > 0.0)); ok++);

    if (! goodnorm)	/* All provided normals are zilch, try and invent one */
    {
	norm.x = 0.0; norm.y = 0.0; norm.z = 0.0;

	for (i = 0; i < 3L; i++)
	{
	    j = (i + 1L) % 3L;
	    norm.x += (V[i]->point.y - V[j]->point.y) * (V[i]->point.z + V[j]->point.z);
	    norm.y += (V[i]->point.z - V[j]->point.z) * (V[i]->point.x + V[j]->point.x);
	    norm.z += (V[i]->point.x - V[j]->point.x) * (V[i]->point.y + V[j]->point.y);
	}
	dist = V3Length( &norm );
	if (dist == 0.0)
	    return;		/* This sucker's hopeless... */

	DIVPT( norm, dist );

	for (i = 0; i < 3; i++)
	{
	    V[i]->normal = norm;
	    V[i]->normLen = dist;
	}
    }
    else	    /* Replace a sick normal with a healthy one nearby */
    {
	for (i = 0; i < 3; i++)
	    if ((i != ok) && (V[i]->normLen == 0.0))
		V[i]->normal = V[ok]->normal;
    }
    return;
}

/*
 * Normalize the normal in a sample.  If it's degenerate,
 * flag it as such by setting the normLen to 0.0
 */
static void
AdjustNormal( SurfSample * samp )
{
    /* If it's not degenerate, do the normalization now */
    samp->normLen = V3Length( &(samp->normal) );

    if (samp->normLen < EPSILON)
	samp->normLen = 0.0;
    else
	DIVPT( (samp->normal), samp->normLen );
}

/*
 * Compute the normal of a corner point of a mesh.  The
 * base is the value of the point at the corner, indU and indV
 * are the mesh indices of that point (either 0 or numU|numV).
 */
static void
GetNormal( NurbSurface * n, long indV, long indU )
{
    Point3 tmpL, tmpR;	/* "Left" and "Right" of the base point */
    SurfSample * crnr;

    if ( (indU && indV) || ((! indU) && (!indV)) )
    {
	if (indU)
	    crnr = &(n->cnn);
	else
	    crnr = &(n->c00);
	DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), &tmpL );
	DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), &tmpR );
    }
    else
    {
	if (indU)
	    crnr = &(n->c0n);
	else
	    crnr = &(n->cn0);
	DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), &tmpR );
	DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), &tmpL );
    }

    (void) V3Sub( &tmpL, &(crnr->point), &tmpL );
    (void) V3Sub( &tmpR, &(crnr->point), &tmpR );
    (void) V3Cross( &tmpL, &tmpR, &(crnr->normal) );
    AdjustNormal( crnr );
}

/*
 * Build the new corners in the two new surfaces, computing both
 * point on the surface along with the normal.	Prevent cracks that may occur.
 */
static void
MakeNewCorners( NurbSurface * parent,
		NurbSurface * kid0,
		NurbSurface * kid1,
		Boolean dirflag )
{
    DIVW( &(kid0->points[maxV(kid0)][maxU(kid0)]), &(kid0->cnn.point) );
    GetNormal( kid0, maxV(kid0), maxU(kid0) );

    if (dirflag)
    {
	kid0->strUn = FALSE;	/* Must re-test new edge straightness */

	DIVW( &(kid0->points[0L][maxU(kid0)]), &(kid0->c0n.point) );
	GetNormal( kid0, 0L, maxU(kid0) );
	/*
	 * Normals must be re-calculated for kid1 in case the surface
	 * was split at a c1 (or c0!) discontinutiy
	 */
	kid1->c00.point = kid0->c0n.point;
	GetNormal( kid1, 0L, 0L );
	kid1->cn0.point = kid0->cnn.point;
	GetNormal( kid1, maxV(kid1), 0L );

	/*
	 * Prevent cracks from forming by forcing the points on the seam to
	 * lie along any straight edges.  (Must do this BEFORE finding normals)
	 */
	if (parent->strV0)
	    ProjectToLine( &(parent->c00.point),
			   &(parent->c0n.point),
			   &(kid0->c0n.point) );
	if (parent->strVn)
	    ProjectToLine( &(parent->cn0.point),
			   &(parent->cnn.point),
			   &(kid0->cnn.point) );

	kid1->c00.point = kid0->c0n.point;
	kid1->cn0.point = kid0->cnn.point;
	kid1->strU0 = FALSE;
    }
    else
    {
	kid0->strVn = FALSE;

	DIVW( &(kid0->points[maxV(kid0)][0]), &(kid0->cn0.point) );
	GetNormal( kid0, maxV(kid0), 0L );
	kid1->c00.point = kid0->cn0.point;
	GetNormal( kid1, 0L, 0L );
	kid1->c0n.point = kid0->cnn.point;
	GetNormal( kid1, 0L, maxU(kid1) );

	if (parent->strU0)
	    ProjectToLine( &(parent->c00.point),
			   &(parent->cn0.point),
			   &(kid0->cn0.point) );
	if (parent->strUn)
	    ProjectToLine( &(parent->c0n.point),
			   &(parent->cnn.point),
			   &(kid0->cnn.point) );

	kid1->c00.point = kid0->cn0.point;
	kid1->c0n.point = kid0->cnn.point;
	kid1->strV0 = FALSE;
    }
}

/*
 * Split a surface into two halves.  First inserts multiplicity k knots
 * in the center of the parametric range.  After refinement, the two
 * resulting surfaces are copied into separate data structures.	 If the
 * parent surface had straight edges, the points of the children are
 * projected onto those edges.
 */
static void
SplitSurface( NurbSurface * parent,
	      NurbSurface * kid0, NurbSurface * kid1,
	      Boolean dirflag )	    /* If true subdivided in U, else in V */
{
    NurbSurface tmp;
    double * newkv;
    long i, j, splitPt;

    /*
     * Add a multiplicty k knot to the knot vector in the direction
     * specified by dirflag, and refine the surface.  This creates two
     * adjacent surfaces with c0 discontinuity at the seam.
     */

    tmp = *parent;	/* Copy order, # of points, etc. */
    if (dirflag)
    {
	tmp.numU = parent->numU + SplitKV( parent->kvU,
					    &newkv,
					    &splitPt,
					    maxU(parent),