/* */
/* 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 02110-1301, 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;
}