/*
======================================================================
libply.c

Alexander Enzmann  01 Nov 94
Ernie Wright  01 Jun 05

A library of routines for generating polygonal object primitives.
Modified for use in a LightWave Modeler plug-in.  The public routines
for object primitives are

   lib_output_polygon_cylcone()
   lib_output_polygon_disc()
   lib_output_polygon_sphere()
   lib_output_polygon_height()
   lib_output_polygon_torus()
   lib_output_polygon_box()

All of these rely on one or both of the following to actually generate
polygons,

   lib_output_polygon()
   lib_output_polypatch()

so those are the only places where calls to LightWave Modeler geometry
functions are made.
====================================================================== */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "lwmod.h"
#include "lib.h"


/*-----------------------------------------------------------------*/
/* defines/constants section */

#define VERT(i, a) ((verts[gPoly_vbuffer[i]])[a])


/*-----------------------------------------------------------------*/
/* Keep track of how many vertices/faces have been emitted */
unsigned long gVertex_count = 0; /* Vertex coordinates */
unsigned long gNormal_count = 0; /* Vertex normals */
unsigned long gFace_count = 0;

/* Storage for polygon indices, allocated in libini.c */
unsigned int *gPoly_vbuffer = NULL;
int *gPoly_end = NULL;

/* Globals to determine which axes can be used to split the polygon */
int gPoly_Axis1 = 0;
int gPoly_Axis2 = 1;


/*
======================================================================
lib_output_polygon_cylcone()

Create polygonal cylinders and cones.  Polygons are output using
lib_output_polypatch().
====================================================================== */

void lib_output_polygon_cylcone( COORD4 base_pt, COORD4 apex_pt )
{
   double angle, delta_angle, height, divisor, ba;
   COORD3 axis, dir, norm_axis, start_dir, start_norm;
   COORD3 norm[4], vert[4], start_radius[4];
   MATRIX nmx, mx;
   int    i;

   SUB3_COORD3(axis, apex_pt, base_pt);
   COPY_COORD3(norm_axis, axis);
   height = lib_normalize_vector(norm_axis);

   SET_COORD3(dir, 0.0, 0.0, 1.0);
   CROSS(start_dir, axis, dir);
   divisor = lib_normalize_vector(start_dir);

   if (ABS(divisor) < EPSILON2) {
      SET_COORD3(dir, 1.0, 0.0, 0.0);
      CROSS(start_dir, axis, dir);
      lib_normalize_vector(start_dir);
   }

   start_radius[0][X] = start_dir[X] * base_pt[W];
   start_radius[0][Y] = start_dir[Y] * base_pt[W];
   start_radius[0][Z] = start_dir[Z] * base_pt[W];
   ADD3_COORD3(vert[2], base_pt, start_radius[0]);

   start_radius[1][X] = start_dir[X] * apex_pt[W];
   start_radius[1][Y] = start_dir[Y] * apex_pt[W];
   start_radius[1][Z] = start_dir[Z] * apex_pt[W];
   ADD3_COORD3(vert[1], apex_pt, start_radius[1]);

   if ( base_pt[W] == apex_pt[W] ) {
      /* it's a cylinder, so simply copy dir to norm */
      COPY_COORD3( start_norm, start_dir ) ;
   } else {
      /* it's a cone, so compute true normal here */
      ba = base_pt[W] - apex_pt[W] ;
      start_norm[X] = start_dir[X] * height + norm_axis[X] * ba;
      start_norm[Y] = start_dir[Y] * height + norm_axis[Y] * ba;
      start_norm[Z] = start_dir[Z] * height + norm_axis[Z] * ba;
      lib_normalize_vector(start_norm);
   }
   COPY_COORD3(norm[2], start_norm);
   COPY_COORD3(norm[1], start_norm);

   delta_angle = 2.0 * PI / (double)(4*gU_resolution);
   for (i=1,angle=delta_angle;i<=4*gU_resolution;++i,angle+=delta_angle) {
      lib_create_axis_rotate_matrix(mx, norm_axis, angle);
      lib_invert_matrix(nmx, mx);
      lib_transform_point(vert[0], start_radius[1], mx);
      ADD2_COORD3(vert[0], apex_pt);
      lib_transform_normal(norm[0], start_norm, nmx);
      lib_output_polypatch(3, vert, norm);
      COPY_COORD3(vert[1], vert[0]);
      COPY_COORD3(norm[1], norm[0]);
      lib_transform_point(vert[0], start_radius[0], mx);
      ADD2_COORD3(vert[0], base_pt);
      lib_output_polypatch(3, vert, norm);

      COPY_COORD3(vert[2], vert[0]);
      COPY_COORD3(norm[2], norm[0]);
   }
}


