GEOM.C

/**********************************Module**********************************\ 
* 
* geom.c 
* 
* 3D FlowerBox screen saver 
* Geometry routines 
* 
* Copyright 1995 - 1998 Microsoft Corporation 
* 
\**************************************************************************/ 
 
#include "precomp.h" 
#pragma hdrstop 
 
// Sphere radius 
#define RADIUS 1 
 
// Values to map a 2D point onto a 3D plane 
// Base point and axes to map X and Y coordinates onto 
typedef struct _PLANE_MAP 
{ 
    PT3 base, x_axis, y_axis; 
} PLANE_MAP; 
 
// Data area used by the current geometry 
// Base points and generated points 
PT3 pts[MAXPTS], npts[MAXPTS]; 
// Scaling factor for spherical projection 
FLT vlen[MAXPTS]; 
// Normals 
PT3 normals[MAXPTS]; 
// Vertex data indices 
int index[MAXPTS*2]; 
// Triangle strip sizes 
int strip_size[MAXSIDES*MAXSUBDIV]; 
 
void InitCube(GEOMETRY *geom); 
void InitTetra(GEOMETRY *geom); 
void InitPyramids(GEOMETRY *geom); 
void InitCylinder(GEOMETRY *geom); 
void InitSpring(GEOMETRY *geom); 
 
GEOMETRY cube_geom = {InitCube}; 
GEOMETRY tetra_geom = {InitTetra}; 
GEOMETRY pyramids_geom = {InitPyramids}; 
GEOMETRY cylinder_geom = {InitCylinder}; 
GEOMETRY spring_geom = {InitSpring}; 
 
GEOMETRY *geom_table[] = 
{ 
    &cube_geom, 
    &tetra_geom, 
    &pyramids_geom, 
    &cylinder_geom, 
    &spring_geom 
}; 
 
extern BOOL bOgl11; 
extern BOOL bCheckerOn; 
 
/******************************Public*Routine******************************\ 
* 
* InitVlen 
* 
* Precomputes scaling factor for spherical projection 
* 
\**************************************************************************/ 
 
void InitVlen(GEOMETRY *geom, int npts, PT3 *pts) 
{ 
    FLT d; 
    FLT *vl; 
 
    vl = vlen; 
    while (npts-- > 0) 
    { 
        d = V3Len(pts); 
 
        // Don't allow really close points because this leads to 
        // numeric instability and really large objects 
        assert(d > 0.01f); 
 
        // Geometries are created with size one, filling the area 
        // from -.5 to .5.  This leads to distances generally less 
        // than one, which leaves off half of the interesting morphing 
        // effects due to the projection 
        // Scaling up the scaling factor allows the values to 
        // be both above and below one 
        d *= geom->init_sf; 
         
        assert(d > 0.0001f); 
         
        *vl++ = (RADIUS-d)/d; 
         
#if 0 
        dprintf(("Distance is %f, vl %f\n", d, *(vl-1))); 
#endif 
 
        pts++; 
    } 
} 
 
/******************************Public*Routine******************************\ 
* 
* MapToSide 
* 
* Takes x,y coordinates in the range 0-1 and maps them onto the given 
* side plane for the current geometry 
* 
\**************************************************************************/ 
 
void MapToSide(PLANE_MAP *map, FLT x, FLT y, PT3 *pt) 
{ 
    pt->x = x*map->x_axis.x+y*map->y_axis.x+map->base.x; 
    pt->y = x*map->x_axis.y+y*map->y_axis.y+map->base.y; 
    pt->z = x*map->x_axis.z+y*map->y_axis.z+map->base.z; 
 
} 
 
