polypartition.cpp
Go to the documentation of this file.
00001 //Copyright (C) 2011 by Ivan Fratric
00002 //
00003 //Permission is hereby granted, free of charge, to any person obtaining a copy
00004 //of this software and associated documentation files (the "Software"), to deal
00005 //in the Software without restriction, including without limitation the rights
00006 //to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
00007 //copies of the Software, and to permit persons to whom the Software is
00008 //furnished to do so, subject to the following conditions:
00009 //
00010 //The above copyright notice and this permission notice shall be included in
00011 //all copies or substantial portions of the Software.
00012 //
00013 //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
00014 //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00015 //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
00016 //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
00017 //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
00018 //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
00019 //THE SOFTWARE.
00020 
00021 
00022 #include <stdio.h>
00023 #include <string.h>
00024 #include <math.h>
00025 #include <list>
00026 #include <algorithm>
00027 #include <set>
00028 
00029 using namespace std;
00030 
00031 #include <polypartition/polypartition.h>
00032 
00033 #define TPPL_VERTEXTYPE_REGULAR 0
00034 #define TPPL_VERTEXTYPE_START 1
00035 #define TPPL_VERTEXTYPE_END 2
00036 #define TPPL_VERTEXTYPE_SPLIT 3
00037 #define TPPL_VERTEXTYPE_MERGE 4
00038 
00039 TPPLPoly::TPPLPoly() { 
00040         hole = false;
00041         numpoints = 0;
00042         points = NULL;
00043 }
00044 
00045 TPPLPoly::~TPPLPoly() {
00046         if(points) delete [] points;
00047 }
00048 
00049 void TPPLPoly::Clear() {
00050         if(points) delete [] points;
00051         hole = false;
00052         numpoints = 0;
00053         points = NULL;
00054 }
00055 
00056 void TPPLPoly::Init(long numpoints) {
00057         Clear();
00058         this->numpoints = numpoints;
00059         points = new TPPLPoint[numpoints];
00060 }
00061 
00062 void TPPLPoly::Triangle(TPPLPoint &p1, TPPLPoint &p2, TPPLPoint &p3) {
00063         Init(3);
00064         points[0] = p1;
00065         points[1] = p2;
00066         points[2] = p3;
00067 }
00068 
00069 TPPLPoly::TPPLPoly(const TPPLPoly &src) {
00070         hole = src.hole;
00071         numpoints = src.numpoints;
00072         points = new TPPLPoint[numpoints];
00073         memcpy(points, src.points, numpoints*sizeof(TPPLPoint));
00074 }
00075 
00076 TPPLPoly& TPPLPoly::operator=(const TPPLPoly &src) {
00077         Clear();
00078         hole = src.hole;
00079         numpoints = src.numpoints;
00080         points = new TPPLPoint[numpoints];
00081         memcpy(points, src.points, numpoints*sizeof(TPPLPoint));
00082         return *this;
00083 }
00084 
00085 int TPPLPoly::GetOrientation() {
00086         long i1,i2;
00087         tppl_float area = 0;
00088         for(i1=0; i1<numpoints; i1++) {
00089                 i2 = i1+1;
00090                 if(i2 == numpoints) i2 = 0;
00091                 area += points[i1].x * points[i2].y - points[i1].y * points[i2].x;
00092         }
00093         if(area>0) return TPPL_CCW;
00094         if(area<0) return TPPL_CW;
00095         return 0;
00096 }
00097 
00098 void TPPLPoly::SetOrientation(int orientation) {
00099         int polyorientation = GetOrientation();
00100         if(polyorientation&&(polyorientation!=orientation)) {
00101                 Invert();
00102         }
00103 }
00104 
00105 void TPPLPoly::Invert() {
00106         long i;
00107         TPPLPoint *invpoints;
00108 
00109         invpoints = new TPPLPoint[numpoints];
00110         for(i=0;i<numpoints;i++) {
00111                 invpoints[i] = points[numpoints-i-1];
00112         }
00113 
00114         delete [] points;
00115         points = invpoints;
00116 }
00117 
00118 TPPLPoint TPPLPartition::Normalize(const TPPLPoint &p) {
00119         TPPLPoint r;
00120         tppl_float n = sqrt(p.x*p.x + p.y*p.y);
00121         if(n!=0) {
00122                 r = p/n;
00123         } else {
00124                 r.x = 0;
00125                 r.y = 0;
00126         }
00127         return r;
00128 }
00129 
00130 tppl_float TPPLPartition::Distance(const TPPLPoint &p1, const TPPLPoint &p2) {
00131         tppl_float dx,dy;
00132         dx = p2.x - p1.x;
00133         dy = p2.y - p1.y;
00134         return(sqrt(dx*dx + dy*dy));
00135 }
00136 
00137 //checks if two lines intersect
00138 int TPPLPartition::Intersects(TPPLPoint &p11, TPPLPoint &p12, TPPLPoint &p21, TPPLPoint &p22) {
00139         if((p11.x == p21.x)&&(p11.y == p21.y)) return 0;
00140         if((p11.x == p22.x)&&(p11.y == p22.y)) return 0;
00141         if((p12.x == p21.x)&&(p12.y == p21.y)) return 0;
00142         if((p12.x == p22.x)&&(p12.y == p22.y)) return 0;
00143 
00144         TPPLPoint v1ort,v2ort,v;
00145         tppl_float dot11,dot12,dot21,dot22;
00146 
00147         v1ort.x = p12.y-p11.y;
00148         v1ort.y = p11.x-p12.x;
00149 
00150         v2ort.x = p22.y-p21.y;
00151         v2ort.y = p21.x-p22.x;
00152 
00153         v = p21-p11;
00154         dot21 = v.x*v1ort.x + v.y*v1ort.y;
00155         v = p22-p11;
00156         dot22 = v.x*v1ort.x + v.y*v1ort.y;
00157 
00158         v = p11-p21;
00159         dot11 = v.x*v2ort.x + v.y*v2ort.y;
00160         v = p12-p21;
00161         dot12 = v.x*v2ort.x + v.y*v2ort.y;
00162 
00163         if(dot11*dot12>0) return 0;
00164         if(dot21*dot22>0) return 0;
00165 
00166         return 1;
00167 }
00168 
00169 //removes holes from inpolys by merging them with non-holes
00170 int TPPLPartition::RemoveHoles(list<TPPLPoly> *inpolys, list<TPPLPoly> *outpolys) {
00171         list<TPPLPoly> polys;
00172         list<TPPLPoly>::iterator holeiter,polyiter,iter,iter2;
00173         long i,i2,holepointindex,polypointindex;
00174         TPPLPoint holepoint,polypoint,bestpolypoint;
00175         TPPLPoint linep1,linep2;
00176         TPPLPoint v1,v2;
00177         TPPLPoly newpoly;
00178         bool hasholes;
00179         bool pointvisible;
00180         bool pointfound;
00181         
00182         //check for trivial case (no holes)
00183         hasholes = false;
00184         for(iter = inpolys->begin(); iter!=inpolys->end(); iter++) {
00185                 if(iter->IsHole()) {
00186                         hasholes = true;
00187                         break;
00188                 }
00189         }
00190         if(!hasholes) {
00191                 for(iter = inpolys->begin(); iter!=inpolys->end(); iter++) {
00192                         outpolys->push_back(*iter);
00193                 }
00194                 return 1;
00195         }
00196 
00197         polys = *inpolys;
00198 
00199         while(1) {
00200                 //find the hole point with the largest x
00201                 hasholes = false;
00202                 for(iter = polys.begin(); iter!=polys.end(); iter++) {
00203                         if(!iter->IsHole()) continue;
00204 
00205                         if(!hasholes) {
00206                                 hasholes = true;
00207                                 holeiter = iter;
00208                                 holepointindex = 0;
00209                         }
00210 
00211                         for(i=0; i < iter->GetNumPoints(); i++) {
00212                                 if(iter->GetPoint(i).x > holeiter->GetPoint(holepointindex).x) {
00213                                         holeiter = iter;
00214                                         holepointindex = i;
00215                                 }
00216                         }
00217                 }
00218                 if(!hasholes) break;
00219                 holepoint = holeiter->GetPoint(holepointindex);
00220                 
00221                 pointfound = false;
00222                 for(iter = polys.begin(); iter!=polys.end(); iter++) {
00223                         if(iter->IsHole()) continue;
00224                         for(i=0; i < iter->GetNumPoints(); i++) {
00225                                 if(iter->GetPoint(i).x <= holepoint.x) continue;
00226                                 if(!InCone(iter->GetPoint((i+iter->GetNumPoints()-1)%(iter->GetNumPoints())),
00227                                         iter->GetPoint(i),
00228                                         iter->GetPoint((i+1)%(iter->GetNumPoints())),
00229                                         holepoint)) 
00230                                         continue;
00231                                 polypoint = iter->GetPoint(i);
00232                                 if(pointfound) {
00233                                         v1 = Normalize(polypoint-holepoint);
00234                                         v2 = Normalize(bestpolypoint-holepoint);
00235                                         if(v2.x > v1.x) continue;                               
00236                                 }
00237                                 pointvisible = true;
00238                                 for(iter2 = polys.begin(); iter2!=polys.end(); iter2++) {
00239                                         if(iter2->IsHole()) continue;
00240                                         for(i2=0; i2 < iter2->GetNumPoints(); i2++) {
00241                                                 linep1 = iter2->GetPoint(i2);
00242                                                 linep2 = iter2->GetPoint((i2+1)%(iter2->GetNumPoints()));
00243                                                 if(Intersects(holepoint,polypoint,linep1,linep2)) {
00244                                                         pointvisible = false;
00245                                                         break;
00246                                                 }
00247                                         }
00248                                         if(!pointvisible) break;
00249                                 }
00250                                 if(pointvisible) {
00251                                         pointfound = true;
00252                                         bestpolypoint = polypoint;
00253                                         polyiter = iter;
00254                                         polypointindex = i;
00255                                 }
00256                         }
00257                 }
00258 
00259                 if(!pointfound) return 0;
00260 
00261                 newpoly.Init(holeiter->GetNumPoints() + polyiter->GetNumPoints() + 2);
00262                 i2 = 0;
00263                 for(i=0;i<=polypointindex;i++) {
00264                         newpoly[i2] = polyiter->GetPoint(i);
00265                         i2++;
00266                 }
00267                 for(i=0;i<=holeiter->GetNumPoints();i++) {
00268                         newpoly[i2] = holeiter->GetPoint((i+holepointindex)%holeiter->GetNumPoints());
00269                         i2++;
00270                 }
00271                 for(i=polypointindex;i<polyiter->GetNumPoints();i++) {
00272                         newpoly[i2] = polyiter->GetPoint(i);
00273                         i2++;
00274                 }
00275                 
00276                 polys.erase(holeiter);
00277                 polys.erase(polyiter);
00278                 polys.push_back(newpoly);
00279         }
00280 
00281         for(iter = polys.begin(); iter!=polys.end(); iter++) {
00282                 outpolys->push_back(*iter);
00283         }
00284         
00285         return 1;
00286 }
00287 
00288 bool TPPLPartition::IsConvex(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3) {
00289         tppl_float tmp;
00290         tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
00291         if(tmp>0) return 1;
00292         else return 0;
00293 }
00294 
00295 bool TPPLPartition::IsReflex(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3) {
00296         tppl_float tmp;
00297         tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
00298         if(tmp<0) return 1;
00299         else return 0;
00300 }
00301 
00302 bool TPPLPartition::IsInside(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3, TPPLPoint &p) {
00303         if(IsConvex(p1,p,p2)) return false;
00304         if(IsConvex(p2,p,p3)) return false;
00305         if(IsConvex(p3,p,p1)) return false;
00306         return true;
00307 }
00308 
00309 bool TPPLPartition::InCone(TPPLPoint &p1, TPPLPoint &p2, TPPLPoint &p3, TPPLPoint &p) {
00310         bool convex;
00311 
00312         convex = IsConvex(p1,p2,p3);
00313 
00314         if(convex) {
00315                 if(!IsConvex(p1,p2,p)) return false;
00316                 if(!IsConvex(p2,p3,p)) return false;
00317                 return true;
00318         } else {
00319                 if(IsConvex(p1,p2,p)) return true;
00320                 if(IsConvex(p2,p3,p)) return true;
00321                 return false;
00322         }
00323 }
00324 
00325 bool TPPLPartition::InCone(PartitionVertex *v, TPPLPoint &p) {
00326         TPPLPoint p1,p2,p3;
00327 
00328         p1 = v->previous->p;
00329         p2 = v->p;
00330         p3 = v->next->p;
00331 
00332         return InCone(p1,p2,p3,p);
00333 }
00334 
00335 void TPPLPartition::UpdateVertexReflexity(PartitionVertex *v) {
00336         PartitionVertex *v1,*v3;
00337         v1 = v->previous;
00338         v3 = v->next;
00339         v->isConvex = !IsReflex(v1->p,v->p,v3->p);      
00340 }
00341 
00342 void TPPLPartition::UpdateVertex(PartitionVertex *v, PartitionVertex *vertices, long numvertices) {
00343         long i;
00344         PartitionVertex *v1,*v3;
00345         TPPLPoint vec1,vec3;
00346 
00347         v1 = v->previous;
00348         v3 = v->next;
00349 
00350         v->isConvex = IsConvex(v1->p,v->p,v3->p);
00351 
00352         vec1 = Normalize(v1->p - v->p);
00353         vec3 = Normalize(v3->p - v->p);
00354         v->angle = vec1.x*vec3.x + vec1.y*vec3.y;
00355 
00356         if(v->isConvex) {
00357                 v->isEar = true;
00358                 for(i=0;i<numvertices;i++) {
00359                         if((vertices[i].p.x==v->p.x)&&(vertices[i].p.y==v->p.y)) continue;
00360                         if((vertices[i].p.x==v1->p.x)&&(vertices[i].p.y==v1->p.y)) continue;
00361                         if((vertices[i].p.x==v3->p.x)&&(vertices[i].p.y==v3->p.y)) continue;
00362                         if(IsInside(v1->p,v->p,v3->p,vertices[i].p)) {
00363                                 v->isEar = false;
00364                                 break;
00365                         }
00366                 }
00367         } else {
00368                 v->isEar = false;
00369         }
00370 }
00371 
00372 //triangulation by ear removal
00373 int TPPLPartition::Triangulate_EC(TPPLPoly *poly, list<TPPLPoly> *triangles) {
00374         long numvertices;
00375         PartitionVertex *vertices;
00376         PartitionVertex *ear;
00377         TPPLPoly triangle;
00378         long i,j;
00379         bool earfound;
00380 
00381         if(poly->GetNumPoints() < 3) return 0;
00382         if(poly->GetNumPoints() == 3) {
00383                 triangles->push_back(*poly);
00384                 return 1;
00385         }
00386 
00387         numvertices = poly->GetNumPoints();
00388 
00389         vertices = new PartitionVertex[numvertices];
00390         for(i=0;i<numvertices;i++) {
00391                 vertices[i].isActive = true;
00392                 vertices[i].p = poly->GetPoint(i);
00393                 if(i==(numvertices-1)) vertices[i].next=&(vertices[0]);
00394                 else vertices[i].next=&(vertices[i+1]);
00395                 if(i==0) vertices[i].previous = &(vertices[numvertices-1]);
00396                 else vertices[i].previous = &(vertices[i-1]);
00397         }
00398         for(i=0;i<numvertices;i++) {
00399                 UpdateVertex(&vertices[i],vertices,numvertices);
00400         }
00401 
00402         for(i=0;i<numvertices-3;i++) {
00403                 earfound = false;
00404                 //find the most extruded ear
00405                 for(j=0;j<numvertices;j++) {
00406                         if(!vertices[j].isActive) continue;
00407                         if(!vertices[j].isEar) continue;
00408                         if(!earfound) {
00409                                 earfound = true;
00410                                 ear = &(vertices[j]);
00411                         } else {
00412                                 if(vertices[j].angle > ear->angle) {
00413                                         ear = &(vertices[j]);                           
00414                                 }
00415                         }
00416                 }
00417                 if(!earfound) {
00418                         delete [] vertices;
00419                         return 0;
00420                 }
00421 
00422                 triangle.Triangle(ear->previous->p,ear->p,ear->next->p);
00423                 triangles->push_back(triangle);
00424 
00425                 ear->isActive = false;
00426                 ear->previous->next = ear->next;
00427                 ear->next->previous = ear->previous;
00428 
00429                 if(i==numvertices-4) break;
00430 
00431                 UpdateVertex(ear->previous,vertices,numvertices);
00432                 UpdateVertex(ear->next,vertices,numvertices);
00433         }
00434         for(i=0;i<numvertices;i++) {
00435                 if(vertices[i].isActive) {
00436                         triangle.Triangle(vertices[i].previous->p,vertices[i].p,vertices[i].next->p);
00437                         triangles->push_back(triangle);
00438                         break;
00439                 }
00440         }
00441 
00442         delete [] vertices;
00443 
00444         return 1;
00445 }
00446 
00447 int TPPLPartition::Triangulate_EC(list<TPPLPoly> *inpolys, list<TPPLPoly> *triangles) {
00448         list<TPPLPoly> outpolys;
00449         list<TPPLPoly>::iterator iter;
00450         
00451         if(!RemoveHoles(inpolys,&outpolys)) return 0;
00452         for(iter=outpolys.begin();iter!=outpolys.end();iter++) {
00453                 if(!Triangulate_EC(&(*iter),triangles)) return 0;
00454         }
00455         return 1;
00456 }
00457 
00458 int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, list<TPPLPoly> *parts) {
00459         list<TPPLPoly> triangles;
00460         list<TPPLPoly>::iterator iter1,iter2;
00461         TPPLPoly *poly1,*poly2;
00462         TPPLPoly newpoly;
00463         TPPLPoint d1,d2,p1,p2,p3;
00464         long i11,i12,i21,i22,i13,i23,j,k;
00465         bool isdiagonal;
00466         long numreflex;
00467 
00468         //check if the poly is already convex
00469         numreflex = 0;
00470         for(i11=0;i11<poly->GetNumPoints();i11++) {
00471                 if(i11==0) i12 = poly->GetNumPoints()-1;
00472                 else i12=i11-1;
00473                 if(i11==(poly->GetNumPoints()-1)) i13=0;
00474                 else i13=i11+1;
00475                 if(IsReflex(poly->GetPoint(i12),poly->GetPoint(i11),poly->GetPoint(i13))) {
00476                         numreflex = 1;
00477                         break;
00478                 }
00479         }
00480         if(numreflex == 0) {
00481                 parts->push_back(*poly);
00482                 return 1;
00483         }
00484 
00485         if(!Triangulate_EC(poly,&triangles)) return 0;
00486 
00487         for(iter1 = triangles.begin(); iter1 != triangles.end(); iter1++) {
00488                 poly1 = &(*iter1);
00489                 for(i11=0;i11<poly1->GetNumPoints();i11++) {
00490                         d1 = poly1->GetPoint(i11);
00491                         i12 = (i11+1)%(poly1->GetNumPoints());
00492                         d2 = poly1->GetPoint(i12);
00493 
00494                         isdiagonal = false;
00495                         for(iter2 = iter1; iter2 != triangles.end(); iter2++) {
00496                                 if(iter1 == iter2) continue;
00497                                 poly2 = &(*iter2);
00498 
00499                                 for(i21=0;i21<poly2->GetNumPoints();i21++) {
00500                                         if((d2.x != poly2->GetPoint(i21).x)||(d2.y != poly2->GetPoint(i21).y)) continue;
00501                                         i22 = (i21+1)%(poly2->GetNumPoints());
00502                                         if((d1.x != poly2->GetPoint(i22).x)||(d1.y != poly2->GetPoint(i22).y)) continue;
00503                                         isdiagonal = true;
00504                                         break;
00505                                 }
00506                                 if(isdiagonal) break;
00507                         }
00508 
00509                         if(!isdiagonal) continue;
00510 
00511                         p2 = poly1->GetPoint(i11);
00512                         if(i11 == 0) i13 = poly1->GetNumPoints()-1;
00513                         else i13 = i11-1;
00514                         p1 = poly1->GetPoint(i13);
00515                         if(i22 == (poly2->GetNumPoints()-1)) i23 = 0;
00516                         else i23 = i22+1;
00517                         p3 = poly2->GetPoint(i23);
00518 
00519                         if(!IsConvex(p1,p2,p3)) continue;
00520                         
00521                         p2 = poly1->GetPoint(i12);
00522                         if(i12 == (poly1->GetNumPoints()-1)) i13 = 0;
00523                         else i13 = i12+1;
00524                         p3 = poly1->GetPoint(i13);
00525                         if(i21 == 0) i23 = poly2->GetNumPoints()-1;
00526                         else i23 = i21-1;
00527                         p1 = poly2->GetPoint(i23);
00528                         
00529                         if(!IsConvex(p1,p2,p3)) continue;
00530 
00531                         newpoly.Init(poly1->GetNumPoints()+poly2->GetNumPoints()-2);
00532                         k = 0;
00533                         for(j=i12;j!=i11;j=(j+1)%(poly1->GetNumPoints())) {
00534                                 newpoly[k] = poly1->GetPoint(j);
00535                                 k++;
00536                         }
00537                         for(j=i22;j!=i21;j=(j+1)%(poly2->GetNumPoints())) {
00538                                 newpoly[k] = poly2->GetPoint(j);
00539                                 k++;
00540                         }
00541 
00542                         triangles.erase(iter2);
00543                         *iter1 = newpoly;
00544                         poly1 = &(*iter1);
00545                         i11 = -1;
00546 
00547                         continue;
00548                 }
00549         }
00550 
00551         for(iter1 = triangles.begin(); iter1 != triangles.end(); iter1++) {
00552                 parts->push_back(*iter1);
00553         }
00554 
00555         return 1;
00556 }
00557 
00558 int TPPLPartition::ConvexPartition_HM(list<TPPLPoly> *inpolys, list<TPPLPoly> *parts) {
00559         list<TPPLPoly> outpolys;
00560         list<TPPLPoly>::iterator iter;
00561         
00562         if(!RemoveHoles(inpolys,&outpolys)) return 0;
00563         for(iter=outpolys.begin();iter!=outpolys.end();iter++) {
00564                 if(!ConvexPartition_HM(&(*iter),parts)) return 0;
00565         }
00566         return 1;
00567 }
00568 
00569 //minimum-weight polygon triangulation by dynamic programming
00570 //O(n^3) time complexity
00571 //O(n^2) space complexity
00572 int TPPLPartition::Triangulate_OPT(TPPLPoly *poly, list<TPPLPoly> *triangles) {
00573         long i,j,k,gap,n;
00574         DPState **dpstates;
00575         TPPLPoint p1,p2,p3,p4;
00576         long bestvertex;
00577         tppl_float weight,minweight,d1,d2;
00578         Diagonal diagonal,newdiagonal;
00579         list<Diagonal> diagonals;
00580         TPPLPoly triangle;
00581         int ret = 1;
00582 
00583         n = poly->GetNumPoints();
00584         dpstates = new DPState *[n];
00585         for(i=1;i<n;i++) {
00586                 dpstates[i] = new DPState[i];
00587         }
00588 
00589         //init states and visibility
00590         for(i=0;i<(n-1);i++) {
00591                 p1 = poly->GetPoint(i);
00592                 for(j=i+1;j<n;j++) {
00593                         dpstates[j][i].visible = true;
00594                         dpstates[j][i].weight = 0;
00595                         dpstates[j][i].bestvertex = -1;
00596                         if(j!=(i+1)) {
00597                                 p2 = poly->GetPoint(j);
00598                                 
00599                                 //visibility check
00600                                 if(i==0) p3 = poly->GetPoint(n-1);
00601                                 else p3 = poly->GetPoint(i-1);
00602                                 if(i==(n-1)) p4 = poly->GetPoint(0);
00603                                 else p4 = poly->GetPoint(i+1);
00604                                 if(!InCone(p3,p1,p4,p2)) {
00605                                         dpstates[j][i].visible = false;
00606                                         continue;
00607                                 }
00608 
00609                                 if(j==0) p3 = poly->GetPoint(n-1);
00610                                 else p3 = poly->GetPoint(j-1);
00611                                 if(j==(n-1)) p4 = poly->GetPoint(0);
00612                                 else p4 = poly->GetPoint(j+1);
00613                                 if(!InCone(p3,p2,p4,p1)) {
00614                                         dpstates[j][i].visible = false;
00615                                         continue;
00616                                 }
00617 
00618                                 for(k=0;k<n;k++) {
00619                                         p3 = poly->GetPoint(k);
00620                                         if(k==(n-1)) p4 = poly->GetPoint(0);
00621                                         else p4 = poly->GetPoint(k+1);
00622                                         if(Intersects(p1,p2,p3,p4)) {
00623                                                 dpstates[j][i].visible = false;
00624                                                 break;
00625                                         }
00626                                 }
00627                         }
00628                 }
00629         }
00630         dpstates[n-1][0].visible = true;
00631         dpstates[n-1][0].weight = 0;
00632         dpstates[n-1][0].bestvertex = -1;
00633 
00634         for(gap = 2; gap<n; gap++) {
00635                 for(i=0; i<(n-gap); i++) {
00636                         j = i+gap;
00637                         if(!dpstates[j][i].visible) continue; 
00638                         bestvertex = -1;
00639                         for(k=(i+1);k<j;k++) {
00640                                 if(!dpstates[k][i].visible) continue; 
00641                                 if(!dpstates[j][k].visible) continue;
00642 
00643                                 if(k<=(i+1)) d1=0;
00644                                 else d1 = Distance(poly->GetPoint(i),poly->GetPoint(k));
00645                                 if(j<=(k+1)) d2=0;
00646                                 else d2 = Distance(poly->GetPoint(k),poly->GetPoint(j));
00647 
00648                                 weight = dpstates[k][i].weight + dpstates[j][k].weight + d1 + d2;
00649 
00650                                 if((bestvertex == -1)||(weight<minweight)) {
00651                                         bestvertex = k;
00652                                         minweight = weight;
00653                                 }
00654                         }
00655                         if(bestvertex == -1) {
00656                                 for(i=1;i<n;i++) {
00657                                         delete [] dpstates[i];
00658                                 }
00659                                 delete [] dpstates;
00660 
00661                                 return 0;
00662                         }
00663                         
00664                         dpstates[j][i].bestvertex = bestvertex;
00665                         dpstates[j][i].weight = minweight;
00666                 }
00667         }
00668 
00669         newdiagonal.index1 = 0;
00670         newdiagonal.index2 = n-1;
00671         diagonals.push_back(newdiagonal);
00672         while(!diagonals.empty()) {
00673                 diagonal = *(diagonals.begin());
00674                 diagonals.pop_front();
00675                 bestvertex = dpstates[diagonal.index2][diagonal.index1].bestvertex;
00676                 if(bestvertex == -1) {
00677                         ret = 0;                
00678                         break;
00679                 }
00680                 triangle.Triangle(poly->GetPoint(diagonal.index1),poly->GetPoint(bestvertex),poly->GetPoint(diagonal.index2));
00681                 triangles->push_back(triangle);
00682                 if(bestvertex > (diagonal.index1+1)) {
00683                         newdiagonal.index1 = diagonal.index1;
00684                         newdiagonal.index2 = bestvertex;
00685                         diagonals.push_back(newdiagonal);
00686                 }
00687                 if(diagonal.index2 > (bestvertex+1)) {
00688                         newdiagonal.index1 = bestvertex;
00689                         newdiagonal.index2 = diagonal.index2;
00690                         diagonals.push_back(newdiagonal);
00691                 }
00692         }
00693 
00694         for(i=1;i<n;i++) {
00695                 delete [] dpstates[i];
00696         }
00697         delete [] dpstates;
00698 
00699         return ret;
00700 }
00701 
00702 void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 **dpstates) {
00703         Diagonal newdiagonal;
00704         list<Diagonal> *pairs;
00705         long w2;
00706 
00707         w2 = dpstates[a][b].weight;
00708         if(w>w2) return;
00709         
00710         pairs = &(dpstates[a][b].pairs);
00711         newdiagonal.index1 = i;
00712         newdiagonal.index2 = j;
00713 
00714         if(w<w2) {
00715                 pairs->clear();
00716                 pairs->push_front(newdiagonal);
00717                 dpstates[a][b].weight = w;
00718         } else {
00719                 if((!pairs->empty())&&(i <= pairs->begin()->index1)) return;
00720                 while((!pairs->empty())&&(pairs->begin()->index2 >= j)) pairs->pop_front();
00721                 pairs->push_front(newdiagonal);
00722         }
00723 }
00724 
00725 void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) {
00726         list<Diagonal> *pairs;
00727         list<Diagonal>::iterator iter,lastiter;
00728         long top;
00729         long w;
00730 
00731         if(!dpstates[i][j].visible) return;
00732         top = j;
00733         w = dpstates[i][j].weight;
00734     if(k-j > 1) {
00735                 if (!dpstates[j][k].visible) return;
00736                 w += dpstates[j][k].weight + 1;
00737     }
00738         if(j-i > 1) {
00739                 pairs = &(dpstates[i][j].pairs);
00740                 iter = pairs->end();
00741                 lastiter = pairs->end();
00742                 while(iter!=pairs->begin()) {
00743                         iter--;
00744                         if(!IsReflex(vertices[iter->index2].p,vertices[j].p,vertices[k].p)) lastiter = iter;
00745                         else break;
00746                 }
00747                 if(lastiter == pairs->end()) w++;
00748                 else {
00749                         if(IsReflex(vertices[k].p,vertices[i].p,vertices[lastiter->index1].p)) w++;
00750                         else top = lastiter->index1;
00751                 }
00752         }
00753         UpdateState(i,k,w,top,j,dpstates);
00754 }
00755 
00756 void TPPLPartition::TypeB(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) {
00757         list<Diagonal> *pairs;
00758         list<Diagonal>::iterator iter,lastiter;
00759         long top;
00760         long w;
00761 
00762         if(!dpstates[j][k].visible) return;
00763         top = j;
00764         w = dpstates[j][k].weight;
00765         
00766         if (j-i > 1) {
00767                 if (!dpstates[i][j].visible) return;
00768                 w += dpstates[i][j].weight + 1;
00769         }
00770         if (k-j > 1) {
00771                 pairs = &(dpstates[j][k].pairs);
00772 
00773                 iter = pairs->begin();
00774                 if((!pairs->empty())&&(!IsReflex(vertices[i].p,vertices[j].p,vertices[iter->index1].p))) {
00775                         lastiter = iter;
00776                         while(iter!=pairs->end()) {
00777                                 if(!IsReflex(vertices[i].p,vertices[j].p,vertices[iter->index1].p)) {
00778                                         lastiter = iter;
00779                                         iter++;
00780                                 }
00781                                 else break;
00782                         }
00783                         if(IsReflex(vertices[lastiter->index2].p,vertices[k].p,vertices[i].p)) w++;
00784                         else top = lastiter->index2;
00785                 } else w++;
00786         }
00787         UpdateState(i,k,w,j,top,dpstates);
00788 }
00789 
00790 int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, list<TPPLPoly> *parts) {
00791         TPPLPoint p1,p2,p3,p4;
00792         PartitionVertex *vertices;
00793         DPState2 **dpstates;
00794         long i,j,k,n,gap;
00795         list<Diagonal> diagonals,diagonals2;
00796         Diagonal diagonal,newdiagonal;
00797         list<Diagonal> *pairs,*pairs2;
00798         list<Diagonal>::iterator iter,iter2;
00799         int ret;
00800         TPPLPoly newpoly;
00801         list<long> indices;
00802         list<long>::iterator iiter;
00803         bool ijreal,jkreal;
00804 
00805         n = poly->GetNumPoints();
00806         vertices = new PartitionVertex[n];
00807 
00808         dpstates = new DPState2 *[n];
00809         for(i=0;i<n;i++) {
00810                 dpstates[i] = new DPState2[n];
00811         }
00812 
00813         //init vertex information
00814         for(i=0;i<n;i++) {
00815                 vertices[i].p = poly->GetPoint(i);
00816                 vertices[i].isActive = true;
00817                 if(i==0) vertices[i].previous = &(vertices[n-1]);
00818                 else vertices[i].previous = &(vertices[i-1]);
00819                 if(i==(poly->GetNumPoints()-1)) vertices[i].next = &(vertices[0]);
00820                 else vertices[i].next = &(vertices[i+1]);
00821         }
00822         for(i=1;i<n;i++) {
00823                 UpdateVertexReflexity(&(vertices[i]));
00824         }
00825 
00826         //init states and visibility
00827         for(i=0;i<(n-1);i++) {
00828                 p1 = poly->GetPoint(i);
00829                 for(j=i+1;j<n;j++) {
00830                         dpstates[i][j].visible = true;
00831                         if(j==i+1) {
00832                                 dpstates[i][j].weight = 0;
00833                         } else {
00834                                 dpstates[i][j].weight = 2147483647;
00835                         }
00836                         if(j!=(i+1)) {
00837                                 p2 = poly->GetPoint(j);
00838                                 
00839                                 //visibility check
00840                                 if(!InCone(&vertices[i],p2)) {
00841                                         dpstates[i][j].visible = false;
00842                                         continue;
00843                                 }
00844                                 if(!InCone(&vertices[j],p1)) {
00845                                         dpstates[i][j].visible = false;
00846                                         continue;
00847                                 }
00848 
00849                                 for(k=0;k<n;k++) {
00850                                         p3 = poly->GetPoint(k);
00851                                         if(k==(n-1)) p4 = poly->GetPoint(0);
00852                                         else p4 = poly->GetPoint(k+1);
00853                                         if(Intersects(p1,p2,p3,p4)) {
00854                                                 dpstates[i][j].visible = false;
00855                                                 break;
00856                                         }
00857                                 }
00858                         }
00859                 }
00860         }
00861         for(i=0;i<(n-2);i++) {
00862                 j = i+2;
00863                 if(dpstates[i][j].visible) {
00864                         dpstates[i][j].weight = 0;
00865                         newdiagonal.index1 = i+1;
00866                         newdiagonal.index2 = i+1;
00867                         dpstates[i][j].pairs.push_back(newdiagonal);
00868                 }
00869         }
00870 
00871         dpstates[0][n-1].visible = true;
00872         vertices[0].isConvex = false; //by convention
00873 
00874         for(gap=3; gap<n; gap++) {
00875                 for(i=0;i<n-gap;i++) {
00876                         if(vertices[i].isConvex) continue;
00877                         k = i+gap;
00878                         if(dpstates[i][k].visible) {
00879                                 if(!vertices[k].isConvex) {
00880                                         for(j=i+1;j<k;j++) TypeA(i,j,k,vertices,dpstates);
00881                                 } else {
00882                                         for(j=i+1;j<(k-1);j++) {
00883                                                 if(vertices[j].isConvex) continue;                              
00884                                                 TypeA(i,j,k,vertices,dpstates);
00885                                         }
00886                                         TypeA(i,k-1,k,vertices,dpstates);
00887                                 }
00888                         }
00889                 }
00890                 for(k=gap;k<n;k++) {
00891                         if(vertices[k].isConvex) continue;
00892                         i = k-gap;
00893                         if((vertices[i].isConvex)&&(dpstates[i][k].visible)) {
00894                                 TypeB(i,i+1,k,vertices,dpstates);
00895                                 for(j=i+2;j<k;j++) {
00896                                         if(vertices[j].isConvex) continue;                              
00897                                         TypeB(i,j,k,vertices,dpstates);
00898                                 }
00899                         }
00900                 }
00901         }
00902 
00903 
00904         //recover solution
00905         ret = 1;
00906         newdiagonal.index1 = 0;
00907         newdiagonal.index2 = n-1;
00908         diagonals.push_front(newdiagonal);
00909         while(!diagonals.empty()) {
00910                 diagonal = *(diagonals.begin());
00911                 diagonals.pop_front();
00912                 if((diagonal.index2 - diagonal.index1) <=1) continue;
00913                 pairs = &(dpstates[diagonal.index1][diagonal.index2].pairs);
00914                 if(pairs->empty()) {
00915                         ret = 0;
00916                         break;
00917                 }
00918                 if(!vertices[diagonal.index1].isConvex) {
00919                         iter = pairs->end();
00920                         iter--;
00921                         j = iter->index2;
00922                         newdiagonal.index1 = j;
00923                         newdiagonal.index2 = diagonal.index2;
00924                         diagonals.push_front(newdiagonal);
00925                         if((j - diagonal.index1)>1) {
00926                                 if(iter->index1 != iter->index2) {
00927                                         pairs2 = &(dpstates[diagonal.index1][j].pairs);
00928                                         while(1) {
00929                                                 if(pairs2->empty()) {
00930                                                         ret = 0;
00931                                                         break;
00932                                                 }
00933                                                 iter2 = pairs2->end();
00934                                                 iter2--;
00935                                                 if(iter->index1 != iter2->index1) pairs2->pop_back();
00936                                                 else break;
00937                                         }
00938                                         if(ret == 0) break;
00939                                 }
00940                                 newdiagonal.index1 = diagonal.index1;
00941                                 newdiagonal.index2 = j;
00942                                 diagonals.push_front(newdiagonal);
00943                         }
00944                 } else {
00945                         iter = pairs->begin();
00946                         j = iter->index1;
00947                         newdiagonal.index1 = diagonal.index1;
00948                         newdiagonal.index2 = j;
00949                         diagonals.push_front(newdiagonal);
00950                         if((diagonal.index2 - j) > 1) {
00951                                 if(iter->index1 != iter->index2) {
00952                                         pairs2 = &(dpstates[j][diagonal.index2].pairs);
00953                                         while(1) {
00954                                                 if(pairs2->empty()) {
00955                                                         ret = 0;
00956                                                         break;
00957                                                 }
00958                                                 iter2 = pairs2->begin();
00959                                                 if(iter->index2 != iter2->index2) pairs2->pop_front();
00960                                                 else break;
00961                                         }
00962                                         if(ret == 0) break;
00963                                 }
00964                                 newdiagonal.index1 = j;
00965                                 newdiagonal.index2 = diagonal.index2;
00966                                 diagonals.push_front(newdiagonal);
00967                         }
00968                 }
00969         }
00970 
00971         if(ret == 0) {
00972                 for(i=0;i<n;i++) {
00973                         delete [] dpstates[i];
00974                 }
00975                 delete [] dpstates;
00976                 delete [] vertices;
00977 
00978                 return ret;
00979         }
00980 
00981         newdiagonal.index1 = 0;
00982         newdiagonal.index2 = n-1;
00983         diagonals.push_front(newdiagonal);
00984         while(!diagonals.empty()) {
00985                 diagonal = *(diagonals.begin());
00986                 diagonals.pop_front();
00987                 if((diagonal.index2 - diagonal.index1) <= 1) continue;
00988                 
00989                 indices.clear();
00990                 diagonals2.clear();
00991                 indices.push_back(diagonal.index1);
00992                 indices.push_back(diagonal.index2);
00993                 diagonals2.push_front(diagonal);
00994 
00995                 while(!diagonals2.empty()) {
00996                         diagonal = *(diagonals2.begin());
00997                         diagonals2.pop_front();
00998                         if((diagonal.index2 - diagonal.index1) <= 1) continue;
00999                         ijreal = true;
01000                         jkreal = true;
01001                         pairs = &(dpstates[diagonal.index1][diagonal.index2].pairs);
01002                         if(!vertices[diagonal.index1].isConvex) {
01003                                 iter = pairs->end();
01004                                 iter--;
01005                                 j = iter->index2;
01006                                 if(iter->index1 != iter->index2) ijreal = false;
01007                         } else {
01008                                 iter = pairs->begin();
01009                                 j = iter->index1;
01010                                 if(iter->index1 != iter->index2) jkreal = false;
01011                         }
01012 
01013                         newdiagonal.index1 = diagonal.index1;
01014                         newdiagonal.index2 = j;
01015                         if(ijreal) {
01016                                 diagonals.push_back(newdiagonal);
01017                         } else {
01018                                 diagonals2.push_back(newdiagonal);
01019                         }
01020 
01021                         newdiagonal.index1 = j;
01022                         newdiagonal.index2 = diagonal.index2;
01023                         if(jkreal) {
01024                                 diagonals.push_back(newdiagonal);
01025                         } else {
01026                                 diagonals2.push_back(newdiagonal);
01027                         }
01028 
01029                         indices.push_back(j);
01030                 }
01031 
01032                 indices.sort();
01033                 newpoly.Init((long)indices.size());
01034                 k=0;
01035                 for(iiter = indices.begin();iiter!=indices.end();iiter++) {
01036                         newpoly[k] = vertices[*iiter].p;
01037                         k++;
01038                 }
01039                 parts->push_back(newpoly);
01040         }
01041 
01042         for(i=0;i<n;i++) {
01043                 delete [] dpstates[i];
01044         }
01045         delete [] dpstates;
01046         delete [] vertices;
01047 
01048         return ret;
01049 }
01050 
01051 //triangulates a set of polygons by first partitioning them into monotone polygons
01052 //O(n*log(n)) time complexity, O(n) space complexity
01053 //the algorithm used here is outlined in the book
01054 //"Computational Geometry: Algorithms and Applications" 
01055 //by Mark de Berg, Otfried Cheong, Marc van Kreveld and Mark Overmars
01056 int TPPLPartition::MonotonePartition(list<TPPLPoly> *inpolys, list<TPPLPoly> *monotonePolys) {
01057         list<TPPLPoly>::iterator iter;
01058         MonotoneVertex *vertices;
01059         long i,numvertices,vindex,vindex2,newnumvertices,maxnumvertices;
01060         long polystartindex, polyendindex;
01061         TPPLPoly *poly;
01062         MonotoneVertex *v,*v2,*vprev,*vnext;
01063         ScanLineEdge newedge;
01064         bool error = false;
01065 
01066         numvertices = 0;
01067         for(iter = inpolys->begin(); iter != inpolys->end(); iter++) {
01068                 numvertices += iter->GetNumPoints();
01069         }
01070 
01071         maxnumvertices = numvertices*3;
01072         vertices = new MonotoneVertex[maxnumvertices];
01073         newnumvertices = numvertices;
01074 
01075         polystartindex = 0;
01076         for(iter = inpolys->begin(); iter != inpolys->end(); iter++) {
01077                 poly = &(*iter);
01078                 polyendindex = polystartindex + poly->GetNumPoints()-1;
01079                 for(i=0;i<poly->GetNumPoints();i++) {
01080                         vertices[i+polystartindex].p = poly->GetPoint(i);
01081                         if(i==0) vertices[i+polystartindex].previous = polyendindex;
01082                         else vertices[i+polystartindex].previous = i+polystartindex-1;
01083                         if(i==(poly->GetNumPoints()-1)) vertices[i+polystartindex].next = polystartindex;
01084                         else vertices[i+polystartindex].next = i+polystartindex+1;
01085                 }
01086                 polystartindex = polyendindex+1;
01087         }
01088 
01089         //construct the priority queue
01090         long *priority = new long [numvertices];
01091         for(i=0;i<numvertices;i++) priority[i] = i;
01092         std::sort(priority,&(priority[numvertices]),VertexSorter(vertices));
01093 
01094         //determine vertex types
01095         char *vertextypes = new char[maxnumvertices];
01096         for(i=0;i<numvertices;i++) {
01097                 v = &(vertices[i]);
01098                 vprev = &(vertices[v->previous]);
01099                 vnext = &(vertices[v->next]);
01100 
01101                 if(Below(vprev->p,v->p)&&Below(vnext->p,v->p)) {
01102                         if(IsConvex(vnext->p,vprev->p,v->p)) {
01103                                 vertextypes[i] = TPPL_VERTEXTYPE_START;
01104                         } else {
01105                                 vertextypes[i] = TPPL_VERTEXTYPE_SPLIT;
01106                         }
01107                 } else if(Below(v->p,vprev->p)&&Below(v->p,vnext->p)) {
01108                         if(IsConvex(vnext->p,vprev->p,v->p))
01109                         {
01110                                 vertextypes[i] = TPPL_VERTEXTYPE_END;
01111                         } else {
01112                                 vertextypes[i] = TPPL_VERTEXTYPE_MERGE;
01113                         }
01114                 } else {
01115                         vertextypes[i] = TPPL_VERTEXTYPE_REGULAR;
01116                 }
01117         }
01118 
01119         //helpers
01120         long *helpers = new long[maxnumvertices];
01121 
01122         //binary search tree that holds edges intersecting the scanline
01123         //note that while set doesn't actually have to be implemented as a tree
01124         //complexity requirements for operations are the same as for the balanced binary search tree
01125         set<ScanLineEdge> edgeTree;
01126         //store iterators to the edge tree elements
01127         //this makes deleting existing edges much faster
01128         set<ScanLineEdge>::iterator *edgeTreeIterators,edgeIter;
01129         edgeTreeIterators = new set<ScanLineEdge>::iterator[maxnumvertices];
01130         pair<set<ScanLineEdge>::iterator,bool> edgeTreeRet;
01131 
01132         //for each vertex
01133         for(i=0;i<numvertices;i++) {
01134                 vindex = priority[i];
01135                 v = &(vertices[vindex]);
01136                 vindex2 = vindex;
01137                 v2 = v;
01138 
01139                 //depending on the vertex type, do the appropriate action
01140                 //comments in the following sections are copied from "Computational Geometry: Algorithms and Applications"
01141                 switch(vertextypes[vindex]) {
01142                         case TPPL_VERTEXTYPE_START:
01143                                 //Insert ei in T and set helper(ei) to vi.
01144                                 newedge.p1 = v->p;
01145                                 newedge.p2 = vertices[v->next].p;
01146                                 newedge.index = vindex;
01147                                 edgeTreeRet = edgeTree.insert(newedge);
01148                                 edgeTreeIterators[vindex] = edgeTreeRet.first;
01149                                 helpers[vindex] = vindex;
01150                                 break;
01151 
01152                         case TPPL_VERTEXTYPE_END:
01153                                 //if helper(ei-1) is a merge vertex
01154                                 if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
01155                                         //Insert the diagonal connecting vi to helper(ei-1) in D.
01156                                         AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous]);
01157                                         vertextypes[newnumvertices-2] = vertextypes[vindex];
01158                                         edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex];
01159                                         helpers[newnumvertices-2] = helpers[vindex];
01160                                         vertextypes[newnumvertices-1] = vertextypes[helpers[v->previous]];
01161                                         edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[v->previous]];
01162                                         helpers[newnumvertices-1] = helpers[helpers[v->previous]];
01163                                 }
01164                                 //Delete ei-1 from T
01165                                 edgeTree.erase(edgeTreeIterators[v->previous]);
01166                                 break;
01167 
01168                         case TPPL_VERTEXTYPE_SPLIT:
01169                                 //Search in T to find the edge e j directly left of vi.
01170                                 newedge.p1 = v->p;
01171                                 newedge.p2 = v->p;
01172                                 edgeIter = edgeTree.lower_bound(newedge);
01173                                 if(edgeIter == edgeTree.begin()) {
01174                                         error = true;
01175                                         break;
01176                                 }
01177                                 edgeIter--;
01178                                 //Insert the diagonal connecting vi to helper(ej) in D.
01179                                 AddDiagonal(vertices,&newnumvertices,vindex,helpers[edgeIter->index]);
01180                                 vertextypes[newnumvertices-2] = vertextypes[vindex];
01181                                 edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex];
01182                                 helpers[newnumvertices-2] = helpers[vindex];
01183                                 vertextypes[newnumvertices-1] = vertextypes[helpers[edgeIter->index]];
01184                                 edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[edgeIter->index]];
01185                                 helpers[newnumvertices-1] = helpers[helpers[edgeIter->index]];
01186                                 vindex2 = newnumvertices-2;
01187                                 v2 = &(vertices[vindex2]);
01188                                 //helper(e j)�vi
01189                                 helpers[edgeIter->index] = vindex;
01190                                 //Insert ei in T and set helper(ei) to vi.
01191                                 newedge.p1 = v2->p;
01192                                 newedge.p2 = vertices[v2->next].p;
01193                                 newedge.index = vindex2;
01194                                 edgeTreeRet = edgeTree.insert(newedge);
01195                                 edgeTreeIterators[vindex2] = edgeTreeRet.first;
01196                                 helpers[vindex2] = vindex2;
01197                                 break;
01198 
01199                         case TPPL_VERTEXTYPE_MERGE:
01200                                 //if helper(ei-1) is a merge vertex
01201                                 if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
01202                                         //Insert the diagonal connecting vi to helper(ei-1) in D.
01203                                         AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous]);
01204                                         vertextypes[newnumvertices-2] = vertextypes[vindex];
01205                                         edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex];
01206                                         helpers[newnumvertices-2] = helpers[vindex];
01207                                         vertextypes[newnumvertices-1] = vertextypes[helpers[v->previous]];
01208                                         edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[v->previous]];
01209                                         helpers[newnumvertices-1] = helpers[helpers[v->previous]];
01210                                         vindex2 = newnumvertices-2;
01211                                         v2 = &(vertices[vindex2]);
01212                                 }
01213                                 //Delete ei-1 from T.
01214                                 edgeTree.erase(edgeTreeIterators[v->previous]);
01215                                 //Search in T to find the edge e j directly left of vi.
01216                                 newedge.p1 = v->p;
01217                                 newedge.p2 = v->p;
01218                                 edgeIter = edgeTree.lower_bound(newedge);
01219                                 if(edgeIter == edgeTree.begin()) {
01220                                         error = true;
01221                                         break;
01222                                 }
01223                                 edgeIter--;
01224                                 //if helper(ej) is a merge vertex
01225                                 if(vertextypes[helpers[edgeIter->index]]==TPPL_VERTEXTYPE_MERGE) {
01226                                         //Insert the diagonal connecting vi to helper(e j) in D.
01227                                         AddDiagonal(vertices,&newnumvertices,vindex2,helpers[edgeIter->index]);
01228                                         vertextypes[newnumvertices-2] = vertextypes[vindex2];
01229                                         edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex2];
01230                                         helpers[newnumvertices-2] = helpers[vindex2];
01231                                         vertextypes[newnumvertices-1] = vertextypes[helpers[edgeIter->index]];
01232                                         edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[edgeIter->index]];
01233                                         helpers[newnumvertices-1] = helpers[helpers[edgeIter->index]];
01234                                 }
01235                                 //helper(e j)�vi
01236                                 helpers[edgeIter->index] = vindex2;
01237                                 break;
01238 
01239                         case TPPL_VERTEXTYPE_REGULAR:
01240                                 //if the interior of P lies to the right of vi
01241                                 if(Below(v->p,vertices[v->previous].p)) {
01242                                         //if helper(ei-1) is a merge vertex
01243                                         if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
01244                                                 //Insert the diagonal connecting vi to helper(ei-1) in D.
01245                                                 AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous]);
01246                                                 vertextypes[newnumvertices-2] = vertextypes[vindex];
01247                                                 edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex];
01248                                                 helpers[newnumvertices-2] = helpers[vindex];
01249                                                 vertextypes[newnumvertices-1] = vertextypes[helpers[v->previous]];
01250                                                 edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[v->previous]];
01251                                                 helpers[newnumvertices-1] = helpers[helpers[v->previous]];
01252                                                 vindex2 = newnumvertices-2;
01253                                                 v2 = &(vertices[vindex2]);
01254                                         }
01255                                         //Delete ei-1 from T.
01256                                         edgeTree.erase(edgeTreeIterators[v->previous]);
01257                                         //Insert ei in T and set helper(ei) to vi.
01258                                         newedge.p1 = v2->p;
01259                                         newedge.p2 = vertices[v2->next].p;
01260                                         newedge.index = vindex2;
01261                                         edgeTreeRet = edgeTree.insert(newedge);
01262                                         edgeTreeIterators[vindex2] = edgeTreeRet.first;
01263                                         helpers[vindex2] = vindex;
01264                                 } else {
01265                                         //Search in T to find the edge ej directly left of vi.
01266                                         newedge.p1 = v->p;
01267                                         newedge.p2 = v->p;
01268                                         edgeIter = edgeTree.lower_bound(newedge);
01269                                         if(edgeIter == edgeTree.begin()) {
01270                                                 error = true;
01271                                                 break;
01272                                         }
01273                                         edgeIter--;
01274                                         //if helper(ej) is a merge vertex
01275                                         if(vertextypes[helpers[edgeIter->index]]==TPPL_VERTEXTYPE_MERGE) {
01276                                                 //Insert the diagonal connecting vi to helper(e j) in D.
01277                                                 AddDiagonal(vertices,&newnumvertices,vindex,helpers[edgeIter->index]);
01278                                                 vertextypes[newnumvertices-2] = vertextypes[vindex];
01279                                                 edgeTreeIterators[newnumvertices-2] = edgeTreeIterators[vindex];
01280                                                 helpers[newnumvertices-2] = helpers[vindex];
01281                                                 vertextypes[newnumvertices-1] = vertextypes[helpers[edgeIter->index]];
01282                                                 edgeTreeIterators[newnumvertices-1] = edgeTreeIterators[helpers[edgeIter->index]];
01283                                                 helpers[newnumvertices-1] = helpers[helpers[edgeIter->index]];
01284                                         }
01285                                         //helper(e j)�vi
01286                                         helpers[edgeIter->index] = vindex;
01287                                 }
01288                                 break;
01289                 }
01290 
01291                 if(error) break;
01292         }
01293 
01294         char *used = new char[newnumvertices];
01295         memset(used,0,newnumvertices*sizeof(char));
01296 
01297         if(!error) {
01298                 //return result
01299                 long size;
01300                 TPPLPoly mpoly;
01301                 for(i=0;i<newnumvertices;i++) {
01302                         if(used[i]) continue;
01303                         v = &(vertices[i]);
01304                         vnext = &(vertices[v->next]);
01305                         size = 1;
01306                         while(vnext!=v) {
01307                                 vnext = &(vertices[vnext->next]);
01308                                 size++;
01309                         }
01310                         mpoly.Init(size);
01311                         v = &(vertices[i]);
01312                         mpoly[0] = v->p;
01313                         vnext = &(vertices[v->next]);
01314                         size = 1;
01315                         used[i] = 1;
01316                         used[v->next] = 1;
01317                         while(vnext!=v) {
01318                                 mpoly[size] = vnext->p;
01319                                 used[vnext->next] = 1;
01320                                 vnext = &(vertices[vnext->next]);
01321                                 size++;
01322                         }
01323                         monotonePolys->push_back(mpoly);
01324                 }
01325         }
01326 
01327         //cleanup
01328         delete [] vertices;
01329         delete [] priority;
01330         delete [] vertextypes;
01331         delete [] edgeTreeIterators;
01332         delete [] helpers;
01333         delete [] used;
01334 
01335         if(error) {
01336                 return 0;
01337         } else {
01338                 return 1;
01339         }
01340 }
01341 
01342 //adds a diagonal to the doubly-connected list of vertices
01343 void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, long index1, long index2) {
01344         long newindex1,newindex2;
01345 
01346         newindex1 = *numvertices;
01347         (*numvertices)++;
01348         newindex2 = *numvertices;
01349         (*numvertices)++;
01350 
01351         vertices[newindex1].p = vertices[index1].p;
01352         vertices[newindex2].p = vertices[index2].p;
01353 
01354         vertices[newindex2].next = vertices[index2].next;
01355         vertices[newindex1].next = vertices[index1].next;
01356 
01357         vertices[vertices[index2].next].previous = newindex2;
01358         vertices[vertices[index1].next].previous = newindex1;
01359 
01360         vertices[index1].next = newindex2;
01361         vertices[newindex2].previous = index1;
01362 
01363         vertices[index2].next = newindex1;
01364         vertices[newindex1].previous = index2;
01365 }
01366 
01367 bool TPPLPartition::Below(TPPLPoint &p1, TPPLPoint &p2) {
01368         if(p1.y < p2.y) return true;
01369         else if(p1.y == p2.y) {
01370                 if(p1.x < p2.x) return true;
01371         }
01372         return false;
01373 }
01374 
01375 //sorts in the falling order of y values, if y is equal, x is used instead
01376 bool TPPLPartition::VertexSorter::operator() (long index1, long index2) {
01377         if(vertices[index1].p.y > vertices[index2].p.y) return true;
01378         else if(vertices[index1].p.y == vertices[index2].p.y) {
01379                 if(vertices[index1].p.x > vertices[index2].p.x) return true;
01380         }
01381         return false;
01382 }
01383 
01384 bool TPPLPartition::ScanLineEdge::IsConvex(const TPPLPoint& p1, const TPPLPoint& p2, const TPPLPoint& p3) const {
01385         tppl_float tmp;
01386         tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
01387         if(tmp>0) return 1;
01388         else return 0;
01389 }
01390 
01391 bool TPPLPartition::ScanLineEdge::operator < (const ScanLineEdge & other) const {
01392         if(other.p1.y == other.p2.y) {
01393                 if(p1.y == p2.y) {
01394                         if(p1.y < other.p1.y) return true;
01395                         else return false;
01396                 }
01397                 if(IsConvex(p1,p2,other.p1)) return true;
01398                 else return false;
01399         } else if(p1.y == p2.y) {
01400                 if(IsConvex(other.p1,other.p2,p1)) return false;
01401                 else return true;       
01402         } else if(p1.y < other.p1.y) {
01403                 if(IsConvex(other.p1,other.p2,p1)) return false;
01404                 else return true;                       
01405         } else {
01406                 if(IsConvex(p1,p2,other.p1)) return true;
01407                 else return false;
01408         }
01409 }
01410 
01411 //triangulates monotone polygon
01412 //O(n) time, O(n) space complexity
01413 int TPPLPartition::TriangulateMonotone(TPPLPoly *inPoly, list<TPPLPoly> *triangles) {
01414         long i,i2,j,topindex,bottomindex,leftindex,rightindex,vindex;
01415         TPPLPoint *points;
01416         long numpoints;
01417         TPPLPoly triangle;
01418 
01419         numpoints = inPoly->GetNumPoints();
01420         points = inPoly->GetPoints();
01421 
01422         //trivial calses
01423         if(numpoints < 3) return 0;
01424         if(numpoints == 3) {
01425                 triangles->push_back(*inPoly);
01426         }
01427 
01428         topindex = 0; bottomindex=0;
01429         for(i=1;i<numpoints;i++) {
01430                 if(Below(points[i],points[bottomindex])) bottomindex = i;
01431                 if(Below(points[topindex],points[i])) topindex = i;
01432         }
01433 
01434         //check if the poly is really monotone
01435         i = topindex;
01436         while(i!=bottomindex) {
01437                 i2 = i+1; if(i2>=numpoints) i2 = 0;
01438                 if(!Below(points[i2],points[i])) return 0;
01439                 i = i2;
01440         }
01441         i = bottomindex;
01442         while(i!=topindex) {
01443                 i2 = i+1; if(i2>=numpoints) i2 = 0;
01444                 if(!Below(points[i],points[i2])) return 0;
01445                 i = i2;
01446         }
01447 
01448         char *vertextypes = new char[numpoints];
01449         long *priority = new long[numpoints];
01450 
01451         //merge left and right vertex chains
01452         priority[0] = topindex;
01453         vertextypes[topindex] = 0;
01454         leftindex = topindex+1; if(leftindex>=numpoints) leftindex = 0;
01455         rightindex = topindex-1; if(rightindex<0) rightindex = numpoints-1;
01456         for(i=1;i<(numpoints-1);i++) {
01457                 if(leftindex==bottomindex) {
01458                         priority[i] = rightindex;
01459                         rightindex--; if(rightindex<0) rightindex = numpoints-1;
01460                         vertextypes[priority[i]] = -1;
01461                 } else if(rightindex==bottomindex) {
01462                         priority[i] = leftindex;
01463                         leftindex++;  if(leftindex>=numpoints) leftindex = 0;
01464                         vertextypes[priority[i]] = 1;
01465                 } else {
01466                         if(Below(points[leftindex],points[rightindex])) {
01467                                 priority[i] = rightindex;
01468                                 rightindex--; if(rightindex<0) rightindex = numpoints-1;
01469                                 vertextypes[priority[i]] = -1;
01470                         } else {
01471                                 priority[i] = leftindex;
01472                                 leftindex++;  if(leftindex>=numpoints) leftindex = 0;
01473                                 vertextypes[priority[i]] = 1;                   
01474                         }
01475                 }
01476         }
01477         priority[i] = bottomindex;
01478         vertextypes[bottomindex] = 0;
01479 
01480         long *stack = new long[numpoints];
01481         long stackptr = 0;
01482 
01483         stack[0] = priority[0];
01484         stack[1] = priority[1];
01485         stackptr = 2;
01486 
01487         //for each vertex from top to bottom trim as many triangles as possible
01488         for(i=2;i<(numpoints-1);i++) {
01489                 vindex = priority[i];
01490                 if(vertextypes[vindex]!=vertextypes[stack[stackptr-1]]) {
01491                         for(j=0;j<(stackptr-1);j++) {
01492                                 if(vertextypes[vindex]==1) {
01493                                         triangle.Triangle(points[stack[j+1]],points[stack[j]],points[vindex]);
01494                                 } else {
01495                                         triangle.Triangle(points[stack[j]],points[stack[j+1]],points[vindex]);
01496                                 }
01497                                 triangles->push_back(triangle);
01498                         }
01499                         stack[0] = priority[i-1];
01500                         stack[1] = priority[i];
01501                         stackptr = 2;
01502                 } else {
01503                         stackptr--;
01504                         while(stackptr>0) {
01505                                 if(vertextypes[vindex]==1) {
01506                                         if(IsConvex(points[vindex],points[stack[stackptr-1]],points[stack[stackptr]])) {
01507                                                 triangle.Triangle(points[vindex],points[stack[stackptr-1]],points[stack[stackptr]]);
01508                                                 triangles->push_back(triangle);
01509                                                 stackptr--;
01510                                         } else {
01511                                                 break;
01512                                         }
01513                                 } else {
01514                                         if(IsConvex(points[vindex],points[stack[stackptr]],points[stack[stackptr-1]])) {
01515                                                 triangle.Triangle(points[vindex],points[stack[stackptr]],points[stack[stackptr-1]]);
01516                                                 triangles->push_back(triangle);
01517                                                 stackptr--;
01518                                         } else {
01519                                                 break;
01520                                         }
01521                                 }
01522                         }
01523                         stackptr++;
01524                         stack[stackptr] = vindex;
01525                         stackptr++;
01526                 }
01527         }
01528         vindex = priority[i];
01529         for(j=0;j<(stackptr-1);j++) {
01530                 if(vertextypes[stack[j+1]]==1) {
01531                         triangle.Triangle(points[stack[j]],points[stack[j+1]],points[vindex]);
01532                 } else {
01533                         triangle.Triangle(points[stack[j+1]],points[stack[j]],points[vindex]);
01534                 }
01535                 triangles->push_back(triangle);
01536         }
01537 
01538         delete [] priority;
01539         delete [] vertextypes;
01540         delete [] stack;
01541 
01542         return 1;
01543 }
01544 
01545 int TPPLPartition::Triangulate_MONO(list<TPPLPoly> *inpolys, list<TPPLPoly> *triangles) {
01546         list<TPPLPoly> monotone;
01547         list<TPPLPoly>::iterator iter;
01548 
01549         if(!MonotonePartition(inpolys,&monotone)) return 0;
01550         for(iter = monotone.begin(); iter!=monotone.end();iter++) {
01551                 if(!TriangulateMonotone(&(*iter),triangles)) return 0;
01552         }
01553         return 1;
01554 }
01555 
01556 int TPPLPartition::Triangulate_MONO(TPPLPoly *poly, list<TPPLPoly> *triangles) {
01557         list<TPPLPoly> polys;
01558         polys.push_back(*poly);
01559 
01560         return Triangulate_MONO(&polys, triangles);
01561 }


polypartition
Author(s): Georg Arbeiter
autogenerated on Wed Aug 26 2015 11:02:14