/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>


/*********************************************************************/
/*

  Author: Andre GUEZIEC, March/April 1997

  */
/*********************************************************************/


#ifndef _simplifysurface_
#define _simplifysurface_

/* SimplifySurface

   is a module that approximates a triangulated surface and resamples data

 */

/* partial bug log 
   ~~~~~~~~~~~~~~~

   this is a partial log of the bugs I discovered when developing this module.
   it is there mainly to help me, and other potential developers working 
   on this module spend less time hunting bugs, especially if we want to
   change the code and redo the same errors.


   3/29/97   DXCull, DXInvalidateConnections, DXInvalidateUnreferencedPositions

   I was trying to run something like:

   obj1 = DXInvalidateUnreferencedPositions(obj)
   obj2 = DXCull(obj1)

   That corresponded to what the User manual suggested
   but I was getting an error at DXGetObjectClass with the result of DXCull

   I decided to code things differently:

   new_in[0] = DXInvalidateConnections(in[0]);
   if (!DXInvalidateUnreferencedPositions(new_in[0])) goto error;
   if (!DXCull(new_in[0]))  goto error;
 
   
   ...........................................................................
   

   3/31/97   DXNewArrayV must be used in place of DXNewArray if a "shape" int * vector
   is passed to the routine rather than an int. Typically if rank > 1 but also if
   rank = 1 and shape is obtained from DXArrayInfo, which expects an int *.
   
   The problem with this bug is that the error message appeared in the end of the module
   something like Out of Memory Out of Memory, and I had no obvious way of tracking it down to the
   source.


   ...........................................................................


   3/31/97      _dxfVertexNormalsfromTriangleNormals(nor, new_nV, new_nT, new_connections,
		  				    face_normals, new_face_areas);
						   
   was executed *after* DXFree((Pointer)new_connections); which is a stupid thing to do!
   I was freeing after 
		      DXAddArrayData(con,  0, new_nT, new_connections);
   and when later on I added the vertex normals from triangle normals routine,
   I forgot that this array was freed. Part of the problem is that I wasnt used to the
   style of allocating and freeing with DXAllocate and DXFree whereby you do a goto error 
   if something bad happens. I was doing things differently:
 
   return code = 0
   buffer = (malloc)

   if (buffer)  { do the work;  free(buffer); set return code to 1}
   else         { write error message }
   
   return(return code)


   ...........................................................................

             

   3/31/97 DXConvertArray() sets an error code if the array is already of the same type that requires
   conversion. Same thing: it is not documented and it took me a while to realize what was going wrong
   since the error was showing up in the end of the module execution, and was not informative:

   Out of Memory, or Null Object, or something similar


   ...........................................................................



   4/2/97 MACRO DEFINITIONS: I had a hard time compiling (2 hours) because of syntax errors
   in the macro definitions that are not picked up by the compiler. They are detected much
   later and point so seemingly incomprehensible errors in the c-code.


   ...........................................................................


   4/2/97 SurfEdgeLabel: when building the "e_table" hash table of edges, I was forgetting 
   to initialize the label of the edge to 1 

   SurfEdgeLabel(edge) = 1;

   Accordingly, edges that were boundary edges had a label of 28800 or something like that
   (an arbitrary value) and were misclassified as interior edges.

   then, the routine that was looping around edges was crashing because it expected boundary
   edges to be interior edges.

   This bug was easy to locate in dbx: I started from the location where the program crashed,
   I printed the edge array and noticed right away that
   the labels were funny. I was leaded quickly to the routine that generated the edges



   ...........................................................................


   4/2/97 DXAllocate the valence instead of DXAllocateZero: stupid error. I had only two
   DXAllocate instead of DXAllocateZero in the entire program. 

   I noticed it in dbx because the module was simplifying surface so little.
   by stopping at the right line, I notices way too many topologically infeasible collapses.
   That made me want to check the valence

   simp_data->boundary_vert was also DXAllocated and that was probably another bug, fixed
   by hitting two birds with a stone.




   ...........................................................................



   4/3/97 rank and shape of a flag input. By default, the module builder was using the test

   if (type != TYPE_INT || category != CATEGORY_REAL || rank != 1 || shape != 2) 
      DXSetError(ERROR_INVALID_DATA, "input \"preserve_volume\"");
   
   the right test to use is 

   if (type != TYPE_INT || category != CATEGORY_REAL || !((rank == 0) || ((rank == 1)&&(shape == 1))))



   ...........................................................................



   4/3/97 DXGetComponentValue(field, "positional_error")

   the "positional_error" component was ignored by the simplification because
   of a faulty test for the "dep" attribute
   if was using 

   if (strcmp("positions", DXGetString((String)attr)))

   instead of 

   if (!strcmp("positions", DXGetString((String)attr)))

   which is the "c" translation of : are "positions" and DXGetString((String)attr) the same

   the (!) will presumably get me many more times in the future. This is the advantage of
   documenting the bugs. Hopefully I'll remember this one.



   ...........................................................................


   4/3/97 bug in the routine _dxfFloatDataBoundingBoxDiagonal. I noticed it when running
   an example that fed several isosurfaces to the module. The bounding box value was inconsistent.

   the code was as follows:

   for(j=0;j<data_dim;j++)
    {
      bbox[j] = bbox[j+data_dim] = data[j];
    }
    

  cur_data = data + data_dim;

  for (i=0;i<n_data;i++,cur_data+=data_dim)

    for(j=0;j<data_dim;j++)
      {
	if (cur_data[j] < bbox[j]) bbox[j] = cur_data[j];

	else if (cur_data[j] > bbox[j+data_dim]) bbox[j+data_dim] = cur_data[j];
      }

      The problem was that the first set of data values was used to initialize the bounding box
      and *then* after incrementing cur_data, I was examining n_data values.
      Hence I was peeking outside of the data array and picking any ridiculous value there or
      risking a core dump

      one solution is to do the loop as

      for (i=1;i<n_data;i++,cur_data+=data_dim)


   ...........................................................................


   4/3/97 GET_POSITIONAL_ERROR the existing positional error was not used. I noticed it because
   the module was apparently simplifying too much on a second simplification.

   what happened is that the portion of the code 

   starting with

   array = (Array)DXGetComponentValue(field, "positional_error");

   was just before a goto label GET_OTHER_INPUTS, and the code jumped to that goto in case
   no data "dep" positions was available, which concerns most of the cases. 

   so I changed the name of the goto label to  GET_POSITIONAL_ERROR and placed it before the
   relevant code.

   ...........................................................................


   4/3/97 this was an important omission rather than a bug.

   when running the debugged module on non-manifold input, I realize that if
   some faces in the input were degenerate, (if the three vertices of the face
   weren't all different), the conversion to manifold would not eliminate such faces.

   Curiously, this elimination of degenerate faces was done in the first version
   of the conversion to manifold, but when re-programming it I decided that it
   was not necessary. I think that it *is* necessary, because it perturbs 
   the joining of the corners. I am supposed to join corner c1 and c2 and c3 and c4
   and all the corners must be different. If a face is degenerate, then for
   instance c1 and c3 are the same and I don't know what will happen.

   This is a nasty omission, because a lot of code was depending on the number of
   faces not being changed in the conversion to manifold. Now there is a "face_lut"
   and a lot of code needed to be updated. This took approximately 3 hours.
   
   However, after incorporating this capability, problems still subsist,
   leading us to the next bug:
   

   ...........................................................................


   4/3/97 - 4/4/97 dealing with erroneous input
  
   in order to convert non manifold input surfaces to manifold surfaces,
   it is necessary to identify the manifold edges and to join triangles
   that share a manifold edge. If the edge is not manifold, no triangles can
   be joined there.

   ...........................................................................

   4/4/97 with new method for fixing non-manifold input, the edge count
   was not correct. I was counting too many edges and hence there were some
   additional spurious edges.

   In fact, in the routine _dxfLabelEdges() it is absolutely *necessary*
   to compute the number of edges again, as some triangles might be
   implicitely joined and there may be fewer edges than we previously thought.

   This bug was very difficult to debug because the -g version was running fine
   in the debugger, but I was convinced that it was a problem of memory being overwritten.
   In practice, no memory was overwritten but there were fewer valid elements in the
   edges array than the counter predicted leading in any catastrophic behavior when
   such edges were processed.

   ...........................................................................

   4/4/97 SIGH! as if it wasn't enough for today: I fixed the counting problem for
   the number of edges nE, and thus passed the pointer &nE to a bunch of routines

   I forgot that _dxfInitHashTable(e_table, _dxfFlog2(nE) + 1,...
   requires an estimate of the number of edges ! nE.
   It was receiving a NE of 0!!!, resulting in all the hash entries being added to the
   same list. 
   
   I noticed it because the program was soooo slooooow! 

   ...........................................................................

   4/6/97 definitions of uint, u_int, ushort, u_short ...

	the problem is that if a compiler knows about them on a particular
	platform, another compiler on a different platform might not know!

	below is a series of directives ifdef and include that seem to
	work

   ...........................................................................

   4/7/97 DXArrayConvert again
   
   I was getting an error ERROR: SimplifySurface: Invalid data: data ranks dont match
   that after one hour of debugging, I attributed to DXArrayConvert again
   so now I call instead of DXArrayConvert(array, TYPE_FLOAT, CATEGORY_REAL, 0)
   
   data_array = DXArrayConvert(array, TYPE_FLOAT, CATEGORY_REAL, rank);

   ...........................................................................

    4/7/97 

	line 6226 of the _simplifysurface.c file
    
    the line was

    s_comp[i] = s_comp[star1[i1]];

    instead of:

    s_comp[i] = t_comp[star1[i1]];

   resulting in reading the memory in some strange place
   it took me 3 hours to hunt this bug which was causing a segmentation
   fault on a dec machine (but not on other machines apparently)
   the debugger dbx was very useful for that, as I could notice that
   the array s_comp was not filled properly. It was surprisingly hard
   to finally find this typo, though, who easily survived compiling.

   ...........................................................................

    4/7/97 

     *new_edges =  (EdgeS *) DXReAllocate(*new_edges, nE * sizeof (EdgeS));
     line 760 of _simplifysurface.c

    the problem is that I was reallocating *after* filling up the h-table
    of edges. The edges were copied to another place and the edge entries
    in the h-table were pointing to nothing=or to a completely different
    piece of memory.

    This bug only appeared in the purified version of the dxexec because
    the reallocate was a realloc with a smaller size, which results in
    leaving the memory where it is with most implementations, except
    the implementation of Paula for the memorystubs.o which mallocs,
    copies and frees (places the information in a completely different memory
    location.

    The most annoying part of this bug is that I am sure to have encountered
    it before. I should have diagnosed it much sooner.


   ...........................................................................

   4/8/97

   SQRT: DOMAIN ERROR observed on hp and on sun. What I did was to
   change most of the sqrt calls, except those embedded into macros,
   and make sure that do not call sqrt on 0 or negative input.

   ...........................................................................
   
   4/9/97

   The "positional_error" component was probably resampled like
   the other components. Then there should have been two components
   called  "positional_error" dep on "positions": the newly created component
   + the resampled component. Apparently there was only one??

   What I did was to add a test to make sure that "positional_error"
   would not be resampled.

   ...........................................................................

   4/10/97

   function 
   Nodetype *_dxfnewnode()

   I wasnt checking that what was returned was different from NULL
   before doing

   Nodetype *retval = _dxfgetnode(table);

   retval->code     = code;
   retval->rec      = rec;

   now it reads 
   if (retval)
    {
      retval->code     = code;
      retval->rec      = rec;
    }
    return(retval);

    I hope that it will get rid of the "Memory Corrupt" signals 
  
   ...........................................................................

   4/11/97 improvement rather than bug fix

   I added a "data" switch in order to turn on or off the use of
   vertex related data to decide whether to collapse edges
   ...........................................................................

   4/12/97-4/13/97 improvement rather than bug fix

   I added a "stats" switch in order to provide information on
   how many vertices and triangles the 
   module starts with and how many vertices and triangles the module ends up
   with


   ...........................................................................


   4/15/97 Lloydt reported that after SimplifySurface, various attributes
   of dx fields were lost

   --> call DXCopyAttributes() after the components are created

   when printing the objects resulting from SimplifySurface, I noticed
   that sometimes the Bounding Box and Neighbors were created and sometimes
   they werent.

   --> call DXEndField() after the field is created (after the worker 
   routine has completed


   
   ...........................................................................


   4/30/97 Lloydt noted that SimplifySurface passes "quads" through
   without processing them

   --> implemented a _dxfQuads2Triangles internal conversion routine
   from quads to triangles. The same connection dependent
   data mapping can be used as before using the "face_lut" array.

   

   ...........................................................................

   4/30/97 treat the following classes in the doLeaf() handler

    case CLASS_PRIVATE:
    case CLASS_LIGHT:
    case CLASS_CAMERA:

    (before, I would have set an error "ERROR_BAD_CLASS", "encountered in object traversal");

   ...........................................................................
 
   4/30/97 series are handled better

   --> set is_series = (DXGetGroupClass((Group)in[0]) == CLASS_SERIES);


   --> use

   new_in[0] = DXGetSeriesMember((Series)in[0], i, &position);
 
   DXSetSeriesMember((Series)new_group, i, (double) position, new_out[0]);

   ...........................................................................

   note: most of the 4/29/96 changes were done during the vacation period
   4/15--4/28. They were checked on 4/30/96

   ...........................................................................

   5/1/97 

   I discovered that the "memory corrupt" signals are set by the DXFree routine
   when it frees a portion of memory that was already freed.

   in the routine _dxfBuildOrientedManifold I had the following piece of code:

   if (*vlut) 
    {DXFree((Pointer)(*vlut));}
 
    *vlut = new_vlut; 

    be careful, here, in case there is an error, new_vlut will be freed


   I had to add the following after the error: flag

   if (new_vlut) 
    {
      DXFree((Pointer)new_vlut); 
      new_vlut was computed and vlut freed and replaced with new_vlut,
      so we must set *vlut = NULL 
      *vlut = NULL;}


      also I redid the routines _dxfCreateSimplificationDataStructure
      _dxfFreeSimplificationDataStructure 
      _dxfPartialFreeSimplificationDataStructure to make sure that pointers were set
     to NULL after being freed.

     _dxfPartialFreeSimplificationDataStructure was created to make space and to avoid running
     out of memory in _dxfBuildSimplifiedSurface

     also. In
     _dxfBuildSimplifiedSurface I made sure pointers were set to NULL after being
     freed, to avoid freeing them twice in the error: part of the code

     the building of new surface vertices and new surface triangles was separated
     for clarity and for freeing un-used memory before allocating new memory

   ...........................................................................

      5/7/97 bug of Floating Point exception noticed on DEC

      in the function  _dxfSimplifiedBoundaryVertexLocation()

      around lines 3700-3800 of _simplifysurface.c

	   A special case is encountered if c == 0,
	   then we also have x = 0 and the formula is simply ye=a
	   this case is handled with the general case without particular
	   attention, except that we can't divide by c and
	   also we should not execute the  portion of code relative to 
	   t == 0 or t == 1
 
   ...........................................................................
   */ 


