nurbsubdiv.c

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

					    parent->orderU );
	AllocNurb( &tmp, newkv, NULL );
	for (i = 0L; i < tmp.numV + tmp.orderV; i++)
	    tmp.kvV[i] = parent->kvV[i];
    }
    else
    {
	tmp.numV = parent->numV + SplitKV( parent->kvV,
					    &newkv,
					    &splitPt,
					    maxV(parent),
					    parent->orderV );
	AllocNurb( &tmp, NULL, newkv );
	for (i = 0L; i < tmp.numU + tmp.orderU; i++)
	    tmp.kvU[i] = parent->kvU[i];
    }
    RefineSurface( parent, &tmp, dirflag );

    /*
     * Build the two child surfaces, and copy the data from the refined
     * version of the parent (tmp) into the two children
     */

    /* First half */

    *kid0 = *parent;	/* copy various edge flags and orders */

    kid0->numU = dirflag ? splitPt+1L : parent->numU;
    kid0->numV = dirflag ? parent->numV : splitPt+1L;
    kid0->kvU = kid0->kvV = NULL;
    kid0->points = NULL;
    AllocNurb( kid0, NULL, NULL );

    for (i = 0L; i < kid0->numV; i++)	/* Copy the point and kv data */
	for (j = 0L; j < kid0->numU; j++)
	    kid0->points[i][j] = tmp.points[i][j];
    for (i = 0L; i < kid0->orderU + kid0->numU; i++)
	kid0->kvU[i] = tmp.kvU[i];
    for (i = 0L; i < kid0->orderV + kid0->numV; i++)
	kid0->kvV[i] = tmp.kvV[i];

    /* Second half */

    splitPt++;
    *kid1 = *parent;

    kid1->numU = dirflag ? tmp.numU - splitPt : parent->numU;
    kid1->numV = dirflag ? parent->numV : tmp.numV - splitPt;
    kid1->kvU = kid1->kvV = NULL;
    kid1->points = NULL;
    AllocNurb( kid1, NULL, NULL );

    for (i = 0L; i < kid1->numV; i++)	/* Copy the point and kv data */
	for (j = 0L; j < kid1->numU; j++)
	    kid1->points[i][j]
		= tmp.points[dirflag ? i: (i + splitPt) ][dirflag ? (j + splitPt) : j];
    for (i = 0L; i < kid1->orderU + kid1->numU; i++)
	kid1->kvU[i] = tmp.kvU[dirflag ? (i + splitPt) : i];
    for (i = 0L; i < kid1->orderV + kid1->numV; i++)
	kid1->kvV[i] = tmp.kvV[dirflag ? i : (i + splitPt)];

    /* Construct new corners on the boundry between the two kids */
    MakeNewCorners( parent, kid0, kid1, dirflag );

    FreeNurb( &tmp );	    /* Get rid of refined parent */
}

/*
 * Test if a particular row or column of control points in a mesh
 * is "straight" with respect to a particular tolerance.  Returns true
 * if it is.
 */

#define GETPT( i )  (( dirflag ? &(n->points[crvInd][i]) : &(n->points[i][crvInd]) ))

static Boolean
IsCurveStraight( NurbSurface * n,
		 double tolerance,
		 long crvInd,
		 Boolean dirflag )  /* If true, test in U direction, else test in V */
{
    Point3 p, vec, prod;
    Point3 cp, e0;
    long i, last;
    double linelen, dist;

    /* Special case: lines are automatically straight. */
    if ((dirflag ? n->numU : n->numV) == 2L)
	return( TRUE );

    last = (dirflag ? n->numU : n->numV) - 1L;
    ScreenProject( GETPT( 0L ), &e0 );

    /* Form an initial line to test the other points against (skiping degen lines) */

    linelen = 0.0;
    for (i = last; (i > 0L) && (linelen < EPSILON); i--)
    {
	ScreenProject( GETPT( i ), &cp );
	(void) V3Sub( &cp, &e0, &vec );

	linelen = sqrt( V3SquaredLength( &vec ) );
    }

    DIVPT( vec, linelen );

    if (linelen > EPSILON)	/* If no non-degenerate lines found, it's all degen */
	for (i = 1L; i <= last; i++)
	{
	    /* The cross product of the vector defining the
	     * initial line with the vector of the current point
	     * gives the distance to the line. */
	    ScreenProject( GETPT( i ), &cp );
	    (void) V3Sub( &cp,&e0,&p );

	    (void) V3Cross( &p, &vec, &prod );
	    dist = V3Length( &prod );

	    if (dist > tolerance)
		return( FALSE );
	}

    return( TRUE );
}

