/****************************************************************************\

  Copyright 1995 The University of North Carolina at Chapel Hill.
  All Rights Reserved.

  Permission to use, copy, modify and distribute this software and its
  documentation for educational, research and non-profit purposes,
  without fee, and without a written agreement is hereby granted,
  provided that the above copyright notice and the following three
  paragraphs appear in all copies.

  IN NO EVENT SHALL THE UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL BE
  LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR
  CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE
  USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY
  OF NORTH CAROLINA HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH
  DAMAGES.


  Permission to use, copy, modify and distribute this software and its
  documentation for educational, research and non-profit purposes,
  without fee, and without a written agreement is hereby granted,
  provided that the above copyright notice and the following three
  paragraphs appear in all copies.

  THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY
  WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE
  PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
  NORTH CAROLINA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

  The authors may be contacted via:

  US Mail:  Jonathan Cohen                      Amitabh Varshney
            Department of Computer Science      Department of Computer Science 
            Sitterson Hall, CB #3175            State University of New York
            University of N. Carolina           Stony Brook, NY 11794-4400, USA 
            Chapel Hill, NC 27599-3175
	    
  Phone:    (919)962-1749                       Phone: (516)632-8446 
	    
  EMail:    cohenj@cs.unc.edu                   varshney@cs.sunysb.edu

\****************************************************************************/
/*****************************************************************************\
  tubes.c
  --
  Description : Functions for building border tubes around mesh borders.

  ----------------------------------------------------------------------------
  $Source: /cvs/RenderPark/SE/Simplify/tubes.c,v $
  $Revision: 1.1.1.1 $
  $Date: 2000/04/06 15:35:32 $
  $Author: philippe $
  $Locker:  $
\*****************************************************************************/


/*----------------------------- Local Includes -----------------------------*/

#include <stdio.h>
#include <math.h>
#include <simplify.h>
#include <geometry.h>
#include <intersect.h>
#include <offset.h>
#include <tubes.h>
#include <ply.h>

/*----------------------------- Local Constants -----------------------------*/


/*------------------------------ Local Macros -------------------------------*/


/*------------------------------- Local Types -------------------------------*/


/*------------------------ Local Function Prototypes ------------------------*/

static void construct_border_polylines(Surface *model, PolyLine **polylines,
				       int *num_polylines);
static void construct_border_tubes(PolyLine *polylines, int num_polylines,
				   int num_sides, double width,
				   Surface *tubes);
static void construct_tube(PolyLine *polyline, int num_sides,
			    double width, Surface *tube);
static void append_tube(Surface *one_tube, Surface *tubes);

/*------------------------------ Local Globals ------------------------------*/


/*---------------------------------Functions-------------------------------- */

/*****************************************************************************\
 @ init_tubes()
 -----------------------------------------------------------------------------
 description : Create a set of polygonal tubes around mesh border edges and
               vertices.  These tubes are typically constructed very narrow
	       and then offset later to their full, non-self-intersecting
	       width.
 input       : Surface mesh, tube width, number of sides in tesselation around
               the tube.
 output      : Tube surface mesh (if there are multiple borders, this will
               contain one closed, connected mesh component for each cycle of
	       border edges).
 notes       :
\*****************************************************************************/
void init_tubes(Surface *model, Surface *tubes, double width, int num_sides)
{
    PolyLine *polylines;
    int       num_polylines;
    
    fprintf(stderr, "building polylines...");
    fflush(stderr);
    
    construct_border_polylines(model, &polylines, &num_polylines);

    fprintf(stderr, "done\n");
    fprintf(stderr, "building tubes...");
    fflush(stderr);

    construct_border_tubes(polylines, num_polylines, num_sides, width, tubes);

    fprintf(stderr, "done\n");
    fflush(stderr);
    
    return;
} /** End of init_tubes() **/