#include <stdio.h>
#include <math.h>
#ifdef sgi
#include <values.h>
#endif

#ifdef sgi
#include  <sys/bsd_types.h>
#endif

#include  <sys/types.h>

/* because they are not standard accross platforms
   (unless the include files were not completely right)
   I have had to re-typedef the following standard definitions at some point */
/*
typedef unsigned int   uint;
typedef unsigned int   u_int;
typedef unsigned short ushort;
typedef unsigned short u_short;
typedef unsigned char  u_char;
*/

/* define minimum and maximum of two numbers */

#undef  MIN
#define MIN(a,b) ((a) < (b) ? (a) : (b))

#undef  MAX
#define MAX(a,b) ((a) > (b) ? (a) : (b))


/* some vector operations: */

#undef  SQUARED_DISTANCE3
#define SQUARED_DISTANCE3(u,v) (((u)[0]-(v)[0])*((u)[0]-(v)[0]) + ((u)[1]-(v)[1])*((u)[1]-(v)[1]) +\
                                ((u)[2]-(v)[2])*((u)[2]-(v)[2]))


     /* At some point, I might wast to make sure that sqrt is called on non-zero input.
	I will have to replace the macro with a function call */
	
#undef  DIST3
#define DIST3(u,v)             ( sqrt ( SQUARED_DISTANCE3((u),(v))))

