Subversion Repositories Scribus

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2454 craig 1
//
2
//  Little cms - profiler construction set
3
//  Copyright (C) 1998-2001 Marti Maria
4
//
5
// THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,
6
// EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
7
// WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
8
//
9
// IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL,
10
// INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
11
// OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
12
// WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF
13
// LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
14
// OF THIS SOFTWARE.
15
//
16
// This file is free software; you can redistribute it and/or modify it
17
// under the terms of the GNU General Public License as published by
18
// the Free Software Foundation; either version 2 of the License, or
19
// (at your option) any later version.
20
//
21
// This program is distributed in the hope that it will be useful, but
22
// WITHOUT ANY WARRANTY; without even the implied warranty of
23
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24
// General Public License for more details.
25
//
26
// You should have received a copy of the GNU General Public License
27
// along with this program; if not, write to the Free Software
28
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
29
//
30
// As a special exception to the GNU General Public License, if you
31
// distribute this file as part of a program that contains a
32
// configuration script generated by Autoconf, you may include it under
33
// the same distribution terms that you use for the rest of that program.
34
//
35
// Version 1.09a
36
 
37
 
38
#include "lcmsprf.h"
39
 
40
// Convex hull management
41
 
42
LCMSHANDLE cdecl cmsxHullInit(void);
43
void       cdecl cmsxHullDone(LCMSHANDLE hHull);
44
BOOL       cdecl cmsxHullAddPoint(LCMSHANDLE hHull, int x, int y, int z);
45
BOOL       cdecl cmsxHullComputeHull(LCMSHANDLE hHull);
46
char       cdecl cmsxHullCheckpoint(LCMSHANDLE hHull, int x, int y, int z);
47
BOOL       cdecl cmsxHullDumpVRML(LCMSHANDLE hHull, const char* fname);
48
 
49
// ---------------------------------------------------------------------
50
 
51
 
52
 
53
// This method is described in "Computational Geometry in C" Chapter 4.  
54
// 
55
// --------------------------------------------------------------------
56
// This code is Copyright 1998 by Joseph O'Rourke.  It may be freely 
57
// redistributed in its entirety provided that this copyright notice is
58
// not removed.
59
// --------------------------------------------------------------------
60
 
61
#define SWAP(t,x,y)     { t = x; x = y; y = t; }
62
 
63
#define XFREE(p)        
64
// if (p) { free ((char *) p); p = NULL; }
65
 
66
 
67
#define ADD( head, p )  if ( head )  { \
68
                                p->Next = head; \
69
                                p->Prev = head->Prev; \
70
                                head->Prev = p; \
71
                                p->Prev->Next = p; \
72
                        } \
73
                        else { \
74
                                head = p; \
75
                                head->Next = head->Prev = p; \
76
                        }
77
 
78
#define XDELETE( head, p ) if ( head )  { \
79
                                if ( head == head->Next ) \
80
                                        head = NULL;  \
81
                                else if ( p == head ) \
82
                                        head = head->Next; \
83
                                p->Next->Prev = p->Prev;  \
84
                                p->Prev->Next = p->Next;  \
85
                                XFREE( p ); \
86
                        } 
87
 
88
// Define Vertex indices. 
89
#define X   0
90
#define Y   1
91
#define Z   2
92
 
93
// Define structures for vertices, edges and faces 
94
 
95
typedef struct _vertex_struct VERTEX,FAR *LPVERTEX;
96
typedef struct _edge_struct EDGE, FAR *LPEDGE;
97
typedef struct _face_struct FACE, FAR *LPFACE;
98
 
99
 
100
struct _edge_struct {
101
 
102
   LPFACE    AdjFace[2];
103
   LPVERTEX  EndPts[2];
104
   LPFACE    NewFace;                       /* pointer to incident cone face. */
105
   BOOL      DoDelete;                          /* T iff Edge should be delete. */
106
 
107
   LPEDGE    Next, Prev;
108
};
109
 
110
struct _face_struct {
111
 
112
   LPEDGE    Edge[3];
113
   LPVERTEX  Vertex[3];
114
   BOOL      Visible;            /* T iff face Visible from new point. */
115
 
116
   LPFACE    Next, Prev;
117
};
118
 
119
struct _vertex_struct {
120
 
121
   int      v[3];
122
   int      vnum;
123
   LPEDGE   duplicate;          /* pointer to incident cone Edge (or NULL) */
124
   BOOL     onhull;                             /* T iff point on hull. */
125
   BOOL     mark;                               /* T iff point already processed. */
126
 
127
   LPVERTEX  Next, Prev;
128
};
129
 
130
// Define flags
131
 
132
#define ONHULL          TRUE
133
#define REMOVED         TRUE
134
#define VISIBLE         TRUE
135
#define PROCESSED       TRUE
136
#define SAFE            1000000         /* Range of safe coord values. */
137
 
138
#define DIM 3                  /* Dimension of points */
139
typedef int    VEC3I[DIM];   /* Type integer point */
140
 
141
#define PMAX 10000       /* Max # of pts */
142
 
143
typedef struct {
144
 
145
        /* Global variable definitions */
146
        LPVERTEX vertices;
147
        LPEDGE edges;
148
        LPFACE faces;
149
 
150
        VEC3I Vertices[PMAX];        /* All the points */
151
        VEC3I Faces[PMAX];           /* Each triangle face is 3 indices */
152
        VEC3I Box[PMAX][2];          /* Box around each face */
153
 
154
 
155
        VEC3I bmin, bmax;
156
        int radius;
157
        int vnumCounter;
158
 
159
        int nfaces;
160
        int nvertex;
161
 
162
} HULL, FAR* LPHULL;
163
 
164
// static HULL Global;
165
 
166
 
167
/*---------------------------------------------------------------------
168
MakeNullVertex: Makes a Vertex, nulls out fields.
169
---------------------------------------------------------------------*/
170
 
171
static
172
LPVERTEX MakeNullVertex(LPHULL hull)
173
{
174
   LPVERTEX  v;
175
 
176
   v = (LPVERTEX) malloc(sizeof(VERTEX));
177
   if (!v) return NULL;
178
 
179
   v->duplicate = NULL;
180
   v->onhull = !ONHULL;
181
   v->mark = !PROCESSED;
182
   ADD( hull->vertices, v );
183
 
184
   return v;
185
}
186
 
187
 
188
 
189
/*---------------------------------------------------------------------
190
MakeNullEdge creates a new cell and initializes all pointers to NULL
191
and sets all flags to off.  It returns a pointer to the empty cell.
192
---------------------------------------------------------------------*/
193
static
194
LPEDGE MakeNullEdge(LPHULL hull)
195
{
196
   LPEDGE  e;
197
 
198
   e = (LPEDGE) malloc(sizeof(EDGE));
199
   if (!e) return NULL;
200
 
201
   e->AdjFace[0] = e->AdjFace[1] = e->NewFace = NULL;
202
   e->EndPts[0] = e->EndPts[1] = NULL;
203
   e->DoDelete = !REMOVED;
204
   ADD( hull->edges, e );
205
   return e;
206
}
207
 