/*****************************************************************************\
 @ construct_border_polylines()
 -----------------------------------------------------------------------------
 description : Find all the borders on a surface mesh and represent each as a
               cyclic polyine.
 input       : Surface mesh
 output      : A set of polylines and how many polylines there are.
 notes       :
\*****************************************************************************/
static void construct_border_polylines(Surface *model, PolyLine **polylines,
				       int *num_polylines)
{
    int            i, j;
    unsigned char *bvert_used;
    Vertex        *vert;
    PolyLine      *polyline;
    Vertex        *start_vert, *adjacent_vert;
    Edge          *edge;
    Vertex        *prev_vert;
    
    ALLOCN(bvert_used, unsigned char, model->num_verts);
    for (i=0; i<model->num_verts; i++)
    {
	vert = &(model->verts[i]);
	if (vert->num_tris == vert->num_edges) /* non-border */
	    bvert_used[i] = TRUE;
	else if (vert->num_tris == (vert->num_edges-1)) /* border */
	    bvert_used[i] = FALSE;
	else /* vertex on multiple borders -- non-manifold */
	{
	    fprintf(stderr,
		    "Vertex #%d is non-manifold (on multiple borders)\n", i);
	    exit(1);
	}
    }

    for (i=0, *num_polylines = 0; i<model->num_verts; i++)
    {
	if (bvert_used[i] == TRUE)
	    continue;

	REALLOCN(*polylines, PolyLine , *num_polylines, *num_polylines + 1);
	polyline = &((*polylines)[*num_polylines]);

	start_vert = &(model->verts[i]);
	ALLOCN(polyline->verts, Vertex *, 1);
	polyline->verts[0] = start_vert;
	polyline->num_verts = 1;
	bvert_used[start_vert->id] = TRUE;
	for (j=0; j<start_vert->num_edges; j++)
	{
	    edge = start_vert->edges[j];
	    adjacent_vert =
		(edge->verts[0] == start_vert) ?
		    edge->verts[1] : edge->verts[0];
	    if (bvert_used[adjacent_vert->id] == FALSE) /* border vert */
		break;
	}
	if (j >= start_vert->num_edges)
	{
	    fprintf(stderr, "Couldn't find second border vert\n");
	    exit(1);
	}
	vert = adjacent_vert;
	prev_vert = start_vert;
	while (vert != start_vert)
	{
	    REALLOCN(polyline->verts, Vertex *,
		     polyline->num_verts, polyline->num_verts+1);
	    polyline->verts[polyline->num_verts++] = vert;
	    bvert_used[vert->id] = TRUE;
	    
	    for (j=0; j<vert->num_edges; j++)
	    {
		edge = vert->edges[j];
		if (edge->tris[1] != NULL)
		    continue;
		adjacent_vert =
		    (edge->verts[0] == vert) ?
			edge->verts[1] : edge->verts[0];
		if (adjacent_vert == prev_vert)
		    continue;
		break;
	    }
	    if (j >= vert->num_edges)
	    {
		fprintf(stderr, "Couldn't follow border\n");
		exit(1);
	    }
	    
	    prev_vert = vert;
	    vert = adjacent_vert;
	}
	(*num_polylines)++;
    }
    return;
} /** End of construct_border_polylines() **/

