voronoi.h

[Game.Programming].Academic - Graphics Gems (6 books source code)

                for(register int i=0; i<D+1; i++) this->b[i]=b[i];
                build(l);
        }
        void operator()(void (*f)(VERTEX<D>* v));       // TRAVERSE VERTICES
        ~VORONOI() {(*this)(free);}                     // TRAVERSE AND DELETE
};

template <int D> void VORONOI<D>::bound(list<VECTOR<D>*>* l) {
        register int i, j;

//      NORMAL VECTORS FOR FACES OF BOUNDING SIMPLEX

        double a[D+1][D]; normals(D, a); 
        VECTOR<D> n[D+1]; 
        for(i=0; i<D+1; i++) n[i]=VECTOR<D>(a[i])/sqrt(ll(VECTOR<D>(a[i])));

//      MAXIMAL DISTANCES IN DIRECTION OF NORMALS

        double d, dmax[D+1], dmin[D+1]; register int first=1;
        for(VECTOR<D>* p=l->first(); p; p=l->next()) {
                for(i=0; i<D+1; i++) {
                        d=n[i]*(*p);
                        if(first || d>dmax[i]) dmax[i]=d;
                        if(first || d<dmin[i]) dmin[i]=d;
                }
                first=0;
        }

//      VERTICES OF BOUNDING SIMPLEX (INTERSECT FACES CYCLICALLY)

        for(i=0; i<D+1; i++) dmax[i]+=(dmax[i]-dmin[i])*.5;     // INACCURACY
        VECTOR<D> A[D]; double t[D];
        for(i=0; i<D+1; i++) {
                for(j=0; j<D; j++) {
                        A[j]=n[(i+j)%(D+1)];
                        t[j]=dmax[(i+j)%(D+1)];
                }
                (void)MATRIX<D>(A)(*b[i],VECTOR<D>(t)); // EQUATION A*b[i]=t
        }

//      CENTRAL VERTEX AND CENTROID

        c=new VERTEX<D>(b);
        for(i=0; i<D+1; i++) C=C+(*b[i]);
        C=C*(1./(double)(D+1));
}

template <int D> void VORONOI<D>::find() {
        register int i, j;
        double P[D+1][D+1]; for(j=0; j<D+1; j++) P[D][j]=1.;
        double q[D+1]; for(j=0; j<D; j++) q[j]=(*this->q)[j]; q[D]=1.;
        VERTEX<D> *v=c;
        VECTOR<D+1> a;                          // BARICENTRIC COORDINATES OF q
        for(;;) {
                for(i=0; i<D; i++) for(j=0; j<D+1; j++) P[i][j]=(*v->p[j])[i];
                (void)MATRIX<D+1>(P)(a,VECTOR<D+1>(q));         // SOLVE P*a=q
                double aminus=0.;
                for(j=0; j<D+1; j++) if(a[j]<aminus) {aminus=a[j]; i=j;}
                if(aminus<0.) {v=v->v[i]; continue;}
                break;                                          // q INSIDE
        }
        ld=new list<VERTEX<D>*>; *ld+=v;
}

template <int D> void VORONOI<D>::search() {
        VERTEX<D> *vstart=ld->first(), *v=vstart; v->b=-1;
        register int back=0;
        for(;;) {
                register int go=0, i; VERTEX<D> *n;
                do {
                        if(back) {                      // STEP BACKWARDS
                                *ld+=v; 
                                i=v->b; v->b=-1; v=v->v[i]; back=0; continue;
                        }
                        if(v->i==v->b) continue;
                        if(v->v[v->i]==(VERTEX<D>*)0) continue;
                        n=v->v[v->i];                   // NEIGHBOR
                        if(n->b>=0) continue;           // ALREADY TRAVERSED
                        if((*ld)[n]) continue;          // ALREADY ON LIST
                        if(dd(*q,n->c)<n->rr) go=1;     // GO IF q IN SPHERE
                        if(go) break;
                } while((v->i=(v->i+1)%(D+1))!=0);
                if(go) {                                // STEP FORWARDS
                        for(i=0; i<D+1; i++)            // COMPUTE BACK INDEX
                                if(v==n->v[i]) {n->b=i; break;}
                        v=n; continue;
                }
                if(v==vstart) break;
                back=1;
        }    
}