static void disc_evaluator( MATRIX trans, double theta, double v, double r,
   COORD3 vert )
{
   COORD3 tvert;

   /* Compute the position of the point */
   SET_COORD3(tvert, (r + v) * cos(theta), (r + v) * sin(theta), 0.0);
   lib_transform_point(vert, tvert, trans);
}


/*
======================================================================
lib_output_polygon_disc()

Create a disk (a circle with a hole in the center).  Polygons are
output using lib_output_polygon().
====================================================================== */

void lib_output_polygon_disc( COORD3 center, COORD3 normal, double iradius,
   double oradius )
{
   double u, v, delta_u, delta_v;
   MATRIX mx, imx;
   int i;
   COORD3 norm, vert[4];

   COPY_COORD3(norm, normal);
   if ( lib_normalize_vector(norm) < EPSILON2) {
      /* bad disc normal */
      return;
   }
   lib_create_canonical_matrix(mx, imx, center, norm);
   delta_u = 2.0 * PI / (double)(4 * gU_resolution);

   /* Dump out polygons */
   for (i=0,u=0.0;i<4*gU_resolution;i++,u+=delta_u) {
      v = 0.0;
      delta_v = oradius-iradius;
      disc_evaluator(imx, u, v, iradius, vert[3]);
      disc_evaluator(imx, u+delta_u, v, iradius, vert[2]);
      disc_evaluator(imx, u+delta_u, v+delta_v, iradius, vert[1]);
      disc_evaluator(imx, u, v+delta_v, iradius, vert[0]);
      lib_output_polygon(4, vert);
   }
}


/*
======================================================================
lib_output_polygon_sphere()

Create a sphere.  Polygons are output using lib_output_polypatch().
====================================================================== */