/*****************************************************************************\
 @ calculate_tube_verts()
 -----------------------------------------------------------------------------
 description : Calculate the positions of the tube vertices around the
               polyline.  For each polyline vertex, an appropriate plane is
	       chosen, and num_sides vertices are placed on this plane in a
	       circular fashion.
 input       : A cyclic polyline, number of sides around tube, width of tube 
 output      : A list of vertices and the list length
 notes       : The tricky bit in all this is placing the plane and choosing an
               appropriate coordinate system in the plane to encourage smooth,
	       non-self-intersecting tubes.  This is not guaranteed, but it
	       seems to work even in some pretty tough cases.  If the border
	       tube DOES self intersect, it will probably just prevent any
	       simplification to occur along the mesh borders in those
	       regions.
\*****************************************************************************/
void calculate_tube_verts(PolyLine *polyline, Vertex **verts,
			  int *num_verts, int num_sides, double width)
{
    Point *plane_points;
    int    point_count;
    int vert_count;
    Vertex *vert, *next_vert, *prev_vert;
    Vector forward_vec, backward_vec, binormal;
    Vector tangent;
    Vector plane_x, plane_y;
    Vector edge_vec, bisector, bisector_perp;
    Matrix mat;
    Vector normal_component;
    double normal_scale;
    Vector edge_bisector, prev_bisector, orig_bisector, forward_bisector;

    *num_verts = 0;
    ALLOCN(*verts, Vertex, num_sides*polyline->num_verts);

    /* build a circle in the plane with num_sides vertices */
    ALLOCN(plane_points, Point, num_sides);
    for (point_count = 0; point_count < num_sides; point_count++)
    {
	plane_points[point_count][X] =
	    width * cos(2.0*M_PI* (double)point_count/(double)num_sides);
	plane_points[point_count][Y] =
	    width*sin(2.0*M_PI * (double)point_count/(double)num_sides);
	plane_points[point_count][Z] = 0.0;
    }

    for (vert_count = 0; vert_count < polyline->num_verts; vert_count++)
    {
	vert = polyline->verts[vert_count];
	next_vert = polyline->verts[(vert_count+1)%polyline->num_verts];
	prev_vert = polyline->verts[(vert_count+polyline->num_verts-1) %
				      polyline->num_verts];

	VEC3_V_OP_V(forward_vec, next_vert->coord, -, vert->coord);
	VEC3_V_OP_V(backward_vec, prev_vert->coord, -, vert->coord);
	NORMALIZE3(forward_vec);
	NORMALIZE3(backward_vec);

	/* calculate a tangent vector at the polyline vertex, which will be
	   used as the normal to the plane on which the tube vertices lie */
	if (DOTPROD3(forward_vec, backward_vec) < 0) /* angle > 90
							      degrees */
	{
	    VEC3_V_OP_V(tangent, next_vert->coord, -, prev_vert->coord);
	    NORMALIZE3(tangent);
	}
	else
	{
	    VEC3_VOPV_OP_S(bisector, forward_vec, +, backward_vec, *, 0.5);
	    CROSSPROD3(binormal, forward_vec, backward_vec);
	    CROSSPROD3(tangent, bisector, binormal);
	    NORMALIZE3(tangent);
	}

	/* Now compute plane_x and plane_y, which determine how we project our
	   circle of vertices onto the plane.  These should be chosen to
	   minimize any twisting of the tube. */
	if (vert_count == 0)
	{
	    /* an arbitrary plane_x and plane_y in the normal plane
	       (perpendicular to the tangent) should do */
	    VEC3_V_OP_V(edge_vec, next_vert->coord, -, vert->coord);
	    CROSSPROD3(plane_x, tangent, edge_vec);
	    CROSSPROD3(plane_y, plane_x, tangent);
	    NORMALIZE3(plane_x);
	    NORMALIZE3(plane_y);
	    VEC3_VOPV_OP_S(bisector, plane_x, +, plane_y, *, 0.5);
	    NORMALIZE3(bisector);
	    VEC3_ASN_OP(orig_bisector, =, bisector);
	}
	else if (vert_count != (polyline->num_verts - 1))
	{
	    /* Use a two step projection process to project the coordinate
	       system from the previous polyline vertex to the current one.
	       First project onto the plane normal to the intervening edge,
	       then project from that plane to the normal plane of the current
	       polyline vertex.  This two step process helps reduce twisting
	       and is more stable than a direct one-step projection in the
	       presence of sharp polyline angles. */
	    
	    VEC3_V_OP_V(edge_vec, vert->coord, -, prev_vert->coord);
	    NORMALIZE3(edge_vec);
	    
	    /* project x- and y-axis bisector from previous normal plane
	       onto the plane normal of the previous edge */
	    normal_scale = DOTPROD3(prev_bisector, edge_vec);
	    VEC3_V_OP_S(normal_component, edge_vec, *, normal_scale);
	    VEC3_V_OP_V(edge_bisector, prev_bisector, -, normal_component);

	    /* project from the previous edge's normal plane onto the current
	       normal plane */
	    normal_scale = DOTPROD3(edge_bisector, tangent);
	    VEC3_V_OP_S(normal_component, tangent, *, normal_scale);
	    VEC3_V_OP_V(bisector, edge_bisector, -, normal_component);
	    NORMALIZE3(bisector);

	    /* construct perpendicular x- and y-axes from the
	       projected axis bisector */

	    CROSSPROD3(bisector_perp, tangent, bisector);
	    NORMALIZE3(bisector_perp);

	    VEC3_VOPV_OP_S(plane_x, bisector, +, bisector_perp, *, 0.5);
	    NORMALIZE3(plane_x);
	    
	    CROSSPROD3(plane_y, plane_x, tangent);
	    NORMALIZE3(plane_y);
	}
	else
	{
	    /* For the final polyline vertex, average the projected coordinate
	       systems of the previous vertex and the initial vertex.  Even if
	       the tube twists gradually, the starting and ending sections may
	       have as much as a 180 degree twist between them.  This
	       averaging reduces the maximum possible twist to 90 degrees. */


	       
	    /* project x- and y-axis bisector from previous normal plane
	       onto the plane normal of the previous edge */

	    VEC3_V_OP_V(edge_vec, vert->coord, -, prev_vert->coord);
	    NORMALIZE3(edge_vec);
	    normal_scale = DOTPROD3(prev_bisector, edge_vec);
	    VEC3_V_OP_S(normal_component, edge_vec, *, normal_scale);
	    VEC3_V_OP_V(edge_bisector, prev_bisector, -, normal_component);

	    /* project from the previous edge's normal plane onto the current
	       normal plane */
	    normal_scale = DOTPROD3(edge_bisector, tangent);
	    VEC3_V_OP_S(normal_component, tangent, *, normal_scale);
	    VEC3_V_OP_V(bisector, edge_bisector, -, normal_component);	
	    NORMALIZE3(bisector);

	    /* project x- and y-axis bisector from initial normal plane
	       onto the plane normal of the final edge */

	    VEC3_V_OP_V(edge_vec, next_vert->coord, -, vert->coord);
	    NORMALIZE3(edge_vec);
	    normal_scale = DOTPROD3(orig_bisector, edge_vec);
	    VEC3_V_OP_S(normal_component, edge_vec, *, normal_scale);
	    VEC3_V_OP_V(edge_bisector, orig_bisector, -, normal_component);

	    /* project from the final edge's normal plane onto the current
	       normal plane */
	    normal_scale = DOTPROD3(edge_bisector, tangent);
	    VEC3_V_OP_S(normal_component, tangent, *, normal_scale);
	    VEC3_V_OP_V(forward_bisector, edge_bisector, -, normal_component);	
	    NORMALIZE3(forward_bisector);

	    /* bisect the bisector from the previous plane and the initial
	       plane to reduce the error of final twisting effects */
	    VEC3_VOPV_OP_S(bisector, bisector, +, forward_bisector, *, 0.5);
	    
	    
	    /* construct perpendicular x- and y-axes from the
	       projected axes */

	    CROSSPROD3(bisector_perp, tangent, bisector);
	    NORMALIZE3(bisector_perp);

	    VEC3_VOPV_OP_S(plane_x, bisector, +, bisector_perp, *, 0.5);
	    NORMALIZE3(plane_x);
	    
	    CROSSPROD3(plane_y, plane_x, tangent);
	    NORMALIZE3(plane_y);
	}
	
	
	VEC3_ASN_OP(prev_bisector, =, bisector);


	/* Construct the matrix for transforming a circle of points around the
	   origin on the xy-plane onto the plane through the polyline vertex
	   perpendicular to its tangent.  Rotate to follow the plane_x and
	   plane_y coordinate system on this normal plane. */
	
	/* x-axis rotates into plane_x */
	mat[0][0] = plane_x[X];
	mat[1][0] = plane_x[Y];
	mat[2][0] = plane_x[Z];

	/* y-axis rotates into plane_y */
	mat[0][1] = plane_y[X];
	mat[1][1] = plane_y[Y];
	mat[2][1] = plane_y[Z];

	/* z-axis rotates into tangent */
	mat[0][2] = tangent[X];
	mat[1][2] = tangent[Y];
	mat[2][2] = tangent[Z];

	/* translate to current vertex */
	mat[0][3] = polyline->verts[vert_count]->coord[X];
	mat[1][3] = polyline->verts[vert_count]->coord[Y];
	mat[2][3] = polyline->verts[vert_count]->coord[Z];
	
	for (point_count = 0; point_count < num_sides; point_count++)
	{
	    TRANSFORM_POINT((*verts)[*num_verts].coord,
			    mat, plane_points[point_count]);
	    TRANSFORM_VECTOR((*verts)[*num_verts].normal,
			     mat, plane_points[point_count]);
	    NORMALIZE3((*verts)[*num_verts].normal);
#ifdef VERTEX_EPSILONS
                (*verts)[*num_verts].epsilon = vert->epsilon;
#endif
	    (*num_verts)++;
	}
    }
    return;
} /** End of calculate_tube_verts() **/