208
/*--------------------------------------------------------------------
209
MakeNullFace creates a new face structure and initializes all of its
210
flags to NULL and sets all the flags to off.  It returns a pointer
211
to the empty cell.
212
---------------------------------------------------------------------*/
213
static
214
LPFACE MakeNullFace(LPHULL hull)
215
{
216
   LPFACE  f;
217
   int    i;
218
 
219
   f = (LPFACE) malloc(sizeof(FACE));
220
   if (!f) return NULL;
221
 
222
   for ( i=0; i < 3; ++i ) {
223
      f->Edge[i] = NULL;
224
      f->Vertex[i] = NULL;
225
   }
226
   f->Visible = !VISIBLE;
227
   ADD( hull->faces, f );
228
   return f;
229
}
230
 
231
 
232
 
233
/*---------------------------------------------------------------------
234
MakeFace creates a new face structure from three vertices (in ccw
235
order).  It returns a pointer to the face.
236
---------------------------------------------------------------------*/
237
static
238
LPFACE MakeFace(LPHULL hull, LPVERTEX v0, LPVERTEX v1, LPVERTEX v2, LPFACE fold)
239
{
240
   LPFACE  f;
241
   LPEDGE  e0, e1, e2;
242
 
243
   /* Create edges of the initial triangle. */
244
   if( !fold ) {
245
     e0 = MakeNullEdge(hull);
246
     e1 = MakeNullEdge(hull);
247
     e2 = MakeNullEdge(hull);
248
   }
249
   else { /* Copy from fold, in reverse order. */
250
     e0 = fold->Edge[2];
251
     e1 = fold->Edge[1];
252
     e2 = fold->Edge[0];
253
   }
254
   e0->EndPts[0] = v0; e0->EndPts[1] = v1;
255
   e1->EndPts[0] = v1; e1->EndPts[1] = v2;
256
   e2->EndPts[0] = v2; e2->EndPts[1] = v0;
257
 
258
   /* Create face for triangle. */
259
   f = MakeNullFace(hull);
260
   f->Edge[0]   = e0;  f->Edge[1]   = e1; f->Edge[2]   = e2;
261
   f->Vertex[0] = v0;  f->Vertex[1] = v1; f->Vertex[2] = v2;
262
 
263
   /* Link edges to face. */
264
   e0->AdjFace[0] = e1->AdjFace[0] = e2->AdjFace[0] = f;
265
 
266
   return f;
267
}
268
 
269
/*---------------------------------------------------------------------
270
Collinear checks to see if the three points given are collinear,
271
by checking to see if each element of the cross product is zero.
272
---------------------------------------------------------------------*/
273
static
274
BOOL Collinear( LPVERTEX a, LPVERTEX b, LPVERTEX c )
275
{
276
   return
277
         ( c->v[Z] - a->v[Z] ) * ( b->v[Y] - a->v[Y] ) -
278
         ( b->v[Z] - a->v[Z] ) * ( c->v[Y] - a->v[Y] ) == 0
279
      && ( b->v[Z] - a->v[Z] ) * ( c->v[X] - a->v[X] ) -
280
         ( b->v[X] - a->v[X] ) * ( c->v[Z] - a->v[Z] ) == 0
281
      && ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) -
282
         ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] ) == 0  ;
283
}
284
 
285
/*---------------------------------------------------------------------
286
VolumeSign returns the sign of the volume of the tetrahedron determined by f
287
and p.  VolumeSign is +1 iff p is on the negative side of f,
288
where the positive side is determined by the rh-rule.  So the volume
289
is positive if the ccw normal to f points outside the tetrahedron.
290
The final fewer-multiplications form is due to Bob Williamson.
291
---------------------------------------------------------------------*/
292
int  VolumeSign( LPFACE f, LPVERTEX p )
293
{
294
   double  vol;  
295
   double  ax, ay, az, bx, by, bz, cx, cy, cz;
296
 
297
   ax = f->Vertex[0]->v[X] - p->v[X];
298
   ay = f->Vertex[0]->v[Y] - p->v[Y];
299
   az = f->Vertex[0]->v[Z] - p->v[Z];
300
   bx = f->Vertex[1]->v[X] - p->v[X];
301
   by = f->Vertex[1]->v[Y] - p->v[Y];
302
   bz = f->Vertex[1]->v[Z] - p->v[Z];
303
   cx = f->Vertex[2]->v[X] - p->v[X];
304
   cy = f->Vertex[2]->v[Y] - p->v[Y];
305
   cz = f->Vertex[2]->v[Z] - p->v[Z];
306
 
307
   vol =   ax * (by*cz - bz*cy)
308
         + ay * (bz*cx - bx*cz)
309
         + az * (bx*cy - by*cx);
310
 
311
 
312
   /* The volume should be an integer. */
313
   if      ( vol >  0.5 )  return  1;
314
   else if ( vol < -0.5 )  return -1;
315
   else                    return  0;
316
}
317
 
318
 
319
 
320
/*---------------------------------------------------------------------
321
CleanEdges runs through the Edge list and cleans up the structure.
322
If there is a NewFace then it will put that face in place of the
323
Visible face and NULL out NewFace. It also deletes so marked edges.
324
---------------------------------------------------------------------*/
325
static
326
void CleanEdges(LPHULL hull)
327
{
328
   LPEDGE  e;   /* Primary index into Edge list. */
329
   LPEDGE  t;   /* Temporary Edge pointer. */
330
 
331
   /* Integrate the NewFace's into the data structure. */
332
   /* Check every Edge. */
333
 
334
   e = hull ->edges;
335
   do {
336
                if ( e->NewFace ) {
337
 
338
                        if ( e->AdjFace[0]->Visible )
339
                                    e->AdjFace[0] = e->NewFace;
340
                        else   
341
                                        e->AdjFace[1] = e->NewFace;
342
 
343
                        e->NewFace = NULL;
344
                }
345
 
346
      e = e->Next;
347
 
348
   } while ( e != hull ->edges );
349
 
350
   /* Delete any edges marked for deletion. */
351
   while ( hull ->edges && hull ->edges->DoDelete ) {
352
 
353
      e = hull ->edges;
354
 
355
      XDELETE( hull ->edges, e );
356
   }
357
 
358
   e = hull ->edges->Next;
359
 
360
   do {
361
                if ( e->DoDelete ) {
362
 
363
                        t = e;
364
                        e = e->Next;
365
                        XDELETE( hull ->edges, t );
366
                }
367
                else e = e->Next;
368
 
369
   } while ( e != hull ->edges );
370
}
371
 
372
/*---------------------------------------------------------------------
373
CleanFaces runs through the face list and deletes any face marked Visible.
374
---------------------------------------------------------------------*/
375
static
376
void CleanFaces(LPHULL hull)
377
{
378
   LPFACE  f;   /* Primary pointer into face list. */
379
   LPFACE  t;   /* Temporary pointer, for deleting. */
380
 
381
 
382
   while ( hull ->faces && hull ->faces->Visible ) {
383
 
384
                f = hull ->faces;
385
                XDELETE( hull ->faces, f );
386
   }
387
 
388
   f = hull ->faces->Next;
389
 
390
   do {
391
      if ( f->Visible ) {
392
 
393
                        t = f;
394
                        f = f->Next;
395
                        XDELETE( hull ->faces, t );
396
      }
397
      else f = f->Next;
398
 
399
   } while ( f != hull ->faces );
400
}
401
 