void DrawWithVArrays (GEOMETRY *geom) 
{ 
    int side, ss, idc, k; 
    unsigned int *idx; 
 
    for (side = 0; side < geom->nsides; side++) { 
        geom->sides[side].dBuf = (GLuint *) LocalAlloc (LMEM_FIXED,  
                                                        sizeof (GLuint) *  
                                                        3 * MAXPTS * 2); 
        k = 0; 
        idx = geom->sides[side].strip_index; 
        for (ss = 0; ss < geom->sides[side].nstrips; ss++) { 
            if (geom->sides[side].strip_size[ss] < 3) continue; 
            for (idc = 2; idc < geom->sides[side].strip_size[ss]; idc++) { 
                if (!(idc % 2)) { //even 
                    geom->sides[side].dBuf[k++] = *(idx+idc-2);  
                    geom->sides[side].dBuf[k++] = *(idx+idc-1);  
                } else { 
                    geom->sides[side].dBuf[k++] = *(idx+idc-1);  
                    geom->sides[side].dBuf[k++] = *(idx+idc-2);  
                } 
                geom->sides[side].dBuf[k++] = *(idx+idc);  
            } 
            idx += geom->sides[side].strip_size[ss]; 
        } 
        geom->sides[side].num_eles = k; 
    } 
    glNormalPointer (GL_FLOAT, sizeof (PT3),  
                     (GLfloat *)&(geom->normals[0].x)); 
    glVertexPointer (3, GL_FLOAT, sizeof (PT3),  
                     (GLfloat *)&(geom->npts[0].x)); 
    glEnableClientState (GL_VERTEX_ARRAY); 
    glEnableClientState (GL_NORMAL_ARRAY); 
     
    glDisableClientState (GL_COLOR_ARRAY); 
    glDisableClientState (GL_INDEX_ARRAY); 
    glDisableClientState (GL_EDGE_FLAG_ARRAY); 
    glDisableClientState (GL_TEXTURE_COORD_ARRAY); 
} 
/******************************Public*Routine******************************\ 
* 
* InitCube 
* 
* Initialize the cube's geometry 
* 
\**************************************************************************/ 
 
#define CUBE_SIDES 6 
 
PLANE_MAP cube_planes[CUBE_SIDES] = 
{ 
    -0.5f, -0.5f,  0.5f,  1.0f,  0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 
     0.5f, -0.5f, -0.5f, -1.0f,  0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 
     0.5f,  0.5f, -0.5f, -1.0f,  0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 
    -0.5f, -0.5f, -0.5f,  1.0f,  0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 
     0.5f, -0.5f, -0.5f,  0.0f,  1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 
    -0.5f,  0.5f, -0.5f,  0.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f 
}; 
 
#define CUBE_IDX(side, x, y) ((side)*side_pts+(x)*(config.subdiv+1)+(y)) 
 
void InitCube(GEOMETRY *geom) 
{ 
    int side, x, y; 
    PT3 *pt; 
    unsigned int *sz, *idx; 
    int side_pts; 
 
    side_pts = (config.subdiv+1)*(config.subdiv+1); 
     
    geom->nsides = CUBE_SIDES; 
    geom->pts = &pts[0]; 
    geom->npts = &npts[0]; 
    geom->normals = &normals[0]; 
 
    geom->min_sf = -1.1f; 
    geom->max_sf = 5.1f; 
    geom->sf_inc = 0.05f; 
    geom->init_sf = 2.0f; 
     
    // Generate triangle strip data 
    sz = &strip_size[0]; 
    idx = &index[0]; 
    for (side = 0; side < geom->nsides; side++) 
    { 
        geom->sides[side].nstrips = config.subdiv; 
        geom->sides[side].strip_size = sz; 
        geom->sides[side].strip_index = idx; 
         
        for (x = 0; x < config.subdiv; x++) 
        { 
            *sz++ = (config.subdiv+1)*2; 
 
            for (y = 0; y < config.subdiv+1; y++) 
            { 
                *idx++ = CUBE_IDX(side, x, y); 
                *idx++ = CUBE_IDX(side, x+1, y); 
            } 
        } 
    } 
 
    assert(sz-strip_size <= DIMA(strip_size)); 
    assert(idx-index <= DIMA(index)); 
 
    
    // Generate base vertices 
    pt = geom->pts; 
    for (side = 0; side < geom->nsides; side++) 
    { 
        for (x = 0; x < config.subdiv+1; x++) 
        { 
            for (y = 0; y < config.subdiv+1; y++) 
            { 
                MapToSide(&cube_planes[side], 
                          (FLT)x/config.subdiv, (FLT)y/config.subdiv, 
                          pt); 
                pt++; 
            } 
        } 
    } 
 
    assert(pt-pts <= DIMA(pts)); 
 
    geom->total_pts = geom->nsides*side_pts; 
} 
 
 
/******************************Public*Routine******************************\ 
* 
* InitTetra 
* 
* Initialize the tetrahedron's geometry 
* 
\**************************************************************************/ 
 
