#pragma ident "@(#)cgpstream.c	1.25 98/08/25 SMI"

/*
 * Copyright (c) 1997 by Sun Microsystems, Inc.
 * All Rights Reserved
 */


/******************************************************************************
 * cgpstream.c
 *
 * This file contains all the routines to convert a list of triangle mesh
 * primitives to a command stream.
 *
 *  List of functions in this file:
 *      cgpStreamChangeQuant
 *      cgpStreamMBR
 *	cgpStreamSetColor
 *      cgpStreamSetNormal
 *	cgpStreamSetState
 *      cgpStreamData
 *	cgpStreamVertex
 *	
 *
 *****************************************************************************/

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

#include "cgpi.h"
#include "cgpstream.h"
#include "zdebug.h"

#define MAX_Y_ANG 0.615479709

static CGvoid cgpCompInvSinTables();
static CGvoid cgpCompNormals();

int first_norm;
int histo_norm[16];
int num_normal_bits = 0;

/* Bit-field masks: BMASK[i] = (1<<i)-1, but works for 32 as well */
static CGuint BMASK[33] = {
    0x0,        0x1,        0x3,        0x7,
    0xF,        0x1F,       0x3F,       0x7F,
    0xFF,       0x1FF,      0x3FF,      0x7FF,
    0xFFF,      0x1FFF,     0x3FFF,     0x7FFF,
    0xFFFF,     0x1FFFF,    0x3FFFF,    0x7FFFF,
    0xFFFFF,    0x1FFFFF,   0x3FFFFF,   0x7FFFFF,
    0xFFFFFF,   0x1FFFFFF,  0x3FFFFFF,  0x7FFFFFF,
    0xFFFFFFF,  0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF,
    0xFFFFFFFF,
};

/* used to zero out the quantizations bits */
static CGuint QMASK[17] = {
    0xFFFF0000,     0xFFFF8000,     0xFFFFC000,     0xFFFFE000,
    0xFFFFF000,     0xFFFFF800,     0xFFFFFC00,     0xFFFFFE00,
    0xFFFFFF00,     0xFFFFFF80,     0xFFFFFFC0,     0xFFFFFFE0,
    0xFFFFFFF0,     0xFFFFFFF8,     0xFFFFFFFC,     0xFFFFFFFE,
    0xFFFFFFFF
};

/* different values of "close to -1" for different quantization levels */
static CGuint MinusOne[17] = {
    0xFFFF0000,     0xFFFF0000,     0xFFFFC000,     0xFFFFA000,
    0xFFFF9000,     0xFFFF8800,     0xFFFF8400,     0xFFFF8200,
    0xFFFF8100,     0xFFFF8080,     0xFFFF8040,     0xFFFF8020,
    0xFFFF8010,     0xFFFF8008,     0xFFFF8004,     0xFFFF8002,
    0xFFFF8001
};

/* functions used only in this routine */
static CGvoid cgpStreamQuant(CGPcommandstream *, CGshort);
static CGvoid cgpStreamMBR(CGPpoint *);
static CGvoid cgpStreamSetNormal(CGPcommandstream *, CGfloat, CGfloat,
                                 CGfloat, CGint, CGboolean);
static CGvoid cgpStreamSetColor(CGPcommandstream *, CGfloat, CGfloat,
                                CGfloat, CGfloat, CGboolean);
static CGvoid cgpStreamSetState(CGPcommandstream *, CGuint);
static CGvoid cgpStreamVertex(CGPcommandstream *, CGPpoint *, CGviewport *);

/* globals to this routine */
static CGuint            curState;     /* pt type, alpha */
static CGuint            cmdPos;       /* free posn in stream to write to */
static CGPcommandstream *cmdHdr;       /* always the beginning of the stream */
static CGuint            cmdBuffSize;  /* current size of buffer */

/* absolute (instead of relative) value of previous vertex */
static CGint lastX;
static CGint lastY;
static CGint lastZ;
static CGint lastR;
static CGint lastG;
static CGint lastB;
static CGint lastA;

/* quantization parameters */
static CGint GeomQuant;
static CGint ColorQuant;
static CGint NormalQuant;

/* the mesh buffer.... buffer */
static CGPmeshbuff meshBuff[16];
static CGint       meshIdx;

static void cgpCompNormals();
int CurrBuffRefs[10000];
static int CurrBuffIndex = 0;