402
 
403
 
404
/*---------------------------------------------------------------------
405
CleanVertices runs through the Vertex list and deletes the
406
vertices that are marked as processed but are not incident to any
407
undeleted edges.
408
---------------------------------------------------------------------*/
409
static
410
void CleanVertices(LPHULL hull)
411
{
412
   LPEDGE    e;
413
   LPVERTEX  v, t;
414
 
415
   /* Mark all vertices incident to some undeleted Edge as on the hull. */
416
 
417
   e = hull ->edges;
418
   do {
419
                e->EndPts[0]->onhull = e->EndPts[1]->onhull = ONHULL;
420
                e = e->Next;
421
 
422
   } while (e != hull ->edges);
423
 
424
 
425
   /* Delete all vertices that have been processed but
426
      are not on the hull. */
427
 
428
   while ( hull ->vertices && hull->vertices->mark && !hull ->vertices->onhull ) {
429
 
430
                v = hull ->vertices;
431
                XDELETE(hull ->vertices, v );
432
   }
433
 
434
 
435
   v = hull ->vertices->Next;
436
   do {
437
                if (v->mark && !v->onhull ) {    
438
                                t = v;
439
                                v = v->Next;
440
                                XDELETE(hull ->vertices, t )
441
                }
442
                else
443
                                v = v->Next;
444
 
445
   } while ( v != hull ->vertices );
446
 
447
 
448
   /* Reset flags. */
449
 
450
   v = hull ->vertices;
451
   do {
452
                        v->duplicate = NULL;
453
                        v->onhull = !ONHULL;
454
                        v = v->Next;
455
 
456
   } while (v != hull->vertices );
457
 
458
}
459
 
460
 
461
 
462
 
463
/*---------------------------------------------------------------------
464
MakeCcw puts the vertices in the face structure in counterclock wise
465
order.  We want to store the vertices in the same
466
order as in the Visible face.  The third Vertex is always p.
467
 
468
Although no specific ordering of the edges of a face are used
469
by the code, the following condition is maintained for each face f:
470
one of the two endpoints of f->Edge[i] matches f->Vertex[i].
471
But note that this does not imply that f->Edge[i] is between
472
f->Vertex[i] and f->Vertex[(i+1)%3].  (Thanks to Bob Williamson.)
473
---------------------------------------------------------------------*/
474
 
475
static
476
void MakeCcw(LPFACE f, LPEDGE e, LPVERTEX p)
477
{
478
   LPFACE  fv;   /* The Visible face adjacent to e */
479
   int    i;    /* Index of e->endpoint[0] in fv. */
480
   LPEDGE  s;   /* Temporary, for swapping */
481
 
482
   if  (e->AdjFace[0]->Visible)      
483
 
484
        fv = e->AdjFace[0];
485
   else
486
                fv = e->AdjFace[1];
487
 
488
   /* Set Vertex[0] & [1] of f to have the same orientation
489
      as do the corresponding vertices of fv. */
490
 
491
   for ( i=0; fv->Vertex[i] != e->EndPts[0]; ++i )
492
      ;
493
 
494
   /* Orient f the same as fv. */
495
 
496
   if ( fv->Vertex[ (i+1) % 3 ] != e->EndPts[1] ) {
497
 
498
                f->Vertex[0] = e->EndPts[1];  
499
                f->Vertex[1] = e->EndPts[0];    
500
   }
501
   else {                              
502
                f->Vertex[0] = e->EndPts[0];  
503
                f->Vertex[1] = e->EndPts[1];      
504
                SWAP( s, f->Edge[1], f->Edge[2] );
505
   }
506
 
507
   /* This swap is tricky. e is Edge[0]. Edge[1] is based on endpt[0],
508
      Edge[2] on endpt[1].  So if e is oriented "forwards," we
509
      need to move Edge[1] to follow [0], because it precedes. */
510
 
511
   f->Vertex[2] = p;
512
}
513
 
514
/*---------------------------------------------------------------------
515
MakeConeFace makes a new face and two new edges between the
516
Edge and the point that are passed to it. It returns a pointer to
517
the new face.
518
---------------------------------------------------------------------*/
519
 
520
static
521
LPFACE MakeConeFace(LPHULL hull, LPEDGE e, LPVERTEX p)
522
{
523
   LPEDGE  new_edge[2];
524
   LPFACE  new_face;
525
   int    i, j;
526
 
527
   /* Make two new edges (if don't already exist). */
528
 
529
   for ( i=0; i < 2; ++i )
530
      /* If the Edge exists, copy it into new_edge. */
531
      if ( !( new_edge[i] = e->EndPts[i]->duplicate) ) {
532
 
533
                /* Otherwise (duplicate is NULL), MakeNullEdge. */
534
                new_edge[i] = MakeNullEdge(hull);
535
                new_edge[i]->EndPts[0] = e->EndPts[i];
536
                new_edge[i]->EndPts[1] = p;
537
                e->EndPts[i]->duplicate = new_edge[i];
538
      }
539
 
540
        /* Make the new face. */
541
        new_face = MakeNullFace(hull);  
542
        new_face->Edge[0] = e;
543
        new_face->Edge[1] = new_edge[0];
544
        new_face->Edge[2] = new_edge[1];
545
        MakeCcw( new_face, e, p );
546
 
547
   /* Set the adjacent face pointers. */
548
        for ( i=0; i < 2; ++i )
549
                for ( j=0; j < 2; ++j )  
550
                        /* Only one NULL link should be set to new_face. */
551
                        if ( !new_edge[i]->AdjFace[j] ) {
552
                                        new_edge[i]->AdjFace[j] = new_face;
553
                                        break;
554
                        }
555
 
556
   return new_face;
557
}
558
 
559
 