#define TETRA_SIDES 4 
 
#define SQRT3 1.73205f 
#define SQRT3_2 (SQRT3/2.0f) 
#define SQRT3_3 (SQRT3/3.0f) 
#define SQRT3_6 (SQRT3/6.0f) 
#define SQRT3_12 (SQRT3/12.0f) 
 
#define TETRA_BASE (-SQRT3/8.0f) 
 
PLANE_MAP tetra_planes[TETRA_SIDES] = 
{ 
    -0.5f, TETRA_BASE, SQRT3_6, 
    1.0f, 0.0f, 0.0f, 0.0f, SQRT3_2, -SQRT3_6, 
     
    0.0f, TETRA_BASE, -SQRT3_3, 
    -0.5f, 0.0f, SQRT3_2, 0.25f, SQRT3_2, SQRT3_12, 
     
    0.5f, TETRA_BASE, SQRT3_6, 
    -0.5f, 0.0f, -SQRT3_2, -0.25f, SQRT3_2, SQRT3_12, 
     
    0.5f, TETRA_BASE, SQRT3_6, 
    -1.0f, 0.0f, 0.0f, 0.0f, 0.0f, -SQRT3_2 
}; 
 
void InitTetra(GEOMETRY *geom) 
{ 
    int side, x, y; 
    PT3 *pt; 
    int *sz, *idx; 
    int side_pts; 
    int base_pt; 
    int row_pts; 
    FLT fx; 
 
    side_pts = (config.subdiv+2)*(config.subdiv+1)/2; 
     
    geom->nsides = TETRA_SIDES; 
    geom->pts = &pts[0]; 
    geom->npts = &npts[0]; 
    geom->normals = &normals[0]; 
 
    geom->min_sf = -1.1f; 
    geom->max_sf = 5.2f; 
    geom->sf_inc = 0.05f; 
    geom->init_sf = 3.75f; 
 
    // Generate triangle strip data 
    sz = &strip_size[0]; 
    idx = &index[0]; 
    base_pt = 0; 
    for (side = 0; side < geom->nsides; side++) 
    { 
        geom->sides[side].nstrips = config.subdiv; 
        geom->sides[side].strip_size = sz; 
        geom->sides[side].strip_index = idx; 
 
        for (x = 0; x < config.subdiv; x++) 
        { 
            row_pts = config.subdiv-x+1; 
            *sz++ = row_pts*2-1; 
 
            *idx++ = base_pt; 
            for (y = 0; y < row_pts-1; y++) 
            { 
                *idx++ = base_pt+row_pts+y; 
                *idx++ = base_pt+1+y; 
            } 
 
            base_pt += row_pts; 
        } 
 
        base_pt++; 
    } 
 
    assert(sz-strip_size <= DIMA(strip_size)); 
    assert(idx-index <= DIMA(index)); 
 
    // Generate base vertices 
    pt = geom->pts; 
    for (side = 0; side < geom->nsides; side++) 
    { 
        for (x = 0; x < config.subdiv+1; x++) 
        { 
            fx = (FLT)x/config.subdiv; 
            for (y = 0; y < config.subdiv-x+1; y++) 
            { 
                MapToSide(&tetra_planes[side], 
                          fx+(FLT)y/(config.subdiv*2), 
                          (FLT)y/config.subdiv, 
                          pt); 
                pt++; 
            } 
        } 
    } 
 
    assert(pt-pts <= DIMA(pts)); 
 
    geom->total_pts = geom->nsides*side_pts; 
} 
 