#undef  SCALPROD3
#define SCALPROD3(u,v)         ( (u)[0]*(v)[0] + (u)[1]*(v)[1] + (u)[2]*(v)[2] )

#undef  NORM3
#define NORM3(u)               ( sqrt ( SCALPROD3((u),(u)) ) )

#undef  MakeVec
#define MakeVec(u,v,w) {(w)[0] = (v)[0] - (u)[0]; (w)[1] = (v)[1] - (u)[1]; (w)[2] = (v)[2] - (u)[2];}

#undef  AddVec
#define AddVec(u,v)    {/* u = u + v */ (u)[0] += (v)[0]; (u)[1] += (v)[1]; (u)[2] += (v)[2];}

#undef  VecCopy
#define VecCopy(u,v)   {(v)[0] = (u)[0]; (v)[1] = (u)[1]; (v)[2] = (u)[2];}

#undef  VecProd3
#define VecProd3(u,v,w) { (w)[0] = (u)[1]*(v)[2]-(u)[2]*(v)[1]; (w)[1] = (u)[2]*(v)[0]-(u)[0]*(v)[2]; \
                          (w)[2] = (u)[0]*(v)[1]-(u)[1]*(v)[0];}

#undef  EdgeBarycentricCoordinates
#define EdgeBarycentricCoordinates(u,v,a,b) { a = SCALPROD3(u,v); b = SCALPROD3(u,u); }

/* prototype of the driver routine for simplification */

extern int _dxfSimplifySurfaceDriver(int nV, float *v, int nT, int *t,
				     int dim_data,
				     float *vertex_data,
				     float tolerance, 
				     float data_tolerance,
				     /* to be able to resample data defined on a non-manifold
					surface, we need look-up tables between vertices before
					and after manifold conversion and between triangles before 
					and after manifold conversion */
				     int *old_nV_after_manifold_conversion,
				     int *old_nT_after_manifold_conversion,
				     int *new_nV,             /* new number of vertices after simplification*/
				     float **new_v,           /* new vertex array (positions) after 
								 simplification */
				     float **new_vertex_data, /* new data values after simplification */
				     int *new_nT,             /* new number of triangles after simplification*/
				     int   **new_t,           /* new triangle array (connections) after
								 simplification */
				     int   **vertex_parents,  /* relationship between simplified vertices
								 and manifold converted vertices */
				     int   **face_parents,    /* relationship between simplified triangles
								 and manifold converted triangles */
				     int   **vertex_lut,      /* for conversion from non-manifold */
				     int   **face_lut,        /* to manifold surfaces */

				     float *old_positional_error, 

     /* in case we resimplify a surface that was already simplified,
	we use the old positional error as a starting point. The old positional
	error may also have been provided by uncertainty measurements */
								     
				     float **new_positional_error, 
				     float **face_normals,   /* normals of the simplified faces */
				     float **old_face_areas, /* face areas after conversion to manifold */
				     float **new_face_areas, /* face areas after simplification */
				     int preserve_volume,    /* flag = 1 is the volume should be preserved */
				     int simplify_boundary,  /* flag = 1 if the boundary should be
								simplified */
				     /* flag = 1 if the boundary length should be preserved, in case
					the boundary is simplified (simplify_boundary = 1) */
				     int preserve_boundary_length);


				    
extern float _dxfFloatDataBoundingBoxDiagonal(float *data, int n_data, int data_dim);
/* compute the bounding box of any multidimensional array of data */
    

extern int   _dxfVertexNormalsfromTriangleNormals(Array nor, int nV, int nT, int *triangles,
						  float *face_normals, float *face_areas);
/* convert the triangle normals to vertex normals */

/* resample components dep positions and dep connections after simplification 
   "resample" means: 1) create a new component with the same name for the simplified surface
                     2) compute new value at each new vertex or at each new triangle
		     by averaging the values at the corresponding old vertices and the corresponding
		     old triangles:
		     
		     the correspondences between old elements and new elements are obtained by
		     querying the following arrays:  *vertex_parents, *face_parents, *vertex_lut, 
		     *face_lut

		     the array "vertex_parents" gives for each old vertex after conversion to manifold
		     the new vertex to which it maps
		     the array "vertex_lut" gives for each vertex after conversion the corresponding
		     vertex before conversion.
 
		     the array "face_parents" gives for each old triangle after conversion to manifold
		     the new triangle to which it maps
		     the array "face_lut" gives for each triangle after conversion the corresponding
		     triangle before conversion.

		     When averaging vertex attributes, each vertex has weight 1
		     When averaging triangle attributes, each triangle has a weight equal to its area

		     the following components are not resampled
		     POSITIONS, CONNECTIONS, INVALID_POSITIONS INVALID_CONNECTIONS
		     NORMALS, NEIGHBORS, POSITIONAL_ERROR
		     and DATA if it is used by SimplifySurface.

   */