560
/*---------------------------------------------------------------------
561
AddOne is passed a Vertex.  It first determines all faces Visible from
562
that point.  If none are Visible then the point is marked as not
563
onhull.  Next is a loop over edges.  If both faces adjacent to an Edge
564
are Visible, then the Edge is marked for deletion.  If just one of the
565
adjacent faces is Visible then a new face is constructed.
566
---------------------------------------------------------------------*/
567
static
568
BOOL AddOne(LPHULL hull, LPVERTEX p)
569
{
570
   LPFACE  f;
571
   LPEDGE  e, temp;
572
   int    vol;
573
   BOOL   vis = FALSE;
574
 
575
 
576
   /* Mark faces Visible from p. */
577
   f = hull -> faces;
578
 
579
   do {
580
 
581
           vol = VolumeSign(f, p);
582
 
583
      if ( vol < 0 ) {
584
                        f->Visible = VISIBLE;  
585
                        vis = TRUE;                      
586
      }
587
 
588
          f = f->Next;
589
 
590
   } while ( f != hull ->faces );
591
 
592
   /* If no faces are Visible from p, then p is inside the hull. */
593
 
594
   if ( !vis ) {
595
 
596
                        p->onhull = !ONHULL;  
597
                        return FALSE;
598
   }
599
 
600
   /* Mark edges in interior of Visible region for deletion.
601
      Erect a NewFace based on each border Edge. */
602
 
603
   e = hull ->edges;
604
 
605
   do {
606
 
607
      temp = e->Next;
608
 
609
      if ( e->AdjFace[0]->Visible && e->AdjFace[1]->Visible )
610
                                /* e interior: mark for deletion. */
611
                                e->DoDelete = REMOVED;
612
 
613
      else
614
                  if ( e->AdjFace[0]->Visible || e->AdjFace[1]->Visible )
615
                                /* e border: make a new face. */
616
                                e->NewFace = MakeConeFace(hull, e, p );
617
 
618
                e = temp;
619
 
620
   } while ( e != hull ->edges );
621
 
622
   return TRUE;
623
}
624
 
625
 
626
/*---------------------------------------------------------------------
627
 DoubleTriangle builds the initial double triangle.  It first finds 3
628
 noncollinear points and makes two faces out of them, in opposite order.
629
 It then finds a fourth point that is not coplanar with that face.  The  
630
 vertices are stored in the face structure in counterclockwise order so
631
 that the volume between the face and the point is negative. Lastly, the
632
 3 newfaces to the fourth point are constructed and the data structures
633
 are cleaned up.
634
---------------------------------------------------------------------*/
635
 
636
static
637
BOOL DoubleTriangle(LPHULL hull)
638
{
639
   LPVERTEX  v0, v1, v2, v3;
640
   LPFACE    f0, f1 = NULL;
641
   int      vol;
642
 
643
   /* Find 3 noncollinear points. */
644
   v0 = hull ->vertices;
645
   while ( Collinear( v0, v0->Next, v0->Next->Next ) )
646
      if ( ( v0 = v0->Next ) == hull->vertices )
647
                                return FALSE; // All points are Collinear!
648
 
649
   v1 = v0->Next;
650
   v2 = v1->Next;
651
 
652
   /* Mark the vertices as processed. */
653
   v0->mark = PROCESSED;
654
   v1->mark = PROCESSED;
655
   v2->mark = PROCESSED;
656
 
657
   /* Create the two "twin" faces. */
658
   f0 = MakeFace(hull, v0, v1, v2, f1 );
659
   f1 = MakeFace(hull, v2, v1, v0, f0 );
660
 
661
   /* Link adjacent face fields. */
662
   f0->Edge[0]->AdjFace[1] = f1;
663
   f0->Edge[1]->AdjFace[1] = f1;
664
   f0->Edge[2]->AdjFace[1] = f1;
665
   f1->Edge[0]->AdjFace[1] = f0;
666
   f1->Edge[1]->AdjFace[1] = f0;
667
   f1->Edge[2]->AdjFace[1] = f0;
668
 
669
   /* Find a fourth, noncoplanar point to form tetrahedron. */
670
   v3 = v2->Next;
671
   vol = VolumeSign( f0, v3 );
672
 
673
   while ( !vol )   {
674
 
675
      if ( ( v3 = v3->Next ) == v0 )
676
                        return FALSE; // All points are coplanar!
677
 
678
      vol = VolumeSign( f0, v3 );
679
   }
680
 
681
   /* Insure that v3 will be the first added. */
682
   hull ->vertices = v3;
683
   return TRUE;
684
}
685
 
686
 
687
 
688
/*---------------------------------------------------------------------
689
ConstructHull adds the vertices to the hull one at a time.  The hull
690
vertices are those in the list marked as onhull.
691
---------------------------------------------------------------------*/
692
static
693
void ConstructHull(LPHULL hull)
694
{
695
 LPVERTEX  v, vnext;
696
 BOOL     changed;    /* T if addition changes hull; not used. */
697
 
698
 v = hull->vertices;
699
 
700
 do {
701
                vnext = v->Next;
702
 
703
        changed = FALSE;
704
 
705
                if (!v->mark ) {
706
 
707
                v->mark = PROCESSED;
708
                changed = AddOne(hull,  v );
709
 
710
                                CleanEdges(hull);
711
                                CleanFaces(hull);
712
                                CleanVertices(hull);
713
    }
714
 
715
    v = vnext;
716
 
717
 } while (v != hull->vertices );
718
 
719
}
720
 
721
 
722
 
723
/*-------------------------------------------------------------------*/
724
 
725
 
726
static
727
void AddVec( VEC3I q, VEC3I ray )
728
{
729
  int i;
730
 
731
  for( i = 0; i < DIM; i++ )
732
    ray[i] = q[i] + ray[i];
733
}
734
 
735
/*---------------------------------------------------------------------
736
a - b ==> c.
737
---------------------------------------------------------------------*/
738
static
739
void  SubVec( VEC3I a, VEC3I b, VEC3I c )
740
{
741
   int i;
742
 
743
   for( i = 0; i < DIM; i++ )
744
      c[i] = a[i] - b[i];
745
}
746
 
747
 
748
/*---------------------------------------------------------------------
749
Returns the dot product of the two input vectors.
750
---------------------------------------------------------------------*/
751
static
752
double  Dot( VEC3I a, LPVEC3 b )
753
{
754
    int i;
755
    double sum = 0.0;
756
 
757
    for( i = 0; i < DIM; i++ )
758
       sum += a[i] * b->n[i];
759
 
760
    return  sum;
761
}
762
 
763
/*---------------------------------------------------------------------
764
Compute the cross product of (b-a)x(c-a) and place into N.
765
---------------------------------------------------------------------*/
766
static
767
void    NormalVec( VEC3I a, VEC3I b, VEC3I c, LPVEC3 N )
768
{
769
    N->n[X] = ( c[Z] - a[Z] ) * ( b[Y] - a[Y] ) -
770
           ( b[Z] - a[Z] ) * ( c[Y] - a[Y] );
771
    N->n[Y] = ( b[Z] - a[Z] ) * ( c[X] - a[X] ) -
772
           ( b[X] - a[X] ) * ( c[Z] - a[Z] );
773
    N->n[Z] = ( b[X] - a[X] ) * ( c[Y] - a[Y] ) -
774
           ( b[Y] - a[Y] ) * ( c[X] - a[X] );
775
}
776
 
777
 
778
 
779
 
780
static
781
int InBox( VEC3I q, VEC3I bmin, VEC3I bmax )
782
{
783
 
784
  if( ( bmin[X] <= q[X] ) && ( q[X] <= bmax[X] ) &&
785
      ( bmin[Y] <= q[Y] ) && ( q[Y] <= bmax[Y] ) &&
786
      ( bmin[Z] <= q[Z] ) && ( q[Z] <= bmax[Z] ) )
787
    return TRUE;
788
 
789
  return FALSE;
790
}
791
 
792
 
793
 