/******************************Public*Routine******************************\ 
* 
* InitPyramids 
* 
* Initializes double pyramid geometry 
* 
\**************************************************************************/ 
 
#define PYRAMIDS_SIDES 8 
 
PLANE_MAP pyramids_planes[PYRAMIDS_SIDES] = 
{ 
    -0.5f, 0.0f,  0.5f,  1.0f, 0.0f,  0.0f,  0.0f,  0.5f, -0.5f, 
     0.5f, 0.0f,  0.5f, -1.0f, 0.0f,  0.0f,  0.0f, -0.5f, -0.5f, 
     0.5f, 0.0f,  0.5f,  0.0f, 0.0f, -1.0f, -0.5f,  0.5f,  0.0f, 
     0.5f, 0.0f, -0.5f,  0.0f, 0.0f,  1.0f, -0.5f, -0.5f,  0.0f, 
     0.5f, 0.0f, -0.5f, -1.0f, 0.0f,  0.0f,  0.0f,  0.5f,  0.5f, 
    -0.5f, 0.0f, -0.5f,  1.0f, 0.0f,  0.0f,  0.0f, -0.5f,  0.5f, 
    -0.5f, 0.0f, -0.5f,  0.0f, 0.0f,  1.0f,  0.5f,  0.5f,  0.0f, 
    -0.5f, 0.0f,  0.5f,  0.0f, 0.0f, -1.0f,  0.5f, -0.5f,  0.0f 
}; 
 
void InitPyramids(GEOMETRY *geom) 
{ 
    int side, x, y; 
    PT3 *pt; 
    int *sz, *idx; 
    int side_pts; 
    int base_pt; 
    int row_pts; 
    FLT fx; 
 
    side_pts = (config.subdiv+2)*(config.subdiv+1)/2; 
     
    geom->nsides = PYRAMIDS_SIDES; 
    geom->pts = &pts[0]; 
    geom->npts = &npts[0]; 
    geom->normals = &normals[0]; 
 
    geom->min_sf = -1.1f; 
    geom->max_sf = 5.2f; 
    geom->sf_inc = 0.05f; 
    geom->init_sf = 3.0f; 
 
    // Generate triangle strip data 
    sz = &strip_size[0]; 
    idx = &index[0]; 
    base_pt = 0; 
 
    for (side = 0; side < geom->nsides; side++) { 
        geom->sides[side].nstrips = config.subdiv; 
        geom->sides[side].strip_size = sz; 
        geom->sides[side].strip_index = idx; 
             
        for (x = 0; x < config.subdiv; x++) { 
            row_pts = config.subdiv-x+1; 
            *sz++ = row_pts*2-1; 
 
            *idx++ = base_pt; 
            for (y = 0; y < row_pts-1; y++) { 
                *idx++ = base_pt+row_pts+y; 
                *idx++ = base_pt+1+y; 
            } 
 
            base_pt += row_pts; 
        } 
             
        base_pt++; 
    } 
 
    assert(sz-strip_size <= DIMA(strip_size)); 
    assert(idx-index <= DIMA(index)); 
 
    // Generate base vertices 
    pt = geom->pts; 
    for (side = 0; side < geom->nsides; side++) 
    { 
        for (x = 0; x < config.subdiv+1; x++) 
        { 
            fx = (FLT)x/config.subdiv; 
            for (y = 0; y < config.subdiv-x+1; y++) 
            { 
                MapToSide(&pyramids_planes[side], 
                          fx+(FLT)y/(config.subdiv*2), 
                          (FLT)y/config.subdiv, 
                          pt); 
                pt++; 
            } 
        } 
    } 
 
    assert(pt-pts <= DIMA(pts)); 
 
    geom->total_pts = geom->nsides*side_pts; 
} 
 
