From adfb4482e2c02090e79ccb4295cfa586294ded49 Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Sat, 1 Sep 2012 17:23:35 -0400 Subject: [PATCH] experimental - use cgal for triangulation try breaking up construct_bucket into stages for better tile matching --- src/BuildTiles/Main/construct.cxx | 30 ++++-- src/BuildTiles/Main/construct.hxx | 5 +- src/BuildTiles/Main/main.cxx | 8 +- src/Lib/Geometry/CMakeLists.txt | 1 + src/Lib/Geometry/poly_cgal.cxx | 151 ++++++++++++++++++++++++++++++ src/Lib/Geometry/poly_support.hxx | 3 + 6 files changed, 182 insertions(+), 16 deletions(-) create mode 100644 src/Lib/Geometry/poly_cgal.cxx diff --git a/src/BuildTiles/Main/construct.cxx b/src/BuildTiles/Main/construct.cxx index 3c2a5350..607bd6e0 100644 --- a/src/BuildTiles/Main/construct.cxx +++ b/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) - // Interpolate elevations, and flatten stuff + // STEP 9) + // Interpolate elevations, and flatten stuff CalcElevations(); - - // STEP 10) - // Generate face_connected list - + + // STEP 10) + // 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(); diff --git a/src/BuildTiles/Main/construct.hxx b/src/BuildTiles/Main/construct.hxx index 54b39cbe..dbf9cad6 100644 --- a/src/BuildTiles/Main/construct.hxx +++ b/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 ); diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 70eef44e..9f018614 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/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; } diff --git a/src/Lib/Geometry/CMakeLists.txt b/src/Lib/Geometry/CMakeLists.txt index e44ccb30..501438e1 100644 --- a/src/Lib/Geometry/CMakeLists.txt +++ b/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 diff --git a/src/Lib/Geometry/poly_cgal.cxx b/src/Lib/Geometry/poly_cgal.cxx new file mode 100644 index 00000000..c20b822b --- /dev/null +++ b/src/Lib/Geometry/poly_cgal.cxx @@ -0,0 +1,151 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include "poly_support.hxx" + +#include +#include +#include +#include +#include +#include + +/* 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 Vb; +typedef CGAL::Triangulation_face_base_with_info_2 Fbb; +typedef CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 TDS; +typedef CGAL::Exact_predicates_tag Itag; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +typedef CDT::Point Point; +typedef CGAL::Polygon_2 Polygon_2; +typedef CGAL::Triangle_2 Triangle_2; + +void mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list& border ) +{ + if(start->info().nesting_level != -1) { + return; + } + + std::list 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 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; +} \ No newline at end of file diff --git a/src/Lib/Geometry/poly_support.hxx b/src/Lib/Geometry/poly_support.hxx index 2cd9a52d..b2fb1715 100644 --- a/src/Lib/Geometry/poly_support.hxx +++ b/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 );