794
/*
795
  This function returns a char:
796
    '0': the segment [ab] does not intersect (completely misses) the
797
         bounding box surrounding the n-th triangle T.  It lies
798
         strictly to one side of one of the six supporting planes.
799
    '?': status unknown: the segment may or may not intersect T.
800
*/
801
 
802
static
803
char BoxTest(LPHULL hull, int n, VEC3I a, VEC3I b)
804
{
805
   int i; /* Coordinate index */
806
   int w;
807
 
808
   for ( i=0; i < DIM; i++ ) {
809
 
810
       w = hull ->Box[ n ][0][i]; /* min: lower left */
811
 
812
           if ( (a[i] < w) && (b[i] < w) ) return '0';
813
 
814
       w = hull ->Box[ n ][1][i]; /* max: upper right */
815
 
816
       if ( (a[i] > w) && (b[i] > w) ) return '0';
817
   }
818
 
819
   return '?';
820
}
821
 
822
 
823
 
824
/* Return a random ray endpoint */
825
 
826
static
827
void RandomRay( VEC3I ray, int radius )
828
{
829
  double x, y, z, w, t;
830
 
831
  /* Generate a random point on a sphere of radius 1. */
832
  /* the sphere is sliced at z, and a random point at angle t
833
     generated on the circle of intersection. */
834
 
835
  z = 2.0 * (double) rand() / RAND_MAX - 1.0;
836
  t = 2.0 * M_PI * (double) rand() / RAND_MAX;
837
  w = sqrt( 1 - z*z );
838
  x = w * cos( t );
839
  y = w * sin( t );
840
 
841
  ray[X] = (int) ( radius * x );
842
  ray[Y] = (int) ( radius * y );
843
  ray[Z] = (int) ( radius * z );
844
}
845
 
846
 
847
 
848
static
849
int ComputeBox(LPHULL hull, int F, VEC3I bmin, VEC3I bmax )
850
{
851
  int i, j;
852
  double radius;
853
 
854
  for( i = 0; i < F; i++ )
855
    for( j = 0; j < DIM; j++ ) {
856
 
857
      if( hull ->Vertices[i][j] < bmin[j] )
858
                        bmin[j] = hull ->Vertices[i][j];
859
 
860
      if( hull ->Vertices[i][j] > bmax[j] )
861
                                bmax[j] = hull ->Vertices[i][j];
862
    }
863
 
864
  radius = sqrt( pow( (double)(bmax[X] - bmin[X]), 2.0 ) +
865
                 pow( (double)(bmax[Y] - bmin[Y]), 2.0 ) +
866
                 pow( (double)(bmax[Z] - bmin[Z]), 2.0 ) );
867
 
868
  return (int)( radius +1 ) + 1;
869
}
870
 
871
 
872
/*---------------------------------------------------------------------
873
Computes N & D and returns index m of largest component.
874
---------------------------------------------------------------------*/
875
static
876
int     PlaneCoeff(LPHULL hull,  VEC3I T, LPVEC3 N, double *D )
877
{
878
    int i;
879
    double t;              /* Temp storage */
880
    double biggest = 0.0;  /* Largest component of normal vector. */
881
    int m = 0;             /* Index of largest component. */
882
 
883
    NormalVec(hull ->Vertices[T[0]], hull ->Vertices[T[1]], hull ->Vertices[T[2]], N );
884
    *D = Dot( hull ->Vertices[T[0]], N );
885
 
886
    /* Find the largest component of N. */
887
    for ( i = 0; i < DIM; i++ ) {
888
      t = fabs( N->n[i] );
889
      if ( t > biggest ) {
890
        biggest = t;
891
        m = i;
892
      }
893
    }
894
    return m;
895
}
896
 
897
/*---------------------------------------------------------------------
898
    'p': The segment lies wholly within the plane.
899
    'q': The q endpoint is on the plane (but not 'p').
900
    'r': The r endpoint is on the plane (but not 'p').
901
    '0': The segment lies strictly to one side or the other of the plane.
902
    '1': The segement intersects the plane, and 'p' does not hold.
903
---------------------------------------------------------------------*/
904
static
905
char SegPlaneInt(LPHULL hull, VEC3I T, VEC3I q, VEC3I r, LPVEC3 p, int *m)
906
{
907
    VEC3 N; double D;
908
    VEC3I rq;
909
    double num, denom, t;
910
    int i;
911
 
912
    *m = PlaneCoeff(hull, T, &N, &D );
913
    num = D - Dot( q, &N );
914
    SubVec( r, q, rq );
915
    denom = Dot( rq, &N );
916
 
917
    if ( denom == 0.0 ) {  /* Segment is parallel to plane. */
918
       if ( num == 0.0 )   /* q is on plane. */
919
           return 'p';
920
       else
921
           return '0';
922
    }
923
    else
924
       t = num / denom;
925
 
926
    for( i = 0; i < DIM; i++ )
927
       p->n[i] = q[i] + t * ( r[i] - q[i] );
928
 
929
    if ( (0.0 < t) && (t < 1.0) )
930
         return '1';
931
    else if ( num == 0.0 )   /* t == 0 */
932
         return 'q';
933
    else if ( num == denom ) /* t == 1 */
934
         return 'r';
935
    else return '0';
936
}
937
 
938
 
939
 
940
static
941
int  AreaSign( VEC3I a, VEC3I b, VEC3I c )  
942
{
943
    double area2;
944
 
945
    area2 = ( b[0] - a[0] ) * (double)( c[1] - a[1] ) -
946
            ( c[0] - a[0] ) * (double)( b[1] - a[1] );
947
 
948
    /* The area should be an integer. */
949
    if      ( area2 >  0.5 ) return  1;
950
    else if ( area2 < -0.5 ) return -1;
951
    else                     return  0;
952
}    
953
 
954
 
955
static
956
char    InTri2D( VEC3I Tp[3], VEC3I pp )
957
{
958
   int area0, area1, area2;
959
 
960
   /* compute three AreaSign() values for pp w.r.t. each Edge of the face in 2D */
961
   area0 = AreaSign( pp, Tp[0], Tp[1] );
962
   area1 = AreaSign( pp, Tp[1], Tp[2] );
963
   area2 = AreaSign( pp, Tp[2], Tp[0] );
964
 
965
   if ( (( area0 == 0 ) && ( area1 > 0 ) && ( area2 > 0 )) ||
966
        (( area1 == 0 ) && ( area0 > 0 ) && ( area2 > 0 )) ||
967
        (( area2 == 0 ) && ( area0 > 0 ) && ( area1 > 0 )) )
968
     return 'E';
969
 
970
   if ( (( area0 == 0 ) && ( area1 < 0 ) && ( area2 < 0 )) ||
971
        (( area1 == 0 ) && ( area0 < 0 ) && ( area2 < 0 )) ||
972
        (( area2 == 0 ) && ( area0 < 0 ) && ( area1 < 0 )))
973
     return 'E';                
974
 
975
   if ( (( area0 >  0 ) && ( area1 > 0 ) && ( area2 > 0 )) ||
976
        (( area0 <  0 ) && ( area1 < 0 ) && ( area2 < 0 )))
977
     return 'F';
978
 
979
   if ( ( area0 == 0 ) && ( area1 == 0 ) && ( area2 == 0 ) )
980
     return '?'; // Error in InTriD
981
 
982
   if ( (( area0 == 0 ) && ( area1 == 0 )) ||
983
        (( area0 == 0 ) && ( area2 == 0 )) ||
984
        (( area1 == 0 ) && ( area2 == 0 )) )
985
     return 'V';
986
 
987
   else  
988
     return '0';  
989
}
990
 
