From bc3ecf8fe903dfd560a4efaf69f15505f78b1ff4 Mon Sep 17 00:00:00 2001 From: PSadrozinski Date: Thu, 22 Mar 2012 10:01:32 +0100 Subject: [PATCH] - Add --texture-line option to ogr-decode to generated polygons with texture parameters - change cout to SG_LOG so we can turn down the verbosity --- src/Lib/Array/array.cxx | 55 ++---- src/Lib/Geometry/trisegs.cxx | 44 ++--- src/Lib/Geometry/util.cxx | 292 ++++++++++++++++++++++++++++++ src/Lib/Geometry/util.hxx | 3 + src/Lib/Polygon/chop-bin.cxx | 271 +++++++++++++++++++++++++-- src/Lib/Polygon/chop.hxx | 6 +- src/Prep/OGRDecode/ogr-decode.cxx | 111 +++++++++++- 7 files changed, 694 insertions(+), 88 deletions(-) diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index 2f0afd96..fb5ac989 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -32,28 +32,25 @@ #include #include +#include #include #include #include -#include #include "array.hxx" using std::string; -using std::cout; -using std::endl; TGArray::TGArray( void ): array_in(NULL), fitted_in(NULL) { - // cout << "class TGArray CONstructor called." << endl; - //in_data = new int[ARRAY_SIZE_1][ARRAY_SIZE_1]; + SG_LOG(SG_GENERAL, SG_DEBUG, "class TGArray CONstructor called." ); in_data = new int*[ARRAY_SIZE_1]; - for (int i = 0; i < ARRAY_SIZE_1; i++) + for (int i = 0; i < ARRAY_SIZE_1; i++) { in_data[i] = new int[ARRAY_SIZE_1]; - // out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1]; + } } @@ -61,11 +58,10 @@ TGArray::TGArray( const string &file ): array_in(NULL), fitted_in(NULL) { - // cout << "class TGArray CONstructor called." << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, "class TGArray CONstructor called." ); in_data = new int* [ARRAY_SIZE_1]; for (int i = 0; i < ARRAY_SIZE_1; i++) in_data[i] = new int[ARRAY_SIZE_1]; - // out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1]; TGArray::open(file); } @@ -79,10 +75,10 @@ bool TGArray::open( const string& file_base ) { string array_name = file_base + ".arr.gz"; array_in = new sg_gzifstream( array_name ); if ( ! array_in->is_open() ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " Cannot open " << array_name ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Cannot open " << array_name ); success = false; } else { - SG_LOG(SG_GENERAL, SG_DEBUG, " Opening array data file: " << array_name ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Opening array data file: " << array_name ); } // open fitted data file @@ -93,7 +89,7 @@ bool TGArray::open( const string& file_base ) { // can do a really stupid/crude fit on the fly, but it will // not be nearly as nice as what the offline terrafit utility // would have produced. - SG_LOG(SG_GENERAL, SG_DEBUG, " Cannot open " << fitted_name ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Cannot open " << fitted_name ); } else { SG_LOG(SG_GENERAL, SG_DEBUG, " Opening fitted data file: " << fitted_name ); } @@ -131,11 +127,6 @@ TGArray::parse( SGBucket& b ) { SG_LOG(SG_GENERAL, SG_DEBUG, " cols = " << cols << " rows = " << rows ); SG_LOG(SG_GENERAL, SG_DEBUG, " col_step = " << col_step << " row_step = " << row_step ); - if ((cols>ARRAY_SIZE_1) || (rows>ARRAY_SIZE_1)) { - cout << "error, ARRAY_SIZE_1=" << ARRAY_SIZE_1 <<" is too small (cols=" << cols << ", rows=" << rows << "), aborting." << endl; - exit(-1); - } - for ( int i = 0; i < cols; i++ ) { for ( int j = 0; j < rows; j++ ) { *array_in >> in_data[i][j]; @@ -178,7 +169,7 @@ TGArray::parse( SGBucket& b ) { for ( int i = 0; i < fitted_size; ++i ) { *fitted_in >> x >> y >> z; fitted_list.push_back( Point3D(x, y, z) ); - // cout << " loading fitted = " << Point3D(x, y, z) << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, " loading fitted = " << Point3D(x, y, z) ); } } @@ -201,7 +192,7 @@ bool TGArray::write( const string root_dir, SGBucket& b ) { // write the file gzFile fp; if ( (fp = gzopen( array_file.c_str(), "wb9" )) == NULL ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "ERROR: cannot open " << array_file << " for writing!" ); + SG_LOG(SG_GENERAL, SG_ALERT, "ERROR: cannot open " << array_file << " for writing!" ); return false; } @@ -309,7 +300,7 @@ void TGArray::add_corner_node( int i, int j, double val ) { double x = (originx + i * col_step) / 3600.0; double y = (originy + j * row_step) / 3600.0; - // cout << "originx = " << originx << " originy = " << originy << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, "originx = " << originx << " originy = " << originy ); SG_LOG(SG_GENERAL, SG_DEBUG, "corner = " << Point3D(x, y, val) ); corner_list.push_back( Point3D(x, y, val) ); } @@ -319,7 +310,7 @@ void TGArray::add_corner_node( int i, int j, double val ) { void TGArray::add_fit_node( int i, int j, double val ) { double x = (originx + i * col_step) / 3600.0; double y = (originy + j * row_step) / 3600.0; - // cout << Point3D(x, y, val) << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, Point3D(x, y, val) ); fitted_list.push_back( Point3D(x, y, val) ); } @@ -557,8 +548,7 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { if ( (xindex < 0) || (xindex + 1 >= cols) || (yindex < 0) || (yindex + 1 >= rows) ) { - // cout << "WARNING: Attempt to interpolate value outside of array!!!" - // << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, "WARNING: Attempt to interpolate value outside of array!!!" ); return -9999; } @@ -567,7 +557,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { if ( dx > dy ) { // lower triangle - // printf(" Lower triangle\n"); x1 = xindex; y1 = yindex; @@ -581,11 +570,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { y3 = yindex + 1; z3 = in_data[x3][y3]; - // printf(" dx = %.2f dy = %.2f\n", dx, dy); - // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); - // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); - // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); - if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { // don't interpolate off a void return closest_nonvoid_elev( lon, lat ); @@ -594,8 +578,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { zA = dx * (z2 - z1) + z1; zB = dx * (z3 - z1) + z1; - // printf(" zA = %.2f zB = %.2f\n", zA, zB); - if ( dx > SG_EPSILON ) { elev = dy * (zB - zA) / dx + zA; } else { @@ -603,7 +585,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { } } else { // upper triangle - // printf(" Upper triangle\n"); x1 = xindex; y1 = yindex; @@ -616,11 +597,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { x3 = xindex + 1; y3 = yindex + 1; z3 = in_data[x3][y3]; - - // printf(" dx = %.2f dy = %.2f\n", dx, dy); - // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); - // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); - // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { // don't interpolate off a void @@ -630,9 +606,6 @@ double TGArray::altitude_from_grid( double lon, double lat ) const { zA = dy * (z2 - z1) + z1; zB = dy * (z3 - z1) + z1; - // printf(" zA = %.2f zB = %.2f\n", zA, zB ); - // printf(" xB - xA = %.2f\n", col_step * dy / row_step); - if ( dy > SG_EPSILON ) { elev = dx * (zB - zA) / dy + zA; } else { @@ -764,11 +737,9 @@ void TGArray::outputmesh_output_nodes( const string& fg_root, SGBucket& p ) TGArray::~TGArray( void ) { - // printf("class TGArray DEstructor called.\n"); for (int i = 0; i < ARRAY_SIZE_1; i++) delete [] in_data[i]; delete [] in_data; - // delete [] out_data; } diff --git a/src/Lib/Geometry/trisegs.cxx b/src/Lib/Geometry/trisegs.cxx index bcaa3b2c..adab2135 100644 --- a/src/Lib/Geometry/trisegs.cxx +++ b/src/Lib/Geometry/trisegs.cxx @@ -22,6 +22,8 @@ #include #include +#include + #include #include @@ -30,10 +32,6 @@ #include "trisegs.hxx" -using std::cout; -using std::endl; - - // Constructor TGTriSegments::TGTriSegments( void ) { } @@ -51,12 +49,11 @@ int TGTriSegments::unique_add( const TGTriSeg& s ) triseg_list_iterator current, last; int counter = 0; - // cout << s.get_n1() << "," << s.get_n2() << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, s.get_n1() << "," << s.get_n2() ); // check if segment has duplicated endpoints if ( s.get_n1() == s.get_n2() ) { - cout << "WARNING: ignoring null segment with the same " - << "point for both endpoints" << endl; + SG_LOG(SG_GENERAL, SG_ALERT, "WARNING: ignoring null segment with the same point for both endpoints" ); return -1; } @@ -65,7 +62,7 @@ int TGTriSegments::unique_add( const TGTriSeg& s ) last = seg_list.end(); for ( ; current != last; ++current ) { if ( s == *current ) { - // cout << "found an existing segment match" << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, "found an existing segment match" ); return counter; } @@ -95,7 +92,7 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, // bool temp = false; // if ( s == TGTriSeg( 170, 206 ) ) { - // cout << "this is it!" << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, "this is it!" ); // temp = true; // } @@ -118,7 +115,7 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, b = p1.y() - m * p1.x(); // if ( temp ) { - // cout << "m = " << m << " b = " << b << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, "m = " << m << " b = " << b ); // } current = nodes.begin(); @@ -129,15 +126,14 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, && (current->x() < (p1.x() - SG_EPSILON)) ) { // if ( temp ) { - // cout << counter << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, counter ); // } y_err = fabs(current->y() - (m * current->x() + b)); if ( y_err < FG_PROXIMITY_EPSILON ) { - //cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; - //cout << p0 << " < " << *current << " < " - // << p1 << endl; + //SG_LOG(SG_GENERAL, SG_DEBUG, "FOUND EXTRA SEGMENT NODE (Y)" ); + //SG_LOG(SG_GENERAL, SG_DEBUG, p0 << " < " << *current << " < " << p1 ); found_extra = true; if ( y_err < y_err_min ) { extra_index = counter; @@ -162,13 +158,13 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, // bool temp = true; // if ( temp ) { - // cout << "xdist = " << xdist << " ydist = " << ydist << endl; - // cout << " p0 = " << p0 << " p1 = " << p1 << endl; - // cout << " m1 = " << m1 << " b1 = " << b1 << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, "xdist = " << xdist << " ydist = " << ydist ); + // SG_LOG(SG_GENERAL, SG_DEBUG, " p0 = " << p0 << " p1 = " << p1 ); + // SG_LOG(SG_GENERAL, SG_DEBUG, " m1 = " << m1 << " b1 = " << b1 ); // } - // cout << " should = 0 = " << fabs(p0.x() - (m1 * p0.y() + b1)) << endl;; - // cout << " should = 0 = " << fabs(p1.x() - (m1 * p1.y() + b1)) << endl;; + // SG_LOG(SG_GENERAL, SG_DEBUG, " should = 0 = " << fabs(p0.x() - (m1 * p0.y() + b1)) ); + // SG_LOG(SG_GENERAL, SG_DEBUG, " should = 0 = " << fabs(p1.x() - (m1 * p1.y() + b1)) ); current = nodes.begin(); last = nodes.end(); @@ -180,13 +176,12 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, x_err = fabs(current->x() - (m1 * current->y() + b1)); // if ( temp ) { - // cout << " (" << counter << ") x_err = " << x_err << endl; + // SG_LOG(SG_GENERAL, SG_DEBUG, " (" << counter << ") x_err = " << x_err ); // } if ( x_err < FG_PROXIMITY_EPSILON ) { - //cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; - //cout << p0 << " < " << *current << " < " - // << p1 << endl; + //SG_LOG(SG_GENERAL, SG_DEBUG, "FOUND EXTRA SEGMENT NODE (X)" ); + //SG_LOG(SG_GENERAL, SG_DEBUG, p0 << " < " << *current << " < " << p1 ); found_extra = true; if ( x_err < x_err_min ) { extra_index = counter; @@ -200,8 +195,7 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, if ( found_extra ) { // recurse with two sub segments - cout << "dividing " << s.get_n1() << " " << extra_index - << " " << s.get_n2() << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, "dividing " << s.get_n1() << " " << extra_index << " " << s.get_n2() ); unique_divide_and_add( nodes, TGTriSeg( s.get_n1(), extra_index, s.get_boundary_marker() ) ); unique_divide_and_add( nodes, TGTriSeg( extra_index, s.get_n2(), diff --git a/src/Lib/Geometry/util.cxx b/src/Lib/Geometry/util.cxx index ad13c47d..6e335da1 100644 --- a/src/Lib/Geometry/util.cxx +++ b/src/Lib/Geometry/util.cxx @@ -11,9 +11,13 @@ #include #include +#include +#include #include +#define MP_STRETCH (0.000001) + using std::vector; namespace tg { @@ -211,6 +215,294 @@ makePolygon (const Line &line, double width, TGPolygon &polygon) polygon.set_hole_flag(0, 0); // mark as solid } +inline double CalculateTheta( Point3D p0, Point3D p1, Point3D p2 ) +{ + Point3D u, v; + double udist, vdist, uv_dot, tmp; + + // u . v = ||u|| * ||v|| * cos(theta) + + u = p1 - p0; + udist = sqrt( u.x() * u.x() + u.y() * u.y() ); + // printf("udist = %.6f\n", udist); + + v = p1 - p2; + vdist = sqrt( v.x() * v.x() + v.y() * v.y() ); + // printf("vdist = %.6f\n", vdist); + + uv_dot = u.x() * v.x() + u.y() * v.y(); + // printf("uv_dot = %.6f\n", uv_dot); + + tmp = uv_dot / (udist * vdist); + // printf("tmp = %.6f\n", tmp); + + return acos(tmp); +} + +Point3D OffsetPointMiddle( Point3D prev, Point3D cur, Point3D next, double offset_by, int& turn_dir ) +{ + double offset_dir; + double next_dir; + double az2; + double dist; + double theta; + double pt_x = 0, pt_y = 0; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find average angle for contour: prev (" << prev << "), " + "cur (" << cur << "), " + "next (" << next << ")" ); + + + // first, find if the line turns left or right ar src + // for this, take the cross product of the vectors from prev to src, and src to next. + // if the cross product is negetive, we've turned to the left + // if the cross product is positive, we've turned to the right + // if the cross product is 0, then we need to use the direction passed in + SGVec3d dir1 = prev.toSGVec3d() - cur.toSGVec3d(); + dir1 = normalize(dir1); + + SGVec3d dir2 = next.toSGVec3d() - cur.toSGVec3d(); + dir2 = normalize(dir2); + + // Now find the average + SGVec3d avg = dir1 + dir2; + avg = normalize(avg); + + // check the turn direction + SGVec3d cp = cross( dir1, dir2 ); + SG_LOG(SG_GENERAL, SG_DEBUG, "\tcross product of dir1: " << dir1 << " and dir2: " << dir2 << " is " << cp ); + + // calculate the angle between cur->prev and cur->next + theta = SGMiscd::rad2deg(CalculateTheta(prev, cur, next)); + + if ( abs(theta - 180.0) < 0.1 ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta close to 180) theta is " << theta ); + + // find the direction to the next point + geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); + + offset_dir = next_dir - 90.0; + while (offset_dir < 0.0) + { + offset_dir += 360.0; + } + + // straight line blows up math - dist should be exactly as given + dist = offset_by; + } + else if ( abs(theta) < 0.1 ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta close to 0) : theta is " << theta ); + + // find the direction to the next point + geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); + + offset_dir = next_dir - 90; + while (offset_dir < 0.0) + { + offset_dir += 360.0; + } + + // straight line blows up math - dist should be exactly as given + dist = offset_by; + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta NOT close to 180) : theta is " << theta ); + + // find the offset angle + geo_inverse_wgs_84( avg.y(), avg.x(), 0.0f, 0.0f, &offset_dir, &az2, &dist); + + // if we turned right, reverse the heading + if (cp.z() < 0.0f) + { + turn_dir = 0; + offset_dir += 180.0; + } + else + { + turn_dir = 1; + } + + while (offset_dir >= 360.0) + { + offset_dir -= 360.0; + } + + // find the direction to the next point + geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); + + // calculate correct distance for the offset point + dist = (offset_by)/sin(SGMiscd::deg2rad(next_dir-offset_dir)); + } + + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << dist ); + + // calculate the point from cur + geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, dist, &pt_y, &pt_x, &az2 ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); + + return Point3D(pt_x, pt_y, 0.0f); +} + +Point3D OffsetPointFirst( Point3D cur, Point3D next, double offset_by ) +{ + double offset_dir; + double az2; + double dist; + double pt_x = 0, pt_y = 0; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at Start : cur (" << cur << "), " + "next (" << next << ")" ); + + // find the offset angle + geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &offset_dir, &az2, &dist); + offset_dir -= 90; + if (offset_dir < 0) + { + offset_dir += 360; + } + + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << offset_by ); + + // calculate the point from cur + geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, offset_by, &pt_y, &pt_x, &az2 ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); + + return Point3D(pt_x, pt_y, 0.0f); +} + +Point3D OffsetPointLast( Point3D prev, Point3D cur, double offset_by ) +{ + double offset_dir; + double az2; + double dist; + double pt_x = 0, pt_y = 0; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at End : prev (" << prev << "), " + "cur (" << cur << ")" ); + + // find the offset angle + geo_inverse_wgs_84( prev.y(), prev.x(), cur.y(), cur.x(), &offset_dir, &az2, &dist); + offset_dir -= 90; + if (offset_dir < 0) + { + offset_dir += 360; + } + + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << offset_by ); + + // calculate the point from cur + geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, offset_by, &pt_y, &pt_x, &az2 ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); + + return Point3D(pt_x, pt_y, 0.0f); +} + +Point3D midpoint( Point3D p0, Point3D p1 ) +{ + return Point3D( (p0.x() + p1.x()) / 2, (p0.y() + p1.y()) / 2, (p0.z() + p1.z()) / 2 ); +} + +void +makePolygonsTP (const Line &line, double width, poly_list& polys, texparams_list &tps) +{ + 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 last_end_v; + double heading, az2, dist; + double pt_x, pt_y; + + TGPolygon poly; + TGTexParams tp; + + // generate poly and texparam lists for each line segment + for (i=0; i #include +#include +#include #include "line.hxx" @@ -91,6 +93,7 @@ void makePolygon (const Point3D &p, double width, TGPolygon &polygon); */ void makePolygon (const Line &line, double width, TGPolygon &polygon); +void makePolygonsTP (const Line &line, double width, poly_list &segments, texparams_list &tps); /** * Make a 2D bounding rectangle for a polygon. diff --git a/src/Lib/Polygon/chop-bin.cxx b/src/Lib/Polygon/chop-bin.cxx index e0a90361..f28cf3e1 100644 --- a/src/Lib/Polygon/chop-bin.cxx +++ b/src/Lib/Polygon/chop-bin.cxx @@ -43,9 +43,6 @@ #include "simple_clip.hxx" #include "chop.hxx" -using std::cout; -using std::endl; - static void clip_and_write_poly( string root, long int p_index, const string &poly_type, SGBucket b, const TGPolygon& shape, @@ -97,24 +94,21 @@ static void clip_and_write_poly( string root, long int p_index, fclose(bfp); */ - cout << "shape contours = " << shape.contours() << " "; + SG_LOG(SG_GENERAL, SG_DEBUG, "shape contours = " << shape.contours() ); for ( int ii = 0; ii < shape.contours(); ii++ ) - cout << "hole = " << shape.get_hole_flag(ii) << " "; - cout << endl; + SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << shape.get_hole_flag(ii) ); result = tgPolygonInt( base, shape ); // write_polygon(shape, "shape"); // write_polygon(result, "result"); - cout << "result contours = " << result.contours() << " "; + SG_LOG(SG_GENERAL, SG_DEBUG, "result contours = " << result.contours() ); for ( int ii = 0; ii < result.contours(); ii++ ) { - cout << "hole = " << result.get_hole_flag(ii) << " "; - //if ( result.get_hole_flag(ii) ) { - // exit(0); - //} + SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << result.get_hole_flag(ii) ); + } - cout << endl; + if ( preserve3d ) result.inherit_elevations( shape ); @@ -157,6 +151,108 @@ static void clip_and_write_poly( string root, long int p_index, } } +static void clip_and_write_poly_tp( string root, long int p_index, + const string &poly_type, + SGBucket b, const TGPolygon& shape, + const TGTexParams& tp, 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, "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 ); + + 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 ) { + 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_TP\n" ); + else + fprintf( rfp, "#2D_TP\n" ); + fprintf( rfp, "%s\n", poly_type.c_str() ); + + fprintf( rfp, "%.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", + tp.get_ref().x(), tp.get_ref().y(), tp.get_width(), tp.get_length(), + tp.get_heading(), + tp.get_minu(), tp.get_maxu(), tp.get_minu(), tp.get_maxu()); + + 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 ); + } +} + + // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) void tgChopNormalPolygon( const string& path, const string& poly_type, @@ -206,7 +302,7 @@ void tgChopNormalPolygon( const string& path, const string& poly_type, int dx, dy; sgBucketDiff(b_min, b_max, &dx, &dy); - SG_LOG( SG_GENERAL, SG_INFO, + SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" ); SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); @@ -303,6 +399,155 @@ void tgChopNormalPolygon( const string& path, const string& poly_type, } } +// 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, + const TGPolygon& shape, const TGTexParams& tp, bool preserve3d ) +{ + Point3D min, max, p; + // point2d min, max; + long int index; + int i, j; + + // bail out immediately if polygon is empty + if ( shape.contours() == 0 ) + return; + + min = Point3D( 200.0 ); + max = Point3D( -200.0 ); + + // 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() ); + } + } + + // get next polygon index + index = poly_index_next(); + + 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_poly_tp( path, index, poly_type, b_min, shape, tp, 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_poly_tp( path, index, poly_type, b_cur, shape, tp, + 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; + + bottom.erase(); + bottom_clip.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) ); + + bottom_clip = tgPolygonInt( bottom, shape ); + + // the texparam should be constant over each clipped poly. + // when they are reassembled, we want the texture map to + // be seamless + tgChopNormalPolygonTP( path, poly_type, bottom_clip, tp, 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; + + top.erase(); + top_clip.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) ); + + top_clip = tgPolygonInt( top, shape ); + + tgChopNormalPolygonTP( path, poly_type, top_clip, tp, preserve3d ); + } +} + // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) This has a front end to a crude clipper // that doesn't handle holes so beware. This routine is appropriate diff --git a/src/Lib/Polygon/chop.hxx b/src/Lib/Polygon/chop.hxx index 52eed6f0..98c7bcd1 100644 --- a/src/Lib/Polygon/chop.hxx +++ b/src/Lib/Polygon/chop.hxx @@ -30,13 +30,17 @@ #include #include "polygon.hxx" - +#include "texparams.hxx" // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) void tgChopNormalPolygon( const std::string& path, const std::string& poly_type, const TGPolygon& shape, 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, + const TGPolygon& shape, const TGTexParams& tp, bool preserve3d ); // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) This has a front end to a crude clipper diff --git a/src/Prep/OGRDecode/ogr-decode.cxx b/src/Prep/OGRDecode/ogr-decode.cxx index 7674020c..620a5bfd 100644 --- a/src/Prep/OGRDecode/ogr-decode.cxx +++ b/src/Prep/OGRDecode/ogr-decode.cxx @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -38,12 +39,16 @@ #include #include #include +#include +#include #include -using std:: cout ; -using std:: string ; -using std:: map ; +/* stretch endpoints to reduce slivers in linear data ~.1 meters */ +#define EP_STRETCH (0.000001) + +using std::string; +using std::map; int line_width=50; string line_width_col; @@ -52,6 +57,7 @@ string point_width_col; string area_type="Default"; string area_type_col; int continue_on_errors=0; +int texture_lines = 0; int max_segment_length=0; // ==0 => don't split int start_record=0; bool use_attribute_query=false; @@ -86,8 +92,7 @@ void processLineString(OGRLineString* poGeometry, TGPolygon shape; if (poGeometry->getNumPoints()<2) { - SG_LOG( SG_GENERAL, SG_WARN, - "Skipping line with less than two points" ); + SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" ); return; } @@ -95,14 +100,91 @@ void processLineString(OGRLineString* poGeometry, for (int i=0;igetNumPoints();i++) { line.addPoint(Point3D(poGeometry->getX(i),poGeometry->getY(i),0)); } + tg::makePolygon(line,width,shape); if ( max_segment_length > 0 ) { shape = tgPolygonSplitLongEdges( shape, max_segment_length ); } + tgChopNormalPolygon(work_dir, area_type, shape, false); } +void processLineStringWithTextureInfo(OGRLineString* poGeometry, + const string& work_dir, + const string& area_type, + int width) +{ + poly_list segments; + TGPolygon segment; + texparams_list tps; + TGTexParams tp; + 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;igetX(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; jgetX(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::makePolygonsTP(line,width,segments,tps); + + for ( i = 0; i < (int)segments.size(); ++i ) { + segment = segments[i]; + tp = tps[i]; + + tgChopNormalPolygonTP(work_dir, area_type, segment, tp, false); + } +} + void processPolygon(OGRPolygon* poGeometry, const string& work_dir, const string& area_type) @@ -321,7 +403,11 @@ void processLayer(OGRLayer* poLayer, if (line_width_field!=-1) { width=poFeature->GetFieldAsInteger(line_width_field); } - processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); + if (texture_lines) { + processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); + } else { + processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); + } break; } case wkbMultiLineString: { @@ -332,7 +418,11 @@ void processLayer(OGRLayer* poLayer, } OGRMultiLineString* multils=(OGRMultiLineString*)poGeometry; for (int i=0;igetNumGeometries();i++) { - processLineString((OGRLineString*)(multils->getGeometryRef(i)),work_dir,area_type_name,width); + if (texture_lines) { + processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); + } else { + processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); + } } break; } @@ -349,6 +439,9 @@ void processLayer(OGRLayer* poLayer, } break; } + default: + /* Ignore unhandled objects */ + break; } } @@ -438,6 +531,10 @@ int main( int argc, char **argv ) { area_type_col=argv[2]; argv+=2; argc-=2; + } else if (!strcmp(argv[1],"--texture-lines")) { + argv++; + argc--; + texture_lines=1; } else if (!strcmp(argv[1],"--continue-on-errors")) { argv++; argc--;