/*******************
 * cgpStreamData
 *
 * Walks through a CGPprimlist and turns everything into a stream of
 * vertex, color, normal and state commands.  All subsequent routines
 * will process this data stream - the primitive structures are no
 * longer needed.
 *
 *  Input:
 *    CGPprimlist*	 prims 	- the list of primitives to convert
 *    CGPcommandstream** cmdPtr - will insert a pointer to the stream here
 *    CGviewport*      viewport - holds scale & offset values to move data
 *                                into [-1.0..1.0] range
 *    CGint           geomQuant - current geometry quantization amount
 *    CGint          colorQuant - current color quantization amount
 *    CGint           normQuant - current normal quantization amount
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamData(CGPprimlist *prims, CGPcommandstream **cmdPtr,
              CGviewport *viewport, CGint geomQuant, CGint colorQuant,
              CGint normQuant)
{
    CGPprimlist  *primlist;
    CGPprimitive *curPrim;
    CGuint numCurBuffPts;  /* number of points processed in current */
                           /*  ptBuffer structure */
    CGuint       numPts;
    CGPptbuffer *curPtBuff;
    CGPpoint    *curPt;
    CGubyte      lastPtType;
    CGint        firstBuffer;  /* needed because first buffer is allocated */
                               /*  in CGPprimitive                         */
    CGbitfield   primFlags;

    /* initialize static values */
    curState = 0;
    meshIdx = 0;
    CurrBuffIndex = 0;

    /* set to large values so first time always use an absolute value */
    lastX = 327680;
    lastY = 327680;
    lastZ = 327680;
    lastR = 327680;
    lastG = 327680;
    lastB = 327680;
    lastA = 327680;

    /* set to default values */
    GeomQuant = geomQuant;
    ColorQuant = colorQuant;
    NormalQuant = normQuant;

    /* set up normal table */
    cgpCompNormals();

    /* allocate first buffer */
    cmdHdr = (CGPcommandstream *)
        malloc(sizeof(CGPcommandstream) * CGP_DATA_STREAM_BUFF_SIZE);
    cmdBuffSize = CGP_DATA_STREAM_BUFF_SIZE;

    primlist = prims;
    lastPtType = -1;

    /* first command is always SET_TABLE and will be filled in by the */
    /*  the compression routine.                                      */
    cmdHdr->dataType = CGP_STREAM_CMD_SET_TABLE;
    cmdHdr->data = NULL;
    cmdPos = 0;

    while (primlist != NULL)
    {
        CGPprimlist *lastPrim = primlist;


        /* at this point all primitives are a triangle mesh */
        curPrim = &primlist->primitive;
        primFlags = curPrim->primFlags;

        /* issues set normal if needed */
        if (primFlags & CGP_GLOBAL_NORMAL_VALID)
        {
            CGP_INC_CMD_POS;
            cgpStreamSetNormal(&cmdHdr[cmdPos], curPrim->globalNormalX,
               curPrim->globalNormalY, curPrim->globalNormalZ, NormalQuant,1);
        }

        /* issues set color if needed */
        if (primFlags & CGP_GLOBAL_COLOR_VALID)
        {
            /* issue a setState */
            if (primFlags & CGP_GLOBAL_ALPHA_VALID)
            {
                CGP_INC_CMD_POS;
                cgpStreamSetState(&cmdHdr[cmdPos],
                                  CGP_PT_TYPE_VTX_COLOR_ALPHA);
            }
            else
            {
                CGP_INC_CMD_POS;
                cgpStreamSetState(&cmdHdr[cmdPos], 0);
            }

            CGP_INC_CMD_POS;
            cgpStreamSetColor(&cmdHdr[cmdPos], curPrim->globalColorR,
               curPrim->globalColorG, curPrim->globalColorB,
               curPrim->globalAlpha, CGP_FALSE);
        }

        /* 0x3007 - all bits needed to process setState command */
        if ((primFlags & 0x3007) != lastPtType)
        {
            CGP_INC_CMD_POS;
            cgpStreamSetState(&cmdHdr[cmdPos], primFlags & 0x3007);

            lastPtType = primFlags & 0x3007;
        }

        curPtBuff = &curPrim->ptBuff;
        curPt = curPtBuff->pts;

        /* the first point buffer is allocated with the primitive structure */
        firstBuffer = 1;

        for (numCurBuffPts = 0, numPts = 0; numPts < curPrim->numPts; numPts++)
        {
            CGP_INC_CMD_POS;
            cgpStreamVertex(&cmdHdr[cmdPos], curPt, viewport);

            if (curState & CGP_CMD_SET_STATE_BNV)
            {
                if (!(curState & CGP_CMD_SET_STATE_MBR_N))
                {
                    CGP_INC_CMD_POS;
                    cgpStreamSetNormal(&cmdHdr[cmdPos], curPt->nX, curPt->nY,
                        curPt->nZ, NormalQuant,
                        (curPt->header & CGP_POINT_HDR_NORMAL_ABSOLUTE ?
                          CGP_TRUE : CGP_FALSE));
                }
                else
                    curState &= ~CGP_CMD_SET_STATE_MBR_N;
            }

            if (curState & CGP_CMD_SET_STATE_BCV)
            {
                if (!(curState & CGP_CMD_SET_STATE_MBR_C))
                {
                    CGP_INC_CMD_POS;
                    cgpStreamSetColor(&cmdHdr[cmdPos], curPt->r, curPt->g,
                        curPt->b, curPt->a,
                        (curPt->header & CGP_POINT_HDR_COLOR_ABSOLUTE ?
                         CGP_TRUE : CGP_FALSE));
                }
                else
                    curState &= ~CGP_CMD_SET_STATE_MBR_C;
            }

            /* see if points continue onto next ptBuffer structure */
            if (++numCurBuffPts == CGP_PT_BUFFER_SIZE)
            {
                CGPptbuffer *tmp = curPtBuff;

                curPtBuff = curPtBuff->nextPtBuffer;
                curPt = &curPtBuff->pts[0];
                numCurBuffPts = 0;

                /* free old pt buffer */
                if (firstBuffer)
                    firstBuffer = 0;
                else
                    free(tmp);
             }
             else
                 curPt++;

        } /* for (numCurBuffPts... */

        if (!firstBuffer)
            free(curPtBuff);

        primlist = primlist->nextPrimitive;

        /* free just processed primitive */
        free(lastPrim);

    } /* while (primlist != NULL) */


    /* add a NULL command to signify the end of the buffer */    
    CGP_INC_CMD_POS;
    cmdHdr[cmdPos].dataType = CGP_STREAM_CMD_NULL;
    cmdHdr[cmdPos].data = NULL;

    *cmdPtr = cmdHdr;
}