/*****************************************************************************\
 @ construct_border_tubes()
 -----------------------------------------------------------------------------
 description : Build a tube around each polyline.
 input       : Set of cyclic polylines, number of sides to use in tesselating
               around the tube, and tube width.
 output      : Surface mesh of the tubes.
 notes       :
\*****************************************************************************/
static void construct_border_tubes(PolyLine *polylines, int num_polylines,
				   int num_sides, double width,
				   Surface *tubes)
{
    int   i;
    Surface tube;
    PolyLine *polyline;
    

    /* allocate space for the tubes */
    tubes->num_verts = tubes->num_edges = tubes->num_tris = 0;
    for (i=0; i<num_polylines; i++)
    {
	polyline = &(polylines[i]);
	tubes->num_verts += num_sides*polyline->num_verts;
	tubes->num_edges += 3*num_sides*polyline->num_verts;
	tubes->num_tris += 2*num_sides*polyline->num_verts;
    }
    ALLOCN(tubes->verts, Vertex, tubes->num_verts);
    ALLOCN(tubes->edges, Edge, tubes->num_edges);
    ALLOCN(tubes->tris, Triangle, tubes->num_tris);
    tubes->num_verts = tubes->num_edges = tubes->num_tris = 0;

    /* construct the tubes */
    for (i=0; i<num_polylines; i++)
    {
	construct_tube(&(polylines[i]), num_sides, width, &tube);
	append_tube(&tube, tubes);
    }

    return;
}


