home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch5_5 / voronoi.h < prev   
C/C++ Source or Header  |  1995-03-04  |  16KB  |  365 lines
  1. //      PERMUTATION TEMPLATE (USED IN GAUSSIAN ELIMINATION)
  2.  
  3. template <int D> class PERMUTATION {
  4.         int n[D];                                               // ELEMENTS
  5. public:
  6.         PERMUTATION() {for(register int i=0; i<D; i++) n[i]=i;}
  7.         int operator[](int i) {return n[i];}
  8.         void operator()(int i, int j) {                         // SWAP
  9.                 if(i==j) return;
  10.                 register int t=n[i]; n[i]=n[j]; n[j]=t;
  11.         }
  12. };
  13.  
  14. //      D-DIMENSIONAL VECTOR TEMPLATE
  15.  
  16. template <int D> class VECTOR {
  17.         friend ostream& operator<<(ostream& o, VECTOR<D>& v);
  18.         double x[D];                                            // COORDINATES
  19. public:
  20.         VECTOR() {for(register int i=0; i<D; i++) x[i]=0.;}
  21.         VECTOR(double x[D]) {for(register int i=0; i<D; i++) this->x[i]=x[i];}
  22.         VECTOR(double x[D], PERMUTATION<D>& p) {
  23.                 for(register int i=0; i<D; i++) this->x[i]=x[p[i]];
  24.         }
  25.         double operator[](int i) {return x[i];}
  26.         void operator-=(VECTOR<D>& v) {
  27.                 for(register int i=0; i<D; i++) x[i]-=v.x[i];
  28.         }
  29.         VECTOR<D> operator*(double d) {
  30.                 VECTOR<D> w; for(register int i=0; i<D; i++) w.x[i]=x[i]*d;
  31.                 return w;
  32.         }
  33.         VECTOR<D> operator/(double d) {
  34.                 VECTOR<D> w; for(register int i=0; i<D; i++) w.x[i]=x[i]/d;
  35.                 return w;
  36.         }
  37.         double operator*(VECTOR<D>& v) {
  38.                 double d=0.; for(register int i=0; i<D; i++) d+=x[i]*v.x[i];
  39.                 return d;
  40.         }
  41.         VECTOR<D> operator+(VECTOR<D>& v) {
  42.                 VECTOR<D> w; for(register int i=0; i<D;i++) w.x[i]=x[i]+v.x[i];
  43.                 return w;
  44.         }
  45.         VECTOR<D> operator-(VECTOR<D>& v) {
  46.                 VECTOR<D> w; for(register int i=0; i<D;i++) w.x[i]=x[i]-v.x[i];
  47.                 return w;
  48.         }
  49. };
  50.  
  51. //      D-DIMENSIONAL SQUARE MATRIX TEMPLATE
  52.  
  53. template <int D> class MATRIX {
  54.         friend ostream& operator<<(ostream& o, MATRIX<D>& A);
  55.         VECTOR<D> a[D];                                         // ROWS
  56. public:
  57.         MATRIX(VECTOR<D> a[D]) {for(register int i=0;i<D;i++) this->a[i]=a[i];}
  58.         MATRIX(double a[D][D]) {
  59.                 for(register int i=0; i<D; i++) this->a[i]=VECTOR<D>(a[i]);
  60.         }
  61.         VECTOR<D> operator*(VECTOR<D>& x) {
  62.                 double y[D];
  63.                 for(register int i=0; i<D; i++) y[i]=a[i]*x;
  64.                 return VECTOR<D>(y);
  65.         }
  66.         int operator()(VECTOR<D>& x, VECTOR<D>& b) {    // SOLVE (*this)x=b
  67. //              GAUSSIAN ELIMINATION METHOD
  68.                 const double EPS=1e-10;
  69.                 VECTOR<D> B[D]; double c[D]; PERMUTATION<D> p; 
  70.                 register int i, j, k;
  71.                 for(i=0; i<D; i++) {B[i]=a[i]; c[i]=b[i];}      // COPY
  72.                 for(i=0; i<D; i++) {                            // THROUGH ROWS
  73.                         double a, amax=0., e, emain;
  74.                         for(j=i; j<D; j++)                      // MAIN ELEMENT
  75.                                 if((a=fabs(e=B[p[j]][i]))>amax)
  76.                                         {emain=e; amax=a; k=j;}
  77.                         if(amax<EPS) return 0;                  // SINGULAR
  78.                         p(i,k);                                 // SWAP
  79.                         for(j=i+1; j<D; j++) {                  // NULL BELOW
  80.                                 double s=B[p[j]][i]/emain;
  81.                                 B[p[j]]-=B[p[i]]*s;
  82.                                 c[p[j]]-=c[p[i]]*s;
  83.                         }
  84.                 }
  85.                 for(i=D-1; i>=0; i--) {                 // BUILD SOLUTION
  86.                         for(j=D-1; j>i; j--) c[p[i]]-=B[p[i]][j]*c[p[j]];
  87.                         c[p[i]]/=B[p[i]][i];
  88.                 }
  89.                 x=VECTOR<D>(c,p); return 1;
  90.         }
  91. };
  92.  
  93. //      VORONOI-VERTEX TEMPLATE
  94.  
  95. template <int D> struct VERTEX {
  96.         VECTOR<D>* p[D+1];                              // FORMING POINTS
  97.         VERTEX<D>* v[D+1];                              // CONTIGUOUS VERTICES
  98.         VECTOR<D> c;                                    // POSITION
  99.         double rr;                                      // RADIUS SQUARE
  100.         int b;                                          // BACK INDEX (WORK)
  101.         int i;                                          // ACTUAL INDEX (WORK)
  102.         long t;                                         // TRAVERSE CODE (WORK)
  103. private:
  104.         void initialize(VECTOR<D>* f[D+1]) {
  105.                 register int i;
  106.                 for(i=0; i<D+1; i++) {p[i]=f[i]; v[i]=(VERTEX<D>*)0;}
  107.                 this->b=-1; this->i=0; this->t=0L;
  108.                 VECTOR<D> A[D]; double b[D];
  109.                 for(i=0; i<D; i++) {
  110.                         A[i]=(*f[i+1])-(*f[i]);
  111.                         b[i]=(((*f[i+1])+(*f[i]))*0.5)*A[i];
  112.                 }
  113.                 if(!MATRIX<D>(A)(c,VECTOR<D>(b))) {     // EQUATION A*c=b
  114.                         rr=-1.;                         // DEGENERATE
  115.                         return;
  116.                 }
  117.                 rr=(*p[0]-c)*(*p[0]-c);
  118.         }
  119. public:
  120.         VERTEX(VECTOR<D>* f[D+1]) {initialize(f);}      // FORMING POINTS
  121.         VERTEX(VECTOR<D> *q, VERTEX<D> *v, int i) {     // POINT q AND RING i
  122.                 VECTOR<D> *f[D+1]; f[0]=q;
  123.                 for(register int j=1; j<D+1; j++) f[j]=v->p[(j+i)%(D+1)];
  124.                 initialize(f);
  125.         }
  126. };
  127.  
  128. //      VORONOI-DIAGRAM TEMPLATE
  129.  
  130. template <int D> class VORONOI {
  131.         friend ostream& operator<<(ostream& o, VORONOI<D>& v);
  132.         VECTOR<D> C;                                    // CENTROID
  133.         VECTOR<D> *b[D+1], B[D+1];                      // BOUNDING SIMPLEX
  134.         VERTEX<D> *c;                                   // CLOSEST TO CENTROID
  135.         double ll(VECTOR<D>& v) {return v*v;}           // LENGTH SQUARE
  136.         double dd(VECTOR<D>& v, VECTOR<D>& w) {         // DISTANCE SQUARE
  137.                 VECTOR<D> d=v-w; return d*d;
  138.         }
  139.         void normals(int d, double n[D+1][D]) { // NORMAL VECTORS FOR bound()
  140.                 register int i, j;
  141.                 if(d==2) {
  142.                         n[0][0]=1.0; n[0][1]=1.0;
  143.                         n[1][0]=-1.0; n[1][1]=1.0;
  144.                         n[2][0]=0.0; n[2][1]=-1.0;
  145.                         return;
  146.                 }
  147.                 normals(d-1, n);
  148.                 for(i=0; i<d; i++) n[i][d-1]=1.0;
  149.                 for(i=0; i<d-1; i++) n[d][i]=0.0;
  150.                 n[d][d-1]=-1.0;
  151.         }
  152.         void bound(list<VECTOR<D>*>* l);        // BUILD BOUNDING SIMPLEX
  153.         VECTOR<D>* q;                           // ACTUAL POINT
  154.         list<VERTEX<D>*> *ld;                   // VERTICES TO DELETE
  155.         void find();                            // FIND A VERTEX TO DELETE
  156.         void search();                          // FIND ALL VERTICES TO DELETE
  157.         list<VERTEX<D>*> *ln;                   // NEW VERTICES
  158.         void create();                          // CREATE NEW VERTICES
  159.         int samering(VERTEX<D>*v, int iv, VERTEX<D>*w, int iw) {
  160.                 for(register int i=(iv+1)%(D+1);i!=iv;i=(i+1)%(D+1)) {
  161.                         VECTOR<D> *p=v->p[i];
  162.                         for(register int j=(iw+1)%(D+1);j!=iw;j=(j+1)%(D+1))
  163.                                 if(w->p[j]==p) {j=-1; break;}
  164.                         if(j>=0) return 0;
  165.                 }
  166.                 return 1;
  167.         }
  168.         void link();                            // LINK NEW VERTICES TO EACH O.
  169.         void build(list<VECTOR<D>*>* l) {       // DISJOINT PARTICLES
  170.                 traverse=0L;
  171.                 bound(l);
  172.                 for(VECTOR<D>* p=l->first(); p; p=l->next())
  173.                         {q=p; find(); search(); create(); link();}
  174.         }
  175.         long traverse;                          // TRAVERSE CODE
  176.         static void free(VERTEX<D>*v){delete v;}// FOR DESTRUCTOR
  177.         static void donothing(VERTEX<D>*v){}    // FOR REINITIALIZE traverse
  178. public:
  179.         VORONOI(list<VECTOR<D>*>* l) {
  180.                 for(register int i=0; i<D+1; i++) b[i]=&B[i];
  181.                 build(l);
  182.         }
  183.         VORONOI(list<VECTOR<D>*>* l, VECTOR<D> *b[D+1]) {
  184.                 for(register int i=0; i<D+1; i++) this->b[i]=b[i];
  185.                 build(l);
  186.         }
  187.         void operator()(void (*f)(VERTEX<D>* v));       // TRAVERSE VERTICES
  188.         ~VORONOI() {(*this)(free);}                     // TRAVERSE AND DELETE
  189. };
  190.  
  191. template <int D> void VORONOI<D>::bound(list<VECTOR<D>*>* l) {
  192.         register int i, j;
  193.  
  194. //      NORMAL VECTORS FOR FACES OF BOUNDING SIMPLEX
  195.  
  196.         double a[D+1][D]; normals(D, a); 
  197.         VECTOR<D> n[D+1]; 
  198.         for(i=0; i<D+1; i++) n[i]=VECTOR<D>(a[i])/sqrt(ll(VECTOR<D>(a[i])));
  199.  
  200. //      MAXIMAL DISTANCES IN DIRECTION OF NORMALS
  201.  
  202.         double d, dmax[D+1], dmin[D+1]; register int first=1;
  203.         for(VECTOR<D>* p=l->first(); p; p=l->next()) {
  204.                 for(i=0; i<D+1; i++) {
  205.                         d=n[i]*(*p);
  206.                         if(first || d>dmax[i]) dmax[i]=d;
  207.                         if(first || d<dmin[i]) dmin[i]=d;
  208.                 }
  209.                 first=0;
  210.         }
  211.  
  212. //      VERTICES OF BOUNDING SIMPLEX (INTERSECT FACES CYCLICALLY)
  213.  
  214.         for(i=0; i<D+1; i++) dmax[i]+=(dmax[i]-dmin[i])*.5;     // INACCURACY
  215.         VECTOR<D> A[D]; double t[D];
  216.         for(i=0; i<D+1; i++) {
  217.                 for(j=0; j<D; j++) {
  218.                         A[j]=n[(i+j)%(D+1)];
  219.                         t[j]=dmax[(i+j)%(D+1)];
  220.                 }
  221.                 (void)MATRIX<D>(A)(*b[i],VECTOR<D>(t)); // EQUATION A*b[i]=t
  222.         }
  223.  
  224. //      CENTRAL VERTEX AND CENTROID
  225.  
  226.         c=new VERTEX<D>(b);
  227.         for(i=0; i<D+1; i++) C=C+(*b[i]);
  228.         C=C*(1./(double)(D+1));
  229. }
  230.  
  231. template <int D> void VORONOI<D>::find() {
  232.         register int i, j;
  233.         double P[D+1][D+1]; for(j=0; j<D+1; j++) P[D][j]=1.;
  234.         double q[D+1]; for(j=0; j<D; j++) q[j]=(*this->q)[j]; q[D]=1.;
  235.         VERTEX<D> *v=c;
  236.         VECTOR<D+1> a;                          // BARICENTRIC COORDINATES OF q
  237.         for(;;) {
  238.                 for(i=0; i<D; i++) for(j=0; j<D+1; j++) P[i][j]=(*v->p[j])[i];
  239.                 (void)MATRIX<D+1>(P)(a,VECTOR<D+1>(q));         // SOLVE P*a=q
  240.                 double aminus=0.;
  241.                 for(j=0; j<D+1; j++) if(a[j]<aminus) {aminus=a[j]; i=j;}
  242.                 if(aminus<0.) {v=v->v[i]; continue;}
  243.                 break;                                          // q INSIDE
  244.         }
  245.         ld=new list<VERTEX<D>*>; *ld+=v;
  246. }
  247.  
  248. template <int D> void VORONOI<D>::search() {
  249.         VERTEX<D> *vstart=ld->first(), *v=vstart; v->b=-1;
  250.         register int back=0;
  251.         for(;;) {
  252.                 register int go=0, i; VERTEX<D> *n;
  253.                 do {
  254.                         if(back) {                      // STEP BACKWARDS
  255.                                 *ld+=v; 
  256.                                 i=v->b; v->b=-1; v=v->v[i]; back=0; continue;
  257.                         }
  258.                         if(v->i==v->b) continue;
  259.                         if(v->v[v->i]==(VERTEX<D>*)0) continue;
  260.                         n=v->v[v->i];                   // NEIGHBOR
  261.                         if(n->b>=0) continue;           // ALREADY TRAVERSED
  262.                         if((*ld)[n]) continue;          // ALREADY ON LIST
  263.                         if(dd(*q,n->c)<n->rr) go=1;     // GO IF q IN SPHERE
  264.                         if(go) break;
  265.                 } while((v->i=(v->i+1)%(D+1))!=0);
  266.                 if(go) {                                // STEP FORWARDS
  267.                         for(i=0; i<D+1; i++)            // COMPUTE BACK INDEX
  268.                                 if(v==n->v[i]) {n->b=i; break;}
  269.                         v=n; continue;
  270.                 }
  271.                 if(v==vstart) break;
  272.                 back=1;
  273.         }    
  274. }
  275.  
  276. template <int D> void VORONOI<D>::create() {
  277.         VERTEX<D> *v;
  278.         ln=new list<VERTEX<D>*>;
  279.         for(v=ld->first(); v; v=ld->next()) {               // TAKE VERTICES
  280.                 for(register int i=0; i<D+1; i++) {         // TAKE NEIGHBORS
  281.                         VERTEX<D> *n=v->v[i];
  282.                         if((*ld)[n]) continue;              // ALSO DELETED
  283.                         VERTEX<D> *m=new VERTEX<D>(q,v,i);  // POINT q + RING i
  284.                         if(m->rr<0.) 
  285.                                 {delete m;*ld+=n;continue;} // DEGENERACY
  286.                         *ln+=m;                             // STORE
  287.                         register int j;
  288.                         for(j=0; j<D+1; j++) 
  289.                                 if(m->p[j]==q) m->v[j]=n;   // OUTER LINK
  290.                         if(n==(VERTEX<D>*)0) continue;      // NO REAL NEIGHBOR
  291.                         for(j=0; j<D+1; j++) 
  292.                                 if(n->v[j]==v) n->v[j]=m;   // BACK LINK
  293.                 }
  294.         }
  295.         if((*ld)[c]) {                                      // NEW c NEEDED
  296.                 double d, ddmin; register int first=1;
  297.                 for(v=ln->first(); v; v=ln->next()) {
  298.                         d=dd(v->c,C);
  299.                         if(first || d<ddmin) {c=v; ddmin=d;}
  300.                         first=0;
  301.                 }
  302.         }
  303.         for(v=ld->first(); v; v=ld->next()) {delete v;} // DELETE VERTICES
  304.         delete ld;
  305. }
  306.  
  307. template <int D> void VORONOI<D>::link() {
  308.         register int i, j, n, iv, iw;
  309.         VERTEX<D> *v, *w;
  310.         n=0; for(v=ln->first(); v; v=ln->next()) n++;
  311.         VERTEX<D> *N[n];
  312.         n=0; for(v=ln->first(); v; v=ln->next()) N[n++]=v;
  313.         for(i=0; i<n-1; i++) {
  314.                 v=N[i];
  315.                 for(j=i+1; j<n; j++) {
  316.                         w=N[j];
  317.                         for(iv=0; iv<D+1; iv++) {
  318.                                 for(iw=0; iw<D+1; iw++) {
  319.                                         if(samering(v,iv,w,iw))
  320.                                                 {v->v[iv]=w; w->v[iw]=v;}
  321.                                 }
  322.                         }
  323.                 }
  324.         }
  325.         delete ln;
  326. }
  327.  
  328. template <int D> void VORONOI<D>::operator()(void (*f)(VERTEX<D>* v)) {
  329.         traverse++;
  330.         if(traverse==-1L) {                             // OVERFLOW
  331.                 traverse=-3L; (*this)(donothing);       // REINITIALIZE
  332.                 traverse=0L;
  333.         } 
  334.         VERTEX<D> *v=c; v->b=D+1;                       // PARTICULAR CASE
  335.         register int back=0;
  336.         for(;;) {                                       // ITERATIVE TRAVERSE
  337.                 v->t=traverse;                          // MARK AS TRAVERSED
  338.                 register int go=0;                      // DON'T GO YET
  339.                 VERTEX<D> *n;                           // ACTUAL NEIGHBOR
  340.                 do {                                    // ACTION ON ACTUAL v
  341.                         if(back) {                      // STEP BACKWARDS
  342.                                 VERTEX<D>*n=v->v[v->b]; // FROM WHERE WE CAME
  343.                                 v->b=-1;                // FOR NEXT USAGE
  344.                                 f(v);                   // PERFORM ACTION
  345.                                 v=n; back=0; continue;  // TAKE BACK VERTEX
  346.                         }
  347.                         if(v->i==v->b) continue;        // DON'T STEP BACK YET
  348.                         if(v->v[v->i]==(VERTEX<D>*)0) 
  349.                                 continue;               // NO REAL NEIGHBOR
  350.                         n=v->v[v->i];                   // WHERE WE SHOULD GO
  351.                         if(n->t==traverse) continue;    // ALREADY TRAVERSED
  352.                         go=1; break;                    // GO ON
  353.                 } while((v->i=(v->i+1)%(D+1))!=0);      // UNTIL NOT ALL DONE
  354.                 if(go) {                                // STEP FORWARDS
  355.                         for(register int i=0;i<D+1;i++) // FIND BACK LINK
  356.                                 if(v==n->v[i])
  357.                                         {n->b=i;break;} // BOOK
  358.                         v=n; continue;                  // LET'S GO
  359.                 }
  360.                 if(v==c) break;                         // RETURNED
  361.                 back=1;                                 // GO BACK IF NO BETTER
  362.         }
  363.         f(c); c->b=-1; c->t=traverse;                   // THE LAST ONE
  364. }
  365.