/******************************Public*Routine******************************\ 
* 
* InitCylinder 
* 
* Initializes the cylinder geometry 
* 
\**************************************************************************/ 
 
void InitCylinder(GEOMETRY *geom) 
{ 
    int side, x, y; 
    PT3 *pt; 
    int *sz, *idx; 
    int base_pt; 
    int row_pts; 
    FLT fx, fz; 
    double ang; 
 
    geom->nsides = 1; 
    geom->pts = &pts[0]; 
    geom->npts = &npts[0]; 
    geom->normals = &normals[0]; 
 
    geom->min_sf = -2.5f; 
    geom->max_sf = 8.5f; 
    geom->sf_inc = 0.05f; 
    geom->init_sf = 2.1f; 
 
    // Generate triangle strip data 
    // If version 1.1 then allocate the index buffer for glDrawElements 
    sz = &strip_size[0]; 
    idx = &index[0]; 
    side = 0; 
    geom->sides[side].nstrips = config.subdiv; 
    geom->sides[side].strip_size = sz; 
    geom->sides[side].strip_index = idx; 
     
    row_pts = config.subdiv+1; 
    base_pt = 0; 
    for (x = 0; x < config.subdiv; x++) { 
        *sz++ = row_pts*2; 
         
        for (y = 0; y < row_pts; y++) { 
            // Wrap around at the edge so the cylinder normals 
            // are properly averaged 
            if (x == config.subdiv-1) { 
                *idx++ = y; 
            } 
            else { 
                *idx++ = base_pt+row_pts+y; 
            } 
            *idx++ = base_pt+y; 
        } 
 
        base_pt += row_pts; 
    } 
     
    assert(sz-strip_size <= DIMA(strip_size)); 
    assert(idx-index <= DIMA(index)); 
 
    // Generate base vertices 
    pt = geom->pts; 
    ang = 0; 
    for (x = 0; x < config.subdiv; x++) 
    { 
        fx = (FLT)cos(ang)*0.5f; 
        fz = (FLT)sin(ang)*0.5f; 
        for (y = 0; y < config.subdiv+1; y++) 
        { 
            pt->x = fx; 
            pt->y = (FLT)y/config.subdiv-0.5f; 
            pt->z = fz; 
            pt++; 
        } 
        ang += (2*PI)/config.subdiv; 
    } 
 
    assert(pt-pts <= DIMA(pts)); 
 
    geom->total_pts = geom->nsides*(config.subdiv+1)*config.subdiv; 
} 
 
/******************************Public*Routine******************************\ 
* 
* InitSpring 
* 
* Initializes the spring geometry 
* 
\**************************************************************************/ 
 
#define SPRING_RADIUS 0.1f 
#define SPRING_CENTER (0.5f-SPRING_RADIUS) 
 
