#ifndef lint static char rcsid [] = "$RCSfile$ $Revision$ $State$"; #endif /**************************************************************************** * File: rtmp_dlist.c * * * * Copyright 1994 by Loral Advanced Distributed Simulation, Inc. * * * * Loral Advanced Distributed Simulation, Inc. * * 10 Moulton Street * * Cambridge, MA 02238 * * 617-873-1850 * * * * This software was developed by Loral under U. S. Government contracts* * and may be reproduced by or for the U. S. Government pursuant to * * the copyright license under the clause at DFARS 252.227-7013 * * (OCT 1988). * * * * Contents: dlist utilities * * Created: * * Author: oded * * Remarks: * * * ****************************************************************************/ #include "rtmp_tree.h" #include "planner.h" #ifndef _MAXFLOAT #define _MAXFLOAT #define MAXFLOAT ((float)3.40282347e+38) #endif static int debug = 0; static int *count_array; static int Space_Dimension = -1; /************************************/ /********* Dlist Utils *************/ /************************************/ void set_dlist_to (DLIST *dlist, DLIST *other_dlist) { DNODE *node; char *next_elt; int i, length; /* empty the dlist */ dlist_free (dlist); /* set it to the new one */ length = dlist_length (other_dlist); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (other_dlist, &node); next_elt = node->value; dlist_add (dlist, next_elt); } } } /************************************/ /********* FT Utils *************/ /************************************/ /* scale v into result */ void ft_scale (float scale, float *v, float *result) { int i; for (i = 0; i < Space_Dimension; i++) result[i] = scale * v[i]; } void ft_copy (float *v1, float *v2) { int i; for (i = 0; i < Space_Dimension; i++) v2[i] = v1[i]; } /* subtract into v3 */ void ft_sub (float *v1, float *v2, float *v3) { int i; for (i = 0; i < Space_Dimension; i++) v3[i] = v1[i] - v2[i]; } float ft_dot (float *v1, float *v2) { int i; float sum; sum = 0.0f; for (i = 0; i < Space_Dimension; i++) sum += (v1[i] * v2[i]); return (sum); } float ft_mag (float *v1) { float sum, mag; double dsum, dmag; sum = ft_dot (v1, v1); dsum = (double)sum; dmag = sqrt (dsum); mag = (float)dmag; return (mag); } /* normalize into result */ void ft_norm (float *v1, float *result) { float mag, scale; mag = ft_mag (v1); scale = 1.0f / mag; ft_scale (scale, v1, result); } /* return 1 for near */ ft_near (float *v1, float *v2, float tolerance) { float *diff, diff_mag; int flag; diff = (float *)calloc (Space_Dimension, sizeof (float)); ft_sub (v1, v2, diff); diff_mag = ft_mag (diff); flag = (diff_mag < tolerance); free (diff); return (flag); } void ft_read (FILE *fp, float *v1) { int i; for (i = 0; i < Space_Dimension; i++) if (fscanf (fp, "%f", v1 + i) != 1) { fprintf (stderr, "ft_read can't get float\n"); exit (0); } } void ft_dump (float *v1) { int i; for (i = 0; i < Space_Dimension; i++) printf ("%6.2f ", v1[i]); printf ("\n"); } void ft_dump_to_file (FILE *fp, float *v1) { int i; for (i = 0; i < Space_Dimension; i++) fprintf (fp, "%f ", v1[i]); fprintf (fp, "\n"); } /**********************************************/ /************* Marker Utils ***************/ /**********************************************/ /* marked, unmarked, and next_unmarked are revolving marker definitions. except when setting all nodes to next_unmarked, the tree is guaranteed to contain only marked and unmarked nodes. after setting all nodes to next_unmarked: unmarked <= next_unmarked marked <= unmarked next_unmarked <= marked nodes are born unmarked. */ static int Marked; static int Unmarked; static int Next_unmarked; static int Mark_semaphore; void init_markers () { Marked = 0; Unmarked = 1; Next_unmarked = 2; Mark_semaphore = AVAILABLE; } void update_markers () { int hold; hold = Unmarked; Unmarked = Next_unmarked; Next_unmarked = Marked; Marked = hold; } /***** clear markers *****/ /* ensures that face and all ascendants are set to Next_unmarked */ void unmark_all_loop (Nface *face) { DNODE *node; Nface *next_face; int i, length; length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->marker != Next_unmarked) unmark_all_loop (next_face); } } face->marker = Next_unmarked; if (face->new_ascendant != NULL) face->new_ascendant->marker = Next_unmarked; } /* unmarks entire tree .. to ensure marker consistency, input must be top */ unmark_all (Nface *hull) { /* check that input is top */ assert (hull->dimension == -1); /* take mark semaphore */ assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; /* set all to Next_unmarked */ unmark_all_loop (hull); /* redefine markers */ update_markers (); Mark_semaphore = AVAILABLE; return (1); } /******************************************/ /********** Nstats Utils **************/ /******************************************/ void Nstats_init (Nstats *stats, int dimension) { int i; stats->basis = NULL; if (dimension > 0) { stats->basis = (float **)calloc (dimension, sizeof (float *)); for (i = 0; i < dimension; i++) stats->basis[i] = (float *)calloc (Space_Dimension, sizeof (float)); } stats->coords = (float *)calloc (Space_Dimension, sizeof (float)); stats->inside_vec = (float *)calloc (Space_Dimension, sizeof (float)); } void Nstats_copy (Nstats *stats, Nstats *copy, int dimension) { int i; for (i = 0; i < dimension; i++) ft_copy (stats->basis[i], copy->basis[i]); ft_copy (stats->coords, copy->coords); ft_copy (stats->inside_vec, copy->inside_vec); } void Nstats_free (Nstats *stats, int dimension) { int i; if (dimension > 0) { for (i = 0; i < dimension; i++) free (stats->basis[i]); free (stats->basis); } free (stats->coords); free (stats->inside_vec); } /******************************************/ /********** Nface Utils **************/ /******************************************/ void Nface_init (Nface *face, int dimension, Nstats *stats) { int i; face->dimension = dimension; face->stats = stats; /* default cmp_func, null free_func .. freeing done elsewhere */ dlist_init (&face->ascendants, NULL, NULL); if (dimension >= 0) face->descendants = (Nface **)calloc (dimension + 1, sizeof (Nface *)); for (i = 0; i < (dimension + 1); i++) face->descendants[i] = NULL; face->new_ascendant = NULL; dlist_init (&face->facets, NULL, NULL); dlist_init (&face->points, NULL, NULL); face->class = CLASS_UNKNOWN; face->marker = Unmarked; if (dimension >= 0) face->ident = count_array[dimension]++; else face->ident = 0; if (debug > 0) { printf ("creating face %d-%d: \n", dimension, face->ident); if (dimension == 0) ft_dump (stats->coords); } } /* copies all info *except* connections to other faces */ void Nface_copy (Nface *face, Nface *copy) { int i; copy->dimension = face->dimension; copy->stats = (Nstats *)malloc (sizeof (Nstats)); Nstats_init (copy->stats, copy->dimension); Nstats_copy (face->stats, copy->stats, copy->dimension); dlist_init (©->ascendants, NULL, NULL); if (copy->dimension >= 0) copy->descendants = (Nface **)calloc (copy->dimension + 1, sizeof (Nface *)); for (i = 0; i < (copy->dimension + 1); i++) copy->descendants[i] = NULL; copy->new_ascendant = NULL; dlist_init (©->facets, NULL, NULL); dlist_init (©->points, NULL, NULL); copy->class = face->class; copy->marker = face->marker; copy->ident = face->ident; if (debug > 0) { printf ("copying face %d-%d: \n", face->dimension, face->ident); if (copy->dimension == 0) ft_dump (copy->stats->coords); } } /* note that all other Nfaces have to be freed explicitly .. this function just kills the stats and the pointers */ void Nface_free (Nface *face) { if (face->stats != NULL) { Nstats_free (face->stats, face->dimension); free (face->stats); } dlist_free (&face->ascendants); if (face->dimension >= 0) free (face->descendants); dlist_free (&face->facets); dlist_free (&face->points); } void Nface_delete_ascendant (Nface *face, Nface *ascendant) { if (ascendant == face->new_ascendant) { face->new_ascendant = NULL; if (debug > 0) printf ("face %d-%d removes new_ascendant connection to %d-%d\n", face->dimension, face->ident, ascendant->dimension, ascendant->ident); } else { dlist_delete (&face->ascendants, (char *)ascendant); if (debug > 0) printf ("face %d-%d removes ascendant connection to %d-%d\n", face->dimension, face->ident, ascendant->dimension, ascendant->ident); } } /* ensures that ascendants have all been killed, removes self from descendants lists, then frees self */ void Nface_kill (Nface *face) { Nface *descendant; int i, length; /* check that face exists */ if (face == NULL) return; /* ascendants should already be gone */ length = dlist_length (&face->ascendants); assert (length == 0); assert (face->new_ascendant == NULL); /* remove self from ascendants list of all descendants */ for (i = 0; i < (face->dimension + 1); i++) { descendant = face->descendants[i]; if (descendant != NULL) Nface_delete_ascendant (descendant, face); } if (debug > 0) printf ("killing face %d-%d\n", face->dimension, face->ident); /* free self */ Nface_free (face); free (face); } /* a point fails the affine test (returns 0) if the point coords minus the dot of these coords with all basis vectors results in something close to the min-dist-to-origin point of the original face .. the original face has dimension basis vectors .. if the point passes the affine test, the function returns a new basis vector and new min-dist-to-origin point so that a new facet can be constructed */ Nface_perpendicular (Nface *face, Nstats *point_stats, float *basis, float *coords) { int i; float *vec, *result, dot, diff; float **face_basis, *face_coords, *point_coords; /* affine set for top is empty.. can't compute new basis */ assert (face->dimension >= 0); /* init */ vec = (float *)calloc (Space_Dimension, sizeof (float)); result = (float *)calloc (Space_Dimension, sizeof (float)); face_basis = face->stats->basis; face_coords = face->stats->coords; point_coords = point_stats->coords; /* results starts with point coordinates */ ft_copy (point_coords, result); /* for each basis vector, subtract from result that vector scaled by the dot product of the vector and the point coords */ for (i = 0; i < face->dimension; i++) { dot = ft_dot (point_coords, face_basis[i]); ft_scale (dot, face_basis[i], vec); ft_sub (result, vec, result); } /* get remaining vector */ ft_sub (result, face_coords, vec); diff = ft_mag (vec); /* if it is too small, the point fails the affine test */ if (diff < FT_AFFINE_EPS) { free (vec); free (result); return (0); } /* otherwise, normalize it into basis */ ft_norm (vec, basis); /* subtract out basis component of results to get coords */ dot = ft_dot (point_coords, basis); ft_scale (dot, basis, vec); ft_sub (result, vec, coords); free (vec); free (result); return (1); } /* variation that just returns the perpendicular distance to the point given */ Nface_distance (Nface *face, Nstats *point_stats, float *diff) { int i; float *vec, *result, dot; float **face_basis, *face_coords, *point_coords; /* affine set for top is empty.. can't compute new basis */ assert (face->dimension >= 0); /* init */ vec = (float *)calloc (Space_Dimension, sizeof (float)); result = (float *)calloc (Space_Dimension, sizeof (float)); face_basis = face->stats->basis; face_coords = face->stats->coords; point_coords = point_stats->coords; /* results starts with point coordinates */ ft_copy (point_coords, result); /* for each basis vector, subtract from result that vector scaled by the dot product of the vector and the point coords */ for (i = 0; i < face->dimension; i++) { dot = ft_dot (point_coords, face_basis[i]); ft_scale (dot, face_basis[i], vec); ft_sub (result, vec, result); } /* get remaining vector */ ft_sub (result, face_coords, vec); *diff = ft_mag (vec); free (vec); free (result); return (1); } /**********************************************/ /********* Nface_Stack Utils *************/ /**********************************************/ Nface_Stack * new_Nface_Stack () { Nface_Stack *st; st = (Nface_Stack *)malloc (sizeof (Nface_Stack)); st->head = NULL; st->next = NULL; st->bottom = 1; return (st); } /* careful! this frees the stacked Nfaces */ void free_Nface_Stack (Nface_Stack *stack) { if (stack->next != NULL) free_Nface_Stack (stack->next); if (stack->head != NULL) { Nface_kill (stack->head); free (stack->head); } free (stack); } void push_Nface_Stack (Nface_Stack *stack, Nface *value) { Nface_Stack *st = new_Nface_Stack (); st->head = stack->head; stack->head = value; st->next = stack->next; stack->next = st; st->bottom = stack->bottom; stack->bottom = NULL; } /* if stack has bottomed out, pop call returns 1 */ pop_Nface_Stack (Nface_Stack *stack, Nface **value) { int bottom_flag = stack->bottom; *value = stack->head; /* if stack is not empty, do the pop */ if (stack->next != NULL) { stack->head = stack->next->head; stack->bottom = stack->next->bottom; stack->next = stack->next->next; } return (bottom_flag); } read_Nface_Stack (Nface_Stack *stack, Nface **value) { if (value != NULL) *value = stack->head; return (stack->bottom); } /*********************************************/ /********** Classification **************/ /*********************************************/ /***** set all to value *****/ void classify_all_loop (Nface *face, int class) { DNODE *node; Nface *next_face; int i, length; /* check whether we've been there */ if (face->marker == Marked) return; /* mark that we're there now.. update class */ face->marker = Marked; face->class = class; /* update all ascendants */ length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; classify_all_loop (next_face, class); } } } void classify_all (Nface *hull, int class) { unmark_all (hull); /* take mark semaphore */ assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; classify_all_loop (hull, class); Mark_semaphore = AVAILABLE; } /***** find facet class from inside vec *****/ /* is the point beneath or beyond the facet? find vec to point, dot with inside_vec */ void classify_facet (Nface *facet, Nstats *point_stats) { float *vec_to_point, *coords, dot; /* init */ vec_to_point = (float *)calloc (Space_Dimension, sizeof (float)); coords = (float *)calloc (Space_Dimension, sizeof (float)); /* get perpendicular vec in direction of point */ Nface_perpendicular (facet, point_stats, vec_to_point, coords); /* check dot product with facet inside vec */ dot = ft_dot (vec_to_point, facet->stats->inside_vec); if (dot > 0.0f) { facet->class = BENEATH; if (debug > 0) printf ("facet %d-%d is class %d\n", facet->dimension, facet->ident, BENEATH); } else { facet->class = BEYOND; if (debug > 0) printf ("facet %d-%d is class %d\n", facet->dimension, facet->ident, BEYOND); } free (vec_to_point); free (coords); } /***** push classification from facets up the tree *****/ void update_face_class (Nface *face, int class) { switch (face->class) { case CLASS_UNKNOWN: face->class = class; break; case BENEATH: if (class != BENEATH) face->class = MIXED; break; case BEYOND: if (class != BEYOND) face->class = MIXED; break; case MIXED: break; default: fprintf (stderr, "unknown class: %d\n", face->class); exit (0); } } void push_classification_loop (Nface *face, int class) { int i; /* check whether we've been there */ if (face->marker == Marked) return; /* mark that we're there now */ face->marker = Marked; /* update face classification */ update_face_class (face, class); /* push facet up through immediate descendants */ for (i = 0; i < (face->dimension + 1); i++) push_classification_loop (face->descendants[i], class); } void push_classification (Nface *hull, Nface *facet) { int i; unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; for (i = 0; i < (facet->dimension + 1); i++) push_classification_loop (facet->descendants[i], facet->class); Mark_semaphore = AVAILABLE; } /***** dump classification of entire tree ****/ void dump_class_loop (Nface *face) { DNODE *node; Nface *next_face; int i, length; if (face->marker == Marked) return; face->marker = Marked; printf ("face %d-%d is class %d\n", face->dimension, face->ident, face->class); /* loop through direct ascendants */ length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; dump_class_loop (next_face); } } } void dump_classes (Nface *hull) { unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; dump_class_loop (hull); Mark_semaphore = AVAILABLE; } /***** top level classify tree *****/ void classify_hull (Nface *hull, DLIST *facets_dlist, Nstats *point_stats) { DNODE *node; Nface *facet; int i, length; /* initialize class to unknown */ classify_all (hull, CLASS_UNKNOWN); /* loop through facets in dlist */ length = dlist_length (facets_dlist); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { /* find and classify next facet */ node = dlist_next_dnode (facets_dlist, &node); facet = (Nface *)node->value; classify_facet (facet, point_stats); push_classification (hull, facet); } } if (debug > 0) dump_classes (hull); } /*********************************************************/ /********** Attach and Remove Commands **************/ /*********************************************************/ void make_new_ascendant (Nface *face, Nface *ascendant) { face->new_ascendant = ascendant; ascendant->descendants[face->dimension + 1] = face; if (debug > 0) { printf ("%d-%d gets new_ascendant %d-%d\n", face->dimension, face->ident, ascendant->dimension, ascendant->ident); printf ("%d-%d is descendant %d of %d-%d\n", face->dimension, face->ident, face->dimension + 1, ascendant->dimension, ascendant->ident); } } void remove_new_ascendant (Nface *face) { if (face->new_ascendant != NULL) Nface_kill (face->new_ascendant); face->new_ascendant = NULL; } void copy_connection (Nface *face, Nface *ascendant, Nface *copy, Nface *copy_ascendant) { int i; dlist_add (©->ascendants, (char *)copy_ascendant); if (debug > 0) printf ("%d-%d gets ascendant %d-%d\n", copy->dimension, copy->ident, copy_ascendant->dimension, copy_ascendant->ident); for (i = 0; i < (ascendant->dimension + 1); i++) if (ascendant->descendants[i] == face) { copy_ascendant->descendants[i] = copy; if (debug > 0) printf ("%d-%d is descendant %d of %d-%d\n", copy->dimension, copy->ident, i, copy_ascendant->dimension, copy_ascendant->ident); } } /***************************************************/ /********** Build Composite Hull **************/ /***************************************************/ /* fills in already allocated stats.. also does the affine test */ get_ascendant_stats (Nface *face, Nstats *point_stats, Nstats *new_stats) { int i; /* get basis vector.. fails affine test if not found */ if (Nface_perpendicular (face, point_stats, new_stats->basis[face->dimension], new_stats->coords) == NULL) return (0); /* fill in basis vectors for new_stats */ for (i = 0; i < face->dimension; i++) ft_copy (face->stats->basis[i], new_stats->basis[i]); return (1); } build_composite_subtree (Nface *hull, Nface *copy, Nstats *point_stats, int level) { Nface *face, *new_face; Nstats *stats; DNODE *node; int length, i; /* never construct ascendants if at last level */ if (level == (Space_Dimension - 1)) return (1); /* deal with direct ascendants of hull */ length = dlist_length (&hull->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { /* find current face in tree */ node = dlist_next_dnode (&hull->ascendants, &node); face = (Nface *)node->value; /* if it hasn't been processed, process it and its subtree */ if (face->new_ascendant == NULL) { /* get new stats and perform affine test */ stats = (Nstats *)malloc (sizeof (Nstats)); Nstats_init (stats, level + 1); if (get_ascendant_stats (face, point_stats, stats) == NULL) { fprintf (stderr, "point failed affine test\n"); return (0); } new_face = (Nface *)malloc (sizeof (Nface)); Nface_init (new_face, level + 1, stats); /* make new_ascendant connection */ make_new_ascendant (face, new_face); /* continue depth-first */ if (build_composite_subtree (face, new_face, point_stats, level + 1) == NULL) return (0); } /* always copy main tree connection */ copy_connection (hull, face, copy, face->new_ascendant); } } return (1); } build_composite_hull (Nface *hull, Nface *point) { /* duplicate top */ make_new_ascendant (hull, point); /* deal with ascendants of top */ return (build_composite_subtree (hull, point, point->stats, 0)); } /***************************************************/ /********** Prune Composite Hull **************/ /***************************************************/ /***** concave loses new_ascendants *****/ void prune_concave_face (Nface *face) { DNODE *node; Nface *next_face; int i, length; /* deal with direct ascendants */ length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->new_ascendant != NULL) prune_concave_face (next_face); } } /* deal with current face */ remove_new_ascendant (face); } /***** tree loses reflex and ascendants *****/ void prune_reflex_face (Nface *face, DLIST *facets_dlist) { DNODE *node; Nface *next_face; int i, length, flag; /* remove ascendants */ flag = 1; while (flag) { length = dlist_length (&face->ascendants); flag = 0; node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->marker != Marked) { prune_reflex_face (next_face, facets_dlist); flag = 1; break; } } } /* remove current point */ remove_new_ascendant (face); if (face->dimension == (Space_Dimension - 1)) { /* remove facet from facets list of all points in facet */ length = dlist_length (&face->points); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->points, &node); next_face = (Nface *)node->value; dlist_delete (&next_face->facets, (char *)face); if (debug > 0) printf ("removing facet %d-%d from list of point %d-%d\n", face->dimension, face->ident, next_face->dimension, next_face->ident); } } /* remove from facets_dlist */ dlist_delete (facets_dlist, (char *)face); if (debug > 0) printf ("removing facet %d-%d from facets_dlist\n", face->dimension, face->ident); } Nface_kill (face); } /* if p is classified as completely BEYOND a point, that point and all its ascendants are removed */ void prune_hull_loop (Nface *face, DLIST *facets_dlist) { DNODE *node; Nface *next_face; int i, length, flag; /* check whether we've been there */ if (face->marker == Marked) return; /* mark that we're there now */ face->marker = Marked; /* process face according to class */ switch (face->class) { case BEYOND: prune_reflex_face (face, facets_dlist); break; case BENEATH: prune_concave_face (face); break; case MIXED: flag = 1; while (flag) { length = dlist_length (&face->ascendants); flag = 0; node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->marker != Marked) { prune_hull_loop (next_face, facets_dlist); flag = 1; break; } } } break; default: fprintf (stderr, "bad class: %d\n", face->class); exit (0); } } void prune_hull (Nface *hull, DLIST *facets_dlist, Nstats *point_stats) { int length; /* check that we're built up enough to be pruning */ length = dlist_length (facets_dlist); assert (length > Space_Dimension); /* classify all faces in hull */ classify_hull (hull, facets_dlist, point_stats); /* trace through tree, treating faces by class */ unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; prune_hull_loop (hull, facets_dlist); Mark_semaphore = AVAILABLE; } /***************************************************/ /********** Merge Composite Hull **************/ /***************************************************/ /***** upgrade new_ascendants *****/ void upgrade_new_ascendant (Nface *face, DLIST *new_facets) { Nface *new_ascendant; /* init */ new_ascendant = face->new_ascendant; /* add new_ascendant to ascendant list of face */ dlist_add_first (&face->ascendants, (char *)new_ascendant); face->new_ascendant = NULL; if (debug > 0) printf ("upgrading new_ascendant %d-%d of face %d-%d\n", new_ascendant->dimension, new_ascendant->ident, face->dimension, face->ident); /* new_ascendant comes in marked */ new_ascendant->marker = Marked; /* if facet status, add to new_facets list */ if (new_ascendant->dimension == (Space_Dimension - 1)) dlist_add_first (new_facets, (char *)new_ascendant); } void update_composite_hull_loop (Nface *face, DLIST *new_facets) { DNODE *node; Nface *next_face; int i, length; /* deal with ascendants */ length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->marker != Marked) { next_face->marker = Marked; update_composite_hull_loop (next_face, new_facets); } } } /* deal with current point */ if (face->new_ascendant != NULL) upgrade_new_ascendant (face, new_facets); } process_facet_inside_vec (Nface *facet, Nface *hull) { DNODE *node; Nface *point, *best_point; float *coords, dist, max_dist; int i, length; /* loop over points */ length = dlist_length (&hull->ascendants); assert (length > 0); node = NULL; best_point = NULL; max_dist = -1.0f; for (i = 0; i < length; i++) { node = dlist_next_dnode (&hull->ascendants, &node); point = (Nface *)node->value; if (dlist_find (&point->facets, (char *)facet) == NULL) { Nface_distance (facet, point->stats, &dist); if (dist > max_dist) { max_dist = dist; best_point = point; } } } /* check that a point was found */ if (best_point == NULL) { fprintf (stderr, "no good points found for inside_vec computation\n"); exit (0); } if (debug > 0) printf ("facet %d-%d external point found: %d-%d\n", facet->dimension, facet->ident, point->dimension, point->ident); /* compute inside_vec using this point */ coords = (float *)calloc (Space_Dimension, sizeof (float)); if (Nface_perpendicular (facet, best_point->stats, facet->stats->inside_vec, coords) == NULL) { fprintf (stderr, "external point no good\n"); free (coords); exit (0); } free (coords); return (1); } /***** process new facets *****/ void update_point (Nface *point, Nface *facet, int add, DLIST *points) { if (add == 1) { dlist_add_first (&point->facets, (char *)facet); dlist_add (points, (char *)point); if (debug > 0) printf ("adding facet %d-%d to list of point %d-%d\n", facet->dimension, facet->ident, point->dimension, point->ident); } else { dlist_delete (&point->facets, (char *)facet); dlist_add (points, (char *)point); if (debug > 0) printf ("removing facet %d-%d from list of point %d-%d\n", facet->dimension, facet->ident, point->dimension, point->ident); } } /* finds point descendants of face and adds facet to their facet list */ void update_constituent_points_loop (Nface *face, Nface *facet, int add, DLIST *points) { int i; /* check whether we've been there */ if (face->marker == Marked) return; /* mark that we're there now */ face->marker = Marked; /* check whether we've reached a point */ if (face->dimension == 0) update_point (face, facet, add, points); else /* push facet up through immediate descendants */ for (i = 0; i < (face->dimension + 1); i++) update_constituent_points_loop (face->descendants[i], facet, add, points); } /* if add is 1, adds facet to point list.. else removes them */ void update_constituent_points (Nface *facet, Nface *hull, int add, DLIST *points) { unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; update_constituent_points_loop (facet, facet, add, points); Mark_semaphore = AVAILABLE; } /* trace through hull, promoting new_ascendants to ascendant status.. if new_ascendant is of facet dimension: add to facets_dlist update facets list of points compute inside_vec from point not on facet */ update_composite_hull (Nface *hull, DLIST *facets_dlist) { DLIST new_facets, new_points; DNODE *node; Nface *next_facet; int i, prev_length, new_length; /* upgrade new_ascendants and get new_facets list */ dlist_init (&new_facets, NULL, NULL); dlist_init (&new_points, NULL, NULL); unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; update_composite_hull_loop (hull, &new_facets); Mark_semaphore = AVAILABLE; /* process new facets */ prev_length = dlist_length (facets_dlist); new_length = dlist_length (&new_facets); if (prev_length == 0) { if (new_length == 1) { /* if only one facet total, cannot get its inside_vec */ next_facet = (Nface *)dlist_first (&new_facets); dlist_add (facets_dlist, (char *)next_facet); update_constituent_points (next_facet, hull, 1, &new_points); set_dlist_to (&next_facet->points, &new_points); if (debug > 0) printf ("add to facets_dlist: %d-%d\n", next_facet->dimension, next_facet->ident); } else assert (new_length == 0); } else { if (prev_length == 1) { /* finish processing old facet */ assert (new_length == Space_Dimension); next_facet = (Nface *)dlist_first (facets_dlist); if (process_facet_inside_vec (next_facet, hull) == NULL) return (0); } /* process new facets */ node = NULL; for (i = 0; i < new_length; i++) { node = dlist_next_dnode (&new_facets, &node); next_facet = (Nface *)node->value; dlist_add (facets_dlist, (char *)next_facet); dlist_free (&new_points); update_constituent_points (next_facet, hull, 1, &new_points); set_dlist_to (&next_facet->points, &new_points); if (process_facet_inside_vec (next_facet, hull) == NULL) return (0); if (debug > 0) printf ("add to facets_dlist: %d-%d\n", next_facet->dimension, next_facet->ident); } } dlist_free (&new_facets); dlist_free (&new_points); return (1); } /*************************************/ /********** Output **************/ /*************************************/ /***** store faces by level *****/ void get_hull_level_array_loop (Nface *face, DLIST *index_array, int level) { DNODE *node; Nface *next_face; int i, length; /* deal with direct ascendants */ length = dlist_length (&face->ascendants); if (length > 0) { node = NULL; for (i = 0; i < length; i++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; if (next_face->marker != Marked) { next_face->marker = Marked; dlist_add (index_array + level, (char *)next_face); get_hull_level_array_loop (next_face, index_array, level + 1); } } } } get_hull_level_array (Nface *hull, DLIST *index_array) { if (Mark_semaphore == TAKEN) { fprintf (stderr, "can't dump hull now\n"); return (0); } unmark_all (hull); assert (Mark_semaphore == AVAILABLE); Mark_semaphore = TAKEN; get_hull_level_array_loop (hull, index_array, 0); Mark_semaphore = AVAILABLE; return (1); } /***** dump parms *****/ void dump_point_parms (Nface *point) { printf ("%d-%d: ", point->dimension, point->ident); ft_dump (point->stats->coords); } void dump_facet_point_list (Nface *facet) { DNODE *node; Nface *next_face; int j, length; printf ("%d-%d: ", facet->dimension, facet->ident); /* print point list */ length = dlist_length (&facet->points); if (length > 0) { printf ("v:("); node = NULL; for (j = 0; j < (length - 1); j++) { node = dlist_next_dnode (&facet->points, &node); next_face = (Nface *)node->value; printf ("%d-%d, ", next_face->dimension, next_face->ident); } node = dlist_next_dnode (&facet->points, &node); next_face = (Nface *)node->value; printf ("%d-%d)\n", next_face->dimension, next_face->ident); } } void dump_facet_inside_vec (Nface *facet) { printf ("%d-%d: ", facet->dimension, facet->ident); ft_dump (facet->stats->inside_vec); } void dump_array_parms (DLIST *index_array) { DNODE *node; Nface *face; int j, length; /* print points and their coords */ printf ("\n\n%dfaces:\n", 0); length = dlist_length (index_array); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array, &node); face = (Nface *)node->value; dump_point_parms (face); } } /* print facet point lists */ printf ("\n\n%dfaces:\n", Space_Dimension - 1); length = dlist_length (index_array + Space_Dimension - 1); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array + Space_Dimension - 1, &node); face = (Nface *)node->value; dump_facet_point_list (face); } } /* print facet inside_vec */ printf ("\n\n%dfaces:\n", Space_Dimension - 1); length = dlist_length (index_array + Space_Dimension - 1); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array + Space_Dimension - 1, &node); face = (Nface *)node->value; dump_facet_inside_vec (face); } } } void dump_array_parms_into_poly (DLIST *index_array, polygonal *poly) { DNODE *node; Nface *face; int j, length; /* print points and their coords */ /* printf ("\n\n%dfaces:\n", 0); */ length = dlist_length (index_array); if (length > 0) { node = NULL; poly->num_vertices=length; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array, &node); face = (Nface *)node->value; /* dump_point_parms (face); */ poly->vertex[j].x=face->stats->coords[0]; poly->vertex[j].y=face->stats->coords[1]; } } /* print facet point lists */ /* printf ("\n\n%dfaces:\n", Space_Dimension - 1); length = dlist_length (index_array + Space_Dimension - 1); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array + Space_Dimension - 1, &node); face = (Nface *)node->value; dump_facet_point_list (face); } } */ /* print facet inside_vec */ /* printf ("\n\n%dfaces:\n", Space_Dimension - 1); length = dlist_length (index_array + Space_Dimension - 1); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array + Space_Dimension - 1, &node); face = (Nface *)node->value; dump_facet_inside_vec (face); } } */ } /* will not dump any new_ascendants */ dump_hull_parms (Nface *hull) { DLIST *index_array; int i; /* fill index array */ index_array = (DLIST *)calloc (Space_Dimension, sizeof (DLIST)); for (i = 0; i < Space_Dimension; i++) dlist_init (index_array + i, NULL, NULL); if (get_hull_level_array (hull, index_array) == NULL) { free (index_array); return (0); } /* dump index array */ dump_array_parms (index_array); /* cleanup */ for (i = 0; i < Space_Dimension; i++) dlist_free (index_array + i); free (index_array); return (1); } /* will not dump any new_ascendants */ dump_hull_parms_into_poly (Nface *hull,polygonal *poly) { DLIST *index_array; int i; /* fill index array */ index_array = (DLIST *)calloc (Space_Dimension, sizeof (DLIST)); for (i = 0; i < Space_Dimension; i++) dlist_init (index_array + i, NULL, NULL); if (get_hull_level_array (hull, index_array) == NULL) { free (index_array); return (0); } /* dump index array */ dump_array_parms_into_poly (index_array,poly); /* cleanup */ for (i = 0; i < Space_Dimension; i++) dlist_free (index_array + i); free (index_array); return (1); } /***** dump structure, connections *****/ void dump_face_structure (Nface *face, int level) { DNODE *node; Nface *next_face; int i, j, length; printf ("%d-%d: ", face->dimension, face->ident); if (level > 0) { /* print descendant indices */ printf ("d:("); for (i = 0; i < level; i++) printf ("%d-%d, ", face->descendants[i]->dimension, face->descendants[i]->ident); printf ("%d-%d) ", face->descendants[level]->dimension, face->descendants[level]->ident); } /* print ascendant list */ length = dlist_length (&face->ascendants); if (length > 0) { printf ("a:("); node = NULL; for (j = 0; j < (length - 1); j++) { node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; printf ("%d-%d, ", next_face->dimension, next_face->ident); } node = dlist_next_dnode (&face->ascendants, &node); next_face = (Nface *)node->value; printf ("%d-%d)", next_face->dimension, next_face->ident); } printf ("\n"); } void dump_array_structure (DLIST *index_array) { DNODE *node; Nface *face; int i, j, length; for (i = 0; i < Space_Dimension; i++) { printf ("\n\n%dfaces:\n", i); length = dlist_length (index_array + i); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (index_array + i, &node); face = (Nface *)node->value; dump_face_structure (face, i); } } } } /* will not dump any new_ascendants */ dump_hull_structure (Nface *hull) { DLIST *index_array; int i; /* fill index array */ index_array = (DLIST *)calloc (Space_Dimension, sizeof (DLIST)); for (i = 0; i < Space_Dimension; i++) dlist_init (index_array + i, NULL, NULL); if (get_hull_level_array (hull, index_array) == NULL) { free (index_array); return (0); } /* dump index array */ dump_array_structure (index_array); /* cleanup */ for (i = 0; i < Space_Dimension; i++) dlist_free (index_array + i); free (index_array); return (1); } /***** dump brief description of hull *****/ dump_convex_hull (Nface *hull) { DLIST *index_array; int i, length; /* fill index array */ index_array = (DLIST *)calloc (Space_Dimension, sizeof (DLIST)); for (i = 0; i < Space_Dimension; i++) dlist_init (index_array + i, NULL, NULL); if (get_hull_level_array (hull, index_array) == NULL) { free (index_array); return (0); } /* print description */ /* printf ("\nHull contains:\n"); for (i = 0; i < Space_Dimension; i++) { length = dlist_length (index_array + i); printf (" %dfaces: %d\n", i, length); } printf ("\n"); */ /* cleanup */ for (i = 0; i < Space_Dimension; i++) dlist_free (index_array + i); free (index_array); return (1); } /**************************************/ /********** Testing **************/ /**************************************/ /* check whether point is part of the facet point list */ point_on_facet (Nface *point, Nface *facet) { DNODE *node; Nface *face; int j, length; length = dlist_length (&facet->points); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (&facet->points, &node); face = (Nface *)node->value; if (face->ident == point->ident) return (1); } } return (0); } float dist_coords_from_facet (float *coords, Nface *facet) { float *vec, dot; /* init */ vec = (float *)calloc (Space_Dimension, sizeof (float)); /* vec goes from facet to point */ ft_sub (coords, facet->stats->coords, vec); /* vec should be in inside_vec direction */ dot = ft_dot (vec, facet->stats->inside_vec); free (vec); return (dot); } /* positive for points beneath facet */ float dist_point_from_facet (Nface *point, Nface *facet) { float dot; dot = dist_coords_from_facet (point->stats->coords, facet); return (dot); } point_beneath_facet (Nface *point, Nface *facet) { float dist; int flag = 1; if (debug > 0) printf (" facet %d-%d: ", facet->dimension, facet->ident); /* if point is part of facet, we're ok */ if (point_on_facet (point, facet) != NULL) { if (debug > 0) printf ("point on facet\n"); return (1); } /* get distance and check that it's close enough to positive */ dist = dist_point_from_facet (point, facet); if ((dist + FT_HULL_EPS) < 0.0f) { fprintf (stderr, "point %d-%d outside facet %d-%d by %f\n", point->dimension, point->ident, facet->dimension, facet->ident, -1.0f * dist); flag = 0; } if (debug > 0) printf ("dist = %f\n", dist); return (flag); } point_inside_hull (Nface *point, DLIST *facets_dlist) { DNODE *node; Nface *facet; int j, length, flag = 1; length = dlist_length (facets_dlist); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (facets_dlist, &node); facet = (Nface *)node->value; if (point_beneath_facet (point, facet) == NULL) flag = 0; } } return (flag); } /* returns 1 if hull, 0 if not */ check_convex_hull (DLIST *facets_dlist, Nface_Stack *point_stack) { Nface *point; int flag = 1; /* work through the stack */ printf ("\n\n"); for (;;) { /* try to pop point_stack */ if (pop_Nface_Stack (point_stack, &point) != NULL) { if (flag) printf ("\nCaptain, it's a convex hull!\n"); return (flag); } /* output point */ printf ("checking %d-%d: ", point->dimension, point->ident); ft_dump (point->stats->coords); /* check point against dlist */ if (point_inside_hull (point, facets_dlist) == NULL) flag = 0; Nface_kill (point); } } int return_convex_hull (DLIST *facets_dlist, Nface_Stack *point_stack, polygonal *the_hull) { Nface *point; int flag = 1; int i=0; for (;;) { /* try to pop point_stack */ if (pop_Nface_Stack (point_stack, &point) != NULL) { the_hull->num_vertices=i; return (0); } /* output point */ the_hull->vertex[i].x=point->stats->coords[0]; the_hull->vertex[i].y=point->stats->coords[1]; /* check point against dlist */ if (point_inside_hull (point, facets_dlist) == NULL) flag = 0; Nface_kill (point); i++; } } /***********************************************/ /********** Dist From Origin **************/ /***********************************************/ float hull_dist_from_origin (DLIST *facets_dlist) { DNODE *node; Nface *facet; int j, length; float *origin; float dist, min_dist = MAXFLOAT; /* init origin */ origin = (float *)calloc (Space_Dimension, sizeof (float)); for (j = 0; j < Space_Dimension; j++) origin[j] = 0.0f; /* get min dist origin from facets */ length = dlist_length (facets_dlist); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (facets_dlist, &node); facet = (Nface *)node->value; dist = dist_coords_from_facet (origin, facet); if (debug > 0) printf ("facet %d-%d is %f from origin\n", facet->dimension, facet->ident, dist); if (dist < min_dist) min_dist = dist; } } free (origin); return (min_dist); } /****************************************/ /********** Top Level **************/ /****************************************/ /***** make the hull *****/ void find_hull_bottom (Nface *face, Nface **bottom) { Nface *next_face; int length; length = dlist_length (&face->ascendants); if (length > 0) { next_face = (Nface *)dlist_first (&face->ascendants); find_hull_bottom (next_face, bottom); } else *bottom = face; } void find_farthest_point (Nface *face, DLIST *points_dlist, Nface **point) { DNODE *node; Nface *next_point; float dist, max_dist; int j, length; max_dist = -1.0f; length = dlist_length (points_dlist); assert (length > 0); /* loop over points to find best */ node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (points_dlist, &node); next_point = (Nface *)node->value; Nface_distance (face, next_point->stats, &dist); if (dist > max_dist) { max_dist = dist; *point = next_point; } } } void select_point (Nface *hull, DLIST *points_dlist, int hull_dimension, Nface **point) { Nface *bottom; /* if complete or null hull, just pop off next point */ if ((hull_dimension >= Space_Dimension) || (hull_dimension < 0)) { *point = (Nface *)dlist_pop (points_dlist); return; } /* find bottom of hull */ find_hull_bottom (hull, &bottom); /* find point with greatest dist from affine set of bottom */ find_farthest_point (bottom, points_dlist, point); /* remove from points_dlist */ dlist_delete (points_dlist, (char *)*point); } add_to_hull (Nface *hull, DLIST *points_dlist, int hull_dimension, DLIST *facets_dlist) { Nface *point; int length; /* see if we're done */ length = dlist_length (points_dlist); if (length == 0) return (1); /* printf (". "); fflush (stdout); */ /* choose a point */ select_point (hull, points_dlist, hull_dimension, &point); /* build parallel structure */ if (build_composite_hull (hull, point) == NULL) return (0); /* if not yet full dimension, this is current hull */ if (hull_dimension < Space_Dimension) { if (update_composite_hull (hull, facets_dlist) == NULL) return (0); return (add_to_hull (hull, points_dlist, hull_dimension + 1, facets_dlist)); } /* if full dimension, prune hull first */ prune_hull (hull, facets_dlist, point->stats); if (update_composite_hull (hull, facets_dlist) == NULL) return (0); return (add_to_hull (hull, points_dlist, hull_dimension, facets_dlist)); } build_convex_hull (DLIST *points_dlist, Nface *hull, DLIST *facets_dlist, polygonal *poly) { /* init */ init_markers (); /* init hull with dimension -1 and null stats */ Nface_init (hull, -1, NULL); /* start points loop */ if (add_to_hull (hull, points_dlist, -1, facets_dlist) == NULL) return (0); /* print result */ if (debug > 0) { dump_hull_parms (hull); dump_hull_structure (hull); } dump_hull_parms_into_poly (hull,poly); return (1); } /***** free the hull *****/ void clean_hull (Nface *face) { DNODE *node; Nface *next_face; while (dlist_length (&face->ascendants) > 0) { node = face->ascendants.head; next_face = (Nface *)node->value; clean_hull (next_face); } remove_new_ascendant (face); Nface_kill (face); } /***** startup *****/ /* random float in range 0 to max */ float random_float (float maxim) { float r1; int r1_int; r1_int = rand (); r1 = (float)r1_int; r1 = r1 / RAND_MAX; r1 = r1 * maxim; return (r1); } /* point ident is guaranteed to correspond to index into points array */ void initialize_point_stacks (int count, float **points, DLIST *points_dlist, Nface_Stack *check_stack) { Nface *face, *face_copy; Nstats *stats; float shake, maxim; int i, j; for (i = 0; i < count; i++) { /* Nstats allocated individually .. they are split up in hull */ stats = (Nstats *)malloc (sizeof (Nstats)); Nstats_init (stats, 0); for (j = 0; j < Space_Dimension; j++) { /* to shake out symmetries */ maxim = (float)(FT_SHAKE); shake = random_float (maxim); shake = shake - 0.5f * maxim; stats->coords[j] = points[i][j] + shake; } face = (Nface *)malloc (sizeof (Nface)); face_copy = (Nface *)malloc (sizeof (Nface)); Nface_init (face, 0, stats); Nface_copy (face, face_copy); dlist_add (points_dlist, (char *)face); push_Nface_Stack (check_stack, face_copy); } } void initialize_points (char *fname, int *count, DLIST *points_dlist, Nface_Stack *check_stack) { FILE *fp; float **points; int i, j, dimension; /* open file */ if ((fp = fopen (fname, "r")) == NULL) { fprintf (stderr, "can't open %s\n", fname); exit (0); } /* read in dimension, set static parms */ fscanf (fp, "%d", &dimension); Space_Dimension = dimension; count_array = (int *)calloc (Space_Dimension, sizeof (int)); for (i = 0; i < Space_Dimension; i++) count_array[i] = 0; /* init and read in points */ fscanf (fp, "%d", count); points = (float **)calloc (*count, sizeof (float *)); for (i = 0; i < *count; i++) points[i] = (float *)calloc (Space_Dimension, sizeof (float)); for (i = 0; i < *count; i++) for (j = 0; j < Space_Dimension; j++) fscanf (fp, "%f", points[i] + j); /* build Nfaces for points */ initialize_point_stacks (*count, points, points_dlist, check_stack); /* free original points */ for (i = 0; i < *count; i++) free (points[i]); free (points); fclose (fp); } /****************************** to run from program .. *******************************/ void coords_list_add_facet (Nface *facet, DLIST *coords_list) { DNODE *node; Nface *next_point; float *vec; int j, length; /* loop through facet points */ length = dlist_length (&facet->points); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (&facet->points, &node); next_point = (Nface *)node->value; /* add inside_vec of facet to dlist indexed by point->ident */ vec = (float *)calloc (Space_Dimension, sizeof (float)); ft_copy (facet->stats->inside_vec, vec); dlist_add (coords_list + next_point->ident, (char *)vec); } } } /* for each facet, enter inside_vec into dlist of each point that is a member of the facet (ident corresponds to index into points array) */ void build_results (DLIST *facets_dlist, DLIST *coords_list) { DNODE *node; Nface *next_facet; int j, length; /* loop through facets_dlist */ length = dlist_length (facets_dlist); if (length > 0) { node = NULL; for (j = 0; j < length; j++) { node = dlist_next_dnode (facets_dlist, &node); next_facet = (Nface *)node->value; coords_list_add_facet (next_facet, coords_list); } } } int convex_hull (int dimension, int count, float **points, DLIST *coords_list, polygonal *the_hull) { DLIST facets_dlist, points_dlist; Nface *hull; Nface_Stack *check_stack; int i, flag = 1; float dist; /* init parms */ Space_Dimension = dimension; count_array = (int *)calloc (Space_Dimension, sizeof (int)); for (i = 0; i < Space_Dimension; i++) count_array[i] = 0; hull = (Nface *)malloc (sizeof (Nface)); dlist_init (&facets_dlist, NULL, NULL); dlist_init (&points_dlist, NULL, NULL); check_stack = new_Nface_Stack (); /* read and init points */ initialize_point_stacks (count, points, &points_dlist, check_stack); /* compute and dump hull */ if (build_convex_hull (&points_dlist, hull, &facets_dlist,the_hull) == NULL) { clean_hull (hull); dlist_free (&facets_dlist); dlist_free (&points_dlist); free_Nface_Stack (check_stack); return (0); } #if 0 check_convex_hull (&facets_dlist, check_stack); return_convex_hull (&facets_dlist, check_stack, the_hull); /* check hull */ check_convex_hull (&facets_dlist, check_stack); /* print stats */ dump_convex_hull (hull); /* check contains origin */ dist = hull_dist_from_origin (&facets_dlist); if (dist < 0.0f) { printf ("\nHull does NOT contain origin .. by %f\n", dist); flag = 0; } else { printf ("\nHull is %f from origin\n", dist); build_results (&facets_dlist, coords_list); } #endif /* cleanup */ clean_hull (hull); dlist_free (&facets_dlist); dlist_free (&points_dlist); free_Nface_Stack (check_stack); free (count_array); return (flag); } /****************************** to run alone .. input filenames: test.dat.6d *******************************/ #if 0 main (int argc, char *argv[]) { DLIST facets_dlist, points_dlist; Nface *hull; Nface_Stack *check_stack; int count; float dist; /* check inputs */ if (argc < 2) { fprintf (stderr, "need a filename\n"); exit (0); } /* init parms */ hull = (Nface *)malloc (sizeof (Nface)); dlist_init (&facets_dlist, NULL, NULL); dlist_init (&points_dlist, NULL, NULL); check_stack = new_Nface_Stack (); /* read and init points */ initialize_points (argv[1], &count, &points_dlist, check_stack); /* compute and dump hull */ build_convex_hull (&points_dlist, hull, &facets_dlist); /* check hull */ check_convex_hull (&facets_dlist, check_stack); /* print stats */ dump_convex_hull (hull); /* check contains origin */ dist = hull_dist_from_origin (&facets_dlist); if (dist < 0.0f) printf ("\nHull does NOT contain origin .. by %f\n", dist); else printf ("\nHull is %f from origin\n", dist); /* cleanup */ clean_hull (hull); dlist_free (&facets_dlist); dlist_free (&points_dlist); free_Nface_Stack (check_stack); free (count_array); return (1); } #endif