#ifdef COMMENT
/*****************************************************************************

                        Normal Encoding Parameterization

First, points on a unit radius sphere are parameterized by two angles,
th and psi, using usual spherical coordinates. th is the angle about
the y axis, psi is the inclination to the plane containing the point.
The mapping between rectangular and spherical coordinates is:

    x = cos(th)*cos(psi)
    y = sin(th)*cos(psi)
    z = sin(psi)

Points on sphere are folded first by octant, and then by sort order
of xyz into one of six hexes. All the table encoding takes place in
the positive octant, in the region bounded by the half spaces:

    x >= z
    z >= y
    y >= 0

This triangular shaped patch runs from 0 to 45 degrees in th, and
from 0 to as much as 0.615479709 (MAX_Y_ANG) in psi. The xyz bounds
of the patch is:

    (1, 0, 0)  (1/sqrt(2), 0, 1/sqrt(2))  (1/sqrt(3), 1/sqrt(3), 1/sqrt(3))

When dicing this space up into descrete points, the choice for y is
linear quantization in psi.  This means that if the y range is to be
divided up into n segments, the angle of segment j is:

    psi(j) = MAX_Y_ANG*(j/n)

The y heigt of the patch (in arc length) is
*not* the same as the xz dimension. However, the subdivision quantization
needs to treat xz and y equally. To achieve this, the th angles are
re-parameterized as reflected psi angles.  That is, the i-th point's
th is:

    th(i) = asin(tan(psi(i))) = asin(tan(MAX_Y_ANG*(i/n)))

to go the other direction, the angle th coresponds to the real index r
in the same 0-n range as i):

    r(th) = n*atan(sin(th))/MAX_Y_ANG

rounded to the nearest integer, this gives the closest integer index i
to the xz angle th. Because the triangle has a straight edge on the
line x=z, it is more intuitive to index the xz angles in reverse
order.  Thus the two equations above are replaced by:

    th(i) = asin(tan(psi(i))) = asin(tan(MAX_Y_ANG*((n-i)/n)))

    r(th) = n*(1 - atan(sin(th))/MAX_Y_ANG)


----

Each level of quantization subdivides the triangular patch twice as densly.
The case in which only the three vertices of the triangle are present is
the first logical stage of representation, but because of how the table
is encoded the first usable case starts one level of sub-division later.
This three point level has an n of 2, by our above conventions.

Encode level 0 adds one subdivision in y and xz, and one point. All
even levels only fill in every other point in the latice, odd levels
fill in all rectangular latice points. n only increases every two
levels, thus it is useful to define an auxzilary variable f, which
when 0 indicates a half filled level, and 1 is a fully filled level.
A level number is converted to n and f by:

    n = 1<<(1 + (level>>1))
    f = level & 1

For half filled levels, xz lines of even y have even number xz vertices
filled in, odd have odd numbered xz vertices filled in. (Note that this
is reversed from what would be the case if xz indices were not reversed
to start at 0 at 45 degrees, and going to n at 0 degrees.) So for vertice
<i j> (xz index, y index) (for fixed level, n and f), for even level (f==0):

    vertex <i j> filled? only if (i&1) == (j&1)

----

            Brute force algorithm to find closest latice
            point to a given normal at a given level

First convert the input normal into the hextant
Loop through all valid latice points, doting the quantized normal
    with the input normal, return latice normal with largest dot product.

----
******************************************************************************/
#endif

/*
 *  The inverse sin table for each of the quantization levels
 *
 *  The input to each table will be fixed point ny or nz, shifted down
 *  by a specified number of bits. The size of the table should be as
 *  small as possible, but with the delta sin at the greatest angle
 *  still smaller than the delta sin between the last two angles to be
 *  encoded.
 *
 *  Example: the inverse sin y table has a maximum angle of 0.615479709,
 *  with a delta angle of 0.615479709/128. The difference between the last
 *  two angles to be encoded is:
 *  sin(0.615479709*128.0/128.0) - sin(0.615479709*127.0/128.0) = 0.003932730
 *
 *  Using 8 significent bits below the binary point, fixed point can
 *  represent sin in increments of 0.003906250, just slightly smaller.
 *
 *  However, because the maximum y angle sin is 0.577350269, rather than 256
 *  table entries, only 148 are needed.
 */

#define UNITY_14 16384.0
short   inv_sin_y_00[  5];  /* Shift 11 */
short   inv_sin_y_02[ 10];  /* Shift 10 */
short   inv_sin_y_04[ 19];  /* Shift  9 */
short   inv_sin_y_06[ 37];  /* Shift  8 */
short   inv_sin_y_08[ 74];  /* Shift  7 */
short   inv_sin_y_10[148];  /* Shift  6 */


CGdouble  angles_y[129];
CGshort   y_sin[129];


static  CGdouble  gc_normals[12][65][65][3];

/* Compute gc_normals[n][j][i][3] */

static CGvoid
cgpCompNormals()
{
    CGint level, n, f, i, j;
    CGint inx, iny, inz;
    CGdouble th, psi, qnx, qny, qnz;

    for (level = 0; level < 12; level++) {
        n = 1<<(1 + (level>>1));
        f = level & 1;
        for (j = 0; j <= n; j++) {
            for (i = 0; i <= n; i++) {
                if (i+j > n) continue;
                if (f == 0 && (i&1) != (j&1)) continue;
                psi = MAX_Y_ANG*(j/((double) n));
                th = asin(tan(MAX_Y_ANG*((n-i)/((double) n))));

                qnx = cos(th)*cos(psi);
                qny = sin(psi);
                qnz = sin(th)*cos(psi);

                /*
                 *  Convert the floating point normal to s1.14 bit notation,
                 *  then back again.
                 */
                qnx = qnx*16384.0; inx = qnx; qnx = inx; qnx = qnx/16384.0;
                qny = qny*16384.0; iny = qny; qny = iny; qny = qny/16384.0;
                qnz = qnz*16384.0; inz = qnz; qnz = inz; qnz = qnz/16384.0;
                gc_normals[level][j][i][0] = qnx;
                gc_normals[level][j][i][1] = qny;
                gc_normals[level][j][i][2] = qnz;
            }
        }
    }

    /* Go compute the sine tables also */
    cgpCompInvSinTables();

}  /* end of cgpCompNormals */

/*
 *  Compute the inverse sin tables for y, for levels 0 through 11.
 *
 *  At any level of compression, there is a fixed number of different
 *  y angles (between 0 and MAX_Y_ANG).  We wish to build a fast
 *  inverse table to take an ny and obtain the closest quantized y angle.
 *  We build the inverse table to have slightly more entries as
 *  there are y angles at any particular level; this ensures that the
 *  inverse look-up will get us within one angle of the right one.
 *
 *  Each different (pair of) levels has roughfly half as many entries
 *  as the last, all this is taken into account.
 *
 */