extern int   _dxfResampleComponentValuesAfterSimplification(Field original_surface, 
							    Field *simplified_surface,
							    int old_nV_after_conversion,
							    int old_nT_after_conversion,
							    int new_nV,
							    int new_nT,
							    int *vertex_parents, int *face_parents, 
							    int *vertex_lut, int *face_lut, float *face_areas);

/* routine specialized for resampling position dependent attributes */

extern Array _dxfResampleComponentDepPosAfterSimplification(Array array, 
							    int old_nV_after_conversion,
							    int new_nV, int *vertex_parents, 
							    int *vertex_lut);

/* routine specialized for resampling connection dependent attributes */

extern Array _dxfResampleComponentDepConAfterSimplification(Array array, 
							    int old_nT_after_conversion,
							    int new_nT,
							    int *face_parents, int *face_lut,
							    float *face_areas);





/*+--------------------------------------------------------------------------+
  |                                                                          |
  |   _dxfComputeVertexValence                                               |
  |                                                                          |
  +--------------------------------------------------------------------------+*/

/* using the set of triangles, compute the valence of each vertex */

#define _dxfComputeVertexValence(nT, t, valence)                               \
{                                                                              \
  register int i, i1;                                                          \
  register u_char k;                                                           \
                                                                               \
  /* for each triangle, increment the valence of its three vertices */         \
                                                                               \
  for (i=i1=0; i<(nT); i++, i1+=3) for (k=0; k<3; k++) (valence)[(t)[i1+k]]++; \
										 }


/* 
   following the suggestion of Greg Abram

   I am using macros to implement the resampling functions,
   so that they can be used on all the different scalar data types
   without having to copy the code every time.

   A few other macros are defined later on to reduce the number of 
   function calls on some computation-intensive routines
 */ 



/*+--------------------------------------------------------------------------+
  |                                                                          |
  |   _dxfRESAMPLE_DEP_POS_AFTER_SIMPLIFICATION                              |
  |                                                                          |
  +--------------------------------------------------------------------------+*/

#define _dxfRESAMPLE_DEP_POS_AFTER_SIMPLIFICATION(type)            \
{                                                                  \
  register int j;                                                  \
                                                                   \
  type                                                             \
    *old_data = (type *) DXGetArrayData(array),                    \
    *old_data_current,                                             \
    *new_data = (type *) new_array_data;                           \
                                                                   \
  float                                                            \
    *vertex_data_averages_current;                                 \
                                                                   \
  int                                                              \
    total_new_num_elements = new_nV * tensor_size;                 \
                                                                   \
  /*for each data element (including each tensor element)          \
    we compute the average using the vertex lut and vertex parent  \
                                                                   \
    we need a vertex_weights array allocated to zero beforehand    \
                                                                   \
    we need an array of floats to store temporarily the averages   \
    at each vertex                                                 \
    */                                                             \
                                                                   \
                                                                                       \
 for(i=0; i<old_nV_after_conversion; i++) /* this is the old number of vertices        \
					     after the conversion to manifolds */      \
    {                                                                                  \
      int                                                                              \
                                                                                       \
	old_vertex = vertex_lut[i],                                                    \
	new_vertex = vertex_parents[i];                                                \
                                                                                       \
      /* we add one to the vertex weights of the parents */                            \
                                                                                       \
      vertex_weights[new_vertex] += 1;                                                 \
                                                                                       \
      /* we add the tensor_size values to the average */                               \
                                                                                       \
      /* get the right offset in both old data and new data arrays */                  \
                                                                                       \
      old_data_current = old_data + old_vertex * tensor_size;                          \
                                                                                       \
      vertex_data_averages_current = vertex_data_averages + new_vertex * tensor_size;  \
                                                                                       \
      for (j=0;j<tensor_size;j++)                                                      \
                                                                                       \
	vertex_data_averages_current[j] += (float) old_data_current[j];                \
                                                                                       \
    }                                                                                  \
                                                                                       \
  /* then loop on all the averages, and divide them by the vertex_weights */           \
                                                                                       \
  vertex_data_averages_current = vertex_data_averages;                                 \
                                                                                       \
  for(i=0; i< new_nV; i++,vertex_data_averages_current+=tensor_size)                   \
                                                                                       \
    for (j=0;j<tensor_size;j++)                                                        \
                                                                                       \
      vertex_data_averages_current[j] /= vertex_weights[i];                            \
                                                                                       \
  /* finally convert the vertex_data_averages to the requested scalar type */          \
                                                                                       \
                                                                                       \
  for(i=0; i< total_new_num_elements; i++)                                             \
                                                                                       \
    new_data[i] = (type) vertex_data_averages[i];                                      \
											 }


/*+--------------------------------------------------------------------------+
  |                                                                          |
  |   _dxfRESAMPLE_DEP_CON_AFTER_SIMPLIFICATION                              |
  |                                                                          |
  +--------------------------------------------------------------------------+*/

#define _dxfRESAMPLE_DEP_CON_AFTER_SIMPLIFICATION(type)            \
{                                                                  \
  register int j;                                                  \
                                                                   \
  type                                                             \
    *old_data = (type *) DXGetArrayData(array),                    \
    *old_data_current = old_data,                                  \
    *new_data = (type *) new_array_data;                           \
                                                                   \
  float                                                            \
    *face_data_averages_current;                                   \
                                                                   \
  int                                                              \
    total_new_num_elements = new_nT * tensor_size;                 \
                                                                   \
  for(i=0; i< old_nT_after_conversion; i++, old_data_current += tensor_size)           \
      /* this is the old number of faces after conversion and before simplification */ \
    {                                                                                  \
      int                                                                              \
	old_face = face_lut[i],                                                        \
	new_face = face_parents[i];                                                    \
                                                                                       \
      /* we add the area of the face to the weight of the parent triangle */           \
                                                                                       \
      face_data_weights[new_face] += old_face_areas[old_face];                         \
                                                                                       \
      /* we next update the averages for all associated tensor_size values */          \
                                                                                       \
      /* get the right offset in  both old data and new data arrays */                 \
                                                                                       \
      old_data_current = old_data + old_face * tensor_size;                            \
                                                                                       \
      face_data_averages_current = face_data_averages + new_face * tensor_size;        \
                                                                                       \
      for (j=0;j<tensor_size;j++)                                                      \
                                                                                       \
	face_data_averages_current[j] += old_face_areas[old_face] * (float)            \
							     old_data_current[j];      \
                                                                                       \
    }                                                                                  \
                                                                                       \
  /* then loop on all the averages, and divide them by the face_weights */             \
                                                                                       \
  face_data_averages_current = face_data_averages;                                     \
                                                                                       \
  for(i=0; i< new_nT; i++,face_data_averages_current+=tensor_size)                     \
                                                                                       \
    for (j=0;j<tensor_size;j++)                                                        \
                                                                                       \
      face_data_averages_current[j] /= face_data_weights[i];                           \
                                                                                       \
  /* finally convert the vertex_data_averages to the requested scalar type */          \
                                                                                       \
                                                                                       \
  for(i=0; i< total_new_num_elements; i++)                                             \
                                                                                       \
    new_data[i] = (type) face_data_averages[i];                                        \
											 }


