experimental - use cgal for triangulation

try breaking up construct_bucket into stages for better tile matching

src/BuildTiles/Main/construct.cxx

@@ -52,11 +52,12 @@
 using std::string;
 
 //double gSnap = 0.00000001;      // approx 1 mm
-double gSnap = 0.0000001;      // approx 1 cm
+double gSnap = 0.0000002;      // approx 2 cm
 
 static const double cover_size = 1.0 / 120.0;
 static const double half_cover_size = cover_size * 0.5;
 
+static unsigned int cur_poly_id = 0;
 
 // If we don't offset land use squares by some amount, then we can get
 // land use square boundaries coinciding with tile boundaries.
@@ -242,7 +243,6 @@ bool TGConstruct::load_poly(const string& path) {
     int hole_flag;
     int num_polys;
     double startx, starty, startz, x, y, z, lastx, lasty, lastz;
-    static unsigned int cur_id = 0;
 
     sg_gzifstream in( path );
 
@@ -437,7 +437,7 @@ bool TGConstruct::load_poly(const string& path) {
         // Once the full poly is loaded, build the clip mask
         shape.BuildMask();
         shape.area = area;
-        shape.id = cur_id++;
+        shape.id = cur_poly_id++;
 
         polys_in.add_shape( area, shape );
 
@@ -1844,6 +1844,7 @@ void TGConstruct::TesselatePolys( void )
 
             for ( unsigned int segment = 0; segment < polys_clipped.shape_size(area, shape); segment++ ) {
                 TGPolygon poly = polys_clipped.get_poly(area, shape, segment);
+
                 poly.get_bounding_box(min, max);
                 poly_extra = nodes.get_geod_inside( min, max );
 
@@ -1851,6 +1852,7 @@ void TGConstruct::TesselatePolys( void )
                         shape+1 << "-" << segment << " of " << (int)polys_clipped.area_size(area) << 
                         ": id = " << id );
 
+//              TGPolygon tri = polygon_tesselate_alt_with_extra( poly, poly_extra, false );
                 TGPolygon tri = polygon_tesselate_alt_with_extra( poly, poly_extra, false );
 
                 // ensure all added nodes are accounted for
@@ -2092,7 +2094,7 @@ void TGConstruct::CalcTextureCoordinates( void )
 //        loading, clipping, tesselating, normals, and output
 //        Also, we are still calculating some thing more than one 
 //        (like face area - need to move this into superpoly )
-void TGConstruct::ConstructBucket() {
+void TGConstruct::ConstructBucketStage1() {
 
     /* If we have some debug IDs, create a datasource */
     if ( debug_shapes.size() || debug_all ) {
@@ -2127,11 +2129,14 @@ void TGConstruct::ConstructBucket() {
     // Clean the polys - after this, we shouldn't change their shape (other than slightly for
     // fix T-Junctions - as This is the end of the first pass for multicore design
     CleanClippedPolys();
-    
+
+    // END OF FIRST PASS : SAVE THE TILE DATA
+
     // STEP 5)
     // Merge in Shared data (just add the nodes to the nodelist)
     // When this step is complete, some nodes will have normals (from shared tiles)
     // and some will not
+    // Load Shared Edge Data X,Y coords only
     LoadSharedEdgeData();
 
     // STEP 6) 
@@ -2154,14 +2159,19 @@ void TGConstruct::ConstructBucket() {
     // NOTE: After this point, no new nodes can be added
     LookupNodesPerVertex();
 
-    // STEP 9) 
+    // STEP 9)
-    // Interpolate elevations, and flatten stuff 
+    // Interpolate elevations, and flatten stuff
     CalcElevations();
-    
+
-    // STEP 10) 
+    // STEP 10)
-    // Generate face_connected list - 
+    // Generate face_connected list - shared data contains faces, too - save them somehow
     LookupFacesPerNode();
 
+    // END OF SECOND PASS : SAVE THE TILE DATA
+
+    // load shared edge data (with elevations, and face connected list)
+    // LoadSharedEdgeDataWithElevation();
+
     // STEP 11)
     // Calculate Face Normals
     CalcFaceNormals();

src/BuildTiles/Main/construct.hxx

@@ -367,8 +367,9 @@ public:
 
     void set_bucket( SGBucket b ) { bucket = b; }
     
-    void ConstructBucket();
+    void ConstructBucketStage1();
-
+    void ConstructBucketStage2();
+    
 
     void calc_gc_course_dist( const Point3D& start, const Point3D& dest, 
                               double *course, double *dist );

src/BuildTiles/Main/main.cxx

@@ -329,7 +329,7 @@ int main(int argc, char **argv) {
 
             c->set_debug( debug_dir, debug_defs );
 
-            c->ConstructBucket();
+            c->ConstructBucketStage1();
             delete c;
         } else {
             // build all the tiles in an area
@@ -357,7 +357,7 @@ int main(int argc, char **argv) {
 
                 c->set_debug( debug_dir, debug_defs );
 
-                c->ConstructBucket();
+                c->ConstructBucketStage1();
                 delete c;
             } else {
                 SGBucket b_cur;
@@ -390,7 +390,7 @@ int main(int argc, char **argv) {
 
                             c->set_debug( debug_dir, debug_defs );
 
-                            c->ConstructBucket();
+                            c->ConstructBucketStage1();
                             delete c;
                         } else {
                             SG_LOG(SG_GENERAL, SG_ALERT, "skipping " << b_cur);
@@ -418,7 +418,7 @@ int main(int argc, char **argv) {
 
         c->set_debug( debug_dir, debug_defs );
 
-        c->ConstructBucket();
+        c->ConstructBucketStage1();
         delete c;
     }
 

src/Lib/Geometry/CMakeLists.txt

@@ -7,6 +7,7 @@ add_library(Geometry STATIC
     line.cxx
     poly_support.cxx
     poly_extra.cxx
+    poly_cgal.cxx
     rectangle.cxx
     tg_nodes.cxx
     trinodes.cxx

src/Lib/Geometry/poly_cgal.cxx

@@ -0,0 +1,151 @@
+#include <float.h>
+#include <stdio.h>
+
+#include <simgear/compiler.h>
+#include <simgear/constants.h>
+#include <Geometry/point3d.hxx>
+#include <simgear/math/sg_types.hxx>
+#include <simgear/debug/logstream.hxx>
+#include <simgear/structure/exception.hxx>
+
+#include <Polygon/polygon.hxx>
+#include "poly_support.hxx"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Triangulation_face_base_with_info_2.h>
+#include <CGAL/Polygon_2.h>
+#include <CGAL/Triangle_2.h>
+#include <iostream>
+
+/* determining if a face is within the reulting poly */
+struct FaceInfo2
+{
+  FaceInfo2() {}
+  int nesting_level;
+
+  bool in_domain(){
+    return nesting_level%2 == 1;
+  }
+};
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
+typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
+typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,K>    Fbb;
+typedef CGAL::Constrained_triangulation_face_base_2<K,Fbb>        Fb;
+typedef CGAL::Triangulation_data_structure_2<Vb,Fb>               TDS;
+typedef CGAL::Exact_predicates_tag                                Itag;
+typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag>  CDT;
+typedef CDT::Point                                                Point;
+typedef CGAL::Polygon_2<K>                                        Polygon_2;
+typedef CGAL::Triangle_2<K>                                       Triangle_2;
+
+void mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list<CDT::Edge>& border )
+{
+    if(start->info().nesting_level != -1) {
+        return;
+    }
+
+    std::list<CDT::Face_handle> queue;
+    queue.push_back(start);
+
+    while( !queue.empty() ){
+        CDT::Face_handle fh = queue.front();
+        queue.pop_front();
+        if(fh->info().nesting_level == -1) {
+            fh->info().nesting_level = index;
+            for(int i = 0; i < 3; i++) {
+                CDT::Edge e(fh,i);
+                CDT::Face_handle n = fh->neighbor(i);
+                if(n->info().nesting_level == -1) {
+                    if(ct.is_constrained(e)) border.push_back(e);
+                    else queue.push_back(n);
+                }
+            }
+        }
+    }
+}
+
+//explore set of facets connected with non constrained edges,
+//and attribute to each such set a nesting level.
+//We start from facets incident to the infinite vertex, with a nesting
+//level of 0. Then we recursively consider the non-explored facets incident
+//to constrained edges bounding the former set and increase the nesting level by 1.
+//Facets in the domain are those with an odd nesting level.
+void mark_domains(CDT& cdt)
+{
+    for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){
+        it->info().nesting_level = -1;
+    }
+
+    int index = 0;
+    std::list<CDT::Edge> border;
+    mark_domains(cdt, cdt.infinite_face(), index++, border);
+    while(! border.empty()) {
+        CDT::Edge e = border.front();
+        border.pop_front();
+        CDT::Face_handle n = e.first->neighbor(e.second);
+        if(n->info().nesting_level == -1) {
+            mark_domains(cdt, n, e.first->info().nesting_level+1, border);
+        }
+    }
+}
+
+void insert_polygon(CDT& cdt,const Polygon_2& polygon)
+{
+    if ( polygon.is_empty() ) return;
+
+    CDT::Vertex_handle v_prev=cdt.insert(*CGAL::cpp0x::prev(polygon.vertices_end()));
+    for (Polygon_2::Vertex_iterator vit=polygon.vertices_begin(); vit!=polygon.vertices_end();++vit) {
+        CDT::Vertex_handle vh=cdt.insert(*vit);
+        cdt.insert_constraint(vh,v_prev);
+        v_prev=vh;
+    }
+}
+
+TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const point_list& extra_nodes, bool verbose ) {
+    TGPolygon result;
+    CDT       cdt;
+
+    result.erase();
+
+    // Bail right away if polygon is empty
+    if ( p.contours() == 0 ) {
+        return result;
+    }
+
+    // insert each polygon into the triangulation
+    for (int c = 0; c < p.contours(); c++) {
+        point_list contour = p.get_contour( c );
+        Polygon_2  poly;
+
+        for (unsigned int n = 0; n < contour.size(); n++ ) {
+            Point3D node = contour[n];
+            poly.push_back( Point( node.x(), node.y()) );
+        }
+
+        insert_polygon(cdt, poly);
+    }
+
+    mark_domains( cdt );
+
+    int count=0;
+    for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end();++fit) {
+        if ( fit->info().in_domain() ) {
+            Triangle_2 tri = cdt.triangle(fit);
+
+            Point3D p0 = Point3D( tri.vertex(0).x(), tri.vertex(0).y(), 0.0f );
+            Point3D p1 = Point3D( tri.vertex(1).x(), tri.vertex(1).y(), 0.0f );
+            Point3D p2 = Point3D( tri.vertex(2).x(), tri.vertex(2).y(), 0.0f );
+
+            result.add_node( count, p0 );
+            result.add_node( count, p1 );
+            result.add_node( count, p2 );
+
+            ++count;
+            // create a contour in result with this face
+        }
+    }
+
+    return result;
+}

src/Lib/Geometry/poly_support.hxx

@@ -69,6 +69,9 @@ TGPolygon polygon_tesselate_alt( TGPolygon &p, bool verbose );
 TGPolygon polygon_tesselate_alt_with_extra( TGPolygon &p, 
             const point_list &extra_nodes, bool verbose );
 
+TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p,
+            const point_list &extra_nodes, bool verbose );
+
 // calculate some "arbitrary" point inside each of the polygons contours
 void calc_points_inside( TGPolygon& p );