static void
cgpCompInvSinTables()
{
    CGint         is, iy;
    CGdouble      y, ay;

    /*
     *  Build table of floating point angles.
     *  At maximum resolution, there are 65 descrete angles used,
     *  but we need twice as many for thresholding between angles.
     */
    for (iy = 0; iy < 129; iy++) {
        angles_y[iy] = ay = MAX_Y_ANG*iy/128.0;
    }


    /*
     *  Now that we have the angles, compute the point on the sphere
     *  (in s1.14 fixedpoint) for ... .
     */
    for (iy = 0; iy < 129; iy++) {
        /* Convert to s1.14 bit notation */
        y_sin[iy] = UNITY_14*sin(angles_y[iy]);
    }

    /*
     *  Compute the inverse sin y table for levels 10 and 11.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/128.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-8))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<6)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<1)] - is) > abs(y_sin[(iy<<1)-(1<<1)] - is))
                iy = iy-1;
        }
        if (iy < 64) {
            if (abs(y_sin[(iy<<1)] - is) > abs(y_sin[(iy<<1)+(1<<1)] - is))
                iy = iy+1;
        }
        inv_sin_y_10[is>>6] = iy;
    }

    /*
     *  Compute the inverse sin y table for levels 8 and 9.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/64.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-7))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<5)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<2)] - is) > abs(y_sin[(iy-1)<<2] - is))
                iy = iy-1;
        }
        if (iy < 32) {
            if (abs(y_sin[(iy<<2)] - is) > abs(y_sin[(iy+1)<<2] - is))
                iy = iy+1;
        }
        inv_sin_y_08[is>>7] = iy;
    }

    /*
     *  Compute the inverse sin y table for levels 6 and 7.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/32.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-6))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<4)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<3)] - is) > abs(y_sin[(iy-1)<<3] - is))
                iy = iy-1;
        }
        if (iy < 16) {
            if (abs(y_sin[(iy<<3)] - is) > abs(y_sin[(iy+1)<<3] - is))
                iy = iy+1;
        }
        inv_sin_y_06[is>>8] = iy;
    }

    /*
     *  Compute the inverse sin y table for levels 4 and 5.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/16.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-5))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<3)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<4)] - is) > abs(y_sin[(iy-1)<<4] - is))
                iy = iy-1;
        }
        if (iy < 8) {
            if (abs(y_sin[(iy<<4)] - is) > abs(y_sin[(iy+1)<<4] - is))
                iy = iy+1;
        }
        inv_sin_y_04[is>>9] = iy;
    }

    /*
     *  Compute the inverse sin y table for levels 2 and 3.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/8.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-4))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<2)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<5)] - is) > abs(y_sin[(iy-1)<<5] - is))
                iy = iy-1;
        }
        if (iy < 2) {
            if (abs(y_sin[(iy<<5)] - is) > abs(y_sin[(iy+1)<<5] - is))
                iy = iy+1;
        }
        inv_sin_y_02[is>>10] = iy;
    }

    /*
     *  Compute the inverse sin y table for levels 0 and 1.
     *  The integer sin varies from 0 to y_sin[128], in incs of 1/4.
     */
    for (is = 0; is <= y_sin[128]; is += (1<<(14-3))) {
        y = is/UNITY_14;
        ay = asin(y);  /* ay is the angle who's sign is is */
        iy = ay*(1<<1)/MAX_Y_ANG;
        if (iy > 0) {
            if (abs(y_sin[(iy<<6)] - is) > abs(y_sin[(iy-1)<<6] - is))
                iy = iy-1;
        }
        if (iy < 4) {
            if (abs(y_sin[(iy<<6)] - is) > abs(y_sin[(iy+1)<<6] - is))
                iy = iy+1;
        }
        inv_sin_y_00[is>>11] = iy;
    }

}  /* end of cgpCompInvSinTables */