/* copy an array from one type to another type */

/*+--------------------------------------------------------------------------+
  |                                                                          |
  |   CopyType2Type                                                          |
  |                                                                          |
  +--------------------------------------------------------------------------+*/

#define CopyType2Type(type1, type2, a, b, n)                       \
{                                                                  \
  register int i;                                                  \
                                                                   \
  type1 *a1 = (type1 *) (a);                                       \
                                                                   \
  type2 *b1 = (type2 *) (b);                                       \
                                                                   \
  for (i=0; i<(n) ; i++) b1[i] = (type2) a1[i];}


/* convert type "quad" connections to type "triangles" */

extern int _dxfQuads2Triangles(int *nT, int **connections, int **face_lut);


/* data structures for the surface edges and associated routines */


typedef struct _EdgeS_ {

  short            label;               /* used for a flag or for
					   classifying edges among boundary, regular and singular edges */
  int		    v[2];		/* end points */
  int		    t[2];		/* incident triangles */
} EdgeS;

/*
  surface edges have an orientation:
  v[0] >= v[1] and t[0] is before t[1] in clockwise order
  when we rotate around v0.

                      /|\v1
                     / | \
                    /  |  \              
                   / t0|t1 \             
                   \   |   /
                    \  |  /
		     \ | /
		      \|/v0    v0 >= v1
*/

/* routines for indexing edges in a hash table */

extern int         _dxfEdgeSKey(EdgeS *);
extern int         _dxfEdgeSFilter(EdgeS *, EdgeS *);


#define SurfEdgeLabel(edg)           (edg)->label
#define SurfEdgeGetVertices(edg)     (edg)->v
#define SurfEdgeGetTriangles(edg)    (edg)->t


/* my hash tables:

   they do not have the restriction that no more than
   16 elements should have the same code, that is documented in the DX hash tables: 
   
   if someone gives me a surface with 1,000,000 triangles, it should
   roughly have 1,500,000 edges and I don't know of a simple coding
   function that would never give 16 times the same code to edges (it
   may be possible to design such a function but I don't know how)

   */
   
/*...........................................................................*/

/* Following are includes needed for hash-tables using linear probing.
   In case of a code conflict, data items are chained */

/*...........................................................................*/

typedef char  *RECTYPE ;
typedef char  *KEYTYPE ;

typedef struct _nodetype_{
  int   code;
  RECTYPE rec;
  struct _nodetype_ *next;
} Nodetype;

typedef struct _Block_{
  RECTYPE rec;
  struct _Block_ *next;
} Block;

typedef struct _htable_{
  Nodetype **bucket;
  Block     *block;
  Nodetype  *freenode;

  int nodes_per_block ;
  int size, mask;
  int        search_code;
  Nodetype  *search_node;
  KEYTYPE     search_key;

  int     (*fcode)(KEYTYPE);
  int     (*filter)(RECTYPE, RECTYPE);
  KEYTYPE (*compute_key)(RECTYPE);

} Htable;

/*...........................................................................*/

/* prototypes for H-table procedures 
*/

/*...........................................................................*/

/* add a record to the hash-table */
extern int         _dxfAddRecord(Htable *, RECTYPE, KEYTYPE);
extern KEYTYPE     _dxfIdentifyRecord2Key(RECTYPE rec);


/* initialize the hash-table */                                                            

extern int         _dxfInitHashTable(Htable *, int, int (*fcode)(KEYTYPE), 
				     int (*filter)(RECTYPE,RECTYPE), 
				     KEYTYPE (*compute_key)(RECTYPE));

/* determine whether a particular record is in the hash-table and if yes, find that record 
   meaning, return a pointer to that record*/

extern char        *_dxfFindRecord(Htable *, KEYTYPE);
extern void        _dxfHashTableReset(Htable *);

extern int        _dxfNewBlock(Htable *);

/* once _dxfFindRecord() has been called, retrieve the next record with the same code and
   passing through the same filter() from the hash table */

extern char        *_dxfNextRecord(Htable *);
extern Nodetype    *_dxfgetnode(Htable *);
extern Nodetype    *_dxfnewnode(Htable *,RECTYPE, int);

extern int          _dxfFlog2(int n);/* log base two using integer arithmetic */

/* reset the hash table and free the buckets
   after that, either (1) inithashtable must be called again to refill the hash table with new records or
   completely different data
   
   or (2) free(table) to finish it off
   unless it was defined as Htable table[1] in which case nothing more has to be done */

#define       _dxfFreeHashTable(table) {_dxfHashTableReset(table);if ((table)->bucket) \
  {DXFree((Pointer)(table)->bucket);(table)->bucket = NULL;}}


/* find a particular surface edge in a Hash table indexing edges */

extern EdgeS *_dxfFindEdge(Htable *e, int vA, int vB);



/* convert a surface mesh to an oriented manifold surface mesh
   manifold means that each vertex has a single sheet of triangles touching it
   and each edge has at most two incident triangles */

extern int    _dxfToOrientedManifold(int nV,   float *v,   int nT,   int *t, 
				  int *nVm, float **vm, int *nTm, int **tm, int **vlut, int **flut,
				  int *n_edges, Htable *e_table, EdgeS **e);

/* this is the first operation of conversion to a manifold:
   triangles that are degenerate are eliminated
   a triangle is degenerate if either (1) one or more vertices are less then 0 or more than nV
   or (2) two or three vertices in the triangle are the same
   */

extern int    _dxfEliminateDegenerateTriangles(int nV, int nT, int *t, int *new_nT, int **new_t, int **flut);

/* second operation of conversion to manifold:
   the vertices that are not referenced after eliminating the degenerate triangles are
   eliminated from the list of vertices */

extern int    _dxfEliminateStandaloneVertices(int nV, float *v, int *nVm, float **vm, int nT, int *t, 
					      int **vlut);

/* then build a list of edges that are manifold edges, and join the corners of the triangles
   that share such an edge */

extern int   _dxfJoinTriangleCorners(int nT, int *t, Htable *e_table, int *n_edges, EdgeS **edges, 
				     int *fathers);

extern int    _dxfBuildOrientedManifold(int *nV, float **v, int nT, int *t, int **vlut, 
					int *nE, EdgeS **edges, Htable *e_table, int *fathers);

extern int    _dxfIndexEdges(int *nE, EdgeS **edges, Htable *e_table, int nT, int *t);

extern void   _dxfInitFather(int n, int *father);


/* perform path conpression on the representatives of a forest data structure
   the macro FATHER is more efficient than the _dxfFather(i,f) procedure
   it calls the procedure only if necessary */

#define FATHER(i, f) (((i) == (f)[i]) ? (i) : \
                      ((f)[i] == (f)[(f)[i]]) ? (f)[i] : _dxfFather(i,f))

extern int _dxfFather(int i, int *father);

