home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch5_5 / voronoi.cxx < prev    next >
C/C++ Source or Header  |  1995-03-04  |  9KB  |  181 lines
  1. #include <iostream.h>
  2. #include <math.h>
  3.  
  4. #include "global.h"
  5. #include "voronoi.h"                            // D-DIM. VORONOI TEMPLATE
  6.  
  7. struct cell : public VECTOR<3> {                // 3-DIMENSIONAL VORONOI-CELL
  8.         vector p;                               // POSITION OF PARTICLE
  9.         list<cell*> ln;                         // CONTIGUOUS CELLS
  10.         list<object*> lh;                       // HOST OBJECT(S)
  11.         list<object*> lo;                       // LIST OF INTERSECTING OBJECTS
  12.         long t;                                 // LOCAL TRAVERSE CODE (WORK)
  13.         cell *b;                                // BACK (voronoi::disperse)
  14.         cell() {t=0L;}
  15.         cell(vector& v, object *o=(object*)0) {
  16.                 double x[3]; x[0]=v[0]; x[1]=v[1]; x[2]=v[2];
  17.                 *((VECTOR<3>*)this)=VECTOR<3>(x); p=v; {if(o)lh+=o;} t=0L;
  18.         }
  19.         int operator&(vector& p) {              // CELL CONTAINS POINT p
  20.                 for(cell*n=ln.first();n;n=ln.next())
  21.                         if(!(halfspace(this->p,n->p)&p)) return 0;
  22.                 return 1;
  23.         }
  24.         int operator&(object& o) {              // CELL INTERSECTS OBJECT o
  25.                 for(cell*n=ln.first();n;n=ln.next())
  26.                         if(!(o&halfspace(this->p,n->p))) return 0;
  27.                 return 1;
  28.         }
  29. };
  30.  
  31. //      GLOBAL METHODS (NOT REALLY AN OBJECT ORIENTED APPROACH!)
  32.  
  33. static void initcells(VERTEX<3>*v) {        // VORONOI-CELLS FROM VORONOI<3>
  34.         for(register int i=0; i<3; i++) {
  35.                 cell *ci=(cell*)(v->p[i]);
  36.                 for(register int j=i+1; j<4; j++) {
  37.                         cell *cj=(cell*)(v->p[j]);
  38.                         if(!ci->ln[cj]) ci->ln+=cj;
  39.                         if(!cj->ln[ci]) cj->ln+=ci;
  40.                 }
  41.         }
  42. }
  43.  
  44. static int lexcmp(const void *v1, const void *v2) {
  45.         cell *c1=*(cell**)v1; cell *c2=*(cell**)v2;
  46.         const double EPS=1e-6;
  47.         for(register int i=0; i<3; i++) {
  48.                 double d=(*c1)[i]-(*c2)[i];
  49.                 if(d<-EPS) return -1;
  50.                 if(d>EPS) return 1;
  51.         }
  52.         return 0;
  53. }
  54.  
  55. //      METHODS OF class voronoi
  56.  
  57. void voronoi::disperse(object *o, cell *C) {
  58.         traverse++;
  59.         if(traverse==-1L) {                             // OVERFLOW
  60.                 traverse=-3L; disperse((object*)0,C);   // REINITIALIZE
  61.                 traverse=0L;
  62.         } 
  63.         cell *c=C; c->b=(cell*)0;                       // PARTICULAR CASE
  64.         cell *n=c->ln.first();                          // ACTUAL NEIGHBOR
  65.         register int back=0;                            // DON'T STEP BACK YET
  66.         for(;;) {                                       // ITERATIVE TRAVERSE
  67.                 c->t=traverse;                          // MARK AS TRAVERSED
  68.                 register int go=0;                      // DON'T GO YET
  69.                 do {                                    // ACTION ON ACTUAL c
  70.                         if(back) {                      // STEP BACKWARDS
  71.                                 cell *b=c->b;           // FROM WHERE WE CAME
  72.                                 if(o) c->lo+=o;         // PERFORM ACTION
  73.                                 c=b; back=0; continue;  // TAKE BACK CELL
  74.                         }
  75.                         if(n==c->b) continue;           // DON'T STEP BACK YET
  76.                         if(n->t==traverse) continue;    // ALREADY TRAVERSED
  77.                         if(o && !(*n&(*o))) continue;   // PARTIAL TRAVERSE
  78.                         go=1; break;                    // GO ON
  79.                 } while(n=c->ln.next());                // UNTIL NOT ALL DONE
  80.                 if(go) {                                // STEP FORWARDS
  81.                         n->b=c;                         // ESTABLISH BACK LINK
  82.                         c=n; n=c->ln.first(); continue; // LET'S GO AHEAD
  83.                 }
  84.                 if(c==C) break;                         // RETURNED
  85.                 back=1;                                 // GO BACK IF NO BETTER
  86.         }
  87.         if(o) C->lo+=o; C->t=traverse;                  // C IS THE LAST ONE
  88. }
  89.  
  90. void voronoi::preprocess(list<object*> *lo) {           // PREPROCESSING
  91.         list<cell*> *lc=new list<cell*>;
  92.         register int n=0;
  93.         for(object* o=lo->first(); o; o=lo->next()) {           // OBJECTS
  94.                 list<vector*> *lp=o->particles();
  95.                 for(vector* p=lp->first(); p; p=lp->next()) 
  96.                         {(*lc)+=new cell(*p,o); n++; delete p;}
  97.                 delete lp;
  98.         }
  99.         for(o=lightsources.first(); o; o=lightsources.next()) { // LIGHTSRC.S
  100.                 list<vector*> *lp=o->particles();
  101.                 for(vector* p=lp->first(); p; p=lp->next()) 
  102.                         {(*lc)+=new cell(*p); n++; delete p;}
  103.                 delete lp;
  104.         }
  105.         (*lc)+=new cell(actcamera.getray(0.,0.).o); n++;        // EYE
  106.         cell **C=new cell*[n];
  107.         n=0;
  108.         for(cell* c=lc->first(); c; c=lc->next()) C[n++]=c;
  109.         delete lc;
  110.         qsort((void*)C, (size_t)n, sizeof(cell*), lexcmp);  // LEXICOGRAPHIC
  111.         list<VECTOR<3>*> *lv=new list<VECTOR<3>*>;
  112.         (*lv)+=C[0]; cell* p=C[0];
  113.         vector cm=C[0]->p;                              // CENTER OF MASS
  114.         for(register int i=1; i<n; i++) {
  115.                 cm=cm+C[i]->p;
  116.                 if(lexcmp((const void**)&C[i],          // MULTIPLE PARTICLE
  117.                           (const void**)&p)==0) {
  118.                         for(o=C[i]->lh.first(); o; o=C[i]->lh.next())
  119.                                 p->lh+=o;
  120.                         delete C[i];
  121.                 } else {(*lv)+=C[i]; p=C[i];}           // SINGLE (YET)
  122.         }
  123.         cm=cm/(double)n;                                // CENTER OF MASS
  124.         delete C;
  125.         VECTOR<3>*b[4]; for(i=0;i<4;i++) b[i]=new cell; // BOUNDARY CELLS
  126.         VORONOI<3> V(lv,b);                             // VERTEX STRUCTURE
  127.         for(i=0;i<4;i++) {                              // BOUNDARY CELLS
  128.                 cell* c=(cell*)b[i];
  129.                 c->p=vector((*c)[0],(*c)[1],(*c)[2]);   // UPDATE
  130.         }
  131.         V(initcells);                                   // 3D VORONOI-CELLS
  132.         traverse=0L;                                    // INITIALIZE disperse
  133.         double ddmin; register int first=1;             // INITIALIZE SEARCH
  134.         for(c=(cell*)lv->first();c;c=(cell*)lv->next()) {
  135.                 for(o=c->lh.first();o;o=c->lh.next())   // BUILD OBJECT LISTS
  136.                         disperse(o,c);                  // PARTIAL TRAVERSE
  137.                 vector u=c->p-cm; double dd=u%u;        // DISTANCE^2 FROM cm
  138.                 if(first || dd<ddmin)                   // FIND CLOSEST TO cm
  139.                         {ddmin=dd;this->C=c;}
  140.                 first=0;
  141.         }
  142.         n=1<<(lmax+1);
  143.         s=new cell*[n]; for(i=0;i<n;i++) s[i]=this->C;  // START FOR firstlist
  144.         delete lv; // ???
  145. }
  146.  
  147. void voronoi::step() {                                  // ONE STEP ALONG RAY r
  148.         double tmin=0.; cell* nmin=(cell*)0;            // INITIALIZE SEARCH
  149.         for(cell*n=a->ln.first(); n; n=a->ln.next()) {  // NEIGHBORING CELLS
  150.                 vector u=n->p-a->p;                     // NORMAL OF BISECTOR
  151.                 double d=r->d%u;                        // DENOMINATOR
  152.                 if(fabs(d)<EPS) continue;               // PARALLEL FACE
  153.                 vector v=(n->p+a->p)*.5-r->o;
  154.                 double t=(v%u)/d;                       // t AT INTERSECTION
  155.                 if(t<=this->t) continue;                // WOULD BE A BACK STEP
  156.                 if(t<EPS)                               // r->o ON BISECTOR PL.
  157.                         if(u%r->d>0.) {tmin=EPS; nmin=n; continue;}
  158.                         else continue;
  159.                 if(!tmin || t<tmin) {tmin=t; nmin=n;}   // CLOSER INTERSECTION
  160.         }
  161.         this->t=tmin; a=nmin;
  162. }
  163.  
  164. list<object*>* voronoi::firstlist(ray& r) {
  165.         register int i=r.c>=0?r.c:(1<<(lmax+1))-1;      // INDEX INTO s
  166.         a=s[i];                                         // INITIALIZE step()
  167.         if(!((*a)&r.o)) {                               // COHERENCE FAILED
  168.                 ray R(a->p,norm(r.o-a->p),0,0);         // PATH OF WALK
  169.                 this->r=&R; t=0.;                       // INITIALIZE step()
  170.                 for(step(); a; step())                  // WALK
  171.                         if((*a)&r.o) {s[i]=a; break;}   // SUCCESS
  172.         }
  173.         this->r=&r; t=0.;                               // INITIALIZE step()
  174.         return a?&a->lo:(list<object*>*)0;
  175. }
  176.  
  177. list<object*>* voronoi::nextlist() {
  178.         step();                                         // ONE STEP
  179.         return a?&a->lo:(list<object*>*)0;              // SUCCESS OR FAIL
  180. }
  181.