From 465fccc6a03bc84605f453bdd27cc46b74807165 Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski <pete@loft.sadrohome> Date: Thu, 2 Aug 2012 21:13:33 -0400 Subject: [PATCH] checkpoint - madeira builds fine with individual line segments --- src/Lib/Geometry/util.cxx | 89 ++++++++++ src/Lib/Geometry/util.hxx | 1 + src/Lib/Polygon/chop-bin.cxx | 284 ++++++++++++++++++++++++++++++ src/Lib/Polygon/chop.hxx | 3 + src/Prep/OGRDecode/ogr-decode.cxx | 81 ++++++++- 5 files changed, 456 insertions(+), 2 deletions(-) diff --git a/src/Lib/Geometry/util.cxx b/src/Lib/Geometry/util.cxx index 6b5b36b7..376b59eb 100644 --- a/src/Lib/Geometry/util.cxx +++ b/src/Lib/Geometry/util.cxx @@ -414,6 +414,95 @@ Point3D midpoint( Point3D p0, Point3D p1 ) return Point3D( (p0.x() + p1.x()) / 2, (p0.y() + p1.y()) / 2, (p0.z() + p1.z()) / 2 ); } +void +makePolygons (const Line &line, double width, poly_list& polys) +{ + int nPoints = line.getPointCount(); + int i; + int turn_dir; + + Point3D cur_inner; + Point3D cur_outer; + Point3D prev_inner = Point3D(0.0f, 0.0f, 0.0f); + Point3D prev_outer = Point3D(0.0f, 0.0f, 0.0f); + + double heading = 0.0f; + double az2 = 0.0f; + double dist = 0.0f; + double pt_x = 0.0f; + double pt_y = 0.0f; + + TGPolygon poly; + + // generate poly and texparam lists for each line segment + for (i=0; i<nPoints; i++) + { + turn_dir = 0; + + SG_LOG(SG_GENERAL, SG_DEBUG, "makePolygonsTP: calculating offsets for segment " << i); + + // for each point on the PointsList, generate a quad from + // start to next, offset by 1/2 width from the edge + if (i == 0) + { + // first point on the list - offset heading is 90deg + cur_outer = OffsetPointFirst( line.getPoint(i), line.getPoint(i+1), -width/2.0f ); + cur_inner = OffsetPointFirst( line.getPoint(i), line.getPoint(i+1), width/2.0f ); + } + else if (i == nPoints-1) + { + // last point on the list - offset heading is 90deg + cur_outer = OffsetPointLast( line.getPoint(i-1), line.getPoint(i), -width/2.0f ); + cur_inner = OffsetPointLast( line.getPoint(i-1), line.getPoint(i), width/2.0f ); + } + else + { + // middle section + cur_outer = OffsetPointMiddle( line.getPoint(i-1), line.getPoint(i), line.getPoint(i+1), -width/2.0f, turn_dir ); + cur_inner = OffsetPointMiddle( line.getPoint(i-1), line.getPoint(i), line.getPoint(i+1), width/2.0f, turn_dir ); + } + + if ( (prev_inner.x() != 0.0f) && (prev_inner.y() != 0.0f) ) + { + Point3D prev_mp = midpoint( prev_outer, prev_inner ); + Point3D cur_mp = midpoint( cur_outer, cur_inner ); + geo_inverse_wgs_84( prev_mp.y(), prev_mp.x(), cur_mp.y(), cur_mp.x(), &heading, &az2, &dist); + + poly.erase(); + + poly.add_node( 0, prev_inner ); + poly.add_node( 0, prev_outer ); + + // we need to extend one of the points so we're sure we don't create adjacent edges + if (turn_dir == 0) + { + // turned right - offset outer + geo_inverse_wgs_84( prev_outer.y(), prev_outer.x(), cur_outer.y(), cur_outer.x(), &heading, &az2, &dist); + geo_direct_wgs_84( cur_outer.y(), cur_outer.x(), heading, MP_STRETCH, &pt_y, &pt_x, &az2 ); + + poly.add_node( 0, Point3D( pt_x, pt_y, 0.0f) ); + //poly.add_node( 0, cur_outer ); + poly.add_node( 0, cur_inner ); + } + else + { + // turned left - offset inner + geo_inverse_wgs_84( prev_inner.y(), prev_inner.x(), cur_inner.y(), cur_inner.x(), &heading, &az2, &dist); + geo_direct_wgs_84( cur_inner.y(), cur_inner.x(), heading, MP_STRETCH, &pt_y, &pt_x, &az2 ); + + poly.add_node( 0, cur_outer ); + poly.add_node( 0, Point3D( pt_x, pt_y, 0.0f) ); + //poly.add_node( 0, cur_inner ); + } + + polys.push_back(poly); + } + + prev_outer = cur_outer; + prev_inner = cur_inner; + } +} + void makePolygonsTP (const Line &line, double width, poly_list& polys, texparams_list &tps) { diff --git a/src/Lib/Geometry/util.hxx b/src/Lib/Geometry/util.hxx index 6b219f8c..22ba6df0 100644 --- a/src/Lib/Geometry/util.hxx +++ b/src/Lib/Geometry/util.hxx @@ -93,6 +93,7 @@ void makePolygon (const Point3D &p, double width, TGPolygon &polygon); */ void makePolygon (const Line &line, double width, TGPolygon &polygon); +void makePolygons (const Line &line, double width, poly_list &segments); void makePolygonsTP (const Line &line, double width, poly_list &segments, texparams_list &tps); /** diff --git a/src/Lib/Polygon/chop-bin.cxx b/src/Lib/Polygon/chop-bin.cxx index c9c0180c..7c67f6bb 100644 --- a/src/Lib/Polygon/chop-bin.cxx +++ b/src/Lib/Polygon/chop-bin.cxx @@ -153,6 +153,132 @@ static void clip_and_write_poly( string root, long int p_index, } } +static void clip_and_write_polys_with_mask( string root, long int p_index, + const string &poly_type, + SGBucket b, const poly_list& segments, + bool preserve3d ) +{ + Point3D c, min, max, p; + + c = Point3D( b.get_center_lon(), b.get_center_lat(), 0 ); + double span = b.get_width(); + TGPolygon base, result; + char tile_name[256], poly_index[256]; + + // calculate bucket dimensions + if ( (c.y() >= -89.0) && (c.y() < 89.0) ) { + min.setx( c.x() - span / 2.0 ); + max.setx( c.x() + span / 2.0 ); + min.sety( c.y() - SG_HALF_BUCKET_SPAN ); + max.sety( c.y() + SG_HALF_BUCKET_SPAN ); + } else if ( c.y() < -89.0) { + min.setx( -90.0 ); + max.setx( -89.0 ); + min.sety( -180.0 ); + max.sety( 180.0 ); + } else if ( c.y() >= 89.0) { + min.setx( 89.0 ); + max.setx( 90.0 ); + min.sety( -180.0 ); + max.sety( 180.0 ); + } else { + SG_LOG( SG_GENERAL, SG_ALERT, + "Out of range latitude in clip_and_write_poly() = " << c.y() ); + } + + min.setz( 0.0 ); + max.setz( 0.0 ); + + SG_LOG( SG_GENERAL, SG_DEBUG, " (" << min << ") (" << max << ")" ); + + // set up clipping tile + base.add_node( 0, Point3D(min.x(), min.y(), 0) ); + base.add_node( 0, Point3D(max.x(), min.y(), 0) ); + base.add_node( 0, Point3D(max.x(), max.y(), 0) ); + base.add_node( 0, Point3D(min.x(), max.y(), 0) ); + + // SG_LOG( SG_GENERAL, SG_DEBUG, "base = 4 vertices" ); + + /* + FILE *bfp= fopen("base", "w"); + gpc_write_polygon(bfp, &base); + fclose(bfp); + */ + + poly_list clipped_shapes; + TGPolygon shape; + unsigned int s; + + for ( s = 0; s < segments.size(); s++ ) { + shape = segments[s]; + + // SG_LOG(SG_GENERAL, SG_DEBUG, "shape contours = " << shape.contours() ); + // for ( int ii = 0; ii < shape.contours(); ii++ ) + // SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << shape.get_hole_flag(ii) ); + + result = tgPolygonInt( base, shape ); + + // write_polygon(shape, "shape"); + // write_polygon(result, "result"); + + // SG_LOG(SG_GENERAL, SG_DEBUG, "result contours = " << result.contours() ); + // for ( int ii = 0; ii < result.contours(); ii++ ) { + // SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << result.get_hole_flag(ii) ); + // } + + if ( preserve3d ) + result.inherit_elevations( shape ); + + if ( result.contours() > 0 ) { + clipped_shapes.push_back( result ); + } + } + + // Now we can write the file + long int t_index = b.gen_index(); + string path = root + "/" + b.gen_base_path(); + + SGPath sgp( path ); + sgp.append( "dummy" ); + sgp.create_dir( 0755 ); + + sprintf( tile_name, "%ld", t_index ); + string polyfile = path + "/" + tile_name; + + sprintf( poly_index, "%ld", p_index ); + polyfile += "."; + polyfile += poly_index; + + FILE *rfp = fopen( polyfile.c_str(), "w" ); + + if ( preserve3d ) + fprintf( rfp, "#3D_WITH_MASK\n" ); + else + fprintf( rfp, "#2D_WITH_MASK\n" ); + fprintf( rfp, "%s\n", poly_type.c_str() ); + + fprintf( rfp, "%d\n", (int)clipped_shapes.size() ); + + for (s=0; s<clipped_shapes.size(); s++) { + result = clipped_shapes[s]; + + fprintf( rfp, "%d\n", result.contours() ); + for ( int i = 0; i < result.contours(); ++i ) { + fprintf( rfp, "%d\n", result.contour_size(i) ); + fprintf( rfp, "%d\n", result.get_hole_flag(i) ); + for ( int j = 0; j < result.contour_size(i); ++j ) { + p = result.get_pt( i, j ); + if ( preserve3d ) + fprintf( rfp, "%.15f %.15f %.15f\n", p.x(), p.y(), p.z() ); + else + fprintf( rfp, "%.15f %.15f\n", p.x(), p.y() ); + } + } + } + + fclose( rfp ); +} + static void clip_and_write_poly_tp( string root, long int p_index, const string &poly_type, SGBucket b, const TGPolygon& shape, @@ -401,6 +527,164 @@ void tgChopNormalPolygon( const string& path, const string& poly_type, } } +void tgChopNormalPolygonsWithMask(const std::string& path, const std::string& poly_type, + const poly_list& segments, bool preserve3d ) +{ + Point3D min, max, p; + TGPolygon shape; + long int index; + int i, j; + unsigned int s; + + // bail out immediately if poly_list is empty + if ( segments.size() == 0 ) + return; + + // get next polygon index + index = poly_index_next(); + + min = Point3D( 200.0 ); + max = Point3D( -200.0 ); + + for ( s = 0; s < segments.size(); s++) { + shape = segments[s]; + + // find min/max of polygon + for ( i = 0; i < shape.contours(); i++ ) { + for ( j = 0; j < shape.contour_size(i); j++ ) { + p = shape.get_pt( i, j ); + + if ( p.x() < min.x() ) min.setx( p.x() ); if ( p.y() < min.y() ) min.sety( p.y() ); if ( p.x() > max.x() ) max.setx( p.x() ); if ( p.y() > max.y() ) max.sety( p.y() ); + } + } + } + + SG_LOG( SG_GENERAL, SG_DEBUG, " min = " << min << " max = " << max ); + + // find buckets for min, and max points of convex hull. + // note to self: self, you should think about checking for + // polygons that span the date line + SGBucket b_min( min.x(), min.y() ); + SGBucket b_max( max.x(), max.y() ); + + SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket min = " << b_min ); + SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket max = " << b_max ); + + if ( b_min == b_max ) { + // shape entirely contained in a single bucket, write and bail + clip_and_write_polys_with_mask( path, index, poly_type, b_min, segments, preserve3d ); + return; + } + + SGBucket b_cur; + int dx, dy; + + sgBucketDiff(b_min, b_max, &dx, &dy); + SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" ); + SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); + + if ( (dx > 2880) || (dy > 1440) ) + throw sg_exception("something is really wrong in split_polygon()!!!!"); + + if ( dy <= 1 ) { + // we are down to at most two rows, write each column and then bail + double min_center_lat = b_min.get_center_lat(); + double min_center_lon = b_min.get_center_lon(); + for ( j = 0; j <= dy; ++j ) { + for ( i = 0; i <= dx; ++i ) { + b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j); + clip_and_write_polys_with_mask( path, index, poly_type, b_cur, segments, preserve3d ); + } + } + return; + } + + // we have two or more rows left, split in half (along a + // horizontal dividing line) and recurse with each half + + // find mid point (integer math) + int mid = (dy + 1) / 2 - 1; + + // determine horizontal clip line + SGBucket b_clip = sgBucketOffset(min.x(), min.y(), 0, mid); + double clip_line = b_clip.get_center_lat(); + if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN) + && (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) ) + clip_line += SG_HALF_BUCKET_SPAN; + else if ( clip_line < -89.0 ) + clip_line = -89.0; + else if ( clip_line >= 89.0 ) + clip_line = 90.0; + else { + SG_LOG( SG_GENERAL, SG_ALERT, + "Out of range latitude in clip_and_write_poly() = " + << clip_line ); + } + + { + // + // Crop bottom area (hopefully by putting this in it's own + // scope we can shorten the life of some really large data + // structures to reduce memory use) + // + SG_LOG( SG_GENERAL, SG_DEBUG, + "Generating bottom half (" << min.y() << "-" << + clip_line << ")" ); + + TGPolygon bottom, bottom_clip; + poly_list bottom_clip_list; + + bottom.erase(); + bottom.add_node( 0, Point3D(-180.0, min.y(), 0) ); + bottom.add_node( 0, Point3D(180.0, min.y(), 0) ); + bottom.add_node( 0, Point3D(180.0, clip_line, 0) ); + bottom.add_node( 0, Point3D(-180.0, clip_line, 0) ); + + for (s=0; s<segments.size(); s++) { + bottom_clip.erase(); + bottom_clip = tgPolygonInt( bottom, segments[s] ); + + if ( bottom_clip.contours() > 0 ) { + bottom_clip_list.push_back( bottom_clip ); + } + } + + tgChopNormalPolygonsWithMask( path, poly_type, bottom_clip_list, preserve3d ); + } + + { + // + // Crop top area (hopefully by putting this in it's own scope + // we can shorten the life of some really large data + // structures to reduce memory use) + // + + SG_LOG( SG_GENERAL, SG_DEBUG, + "Generating top half (" << clip_line << "-" << + max.y() << ")" ); + + TGPolygon top, top_clip; + poly_list top_clip_list; + + top.erase(); + top.add_node( 0, Point3D(-180.0, clip_line, 0) ); + top.add_node( 0, Point3D(180.0, clip_line, 0) ); + top.add_node( 0, Point3D(180.0, max.y(), 0) ); + top.add_node( 0, Point3D(-180.0, max.y(), 0) ); + + for (s=0; s<segments.size(); s++) { + top_clip.erase(); + top_clip = tgPolygonInt( top, segments[s] ); + + if ( top_clip.contours() > 0 ) { + top_clip_list.push_back( top_clip ); + } + } + + tgChopNormalPolygonsWithMask( path, poly_type, top_clip_list, preserve3d ); + } +} + // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) void tgChopNormalPolygonTP( const string& path, const string& poly_type, diff --git a/src/Lib/Polygon/chop.hxx b/src/Lib/Polygon/chop.hxx index 3ea2c382..a10e6bd3 100644 --- a/src/Lib/Polygon/chop.hxx +++ b/src/Lib/Polygon/chop.hxx @@ -37,6 +37,9 @@ void tgChopNormalPolygon( const std::string& path, const std::string& poly_type, const TGPolygon& shape, bool preserve3d ); +void tgChopNormalPolygonsWithMask(const std::string& path, const std::string& poly_type, + const poly_list& segments, bool preserve3d ); + // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) void tgChopNormalPolygonTP( const std::string& path, const std::string& poly_type, diff --git a/src/Prep/OGRDecode/ogr-decode.cxx b/src/Prep/OGRDecode/ogr-decode.cxx index 2a0316f9..28c5775d 100644 --- a/src/Prep/OGRDecode/ogr-decode.cxx +++ b/src/Prep/OGRDecode/ogr-decode.cxx @@ -110,6 +110,83 @@ void processLineString(OGRLineString* poGeometry, tgChopNormalPolygon(work_dir, area_type, shape, false); } +void processLineStringWithMask(OGRLineString* poGeometry, + const string& work_dir, + const string& area_type, + int width) +{ + poly_list segments; + TGPolygon segment; + tg::Line line; + Point3D p0, p1; + double heading, dist, az2; + double pt_x = 0.0f, pt_y = 0.0f; + int i, j, numPoints, numSegs; + double max_dist; + + numPoints = poGeometry->getNumPoints(); + if (numPoints < 2) { + SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" ); + return; + } + + max_dist = (double)width * 10.0f; + + // because vector data can generate adjacent polys, lets stretch the two enpoints by a little bit + p0 = Point3D(poGeometry->getX(0),poGeometry->getY(0),0); + p1 = Point3D(poGeometry->getX(1),poGeometry->getY(1),0); + + geo_inverse_wgs_84( p1.y(), p1.x(), p0.y(), p0.x(), &heading, &az2, &dist); + geo_direct_wgs_84( p0.y(), p0.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); + line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); + + // now add the middle points : if they are too far apart, add intermediate nodes + for ( i=1;i<numPoints-1;i++) { + p0 = Point3D(poGeometry->getX(i-1),poGeometry->getY(i-1),0); + p1 = Point3D(poGeometry->getX(i ),poGeometry->getY(i ),0); + geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); + + if (dist > max_dist) + { + numSegs = (dist / max_dist) + 1; + dist = dist / (double)numSegs; + + for (j=0; j<numSegs; j++) + { + geo_direct_wgs_84( p0.y(), p0.x(), heading, dist*(j+1), &pt_y, &pt_x, &az2 ); + line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); + } + } + else + { + line.addPoint(p1); + } + } + + // then stretch the last point + // because vector data can generate adjacent polys, lets stretch the two enpoints by a little bit + p0 = Point3D(poGeometry->getX(numPoints-2),poGeometry->getY(numPoints-2),0); + p1 = Point3D(poGeometry->getX(numPoints-1),poGeometry->getY(numPoints-1),0); + + geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); + geo_direct_wgs_84( p1.y(), p1.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); + line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); + + // make a plygons from the line segments + tg::makePolygons(line,width,segments); + +#if 1 + for ( i = 0; i < (int)segments.size(); ++i ) { + segment = segments[i]; + + tgChopNormalPolygon(work_dir, area_type, segment, false); + } +#else + tgChopNormalPolygonsWithMask(work_dir, area_type, segments, false); +#endif + +} + void processLineStringWithTextureInfo(OGRLineString* poGeometry, const string& work_dir, const string& area_type, @@ -406,7 +483,7 @@ void processLayer(OGRLayer* poLayer, if (texture_lines) { processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); } else { - processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); + processLineStringWithMask((OGRLineString*)poGeometry,work_dir,area_type_name,width); } break; } @@ -421,7 +498,7 @@ void processLayer(OGRLayer* poLayer, if (texture_lines) { processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); } else { - processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); + processLineStringWithMask((OGRLineString*)poGeometry,work_dir,area_type_name,width); } } break;