extern int _dxfJoin(int i, int j, int *father);



/*...........................................................................*/
/* 
	default parameters for simplification
*/
/*...........................................................................*/


#define SIMPLIFY_MAX_ANGLE_NOR      1.5     /* maximum angle between normals of corresponding triangles
					       before and after an edge-collapse */
#define SIMPLIFY_COMPACTNESS_RATIO .2       /* maximum ratio between the minimum compactness
					       in the vertex star after simplification and the
					       minimum compactness in the edge star before simplification */
#define SIMPLIFY_VALENCE_MAX        30      /* maximum valence at a vertex after simplification */
#define SIMPLIFY_VALENCE_MAX1       31
#define SIMPLIFY_VALENCE_MAX4       34

#define FOUR_SQRT_3 6.928203230275508
#define EPS_VERTEX_OUTSIDE_TRIANGLE -1e5   /* threshold on the barycentric coordinate such that if 
					      a barycentric coordinate is less than that the corresponding
					      point is decided to be on the other side of the edge */

/* typedef for some simple geometric entities: */

typedef float  Vertex[3];
typedef float  Plane[4];
typedef double VertexD[3];
typedef int    Face[3];

#define DIRECTION_COUNTER_CLOCKWISE 2
#define DIRECTION_CLOCKWISE         1

/*...........................................................................*/

/*
  data structures for binary heaps
*/

/*...........................................................................*/


typedef struct _Heap_ {
  int		n,last;
  int		*perm;
  int		*invp;
  float		*key;
  
} Heap;

#define LeftSon(i)	2*(i)+1
#define RightSon(i)	2*(i)+2

extern Heap          *_dxfNewHeap(int n);
extern unsigned long  _dxfHeapSize(int n);
extern void           _dxfResetHeap(Heap *heap);
extern int            _dxfHeapAdd(float k, Heap *heap);
extern int            _dxfHeapDelMin(Heap *heap);
extern int            _dxfHeapDelete(int i, Heap *heap);

#define              HeapLength(heap)       (heap)->last
#define              HeapLengthMax(heap)    (heap)->n
#define              HeapFull(heap)         (HeapLength(heap) == HeapLengthMax(heap))
#define              HeapValidIndex(heap, index) \
                           ((index) < HeapLengthMax(heap) && (index) >= 0)
   
     /* you can't remove these elements from the heap  but you can add them in: */

#define              HeapOutsideIndex(heap) HeapLengthMax(heap) + 1
#define              HeapInitialIndex(heap) HeapLengthMax(heap)

  /* you can neither add nor remove these elements */

#define              HeapInvalidIndex(heap) -1


/* internal structure used for simplification: */

typedef struct _SimpData_ {
  
  /* data marked "to be allocated" is managed by the 

     _dxfCreateSimplification....

     _dxfFreeSimplification....

     routines. it is used only internally to the _dxfSimplifySurfaceDriver routine and is not exported

     */

  Heap   
        *edge_heap;                  /* to be allocated */

  int   
        nV, nT, nE, 
        data_dim, 
        *triangles,
        *edge2index,                 /* to be allocated */
        *index2edge,                 /* to be allocated */
        *vertex_father,
        *edge_father,                /* to be allocated */
        *triangle_father,
        *vertex_lut,
        preserve_volume,
	simplify_boundary,
	preserve_boundary_length,
        valence_max,

    /* various counters: */

        num_edg_remaining_4_testing,
        num_edg_collapsed,
        num_edg_weights,
        edg_rejected_4_topology,
        edg_rejected_4_geometry,
        edg_rejected_4_tolerance,
        last_edge_added2heap,     /* edge number for the last edge added to the heap */
        heap_initial_index,
        heap_outside_index;
 
  EdgeS 
        *edge_array;

  Htable 
        *e_table;

  Vertex
        *vert,
        *normal,                   /* to be allocated */
        simplified_vertex;
      
  u_char 
        *boundary_vert;            /* to be allocated */

  u_short 
        *valence;                  /* to be allocated */
  
  float 
        tolerance,
        data_tolerance,
        *err_volume,               /* to be allocated */
        *tol_volume,               /* to be allocated */
        *old_face_areas,
        *area,                     /* to be allocated */
        *compactness,              /* to be allocated */ 
        min_scalprod_nor,
        compactness_ratio,
        *vx_data,                  /* to be allocated */
        *vx_data_error,            /* to be allocated */
        *vx_data_potential_values, /* to be allocated */
                                   /* potential data values at the simplified vertex */
        *vx_old_data,
        *vx_new_data,              /* necessary to compute the discrepancy in data values
				      at the vertices (or corners) of the mapping between
				      original and simplified surfaces */

        vx_data_potential_err;     /* potential data error  at the simplified vertex */

  int   (*get_lowest_weight_edge)(struct _SimpData_ *); /* pointer to a function used
							   for extracting the edge
							   with lowest weight from
							   the heap */

  int   (*add_edge)              (struct _SimpData_ *, int);
  int   (*remove_edge)           (struct _SimpData_ *, int, int);

  /* other buffers used for storing information relative to vertex and edge stars */

  int 
      *vertex_star_buffer;

  Pointer 
      edge_star_buffer_fl,
      edge_star_buffer_vx;
 
} SimpData;



extern int  _dxfRemoveLowestWeightEdgeFromHeap(SimpData *simp_data);

extern int  _dxfAddEdge2Heap(SimpData *simp_data, int the_edge);

extern int  _dxfRemoveEdgeFromHeap(SimpData *simp_data, int index, int the_edge);

extern int  _dxfBuildSimplifiedSurface(SimpData *simp_data, int *new_nV, float **new_v, 
				       int *new_nT, int **new_t,
				       float **new_face_areas, float **face_normals,
				       float **new_positional_error, float **new_vertex_data);

extern int  _dxfSimplifyManifoldSurface(SimpData *simp_data);

extern int  _dxfCreateSimplificationDataStructure(SimpData *simp_data, float *data, float *old_pos_err);

/* free all the simplification data structures: */

extern int  _dxfFreeSimplificationDataStructure(SimpData *simp_data, float *data, float *old_pos_err);

/* free only the simplification data structures that are not exported by the driver routine,
   to avoid running out of memory in _dxfBuildSimplifiedSurface() */

extern int  _dxfPartialFreeSimplificationDataStructure(SimpData *simp_data);

extern int  _dxfTrianglesNormalAreaCompactness( int nT, Face *t, Vertex *v, Vertex *t_normal,
						float *t_area, float *t_compactness);

extern int  _dxfTriangleNormalQR2D(VertexD *tri, double *n);

extern int  _dxfTriangleNormalQR2(Vertex *tri, Vertex n);

extern int  _dxfVectorProductQRD(VertexD x1, VertexD x2, double *n);

extern int  _dxfFlagBoundaryVertices(int nE, EdgeS *edges, u_char *boundary_vert);
 
extern int  _dxfMarkEdgesAdjacent2Boundary(SimpData *simp_data);

extern int  _dxfBuildEdgeHeap(SimpData *simp_data);

extern void _dxfInitArray(char *array, char *data, int n, int size);