void lib_output_polygon_sphere( COORD4 center_pt )
{
   double  angle;
   COORD3  edge_norm[3], edge_pt[3];
   long    num_face, num_edge, num_tri, num_vert;
   COORD3  *x_axis, *y_axis, **pt;
   COORD3  mid_axis;
   MATRIX  rot_mx;
   long    u_pol, v_pol;

   /* Allocate storage for the polygon vertices */
   x_axis = (COORD3 *)malloc((gU_resolution+1) * sizeof(COORD3));
   y_axis = (COORD3 *)malloc((gV_resolution+1) * sizeof(COORD3));
   pt     = (COORD3 **)malloc((gU_resolution+1) * sizeof(COORD3 *));
   if (x_axis == NULL || y_axis == NULL || pt == NULL) {
      /* memory allocation failed */
      return;
    }

   for (num_edge=0;num_edge<gU_resolution+1;num_edge++) {
      pt[num_edge] = (COORD3 *)malloc((gV_resolution+1) * sizeof(COORD3));
      if (pt[num_edge] == NULL) {
         /* memory allocation failed */
         return;
      }
   }

   /* calculate axes used to find grid points */
   for (num_edge=0;num_edge<=gU_resolution;++num_edge) {
      angle = (PI/4.0) * (2.0*(double)num_edge/gU_resolution - 1.0);
      mid_axis[X] = 1.0; mid_axis[Y] = 0.0; mid_axis[Z] = 0.0;
      lib_create_rotate_matrix(rot_mx, Y_AXIS, angle);
      lib_transform_vector(x_axis[num_edge], mid_axis, rot_mx);
   }

   for (num_edge=0;num_edge<=gV_resolution;++num_edge) {
      angle = (PI/4.0) * (2.0*(double)num_edge/gV_resolution - 1.0);
      mid_axis[X] = 0.0; mid_axis[Y] = 1.0; mid_axis[Z] = 0.0;
      lib_create_rotate_matrix(rot_mx, X_AXIS, angle);
      lib_transform_vector(y_axis[num_edge], mid_axis, rot_mx);
   }

   /* set up grid of points on +Z sphere surface */
   for (u_pol=0;u_pol<=gU_resolution;++u_pol) {
      for (v_pol=0;v_pol<=gU_resolution;++v_pol) {
         CROSS(pt[u_pol][v_pol], x_axis[u_pol], y_axis[v_pol]);
         lib_normalize_vector(pt[u_pol][v_pol]);
      }
   }

   for (num_face=0;num_face<6;++num_face) {
      /* transform points to cube face */
      for (u_pol=0;u_pol<=gU_resolution;++u_pol) {
         for (v_pol=0;v_pol<=gV_resolution;++v_pol) {
            lib_rotate_cube_face(pt[u_pol][v_pol], Z_AXIS, num_face);
         }
      }

      /* output grid */
      for (u_pol=0;u_pol<gU_resolution;++u_pol) {
         for (v_pol=0;v_pol<gV_resolution;++v_pol) {
            for (num_tri=0;num_tri<2;++num_tri) {
               for (num_edge=0;num_edge<3;++num_edge) {
                  num_vert = (num_tri*2 + num_edge) % 4;
                  if (num_vert == 0) {
                     COPY_COORD3(edge_pt[num_edge], pt[u_pol][v_pol]);
                  } else if ( num_vert == 1 ) {
                     COPY_COORD3(edge_pt[num_edge], pt[u_pol][v_pol+1]);
                  } else if ( num_vert == 2 ) {
                     COPY_COORD3(edge_pt[num_edge],pt[u_pol+1][v_pol+1]);
                  } else {
                     COPY_COORD3(edge_pt[num_edge], pt[u_pol+1][v_pol]);
                  }
                  COPY_COORD3(edge_norm[num_edge], edge_pt[num_edge]);
                  edge_pt[num_edge][X] =
                     edge_pt[num_edge][X] * center_pt[W] +
                     center_pt[X];
                  edge_pt[num_edge][Y] =
                     edge_pt[num_edge][Y] * center_pt[W] +
                     center_pt[Y];
                  edge_pt[num_edge][Z] =
                     edge_pt[num_edge][Z] * center_pt[W] +
                     center_pt[Z];

               }
               lib_output_polypatch(3, edge_pt, edge_norm);
            }
         }
      }
   }

   /* Release any memory used */
   for (num_edge=0;num_edge<gU_resolution+1;num_edge++)
      free(pt[num_edge]);
   free(pt);
   free(y_axis);
   free(x_axis);
}


/*
======================================================================
lib_output_polygon_height()

Create height field geometry.  Data for the height field is in an
array of pointers to float, treated as a 2D array of float.  Polygons
are output using lib_output_polygon().
====================================================================== */

void lib_output_polygon_height( int height, int width, float **data,
   double x0, double x1, double y0, double y1, double z0, double z1 )
{
   int i, j;
   double xdelta, zdelta;
   COORD3 verts[3];

   xdelta = (x1 - x0) / (double)(width - 1);
   zdelta = (z1 - z0) / (double)(height - 1);
   for (i=0;i<height-1;i++) {
      for (j=0;j<width-1;j++) {
         SET_COORD3(verts[0], x0 + j * xdelta, y0 + data[i][j],
            z0 + i * zdelta);
         SET_COORD3(verts[1], x0 + (j+1) * xdelta, y0 + data[i+1][j+1],
            z0 + (i + 1) * zdelta);
         SET_COORD3(verts[2], x0 + (j+1) * xdelta, y0 + data[i][j+1],
            z0 + i * zdelta);
         lib_output_polygon(3, verts);
         COPY_COORD3(verts[2], verts[1]); /* copy corner from previous */
         SET_COORD3(verts[1], x0 + j * xdelta, y0 + data[i+1][j],
            z0 + (i + 1) * zdelta);
         lib_output_polygon(3, verts);
      }
   }
}