/*****************************************************************************\
 @ construct_tube()
 -----------------------------------------------------------------------------
 description : Build a border tube.  First calculate the coordinates of the
               vertices, then build the tube's vertex, edge, and face
	       structures, complete with nasty, nasty pointers.
 input       : Polyline to build tube around, number of sides to tesselate
               around tube, tube width.
 output      : Surface mesh defining the tube.
 notes       : You _really_ don't want to mess with this code...
               (There MUST be some easier way to do this...)

 I'd better describe the topology of the tube here in one place so I can get
 all the vertex, edge, and face pointers correct.  First, the polyline is
 closed, with the last vertex implicitly connected to the first.  I'll
 construct a torus around this polyline.

 If I'm at a vertex on the polyline looking down an edge towards the "next"
 vertex, I should be surrounded by a counter-clockwise circle of num_sides
 vertices.  The first num_sides vertices are around the first polyline vertex,
 the next num_sides vertices around the second polyline vertex, etc.

 Now for the edges.  There are "circle" edges that connect the vertices around
 each polyline vertex.  There are "lengthwise" edges that are parallel to the
 polylines, and there are "diagonal" edges that triangulate the tube polygons.
 The edges will be numbered as follows.  The first num_sides edges will be
 circle edges surrounding the first polyline vertex.  The next num_sides
 edges will be lengthwise edges connecting the vertices of the first circle to
 the vertices of the second circle (the circle around the second polyline
 vertex).  The next num_sides vertices will be the diagonal edges between the
 first and the second circles.  The first diagonal will connect the second
 vertex of the first circle to the first vertex of the second circle.

 The "top" faces will each share an edge with some circle and share a vertex
 with the _next_ circle.  The "bottom" faces will share an edge with some
 circle and share a vertex with the _previous_ circle.  The faces shall be
 numbered as follows.  The first num_sides faces will be bottom faces between
 the first and second circles.  The next num_sides faces will be top faces,
 also between the first and second circles.  The first top face will sit on
 the first edge of the first circle, and share the first vertex of the second
 circle.  The first bottom face will sit beneath the first edge of the second
 circle and share the second vertex of the first circle.
 
\*****************************************************************************/
static void construct_tube(PolyLine *polyline, int num_sides,
			   double width, Surface *tube)
{
    int       i, j;
    int       edge_index, face_index, vert_index;
    Edge     *edge;
    Triangle *face;
    Vertex   *vert;
    
    /* get the coordinates and normals of the tube vertices */
    calculate_tube_verts(polyline, &(tube->verts), &(tube->num_verts),
			 num_sides, width);
    
    ALLOCN(tube->tris, Triangle, polyline->num_verts * 2 * num_sides);
    tube->num_tris = polyline->num_verts * 2 * num_sides;

    ALLOCN(tube->edges, Edge, polyline->num_verts * 3 * num_sides);
    tube->num_edges = polyline->num_verts * 3 * num_sides; /* f-e+v==0 for torus */

    

    /* You are now entering POINTER HELL!!!!  Turn back while you still can.
       Modify this code at your own risk.  You have been warned!!! */

    
    
    /* add the "circle" edges surrounding each polyline vertex */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    edge_index = i*num_sides*3 + j;
	    edge = &(tube->edges[edge_index]);
	    edge->id = edge_index;
	    edge->verts[0] = &(tube->verts[i*num_sides+j]);
	    edge->verts[1] = &(tube->verts[i*num_sides+(j+1)%num_sides]);
	    edge->tris[0] = &(tube->tris[2*i*num_sides+j]);
	    edge->tris[1] =
		&(tube->tris[((i+polyline->num_verts-1)%polyline->num_verts) *
			     2*num_sides + num_sides + j]);
	}
    }

    /* add the "lengthwise" edges parallel to the polyline */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    edge_index = i*num_sides*3 + num_sides + j;
	    edge = &(tube->edges[edge_index]);
	    edge->id = edge_index;
	    edge->verts[0] = &(tube->verts[i*num_sides+j]);
	    edge->verts[1] = &(tube->verts[((i+1)%polyline->num_verts) *
					   num_sides + j]);
	    edge->tris[0] = &(tube->tris[2*i*num_sides+j]);
	    edge->tris[1] = &(tube->tris[2*i*num_sides + num_sides +
					 (j+num_sides-1)%num_sides]);
	}
    }
    
    /* add the "diagonal" edges that triangulate the sides */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    edge_index = i*num_sides*3 + 2*num_sides + j;
	    edge = &(tube->edges[edge_index]);
	    edge->id = edge_index;
	    edge->verts[0] = &(tube->verts[i*num_sides + (j+1)%num_sides]);
	    edge->verts[1] = &(tube->verts[((i+1)%polyline->num_verts) *
					   num_sides + j]);
	    edge->tris[0] = &(tube->tris[2*i*num_sides+j]);
	    edge->tris[1] = &(tube->tris[2*i*num_sides+num_sides+j]);
	}
    }
    
    /* add the "bottom" faces */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    face_index = i*num_sides*2 + j;
	    face = &(tube->tris[face_index]);
	    face->id = face_index;
	    face->verts[0] = &(tube->verts[i*num_sides+j]);
	    face->verts[1] = &(tube->verts[((i+1)%polyline->num_verts) *
					   num_sides + j]);
	    face->verts[2] = &(tube->verts[i*num_sides+(j+1)%num_sides]);
	    face->edges[0] = &(tube->edges[i*num_sides*3 + num_sides + j]);
	    face->edges[1] = &(tube->edges[i*num_sides*3 + 2*num_sides + j]);
	    face->edges[2] = &(tube->edges[i*num_sides*3 + j]);
	    find_plane_eq(face->verts[0]->coord,
			  face->verts[1]->coord,
			  face->verts[2]->coord,
			  face->plane_eq);
	}
    }
    
    
    /* add the "top" faces */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    face_index = i*num_sides*2 + num_sides + j;
	    face = &(tube->tris[face_index]);
	    face->id = face_index;
	    face->verts[0] = &(tube->verts[((i+1)%polyline->num_verts) *
					   num_sides + j]);
	    face->verts[1] = &(tube->verts[((i+1)%polyline->num_verts) *
					   num_sides + (j+1)%num_sides]);
	    face->verts[2] = &(tube->verts[i*num_sides+(j+1)%num_sides]);
	    face->edges[0] = &(tube->edges[((i+1)%polyline->num_verts) *
					   num_sides*3 + j]);
	    face->edges[1] = &(tube->edges[i*num_sides*3 + num_sides +
					   (j+1)%num_sides]);
	    face->edges[2] = &(tube->edges[i*num_sides*3 + 2*num_sides + j]);
	    find_plane_eq(face->verts[0]->coord,
			  face->verts[1]->coord,
			  face->verts[2]->coord,
			  face->plane_eq);
	}
    }
    
    /* add the edge and face indices to vertices */
    for (i=0; i<polyline->num_verts; i++)
    {
	for (j=0; j<num_sides; j++)
	{
	    vert_index = i*num_sides + j;
	    vert = &(tube->verts[vert_index]);
	    vert->id = vert_index;
	    ALLOCN(vert->edges, Edge *, 6);
	    ALLOCN(vert->tris, Triangle *, 6);
	    vert->num_tris = vert->num_edges = 6;
	    vert->edges[0] = &(tube->edges[i*num_sides*3 +
					   (j+num_sides-1) % num_sides]);
	    vert->edges[1] = &(tube->edges[i*num_sides*3 + 2*num_sides +
					   (j+num_sides-1) % num_sides]);
	    vert->edges[2] = &(tube->edges[i*num_sides*3 + num_sides + j]);
	    vert->edges[3] = &(tube->edges[i*num_sides*3 + j]);
	    vert->edges[4] = &(tube->edges[((i+polyline->num_verts-1) %
					   polyline->num_verts)*num_sides*3 +
					   2*num_sides + j]);
	    vert->edges[5] = &(tube->edges[((i+polyline->num_verts-1) %
					   polyline->num_verts)*num_sides*3 +
					   num_sides + j]);
	    vert->tris[0] = &(tube->tris[i*num_sides*2 +
					 (j+num_sides-1)%num_sides]);
	    vert->tris[1] = &(tube->tris[i*num_sides*2 + num_sides +
					 (j+num_sides-1)%num_sides]);
	    vert->tris[2] = &(tube->tris[i*num_sides*2 + j]);
	    vert->tris[3] = &(tube->tris[((i+polyline->num_verts-1) %
					  polyline->num_verts)*num_sides*2 +
					 num_sides + j]);
	    vert->tris[4] = &(tube->tris[((i+polyline->num_verts-1) %
					  polyline->num_verts)*num_sides*2 +
					 j]);
	    vert->tris[5] = &(tube->tris[((i+polyline->num_verts-1) %
					  polyline->num_verts)*num_sides*2 +
					 num_sides +
					 (j+num_sides-1)%num_sides]);
	}
    }
    return;
} /** End of construct_tube() **/