extern int  _dxfCollapsibilityTestsBoundaryEdge2ndKind(SimpData *simp_data, 
						       int edge_num, int v0, int v1, int val0, int val1, 
						       int move_simplified_vertex);

extern int  _dxfDirectedParentVertexStar(int parent_vertex, int first_triangle,
					 int valence, int *star, SimpData *simp_data, int direction);

extern int  _dxfRotateParentTriangle(int  parent_vertex, int *vertex_fathers, 
				     int *tri_v, int *tri_v_parents, int *tri_v_rotated);

extern int  _dxfRotateParentEdge(int parent_vertex, int *vertex_fathers,
				 int *edge_v, int *edge_t, int *edge_v_parents, 
				 int *edge_v_rotated, int *edge_t_rotated);

/* macros used for finding the next edge when rotating around a vertex in
   on a simplified surface */



#define _dxfNextParentEdge(oriented_tri, directionCW)                                                     \
	 _dxfFather(                                                                                      \
	  (int) (((EdgeS *)_dxfFindEdge(simp_data->e_table,(oriented_tri)[0],(oriented_tri)[directionCW]))\
		 -simp_data->edge_array), simp_data->edge_father)



#define _dxfNextParentTriangle(oriented_edge_t, directionCW)                                \
             ((directionCW) == 1)? FATHER((oriented_edge_t)[1], simp_data->triangle_father):\
				   FATHER((oriented_edge_t)[0], simp_data->triangle_father)

extern int   _dxfManifoldLinkTest(int *link0, int *link1, int val0, int val1);

extern int   _dxfCmpIntSmallFirst(char *s1, char *s2);

extern float _dxfClosestPointOnEdge(Vertex X, Vertex XP, Vertex A, Vertex B, 
				    float *t);/* barycentric coordinates: t, 1-t */

extern void  _dxfMakeEdgeBarycenter(Vertex v0, Vertex v1, float alpha0, Vertex v);

extern float _dxfNormalize3Vector(Vertex v);

extern int   _dxfSolveQuadraticEquation(double A, double B, double C, 
		 	 	        double *sol1, double *sol2, double eps);

extern int   _dxfBoundaryCollapse2KndGeometricallyOk(SimpData *simp_data, int v0, int v1,
					int *star0, int *vstar0, int val0, Vertex *s_normal0, float *s_area0,
					int direction, float min_scalprod_nor, float compactness_ratio);

extern float _dxfFastCompactness3DTriangle2(Vertex *tri, float *triangle_area);



/*+--------------------------------------------------------------------------+
  |                                                                          |
  |   _dxfTriangleBasisVectors                                               |
  |                                                                          |
  +--------------------------------------------------------------------------+*/



#define  _dxfTriangleBasisVectors(tri, x1, x2, origin)                         \
{                                                                              \
  /* compute two basis vectors of a triangle, such that the first basis vector \
     will correspond to the shortest edge,                                     \
     return the number of the vertex that is chosen as the triangle origin */  \
                                                                               \
                                                                               \
  int the_origin = 0;                                                          \
                                                                               \
  /* I- find the shortest edge of the triangle                                 \
                   2                                                           \
                  /\                                                           \
          edg1   /  \ edg0                                                     \
                /____\                                                         \
               0 edg2 1 */                                                     \
                                                                               \
  register u_char j;                                                           \
  double min_len = SQUARED_DISTANCE3((tri)[0],(tri)[2]);                       \
  u_char shortest_edge = 1;                                                    \
                                                                               \
  for (j=0;j<2;j++)                                                            \
    {                                                                          \
      double len = SQUARED_DISTANCE3((tri)[j],(tri)[j+1]);                     \
                                                                               \
      if (len < min_len)                                                       \
	{                                                                      \
	  min_len = len;                                                       \
	  shortest_edge = (j+2)%3;                                             \
	}                                                                      \
      }                                                                        \
                                                                               \
  the_origin = (shortest_edge + 1) %3;                                         \
                                                                               \
  (origin) = the_origin;                                                       \
                                                                               \
  /* II- assign the smallest edge to one of the vectors x1 and x2 */           \
  {                                                                            \
    u_char                                                                     \
      dest1 =  (shortest_edge+2)%3,                                            \
      dest2 =   shortest_edge;                                                 \
                                                                               \
    register u_char k;                                                         \
                                                                               \
    for(k=0;k<3;k++)                                                           \
      {                                                                        \
	(x1)[k] = (tri)[dest1][k] - (tri)[the_origin][k];                      \
	(x2)[k] = (tri)[dest2][k] - (tri)[the_origin][k];                      \
      }                                                                        \
  }                                                                            \
										 }

extern int _dxfErrorWithinToleranceVBoundary2ndKnd(SimpData *simp_data, int v0, int v1, 
						   int *vstar0, int *vstar1, int val0, int val1);


#define ExtractTriangleFromVertexStar( vstar, i) a_triangle; {                 \
  memcpy(a_triangle[1], simp_data->vert[(vstar)[(i)]],   sizeof(Vertex));      \
  memcpy(a_triangle[2], simp_data->vert[(vstar)[(i)+1]], sizeof(Vertex));}

/* macro for taking the minimum of the potential errors at the pivot vertex,
   also called the *simplified_vertex in my article */



#define UPDATE_ERR_PIVOT                                                       \
{	                                                                       \
              if (potential_err_pivot > largest_err_pivot)                     \
		largest_err_pivot = potential_err_pivot;                       \
	       edge_collapsible = (largest_err_pivot <= tol_pivot);}



#define    _dxfErrorTestDataOnSurfaceInit(data)                                                         \
{                                                                                                       \
 /* initialize the potential error value for the data at the potential simplified vertex location */    \
	       if ((data)->vx_data) {(data)->vx_data_potential_err = 0.0;}}



#define    _dxfErrorTestDataOnSurfaceEnd(data,v)                                                        \
{                                                                                                       \
                                                                                                        \
 /* replace the data error value at the simplified vertex with the potential data error value */        \
                                                                                                        \
    if ((data)->vx_data)                                                                                \
        {(data)->vx_data_error[(v)] = (data)->vx_data_potential_err;                                    \
                                                                                                        \
 /* replace the data value with the potential data value */                                             \
         memcpy((data)->vx_data + (v) * (data)->data_dim, (data)->vx_data_potential_values,             \
                (data)->data_dim * sizeof (float));                                                     \
													  }}


#define _dxfErrorTestDataOnSurfaceInterpolate2(data, v0, v1, t)                                             \
{                                                                                                           \
  if ((data)->vx_data)                                                                                      \
	{                                                                                                   \
	  /* interpolate the data value at the potential simplified vertex */                               \
                                                                                                            \
	  int coordinate = 0;                                                                               \
                                                                                                            \
	  float                                                                                             \
	    *data0 = (data)->vx_data + (v0) * (data)->data_dim,                                             \
	    *data1 = (data)->vx_data + (v1) * (data)->data_dim,                                             \
	    t1     = 1. -(t);                                                                               \
                                                                                                            \
	  for (coordinate = 0; coordinate < (data)->data_dim; coordinate++)                                 \
                                                                                                            \
	    (data)->vx_data_potential_values[coordinate] = data0[coordinate] * (t) + data1[coordinate] * t1;\
	}                                                                                                   \
													      }