991
/* Assumption: p lies in the plane containing T.
992
    Returns a char:
993
     'V': the query point p coincides with a Vertex of triangle T.
994
     'E': the query point p is in the relative interior of an Edge of triangle T.
995
     'F': the query point p is in the relative interior of a Face of triangle T.
996
     '0': the query point p does not intersect (misses) triangle T.
997
*/
998
 
999
static
1000
char    InTri3D(LPHULL hull,  VEC3I T, int m, VEC3I p )
1001
{
1002
   int i;           /* Index for X,Y,Z           */
1003
   int j;           /* Index for X,Y             */
1004
   int k;           /* Index for triangle Vertex */
1005
   VEC3I pp;      /* projected p */
1006
   VEC3I Tp[3];   /* projected T: three new vertices */
1007
 
1008
   /* Project out coordinate m in both p and the triangular face */
1009
   j = 0;
1010
   for ( i = 0; i < DIM; i++ ) {
1011
     if ( i != m ) {    /* skip largest coordinate */
1012
       pp[j] = p[i];
1013
       for ( k = 0; k < 3; k++ )
1014
                        Tp[k][j] = hull->Vertices[T[k]][i];
1015
       j++;
1016
     }
1017
   }
1018
   return( InTri2D( Tp, pp ) );
1019
}
1020
 
1021
 
1022
 
1023
static        
1024
int     VolumeSign2( VEC3I a, VEC3I b, VEC3I c, VEC3I d )
1025
{
1026
   double vol;
1027
   double ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz;
1028
   double bxdx, bydy, bzdz, cxdx, cydy, czdz;
1029
 
1030
   ax = a[X];
1031
   ay = a[Y];
1032
   az = a[Z];
1033
   bx = b[X];
1034
   by = b[Y];
1035
   bz = b[Z];
1036
   cx = c[X];
1037
   cy = c[Y];
1038
   cz = c[Z];
1039
   dx = d[X];
1040
   dy = d[Y];
1041
   dz = d[Z];
1042
 
1043
   bxdx=bx-dx;
1044
   bydy=by-dy;
1045
   bzdz=bz-dz;
1046
   cxdx=cx-dx;
1047
   cydy=cy-dy;
1048
   czdz=cz-dz;
1049
   vol =   (az-dz) * (bxdx*cydy - bydy*cxdx)
1050
         + (ay-dy) * (bzdz*cxdx - bxdx*czdz)
1051
         + (ax-dx) * (bydy*czdz - bzdz*cydy);
1052
 
1053
 
1054
   /* The volume should be an integer. */
1055
   if      ( vol > 0.5 )   return  1;
1056
   else if ( vol < -0.5 )  return -1;
1057
   else                    return  0;
1058
}
1059
 
1060
 
1061
 
1062
 
1063
/*---------------------------------------------------------------------
1064
The signed volumes of three tetrahedra are computed, determined
1065
by the segment qr, and each Edge of the triangle.  
1066
Returns a char:
1067
   'v': the open segment includes a Vertex of T.
1068
   'e': the open segment includes a point in the relative interior of an Edge
1069
   of T.
1070
   'f': the open segment includes a point in the relative interior of a face
1071
   of T.
1072
   '0': the open segment does not intersect triangle T.
1073
---------------------------------------------------------------------*/
1074
 
1075
static
1076
char SegTriCross(LPHULL hull,  VEC3I T, VEC3I q, VEC3I r )
1077
{
1078
   int vol0, vol1, vol2;
1079
 
1080
   vol0 = VolumeSign2( q, hull->Vertices[ T[0] ], hull->Vertices[ T[1] ], r );
1081
   vol1 = VolumeSign2( q, hull->Vertices[ T[1] ], hull->Vertices[ T[2] ], r );
1082
   vol2 = VolumeSign2( q, hull->Vertices[ T[2] ], hull->Vertices[ T[0] ], r );
1083
 
1084
 
1085
   /* Same sign: segment intersects interior of triangle. */
1086
   if ( ( ( vol0 > 0 ) && ( vol1 > 0 ) && ( vol2 > 0 ) ) ||
1087
        ( ( vol0 < 0 ) && ( vol1 < 0 ) && ( vol2 < 0 ) ) )
1088
      return 'f';
1089
 
1090
   /* Opposite sign: no intersection between segment and triangle */
1091
   if ( ( ( vol0 > 0 ) || ( vol1 > 0 ) || ( vol2 > 0 ) ) &&
1092
        ( ( vol0 < 0 ) || ( vol1 < 0 ) || ( vol2 < 0 ) ) )
1093
      return '0';
1094
 
1095
   else if ( ( vol0 == 0 ) && ( vol1 == 0 ) && ( vol2 == 0 ) )
1096
     return '?'; // Error 1 in SegTriCross
1097
 
1098
   /* Two zeros: segment intersects Vertex. */
1099
   else if ( ( ( vol0 == 0 ) && ( vol1 == 0 ) ) ||
1100
             ( ( vol0 == 0 ) && ( vol2 == 0 ) ) ||
1101
             ( ( vol1 == 0 ) && ( vol2 == 0 ) ) )
1102
      return 'v';
1103
 
1104
   /* One zero: segment intersects Edge. */
1105
   else if ( ( vol0 == 0 ) || ( vol1 == 0 ) || ( vol2 == 0 ) )
1106
      return 'e';
1107
 
1108
   else
1109
     return '?'; // Error 2 in SegTriCross
1110
}
1111
 
1112
 
1113
 
1114
static
1115
char    SegTriInt(LPHULL hull, VEC3I T, VEC3I q, VEC3I r, LPVEC3 p )
1116
{
1117
    int code;
1118
    int m = -1;
1119
 
1120
    code = SegPlaneInt(hull, T, q, r, p, &m );
1121
 
1122
    if ( code == '0')        return '0';
1123
    else if ( code == 'q')   return InTri3D(hull, T, m, q );
1124
    else if ( code == 'r')   return InTri3D(hull,  T, m, r );
1125
    else if ( code == 'p')   return 'p';
1126
    else if ( code == '1' )  return SegTriCross(hull, T, q, r );
1127
    else
1128
       return code;             // Error
1129
}
1130
 
1131
 
1132
 
1133
 
