diff --git a/src/BuildTiles/Main/construct.cxx b/src/BuildTiles/Main/construct.cxx index 4b9c1520..2c95968e 100644 --- a/src/BuildTiles/Main/construct.cxx +++ b/src/BuildTiles/Main/construct.cxx @@ -53,6 +53,8 @@ double gSnap = 0.00000001; // approx 1 mm static const double cover_size = 1.0 / 120.0; static const double half_cover_size = cover_size * 0.5; +#define USE_CLIPPER (0) + // If we don't offset land use squares by some amount, then we can get // land use square boundaries coinciding with tile boundaries. // @@ -63,6 +65,15 @@ static const double half_cover_size = cover_size * 0.5; // problem. static const double quarter_cover_size = cover_size * 0.25; + +// For debug: +void* ds_id = NULL; // If we are going to build shapefiles +void* l_id = NULL; // datasource and layer IDs +char ds_name[128]; // need to clean this up and add cmdline options +char layer_name[128]; +char feature_name[128]; + + // Constructor TGConstruct::TGConstruct(): useUKGrid(false), @@ -80,8 +91,6 @@ bool TGConstruct::load_array() { string base = bucket.gen_base_path(); int i; - SG_LOG(SG_GENERAL, SG_ALERT, "load array"); - for ( i = 0; i < (int)load_dirs.size(); ++i ) { string array_path = get_work_base() + "/" + load_dirs[i] + "/" + base + "/" + bucket.gen_index_str(); @@ -92,9 +101,6 @@ bool TGConstruct::load_array() { } } - - SG_LOG(SG_GENERAL, SG_ALERT, "parse array"); - array.parse( bucket ); array.remove_voids( ); @@ -353,9 +359,7 @@ int TGConstruct::load_polys( ) { poly_path = get_work_base() + "/" + load_dirs[i] + '/' + base; SG_LOG(SG_GENERAL, SG_ALERT, "poly_path = " << poly_path); -// count += actual_load_polys( poly_path ); string tile_str = bucket.gen_index_str(); - simgear::Dir d(poly_path); if (!d.exists()) { SG_LOG(SG_GENERAL, SG_ALERT, "directory not found:" << poly_path); @@ -396,7 +400,11 @@ int TGConstruct::load_polys( ) { // to reduce the number of separate polygons. void TGConstruct::add_to_polys ( TGPolygon &accum, const TGPolygon &poly) { if ( accum.contours() > 0 ) { +#if USE_CLIPPER + accum = tgPolygonUnionClipper( accum, poly ); +#else accum = tgPolygonUnion( accum, poly ); +#endif } else { accum = poly; } @@ -449,7 +457,7 @@ double TGConstruct::measure_roughness( TGPolygon &poly ) { z = array.closest_nonvoid_elev( points[j].x() * 3600.0, points[j].y() * 3600.0 ); } - // cout << "elevation = " << z << endl; + if ( z < min_z ) { min_z = z; } @@ -579,7 +587,7 @@ void TGConstruct::add_intermediate_nodes( ) { after = current.total_size(); if (before != after) { - SG_LOG( SG_CLIPPER, SG_INFO, "Fixed t-juntions in " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() << " nodes increased from " << before << " to " << after ); + SG_LOG( SG_CLIPPER, SG_INFO, "Fixed t-juntions in " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() << " nodes increased from " << before << " to " << after ); } /* Save it back */ @@ -618,6 +626,11 @@ double TGConstruct::distanceSphere( const Point3D p1, const Point3D p2 ) { // hopefully, this will get better when we have the area lookup via superpoly... void TGConstruct::fix_point_heights() { + TGPolygon tri_poly; + double e1, e2, e3, min; + int n1, n2, n3; + Point3D p; + SG_LOG(SG_GENERAL, SG_ALERT, "fixing node heights"); for (int i = 0; i < (int)nodes.size(); ++i) { @@ -630,45 +643,50 @@ void TGConstruct::fix_point_heights() } } - // now flatten some stuuf + // now flatten some stuff for (int i = 0; i < TG_MAX_AREA_TYPES; i++) { if ( is_lake_area( (AreaType)i ) ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { - SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); + + tri_poly.erase(); + tri_poly = polys_clipped.superpolys[i][j].get_tris(); - TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); for (int k=0; k< tri_poly.contours(); k++) { if (tri_poly.contour_size(k) != 3) { SG_LOG(SG_GENERAL, SG_ALERT, "triangle doesnt have 3 nodes" << tri_poly.contour_size(k) ); exit(0); } - double e1, e2, e3, min; - int n1, n2, n3; + n1 = nodes.find( tri_poly.get_pt( k, 0 ) ); + e1 = nodes.get_node(n1).GetPosition().z(); - e1 = tri_poly.get_pt( k, 0 ).z(); - n1 = get_nodes()->find( tri_poly.get_pt( k, 0 ) ); + n2 = nodes.find( tri_poly.get_pt( k, 1 ) ); + e2 = nodes.get_node(n2).GetPosition().z(); - e2 = tri_poly.get_pt( k, 1 ).z(); - n2 = get_nodes()->find( tri_poly.get_pt( k, 1 ) ); - - e3 = tri_poly.get_pt( k, 2 ).z(); - n3 = get_nodes()->find( tri_poly.get_pt( k, 2 ) ); + n3 = nodes.find( tri_poly.get_pt( k, 2 ) ); + e3 = nodes.get_node(n3).GetPosition().z(); min = e1; if ( e2 < min ) { min = e2; } if ( e3 < min ) { min = e3; } - get_nodes()->SetElevation( n1, min ); - get_nodes()->SetElevation( n2, min ); - get_nodes()->SetElevation( n3, min ); + SG_LOG(SG_GENERAL, SG_ALERT, "FLATTEN LAKE: original elevations are " << + nodes.get_node(n1).GetPosition().z() << "(" << e1 << "), " << + nodes.get_node(n2).GetPosition().z() << "(" << e2 << "), " << + nodes.get_node(n3).GetPosition().z() << "(" << e3 << ")" << + " new elevation is " << min ); + + nodes.SetElevation( n1, min ); + nodes.SetElevation( n2, min ); + nodes.SetElevation( n3, min ); } } } if ( is_stream_area( (AreaType)i ) ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { - SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); for (int k=0; k< tri_poly.contours(); k++) { @@ -678,18 +696,15 @@ void TGConstruct::fix_point_heights() } point_list raw_nodes = nodes.get_geod_nodes(); - double e1, e2, e3, min; - int n1, n2, n3; - Point3D p; - e1 = tri_poly.get_pt( k, 0 ).z(); - n1 = get_nodes()->find( tri_poly.get_pt( k, 0 ) ); + n1 = nodes.find( tri_poly.get_pt( k, 0 ) ); + e1 = nodes.get_node(n1).GetPosition().z(); - e2 = tri_poly.get_pt( k, 1 ).z(); - n2 = get_nodes()->find( tri_poly.get_pt( k, 1 ) ); + n2 = nodes.find( tri_poly.get_pt( k, 1 ) ); + e2 = nodes.get_node(n2).GetPosition().z(); - e3 = tri_poly.get_pt( k, 2 ).z(); - n3 = get_nodes()->find( tri_poly.get_pt( k, 2 ) ); + n3 = nodes.find( tri_poly.get_pt( k, 2 ) ); + e3 = nodes.get_node(n3).GetPosition().z(); min = e1; p = raw_nodes[n1]; @@ -705,16 +720,16 @@ void TGConstruct::fix_point_heights() double max2 = d2 * 0.20 + min; double max3 = d3 * 0.20 + min; - if ( max1 < e1 ) { get_nodes()->SetElevation( n1, max1 ); } - if ( max2 < e2 ) { get_nodes()->SetElevation( n2, max2 ); } - if ( max3 < e3 ) { get_nodes()->SetElevation( n3, max3 ); } + if ( max1 < e1 ) { nodes.SetElevation( n1, max1 ); } + if ( max2 < e2 ) { nodes.SetElevation( n2, max2 ); } + if ( max3 < e3 ) { nodes.SetElevation( n3, max3 ); } } } } if ( is_road_area( (AreaType)i ) ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { - SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); for (int k=0; k< tri_poly.contours(); k++) { @@ -724,18 +739,15 @@ void TGConstruct::fix_point_heights() } point_list raw_nodes = nodes.get_geod_nodes(); - double e1, e2, e3, min; - int n1, n2, n3; - Point3D p; - e1 = tri_poly.get_pt( k, 0 ).z(); - n1 = get_nodes()->find( tri_poly.get_pt( k, 0 ) ); + n1 = nodes.find( tri_poly.get_pt( k, 0 ) ); + e1 = nodes.get_node(n1).GetPosition().z(); - e2 = tri_poly.get_pt( k, 1 ).z(); - n2 = get_nodes()->find( tri_poly.get_pt( k, 1 ) ); + n2 = nodes.find( tri_poly.get_pt( k, 1 ) ); + e2 = nodes.get_node(n2).GetPosition().z(); - e3 = tri_poly.get_pt( k, 2 ).z(); - n3 = get_nodes()->find( tri_poly.get_pt( k, 2 ) ); + n3 = nodes.find( tri_poly.get_pt( k, 2 ) ); + e3 = nodes.get_node(n3).GetPosition().z(); min = e1; p = raw_nodes[n1]; @@ -762,7 +774,7 @@ void TGConstruct::fix_point_heights() for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); - SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Flattening " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); for (int k=0; k< tri_poly.contours(); k++) { if (tri_poly.contour_size(k) != 3) { @@ -770,16 +782,13 @@ void TGConstruct::fix_point_heights() exit(0); } - int n1, n2, n3; + n1 = nodes.find( tri_poly.get_pt( k, 0 ) ); + n2 = nodes.find( tri_poly.get_pt( k, 1 ) ); + n3 = nodes.find( tri_poly.get_pt( k, 2 ) ); - n1 = get_nodes()->find( tri_poly.get_pt( k, 0 ) ); - n2 = get_nodes()->find( tri_poly.get_pt( k, 1 ) ); - n3 = get_nodes()->find( tri_poly.get_pt( k, 2 ) ); - - SG_LOG(SG_GENERAL, SG_ALERT, "Set Ocean Triangle " << n1 << "," << n2 << "," << n3 << " to 0.0" ); - get_nodes()->SetElevation( n1, 0.0 ); - get_nodes()->SetElevation( n2, 0.0 ); - get_nodes()->SetElevation( n3, 0.0 ); + nodes.SetElevation( n1, 0.0 ); + nodes.SetElevation( n2, 0.0 ); + nodes.SetElevation( n3, 0.0 ); } } } @@ -1027,7 +1036,11 @@ void TGConstruct::merge_slivers( TGPolyList& clipped, poly_list& slivers_list ) // cout << " polygon = " << j << endl; poly = clipped.superpolys[area][k].get_poly(); original_contours = poly.contours(); +#if USE_CLIPPER + result = tgPolygonUnionClipper( poly, sliver ); +#else result = tgPolygonUnion( poly, sliver ); +#endif result_contours = result.contours(); if ( original_contours == result_contours ) { @@ -1089,15 +1102,28 @@ bool TGConstruct::clip_all(const point2d& min, const point2d& max) { for ( i = 0; i < TG_MAX_AREA_TYPES; i++ ) { if ( is_landmass_area( i ) && !ignoreLandmass ) { for ( unsigned int j = 0; j < polys_in.superpolys[i].size(); ++j ) { +#if USE_CLIPPER + land_mask = tgPolygonUnionClipper( land_mask, polys_in.superpolys[i][j].get_poly() ); +#else land_mask = tgPolygonUnion( land_mask, polys_in.superpolys[i][j].get_poly() ); +#endif + } } else if ( is_water_area( i ) ) { for (unsigned int j = 0; j < polys_in.superpolys[i].size(); j++) { +#if USE_CLIPPER + water_mask = tgPolygonUnionClipper( water_mask, polys_in.superpolys[i][j].get_poly() ); +#else water_mask = tgPolygonUnion( water_mask, polys_in.superpolys[i][j].get_poly() ); +#endif } } else if ( is_island_area( i ) ) { for (unsigned int j = 0; j < polys_in.superpolys[i].size(); j++) { +#if USE_CLIPPER + island_mask = tgPolygonUnionClipper( island_mask, polys_in.superpolys[i][j].get_poly() ); +#else island_mask = tgPolygonUnion( island_mask, polys_in.superpolys[i][j].get_poly() ); +#endif } } } @@ -1108,33 +1134,55 @@ bool TGConstruct::clip_all(const point2d& min, const point2d& max) { for( j = 0; j < (int)polys_in.superpolys[i].size(); ++j ) { TGPolygon current = polys_in.superpolys[i][j].get_poly(); - SG_LOG( SG_CLIPPER, SG_INFO, "Clipping " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_in.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Clipping " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_in.superpolys[i].size() ); tmp = current; // if not a hole, clip the area to the land_mask if ( !ignoreLandmass && !is_hole_area( i ) ) { +#if USE_CLIPPER + tmp = tgPolygonIntClipper( tmp, land_mask ); +#else tmp = tgPolygonInt( tmp, land_mask ); +#endif } // if a water area, cut out potential islands if ( is_water_area( i ) ) { // clip against island mask +#if USE_CLIPPER + tmp = tgPolygonDiffClipper( tmp, island_mask ); +#else tmp = tgPolygonDiff( tmp, island_mask ); +#endif } - clipped = tgPolygonDiff( tmp, accum ); +#if 0 + if ( (i == 3) && (j == 137) ) { + SG_LOG( SG_CLIPPER, SG_INFO, "Error Poly\n" << tmp ); + sprintf( layer_name, "preclip" ); + l_id = tgShapefileOpenLayer( ds_id, layer_name ); + sprintf( feature_name, "preclip" ); + tgShapefileCreateFeature( ds_id, l_id, tmp, feature_name ); + } +#endif + +#if USE_CLIPPER + clipped = tgPolygonDiffClipper( tmp, accum ); +#else + clipped = tgPolygonDiff( tmp, accum ); +#endif // only add to output list if the clip left us with a polygon if ( clipped.contours() > 0 ) { // move slivers from clipped polygon to slivers polygon tgPolygonFindSlivers( clipped, slivers ); -// // merge any slivers with previously clipped -// // neighboring polygons -// if ( slivers.contours() > 0 ) { -// merge_slivers(polys_clipped, slivers); -// } + // merge any slivers with previously clipped + // neighboring polygons + //if ( slivers.contours() > 0 ) { + // merge_slivers(polys_clipped, slivers); + //} // add the sliverless result polygon to the clipped polys list if ( clipped.contours() > 0 ) { @@ -1147,7 +1195,11 @@ bool TGConstruct::clip_all(const point2d& min, const point2d& max) { } } +#if USE_CLIPPER + accum = tgPolygonUnionClipper( tmp, accum ); +#else accum = tgPolygonUnion( tmp, accum ); +#endif } } @@ -1160,7 +1212,11 @@ bool TGConstruct::clip_all(const point2d& min, const point2d& max) { // remains = new gpc_polygon; // remains->num_contours = 0; // remains->contour = NULL; +#if USE_CLIPPER + remains = tgPolygonDiffClipper( polys_in.safety_base, accum ); +#else remains = tgPolygonDiff( polys_in.safety_base, accum ); +#endif if ( remains.contours() > 0 ) { // cout << "remains contours = " << remains.contours() << endl; @@ -1192,8 +1248,56 @@ bool TGConstruct::clip_all(const point2d& min, const point2d& max) { return true; } +bool TGNodesSortByLon( const TGNode& n1, const TGNode& n2 ) +{ + return ( n1.GetPosition().x() < n2.GetPosition().x() ); +} + +static void dump_nodes( TGNodes* nodes ) { + for (unsigned int i=0; isize(); i++) { + TGNode node = nodes->get_node( i ); + string fixed; + + if ( node.GetFixedPosition() ) { + fixed = " z is fixed elevation "; + } else { + fixed = " z is interpolated elevation "; + } + + SG_LOG(SG_GENERAL, SG_ALERT, "Point[" << i << "] is " << node.GetPosition() << fixed ); + } +} + +static void dump_lat_nodes( TGConstruct& c, double lat ) { + node_list all_nodes = c.get_nodes()->get_node_list(); + node_list sorted_nodes; + for (unsigned int i=0; iunique_add_fixed_elevation(corner_list[i]); - get_nodes()->unique_add(corner_list[i]); + nodes.unique_add(corner_list[i]); } point_list fit_list = array.get_fitted_list(); for (unsigned int i=0; iunique_add_fixed_elevation(fit_list[i]); - get_nodes()->unique_add(fit_list[i]); + nodes.unique_add(fit_list[i]); } SG_LOG(SG_GENERAL, SG_ALERT, "NODE LIST AFTER FITTING" ); @@ -1254,7 +1356,7 @@ void TGConstruct::construct_bucket( SGBucket b ) { // When this step is complete, some nodes will have normals (from shared tiles) // and some will not - need to indicate this in the new node class - SG_LOG(SG_GENERAL, SG_ALERT, "number of geod nodes = before adding adding shared edges" << get_nodes()->size() ); + SG_LOG(SG_GENERAL, SG_ALERT, "number of geod nodes = before adding adding shared edges" << nodes.size() ); TGMatch m; m.load_neighbor_shared( bucket, work_base ); @@ -1270,7 +1372,7 @@ void TGConstruct::construct_bucket( SGBucket b ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { TGPolygon poly = polys_clipped.superpolys[i][j].get_poly(); - SG_LOG( SG_CLIPPER, SG_INFO, "Collecting nodes for " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Collecting nodes for " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); for (int k=0; k< poly.contours(); k++) { for (int l = 0; l < poly.contour_size(k); l++) { @@ -1297,7 +1399,20 @@ void TGConstruct::construct_bucket( SGBucket b ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { TGPolygon poly = polys_clipped.superpolys[i][j].get_poly(); - SG_LOG( SG_CLIPPER, SG_INFO, "Tesselating " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() << " : flag = " << polys_clipped.superpolys[i][j].get_flag()); + SG_LOG( SG_CLIPPER, SG_INFO, "Tesselating " << get_area_name( (AreaType)i ) << "(" << i << ") :" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() << " : flag = " << polys_clipped.superpolys[i][j].get_flag()); + + if ( (i == 3) && (j == 137) ) { + SG_LOG( SG_CLIPPER, SG_INFO, "Error Poly\n" << poly ); + + + sprintf( layer_name, "pretesselate" ); + l_id = tgShapefileOpenLayer( ds_id, layer_name ); + sprintf( feature_name, "pretesselate" ); + tgShapefileCreateFeature( ds_id, l_id, poly, feature_name ); + + // close befoe the crash + tgShapefileCloseDatasource( ds_id ); + } TGPolygon tri = polygon_tesselate_alt_with_extra( poly, extra, false ); TGPolygon tc; @@ -1316,8 +1431,7 @@ void TGConstruct::construct_bucket( SGBucket b ) { } } -#if 1 // I don't think this is necessary - triangulations set to don't add points... - // Add triangulation points + /* Add any points from triangulation */ for (int i = 0; i < TG_MAX_AREA_TYPES; i++) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); @@ -1329,11 +1443,16 @@ void TGConstruct::construct_bucket( SGBucket b ) { } } } -#endif + + SG_LOG(SG_GENERAL, SG_ALERT, "NODE LIST BEFORE FLATTEN" ); + dump_nodes( get_nodes() ); // Step 7) Flatten fix_point_heights(); + SG_LOG(SG_GENERAL, SG_ALERT, "NODE LIST AFTER FLATTEN" ); + dump_nodes( get_nodes() ); + #if 0 SG_LOG(SG_GENERAL, SG_ALERT, "REVERSE ELE LOOKUP "); @@ -1419,7 +1538,7 @@ void TGConstruct::construct_bucket( SGBucket b ) { // only tesselate non holes if ( !is_hole_area( i ) ) { for (int j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { - SG_LOG( SG_CLIPPER, SG_INFO, "Ouput nodes for " << get_area_name( (AreaType)i ) << ":" << j << " of " << (int)polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Ouput nodes for " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << (int)polys_clipped.superpolys[i].size() ); TGPolygon tri_poly = polys_clipped.superpolys[i][j].get_tris(); TGPolygon tri_txs = polys_clipped.superpolys[i][j].get_texcoords(); @@ -1433,7 +1552,7 @@ void TGConstruct::construct_bucket( SGBucket b ) { for (int l = 0; l < tri_poly.contour_size(k); ++l) { p = tri_poly.get_pt( k, l ); - index = get_nodes()->find( p ); + index = nodes.find( p ); if (index < 0) { SG_LOG(SG_GENERAL, SG_ALERT, "NODE NOT FOUND " << p); exit(0); @@ -1566,11 +1685,23 @@ void TGConstruct::clean_clipped_polys() { for ( i = 0; i < TG_MAX_AREA_TYPES; ++i ) { for( j = 0; j < (int)polys_clipped.superpolys[i].size(); ++j ) { TGPolygon poly = polys_clipped.superpolys[i][j].get_poly(); - SG_LOG( SG_CLIPPER, SG_INFO, "Cleaning poly " << get_area_name( (AreaType)i ) << ":" << j << " of " << polys_clipped.superpolys[i].size() ); + SG_LOG( SG_CLIPPER, SG_INFO, "Cleaning poly " << get_area_name( (AreaType)i ) << ":" << j+1 << " of " << polys_clipped.superpolys[i].size() ); - poly = snap(poly, gSnap); - poly = remove_dups( poly ); - poly = remove_bad_contours( poly ); + if ( i == 3 ) { + poly = remove_cycles( poly ); + poly = remove_dups( poly ); + poly = remove_bad_contours( poly ); + poly = tgPolygonSimplify( poly ); + poly = remove_tiny_contours( poly ); + poly = remove_spikes( poly ); + poly = remove_dups( poly ); + poly = remove_bad_contours( poly ); + poly = remove_tiny_contours( poly ); + } else { + poly = snap(poly, gSnap); + poly = remove_dups( poly ); + poly = remove_bad_contours( poly ); + } polys_clipped.superpolys[i][j].set_poly( poly ); } diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 722b7bb2..de188b76 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -326,50 +326,6 @@ static void do_output( TGConstruct& c, TGGenOutput& output ) { } #endif -bool TGNodesSortByLon( const TGNode& n1, const TGNode& n2 ) -{ - return ( n1.GetPosition().x() < n2.GetPosition().x() ); -} - -static void dump_nodes( TGConstruct& c ) { - for (unsigned int i=0; isize(); i++) { - TGNode node = c.get_nodes()->get_node( i ); - string fixed; - - if ( node.GetFixedPosition() ) { - fixed = " z is fixed elevation "; - } else { - fixed = " z is interpolated elevation "; - } - - SG_LOG(SG_GENERAL, SG_ALERT, "Point[" << i << "] is " << node.GetPosition() << fixed ); - } -} - -static void dump_lat_nodes( TGConstruct& c, double lat ) { - node_list all_nodes = c.get_nodes()->get_node_list(); - node_list sorted_nodes; - for (unsigned int i=0; i (const Int128 &val) const { - if (hi > val.hi) return true; - else if (hi < val.hi) return false; - else return ulong64(lo) > ulong64(val.lo); + if (hi != val.hi) + return hi > val.hi; + else + return lo > val.lo; } bool operator < (const Int128 &val) const { - if (hi < val.hi) return true; - else if (hi > val.hi) return false; - else return ulong64(lo) < ulong64(val.lo); + if (hi != val.hi) + return hi < val.hi; + else + return lo < val.lo; } Int128& operator += (const Int128 &rhs) @@ -140,14 +137,28 @@ class Int128 return *this; } + //Int128 operator -() const + //{ + // Int128 result(*this); + // if (result.lo == 0) { + // if (result.hi != 0) result.hi = -1; + // } + // else { + // result.lo = -result.lo; + // result.hi = ~result.hi; + // } + // return result; + //} + Int128 operator - (const Int128 &rhs) const { Int128 result(*this); - result-= rhs; + result -= rhs; return result; } - Int128 operator * (const Int128 &rhs) const { + Int128 operator * (const Int128 &rhs) const + { if ( !(hi == 0 || hi == -1) || !(rhs.hi == 0 || rhs.hi == -1)) throw "Int128 operator*: overflow error"; bool negate = (hi < 0) != (rhs.hi < 0); @@ -225,25 +236,25 @@ class Int128 } //for bug testing ... - std::string AsString() const - { - std::string result; - unsigned char r = 0; - Int128 tmp(0), val(*this); - if (hi < 0) Negate(val); - result.resize(50); - std::string::size_type i = result.size() -1; - while (val.hi != 0 || val.lo != 0) - { - Div10(val, tmp, r); - result[i--] = char('0' + r); - val = tmp; - } - if (hi < 0) result[i--] = '-'; - result.erase(0,i+1); - if (result.size() == 0) result = "0"; - return result; - } + //std::string AsString() const + //{ + // std::string result; + // unsigned char r = 0; + // Int128 tmp(0), val(*this); + // if (hi < 0) Negate(val); + // result.resize(50); + // std::string::size_type i = result.size() -1; + // while (val.hi != 0 || val.lo != 0) + // { + // Div10(val, tmp, r); + // result[i--] = char('0' + r); + // val = tmp; + // } + // if (hi < 0) result[i--] = '-'; + // result.erase(0,i+1); + // if (result.size() == 0) result = "0"; + // return result; + //} private: long64 hi; @@ -251,79 +262,70 @@ private: static void Negate(Int128 &val) { - if (val.lo == 0) - { - if( val.hi == 0) return; - val.lo = ~val.lo; - val.hi = ~val.hi +1; + if (val.lo == 0) { + if (val.hi != 0) val.hi = -val.hi;; } - else - { - val.lo = ~val.lo +1; + else { + val.lo = -val.lo; val.hi = ~val.hi; } } //debugging only ... - void Div10(const Int128 val, Int128& result, unsigned char & remainder) const - { - remainder = 0; - result = 0; - for (int i = 63; i >= 0; --i) - { - if ((val.hi & ((long64)1 << i)) != 0) - remainder = char((remainder * 2) + 1); else - remainder *= char(2); - if (remainder >= 10) - { - result.hi += ((long64)1 << i); - remainder -= char(10); - } - } - for (int i = 63; i >= 0; --i) - { - if ((val.lo & ((long64)1 << i)) != 0) - remainder = char((remainder * 2) + 1); else - remainder *= char(2); - if (remainder >= 10) - { - result.lo += ((long64)1 << i); - remainder -= char(10); - } - } - } + //void Div10(const Int128 val, Int128& result, unsigned char & remainder) const + //{ + // remainder = 0; + // result = 0; + // for (int i = 63; i >= 0; --i) + // { + // if ((val.hi & ((long64)1 << i)) != 0) + // remainder = char((remainder * 2) + 1); else + // remainder *= char(2); + // if (remainder >= 10) + // { + // result.hi += ((long64)1 << i); + // remainder -= char(10); + // } + // } + // for (int i = 63; i >= 0; --i) + // { + // if ((val.lo & ((long64)1 << i)) != 0) + // remainder = char((remainder * 2) + 1); else + // remainder *= char(2); + // if (remainder >= 10) + // { + // result.lo += ((long64)1 << i); + // remainder -= char(10); + // } + // } + //} }; //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ -RangeTest TestRange(const Polygon &pts) +bool FullRangeNeeded(const Polygon &pts) { - RangeTest result = rtLo; + bool result = false; for (Polygon::size_type i = 0; i < pts.size(); ++i) { if (Abs(pts[i].X) > hiRange || Abs(pts[i].Y) > hiRange) - return rtError; + throw "Coordinate exceeds range bounds."; else if (Abs(pts[i].X) > loRange || Abs(pts[i].Y) > loRange) - result = rtHi; + result = true; } return result; } //------------------------------------------------------------------------------ - + bool Orientation(const Polygon &poly) { int highI = (int)poly.size() -1; if (highI < 2) return false; - bool UseFullInt64Range = false; int j = 0, jplus, jminus; for (int i = 0; i <= highI; ++i) { - if (Abs(poly[i].X) > hiRange || Abs(poly[i].Y) > hiRange) - throw "Coordinate exceeds range bounds."; - if (Abs(poly[i].X) > loRange || Abs(poly[i].Y) > loRange) - UseFullInt64Range = true; if (poly[i].Y < poly[j].Y) continue; if ((poly[i].Y > poly[j].Y || poly[i].X < poly[j].X)) j = i; }; @@ -339,16 +341,18 @@ bool Orientation(const Polygon &poly) vec2.X = poly[jplus].X - poly[j].X; vec2.Y = poly[jplus].Y - poly[j].Y; - if (UseFullInt64Range) + if (Abs(vec1.X) > loRange || Abs(vec1.Y) > loRange || + Abs(vec2.X) > loRange || Abs(vec2.Y) > loRange) { + if (Abs(vec1.X) > hiRange || Abs(vec1.Y) > hiRange || + Abs(vec2.X) > hiRange || Abs(vec2.Y) > hiRange) + throw "Coordinate exceeds range bounds."; Int128 cross = Int128(vec1.X) * Int128(vec2.Y) - Int128(vec2.X) * Int128(vec1.Y); return cross > 0; } else - { return (vec1.X * vec2.Y - vec2.X * vec1.Y) > 0; - } } //------------------------------------------------------------------------------ @@ -378,7 +382,7 @@ bool Orientation(OutRec *outRec, bool UseFullInt64Range) //find vertices either side of bottomPt (skipping duplicate points) .... OutPt *opPrev = op->prev; OutPt *opNext = op->next; - while (op != opPrev && PointsEqual(op->pt, opPrev->pt)) + while (op != opPrev && PointsEqual(op->pt, opPrev->pt)) opPrev = opPrev->prev; while (op != opNext && PointsEqual(op->pt, opNext->pt)) opNext = opNext->next; @@ -400,21 +404,9 @@ double Area(const Polygon &poly) { int highI = (int)poly.size() -1; if (highI < 2) return 0; - bool UseFullInt64Range; - RangeTest rt = TestRange(poly); - switch (rt) { - case rtLo: - UseFullInt64Range = false; - break; - case rtHi: - UseFullInt64Range = true; - break; - default: - throw "Coordinate exceeds range bounds."; - } - if (UseFullInt64Range) { - Int128 a(0); + if (FullRangeNeeded(poly)) { + Int128 a; a = (Int128(poly[highI].X) * Int128(poly[0].Y)) - Int128(poly[0].X) * Int128(poly[highI].Y); for (int i = 0; i < highI; ++i) @@ -433,6 +425,30 @@ double Area(const Polygon &poly) } //------------------------------------------------------------------------------ +double Area(const OutRec &outRec, bool UseFullInt64Range) +{ + OutPt *op = outRec.pts; + if (UseFullInt64Range) { + Int128 a(0); + do { + a += (Int128(op->prev->pt.X) * Int128(op->pt.Y)) - + Int128(op->pt.X) * Int128(op->prev->pt.Y); + op = op->next; + } while (op != outRec.pts); + return a.AsDouble() / 2; + } + else + { + double a = 0; + do { + a += (op->prev->pt.X * op->pt.Y) - (op->pt.X * op->prev->pt.Y); + op = op->next; + } while (op != outRec.pts); + return a/2; + } +} +//------------------------------------------------------------------------------ + bool PointIsVertex(const IntPoint &pt, OutPt *pp) { OutPt *pp2 = pp; @@ -510,17 +526,15 @@ bool SlopesEqual(const IntPoint pt1, const IntPoint pt2, double GetDx(const IntPoint pt1, const IntPoint pt2) { - if (pt1.Y == pt2.Y) return HORIZONTAL; - else return - (double)(pt2.X - pt1.X) / (double)(pt2.Y - pt1.Y); + return (pt1.Y == pt2.Y) ? + HORIZONTAL : (double)(pt2.X - pt1.X) / (double)(pt2.Y - pt1.Y); } //--------------------------------------------------------------------------- void SetDx(TEdge &e) { if (e.ybot == e.ytop) e.dx = HORIZONTAL; - else e.dx = - (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot); + else e.dx = (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot); } //--------------------------------------------------------------------------- @@ -542,15 +556,15 @@ void SwapPolyIndexes(TEdge &edge1, TEdge &edge2) inline long64 Round(double val) { - if ((val < 0)) return static_cast(val - 0.5); - else return static_cast(val + 0.5); + return (val < 0) ? + static_cast(val - 0.5) : static_cast(val + 0.5); } //------------------------------------------------------------------------------ long64 TopX(TEdge &edge, const long64 currentY) { - if( currentY == edge.ytop ) return edge.xtop; - return edge.xbot + Round(edge.dx *(currentY - edge.ybot)); + return ( currentY == edge.ytop ) ? + edge.xtop : edge.xbot + Round(edge.dx *(currentY - edge.ybot)); } //------------------------------------------------------------------------------ @@ -710,17 +724,60 @@ bool GetOverlapSegment(IntPoint pt1a, IntPoint pt1b, IntPoint pt2a, } //------------------------------------------------------------------------------ -OutPt* PolygonBottom(OutPt* pp) +bool FirstIsBottomPt(const OutPt* btmPt1, const OutPt* btmPt2) { + OutPt *p = btmPt1->prev; + while (PointsEqual(p->pt, btmPt1->pt) && (p != btmPt1)) p = p->prev; + double dx1p = std::fabs(GetDx(btmPt1->pt, p->pt)); + p = btmPt1->next; + while (PointsEqual(p->pt, btmPt1->pt) && (p != btmPt1)) p = p->next; + double dx1n = std::fabs(GetDx(btmPt1->pt, p->pt)); + + p = btmPt2->prev; + while (PointsEqual(p->pt, btmPt2->pt) && (p != btmPt2)) p = p->prev; + double dx2p = std::fabs(GetDx(btmPt2->pt, p->pt)); + p = btmPt2->next; + while (PointsEqual(p->pt, btmPt2->pt) && (p != btmPt2)) p = p->next; + double dx2n = std::fabs(GetDx(btmPt2->pt, p->pt)); + return (dx1p >= dx2p && dx1p >= dx2n) || (dx1n >= dx2p && dx1n >= dx2n); +} +//------------------------------------------------------------------------------ + +OutPt* GetBottomPt(OutPt *pp) +{ + OutPt* dups = 0; OutPt* p = pp->next; - OutPt* result = pp; while (p != pp) { - if (p->pt.Y > result->pt.Y) result = p; - else if (p->pt.Y == result->pt.Y && p->pt.X < result->pt.X) result = p; + if (p->pt.Y > pp->pt.Y) + { + pp = p; + dups = 0; + } + else if (p->pt.Y == pp->pt.Y && p->pt.X <= pp->pt.X) + { + if (p->pt.X < pp->pt.X) + { + dups = 0; + pp = p; + } else + { + if (p->next != pp && p->prev != pp) dups = p; + } + } p = p->next; } - return result; + if (dups) + { + //there appears to be at least 2 vertices at bottomPt so ... + while (dups != p) + { + if (!FirstIsBottomPt(p, dups)) pp = dups; + dups = dups->next; + while (!PointsEqual(dups->pt, pp->pt)) dups = dups->next; + } + } + return pp; } //------------------------------------------------------------------------------ @@ -806,11 +863,9 @@ bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType) { if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal) { - if (m_UseFullRange) + if (Abs(pg[i].X) > hiRange || Abs(pg[i].Y) > hiRange) throw "Coordinate exceeds range bounds"; maxVal = hiRange; - if (Abs(pg[i].X) > maxVal || Abs(pg[i].Y) > maxVal) - throw "Coordinate exceeds range bounds"; m_UseFullRange = true; } @@ -824,7 +879,7 @@ bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType) if (j < 2) return false; len = j+1; - for (;;) + while (len > 2) { //nb: test for point equality before testing slopes ... if (PointsEqual(p[j], p[0])) j--; @@ -837,9 +892,8 @@ bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType) for (int i = 2; i <= j; ++i) p[i-1] = p[i]; j--; } - //exit loop if nothing is changed or there are too few vertices ... - if (j == len-1 || j < 2) break; - len = j +1; + else break; + len--; } if (len < 3) return false; @@ -962,9 +1016,9 @@ TEdge* ClipperBase::AddBoundsToLML(TEdge *e) bool ClipperBase::AddPolygons(const Polygons &ppg, PolyType polyType) { - bool result = true; + bool result = false; for (Polygons::size_type i = 0; i < ppg.size(); ++i) - if (AddPolygon(ppg[i], polyType)) result = false; + if (AddPolygon(ppg[i], polyType)) result = true; return result; } //------------------------------------------------------------------------------ @@ -1167,7 +1221,7 @@ bool PolySort(OutRec *or1, OutRec *or2) { if (or1->pts != or2->pts) { - if (or1->pts) return true; else return false; + return or1->pts ? true : false; } else return false; } @@ -1181,8 +1235,7 @@ bool PolySort(OutRec *or1, OutRec *or2) int result = i1 - i2; if (result == 0 && (or1->isHole != or2->isHole)) { - if (or1->isHole) return false; - else return true; + return or1->isHole ? false : true; } else return result < 0; } @@ -1252,6 +1305,14 @@ bool Clipper::ExecuteInternal(bool fixHoleLinkages) FixupOutPolygon(*outRec); if (!outRec->pts) continue; if (outRec->isHole && fixHoleLinkages) FixHoleLinkage(outRec); + + if (outRec->bottomPt == outRec->bottomFlag && + (Orientation(outRec, m_UseFullRange) != (Area(*outRec, m_UseFullRange) > 0))) + { + DisposeBottomPt(*outRec); + FixupOutPolygon(*outRec); + }; + if (outRec->isHole == (m_ReverseOutput ^ Orientation(outRec, m_UseFullRange))) ReversePolyPtLinks(*outRec->pts); @@ -1312,10 +1373,10 @@ void Clipper::DisposeAllPolyPts(){ } //------------------------------------------------------------------------------ -void Clipper::DisposeOutRec(PolyOutList::size_type index, bool ignorePts) +void Clipper::DisposeOutRec(PolyOutList::size_type index) { OutRec *outRec = m_PolyOuts[index]; - if (!ignorePts && outRec->pts) DisposeOutPts(outRec->pts); + if (outRec->pts) DisposeOutPts(outRec->pts); delete outRec; m_PolyOuts[index] = 0; } @@ -1481,7 +1542,7 @@ void Clipper::AddLocalMinPoly(TEdge *e1, TEdge *e2, const IntPoint &pt) TEdge *e, *prevE; if( NEAR_EQUAL(e2->dx, HORIZONTAL) || ( e1->dx > e2->dx ) ) { - AddOutPt( e1, e2, pt ); + AddOutPt( e1, pt ); e2->outIdx = e1->outIdx; e1->side = esLeft; e2->side = esRight; @@ -1492,7 +1553,7 @@ void Clipper::AddLocalMinPoly(TEdge *e1, TEdge *e2, const IntPoint &pt) prevE = e->prevInAEL; } else { - AddOutPt( e2, e1, pt ); + AddOutPt( e2, pt ); e1->outIdx = e2->outIdx; e1->side = esRight; e2->side = esLeft; @@ -1511,14 +1572,16 @@ void Clipper::AddLocalMinPoly(TEdge *e1, TEdge *e2, const IntPoint &pt) void Clipper::AddLocalMaxPoly(TEdge *e1, TEdge *e2, const IntPoint &pt) { - AddOutPt( e1, 0, pt ); + AddOutPt( e1, pt ); if( e1->outIdx == e2->outIdx ) { e1->outIdx = -1; e2->outIdx = -1; } - else - AppendPolygon( e1, e2 ); + else if (e1->outIdx < e2->outIdx) + AppendPolygon(e1, e2); + else + AppendPolygon(e2, e1); } //------------------------------------------------------------------------------ @@ -1829,8 +1892,8 @@ void Clipper::IntersectEdges(TEdge *e1, TEdge *e2, AddLocalMinPoly(e1, e2, pt); break; case ctDifference: - if (e1->polyType == ptClip && e1Wc2 > 0 && e2Wc2 > 0 || - e1->polyType == ptSubject && e1Wc2 <= 0 && e2Wc2 <= 0) + if (((e1->polyType == ptClip) && (e1Wc2 > 0) && (e2Wc2 > 0)) || + ((e1->polyType == ptSubject) && (e1Wc2 <= 0) && (e2Wc2 <= 0))) AddLocalMinPoly(e1, e2, pt); break; case ctXor: @@ -1871,24 +1934,6 @@ void Clipper::SetHoleState(TEdge *e, OutRec *outRec) } //------------------------------------------------------------------------------ -bool GetNextNonDupOutPt(OutPt* pp, OutPt*& next) -{ - next = pp->next; - while (next != pp && PointsEqual(pp->pt, next->pt)) - next = next->next; - return next != pp; -} -//------------------------------------------------------------------------------ - -bool GetPrevNonDupOutPt(OutPt* pp, OutPt*& prev) -{ - prev = pp->prev; - while (prev != pp && PointsEqual(pp->pt, prev->pt)) - prev = prev->prev; - return prev != pp; -} -//------------------------------------------------------------------------------ - OutRec* GetLowermostRec(OutRec *outRec1, OutRec *outRec2) { //work out which polygon fragment has the correct hole state ... @@ -1898,21 +1943,10 @@ OutRec* GetLowermostRec(OutRec *outRec1, OutRec *outRec2) else if (outPt1->pt.Y < outPt2->pt.Y) return outRec2; else if (outPt1->pt.X < outPt2->pt.X) return outRec1; else if (outPt1->pt.X > outPt2->pt.X) return outRec2; - else if (outRec1->bottomE2 == 0) return outRec2; - else if (outRec2->bottomE2 == 0) return outRec1; - else - { - long64 y1 = std::max(outRec1->bottomE1->ybot, outRec1->bottomE2->ybot); - long64 y2 = std::max(outRec2->bottomE1->ybot, outRec2->bottomE2->ybot); - if (y2 == y1 || (y1 > outPt1->pt.Y && y2 > outPt1->pt.Y)) - { - double dx1 = std::max(outRec1->bottomE1->dx, outRec1->bottomE2->dx); - double dx2 = std::max(outRec2->bottomE1->dx, outRec2->bottomE2->dx); - if (dx2 > dx1) return outRec2; else return outRec1; - } - else if (y2 > y1) return outRec2; - else return outRec1; - } + else if (outPt1->next == outPt1) return outRec2; + else if (outPt2->next == outPt2) return outRec1; + else if (FirstIsBottomPt(outPt1, outPt2)) return outRec1; + else return outRec2; } //------------------------------------------------------------------------------ @@ -1921,13 +1955,11 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2) //get the start and ends of both output polygons ... OutRec *outRec1 = m_PolyOuts[e1->outIdx]; OutRec *outRec2 = m_PolyOuts[e2->outIdx]; - OutRec *holeStateRec = GetLowermostRec(outRec1, outRec2); - //fixup hole status ... - if (holeStateRec == outRec2) - outRec1->isHole = outRec2->isHole; - else - outRec2->isHole = outRec1->isHole; + OutRec *holeStateRec; + if (outRec1->FirstLeft == outRec2) holeStateRec = outRec2; + else if (outRec2->FirstLeft == outRec1) holeStateRec = outRec1; + else holeStateRec = GetLowermostRec(outRec1, outRec2); OutPt* p1_lft = outRec1->pts; OutPt* p1_rt = p1_lft->prev; @@ -1982,11 +2014,9 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2) { outRec1->bottomPt = outRec2->bottomPt; outRec1->bottomPt->idx = outRec1->idx; - outRec1->bottomE1 = outRec2->bottomE1; - outRec1->bottomE2 = outRec2->bottomE2; - if (outRec2->FirstLeft != outRec1) outRec1->FirstLeft = outRec2->FirstLeft; + outRec1->isHole = outRec2->isHole; } outRec2->pts = 0; outRec2->bottomPt = 0; @@ -2032,11 +2062,26 @@ OutRec* Clipper::CreateOutRec() result->AppendLink = 0; result->pts = 0; result->bottomPt = 0; + result->sides = esNeither; + result->bottomFlag = 0; + return result; } //------------------------------------------------------------------------------ -void Clipper::AddOutPt(TEdge *e, TEdge *altE, const IntPoint &pt) +void Clipper::DisposeBottomPt(OutRec &outRec) +{ + OutPt* next = outRec.bottomPt->next; + OutPt* prev = outRec.bottomPt->prev; + if (outRec.pts == outRec.bottomPt) outRec.pts = next; + delete outRec.bottomPt; + next->prev = prev; + prev->next = next; + outRec.bottomPt = next; +} +//------------------------------------------------------------------------------ + +void Clipper::AddOutPt(TEdge *e, const IntPoint &pt) { bool ToFront = (e->side == esLeft); if( e->outIdx < 0 ) @@ -2047,8 +2092,6 @@ void Clipper::AddOutPt(TEdge *e, TEdge *altE, const IntPoint &pt) e->outIdx = outRec->idx; OutPt* op = new OutPt; outRec->pts = op; - outRec->bottomE1 = e; - outRec->bottomE2 = altE; outRec->bottomPt = op; op->pt = pt; op->idx = outRec->idx; @@ -2061,16 +2104,59 @@ void Clipper::AddOutPt(TEdge *e, TEdge *altE, const IntPoint &pt) OutPt* op = outRec->pts; if ((ToFront && PointsEqual(pt, op->pt)) || (!ToFront && PointsEqual(pt, op->prev->pt))) return; + + if ((e->side | outRec->sides) != outRec->sides) + { + //check for 'rounding' artefacts ... + if (outRec->sides == esNeither && pt.Y == op->pt.Y) + if (ToFront) + { + if (pt.X == op->pt.X +1) return; //ie wrong side of bottomPt + } + else if (pt.X == op->pt.X -1) return; //ie wrong side of bottomPt + + outRec->sides = (EdgeSide)(outRec->sides | e->side); + if (outRec->sides == esBoth) + { + //A vertex from each side has now been added. + //Vertices of one side of an output polygon are quite commonly close to + //or even 'touching' edges of the other side of the output polygon. + //Very occasionally vertices from one side can 'cross' an edge on the + //the other side. The distance 'crossed' is always less that a unit + //and is purely an artefact of coordinate rounding. Nevertheless, this + //results in very tiny self-intersections. Because of the way + //orientation is calculated, even tiny self-intersections can cause + //the Orientation function to return the wrong result. Therefore, it's + //important to ensure that any self-intersections close to BottomPt are + //detected and removed before orientation is assigned. + + OutPt *opBot, *op2; + if (ToFront) + { + opBot = outRec->pts; + op2 = opBot->next; //op2 == right side + if (opBot->pt.Y != op2->pt.Y && opBot->pt.Y != pt.Y && + ((opBot->pt.X - pt.X)/(opBot->pt.Y - pt.Y) < + (opBot->pt.X - op2->pt.X)/(opBot->pt.Y - op2->pt.Y))) + outRec->bottomFlag = opBot; + } else + { + opBot = outRec->pts->prev; + op2 = opBot->prev; //op2 == left side + if (opBot->pt.Y != op2->pt.Y && opBot->pt.Y != pt.Y && + ((opBot->pt.X - pt.X)/(opBot->pt.Y - pt.Y) > + (opBot->pt.X - op2->pt.X)/(opBot->pt.Y - op2->pt.Y))) + outRec->bottomFlag = opBot; + } + } + } + OutPt* op2 = new OutPt; op2->pt = pt; op2->idx = outRec->idx; if (op2->pt.Y == outRec->bottomPt->pt.Y && op2->pt.X < outRec->bottomPt->pt.X) - { - outRec->bottomPt = op2; - outRec->bottomE1 = e; - outRec->bottomE2 = altE; - } + outRec->bottomPt = op2; op2->next = op; op2->prev = op->prev; op2->prev->next = op2; @@ -2225,8 +2311,7 @@ void Clipper::SwapPositionsInSEL(TEdge *edge1, TEdge *edge2) TEdge* GetNextInAEL(TEdge *e, Direction dir) { - if( dir == dLeftToRight ) return e->nextInAEL; - else return e->prevInAEL; + return dir == dLeftToRight ? e->nextInAEL : e->prevInAEL; } //------------------------------------------------------------------------------ @@ -2319,7 +2404,7 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge) if( horzEdge->nextInLML ) { if( horzEdge->outIdx >= 0 ) - AddOutPt( horzEdge, 0, IntPoint(horzEdge->xtop, horzEdge->ytop)); + AddOutPt( horzEdge, IntPoint(horzEdge->xtop, horzEdge->ytop)); UpdateEdgeIntoAEL( horzEdge ); } else @@ -2443,12 +2528,12 @@ bool Process1Before2(IntersectNode &node1, IntersectNode &node2) if (node1.edge1 == node2.edge1 || node1.edge2 == node2.edge1) { result = node2.pt.X > node1.pt.X; - if (node2.edge1->dx > 0) return !result; else return result; + return node2.edge1->dx > 0 ? !result : result; } else if (node1.edge1 == node2.edge2 || node1.edge2 == node2.edge2) { result = node2.pt.X > node1.pt.X; - if (node2.edge2->dx > 0) return !result; else return result; + return node2.edge2->dx > 0 ? !result : result; } else return node2.pt.X > node1.pt.X; } @@ -2542,7 +2627,7 @@ void Clipper::ProcessEdgesAtTopOfScanbeam(const long64 topY) { if (e->outIdx >= 0) { - AddOutPt(e, 0, IntPoint(e->xtop, e->ytop)); + AddOutPt(e, IntPoint(e->xtop, e->ytop)); for (HorzJoinList::size_type i = 0; i < m_HorizJoins.size(); ++i) { @@ -2578,7 +2663,7 @@ void Clipper::ProcessEdgesAtTopOfScanbeam(const long64 topY) { if( IsIntermediate( e, topY ) ) { - if( e->outIdx >= 0 ) AddOutPt(e, 0, IntPoint(e->xtop,e->ytop)); + if( e->outIdx >= 0 ) AddOutPt(e, IntPoint(e->xtop,e->ytop)); UpdateEdgeIntoAEL(e); //if output polygons share an edge, they'll need joining later ... @@ -2588,18 +2673,18 @@ void Clipper::ProcessEdgesAtTopOfScanbeam(const long64 topY) IntPoint(e->xbot,e->ybot), IntPoint(e->prevInAEL->xtop, e->prevInAEL->ytop), m_UseFullRange)) { - AddOutPt(e->prevInAEL, 0, IntPoint(e->xbot, e->ybot)); + AddOutPt(e->prevInAEL, IntPoint(e->xbot, e->ybot)); AddJoin(e, e->prevInAEL); } else if (e->outIdx >= 0 && e->nextInAEL && e->nextInAEL->outIdx >= 0 && e->nextInAEL->ycurr > e->nextInAEL->ytop && - e->nextInAEL->ycurr < e->nextInAEL->ybot && + e->nextInAEL->ycurr <= e->nextInAEL->ybot && e->nextInAEL->xcurr == e->xbot && e->nextInAEL->ycurr == e->ybot && SlopesEqual(IntPoint(e->xbot,e->ybot), IntPoint(e->xtop, e->ytop), IntPoint(e->xbot,e->ybot), IntPoint(e->nextInAEL->xtop, e->nextInAEL->ytop), m_UseFullRange)) { - AddOutPt(e->nextInAEL, 0, IntPoint(e->xbot, e->ybot)); + AddOutPt(e->nextInAEL, IntPoint(e->xbot, e->ybot)); AddJoin(e, e->nextInAEL); } } @@ -2632,13 +2717,7 @@ void Clipper::FixupOutPolygon(OutRec &outRec) lastOK = 0; OutPt *tmp = pp; if (pp == outRec.bottomPt) - { - if (tmp->prev->pt.Y > tmp->next->pt.Y) - outRec.bottomPt = tmp->prev; else - outRec.bottomPt = tmp->next; - outRec.pts = outRec.bottomPt; - outRec.bottomPt->idx = outRec.idx; - } + outRec.bottomPt = 0; //flags need for updating pp->prev->next = pp->next; pp->next->prev = pp->prev; pp = pp->prev; @@ -2651,6 +2730,11 @@ void Clipper::FixupOutPolygon(OutRec &outRec) pp = pp->next; } } + if (!outRec.bottomPt) { + outRec.bottomPt = GetBottomPt(pp); + outRec.bottomPt->idx = outRec.idx; + outRec.pts = outRec.bottomPt; + } } //------------------------------------------------------------------------------ @@ -2775,8 +2859,7 @@ bool Clipper::FixupIntersections() bool E2InsertsBeforeE1(TEdge &e1, TEdge &e2) { - if (e2.xcurr == e1.xcurr) return e2.dx > e1.dx; - else return e2.xcurr < e1.xcurr; + return e2.xcurr == e1.xcurr ? e2.dx > e1.dx : e2.xcurr < e1.xcurr; } //------------------------------------------------------------------------------ @@ -2808,7 +2891,7 @@ void Clipper::InsertEdgeIntoAEL(TEdge *edge) void Clipper::DoEdge1(TEdge *edge1, TEdge *edge2, const IntPoint &pt) { - AddOutPt(edge1, edge2, pt); + AddOutPt(edge1, pt); SwapSides(*edge1, *edge2); SwapPolyIndexes(*edge1, *edge2); } @@ -2816,7 +2899,7 @@ void Clipper::DoEdge1(TEdge *edge1, TEdge *edge2, const IntPoint &pt) void Clipper::DoEdge2(TEdge *edge1, TEdge *edge2, const IntPoint &pt) { - AddOutPt(edge2, edge1, pt); + AddOutPt(edge2, pt); SwapSides(*edge1, *edge2); SwapPolyIndexes(*edge1, *edge2); } @@ -2824,8 +2907,8 @@ void Clipper::DoEdge2(TEdge *edge1, TEdge *edge2, const IntPoint &pt) void Clipper::DoBothEdges(TEdge *edge1, TEdge *edge2, const IntPoint &pt) { - AddOutPt(edge1, edge2, pt); - AddOutPt(edge2, edge1, pt); + AddOutPt(edge1, pt); + AddOutPt(edge2, pt); SwapSides( *edge1 , *edge2 ); SwapPolyIndexes( *edge1 , *edge2 ); } @@ -2929,14 +3012,14 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages) { //instead of joining two polygons, we've just created a new one by //splitting one polygon into two. - outRec1->pts = PolygonBottom(p1); + outRec1->pts = GetBottomPt(p1); outRec1->bottomPt = outRec1->pts; outRec1->bottomPt->idx = outRec1->idx; outRec2 = CreateOutRec(); m_PolyOuts.push_back(outRec2); outRec2->idx = (int)m_PolyOuts.size()-1; j->poly2Idx = outRec2->idx; - outRec2->pts = PolygonBottom(p2); + outRec2->pts = GetBottomPt(p2); outRec2->bottomPt = outRec2->pts; outRec2->bottomPt->idx = outRec2->idx; @@ -2986,20 +3069,14 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages) //make sure any holes contained by outRec2 now link to outRec1 ... if (fixHoleLinkages) CheckHoleLinkages2(outRec1, outRec2); - //since it's possible for outRec2->bottomPt to be cleaned up - //in FixupOutPolygon() below, we need to keep a copy of it ... - IntPoint ip = outRec2->bottomPt->pt; - //now cleanup redundant edges too ... FixupOutPolygon(*outRec1); - if (outRec1->pts) { - //sort out hole vs outer and then recheck orientation ... - if (outRec1->isHole != outRec2->isHole && (ip.Y > outRec1->bottomPt->pt.Y || - (ip.Y == outRec1->bottomPt->pt.Y && ip.X < outRec1->bottomPt->pt.X))) - outRec1->isHole = outRec2->isHole; - if (outRec1->isHole == Orientation(outRec1, m_UseFullRange)) - ReversePolyPtLinks(*outRec1->pts); + if (outRec1->pts) + { + outRec1->isHole = !Orientation(outRec1, m_UseFullRange); + if (outRec1->isHole && !outRec1->FirstLeft) + outRec1->FirstLeft = outRec2->FirstLeft; } //delete the obsolete pointer ... @@ -3048,12 +3125,13 @@ struct DoublePoint Polygon BuildArc(const IntPoint &pt, const double a1, const double a2, const double r) { - int steps = std::max(6, int(std::sqrt(std::fabs(r)) * std::fabs(a2 - a1))); - Polygon result(steps); - int n = steps - 1; - double da = (a2 - a1) / n; + long64 steps = std::max(6, int(std::sqrt(std::fabs(r)) * std::fabs(a2 - a1))); + if (steps > 0x100000) steps = 0x100000; + int n = (unsigned)steps; + Polygon result(n); + double da = (a2 - a1) / (n -1); double a = a1; - for (int i = 0; i <= n; ++i) + for (int i = 0; i < n; ++i) { result[i].X = pt.X + Round(std::cos(a)*r); result[i].Y = pt.Y + Round(std::sin(a)*r); @@ -3206,19 +3284,18 @@ void DoSquare(double mul = 1.0) (long64)Round(m_p[m_i][m_j].Y + normals[m_k].Y * m_delta)); IntPoint pt2 = IntPoint((long64)Round(m_p[m_i][m_j].X + normals[m_j].X * m_delta), (long64)Round(m_p[m_i][m_j].Y + normals[m_j].Y * m_delta)); - if ((normals[m_k].X * normals[m_j].Y - normals[m_j].X * normals[m_k].Y) * m_delta >= 0) + double sinAngle = normals[m_k].X * normals[m_j].Y - normals[m_j].X * normals[m_k].Y; + if (sinAngle * m_delta >= 0) { - double a1 = std::atan2(normals[m_k].Y, normals[m_k].X); - double a2 = std::atan2(-normals[m_j].Y, -normals[m_j].X); - a1 = std::fabs(a2 - a1); - if (a1 > pi) a1 = pi * 2 - a1; - double dx = std::tan((pi - a1)/4) * std::fabs(m_delta * mul); - pt1 = IntPoint((long64)(pt1.X -normals[m_k].Y * dx), - (long64)(pt1.Y + normals[m_k].X * dx)); - AddPoint(pt1); - pt2 = IntPoint((long64)(pt2.X + normals[m_j].Y * dx), - (long64)(pt2.Y -normals[m_j].X * dx)); - AddPoint(pt2); + //occasionally (due to floating point math) sinAngle can be > 1 so ... + if (sinAngle > 1) sinAngle = 1; else if (sinAngle < -1) sinAngle = -1; + double dx = tan((pi - asin(sinAngle))/4) * abs(m_delta*mul); + pt1 = IntPoint((long64)(pt1.X -normals[m_k].Y * dx), + (long64)(pt1.Y + normals[m_k].X * dx)); + AddPoint(pt1); + pt2 = IntPoint((long64)(pt2.X + normals[m_j].Y * dx), + (long64)(pt2.Y -normals[m_j].X * dx)); + AddPoint(pt2); } else { diff --git a/src/Lib/Polygon/clipper.hpp b/src/Lib/Polygon/clipper.hpp index 92b58653..26ad6281 100644 --- a/src/Lib/Polygon/clipper.hpp +++ b/src/Lib/Polygon/clipper.hpp @@ -1,8 +1,8 @@ /******************************************************************************* * * * Author : Angus Johnson * -* Version : 4.7.3 * -* Date : 7 March 2012 * +* Version : 4.8.5 * +* Date : 15 July 2012 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2012 * * * @@ -87,7 +87,7 @@ void ReversePoints(Polygon& p); void ReversePoints(Polygons& p); //used internally ... -enum EdgeSide { esLeft, esRight }; +enum EdgeSide { esNeither = 0, esLeft = 1, esRight = 2, esBoth = 3 }; enum IntersectProtects { ipNone = 0, ipLeft = 1, ipRight = 2, ipBoth = 3 }; struct TEdge { @@ -142,8 +142,8 @@ struct OutRec { OutRec *AppendLink; OutPt *pts; OutPt *bottomPt; - TEdge *bottomE1; - TEdge *bottomE2; + OutPt *bottomFlag; + EdgeSide sides; }; struct OutPt { @@ -259,9 +259,10 @@ private: void IntersectEdges(TEdge *e1, TEdge *e2, const IntPoint &pt, IntersectProtects protects); OutRec* CreateOutRec(); - void AddOutPt(TEdge *e, TEdge *altE, const IntPoint &pt); + void AddOutPt(TEdge *e, const IntPoint &pt); + void DisposeBottomPt(OutRec &outRec); void DisposeAllPolyPts(); - void DisposeOutRec(PolyOutList::size_type index, bool ignorePts = false); + void DisposeOutRec(PolyOutList::size_type index); bool ProcessIntersections(const long64 botY, const long64 topY); void AddIntersectNode(TEdge *e1, TEdge *e2, const IntPoint &pt); void BuildIntersectList(const long64 botY, const long64 topY); diff --git a/src/Lib/Polygon/polygon.cxx b/src/Lib/Polygon/polygon.cxx index 25a0f966..89e8f3ab 100644 --- a/src/Lib/Polygon/polygon.cxx +++ b/src/Lib/Polygon/polygon.cxx @@ -709,6 +709,10 @@ TGPolygon tgPolygonDiffClipper( const TGPolygon& subject, const TGPolygon& clip return polygon_clip_clipper( POLY_DIFF, subject, clip ); } +TGPolygon tgPolygonIntClipper( const TGPolygon& subject, const TGPolygon& clip ) { + return polygon_clip_clipper( POLY_INT, subject, clip ); +} + TGPolygon tgPolygonUnionClipper( const TGPolygon& subject, const TGPolygon& clip ) { return polygon_clip_clipper( POLY_UNION, subject, clip ); } diff --git a/src/Lib/Polygon/polygon.hxx b/src/Lib/Polygon/polygon.hxx index d77d7639..96ea49ea 100644 --- a/src/Lib/Polygon/polygon.hxx +++ b/src/Lib/Polygon/polygon.hxx @@ -271,6 +271,9 @@ TGPolygon tgPolygonUnion( const TGPolygon& subject, const TGPolygon& clip ); // Difference TGPolygon tgPolygonDiffClipper( const TGPolygon& subject, const TGPolygon& clip ); +// Intersection +TGPolygon tgPolygonIntClipper( const TGPolygon& subject, const TGPolygon& clip ); + // Union TGPolygon tgPolygonUnionClipper( const TGPolygon& subject, const TGPolygon& clip );