/*
 * Check to see if a surface is flat.  Tests are only performed on edges and
 * directions that aren't already straight.  If an edge is flagged as straight
 * (from the parent surface) it is assumed it will stay that way.
 */
static Boolean
TestFlat( NurbSurface * n, double tolerance )
{
    long i;
    Boolean straight;
    Point3 cp00, cp0n, cpn0, cpnn, planeEqn;
    double dist,d ;

    /* Check edge straightness */

    if (! n->strU0)
	n->strU0 = IsCurveStraight( n, tolerance, 0L, FALSE );
    if (! n->strUn)
	n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE );
    if (! n->strV0)
	n->strV0 = IsCurveStraight( n, tolerance, 0L, TRUE );
    if (! n->strVn)
	n->strVn = IsCurveStraight( n, tolerance, maxV(n), TRUE );

    /* Test to make sure control points are straight in U and V */

    straight = TRUE;
    if ( (! n->flatU) && (n->strV0) && (n->strVn) )
	for (i = 1L;
	     (i < maxV(n)) && (straight = IsCurveStraight( n, tolerance, i, TRUE ));
	     i++);

    if (straight && n->strV0 && n->strVn)
	n->flatU = TRUE;

    straight = TRUE;
    if ( (! n->flatV) && (n->strU0) && (n->strUn) )
	for (i = 1L;
	     (i < maxU(n)) && (straight = IsCurveStraight( n, tolerance, i, FALSE ));
	     i++);

    if (straight && n->strU0 && n->strUn)
	n->flatV = TRUE;

    if ( (! n->flatV) || (! n->flatU) )
	return( FALSE );

    /* The surface can pass the above tests but still be twisted. */

    ScreenProject( &(n->points[0L][0L]),	    &cp00 );
    ScreenProject( &(n->points[0L][maxU(n)]),	    &cp0n );
    ScreenProject( &(n->points[maxV(n)][0L]),	    &cpn0 );
    ScreenProject( &(n->points[maxV(n)][maxU(n)]),  &cpnn );

    (void) V3Sub( &cp0n, &cp00, &cp0n );    /* Make edges into vectors */

    (void) V3Sub( &cpn0, &cp00, &cpn0 );

    /*
     * Compute the plane equation from two adjacent sides, and
     * measure the distance from the far point to the plane.  If it's
     * larger than tolerance, the surface is twisted.
     */

    (void) V3Cross( &cpn0, &cp0n, &planeEqn );

    (void) V3Normalize( &planeEqn );	/* Normalize to keep adds in sync w/ mults */

    d = V3Dot( &planeEqn, &cp00 );
    dist = fabs( V3Dot( &planeEqn, &cpnn ) - d );

    if ( dist > tolerance ) /* Surface is twisted */
	return( FALSE );
    else
	return( TRUE );
}

/*
 * Turn a sufficiently flat surface into triangles.
 */