#define _dxfErrorTestDataOnSurfaceInterpolate3(data, v0, v1, v2, t0, t1, t2)                              \
{                                                                                                         \
  if ((data)->vx_data)                                                                                    \
	{                                                                                                 \
	  /* interpolate the data value at the potential simplified vertex */                             \
                                                                                                          \
	  int coordinate = 0;                                                                             \
                                                                                                          \
	  float                                                                                           \
	    *data0 = (data)->vx_data + (v0) * (data)->data_dim,                                           \
	    *data1 = (data)->vx_data + (v1) * (data)->data_dim,                                           \
	    *data2 = (data)->vx_data + (v2) * (data)->data_dim;                                           \
                                                                                                          \
	  for (coordinate = 0; coordinate < (data)->data_dim; coordinate++)                               \
                                                                                                          \
	    (data)->vx_data_potential_values[coordinate] =                                                \
                  data0[coordinate] * (t0) + data1[coordinate] * (t1) + data2[coordinate] * (t2);         \
	}                                                                                                 \
													    }

extern int   _dxfErrorTestDataOnSurfaceUpdate(SimpData *simp_data, int new_v, int old_v, 
					      int new_v2, int new_v3,
					      float alpha_sv, float alpha_n_v2, float alpha_n_v3,
					      int old_v1, int old_v2, int old_v3, 
					      float alpha_o_v1, float alpha_o_v2, float alpha_o_v3);

extern float  _dxfSegmentIntersection(Vertex A, Vertex B, Vertex C, Vertex D,
				      float *lambda, float *mu, 
				      int changed_AorB, int changed_C, int changed_D);

extern int    _dxfHouseholderPreMultiplication(float *A, int mA, float *v, int m, int n, float *w );

extern int    _dxfHouseholderPreMultiplicationDbl(double *A, int mA, double *v, int m, int n, double *w );

extern int    _dxfCollapseBoundaryEdge2ndKnd(SimpData *simp_data, int edg_num, int v0, int v1, int val0, 
					     int val1, int *star0, int *star1, int *vstar0, int *vstar1, 
					     int *estar0, int *estar1, Vertex *s_normal0, Vertex *s_normal1, 
					     float *s_area0, float *s_area1);

extern int    _dxfRemoveEdgesFromHeap(int *estar, u_short val, int *edge2index, SimpData *simp_data);

extern int    _dxfReinstateEdgesInHeap(int *estar, u_short val, SimpData *simp_data, 
				       int *num_edg_added);

extern int    _dxfCollapsibilityTestsBoundaryEdge1stKind(SimpData *simp_data, 
							 int edge_num, int v0, int v1, int val0, int val1);

extern int    _dxfManifoldLinkTest1stKnd(int v0, int val0, int *star1, int *link1, int val1, 
					 SimpData *simp_data);

extern int    _dxfBoundaryCollapse1stKndGeometricallyOk(SimpData *simp_data, int v0, int *star0, int *vstar0, 
							int val0, Vertex *s_normal0, float *s_area0,
							float min_scalprod_nor, float compactness_ratio);

extern int    _dxfErrorWithinToleranceVBoundary1stKnd(SimpData *simp_data, int v0, int v1, int *vstar1, 
						      int val1);

extern int    _dxfClosestPointOnTriangle(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *dist);

extern int    _dxfTriangle3DBarycentricCoordinates2(Vertex *tri, Vertex w, Vertex wp, Vertex bary, 
						    float *residual);

extern void   _dxfMakeBarycenter(int nV, Vertex *vv, Vertex bary, float *bary_coord);

extern int    _dxfSolveSymmetric2x2Eqn(double a,  double b, double d, double e, double f,
				       double *x, double *y);


extern int    _dxfCollapseBoundaryEdge1stKnd(SimpData *simp_data, int edg_num, int v0, int v1, int val0, 
					     int val1, int *star1, int *vstar1, int *estar1,
					     Vertex *s_normal1, float *s_area1);

extern int    _dxfCollapseTopologicallyFeasible(SimpData *simp_data, int *vstar0, int *vstar1, 
						u_short val0, u_short val1);

extern void   _dxfBuildParentEdgeStars(EdgeS *edg, int v0, int v1, u_short val0,
				       u_short val1, int *star0, int *star1, SimpData *s);

extern void   _dxfParentVertexStar(int vf, int t0, u_short val, int *star, SimpData *simp_data);

extern void   _dxfApplyTranslation(int nV, Vertex *vv, Vertex T);

extern void   _dxfOppositeVector(Vertex u, Vertex v);

extern void   _dxfCopyFaceNormalsAreasCompactness(Plane *plane, float *s_area, float *s_comp, 
						  int *star0, int *star1, int val0, int val013, 
						  Vertex *t_normal, float *t_area, float *t_comp);

extern float  _dxfMinFloatArray(float *array, int n);

extern void   _dxfSimplifiedVertexLocation(SimpData *simpdata, u_short val0, u_short val1, u_short star_val,
					   Vertex *v, Face   *f, Plane  *p, float *area, float *w, 
					   u_char method);

extern double _dxfStarPlaneEquationsAndVolume(Vertex *star_vert,  Face *star_face, Plane  *star_plane, 
					      float *star_areas, int   star_val);

extern void   _dxfComposeVectors(Vertex u, Vertex v, float lambda, Vertex w);

extern int    _dxfPositionVertexLSConstraints2(Plane *star_plane, u_short star_val, Vertex simplified_vertex);

extern void   _dxfFPlaneForIdenticalVolume(u_short val0, u_short val1, Plane p, Vertex *star_vert, 
					   Face *star_face, float star_volume);

extern int    _dxfCollapseGeometricallyFeasible(Vertex *svert, Face *sface, 
						Plane *splane, float *old_comp, float new_comp_min,
						Vertex *snor0, Vertex *snor1,   
						u_short val0, u_short sval, 
						float min_scalprod, float compactness_ratio);

extern int    _dxfNoFlipOverCheck(u_short first_face, u_short last_face, 
				  Vertex *nor, Plane *plane, float min_scalprod_nor);

extern int    _dxfErrorWithinToleranceV(SimpData *simp_data, Vertex *star_vert, Face *star_face,
					Plane *star_plane, int *star0, int *star1, int val0, int val1);

extern void   _dxfSimplifiedStarNormals(Face *sface, Vertex *svert, Vertex *snor, float *sarea, float *scomp, 
					Vertex simpvert, u_short first_face, u_short last_face);

extern int    _dxfCollapseEdge(SimpData *simp_data, int v0, int v1, int *star0, int *star1, 
			       u_short val0, u_short val1, Vertex *nor0, Vertex *nor1, float *area0, 
			       float *area1, float *comp0, float *comp1);

extern void   _dxfReplaceFaceNormals(int *star, Vertex *new_nor, float *new_area, float *new_comp, 
				     Vertex *nor, float *area, float *comp, u_short val);

#endif  /* _simplifysurface_ */