static void torus_evaluator( MATRIX trans, double theta, double phi,
   double r0, double r1, COORD3 vert, COORD3 norm )
{
   COORD3 v0, v1, tvert, tnorm;

   /* Compute the position of the point */
   SET_COORD3(tvert, (r0 + r1 * sin(theta)) * cos(phi),
      (r0 + r1 * sin(theta)) * sin(phi),
      r1 * cos(theta));
   /* Compute the normal at that point */
   SET_COORD3(v0, r1*cos(theta)*cos(phi),
      r1*cos(theta)*sin(phi),
      -r1*sin(theta));
   SET_COORD3(v1,-(r0+r1*sin(theta))*sin(phi),
      (r0+r1*sin(theta))*cos(phi),
      0.0);
   CROSS(tnorm, v0, v1);
   lib_normalize_vector(tnorm);
   lib_transform_point(vert, tvert, trans);
   lib_transform_vector(norm, tnorm, trans);
}


/*
======================================================================
lib_output_polygon_torus()

Create a torus.  Polygons are output using lib_output_polypatch().
====================================================================== */

void lib_output_polygon_torus( COORD3 center, COORD3 normal, double iradius,
   double oradius )
{
   double u, v, delta_u, delta_v;
   MATRIX mx, imx;
   int i, j;
   COORD3 vert[4], norm[4];

   if ( lib_normalize_vector(normal) < EPSILON2) {
      /* bad torus normal */
      return;
   }
   lib_create_canonical_matrix(mx, imx, center, normal);
   delta_u = 2.0 * PI / (double)(4*gU_resolution);
   delta_v = 2.0 * PI / (double)(4*gV_resolution);

   /* Dump out polygons */
   for (i=0,u=0.0;i<(4*gU_resolution);i++,u+=delta_u) {
      for (j=0,v=0.0;j<(4*gV_resolution);j++,v+=delta_v) {
         torus_evaluator(imx, u, v, iradius, oradius, vert[2], norm[2]);
         torus_evaluator(imx, u, v+delta_v, iradius, oradius,
            vert[1], norm[1]);
         torus_evaluator(imx, u+delta_u, v+delta_v,
            iradius, oradius, vert[0], norm[0]);
         lib_output_polypatch(3, vert, norm);
         COPY_COORD3(vert[1], vert[0]);
         COPY_COORD3(norm[1], norm[0]);
         torus_evaluator(imx, u+delta_u, v, iradius, oradius,
            vert[0], norm[0]);
         lib_output_polypatch(3, vert, norm);
      }
   }
}


/*
======================================================================
lib_output_polygon_box()

Create a box.  Polygons are output using lib_output_polygon().
====================================================================== */