/*******************
 * cgpStreamSetNormal
 *
 * Inserts a SetNormal command into the data stream.
 *
 *  Input:
 *    CGPcommandstream* cmd - current position in command stream
 *    CGfloat           nx
 *    CGfloat           ny
 *    CGfloat           nz
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamSetNormal(CGPcommandstream *cmd, CGfloat nx, CGfloat ny, CGfloat nz,
                   CGint quant, CGboolean absolute)
{
    CGPsetnormal* normalEntry;
    CGint    level;
    CGint    oct, sext, u, v;
    CGint    inx, iny, inz;
    CGint    yIndex;
    CGint n, f, j, i;
    CGint besti = -100;
    CGint bestj = -100;
    CGfloat tmp;
    CGdouble dot;
    CGdouble bestDot = -100.0;
    CGint index, special;

    cmd->dataType = CGP_STREAM_CMD_SET_NORMAL;
    cmd->data = (CGvoid *) malloc(sizeof(CGPsetnormal));
    normalEntry = (CGPsetnormal *) cmd->data;
    normalEntry->flags = 0;

    /* quant must be an even number */
    quant &= ~0x1;

    /* currently can compress down to only 7 bits */
    level = quant - 7;

    if (level < 0)
        level = 0;

    /* normalize the fixed point normal to the positive signed octant */
    oct = 0;
    if (nx < 0.0)
    {
        oct |= 4;
        nx = -nx;
    }

    if (ny < 0.0)
    {
        oct |= 2;
        ny = -ny;
    }

    if (nz < 0.0)
    {
        oct |= 1;
        nz = -nz;
    }

    /* normalize the normal to proper sextant of octant */
    sext = 0;
    if (nx < ny)
    {
        tmp = nx;
        nx = ny;
        ny = tmp;
        sext |= 1;
    }

    if (nz < ny)
    {
        tmp = ny;
        ny = nz;
        nz = tmp;
        sext |= 2;
    }

    if (nx < nz)
    {
        tmp = nx;
        nx = nz;
        nz = tmp;
        sext |= 4;
    }

    /* convert the floating point y normal to s1.14 notation */
    tmp = ny * 16384.0;
    iny = (CGint) tmp;

    /* now index into table (which is sin(ay)) based on iny */
    switch (level)
    {
        case 0:
        case 1:
            yIndex = inv_sin_y_00[iny >> 11];
            break;

        case 2:
        case 3:
            yIndex = inv_sin_y_02[iny >> 10];
            break;

        case 4:
        case 5:
            yIndex = inv_sin_y_04[iny >> 9];
            break;

        case 6:
        case 7:
            yIndex = inv_sin_y_06[iny >> 8];
            break;

        case 8:
        case 9:
            yIndex = inv_sin_y_08[iny >> 7];
            break;

        case 10:
        case 11:
            yIndex = inv_sin_y_10[iny >> 6];
            break;
    }

    /* search the three xz rows near y for the best match */
    n = 1 << (1 + (level >> 1));
    f = level & 1;
    for (j = yIndex - 1; j < yIndex + 1 && j <= n; j++)
    {
        if (j < 0)
            continue;

        for (i = 0; i <= n; i++)
        {
            if (i + j > n)
                continue;

            if (f == 0 && (i & 1) != (j & 1))
                continue;

            dot = nx * gc_normals[level][j][i][0] +
                  ny * gc_normals[level][j][i][1] +
                  nz * gc_normals[level][j][i][2];

            if (dot > bestDot)
            {
                bestDot = dot;
                besti = i;
                bestj = j;
            }
        }
    }  

    /* convert besti and bestj to standard grid form */
    normalEntry->u = u = besti << (5 - (level >> 1));
    normalEntry->v = v = bestj << (5 - (level >> 1));

    /* compute the 18 bit representation, check for special normals first */
    if (u == 64 && v == 0)
    {
        /* six coordinate axis case */
        if (sext == 0 || sext == 2)
            index = CGP_SPECIAL_NORMAL_X_AXIS | ((oct & 4) ? 0x2000 : 0);
        else if (sext == 3 || sext == 1)
            index = CGP_SPECIAL_NORMAL_Y_AXIS | ((oct & 2) ? 0x2000 : 0);
        else if (sext == 5 || sext == 4)
            index = CGP_SPECIAL_NORMAL_Z_AXIS | ((oct & 1) ? 0x2000 : 0);

        normalEntry->index = index;
        normalEntry->flags |= CGP_SET_NORMAL_FLAG_SPECIAL;
    }
    else
    {
        if (u == 0 && v == 64)
        {
            /* eight mid-point case */
            index = CGP_SPECIAL_NORMAL_MID_OCT | (oct << 13);

            normalEntry->index = index;
            normalEntry->flags |= CGP_SET_NORMAL_FLAG_SPECIAL;
        }
        else if (u == 0 && v == 0)
        {
            /* twelve mid-circle case */
            if (sext == 2 || sext == 3)
                index = CGP_SPECIAL_NORMAL_CIR_XY | ((oct & 4) ? 0x2000 : 0) |
                    ((oct & 2) ? 0x1000 : 0);
            else if (sext == 1 || sext == 5)
                index = CGP_SPECIAL_NORMAL_CIR_YZ | ((oct & 2) ? 0x2000 : 0) |
                    ((oct & 1) ? 0x1000 : 0);
            else if (sext == 4 || sext == 0)
                index = CGP_SPECIAL_NORMAL_CIR_XZ | ((oct & 4) ? 0x2000 : 0) |
                    ((oct & 1) ? 0x1000 : 0);

            normalEntry->index = index;
            normalEntry->flags |= CGP_SET_NORMAL_FLAG_SPECIAL;
        }
        else
        {
            /* non-special normal, encode using general rule */
            normalEntry->index = (sext << (15 - (11 - level))) |
                    (oct << (12 - (11 - level))) |
                    ((u >> (5 - (level >> 1))) << (6 - 5 + (level >> 1))) |
                    (v >> (5 - (level >> 1)));
        }
    }
/* TODO - bug such that all normals must be absolute */
absolute = 1;

    if (absolute)
        normalEntry->flags |= CGP_SET_NORMAL_FLAG_ABSOLUTE;

    if (normalEntry->flags & CGP_SET_NORMAL_FLAG_SPECIAL)
    {
        normalEntry->normalLen = 6;
        normalEntry->normalShift = 0;
    }
    else
    {
        /* for absolute, sext/oct always assumed */
        normalEntry->normalShift = 5 - (level >> 1);
        normalEntry->normalLen = (quant - 6) / 2 + normalEntry->normalShift;
    }

    normalEntry->sext = sext;
    normalEntry->oct = oct;
    normalEntry->normalTagLen = 0;
    normalEntry->normalTag = 0;
}