/*****************************************************************************\
 @ append_tube()
 -----------------------------------------------------------------------------
 description : Append a single tube mesh to a mesh containing several tubes
               (mostly so I could focus on the layout of a single tube at a
	       time in other routines).
 input       : A single tube mesh and a mesh containing zero or more tubes.
 output      : The multi-tube mesh now has the single tube appended, and the
               single tube is deallocated.
 notes       : assumes space in tubes is already allocated properly
\*****************************************************************************/
static void append_tube(Surface *one_tube, Surface *tubes)
{
    int       i, j;
    int       start_vert, start_edge, start_tri;
    Vertex   *src_vert, *dest_vert;
    Edge     *src_edge, *dest_edge;
    Triangle *src_tri, *dest_tri;
    Vertex   *vert;
    

    start_vert = tubes->num_verts;
    start_edge = tubes->num_edges;
    start_tri  = tubes->num_tris;
    
    /* copy vertices */
    for (i=0; i<one_tube->num_verts; i++)
    {
	src_vert = &(one_tube->verts[i]);
	dest_vert = &(tubes->verts[tubes->num_verts]);
	dest_vert->id = tubes->num_verts;
	VEC3_ASN_OP(dest_vert->coord, =, src_vert->coord);
	VEC3_ASN_OP(dest_vert->normal, =, src_vert->normal);
	
	/* copy edge pointers */
	ALLOCN(dest_vert->edges, Edge *, src_vert->num_edges);
	for (j=0; j<src_vert->num_edges; j++)
	    dest_vert->edges[j] =
		&(tubes->edges[src_vert->edges[j]->id + start_edge]);
	dest_vert->num_edges = src_vert->num_edges;

	/* copy tri pointers */
	ALLOCN(dest_vert->tris, Triangle *, src_vert->num_tris);
	for (j=0; j<src_vert->num_tris; j++)
	    dest_vert->tris[j] =
		&(tubes->tris[src_vert->tris[j]->id + start_tri]);
	dest_vert->num_tris = src_vert->num_tris;
	
	tubes->num_verts++;
    }
    
    /* copy edges */
    for (i=0; i<one_tube->num_edges; i++)
    {
	src_edge = &(one_tube->edges[i]);
	dest_edge = &(tubes->edges[tubes->num_edges]);
	dest_edge->id = tubes->num_edges;

	/* copy vertex pointers */
	for (j=0; j<2; j++)
	    dest_edge->verts[j] =
		&(tubes->verts[src_edge->verts[j]->id + start_vert]);
	
	/* copy tri pointers */
	for (j=0; j<2; j++)
	    dest_edge->tris[j] =
		&(tubes->tris[src_edge->tris[j]->id + start_tri]);

	tubes->num_edges++;
    }
    
    /* copy tris */
    for (i=0; i<one_tube->num_tris; i++)
    {
	src_tri = &(one_tube->tris[i]);
	dest_tri = &(tubes->tris[tubes->num_tris]);
	dest_tri->id = tubes->num_tris;

	/* copy vertex pointers */
	for (j=0; j<3; j++)
	    dest_tri->verts[j] =
		&(tubes->verts[src_tri->verts[j]->id + start_vert]);

	/* copy edge pointers */
	for (j=0; j<3; j++)
	    dest_tri->edges[j] =
		&(tubes->edges[src_tri->edges[j]->id + start_edge]);

	VEC_ASN_OP(dest_tri->plane_eq, =, src_tri->plane_eq, 4);
	
	tubes->num_tris++;
    }
    
    /* free one_tube */
    for (i=0; i<one_tube->num_verts; i++)
    {
	vert = &(one_tube->verts[i]);
	FREE(vert->edges);
	FREE(vert->tris);
    }
    FREE(one_tube->verts);
    FREE(one_tube->edges);
    FREE(one_tube->tris);
    
    return;
} /** End of append_tube() **/

/*****************************************************************************\
  $Log: tubes.c,v $
  Revision 1.1.1.1  2000/04/06 15:35:32  philippe
  Initial CVS release

  Revision 1.5  1997/04/10 20:10:46  cohenj
  Added Amitabh's name to credits

  Revision 1.4  1996/11/15 20:35:20  cohenj
  modified to use vector/vector/scalar operators with CORRECT
  parenthisization

  Revision 1.3  1996/04/08 19:15:06  cohenj
  added support for per-vertex epsilons, procedure comments,
  and copyright notice

 * Revision 1.2  1995/10/16  16:29:41  cohenj
 * Cleaned up.
 *
 * Revision 1.1  95/09/15  16:28:21  cohenj
 * Initial revision
 * 
 * Revision 1.2  1995/09/12  18:24:44  cohenj
 * fixed calculate_tube_verts() to give the least twisted tubes, which seem to produce tubes 
without
 * self-intersections (at least for the bunny borders)
 *
 * Revision 1.1  1995/08/30  20:57:04  cohenj
 * Initial revision
 *
\*****************************************************************************/