template <int D> void VORONOI<D>::create() {
        VERTEX<D> *v;
        ln=new list<VERTEX<D>*>;
        for(v=ld->first(); v; v=ld->next()) {               // TAKE VERTICES
                for(register int i=0; i<D+1; i++) {         // TAKE NEIGHBORS
                        VERTEX<D> *n=v->v[i];
                        if((*ld)[n]) continue;              // ALSO DELETED
                        VERTEX<D> *m=new VERTEX<D>(q,v,i);  // POINT q + RING i
                        if(m->rr<0.) 
                                {delete m;*ld+=n;continue;} // DEGENERACY
                        *ln+=m;                             // STORE
                        register int j;
                        for(j=0; j<D+1; j++) 
                                if(m->p[j]==q) m->v[j]=n;   // OUTER LINK
                        if(n==(VERTEX<D>*)0) continue;      // NO REAL NEIGHBOR
                        for(j=0; j<D+1; j++) 
                                if(n->v[j]==v) n->v[j]=m;   // BACK LINK
                }
        }
        if((*ld)[c]) {                                      // NEW c NEEDED
                double d, ddmin; register int first=1;
                for(v=ln->first(); v; v=ln->next()) {
                        d=dd(v->c,C);
                        if(first || d<ddmin) {c=v; ddmin=d;}
                        first=0;
                }
        }
        for(v=ld->first(); v; v=ld->next()) {delete v;} // DELETE VERTICES
        delete ld;
}

template <int D> void VORONOI<D>::link() {
        register int i, j, n, iv, iw;
        VERTEX<D> *v, *w;
        n=0; for(v=ln->first(); v; v=ln->next()) n++;
        VERTEX<D> *N[n];
        n=0; for(v=ln->first(); v; v=ln->next()) N[n++]=v;
        for(i=0; i<n-1; i++) {
                v=N[i];
                for(j=i+1; j<n; j++) {
                        w=N[j];
                        for(iv=0; iv<D+1; iv++) {
                                for(iw=0; iw<D+1; iw++) {
                                        if(samering(v,iv,w,iw))
                                                {v->v[iv]=w; w->v[iw]=v;}
                                }
                        }
                }
        }
        delete ln;
}

template <int D> void VORONOI<D>::operator()(void (*f)(VERTEX<D>* v)) {
        traverse++;
        if(traverse==-1L) {                             // OVERFLOW
                traverse=-3L; (*this)(donothing);       // REINITIALIZE
                traverse=0L;
        } 
        VERTEX<D> *v=c; v->b=D+1;                       // PARTICULAR CASE
        register int back=0;
        for(;;) {                                       // ITERATIVE TRAVERSE
                v->t=traverse;                          // MARK AS TRAVERSED
                register int go=0;                      // DON'T GO YET
                VERTEX<D> *n;                           // ACTUAL NEIGHBOR
                do {                                    // ACTION ON ACTUAL v
                        if(back) {                      // STEP BACKWARDS
                                VERTEX<D>*n=v->v[v->b]; // FROM WHERE WE CAME
                                v->b=-1;                // FOR NEXT USAGE
                                f(v);                   // PERFORM ACTION
                                v=n; back=0; continue;  // TAKE BACK VERTEX
                        }
                        if(v->i==v->b) continue;        // DON'T STEP BACK YET
                        if(v->v[v->i]==(VERTEX<D>*)0) 
                                continue;               // NO REAL NEIGHBOR
                        n=v->v[v->i];                   // WHERE WE SHOULD GO
                        if(n->t==traverse) continue;    // ALREADY TRAVERSED
                        go=1; break;                    // GO ON
                } while((v->i=(v->i+1)%(D+1))!=0);      // UNTIL NOT ALL DONE
                if(go) {                                // STEP FORWARDS
                        for(register int i=0;i<D+1;i++) // FIND BACK LINK
                                if(v==n->v[i])
                                        {n->b=i;break;} // BOOK
                        v=n; continue;                  // LET'S GO
                }
                if(v==c) break;                         // RETURNED
                back=1;                                 // GO BACK IF NO BETTER
        }
        f(c); c->b=-1; c->t=traverse;                   // THE LAST ONE
}