/*******************
 * cgpStreamSetColor
 *
 * Inserts a SetColor command into the data stream.
 *
 *  Input:
 *    CGPcommandstream* cmd - current position in command stream
 *    CGfloat           r
 *    CGfloat           g
 *    CGfloat           b
 *    CGfloat           a		- if alpha is not used,
 *                                        this can be garbage
 *    CGboolean         absolute	- CGP_FALSE for a relative delta to be
 *                                        inserted if possible, all else
 *                                        means absolute value desired
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamSetColor(CGPcommandstream *cmd, CGfloat r, CGfloat g, CGfloat b,
                                              CGfloat a, CGboolean absolute)
{
    CGPsetcolor* colorEntry;
    CGint dr, dg, db, da;
    CGint ir, ig, ib, ia;

    cmd->dataType = CGP_STREAM_CMD_SET_COLOR;
    cmd->data = (CGvoid *) malloc(sizeof(CGPsetcolor));
    colorEntry = (CGPsetcolor*) cmd->data;

    ir = r * 32768.0;
    ig = g * 32768.0;
    ib = b * 32768.0;
    ia = a * 32768.0;

    /* handles 1.0 case */
    if (ir == 0x8000)
        ir = 0x7fff;
    if (ig == 0x8000)
        ig = 0x7fff;
    if (ib == 0x8000)
        ib = 0x7fff;
    if (ia == 0x8000)
        ia = 0x7fff;

    /* quantize to specified amount */

    /* check for -1.0 case - must be based on current quant level */
    if (ir == 0xffff8000)
        ir = MinusOne[ColorQuant];
    else
        ir &= QMASK[ColorQuant];

    if (ig == 0xffff8000)
        ig = MinusOne[ColorQuant];
    else
        ig &= QMASK[ColorQuant];

    if (ib == 0xffff8000)
        ib = MinusOne[ColorQuant];
    else
        ib &= QMASK[ColorQuant];

    if (ia == 0xffff8000)
        ia = MinusOne[ColorQuant];
    else
        ia &= QMASK[ColorQuant];

    /* calculate deltas (used for relative points) */
    dr = ir - lastR;
    dg = ig - lastG;
    db = ib - lastB;

    if (curState & CGP_CMD_SET_STATE_CAP)
        da = ia - lastA;
    else
	da = 0;

    /* check for special cases */
    if (dr == 0 && dg == 0 && db == 0 && da == 0)
        absolute = CGP_TRUE;

    /* make sure can specify this delta in 15 bits */
    if (absolute || abs(dr) > CGP_MAX_DELTA || abs(dg) > CGP_MAX_DELTA ||
                       abs(db) > CGP_MAX_DELTA || abs(da) > CGP_MAX_DELTA)
    {
        lastR = ir;
        lastG = ig;
        lastB = ib;

        if (curState & CGP_CMD_SET_STATE_CAP)
            lastA = ia;

        colorEntry->ir = ir;
        colorEntry->ig = ig;
        colorEntry->ib = ib;
        colorEntry->ia = ia;
        colorEntry->colorLen = 16;
        colorEntry->colorShift = 0;
        colorEntry->absolute = 1;
    }
    else
    {
        lastR += dr;
        lastG += dg;
        lastB += db;

        if (curState & CGP_CMD_SET_STATE_CAP)
            lastA += da;

        colorEntry->ir = dr;
        colorEntry->ig = dg;
        colorEntry->ib = db;
        colorEntry->ia = da;
        colorEntry->colorLen = 16;
        colorEntry->colorShift = 0;
        colorEntry->absolute = 0;
    }
}


/*******************
 * cgpStreamSetState
 *
 * Inserts a SetState command into the data stream.
 * SetState controls the point type (sending color and/or normal
 * values with each vertex) and if ALPHA is valid in the color
 *
 *  Input:
 *    CGPcommandstream* cmd - current position in command stream
 *    CGuint            ptType
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamSetState(CGPcommandstream *cmd, CGuint ptType)
{
    CGPsetstate *setState;

    cmd->dataType = CGP_STREAM_CMD_SET_STATE;
    cmd->data = (CGvoid *) malloc(sizeof(CGPsetstate));
    setState = (CGPsetstate*) cmd->data;

    curState &=
        ~(CGP_CMD_SET_STATE_BNV | CGP_CMD_SET_STATE_BCV |
          CGP_CMD_SET_STATE_CAP );

    if (ptType & CGP_PT_TYPE_VTX_NORMAL)
    	curState |= CGP_CMD_SET_STATE_BNV;

    if (ptType & CGP_PT_TYPE_VTX_COLOR)
    	curState |= CGP_CMD_SET_STATE_BCV;

    if (ptType & CGP_PT_TYPE_VTX_COLOR_ALPHA)
    	curState |= CGP_CMD_SET_STATE_CAP;

    setState->state = curState;
}


/*******************
 * cgpStreamVertex
 *
 * Inserts a SetVertex command into the data stream.
 *
 *  Input:
 *    CGPcommandstream* cmd - current position in command stream
 *    CGPpoint*          pt
 *    CGviewport*  viewport - holds scale & offset values to move data
 *                            into [-1.0..1.0] range
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamVertex(CGPcommandstream *cmd, CGPpoint *pt, CGviewport *viewport)
{
    CGPvertex *cpt;
    CGint ix, iy, iz;
    CGint dx, dy, dz;   /* position delta's */
    float oldx, oldy, oldz, backx, backy, backz;
    int shift;

    /* check for quntization change */
    if (pt->header & (CGP_POINT_HDR_QUANT_CHANGE_GEOM   |
                      CGP_POINT_HDR_QUANT_CHANGE_NORMAL |
                      CGP_POINT_HDR_QUANT_CHANGE_COLOR))
    {
        if (pt->header & CGP_POINT_HDR_QUANT_CHANGE_GEOM &&
            pt->geomQuant != GeomQuant)
        {
            cmd->dataType = CGP_STREAM_CMD_SET_QUANT_GEOM;
            cgpStreamQuant(cmd, pt->geomQuant);
            CGP_INC_CMD_POS;
            GeomQuant = pt->geomQuant;
        }
        
        if (pt->header & CGP_POINT_HDR_QUANT_CHANGE_NORMAL &&
            pt->normalQuant != NormalQuant)
        {
            cmd = &cmdHdr[cmdPos];
            cmd->dataType = CGP_STREAM_CMD_SET_QUANT_NORMAL;
            cgpStreamQuant(cmd, pt->normalQuant);
            CGP_INC_CMD_POS;
            NormalQuant = pt->normalQuant;
        }
        
        if (pt->header & CGP_POINT_HDR_QUANT_CHANGE_COLOR &&
            pt->colorQuant != ColorQuant)
        {
            cmd = &cmdHdr[cmdPos];
            cmd->dataType = CGP_STREAM_CMD_SET_QUANT_COLOR;
            cgpStreamQuant(cmd, pt->colorQuant);
            CGP_INC_CMD_POS;
            ColorQuant = pt->colorQuant;
        }

        cmd = &cmdHdr[cmdPos];
    }

    /* normalize the data */
    oldx = pt->x;
    oldy = pt->y;
    oldz = pt->z;
    pt->x = (pt->x - viewport->xo) * viewport->scale;
    pt->y = (pt->y - viewport->yo) * viewport->scale;
    pt->z = (pt->z - viewport->zo) * viewport->scale;

    if (pt->header & CGP_POINT_HDR_REF)
    {
       /* TODO - this should be all that is needed
        cgpStreamMBR(pt);
        return;
    } */
/* TODO - remove once mesh push logic in meshifier is fixed */
    CGint   mbIdx = CGP_GET_POINT_HDR_REF_NUMBER(pt->header);
    CGint   arrayIdx;
   /* adjust index to proper place in fifo */
    arrayIdx = ((meshIdx - 1) - mbIdx) & 0xf;

    if (fabs(meshBuff[arrayIdx].x - pt->x) <= .00001 &&
        fabs(meshBuff[arrayIdx].y - pt->y) <= .00001 &&
        fabs(meshBuff[arrayIdx].z - pt->z) <= .00001)
        {
            cgpStreamMBR(pt);
            return;
        }
        else {
            pt->header &= ~CGP_POINT_HDR_REF;
            /* Don't print this out for now 
            printf("Bad mesh buffer reference\n"); */
        }
    }

    cmd->dataType = CGP_STREAM_CMD_VTX;
    cmd->data = (CGvoid *) malloc(sizeof(CGPvertex));
    cpt = (CGPvertex*) cmd->data;

    cpt->header = pt->header;

    shift = 16 - GeomQuant;

    /* convert floating point to s.15 fixed (all floating point values */
    /*  are between (-1.0, 1.0) )                                      */ 

    ix = pt->x * 32768.0;
    iy = pt->y * 32768.0;
    iz = pt->z * 32768.0;

    /* handles 1.0 case */
    if (ix == 0x8000)
        ix = 0x7fff;
    if (iy == 0x8000)
        iy = 0x7fff;
    if (iz == 0x8000)
        iz = 0x7fff;

    /* quantize data */
    /* check for -1.0 case - must be based on current quant level */
    if (ix == 0xffff8000)
        ix = MinusOne[GeomQuant];
    else
        ix &= QMASK[GeomQuant];

    if (iy == 0xffff8000)
        iy = MinusOne[GeomQuant];
    else
        iy &= QMASK[GeomQuant];

    if (iz == 0xffff8000)
        iz = MinusOne[GeomQuant];
    else
        iz &= QMASK[GeomQuant];

    /* save info if this is a mesh buffer push */
    if (pt->header & CGP_POINT_HDR_PUSH)
    {
/* TODO - remove pushing of vtx info once mesh push is fixed */
        meshBuff[meshIdx].x = pt->x;
        meshBuff[meshIdx].y = pt->y;
        meshBuff[meshIdx].z = pt->z;
        meshBuff[meshIdx].nX = pt->nX;
        meshBuff[meshIdx].nY = pt->nY;
        meshBuff[meshIdx].nZ = pt->nZ;
        meshBuff[meshIdx].r = pt->r;
        meshBuff[meshIdx].g = pt->g;
        meshBuff[meshIdx].b = pt->b;
        meshBuff[meshIdx++].a = pt->a;

        if (meshIdx == 16)
            meshIdx = 0;
    }