void lib_output_polygon_box( COORD3 p1, COORD3 p2 )
{
   COORD3 box_verts[4];

   /* Sides */
   SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
   SET_COORD3(box_verts[1], p1[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
   SET_COORD3(box_verts[3], p1[X], p2[Y], p1[Z]);
   lib_output_polygon(4, box_verts);
   SET_COORD3(box_verts[0], p2[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[1], p2[X], p1[Y], p1[Z]);
   SET_COORD3(box_verts[2], p2[X], p2[Y], p1[Z]);
   SET_COORD3(box_verts[3], p2[X], p2[Y], p2[Z]);
   lib_output_polygon(4, box_verts);

   /* Front/Back */
   SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
   SET_COORD3(box_verts[3], p2[X], p1[Y], p1[Z]);
   SET_COORD3(box_verts[2], p2[X], p2[Y], p1[Z]);
   SET_COORD3(box_verts[1], p1[X], p2[Y], p1[Z]);
   lib_output_polygon(4, box_verts);
   SET_COORD3(box_verts[0], p2[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[3], p1[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
   SET_COORD3(box_verts[1], p2[X], p2[Y], p2[Z]);
   lib_output_polygon(4, box_verts);

   /* Top/Bottom */
   SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
   SET_COORD3(box_verts[3], p1[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[2], p2[X], p1[Y], p2[Z]);
   SET_COORD3(box_verts[1], p2[X], p1[Y], p1[Z]);
   lib_output_polygon(4, box_verts);
   SET_COORD3(box_verts[0], p2[X], p2[Y], p1[Z]);
   SET_COORD3(box_verts[3], p2[X], p2[Y], p2[Z]);
   SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
   SET_COORD3(box_verts[1], p1[X], p2[Y], p1[Z]);
   lib_output_polygon(4, box_verts);
}


/*
======================================================================
find_axes()

Given a polygon defined by vertices in verts, determine which of the
components of the vertex correspond to useful x and y coordinates -
with these we can pretend the polygon is 2D to do our work on it.
====================================================================== */

static void find_axes( COORD3 *verts )
{
   double P1[3], P2[3], x, y, z;

   P1[0] = VERT(1, 0) - VERT(0, 0);
   P1[1] = VERT(1, 1) - VERT(0, 1);
   P1[2] = VERT(1, 2) - VERT(0, 2);

   P2[0] = VERT(2, 0) - VERT(0, 0);
   P2[1] = VERT(2, 1) - VERT(0, 1);
   P2[2] = VERT(2, 2) - VERT(0, 2);

   /* Cross product - don't need to normalize cause we're only interested
   in the size of the components */
   x = fabs(P1[1] * P2[2] - P1[2] * P2[1]);
   y = fabs(P1[2] * P2[0] - P1[0] * P2[2]);
   z = fabs(P1[0] * P2[1] - P1[1] * P2[0]);

   if (x > y && x > z) {
      gPoly_Axis1 = 1;
      gPoly_Axis2 = 2;
   } else if (y > x && y > z) {
      gPoly_Axis1 = 0;
      gPoly_Axis2 = 2;
   } else {
      gPoly_Axis1 = 0;
      gPoly_Axis2 = 1;
   }
}


/*
======================================================================
leftmost_vertex()

Find the left most vertex in the polygon that has vertices m ... n.
====================================================================== */

static int leftmost_vertex( int m, int n, COORD3 *verts )
{
   int l, i;
   double x;

   /* Assume the first vertex is the farthest to the left */
   l = m;
   x = VERT(m, gPoly_Axis1);

   /* Now see if any of the others are farther to the left */
   for (i=m+1;i<=n;i++) {
      if (VERT(i, gPoly_Axis1) < x) {
         l = i;
         x = VERT(i, gPoly_Axis1);
      }
   }
   return l;
}


/*
======================================================================
split_vertex()

Given the leftmost vertex in a polygon, this routine finds another
vertex that can be used to safely split the polygon.
====================================================================== */

static int split_vertex( int l, int la, int lb, int m, int n, COORD3 *verts )
{
   int t, k, lpu, lpl;
   double yu, yl;

   yu = MAX(VERT(l, gPoly_Axis2), MAX(VERT(la, gPoly_Axis2), VERT(lb, gPoly_Axis2)));
   yl = MIN(VERT(l, gPoly_Axis2), MIN(VERT(la, gPoly_Axis2), VERT(lb, gPoly_Axis2)));
   if (VERT(lb, gPoly_Axis2) > VERT(la, gPoly_Axis2)) {
      lpu = lb;
      lpl = la;
   } else {
      lpu = la;
      lpl = lb;
   }
   t = (VERT(lb, gPoly_Axis1) > VERT(la, gPoly_Axis1) ? lb : la);
   for (k=m;k<n;k++) {
      if (k != la && k != l && k != lb) {
         if (VERT(k, gPoly_Axis2) <= yu && VERT(k, gPoly_Axis2) >= yl) {
            if (VERT(k, gPoly_Axis1) < VERT(t, gPoly_Axis1) &&
               ((VERT(k, gPoly_Axis2) - VERT(l, gPoly_Axis2)) *
               (VERT(lpu, gPoly_Axis1) - VERT(l, gPoly_Axis1))) <=
               ((VERT(lpu, gPoly_Axis2) - VERT(l, gPoly_Axis2)) *
               (VERT(k, gPoly_Axis1) - VERT(l, gPoly_Axis1)))) {
               if (((VERT(k, gPoly_Axis2) - VERT(l, gPoly_Axis2)) *
                  (VERT(lpl, gPoly_Axis1) - VERT(l, gPoly_Axis1))) >=
                  ((VERT(lpl, gPoly_Axis2) - VERT(l, gPoly_Axis2)) *
                  (VERT(k, gPoly_Axis1) - VERT(l, gPoly_Axis1)))) {
                  t = k;
               }
            }
         }
      }
   }
   return t;
}


/*
======================================================================
linear_vertices()

Test polygon vertices to see if they are linear.
====================================================================== */

static int linear_vertices( int m, int n, COORD3 *verts )
{
   /* Not doing anything right now */
   return 0;
}


/*
======================================================================
perform_split()

Shift vertex indices around to make two polygons out of one.
====================================================================== */

static void perform_split( int m, int m1, int n, int n1 )
{
   int i, j, k;

   k = n + 3 - m;
   /* Move the new polygon up over the place the current one sits */
   for (j=m1;j<=n1;j++) gPoly_vbuffer[j+k] = gPoly_vbuffer[j];

   /* Move top part of remaining polygon */
   for (j=n;j>=n1;j--) gPoly_vbuffer[j+2] = gPoly_vbuffer[j];

   /* Move bottom part of remaining polygon */
   k = n1 - m1 + 1;
   for (j=m1;j>=m;j--) gPoly_vbuffer[j+k] = gPoly_vbuffer[j];

   /* Copy the new polygon so that it sits before the remaining polygon */
   i = n + 3 - m;
   k = m - m1;
   for (j=m1;j<=n1;j++) gPoly_vbuffer[j+k] = gPoly_vbuffer[j+i];
}


/*
======================================================================
add_new_triangle()

Copy an indirectly referenced triangle into the output triangle
buffer.
====================================================================== */

static void add_new_triangle( int m, COORD3 *verts, COORD3 *norms,
   int *out_cnt, COORD3 **out_verts, COORD3 **out_norms )
{
   if (out_verts != NULL) {
      COPY_COORD3(out_verts[*out_cnt][0], verts[gPoly_vbuffer[m]]);
      COPY_COORD3(out_verts[*out_cnt][1], verts[gPoly_vbuffer[m+1]]);
      COPY_COORD3(out_verts[*out_cnt][2], verts[gPoly_vbuffer[m+2]]);
   }
   if (out_norms != NULL) {
      COPY_COORD3(out_norms[*out_cnt][0], norms[gPoly_vbuffer[m]]);
      COPY_COORD3(out_norms[*out_cnt][1], norms[gPoly_vbuffer[m+1]]);
      COPY_COORD3(out_norms[*out_cnt][2], norms[gPoly_vbuffer[m+2]]);
   }
   *out_cnt += 1;
}


/*
======================================================================
split_buffered_polygon()
====================================================================== */

static void split_buffered_polygon( int cnt, COORD3 *verts, COORD3 *norms,
   int *out_cnt, COORD3 **out_verts, COORD3 **out_norms )
{
   int i, m, m1, n, n1;
   int l, la, lb, ls;

   /* No triangles to start with */
   *out_cnt = 0;

   /* Initialize the polygon splitter */
   gPoly_end[0] = -1;
   gPoly_end[1] = cnt-1;

   /* Split and push polygons until they turn into triangles */
   for (i=1;i>0;) {
      m = gPoly_end[i-1] + 1;
      n = gPoly_end[i];
      if (n - m == 2) {
         if (!linear_vertices(m, n, verts)) {
            add_new_triangle(m, verts, norms, out_cnt,
               out_verts, out_norms);
         }
         i = i - 1;
      } else {
         l = leftmost_vertex(m, n, verts);
         la = (l == n ? m : l + 1);
         lb = (l == m ? n : l - 1);
         ls = split_vertex(l, la, lb, m, n, verts);
         if (ls == la || ls == lb) {
            m1 = (la < lb ? la : lb);
            n1 = (la > lb ? la : lb);
         } else {
            m1 = (l < ls ? l : ls);
            n1 = (l > ls ? l : ls);
         }
         perform_split(m, m1, n, n1);
         gPoly_end[i++] = m + n1 - m1;
         gPoly_end[i] = n + 2;
      }
   }
}


/*
======================================================================
split_polygon()

Split an arbitrary polygon into triangles.
====================================================================== */

static void split_polygon( int n, COORD3 *vert, COORD3 *norm )
{
   LWPntID *ptid;
   double lwvert[ 3 ];
   COORD3 **out_verts, **out_norms;
   MATRIX nmx, txmat;
   int i;
   int t, out_n;

   /* Can't split a NULL vertex list */
   if (vert == NULL) return;
   if (gPoly_vbuffer == NULL) { /* [esp] Added error */
      lib_storage_initialize();
      /* [are] removed error, go and initialize if it hasn't been done. */
   }

   /* Allocate space to hold the intermediate polygon stacks */
   out_verts = (COORD3 **)malloc((n - 2) * sizeof(COORD3 *));
   if (norm != NULL)
      out_norms = (COORD3 **)malloc((n - 2) * sizeof(COORD3 *));
   else
      out_norms = NULL;

   for (i=0;i<n-2;i++) {
      out_verts[i] = (COORD3 *)malloc(3 * sizeof(COORD3));
      if (norm != NULL)
         out_norms[i] = (COORD3 *)malloc(3 * sizeof(COORD3));
   }

   /* Start with a strict identity of vertices in verts and vertices in
   the polygon buffer */
   for (i=0;i<n;i++) gPoly_vbuffer[i] = i;

   /* Make sure we know which axes to look at */
   find_axes(vert);

   out_n = 0;
   split_buffered_polygon(n, vert, norm, &out_n, out_verts, out_norms);

   if (lib_tx_active()) {
      /* Perform transformations of the vertices and normals of
      the polygon(s) */
      lib_get_current_tx(txmat);
      lib_invert_matrix(nmx, txmat);
      for (t=0;t<out_n;t++)
         for (i=0;i<3;i++) {
            lib_transform_point(out_verts[t][i], out_verts[t][i], txmat);
            if (out_norms != NULL)
               lib_transform_normal(out_norms[t][i], out_norms[t][i], nmx);
         }
   }

   /* Now output the triangles that we generated */

   ptid = calloc( 3, sizeof( LWPntID ));

   for ( t = 0; t < out_n; t++ ) {
      for ( i = 0; i < 3; i++ ) {
         lwvert[ 0 ] =  out_verts[ t ][ i ][ X ];
         lwvert[ 2 ] = -out_verts[ t ][ i ][ Y ];
         lwvert[ 1 ] =  out_verts[ t ][ i ][ Z ];
         ptid[ i ] = lwAddPoint( lwvert );
      }
      lwAddFace( 3, ptid );
   }

   /* Clean up intermediate storage */
   free( ptid );
   for (i=0;i<n-2;i++) {
      free(out_verts[i]);
      if (norm != NULL)
         free(out_norms[i]);
   }

   free(out_verts);
   if (out_norms != NULL) free(out_norms);
}


/*
======================================================================
lib_output_polygon()

Emit a polygon.
====================================================================== */

void lib_output_polygon( int tot_vert, COORD3 *vert )
{
   LWPntID *ptid;
   double lwvert[ 3 ];
   int i, j;
   COORD3 x;
   MATRIX txmat;

   /* First let's do a couple of checks to see if this is a valid polygon */
   for (i=0;i<tot_vert;) {
      /* If there are two adjacent coordinates that degenerate then
         collapse them down to one */
      SUB3_COORD3(x, vert[i], vert[(i+1)%tot_vert]);
      if (lib_normalize_vector(x) < EPSILON2) {
         for (j=i;j<tot_vert-1;j++)
            memcpy(&vert[j], &vert[j+1], sizeof(COORD3));
         tot_vert--;
      }
      else {
         i++;
      }
   }

   if (tot_vert < 3)
       /* No such thing as a poly that only has two sides */
       return;

   if (lib_tx_active()) {
      /* Perform transformations of the vertices and normals of
         the polygon(s) */
      lib_get_current_tx(txmat);
      for (i=0;i<tot_vert;i++)
         lib_transform_point(vert[i], vert[i], txmat);
   }

   /* call LightWave Modeler to add the polygon */

   ptid = calloc( tot_vert, sizeof( LWPntID ));

   for ( i = 0; i < tot_vert; i++ ) {
      lwvert[ 0 ] =  vert[ i ][ X ];
      lwvert[ 2 ] = -vert[ i ][ Y ];
      lwvert[ 1 ] =  vert[ i ][ Z ];
      ptid[ i ] = lwAddPoint( lwvert );
   }
   lwAddFace( tot_vert, ptid );
   free( ptid );
}


/*
======================================================================
lib_output_polypatch()

Output polygonal patch.  A patch is defined by a set of vertices and
their normals.  Since polygon patches aren't supported directly, we
call a routine to split the patch into triangles.
====================================================================== */

void lib_output_polypatch( int tot_vert, COORD3 *vert, COORD3 *norm )
{
   split_polygon( tot_vert, vert, norm );
}