1134
/*
1135
  This function returns a char:
1136
    'i': the query point a is strictly interior to polyhedron P.
1137
    'o': the query point a is strictly exterior to( or outside of) polyhedron P.
1138
*/
1139
char InPolyhedron(LPHULL hull, VEC3I q)
1140
{
1141
   int F = hull->nfaces;
1142
   VEC3I Ray;  /* Ray endpoint. */
1143
   VEC3 p;  /* Intersection point; not used. */
1144
   int f, k = 0, crossings = 0;
1145
   char code = '?';
1146
 
1147
 
1148
   /* If query point is outside bounding box, finished. */
1149
   if ( !InBox( q, hull->bmin, hull->bmax ) )
1150
      return 'o';
1151
 
1152
   LOOP:
1153
   while( k++ < F ) {
1154
 
1155
      crossings = 0;
1156
 
1157
      RandomRay(Ray, hull->radius );
1158
 
1159
      AddVec( q, Ray );
1160
 
1161
 
1162
      for ( f = 0; f < F; f++ ) {  /* Begin check each face */
1163
 
1164
         if ( BoxTest(hull,  f, q, Ray ) == '0' ) {
1165
              code = '0';
1166
 
1167
         }
1168
         else code = SegTriInt(hull, hull->Faces[f], q, Ray, &p );
1169
 
1170
 
1171
         /* If ray is degenerate, then goto outer while to generate another. */
1172
         if ( code == 'p' || code == 'v' || code == 'e' ) {
1173
 
1174
            goto LOOP;
1175
         }
1176
 
1177
         /* If ray hits face at interior point, increment crossings. */
1178
         else if ( code == 'f' ) {
1179
            crossings++;
1180
 
1181
         }
1182
 
1183
         /* If query endpoint q sits on a V/E/F, return inside. */
1184
         else if ( code == 'V' || code == 'E' || code == 'F' )
1185
                        return code; // 'i'; // MM2
1186
 
1187
         /* If ray misses triangle, do nothing. */
1188
         else if ( code == '0' )
1189
            ;
1190
 
1191
         else
1192
            return '?'; // Error
1193
 
1194
      }
1195
      /* No degeneracies encountered: ray is generic, so finished. */
1196
      break;
1197
 
1198
   } /* End while loop */
1199
 
1200
 
1201
   /* q strictly interior to polyhedron iff an odd number of crossings. */
1202
   if( ( crossings % 2 ) == 1 )
1203
      return   'i';
1204
 
1205
   else return 'o';
1206
}
1207
 
1208
 
1209
/// ----------------------------------------------------------------------------------
1210
 
1211
 
1212
 
1213
 
1214
static
1215
void StoreResults(LPHULL hull)
1216
{
1217
 
1218
   int   i, w;
1219
   LPVERTEX  v;
1220
   LPFACE    f;
1221
   int  V = 0, F = 0;
1222
   int j, k;
1223
 
1224
   // Vertices
1225
 
1226
   v = hull ->vertices;
1227
   V = 0;
1228
   do {
1229
 
1230
      v -> vnum = V;
1231
          hull ->Vertices[V][X] = v -> v[X];
1232
          hull ->Vertices[V][Y] = v -> v[Y];
1233
          hull ->Vertices[V][Z] = v -> v[Z];
1234
 
1235
      v = v->Next;
1236
      V++;
1237
 
1238
   } while ( v != hull ->vertices );
1239
 
1240
   hull ->nvertex = V;
1241
 
1242
   // Faces
1243
   f = hull ->faces;
1244
   F = 0;
1245
   do {
1246
 
1247
          hull ->Faces[F][0] = f->Vertex[0]->vnum;
1248
          hull ->Faces[F][1] = f->Vertex[1]->vnum;
1249
          hull ->Faces[F][2] = f->Vertex[2]->vnum;
1250
 
1251
          for ( j=0; j < 3; j++ ) {
1252
 
1253
       hull ->Box[F][0][j] = hull ->Vertices[ hull ->Faces[F][0] ][j];
1254
       hull ->Box[F][1][j] = hull ->Vertices[ hull ->Faces[F][0] ][j];
1255
          }
1256
 
1257
          /* Check k=1,2 vertices of face. */
1258
          for ( k=1; k < 3; k++ )
1259
                for ( j=0; j < 3; j++ ) {
1260
 
1261
                        w = hull ->Vertices[ hull ->Faces[F][k] ][j];
1262
                        if ( w < hull ->Box[F][0][j] ) hull ->Box[F][0][j] = w;
1263
                        if ( w > hull ->Box[F][1][j] ) hull ->Box[F][1][j] = w;
1264
                }
1265
 
1266
 
1267
      f = f->Next; F++;
1268
 
1269
   } while ( f != hull ->faces );
1270
 
1271
 
1272
  hull ->nfaces = F;
1273
 
1274
 
1275
    /* Initialize the bounding box */
1276
  for ( i = 0; i < DIM; i++ )
1277
    hull ->bmin[i] = hull ->bmax[i] = hull ->Vertices[0][i];
1278
 
1279
  hull ->radius = ComputeBox(hull, V, hull ->bmin, hull ->bmax );
1280
 
1281
 
1282
}
1283
 
1284
 
1285
LCMSHANDLE cmsxHullInit(void)
1286
{
1287
        LPHULL hull = (LPHULL) malloc(sizeof(HULL));
1288
 
1289
        ZeroMemory(hull, sizeof(HULL));
1290
 
1291
        hull->vnumCounter  = 0;
1292
        hull->vertices   = NULL;
1293
        hull->edges     = NULL;
1294
        hull->faces     = NULL;
1295
        hull->nfaces    = 0;
1296
        hull->nvertex   = 0;
1297
 
1298
        return (LCMSHANDLE) (LPSTR) hull;
1299
}
1300
 
1301
 
1302
void cmsxHullDone(LCMSHANDLE hHull)
1303
{
1304
        LPHULL hull = (LPHULL) (LPSTR) hHull;
1305
 
1306
        if (hull)
1307
                free((LPVOID) hull);
1308
}
1309
 
1310
 
1311
BOOL cmsxHullAddPoint(LCMSHANDLE hHull, int x, int y, int z)
1312
{
1313
         LPVERTEX  v;
1314
         LPHULL hull = (LPHULL) (LPSTR) hHull;
1315
 
1316
 
1317
      v = MakeNullVertex(hull);
1318
      v->v[X] = x;
1319
      v->v[Y] = y;
1320
      v->v[Z] = z;
1321
      v->vnum = hull->vnumCounter++;
1322
 
1323
          return TRUE;
1324
}
1325
 
1326
BOOL cmsxHullComputeHull(LCMSHANDLE hHull)
1327
{
1328
 
1329
  LPHULL hull = (LPHULL) (LPSTR) hHull;
1330
 
1331
  if (!DoubleTriangle(hull)) return FALSE;
1332
 
1333
  ConstructHull(hull);
1334
  StoreResults(hull);
1335
 
1336
  return TRUE;
1337
}
1338
 
1339
 
1340
char cmsxHullCheckpoint(LCMSHANDLE hHull, int x, int y, int z)
1341
{
1342
         VEC3I q;
1343
         LPHULL hull = (LPHULL) (LPSTR) hHull;
1344
 
1345
     q[X] = x; q[Y] = y; q[Z] = z;
1346
 
1347
     return InPolyhedron(hull, q ) ;
1348
}
1349
 
1350
 