void InitSpring(GEOMETRY *geom) 
{ 
    int side, x, y; 
    PT3 *pt; 
    int *sz, *idx; 
    double ang_center, ang_surf; 
    FLT cs, sn; 
    FLT rad; 
    PLANE_MAP plane; 
    int spin_pts; 
    int row_pts; 
 
    geom->nsides = 1; 
    geom->pts = &pts[0]; 
    geom->npts = &npts[0]; 
    geom->normals = &normals[0]; 
 
    geom->min_sf = -2.2f; 
    geom->max_sf = 0.2f; 
    geom->sf_inc = 0.05f; 
    geom->init_sf = 1.0f; 
 
    // Generate triangle strip data 
    // If version 1.1 then allocate the index buffer for glDrawElements 
    sz = &strip_size[0]; 
    idx = &index[0]; 
    side = 0; 
    geom->sides[side].nstrips = config.subdiv; 
    geom->sides[side].strip_size = sz; 
    geom->sides[side].strip_index = idx; 
     
    row_pts = config.subdiv; 
    spin_pts = 4*config.subdiv+1; 
    for (x = 0; x < config.subdiv; x++) { 
        *sz++ = spin_pts*2; 
 
        for (y = 0; y < spin_pts; y++) { 
            *idx++ = x+row_pts*y; 
            // Wrap around at the edge so the cylindrical surface 
            // of the tube is seamless.  Without this the normal 
            // averaging would be incorrect and a seam would be visible 
            if (x == config.subdiv-1) { 
                *idx++ = row_pts*y; 
            } 
            else { 
                *idx++ = x+row_pts*y+1; 
            } 
        } 
    } 
     
    assert(sz-strip_size <= DIMA(strip_size)); 
    assert(idx-index <= DIMA(index)); 
 
    // Generate base vertices 
    pt = geom->pts; 
    ang_center = 0; 
    plane.y_axis.x = 0.0f; 
    plane.y_axis.y = SPRING_RADIUS; 
    plane.y_axis.z = 0.0f; 
    plane.x_axis.y = 0.0f; 
    for (x = 0; x < spin_pts; x++) 
    { 
        cs = (FLT)cos(ang_center); 
        sn = (FLT)sin(ang_center); 
        rad = 0.5f-(FLT)x/(spin_pts-1)*(SPRING_CENTER/2); 
        plane.base.x = cs*rad; 
        plane.base.y = -0.5f+(FLT)x/(spin_pts-1); 
        plane.base.z = sn*rad; 
        plane.x_axis.x = cs*SPRING_RADIUS; 
        plane.x_axis.z = sn*SPRING_RADIUS; 
 
        ang_surf = 0; 
        for (y = 0; y < config.subdiv; y++) 
        { 
            MapToSide(&plane, 
                      (FLT)cos(ang_surf), (FLT)sin(ang_surf), 
                      pt); 
            pt++; 
            ang_surf += (2*PI)/config.subdiv; 
        } 
         
        ang_center += (4*PI)/(spin_pts-1); 
    } 
 
    assert(pt-pts <= DIMA(pts)); 
 
    geom->total_pts = geom->nsides*spin_pts*config.subdiv; 
} 
 
/******************************Public*Routine******************************\ 
* 
* DrawGeom 
* 
* Draw the current geometry 
* 
\**************************************************************************/ 
 
void DrawGeom(GEOMETRY *geom) 
{ 
    int side, strip, k; 
    int *idx, idc, idxv; 
 
    if (config.smooth_colors) 
    { 
        glShadeModel(GL_SMOOTH); 
    } 
    else 
    { 
        glShadeModel(GL_FLAT); 
    } 
 
    if (config.color_pick == ID_COL_SINGLE) 
    { 
        glMaterialfv(config.two_sided, GL_DIFFUSE, solid_cols); 
    } 
 
    if (!(bOgl11 && !bCheckerOn)) { 
        for (side = 0; side < geom->nsides; side++) { 
            if (config.color_pick == ID_COL_PER_SIDE) { 
                glMaterialfv(config.two_sided, GL_DIFFUSE, side_cols[side]); 
            } 
 
            idx = geom->sides[side].strip_index; 
            for (strip = 0; strip < geom->sides[side].nstrips; strip++) { 
                glBegin(GL_TRIANGLE_STRIP); 
 
                for (idc = 0; idc < geom->sides[side].strip_size[strip];  
                     idc++) { 
                    idxv = *idx++; 
 
                    assert(idxv >=0 && idxv < geom->total_pts); 
 
                    if (config.color_pick == ID_COL_CHECKER) { 
                        if (config.triangle_colors) { 
                            glMaterialfv(config.two_sided, GL_DIFFUSE, 
                                         checker_cols[side][(idc+1)/2+strip & 
                                                           1]); 
                        } 
                        else { 
                            glMaterialfv(config.two_sided, GL_DIFFUSE, 
                                         checker_cols[side][idc/2+strip & 1]); 
                        } 
                    } 
                 
                    glNormal3fv((GLfloat *)&geom->normals[idxv]); 
                    glVertex3fv((GLfloat *)&geom->npts[idxv]); 
                } 
 
                glEnd(); 
            } 
        } 
    } else { 
        k = 0; 
        for (side = 0; side < geom->nsides; side++) { 
            if (config.color_pick == ID_COL_PER_SIDE) 
                glMaterialfv(config.two_sided, GL_DIFFUSE, side_cols[side]); 
 
            glDrawElements (GL_TRIANGLES, geom->sides[side].num_eles,  
                            GL_UNSIGNED_INT, &(geom->sides[side].dBuf[0])); 
            k += geom->sides[side].num_eles; 
        } 
    }     
} 
 
