home *** CD-ROM | disk | FTP | other *** search
- //========================================================================
- //
- // Poly.cc
- //
- // Copyright 1999 Emmanuel Lesueur
- //
- // A basic polygon manipulation package. This is experimental.
- // There are probably much better algorythms...
- // There is plenty of room for optimizations.
- //
- //========================================================================
-
- #ifdef __GNUC__
- #pragma implementation
- #endif
-
- // the "export" keyword not yet supported by compilers
- #define export
-
- #include <stdio.h>
- #include <stdlib.h>
- #include "poly.h"
-
- #undef STR
- #define STR(x) #x
- #define MYASSERT(x) {if(!(x)) {printf(__FILE__"/%d\n",__LINE__);throw("Internal Error");};}
-
- #if 0
-
- # define DB(x) x
- template<class T> void showp(const Polygon<T>& p) {
- printf("poly(%d):",p.size());
- for(typename Polygon<T>::const_iterator i(p.begin());i!=p.end();++i)
- printf(" (%g,%g)",double(i->x),double(i->y));
- printf("\n");
- }
- # if 0
- extern bool show_steps;
- void clear();
- void draw(int color,int x0,int y0,const Vector<double>&);
- void draw(int color,int x0,int y0,const Polygon<double>&);
- void draw(int color,int x0,int y0,const PArea<double>&);
- bool wait();
- # define CLEAR {if(show_steps) clear();}
- # define DRAW(c,x) {if(show_steps) draw(c,370,250,(x));}
- # define WAIT {if(show_steps && wait()) return res;}
- # else
- # define CLEAR
- # define DRAW(c,x)
- # define WAIT
- # endif
-
- #else
-
- # define DB(x)
-
- #endif
-
- export template<class T> void array<T>::destroy(iterator p,iterator q) {
- while(q!=p)
- (--q)->~T();
- }
-
- export template<class T> void array<T>::destroy() {
- if(d && --d->count==0) {
- destroy(pbeg,pcur);
- myfree(d);
- }
- }
-
- export template<class T> void array<T>::copy(iterator p,const_iterator q,const_iterator r) {
- iterator s=p;
- try {
- while(q!=r) {
- new (p) T(*q);
- ++p;
- ++q;
- }
- }
- catch(...) {
- destroy(s,p);
- rethrow;
- }
- }
-
- export template<class T> void array<T>::reserve(size_type sz) {
- if(pend-pbeg<sz || (d && d->count>1)) {
- data* d1=static_cast<data*>(myalloc(sizeof(data)-sizeof(T)+sizeof(T)*sz));
- d1->count=1;
- try {
- if(pbeg!=pcur)
- copy(d1->buf,pbeg,pcur);
- }
- catch(...) {
- myfree(d1);
- rethrow;
- }
- destroy();
- d=d1;
- pcur=d1->buf+(pcur-pbeg);
- pbeg=d1->buf;
- pend=pbeg+sz;
- }
- }
-
- export template<class T> void array<T>::resize(size_type sz) {
- size_type s=size();
- size_type t=s/2;
- if(sz>=t)
- t+=sz;
- if(t<10)
- t+=10;
- reserve(s+t);
- }
-
- export template<class T> void array<T>::move_forward(iterator p,iterator q,int n) {
- iterator r=q;
- try {
- while(q!=p) {
- --q;
- new (q+n) T(*q);
- }
- }
- catch(...) {
- destroy(q+n+1,r+n);
- destroy(p,q);
- rethrow;
- }
- }
-
- export template<class T> void array<T>::move_backward(iterator p,iterator q,int n) {
- iterator r=p;
- try {
- while(p!=q) {
- new (p-n) T(*p);
- ++p;
- }
- }
- catch(...) {
- destroy(r-n,p-n);
- destroy(p,q);
- rethrow;
- }
- }
-
- export template<class T> void array<T>::insert(const_iterator p,const T& x) {
- int k=p-pbeg;
- if(pcur==pend || d->count!=1)
- resize(1);
- try {
- move_forward(pbeg+k,pcur,1);
- }
- catch(...) {
- pcur=pbeg+k;
- rethrow;
- }
- ++pcur;
- try {
- new (pbeg+k) T(x);
- }
- catch(...) {
- destroy(pbeg+k+1,pcur);
- pcur=pbeg+k;
- rethrow;
- }
- }
-
- export template<class T> void array<T>::erase(const_iterator p) {
- int k=p-pbeg;
- resize(0);
- pbeg[k].~T();
- try {
- move_backward(pbeg+k+1,pcur,1);
- }
- catch(...) {
- pcur=pbeg+k;
- rethrow;
- }
- --pcur;
- }
-
- /*
- * Class Line.
- */
-
- export template<class T> Vertex<T> Line<T>::intersection(const Line<T>& l) const {
- T d=b*l.a-a*l.b;
- return Vertex<T>((c*l.b-b*l.c)/d,(a*l.c-c*l.a)/d);
- }
-
- export template<class T> int Line<T>::side(const Vertex<T>& v) const {
- T d=a*v.x+b*v.y+c;
- return d<-0.01?-1:(d>0.01?1:0);
- }
-
- /*
- * Class Polygon.
- */
-
- export template<class T> Vertex<T> Polygon<T>::center() const {
- Vertex<T> r(0,0);
- if(!empty()) {
- for(const_iterator i(begin());i!=end();++i) {
- r.x+=i->x;
- r.y+=i->y;
- }
- r.x/=size();
- r.y/=size();
- }
- return r;
- }
-
- export template<class T> bool Polygon<T>::is_rectangle() const {
- return (size()==4 ||
- (size()==5 && (*this)[0].x==(*this)[4].x &&
- (*this)[0].y==(*this)[4].y)) &&
- (((*this)[0].x==(*this)[1].x &&
- (*this)[2].x==(*this)[3].x &&
- (*this)[0].y==(*this)[3].y &&
- (*this)[1].y==(*this)[2].y) ||
- ((*this)[1].x==(*this)[2].x &&
- (*this)[0].x==(*this)[3].x &&
- (*this)[0].y==(*this)[1].y &&
- (*this)[2].y==(*this)[3].y));
- }
-
- export template<class T> void Polygon<T>::reduce() {
- if(!empty() && !compare(begin()->x,back().x) &&
- !compare(begin()->y,back().y))
- pop_back();
- }
-
-
- /*
- * Split a polygon along a line.
- */
-
- template<class T> class splitter {
- friend bool PArea<T>::split(PArea<T>&,PArea<T>&,const Line<T>&) const;
-
- struct vinfo {
- Vertex<T> vertex;
- signed char side;
- signed char real_side;
- };
- typedef array<vinfo> varray;
- typedef array<varray> parray;
- typedef typename array<vinfo>::const_iterator viterator;
- typedef typename array<varray>::const_iterator piterator;
-
- struct iinfo {
- piterator poly;
- viterator vertex;
- int dir; // 1 if going from - to +, -1 otherwise.
- bool inside;
- bool pos_side_done;
- bool neg_side_done;
- };
-
- typedef array<iinfo> iarray;
- typedef typename array<iinfo>::const_iterator iiterator;
-
- void follow(PArea<T>& res,piterator p,viterator q,
- int dir,int side,bool,int,viterator stop);
-
- iarray inter;
- };
-
- template<class T> void splitter<T>::follow(PArea<T>& r,
- piterator p,
- viterator q,
- int dir,
- int side,
- bool odd,
- int inum,
- viterator stop) {
- DB(printf("following until (%g,%g,%d)\n",double(stop->vertex.x),double(stop->vertex.y),stop->side);)
- DB(printf("(%g,%g,%d)\n",double(q->vertex.x),double(q->vertex.y),q->side);)
- Polygon<T> res;
- res.push_back(q->vertex);
- bool is_non_trivial=false;
- viterator first(p->begin());
- viterator last(p->end());
- viterator lastq=stop;
- viterator lastlastq;
- viterator q0=q;
- typename iarray::modifier interm(inter);
- int check=inter.size();
- while(true) {
- lastlastq=lastq;
- lastq=q;
- if(dir==1) {
- ++q;
- if(q==last)
- q=first;
- } else {
- if(q==first)
- q=last;
- --q;
- }
- if(q->real_side!=0)
- is_non_trivial=true;
- else {
- if(lastq->side==0 && lastlastq->side==0) {
- DB(printf("...removing last\n");)
- res.pop_back();
- lastq=lastlastq;
- }
- }
- DB(printf("...(%g,%g,%d)\n",double(q->vertex.x),double(q->vertex.y),q->side);)
- res.push_back(q->vertex);
- if(q->side==0) {
- MYASSERT(check-->=0);
- int i=0;
- while(inter[i].vertex!=q) {
- ++i;
- MYASSERT(i<inter.size());
- }
- DB(printf("this is intersection %d\n",i);)
- if(side==1)
- interm[i].pos_side_done=true;
- else
- interm[i].neg_side_done=true;
- if(q==stop)
- break;
- if(i>inum) {
- int k=inum;
- while(k<i && inter[k].inside)
- ++k;
- if(k==i)
- break;
- }
- if(!(i&1)) {
-
- int check=inter.size();
- do {
- ++i;
- if(i==inter.size()) // might happen for winding rule
- i=0; // and self-intersecting polygons.
- MYASSERT(check-->=0);
- } while((side==1 && inter[i].pos_side_done) ||
- (side==-1 && inter[i].neg_side_done));
- } else {
-
- int check=inter.size();
- do {
- if(i==0) // idem.
- i=(int)inter.size()-1;
- else
- --i;
- MYASSERT(check-->=0);
- } while((side==1 && inter[i].pos_side_done) ||
- (side==-1 && inter[i].neg_side_done));
- }
- if(side==1)
- interm[i].pos_side_done=true;
- else
- interm[i].neg_side_done=true;
- p=inter[i].poly;
- first=p->begin();
- last=p->end();
- lastlastq=lastq;
- lastq=q;
- q=inter[i].vertex;
- dir=inter[i].dir*side;
- DB(printf("next=%d, dir=%d\n",i,dir);)
- if(q==stop)
- break;
- if(q->real_side!=0)
- is_non_trivial=true;
- else {
- if(lastq->side==0 && lastlastq->side==0) {
- DB(printf("...removing last\n");)
- res.pop_back();
- lastq=lastlastq;
- }
- }
- DB(printf("...(%g,%g,%d)\n",double(q->vertex.x),double(q->vertex.y),q->side);)
- res.push_back(q->vertex);
- }
- }
- if(q0->real_side!=0)
- is_non_trivial=true;
- else {
- if(lastq->side==0 && lastlastq->side==0) {
- DB(printf("...removing last\n");)
- res.pop_back();
- }
- }
- if(is_non_trivial) {
- DB(printf("follow: added ");showp(res);)
- r.push_back(res);
- }
- }
-
- export template<class T> bool PArea<T>::split(PArea<T>& pos,PArea<T>& neg,const Line<T>& line) const {
-
- DB(printf("*** splitting along (%g,%g,%g)\n",double(line.a),double(line.b),double(line.c));)
-
- splitter<T> sp;
-
- typedef typename splitter<T>::parray parray;
- typedef typename splitter<T>::varray varray;
- typedef typename splitter<T>::iarray iarray;
- typedef typename splitter<T>::viterator viterator;
- typedef typename splitter<T>::iiterator iiterator;
- typedef typename splitter<T>::vinfo vinfo;
- typedef typename splitter<T>::iinfo iinfo;
-
- parray area;
- array<signed char> side;
- area.reserve(size());
- bool has_true_inter=false;
-
- for(const_iterator p(begin());p!=end();++p) {
- DB(printf("split1: ");showp(*p);)
- if(p->empty())
- continue;
- vinfo lastv;
- lastv.vertex=p->end()[-1];
- lastv.side=(signed char)line.side(lastv.vertex);
- varray poly;
-
- // Get the list of all vertices, with their sides.
- // Add all intersection points.
- bool has_pos_side=false;
- bool has_neg_side=false;
- for(typename Polygon<T>::const_iterator i(p->begin());i!=p->end();++i) {
- vinfo v;
- v.vertex=*i;
- v.real_side=(signed char)line.side(v.vertex);
- v.side=v.real_side;
- if(v.side==1)
- has_pos_side=true;
- else if(v.side==-1)
- has_neg_side=true;
- if(v.side!=0 && v.side+lastv.side==0) {
- vinfo vi;
- vi.vertex=line.intersection(Line<T>(lastv.vertex,v.vertex));
- vi.side=0;
- vi.real_side=0;
- DB(printf("(%g,%g,0)\n",double(vi.vertex.x),double(vi.vertex.y));)
- poly.push_back(vi);
- }
- DB(printf("(%g,%g,%d)\n",double(v.vertex.x),double(v.vertex.y),v.side);)
- poly.push_back(v);
- lastv=v;
- }
- if(has_pos_side && has_neg_side)
- has_true_inter=true;
-
- int k=0;
- size_t n=poly.size();
- while(k<n && poly[k].side==0)
- ++k;
- if(k!=n) {
- DB(printf("normalizing\n");)
- int k0=k;
- // ... + 0 ... 0 - ... -> ... + 0 - ...
- // ... + 0 ... 0 + ... -> ... + 0 0 + ...
- signed char s1=poly[k].side;
- signed char s2=s1;
- {
- typename varray::modifier polym(poly);
- do {
- int l=k;
- ++k;
- if(k==n)
- k=0;
- signed char s=poly[k].side;
- //printf("l=%d k=%d s=%d s1=%d s2=%d\n",l,k,s,s1,s2);
- if(s==0) {
- if(s1==0)
- polym[k].side=(signed char)-s2;
- else
- s2=s1;
- } else {
- if(s1==0 && s==s2) {
- if(poly[l].side==0) {
- // |/ | /
- // + == | |
- // |\ | \
- //printf("-> (%g,%g,%d) !!\n",poly[l].vertex.x,poly[l].vertex.y,poly[l].side);
- //printf("l=%d k=%d s=%d s1=%d s2=%d\n",l,k,s,s1,s2);
- polym[l].side=s2;
- //printf("-> (%g,%g,%d) !!\n",poly[l].vertex.x,poly[l].vertex.y,poly[l].side);
- } else
- polym[l].side=0;
- }
- }
- s1=s;
- } while (k!=k0);
- DB(printf("poly(%d):",poly.size());
- for(viterator i(poly.begin());i!=poly.end();++i)
- printf(" (%g,%g,%d)",double(i->vertex.x),double(i->vertex.y),i->side);
- printf("\n");)
- }
-
- area.push_back(poly);
- // Collect intersection points.
- bool has_inter=false;
- s1=poly[n-1].side;
- for(int i=0;i<n;++i) {
- if(poly[i].side==0) {
- has_inter=true;
- iinfo t;
- t.poly=area.end()-1;
- t.vertex=t.poly->begin()+i;
- if(s1==0) {
- if(i==n-1)
- t.dir=poly[0].side;
- else
- t.dir=poly[i+1].side;
- } else
- t.dir=-s1;
-
- // insert sorted.
- size_t l=sp.inter.size();
- iiterator q(sp.inter.begin());
- while(l) {
- size_t k=l/2;
- iiterator q1(q+k);
- T tst=line.a*(t.vertex->vertex.y-q1->vertex->vertex.y)-
- line.b*(t.vertex->vertex.x-q1->vertex->vertex.x);
- if(tst>0) {
- q=q1+1;
- l-=k+1;
- } else if(tst==0) {
- q=q1;
- break;
- } else
- l=k;
- }
- sp.inter.insert(q,t);
- }
- s1=poly[i].side;
- }
-
- DB(printf("inter(%d):",sp.inter.size());
- for(iiterator i(sp.inter.begin());i!=sp.inter.end();++i)
- printf(" (%g,%g,%d)",double(i->vertex->vertex.x),double(i->vertex->vertex.y),i->dir);
- printf("\n");)
-
-
- if(!has_inter) {
- side.push_back(0);
- area.pop_back();
- if(has_neg_side) {
- DB(printf("adding to - side: ");showp(*p);)
- neg.push_back(*p);
- } else {
- DB(printf("adding to + side: ");showp(*p);)
- pos.push_back(*p);
- }
- } else {
- MYASSERT(!(sp.inter.size()&1));
- side.push_back((signed char)(has_neg_side?-1:1));
- }
- } else {
- // The polygon is degenerated and on the line.
- //pos.push_back(*p);
- side.push_back(0);
- }
- }
-
- if(has_true_inter) {
- int num_segs=sp.inter.size()-1;
- int w=0;
- typename iarray::modifier m(sp.inter);
- for(int k=0;k<num_segs+1;++k) {
- w+=m[k].dir;
- m[k].inside=!(k&1) || (rule()==winding && w);
- m[k].pos_side_done=false;
- m[k].neg_side_done=false;
- }
- for(int k=0;k<num_segs;++k) {
- iinfo& i=m[k];
- if(i.inside) {
- if(!i.pos_side_done) {
- DB(printf("following + side of intersection %d\n",k);)
- i.pos_side_done=true;
- sp.follow(pos,i.poly,i.vertex,
- i.dir,1,k&1,k,sp.inter[k+1].vertex);
- }
- if(!i.neg_side_done) {
- DB(printf("following - side of intersection %d\n",k);)
- i.neg_side_done=true;
- sp.follow(neg,i.poly,i.vertex,
- -i.dir,-1,k&1,k,sp.inter[k+1].vertex);
- }
- }
- }
- return true;
- } else {
- DB(printf("No real intersection.\n");)
- array<signed char>::const_iterator q(side.begin());
- for(const_iterator p(begin());p!=end() && q!=side.end();++p)
- if(!p->empty()) {
- if(*q) {
- if(*q==1) {
- DB(printf("adding to + side: ");showp(*p);)
- pos.push_back(*p);
- } else {
- DB(printf("adding to - side: ");showp(*p);)
- neg.push_back(*p);
- }
- }
- ++q;
- }
- return false;
- }
- }
-
- /*
- * Transform a PArea into a union of interiors of convex polygons.
- */
- // Issue: it does not work correctly for areas bounded by
- // non-convex polygons with empty interior.
-
- export template<class T> PArea<T> PArea<T>::make_convex_union() const {
-
- if(is_convex_union())
- return *this;
-
- PArea<T> res;
- res.rule(rule());
-
- if(!empty()) {
- array<PArea<T> > stack;
-
- stack.push_back(*this);
- do {
- PArea<T> a(stack.back());
- stack.pop_back();
-
- DB(CLEAR;DRAW(1,a);printf("----pop\n");WAIT;)
-
- bool splitted=false;
- for(const_iterator poly(a.begin());!splitted && poly!=a.end();++poly)
- if(!poly->empty()) {
- typename Polygon<T>::const_iterator v(poly->end()-1);
- typename Polygon<T>::const_iterator u(poly->begin());
- while(!splitted && u!=poly->end()) {
- if(compare(u->x,v->x) || compare(u->y,v->y)) {
-
- DB(CLEAR;DRAW(1,a);DRAW(2,Vector<T>(*v,*u));printf("----splitting\n");WAIT;)
-
- Line<T> line(*v,*u);
- PArea<T> pos;
- PArea<T> neg;
- pos.rule(rule());
- neg.rule(rule());
- a.split(pos,neg,line);
- if(!pos.empty() && !neg.empty()) {
- DB(CLEAR;DRAW(1,pos);DRAW(2,neg);printf("----splitted\n");WAIT;)
- stack.push_back(pos);
- stack.push_back(neg);
- splitted=true;
- } else {
- v=u;
- ++u;
- }
- } else {
- v=u;
- ++u;
- }
- }
- }
- if(!splitted) {
- // 'a' is a convex area. Its border is made
- // of exactly one polygon.
- const_iterator p(a.begin());
- while(p!=a.end() && p->empty())
- ++p;
- if(p!=a.end()) { // just in case...
- Polygon<T> x(*p);
- x.is_convex(true);
- res.push_back(x);
- }
- }
- } while(!stack.empty());
- }
- res.is_convex_union(true);
- return res;
- }
-
- export template<class T> PArea<T> operator & (const PArea<T>& a,const PArea<T>& b) {
- PArea<T> a1(a.make_convex_union());
- PArea<T> res;
-
- for(typename PArea<T>::const_iterator i(a1.begin());i!=a1.end();++i) {
- if(!i->empty()) {
- PArea<T> c(b);
- Vertex<T> o(i->center());
- typename Polygon<T>::const_iterator v(i->end()-1);
- typename Polygon<T>::const_iterator u(i->begin());
- while(u!=i->end()) {
- if(compare(u->x,v->x) || compare(u->y,v->y)) {
- Line<T> line(*v,*u);
- PArea<T> pos;
- PArea<T> neg;
- pos.rule(b.rule());
- neg.rule(b.rule());
- c.split(pos,neg,line);
- if(line.side(o)==-1)
- c=neg;
- else
- c=pos;
- }
- v=u;
- ++u;
- }
- for(typename PArea<T>::const_iterator j(c.begin());j!=c.end();++j)
- res.push_back(*j);
- }
- }
- return res;
- }
-
- export template<class T> void PArea<T>::reduce() {
- typename array<Polygon<T> >::modifier m(*this);
- for(int k=0;k<size();++k)
- m[k].reduce();
- }
-
- export template<class T> bool PArea<T>::is_rectangle_union() const {
- for(const_iterator p(begin());p!=end();++p)
- if(!p->is_rectangle())
- return false;
- return true;
- }
-
- // not really implemented...
- export template<class T> bool PArea<T>::is_disjoint_rectangle_union() const {
- return size()==1 && begin()->is_rectangle();
- }
-
-
- // until compilers support the "export" keyword...
-
- template int Line<double>::side(const Vertex<double>&) const;
- template Vertex<double> Line<double>::intersection(const Line<double>&) const;
- template void Polygon<double>::reduce();
- template Vertex<double> Polygon<double>::center() const;
- template bool PArea<double>::split(PArea<double>&,PArea<double>&,const Line<double>&) const;
- template PArea<double> PArea<double>::make_convex_union() const;
- template PArea<double> operator & (const PArea<double>&,const PArea<double>&);
- template void PArea<double>::reduce();
- template void splitter<double>::follow(PArea<double>&,piterator,viterator,int,int,bool,int,viterator);
-
- template void array<Vertex<double> >::reserve(size_type);
- template void array<Vertex<double> >::resize(size_type);
- template void array<Vertex<double> >::destroy();
- template void array<Polygon<double> >::reserve(size_type);
- template void array<Polygon<double> >::resize(size_type);
- template void array<Polygon<double> >::destroy();
- template void array<PArea<double> >::reserve(size_type);
- template void array<PArea<double> >::resize(size_type);
- template void array<PArea<double> >::destroy();
-
- template void array<splitter<double>::iinfo>::reserve(size_type);
- template void array<splitter<double>::iinfo>::resize(size_type);
- template void array<splitter<double>::iinfo>::destroy();
- //template void array<splitter<double>::iinfo>::destroy(iterator,iterator);
- template void array<splitter<double>::vinfo>::reserve(size_type);
- template void array<splitter<double>::vinfo>::resize(size_type);
- template void array<splitter<double>::vinfo>::destroy();
- template void array<array<splitter<double>::vinfo> >::reserve(size_type);
- template void array<array<splitter<double>::vinfo> >::resize(size_type);
- template void array<array<splitter<double>::vinfo> >::destroy();
-
- template void array<Vertex<short> >::reserve(size_type);
- template void array<Vertex<short> >::resize(size_type);
- template void array<Vertex<short> >::destroy();
- template void array<Polygon<short> >::reserve(size_type);
- template void array<Polygon<short> >::resize(size_type);
- template void array<Polygon<short> >::destroy();
- template bool PArea<short>::is_disjoint_rectangle_union() const;
- template bool PArea<short>::is_rectangle_union() const;
- template bool Polygon<short>::is_rectangle() const;
- //template void PArea<short>::reduce();
-
- template void array<signed char>::reserve(size_type);
- template void array<signed char>::resize(size_type);
- template void array<signed char>::destroy();
-
-