static void
EmitTriangles( NurbSurface * n )
{
    Point3 vecnn, vec0n;		/* Diagonal vectors */
    double len2nn, len20n;		/* Diagonal lengths squared */
    double u0, un, v0, vn;		/* Texture coords;

    /*
     * Measure the distance along the two diagonals to decide the best
     * way to cut the rectangle into triangles.
     */

    (void) V3Sub( &n->c00.point, &n->cnn.point, &vecnn );
    (void) V3Sub( &n->c0n.point, &n->cn0.point, &vec0n );

    len2nn = V3SquaredLength( &vecnn ); /* Use these to reject triangles */
    len20n = V3SquaredLength( &vec0n ); /* that are too small to render */

    if (MAX(len2nn, len20n) < EPSILON)
	return;				/* Triangles are too small to render */

    /*
     * Assign the texture coordinates
     */
    u0 = n->kvU[n->orderU-1L];
    un = n->kvU[n->numU];
    v0 = n->kvV[n->orderV-1L];
    vn = n->kvV[n->numV];
    n->c00.u = u0; n->c00.v = v0;
    n->c0n.u = un; n->c0n.v = v0;
    n->cn0.u = u0; n->cn0.v = vn;
    n->cnn.u = un; n->cnn.v = vn;

    /*
     * If any normals are sick, fix them now.
     */
    if ((n->c00.normLen == 0.0) || (n->cnn.normLen == 0.0) || (n->cn0.normLen == 0.0))
	FixNormals( &(n->c00), &(n->cnn), &(n->cn0) );
    if (n->c0n.normLen == 0.0)
	FixNormals( &(n->c00), &(n->c0n), &(n->cnn) );

    if ( len2nn < len20n )
    {
	(*DrawTriangle)( &n->c00, &n->cnn, &n->cn0 );
	(*DrawTriangle)( &n->c00, &n->c0n, &n->cnn );
    }
    else
    {
	(*DrawTriangle)( &n->c0n, &n->cnn, &n->cn0 );
	(*DrawTriangle)( &n->c0n, &n->cn0, &n->c00 );
    }
}

/*
 * The recursive subdivision algorithm.	 Test if the surface is flat.
 * If so, split it into triangles.  Otherwise, split it into two halves,
 * and invoke the procedure on each half.
 */
static void
DoSubdivision( NurbSurface * n, double tolerance, Boolean dirflag, long level )
{
    NurbSurface left, right;	/* ...or top or bottom. Whatever spins your wheels. */

    if (TestFlat( n, tolerance ))
    {
	EmitTriangles( n );
    }
    else
    {
	if ( ((! n->flatV) && (! n->flatU)) || ((n->flatV) && (n->flatU)) )
	    dirflag = ! dirflag;    /* If twisted or curved in both directions, */
	else			    /* then alternate subdivision direction */
	{
	    if (n->flatU)	    /* Only split in directions that aren't flat */
		dirflag = FALSE;
	    else
		dirflag = TRUE;
	}
	SplitSurface( n, &left, &right, dirflag );
	DoSubdivision( &left, tolerance, dirflag, level + 1L );
	DoSubdivision( &right, tolerance, dirflag, level + 1L );
	FreeNurb( &left );
	FreeNurb( &right );	    /* Deallocate surfaces made by SplitSurface */
    }
}

/*
 * Main entry point for subdivision */
void
DrawSubdivision( NurbSurface * surf )
{
    surf->flatV = FALSE;
    surf->flatU = FALSE;
    surf->strU0 = FALSE;
    surf->strUn = FALSE;
    surf->strV0 = FALSE;
    surf->strVn = FALSE;

    /*
     * Initialize the projected corners of the surface
     * and the normals.
     */
    DIVW( &(surf->points[0L][0L]),			 &surf->c00.point );
    DIVW( &(surf->points[0L][surf->numU-1L]),		 &surf->c0n.point );
    DIVW( &(surf->points[surf->numV-1L][0L]),		 &surf->cn0.point );
    DIVW( &(surf->points[surf->numV-1L][surf->numU-1L]), &surf->cnn.point );

    GetNormal( surf, 0L, 0L );
    GetNormal( surf, 0L, maxU(surf) );
    GetNormal( surf, maxV(surf), 0L );
    GetNormal( surf, maxV(surf), maxU(surf) );

    DoSubdivision( surf, SubdivTolerance, TRUE, 0L );
    /* Note surf is deallocated by the subdivision process */
}