diff options
Diffstat (limited to 'digikam/libs/lprof/cmshull.cpp')
-rw-r--r-- | digikam/libs/lprof/cmshull.cpp | 1480 |
1 files changed, 1480 insertions, 0 deletions
diff --git a/digikam/libs/lprof/cmshull.cpp b/digikam/libs/lprof/cmshull.cpp new file mode 100644 index 0000000..e44ba80 --- /dev/null +++ b/digikam/libs/lprof/cmshull.cpp @@ -0,0 +1,1480 @@ +/* */ +/* Little cms - profiler construction set */ +/* Copyright (C) 1998-2001 Marti Maria <marti@littlecms.com> */ +/* */ +/* THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, */ +/* EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY */ +/* WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. */ +/* */ +/* IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL, */ +/* INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, */ +/* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, */ +/* WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF */ +/* LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE */ +/* OF THIS SOFTWARE. */ +/* */ +/* This file is free software; you can redistribute it and/or modify it */ +/* under the terms of the GNU General Public License as published by */ +/* the Free Software Foundation; either version 2 of the License, or */ +/* (at your option) any later version. */ +/* */ +/* This program is distributed in the hope that it will be useful, but */ +/* WITHOUT ANY WARRANTY; without even the implied warranty of */ +/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */ +/* General Public License for more details. */ +/* */ +/* You should have received a copy of the GNU General Public License */ +/* along with this program; if not, write to the Free Software */ +/* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA. */ +/* */ +/* As a special exception to the GNU General Public License, if you */ +/* distribute this file as part of a program that contains a */ +/* configuration script generated by Autoconf, you may include it under */ +/* the same distribution terms that you use for the rest of that program. */ +/* */ +/* Version 1.09a */ + + +#include "lcmsprf.h" + +/* Convex hull management */ + +LCMSHANDLE cdecl cmsxHullInit(void); +void cdecl cmsxHullDone(LCMSHANDLE hHull); +BOOL cdecl cmsxHullAddPoint(LCMSHANDLE hHull, int x, int y, int z); +BOOL cdecl cmsxHullComputeHull(LCMSHANDLE hHull); +char cdecl cmsxHullCheckpoint(LCMSHANDLE hHull, int x, int y, int z); +BOOL cdecl cmsxHullDumpVRML(LCMSHANDLE hHull, const char* fname); + +/* --------------------------------------------------------------------- */ + + + +/* This method is described in "Computational Geometry in C" Chapter 4. */ +/* */ +/* -------------------------------------------------------------------- */ +/* This code is Copyright 1998 by Joseph O'Rourke. It may be freely */ +/* redistributed in its entirety provided that this copyright notice is */ +/* not removed. */ +/* -------------------------------------------------------------------- */ + +#define SWAP(t,x,y) { t = x; x = y; y = t; } + +#define XFREE(p) +/* if (p) { free ((char *) p); p = NULL; } */ + + +#define ADD( head, p ) if ( head ) { \ + p->Next = head; \ + p->Prev = head->Prev; \ + head->Prev = p; \ + p->Prev->Next = p; \ + } \ + else { \ + head = p; \ + head->Next = head->Prev = p; \ + } + +#define XDELETE( head, p ) if ( head ) { \ + if ( head == head->Next ) \ + head = NULL; \ + else if ( p == head ) \ + head = head->Next; \ + p->Next->Prev = p->Prev; \ + p->Prev->Next = p->Next; \ + XFREE( p ); \ + } + +/* Define Vertex indices. */ +#define X 0 +#define Y 1 +#define Z 2 + +/* Define structures for vertices, edges and faces */ + +typedef struct _vertex_struct VERTEX,FAR *LPVERTEX; +typedef struct _edge_struct EDGE, FAR *LPEDGE; +typedef struct _face_struct FACE, FAR *LPFACE; + + +struct _edge_struct { + + LPFACE AdjFace[2]; + LPVERTEX EndPts[2]; + LPFACE NewFace; /* pointer to incident cone face. */ + BOOL DoDelete; /* T iff Edge should be delete. */ + + LPEDGE Next, Prev; +}; + +struct _face_struct { + + LPEDGE Edge[3]; + LPVERTEX Vertex[3]; + BOOL Visible; /* T iff face Visible from new point. */ + + LPFACE Next, Prev; +}; + +struct _vertex_struct { + + int v[3]; + int vnum; + LPEDGE duplicate; /* pointer to incident cone Edge (or NULL) */ + BOOL onhull; /* T iff point on hull. */ + BOOL mark; /* T iff point already processed. */ + + LPVERTEX Next, Prev; +}; + +/* Define flags */ + +#define ONHULL true +#define REMOVED true +#define VISIBLE true +#define PROCESSED true +#define SAFE 1000000 /* Range of safe coord values. */ + +#define DIM 3 /* Dimension of points */ +typedef int VEC3I[DIM]; /* Type integer point */ + +#define PMAX 10000 /* Max # of pts */ + +typedef struct { + + /* Global variable definitions */ + LPVERTEX vertices; + LPEDGE edges; + LPFACE faces; + + VEC3I Vertices[PMAX]; /* All the points */ + VEC3I Faces[PMAX]; /* Each triangle face is 3 indices */ + VEC3I Box[PMAX][2]; /* Box around each face */ + + + VEC3I bmin, bmax; + int radius; + int vnumCounter; + + int nfaces; + int nvertex; + +} HULL, FAR* LPHULL; + +/* static HULL Global; */ + + +/*--------------------------------------------------------------------- +MakeNullVertex: Makes a Vertex, nulls out fields. +---------------------------------------------------------------------*/ + +static +LPVERTEX MakeNullVertex(LPHULL hull) +{ + LPVERTEX v; + + v = (LPVERTEX) malloc(sizeof(VERTEX)); + if (!v) return NULL; + + v->duplicate = NULL; + v->onhull = !ONHULL; + v->mark = !PROCESSED; + ADD( hull->vertices, v ); + + return v; +} + + + +/*--------------------------------------------------------------------- +MakeNullEdge creates a new cell and initializes all pointers to NULL +and sets all flags to off. It returns a pointer to the empty cell. +---------------------------------------------------------------------*/ +static +LPEDGE MakeNullEdge(LPHULL hull) +{ + LPEDGE e; + + e = (LPEDGE) malloc(sizeof(EDGE)); + if (!e) return NULL; + + e->AdjFace[0] = e->AdjFace[1] = e->NewFace = NULL; + e->EndPts[0] = e->EndPts[1] = NULL; + e->DoDelete = !REMOVED; + ADD( hull->edges, e ); + return e; +} + +/*-------------------------------------------------------------------- +MakeNullFace creates a new face structure and initializes all of its +flags to NULL and sets all the flags to off. It returns a pointer +to the empty cell. +---------------------------------------------------------------------*/ +static +LPFACE MakeNullFace(LPHULL hull) +{ + LPFACE f; + int i; + + f = (LPFACE) malloc(sizeof(FACE)); + if (!f) return NULL; + + for ( i=0; i < 3; ++i ) { + f->Edge[i] = NULL; + f->Vertex[i] = NULL; + } + f->Visible = !VISIBLE; + ADD( hull->faces, f ); + return f; +} + + + +/*--------------------------------------------------------------------- +MakeFace creates a new face structure from three vertices (in ccw +order). It returns a pointer to the face. +---------------------------------------------------------------------*/ +static +LPFACE MakeFace(LPHULL hull, LPVERTEX v0, LPVERTEX v1, LPVERTEX v2, LPFACE fold) +{ + LPFACE f; + LPEDGE e0, e1, e2; + + /* Create edges of the initial triangle. */ + if( !fold ) { + e0 = MakeNullEdge(hull); + e1 = MakeNullEdge(hull); + e2 = MakeNullEdge(hull); + } + else { /* Copy from fold, in reverse order. */ + e0 = fold->Edge[2]; + e1 = fold->Edge[1]; + e2 = fold->Edge[0]; + } + e0->EndPts[0] = v0; e0->EndPts[1] = v1; + e1->EndPts[0] = v1; e1->EndPts[1] = v2; + e2->EndPts[0] = v2; e2->EndPts[1] = v0; + + /* Create face for triangle. */ + f = MakeNullFace(hull); + f->Edge[0] = e0; f->Edge[1] = e1; f->Edge[2] = e2; + f->Vertex[0] = v0; f->Vertex[1] = v1; f->Vertex[2] = v2; + + /* Link edges to face. */ + e0->AdjFace[0] = e1->AdjFace[0] = e2->AdjFace[0] = f; + + return f; +} + +/*--------------------------------------------------------------------- +Collinear checks to see if the three points given are collinear, +by checking to see if each element of the cross product is zero. +---------------------------------------------------------------------*/ +static +BOOL Collinear( LPVERTEX a, LPVERTEX b, LPVERTEX c ) +{ + return + ( c->v[Z] - a->v[Z] ) * ( b->v[Y] - a->v[Y] ) - + ( b->v[Z] - a->v[Z] ) * ( c->v[Y] - a->v[Y] ) == 0 + && ( b->v[Z] - a->v[Z] ) * ( c->v[X] - a->v[X] ) - + ( b->v[X] - a->v[X] ) * ( c->v[Z] - a->v[Z] ) == 0 + && ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) - + ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] ) == 0 ; +} + +/*--------------------------------------------------------------------- +VolumeSign returns the sign of the volume of the tetrahedron determined by f +and p. VolumeSign is +1 iff p is on the negative side of f, +where the positive side is determined by the rh-rule. So the volume +is positive if the ccw normal to f points outside the tetrahedron. +The final fewer-multiplications form is due to Bob Williamson. +---------------------------------------------------------------------*/ +int VolumeSign( LPFACE f, LPVERTEX p ) +{ + double vol; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + ax = f->Vertex[0]->v[X] - p->v[X]; + ay = f->Vertex[0]->v[Y] - p->v[Y]; + az = f->Vertex[0]->v[Z] - p->v[Z]; + bx = f->Vertex[1]->v[X] - p->v[X]; + by = f->Vertex[1]->v[Y] - p->v[Y]; + bz = f->Vertex[1]->v[Z] - p->v[Z]; + cx = f->Vertex[2]->v[X] - p->v[X]; + cy = f->Vertex[2]->v[Y] - p->v[Y]; + cz = f->Vertex[2]->v[Z] - p->v[Z]; + + vol = ax * (by*cz - bz*cy) + + ay * (bz*cx - bx*cz) + + az * (bx*cy - by*cx); + + + /* The volume should be an integer. */ + if ( vol > 0.5 ) return 1; + else if ( vol < -0.5 ) return -1; + else return 0; +} + + + +/*--------------------------------------------------------------------- +CleanEdges runs through the Edge list and cleans up the structure. +If there is a NewFace then it will put that face in place of the +Visible face and NULL out NewFace. It also deletes so marked edges. +---------------------------------------------------------------------*/ +static +void CleanEdges(LPHULL hull) +{ + LPEDGE e; /* Primary index into Edge list. */ + LPEDGE t; /* Temporary Edge pointer. */ + + /* Integrate the NewFace's into the data structure. */ + /* Check every Edge. */ + + e = hull ->edges; + do { + if ( e->NewFace ) { + + if ( e->AdjFace[0]->Visible ) + e->AdjFace[0] = e->NewFace; + else + e->AdjFace[1] = e->NewFace; + + e->NewFace = NULL; + } + + e = e->Next; + + } while ( e != hull ->edges ); + + /* Delete any edges marked for deletion. */ + while ( hull ->edges && hull ->edges->DoDelete ) { + + e = hull ->edges; + + XDELETE( hull ->edges, e ); + } + + e = hull ->edges->Next; + + do { + if ( e->DoDelete ) { + + t = e; + e = e->Next; + XDELETE( hull ->edges, t ); + } + else e = e->Next; + + } while ( e != hull ->edges ); +} + +/*--------------------------------------------------------------------- +CleanFaces runs through the face list and deletes any face marked Visible. +---------------------------------------------------------------------*/ +static +void CleanFaces(LPHULL hull) +{ + LPFACE f; /* Primary pointer into face list. */ + LPFACE t; /* Temporary pointer, for deleting. */ + + + while ( hull ->faces && hull ->faces->Visible ) { + + f = hull ->faces; + XDELETE( hull ->faces, f ); + } + + f = hull ->faces->Next; + + do { + if ( f->Visible ) { + + t = f; + f = f->Next; + XDELETE( hull ->faces, t ); + } + else f = f->Next; + + } while ( f != hull ->faces ); +} + + + +/*--------------------------------------------------------------------- +CleanVertices runs through the Vertex list and deletes the +vertices that are marked as processed but are not incident to any +undeleted edges. +---------------------------------------------------------------------*/ +static +void CleanVertices(LPHULL hull) +{ + LPEDGE e; + LPVERTEX v, t; + + /* Mark all vertices incident to some undeleted Edge as on the hull. */ + + e = hull ->edges; + do { + e->EndPts[0]->onhull = e->EndPts[1]->onhull = ONHULL; + e = e->Next; + + } while (e != hull ->edges); + + + /* Delete all vertices that have been processed but + are not on the hull. */ + + while ( hull ->vertices && hull->vertices->mark && !hull ->vertices->onhull ) { + + v = hull ->vertices; + XDELETE(hull ->vertices, v ); + } + + + v = hull ->vertices->Next; + do { + if (v->mark && !v->onhull ) { + t = v; + v = v->Next; + XDELETE(hull ->vertices, t ) + } + else + v = v->Next; + + } while ( v != hull ->vertices ); + + + /* Reset flags. */ + + v = hull ->vertices; + do { + v->duplicate = NULL; + v->onhull = !ONHULL; + v = v->Next; + + } while (v != hull->vertices ); + +} + + + + +/*--------------------------------------------------------------------- +MakeCcw puts the vertices in the face structure in counterclock wise +order. We want to store the vertices in the same +order as in the Visible face. The third Vertex is always p. + +Although no specific ordering of the edges of a face are used +by the code, the following condition is maintained for each face f: +one of the two endpoints of f->Edge[i] matches f->Vertex[i]. +But note that this does not imply that f->Edge[i] is between +f->Vertex[i] and f->Vertex[(i+1)%3]. (Thanks to Bob Williamson.) +---------------------------------------------------------------------*/ + +static +void MakeCcw(LPFACE f, LPEDGE e, LPVERTEX p) +{ + LPFACE fv; /* The Visible face adjacent to e */ + int i; /* Index of e->endpoint[0] in fv. */ + LPEDGE s; /* Temporary, for swapping */ + + if (e->AdjFace[0]->Visible) + + fv = e->AdjFace[0]; + else + fv = e->AdjFace[1]; + + /* Set Vertex[0] & [1] of f to have the same orientation + as do the corresponding vertices of fv. */ + + for ( i=0; fv->Vertex[i] != e->EndPts[0]; ++i ) + ; + + /* Orient f the same as fv. */ + + if ( fv->Vertex[ (i+1) % 3 ] != e->EndPts[1] ) { + + f->Vertex[0] = e->EndPts[1]; + f->Vertex[1] = e->EndPts[0]; + } + else { + f->Vertex[0] = e->EndPts[0]; + f->Vertex[1] = e->EndPts[1]; + SWAP( s, f->Edge[1], f->Edge[2] ); + } + + /* This swap is tricky. e is Edge[0]. Edge[1] is based on endpt[0], + Edge[2] on endpt[1]. So if e is oriented "forwards," we + need to move Edge[1] to follow [0], because it precedes. */ + + f->Vertex[2] = p; +} + +/*--------------------------------------------------------------------- +MakeConeFace makes a new face and two new edges between the +Edge and the point that are passed to it. It returns a pointer to +the new face. +---------------------------------------------------------------------*/ + +static +LPFACE MakeConeFace(LPHULL hull, LPEDGE e, LPVERTEX p) +{ + LPEDGE new_edge[2]; + LPFACE new_face; + int i, j; + + /* Make two new edges (if don't already exist). */ + + for ( i=0; i < 2; ++i ) + /* If the Edge exists, copy it into new_edge. */ + if ( !( new_edge[i] = e->EndPts[i]->duplicate) ) { + + /* Otherwise (duplicate is NULL), MakeNullEdge. */ + new_edge[i] = MakeNullEdge(hull); + new_edge[i]->EndPts[0] = e->EndPts[i]; + new_edge[i]->EndPts[1] = p; + e->EndPts[i]->duplicate = new_edge[i]; + } + + /* Make the new face. */ + new_face = MakeNullFace(hull); + new_face->Edge[0] = e; + new_face->Edge[1] = new_edge[0]; + new_face->Edge[2] = new_edge[1]; + MakeCcw( new_face, e, p ); + + /* Set the adjacent face pointers. */ + for ( i=0; i < 2; ++i ) + for ( j=0; j < 2; ++j ) + /* Only one NULL link should be set to new_face. */ + if ( !new_edge[i]->AdjFace[j] ) { + new_edge[i]->AdjFace[j] = new_face; + break; + } + + return new_face; +} + + +/*--------------------------------------------------------------------- +AddOne is passed a Vertex. It first determines all faces Visible from +that point. If none are Visible then the point is marked as not +onhull. Next is a loop over edges. If both faces adjacent to an Edge +are Visible, then the Edge is marked for deletion. If just one of the +adjacent faces is Visible then a new face is constructed. +---------------------------------------------------------------------*/ +static +BOOL AddOne(LPHULL hull, LPVERTEX p) +{ + LPFACE f; + LPEDGE e, temp; + int vol; + BOOL vis = false; + + + /* Mark faces Visible from p. */ + f = hull -> faces; + + do { + + vol = VolumeSign(f, p); + + if ( vol < 0 ) { + f->Visible = VISIBLE; + vis = true; + } + + f = f->Next; + + } while ( f != hull ->faces ); + + /* If no faces are Visible from p, then p is inside the hull. */ + + if ( !vis ) { + + p->onhull = !ONHULL; + return false; + } + + /* Mark edges in interior of Visible region for deletion. + Erect a NewFace based on each border Edge. */ + + e = hull ->edges; + + do { + + temp = e->Next; + + if ( e->AdjFace[0]->Visible && e->AdjFace[1]->Visible ) + /* e interior: mark for deletion. */ + e->DoDelete = REMOVED; + + else + if ( e->AdjFace[0]->Visible || e->AdjFace[1]->Visible ) + /* e border: make a new face. */ + e->NewFace = MakeConeFace(hull, e, p ); + + e = temp; + + } while ( e != hull ->edges ); + + return true; +} + + +/*--------------------------------------------------------------------- + DoubleTriangle builds the initial double triangle. It first finds 3 + noncollinear points and makes two faces out of them, in opposite order. + It then finds a fourth point that is not coplanar with that face. The + vertices are stored in the face structure in counterclockwise order so + that the volume between the face and the point is negative. Lastly, the + 3 newfaces to the fourth point are constructed and the data structures + are cleaned up. +---------------------------------------------------------------------*/ + +static +BOOL DoubleTriangle(LPHULL hull) +{ + LPVERTEX v0, v1, v2, v3; + LPFACE f0, f1 = NULL; + int vol; + + /* Find 3 noncollinear points. */ + v0 = hull ->vertices; + while ( Collinear( v0, v0->Next, v0->Next->Next ) ) + if ( ( v0 = v0->Next ) == hull->vertices ) + return false; /* All points are Collinear! */ + + v1 = v0->Next; + v2 = v1->Next; + + /* Mark the vertices as processed. */ + v0->mark = PROCESSED; + v1->mark = PROCESSED; + v2->mark = PROCESSED; + + /* Create the two "twin" faces. */ + f0 = MakeFace(hull, v0, v1, v2, f1 ); + f1 = MakeFace(hull, v2, v1, v0, f0 ); + + /* Link adjacent face fields. */ + f0->Edge[0]->AdjFace[1] = f1; + f0->Edge[1]->AdjFace[1] = f1; + f0->Edge[2]->AdjFace[1] = f1; + f1->Edge[0]->AdjFace[1] = f0; + f1->Edge[1]->AdjFace[1] = f0; + f1->Edge[2]->AdjFace[1] = f0; + + /* Find a fourth, noncoplanar point to form tetrahedron. */ + v3 = v2->Next; + vol = VolumeSign( f0, v3 ); + + while ( !vol ) { + + if ( ( v3 = v3->Next ) == v0 ) + return false; /* All points are coplanar! */ + + vol = VolumeSign( f0, v3 ); + } + + /* Insure that v3 will be the first added. */ + hull ->vertices = v3; + return true; +} + + + +/*--------------------------------------------------------------------- +ConstructHull adds the vertices to the hull one at a time. The hull +vertices are those in the list marked as onhull. +---------------------------------------------------------------------*/ +static +void ConstructHull(LPHULL hull) +{ + LPVERTEX v, vnext; + BOOL changed; /* T if addition changes hull; not used. */ + + v = hull->vertices; + + do { + vnext = v->Next; + + changed = false; + + if (!v->mark ) { + + v->mark = PROCESSED; + changed = AddOne(hull, v ); + + CleanEdges(hull); + CleanFaces(hull); + CleanVertices(hull); + } + + v = vnext; + + } while (v != hull->vertices ); + +} + + + +/*-------------------------------------------------------------------*/ + + +static +void AddVec( VEC3I q, VEC3I ray ) +{ + int i; + + for( i = 0; i < DIM; i++ ) + ray[i] = q[i] + ray[i]; +} + +/*--------------------------------------------------------------------- +a - b ==> c. +---------------------------------------------------------------------*/ +static +void SubVec( VEC3I a, VEC3I b, VEC3I c ) +{ + int i; + + for( i = 0; i < DIM; i++ ) + c[i] = a[i] - b[i]; +} + + +/*--------------------------------------------------------------------- +Returns the dot product of the two input vectors. +---------------------------------------------------------------------*/ +static +double Dot( VEC3I a, LPVEC3 b ) +{ + int i; + double sum = 0.0; + + for( i = 0; i < DIM; i++ ) + sum += a[i] * b->n[i]; + + return sum; +} + +/*--------------------------------------------------------------------- +Compute the cross product of (b-a)x(c-a) and place into N. +---------------------------------------------------------------------*/ +static +void NormalVec( VEC3I a, VEC3I b, VEC3I c, LPVEC3 N ) +{ + N->n[X] = ( c[Z] - a[Z] ) * ( b[Y] - a[Y] ) - + ( b[Z] - a[Z] ) * ( c[Y] - a[Y] ); + N->n[Y] = ( b[Z] - a[Z] ) * ( c[X] - a[X] ) - + ( b[X] - a[X] ) * ( c[Z] - a[Z] ); + N->n[Z] = ( b[X] - a[X] ) * ( c[Y] - a[Y] ) - + ( b[Y] - a[Y] ) * ( c[X] - a[X] ); +} + + + + +static +int InBox( VEC3I q, VEC3I bmin, VEC3I bmax ) +{ + + if( ( bmin[X] <= q[X] ) && ( q[X] <= bmax[X] ) && + ( bmin[Y] <= q[Y] ) && ( q[Y] <= bmax[Y] ) && + ( bmin[Z] <= q[Z] ) && ( q[Z] <= bmax[Z] ) ) + return true; + + return false; +} + + + +/* + This function returns a char: + '0': the segment [ab] does not intersect (completely misses) the + bounding box surrounding the n-th triangle T. It lies + strictly to one side of one of the six supporting planes. + '?': status unknown: the segment may or may not intersect T. +*/ + +static +char BoxTest(LPHULL hull, int n, VEC3I a, VEC3I b) +{ + int i; /* Coordinate index */ + int w; + + for ( i=0; i < DIM; i++ ) { + + w = hull ->Box[ n ][0][i]; /* min: lower left */ + + if ( (a[i] < w) && (b[i] < w) ) return '0'; + + w = hull ->Box[ n ][1][i]; /* max: upper right */ + + if ( (a[i] > w) && (b[i] > w) ) return '0'; + } + + return '?'; +} + + + +/* Return a random ray endpoint */ + +static +void RandomRay( VEC3I ray, int radius ) +{ + double x, y, z, w, t; + + /* Generate a random point on a sphere of radius 1. */ + /* the sphere is sliced at z, and a random point at angle t + generated on the circle of intersection. */ + + z = 2.0 * (double) rand() / RAND_MAX - 1.0; + t = 2.0 * M_PI * (double) rand() / RAND_MAX; + w = sqrt( 1 - z*z ); + x = w * cos( t ); + y = w * sin( t ); + + ray[X] = (int) ( radius * x ); + ray[Y] = (int) ( radius * y ); + ray[Z] = (int) ( radius * z ); +} + + + +static +int ComputeBox(LPHULL hull, int F, VEC3I bmin, VEC3I bmax ) +{ + int i, j; + double radius; + + for( i = 0; i < F; i++ ) + for( j = 0; j < DIM; j++ ) { + + if( hull ->Vertices[i][j] < bmin[j] ) + bmin[j] = hull ->Vertices[i][j]; + + if( hull ->Vertices[i][j] > bmax[j] ) + bmax[j] = hull ->Vertices[i][j]; + } + + radius = sqrt( pow( (double)(bmax[X] - bmin[X]), 2.0 ) + + pow( (double)(bmax[Y] - bmin[Y]), 2.0 ) + + pow( (double)(bmax[Z] - bmin[Z]), 2.0 ) ); + + return (int)( radius +1 ) + 1; +} + + +/*--------------------------------------------------------------------- +Computes N & D and returns index m of largest component. +---------------------------------------------------------------------*/ +static +int PlaneCoeff(LPHULL hull, VEC3I T, LPVEC3 N, double *D ) +{ + int i; + double t; /* Temp storage */ + double biggest = 0.0; /* Largest component of normal vector. */ + int m = 0; /* Index of largest component. */ + + NormalVec(hull ->Vertices[T[0]], hull ->Vertices[T[1]], hull ->Vertices[T[2]], N ); + *D = Dot( hull ->Vertices[T[0]], N ); + + /* Find the largest component of N. */ + for ( i = 0; i < DIM; i++ ) { + t = fabs( N->n[i] ); + if ( t > biggest ) { + biggest = t; + m = i; + } + } + return m; +} + +/*--------------------------------------------------------------------- + 'p': The segment lies wholly within the plane. + 'q': The q endpoint is on the plane (but not 'p'). + 'r': The r endpoint is on the plane (but not 'p'). + '0': The segment lies strictly to one side or the other of the plane. + '1': The segement intersects the plane, and 'p' does not hold. +---------------------------------------------------------------------*/ +static +char SegPlaneInt(LPHULL hull, VEC3I T, VEC3I q, VEC3I r, LPVEC3 p, int *m) +{ + VEC3 N; double D; + VEC3I rq; + double num, denom, t; + int i; + + *m = PlaneCoeff(hull, T, &N, &D ); + num = D - Dot( q, &N ); + SubVec( r, q, rq ); + denom = Dot( rq, &N ); + + if ( denom == 0.0 ) { /* Segment is parallel to plane. */ + if ( num == 0.0 ) /* q is on plane. */ + return 'p'; + else + return '0'; + } + else + t = num / denom; + + for( i = 0; i < DIM; i++ ) + p->n[i] = q[i] + t * ( r[i] - q[i] ); + + if ( (0.0 < t) && (t < 1.0) ) + return '1'; + else if ( num == 0.0 ) /* t == 0 */ + return 'q'; + else if ( num == denom ) /* t == 1 */ + return 'r'; + else return '0'; +} + + + +static +int AreaSign( VEC3I a, VEC3I b, VEC3I c ) +{ + double area2; + + area2 = ( b[0] - a[0] ) * (double)( c[1] - a[1] ) - + ( c[0] - a[0] ) * (double)( b[1] - a[1] ); + + /* The area should be an integer. */ + if ( area2 > 0.5 ) return 1; + else if ( area2 < -0.5 ) return -1; + else return 0; +} + + +static +char InTri2D( VEC3I Tp[3], VEC3I pp ) +{ + int area0, area1, area2; + + /* compute three AreaSign() values for pp w.r.t. each Edge of the face in 2D */ + area0 = AreaSign( pp, Tp[0], Tp[1] ); + area1 = AreaSign( pp, Tp[1], Tp[2] ); + area2 = AreaSign( pp, Tp[2], Tp[0] ); + + if ( (( area0 == 0 ) && ( area1 > 0 ) && ( area2 > 0 )) || + (( area1 == 0 ) && ( area0 > 0 ) && ( area2 > 0 )) || + (( area2 == 0 ) && ( area0 > 0 ) && ( area1 > 0 )) ) + return 'E'; + + if ( (( area0 == 0 ) && ( area1 < 0 ) && ( area2 < 0 )) || + (( area1 == 0 ) && ( area0 < 0 ) && ( area2 < 0 )) || + (( area2 == 0 ) && ( area0 < 0 ) && ( area1 < 0 ))) + return 'E'; + + if ( (( area0 > 0 ) && ( area1 > 0 ) && ( area2 > 0 )) || + (( area0 < 0 ) && ( area1 < 0 ) && ( area2 < 0 ))) + return 'F'; + + if ( ( area0 == 0 ) && ( area1 == 0 ) && ( area2 == 0 ) ) + return '?'; /* Error in InTriD */ + + if ( (( area0 == 0 ) && ( area1 == 0 )) || + (( area0 == 0 ) && ( area2 == 0 )) || + (( area1 == 0 ) && ( area2 == 0 )) ) + return 'V'; + + else + return '0'; +} + +/* Assumption: p lies in the plane containing T. + Returns a char: + 'V': the query point p coincides with a Vertex of triangle T. + 'E': the query point p is in the relative interior of an Edge of triangle T. + 'F': the query point p is in the relative interior of a Face of triangle T. + '0': the query point p does not intersect (misses) triangle T. +*/ + +static +char InTri3D(LPHULL hull, VEC3I T, int m, VEC3I p ) +{ + int i; /* Index for X,Y,Z */ + int j; /* Index for X,Y */ + int k; /* Index for triangle Vertex */ + VEC3I pp; /* projected p */ + VEC3I Tp[3]; /* projected T: three new vertices */ + + /* Project out coordinate m in both p and the triangular face */ + j = 0; + for ( i = 0; i < DIM; i++ ) { + if ( i != m ) { /* skip largest coordinate */ + pp[j] = p[i]; + for ( k = 0; k < 3; k++ ) + Tp[k][j] = hull->Vertices[T[k]][i]; + j++; + } + } + return( InTri2D( Tp, pp ) ); +} + + + +static +int VolumeSign2( VEC3I a, VEC3I b, VEC3I c, VEC3I d ) +{ + double vol; + double ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz; + double bxdx, bydy, bzdz, cxdx, cydy, czdz; + + ax = a[X]; + ay = a[Y]; + az = a[Z]; + bx = b[X]; + by = b[Y]; + bz = b[Z]; + cx = c[X]; + cy = c[Y]; + cz = c[Z]; + dx = d[X]; + dy = d[Y]; + dz = d[Z]; + + bxdx=bx-dx; + bydy=by-dy; + bzdz=bz-dz; + cxdx=cx-dx; + cydy=cy-dy; + czdz=cz-dz; + vol = (az-dz) * (bxdx*cydy - bydy*cxdx) + + (ay-dy) * (bzdz*cxdx - bxdx*czdz) + + (ax-dx) * (bydy*czdz - bzdz*cydy); + + + /* The volume should be an integer. */ + if ( vol > 0.5 ) return 1; + else if ( vol < -0.5 ) return -1; + else return 0; +} + + + + +/*--------------------------------------------------------------------- +The signed volumes of three tetrahedra are computed, determined +by the segment qr, and each Edge of the triangle. +Returns a char: + 'v': the open segment includes a Vertex of T. + 'e': the open segment includes a point in the relative interior of an Edge + of T. + 'f': the open segment includes a point in the relative interior of a face + of T. + '0': the open segment does not intersect triangle T. +---------------------------------------------------------------------*/ + +static +char SegTriCross(LPHULL hull, VEC3I T, VEC3I q, VEC3I r ) +{ + int vol0, vol1, vol2; + + vol0 = VolumeSign2( q, hull->Vertices[ T[0] ], hull->Vertices[ T[1] ], r ); + vol1 = VolumeSign2( q, hull->Vertices[ T[1] ], hull->Vertices[ T[2] ], r ); + vol2 = VolumeSign2( q, hull->Vertices[ T[2] ], hull->Vertices[ T[0] ], r ); + + + /* Same sign: segment intersects interior of triangle. */ + if ( ( ( vol0 > 0 ) && ( vol1 > 0 ) && ( vol2 > 0 ) ) || + ( ( vol0 < 0 ) && ( vol1 < 0 ) && ( vol2 < 0 ) ) ) + return 'f'; + + /* Opposite sign: no intersection between segment and triangle */ + if ( ( ( vol0 > 0 ) || ( vol1 > 0 ) || ( vol2 > 0 ) ) && + ( ( vol0 < 0 ) || ( vol1 < 0 ) || ( vol2 < 0 ) ) ) + return '0'; + + else if ( ( vol0 == 0 ) && ( vol1 == 0 ) && ( vol2 == 0 ) ) + return '?'; /* Error 1 in SegTriCross */ + + /* Two zeros: segment intersects Vertex. */ + else if ( ( ( vol0 == 0 ) && ( vol1 == 0 ) ) || + ( ( vol0 == 0 ) && ( vol2 == 0 ) ) || + ( ( vol1 == 0 ) && ( vol2 == 0 ) ) ) + return 'v'; + + /* One zero: segment intersects Edge. */ + else if ( ( vol0 == 0 ) || ( vol1 == 0 ) || ( vol2 == 0 ) ) + return 'e'; + + else + return '?'; /* Error 2 in SegTriCross */ +} + + + +static +char SegTriInt(LPHULL hull, VEC3I T, VEC3I q, VEC3I r, LPVEC3 p ) +{ + int code; + int m = -1; + + code = SegPlaneInt(hull, T, q, r, p, &m ); + + if ( code == '0') return '0'; + else if ( code == 'q') return InTri3D(hull, T, m, q ); + else if ( code == 'r') return InTri3D(hull, T, m, r ); + else if ( code == 'p') return 'p'; + else if ( code == '1' ) return SegTriCross(hull, T, q, r ); + else + return code; /* Error */ +} + + + + +/* + This function returns a char: + 'i': the query point a is strictly interior to polyhedron P. + 'o': the query point a is strictly exterior to( or outside of) polyhedron P. +*/ +char InPolyhedron(LPHULL hull, VEC3I q) +{ + int F = hull->nfaces; + VEC3I Ray; /* Ray endpoint. */ + VEC3 p; /* Intersection point; not used. */ + int f, k = 0, crossings = 0; + char code = '?'; + + + /* If query point is outside bounding box, finished. */ + if ( !InBox( q, hull->bmin, hull->bmax ) ) + return 'o'; + + LOOP: + while( k++ < F ) { + + crossings = 0; + + RandomRay(Ray, hull->radius ); + + AddVec( q, Ray ); + + + for ( f = 0; f < F; f++ ) { /* Begin check each face */ + + if ( BoxTest(hull, f, q, Ray ) == '0' ) { + code = '0'; + + } + else code = SegTriInt(hull, hull->Faces[f], q, Ray, &p ); + + + /* If ray is degenerate, then goto outer while to generate another. */ + if ( code == 'p' || code == 'v' || code == 'e' ) { + + goto LOOP; + } + + /* If ray hits face at interior point, increment crossings. */ + else if ( code == 'f' ) { + crossings++; + + } + + /* If query endpoint q sits on a V/E/F, return inside. */ + else if ( code == 'V' || code == 'E' || code == 'F' ) + return code; /* 'i'; MM2 */ + + /* If ray misses triangle, do nothing. */ + else if ( code == '0' ) + ; + + else + return '?'; /* Error */ + + } + /* No degeneracies encountered: ray is generic, so finished. */ + break; + + } /* End while loop */ + + + /* q strictly interior to polyhedron iff an odd number of crossings. */ + if( ( crossings % 2 ) == 1 ) + return 'i'; + + else return 'o'; +} + + +/*/ ---------------------------------------------------------------------------------- */ + + + + +static +void StoreResults(LPHULL hull) +{ + + int i, w; + LPVERTEX v; + LPFACE f; + int V = 0, F = 0; + int j, k; + + /* Vertices */ + + v = hull ->vertices; + V = 0; + do { + + v -> vnum = V; + hull ->Vertices[V][X] = v -> v[X]; + hull ->Vertices[V][Y] = v -> v[Y]; + hull ->Vertices[V][Z] = v -> v[Z]; + + v = v->Next; + V++; + + } while ( v != hull ->vertices ); + + hull ->nvertex = V; + + /* Faces */ + f = hull ->faces; + F = 0; + do { + + hull ->Faces[F][0] = f->Vertex[0]->vnum; + hull ->Faces[F][1] = f->Vertex[1]->vnum; + hull ->Faces[F][2] = f->Vertex[2]->vnum; + + for ( j=0; j < 3; j++ ) { + + hull ->Box[F][0][j] = hull ->Vertices[ hull ->Faces[F][0] ][j]; + hull ->Box[F][1][j] = hull ->Vertices[ hull ->Faces[F][0] ][j]; + } + + /* Check k=1,2 vertices of face. */ + for ( k=1; k < 3; k++ ) + for ( j=0; j < 3; j++ ) { + + w = hull ->Vertices[ hull ->Faces[F][k] ][j]; + if ( w < hull ->Box[F][0][j] ) hull ->Box[F][0][j] = w; + if ( w > hull ->Box[F][1][j] ) hull ->Box[F][1][j] = w; + } + + + f = f->Next; F++; + + } while ( f != hull ->faces ); + + + hull ->nfaces = F; + + + /* Initialize the bounding box */ + for ( i = 0; i < DIM; i++ ) + hull ->bmin[i] = hull ->bmax[i] = hull ->Vertices[0][i]; + + hull ->radius = ComputeBox(hull, V, hull ->bmin, hull ->bmax ); + + +} + + +LCMSHANDLE cmsxHullInit(void) +{ + LPHULL hull = (LPHULL) malloc(sizeof(HULL)); + + ZeroMemory(hull, sizeof(HULL)); + + hull->vnumCounter = 0; + hull->vertices = NULL; + hull->edges = NULL; + hull->faces = NULL; + hull->nfaces = 0; + hull->nvertex = 0; + + return (LCMSHANDLE) (LPSTR) hull; +} + + +void cmsxHullDone(LCMSHANDLE hHull) +{ + LPHULL hull = (LPHULL) (LPSTR) hHull; + + if (hull) + free((LPVOID) hull); +} + + +BOOL cmsxHullAddPoint(LCMSHANDLE hHull, int x, int y, int z) +{ + LPVERTEX v; + LPHULL hull = (LPHULL) (LPSTR) hHull; + + + v = MakeNullVertex(hull); + v->v[X] = x; + v->v[Y] = y; + v->v[Z] = z; + v->vnum = hull->vnumCounter++; + + return true; +} + +BOOL cmsxHullComputeHull(LCMSHANDLE hHull) +{ + + LPHULL hull = (LPHULL) (LPSTR) hHull; + + if (!DoubleTriangle(hull)) return false; + + ConstructHull(hull); + StoreResults(hull); + + return true; +} + + +char cmsxHullCheckpoint(LCMSHANDLE hHull, int x, int y, int z) +{ + VEC3I q; + LPHULL hull = (LPHULL) (LPSTR) hHull; + + q[X] = x; q[Y] = y; q[Z] = z; + + return InPolyhedron(hull, q ) ; +} + + +BOOL cmsxHullDumpVRML(LCMSHANDLE hHull, const char* fname) +{ + FILE* fp; + int i; + LPHULL hull = (LPHULL) (LPSTR) hHull; + + fp = fopen (fname, "wt"); + if (fp == NULL) + return false; + + fprintf (fp, "#VRML V2.0 utf8\n"); + + /* set the viewing orientation and distance */ + fprintf (fp, "DEF CamTest Group {\n"); + fprintf (fp, "\tchildren [\n"); + fprintf (fp, "\t\tDEF Cameras Group {\n"); + fprintf (fp, "\t\t\tchildren [\n"); + fprintf (fp, "\t\t\t\tDEF DefaultView Viewpoint {\n"); + fprintf (fp, "\t\t\t\t\tposition 0 0 340\n"); + fprintf (fp, "\t\t\t\t\torientation 0 0 1 0\n"); + fprintf (fp, "\t\t\t\t\tdescription \"default view\"\n"); + fprintf (fp, "\t\t\t\t}\n"); + fprintf (fp, "\t\t\t]\n"); + fprintf (fp, "\t\t},\n"); + fprintf (fp, "\t]\n"); + fprintf (fp, "}\n"); + + /* Output the background stuff */ + fprintf (fp, "Background {\n"); + fprintf (fp, "\tskyColor [\n"); + fprintf (fp, "\t\t.5 .5 .5\n"); + fprintf (fp, "\t]\n"); + fprintf (fp, "}\n"); + + /* Output the shape stuff */ + fprintf (fp, "Transform {\n"); + fprintf (fp, "\tscale 8 8 8\n"); + fprintf (fp, "\tchildren [\n"); + + /* Draw the axes as a shape: */ + fprintf (fp, "\t\tShape {\n"); + fprintf (fp, "\t\t\tappearance Appearance {\n"); + fprintf (fp, "\t\t\t\tmaterial Material {\n"); + fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n"); + fprintf (fp, "\t\t\t\t\temissiveColor 1.0 1.0 1.0\n"); + fprintf (fp, "\t\t\t\t\tshininess 0.8\n"); + fprintf (fp, "\t\t\t\t}\n"); + fprintf (fp, "\t\t\t}\n"); + fprintf (fp, "\t\t\tgeometry IndexedLineSet {\n"); + fprintf (fp, "\t\t\t\tcoord Coordinate {\n"); + fprintf (fp, "\t\t\t\t\tpoint [\n"); + fprintf (fp, "\t\t\t\t\t0.0 0.0 0.0,\n"); + fprintf (fp, "\t\t\t\t\t%f 0.0 0.0,\n", 255.0); + fprintf (fp, "\t\t\t\t\t0.0 %f 0.0,\n", 255.0); + fprintf (fp, "\t\t\t\t\t0.0 0.0 %f]\n", 255.0); + fprintf (fp, "\t\t\t\t}\n"); + fprintf (fp, "\t\t\t\tcoordIndex [\n"); + fprintf (fp, "\t\t\t\t\t0, 1, -1\n"); + fprintf (fp, "\t\t\t\t\t0, 2, -1\n"); + fprintf (fp, "\t\t\t\t\t0, 3, -1]\n"); + fprintf (fp, "\t\t\t}\n"); + fprintf (fp, "\t\t}\n"); + + + /* Draw the triangles as a shape: */ + fprintf (fp, "\t\tShape {\n"); + fprintf (fp, "\t\t\tappearance Appearance {\n"); + fprintf (fp, "\t\t\t\tmaterial Material {\n"); + fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n"); + fprintf (fp, "\t\t\t\t\temissiveColor 0 0 0\n"); + fprintf (fp, "\t\t\t\t\tshininess 0.8\n"); + fprintf (fp, "\t\t\t\t}\n"); + fprintf (fp, "\t\t\t}\n"); + fprintf (fp, "\t\t\tgeometry IndexedFaceSet {\n"); + fprintf (fp, "\t\t\t\tsolid false\n"); + + /* fill in the points here */ + fprintf (fp, "\t\t\t\tcoord Coordinate {\n"); + fprintf (fp, "\t\t\t\t\tpoint [\n"); + + for (i = 0; i < hull->nvertex; ++i) + { + fprintf (fp, "\t\t\t\t\t%g %g %g%c\n", + (double) hull->Vertices[i][X], (double) hull->Vertices[i][Y], (double) hull->Vertices[i][Z], + i == hull->nvertex-1? ']': ','); + } + fprintf (fp, "\t\t\t\t}\n"); + + /* fill in the Vertex indices (followed by -1) */ + + + fprintf (fp, "\t\t\t\tcoordIndex [\n"); + for (i = 0; i < hull->nfaces; ++i) + { + fprintf (fp, "\t\t\t\t\t%d, %d, %d, -1\n", + hull->Faces[i][0], hull->Faces[i][1], hull->Faces[i][2]); + + } + fprintf (fp, "]\n"); + + + /* fill in the face colors */ + fprintf (fp, "\t\t\t\tcolor Color {\n"); + fprintf (fp, "\t\t\t\t\tcolor [\n"); + for (i = 0; i < hull->nfaces; ++i) + { + int vx, vy, vz; + double r, g, b; + + vx = hull->Faces[i][0]; vy = hull->Faces[i][1]; vz = hull->Faces[i][2]; + r = (double) (hull->Vertices[vx][X] + hull->Vertices[vy][X] + hull->Vertices[vz][X]) / (3* 255); + g = (double) (hull->Vertices[vx][Y] + hull->Vertices[vy][Y] + hull->Vertices[vz][Y]) / (3* 255); + b = (double) (hull->Vertices[vx][Z] + hull->Vertices[vy][Z] + hull->Vertices[vz][Z]) / (3* 255); + + fprintf (fp, "\t\t\t\t\t%g %g %g%c\n", r, g, b, + i == hull->nfaces-1? ']': ','); + } + fprintf (fp, "\t\t\t}\n"); + + fprintf (fp, "\t\t\tcolorPerVertex false\n"); + + fprintf (fp, "\t\t\t}\n"); + fprintf (fp, "\t\t}\n"); + fprintf (fp, "\t]\n"); + fprintf (fp, "}\n"); + + fclose (fp); + + return true; +} |