/******************************Public*Routine******************************\ 
* 
* ComputeAveragedNormals 
* 
* Compute face-averaged normals for each vertex 
* 
\**************************************************************************/ 
 
void ComputeAveragedNormals(GEOMETRY *geom) 
{ 
    int side, strip; 
    int *sz; 
    int *idx, idx1, idx2, idx3; 
    int tc, idc; 
    PT3 v1, v2, n1; 
     
    memset(geom->normals, 0, sizeof(PT3)*geom->total_pts); 
     
    for (side = 0; side < geom->nsides; side++) 
    { 
        idx = geom->sides[side].strip_index; 
        sz = geom->sides[side].strip_size; 
        for (strip = 0; strip < geom->sides[side].nstrips; strip++) 
        { 
            idx1 = *idx++; 
            idx2 = *idx++; 
 
            assert(idx1 >= 0 && idx1 < geom->total_pts && 
                   idx2 >= 0 && idx2 < geom->total_pts); 
             
            tc = (*sz++)-2; 
            for (idc = 0; idc < tc; idc++) 
            { 
                idx3 = *idx++; 
 
                assert(idx3 >= 0 && idx3 < geom->total_pts); 
                 
                V3Sub(&geom->npts[idx3], &geom->npts[idx1], &v1); 
                V3Sub(&geom->npts[idx2], &geom->npts[idx1], &v2); 
                V3Cross(&v1, &v2, &n1); 
                // Triangle strip ordering causes half of the triangles 
                // to be oriented oppositely from the others 
                // Those triangles need to have their normals flipped 
                // so the whole triangle strip has consistent normals 
                if ((idc & 1) == 0) 
                { 
                    n1.x = -n1.x; 
                    n1.y = -n1.y; 
                    n1.z = -n1.z; 
                } 
                 
                V3Add(&geom->normals[idx1], &n1, &geom->normals[idx1]); 
                V3Add(&geom->normals[idx2], &n1, &geom->normals[idx2]); 
                V3Add(&geom->normals[idx3], &n1, &geom->normals[idx3]); 
 
                idx1 = idx2; 
                idx2 = idx3; 
            } 
        } 
    } 
} 
 
/******************************Public*Routine******************************\ 
* 
* UpdatePts 
* 
* Project the point array through a sphere according to the given scale factor 
* 
\**************************************************************************/ 
 
void UpdatePts(GEOMETRY *geom, FLT sf) 
{ 
    int pt; 
    FLT f, *vl; 
    PT3 *v; 
    PT3 *p; 
 
    vl = vlen; 
    p = &geom->pts[0]; 
    v = &geom->npts[0]; 
    for (pt = 0; pt < geom->total_pts; pt++) 
    { 
        f = (*vl++)*sf+1; 
        v->x = p->x*f; 
        v->y = p->y*f; 
        v->z = p->z*f; 
        p++; 
        v++; 
    } 
 
    ComputeAveragedNormals(geom); 
}