#if 0
    backx = ((float) ix) / 32768.0;
    backy = ((float) iy) / 32768.0;
    backz = ((float) iz) / 32768.0;

    backx = backx / viewport->scale + viewport->xo;
    backy = backy / viewport->scale + viewport->yo;
    backz = backz / viewport->scale + viewport->zo;

    zdo6 printf("oldpt %f %f %f, back %f %f %f\n",
		 oldx, oldy, oldz,
		 backx, backy, backz);

    if (fabs(backx - pt->qx) <= .001 &&
	fabs(backy - pt->qy) <= .001 &&
	fabs(backz - pt->qz) <= .001){
    }
    else {
      zdo printf("cgp quantized pt %f %f %f doesn't match precg %f %f %f\n",
		 backx, backy, backz, 
		 pt->qx, pt->qy, pt->qz);
    }

#endif

    dx = ix - lastX;
    dy = iy - lastY;
    dz = iz - lastZ;

    zdo6 printf("shift %d oldpt %f %f %f, cubed %f %f %f\n",
	       shift, oldx, oldy, oldz,
	       pt->x, pt->y, pt->z);
    zdo6 printf("ixyz %d %d %d, dxyz %d %d %d last %d %d %d \n",
		ix, iy, iz, dx, dy, dz, lastX, lastY, lastZ);
		
    if (dx == 0 && dy == 0 && dz == 0)
        cpt->header |= CGP_POINT_HDR_GEOMETRY_ABSOLUTE;

    /* make sure can specify this delta in 15 bits */
    if ((cpt->header & CGP_POINT_HDR_GEOMETRY_ABSOLUTE ||
	      abs(dx) > CGP_MAX_DELTA || abs(dy) > CGP_MAX_DELTA ||
	      abs(dz) > CGP_MAX_DELTA))

    {
      zdo6 printf("max delta exceeded, header %x\n",
		 cpt->header);
      zdo6 printf("pt %2.2f %2.2f %2.2f, ixyz %d %d %d, dxyz %d %d %d last %d %d %d \n",
		 pt->x, pt->y, pt->z, ix, iy, iz, dx, dy, dz, lastX, lastY, lastZ);
        /* check for -1.0 case - must be based on current quant level */
        cpt->ix = ix & 0xffff;
        cpt->iy = iy & 0xffff;
        cpt->iz = iz & 0xffff;
        cpt->geomLen = 16;
        cpt->geomShift = 0;

        /* set this in case we got here due to too large a delta */
        cpt->header |= CGP_POINT_HDR_GEOMETRY_ABSOLUTE;

        lastX = ix;
        lastY = iy;
        lastZ = iz;
    }
    else
    {
        /* will be making this a relative offset */
        cpt->ix = dx & 0xffff;
        cpt->iy = dy & 0xffff;
        cpt->iz = dz & 0xffff;
        cpt->geomLen = 16;
        cpt->geomShift = 0;

        lastX += dx;
        lastY += dy;
        lastZ += dz;
    }