1351
BOOL cmsxHullDumpVRML(LCMSHANDLE hHull, const char* fname)
1352
{
1353
        FILE*           fp;
1354
        int             i;
1355
        LPHULL hull = (LPHULL) (LPSTR) hHull;
1356
 
1357
        fp = fopen (fname, "wt");
1358
        if (fp == NULL)
1359
                return FALSE;
1360
 
1361
        fprintf (fp, "#VRML V2.0 utf8\n");
1362
 
1363
        /* set the viewing orientation and distance */
1364
        fprintf (fp, "DEF CamTest Group {\n");
1365
        fprintf (fp, "\tchildren [\n");
1366
        fprintf (fp, "\t\tDEF Cameras Group {\n");
1367
        fprintf (fp, "\t\t\tchildren [\n");
1368
        fprintf (fp, "\t\t\t\tDEF DefaultView Viewpoint {\n");
1369
        fprintf (fp, "\t\t\t\t\tposition 0 0 340\n");
1370
        fprintf (fp, "\t\t\t\t\torientation 0 0 1 0\n");
1371
        fprintf (fp, "\t\t\t\t\tdescription \"default view\"\n");
1372
        fprintf (fp, "\t\t\t\t}\n");
1373
        fprintf (fp, "\t\t\t]\n");
1374
        fprintf (fp, "\t\t},\n");
1375
        fprintf (fp, "\t]\n");
1376
        fprintf (fp, "}\n");
1377
 
1378
        /* Output the background stuff */
1379
        fprintf (fp, "Background {\n");
1380
        fprintf (fp, "\tskyColor [\n");
1381
        fprintf (fp, "\t\t.5 .5 .5\n");
1382
        fprintf (fp, "\t]\n");
1383
        fprintf (fp, "}\n");
1384
 
1385
        /* Output the shape stuff */
1386
        fprintf (fp, "Transform {\n");
1387
        fprintf (fp, "\tscale 8 8 8\n");
1388
        fprintf (fp, "\tchildren [\n");
1389
 
1390
        /* Draw the axes as a shape: */
1391
        fprintf (fp, "\t\tShape {\n");
1392
        fprintf (fp, "\t\t\tappearance Appearance {\n");
1393
        fprintf (fp, "\t\t\t\tmaterial Material {\n");
1394
        fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
1395
        fprintf (fp, "\t\t\t\t\temissiveColor 1.0 1.0 1.0\n");
1396
        fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
1397
        fprintf (fp, "\t\t\t\t}\n");
1398
        fprintf (fp, "\t\t\t}\n");
1399
        fprintf (fp, "\t\t\tgeometry IndexedLineSet {\n");
1400
        fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
1401
        fprintf (fp, "\t\t\t\t\tpoint [\n");
1402
        fprintf (fp, "\t\t\t\t\t0.0 0.0 0.0,\n");
1403
        fprintf (fp, "\t\t\t\t\t%f 0.0 0.0,\n",  255.0);
1404
        fprintf (fp, "\t\t\t\t\t0.0 %f 0.0,\n",  255.0);
1405
        fprintf (fp, "\t\t\t\t\t0.0 0.0 %f]\n",  255.0);
1406
        fprintf (fp, "\t\t\t\t}\n");
1407
        fprintf (fp, "\t\t\t\tcoordIndex [\n");
1408
        fprintf (fp, "\t\t\t\t\t0, 1, -1\n");
1409
        fprintf (fp, "\t\t\t\t\t0, 2, -1\n");
1410
        fprintf (fp, "\t\t\t\t\t0, 3, -1]\n");
1411
        fprintf (fp, "\t\t\t}\n");
1412
        fprintf (fp, "\t\t}\n");
1413
 
1414
 
1415
        /* Draw the triangles as a shape: */
1416
        fprintf (fp, "\t\tShape {\n");
1417
        fprintf (fp, "\t\t\tappearance Appearance {\n");
1418
        fprintf (fp, "\t\t\t\tmaterial Material {\n");
1419
        fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
1420
        fprintf (fp, "\t\t\t\t\temissiveColor 0 0 0\n");
1421
        fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
1422
        fprintf (fp, "\t\t\t\t}\n");
1423
        fprintf (fp, "\t\t\t}\n");
1424
        fprintf (fp, "\t\t\tgeometry IndexedFaceSet {\n");
1425
        fprintf (fp, "\t\t\t\tsolid FALSE\n");
1426
 
1427
        /* fill in the points here */
1428
        fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
1429
        fprintf (fp, "\t\t\t\t\tpoint [\n");
1430
 
1431
        for (i = 0; i < hull->nvertex; ++i)
1432
        {
1433
                fprintf (fp, "\t\t\t\t\t%g %g %g%c\n",
1434
                        (double) hull->Vertices[i][X], (double) hull->Vertices[i][Y], (double) hull->Vertices[i][Z],
1435
                        i == hull->nvertex-1? ']': ',');
1436
        }
1437
        fprintf (fp, "\t\t\t\t}\n");
1438
 
1439
        /* fill in the Vertex indices (followed by -1) */
1440
 
1441
 
1442
        fprintf (fp, "\t\t\t\tcoordIndex [\n");
1443
        for (i = 0; i < hull->nfaces; ++i)
1444
        {
1445
                fprintf (fp, "\t\t\t\t\t%d, %d, %d, -1\n",
1446
                        hull->Faces[i][0], hull->Faces[i][1], hull->Faces[i][2]);
1447
 
1448
        }
1449
        fprintf (fp, "]\n");
1450
 
1451
 
1452
        /* fill in the face colors */
1453
        fprintf (fp, "\t\t\t\tcolor Color {\n");
1454
        fprintf (fp, "\t\t\t\t\tcolor [\n");
1455
        for (i = 0; i < hull->nfaces; ++i)
1456
        {
1457
                int vx, vy, vz;
1458
                double r, g, b;
1459
 
1460
                vx = hull->Faces[i][0]; vy = hull->Faces[i][1]; vz = hull->Faces[i][2];
1461
                r = (double) (hull->Vertices[vx][X] + hull->Vertices[vy][X] + hull->Vertices[vz][X]) / (3* 255);
1462
                g = (double) (hull->Vertices[vx][Y] + hull->Vertices[vy][Y] + hull->Vertices[vz][Y]) / (3* 255);
1463
                b = (double) (hull->Vertices[vx][Z] + hull->Vertices[vy][Z] + hull->Vertices[vz][Z]) / (3* 255);
1464
 
1465
                fprintf (fp, "\t\t\t\t\t%g %g %g%c\n", r, g, b,                
1466
                        i == hull->nfaces-1? ']': ',');
1467
        }
1468
        fprintf (fp, "\t\t\t}\n");
1469
 
1470
        fprintf (fp, "\t\t\tcolorPerVertex FALSE\n");
1471
 
1472
        fprintf (fp, "\t\t\t}\n");
1473
        fprintf (fp, "\t\t}\n");
1474
        fprintf (fp, "\t]\n");
1475
        fprintf (fp, "}\n");
1476
 
1477
        fclose (fp);
1478
 
1479
        return TRUE;
1480
}