/* TODO return error amount */
}

/*******************
 * cgpStreamMBR
 *
 * Inserts a SetMBR (Mesh Buffer Reference) command into the data stream.
 *
 *  Input:
 *    CGPpoint* vtx	- the header will contain the referenced vertex
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamMBR(CGPpoint *vtx)
{
    CGPmbr *mbr;
    CGint   mbIdx = CGP_GET_POINT_HDR_REF_NUMBER(vtx->header);
    CGint   arrayIdx;

    /* adjust index to proper place in fifo */
    arrayIdx = ((meshIdx - 1) - mbIdx) & 0xf;

    /* if color/normal different from point in buffer, issue set states */
    if (curState & CGP_CMD_SET_STATE_BNV)
    {
        if (0 || !((fabs(vtx->nX - meshBuff[arrayIdx].nX) < CGP_SMALL_FLOAT) &&
	      (fabs(vtx->nY - meshBuff[arrayIdx].nY) < CGP_SMALL_FLOAT) &&
	      (fabs(vtx->nZ - meshBuff[arrayIdx].nZ) < CGP_SMALL_FLOAT)))
        {
	    zdo printf("normals not same as mbr, setting new one\n");
	    /* insert a set normal into the data stream */
            /* cmdPos is already open from call to cgVertex */
            cgpStreamSetNormal(&cmdHdr[cmdPos], vtx->nX, vtx->nY, vtx->nZ,
            NormalQuant, CGP_FALSE);
            CGP_INC_CMD_POS;
        }

        /* don't send down normal from standard vertex processing */
        /* the mesh buffer or explicit set normal will cover it */
        curState |= CGP_CMD_SET_STATE_MBR_N;
	
    }

    if (curState & CGP_CMD_SET_STATE_BCV)
    {
        if (curState & CGP_CMD_SET_STATE_CAP)
        {
            if (!(vtx->r == meshBuff[arrayIdx].r &&
                vtx->g == meshBuff[arrayIdx].g &&
                vtx->b == meshBuff[arrayIdx].b &&
                vtx->a == meshBuff[arrayIdx].a))
            {
                /* insert a set color into the data stream */
                /* cmdPos is already open from call to cgVertex */
                cgpStreamSetColor(&cmdHdr[cmdPos], vtx->r, vtx->g, vtx->b,
                   vtx->a, CGP_FALSE);
                CGP_INC_CMD_POS;
            }

            /* don't send down color from standard vertex processing */
            curState |= CGP_CMD_SET_STATE_MBR_C;
        }
        else
        {
            if (!(vtx->r == meshBuff[arrayIdx].r &&
                vtx->g == meshBuff[arrayIdx].g &&
                vtx->b == meshBuff[arrayIdx].b))
            {
                /* insert a set color into the data stream */
                /* cmdPos is already open from call to cgVertex */
                cgpStreamSetColor(&cmdHdr[cmdPos], vtx->r, vtx->g, vtx->b,
                  vtx->a, 0);
                CGP_INC_CMD_POS;
            }

            /* don't send down color from standard vertex processing */
            curState |= CGP_CMD_SET_STATE_MBR_C;
        }
    }

    cmdHdr[cmdPos].dataType = CGP_STREAM_CMD_MBR;
    cmdHdr[cmdPos].data = (CGvoid *) malloc(sizeof(CGPmbr));
    mbr = (CGPmbr *) cmdHdr[cmdPos].data;

    mbr->mbRef = vtx->header;
    CGP_SET_POINT_HDR_REF_NUMBER(mbr->mbRef, mbIdx);

    zdo6 printf("mbref %d %f %f %f\n",
	       arrayIdx,
	       vtx->x, 
	       vtx->y,
	       vtx->z);

    CurrBuffRefs[CurrBuffIndex++] = mbIdx;
    zdo6 printf("header index %d \n", CGP_GET_POINT_HDR_REF_NUMBER(vtx->header));

    lastX = vtx->x * 32768.0;
    lastY = vtx->y * 32768.0;
    lastZ = vtx->z * 32768.0;

    /* handles 1.0 case */
    if (lastX == 0x8000)
        lastX = 0x7fff;
    if (lastY == 0x8000)
        lastY = 0x7fff;
    if (lastZ == 0x8000)
        lastZ = 0x7fff;

    /* check for -1.0 case - must be based on current quant level */
    if (lastX == 0xffff8000)
        lastX = MinusOne[GeomQuant];
    else
        lastX &= QMASK[GeomQuant];

    if (lastY == 0xffff8000)
        lastY = MinusOne[GeomQuant];
    else
        lastY &= QMASK[GeomQuant];

    if (lastZ == 0xffff8000)
        lastZ = MinusOne[GeomQuant];
    else
        lastZ &= QMASK[GeomQuant];

    zdo6 printf("last %d %d %d\n",
	       lastX, lastY, lastZ);

}


/*******************
 * cgpStreamQuant
 *
 * Inserts a SetQuant command into the data stream.  This controls
 * how much quantization of the data is allowed.
 *
 *  Input:
 *    CGPcommandstream* cmd - current position in command stream
 *    CGfloat           nx
 *    CGfloat           ny
 *    CGfloat           nz
 *
 *  Output:
 *    None.
 *
 */

CGvoid
cgpStreamQuant(CGPcommandstream *cmd, CGshort quant)
{
    CGPsetquant *setQuant;

    cmd->data = (CGvoid *) malloc(sizeof(CGPsetquant));
    setQuant = (CGPsetquant *) cmd->data;

    setQuant->quant = quant;
}
