From c0555a1270685a9b0008e1bc251a9112eb195f2d Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski <psadrozinski@gmail.com> Date: Thu, 29 Nov 2012 19:24:31 -0500 Subject: [PATCH] convert tg_nodes to SGGeod and SGVecd/f --- src/BuildTiles/Main/tgconstruct.hxx | 8 +- src/BuildTiles/Main/tgconstruct_cleanup.cxx | 2 +- src/BuildTiles/Main/tgconstruct_clip.cxx | 10 +- src/BuildTiles/Main/tgconstruct_elevation.cxx | 58 ++-- src/BuildTiles/Main/tgconstruct_lookup.cxx | 4 +- src/BuildTiles/Main/tgconstruct_math.cxx | 37 ++- src/BuildTiles/Main/tgconstruct_output.cxx | 5 +- src/BuildTiles/Main/tgconstruct_poly.cxx | 12 +- src/BuildTiles/Main/tgconstruct_shared.cxx | 59 ++-- src/BuildTiles/Main/tgconstruct_tesselate.cxx | 6 +- src/Lib/Geometry/poly_cgal.cxx | 8 +- src/Lib/Geometry/poly_extra.cxx | 30 +- src/Lib/Geometry/poly_extra.hxx | 8 + src/Lib/Geometry/poly_support.cxx | 120 +++++++ src/Lib/Geometry/poly_support.hxx | 17 +- src/Lib/Geometry/tg_nodes.cxx | 294 ++++-------------- src/Lib/Geometry/tg_nodes.hxx | 98 +++--- 17 files changed, 376 insertions(+), 400 deletions(-) diff --git a/src/BuildTiles/Main/tgconstruct.hxx b/src/BuildTiles/Main/tgconstruct.hxx index 4435af05..04cae17e 100644 --- a/src/BuildTiles/Main/tgconstruct.hxx +++ b/src/BuildTiles/Main/tgconstruct.hxx @@ -164,7 +164,7 @@ private: void CalcPointNormals( void ); void CalcTextureCoordinates( void ); // Helpers - SGVec3d calc_normal( double area, Point3D p1, Point3D p2, Point3D p3 ); + SGVec3d calc_normal( double area, const SGVec3d& p1, const SGVec3d& p2, const SGVec3d& p3 ); TGPolygon linear_tex_coords( const TGPolygon& tri, const TGTexParams& tp ); TGPolygon area_tex_coords( const TGPolygon& tri ); @@ -173,7 +173,7 @@ private: void AddCustomObjects( void ); // Misc - void calc_normals( point_list& wgs84_nodes, TGSuperPoly& sp ); + void calc_normals( std::vector<SGVec3d>& wgs84_nodes, TGSuperPoly& sp ); double calc_tri_area( int_list& triangle_nodes ); // debug @@ -244,11 +244,11 @@ public: inline TGNodes* get_nodes() { return &nodes; } // node list in geodetic coords (with fixed elevation) - inline point_list get_geod_nodes() const { return nodes.get_geod_nodes(); } + inline void get_geod_nodes( std::vector<SGGeod>& points ) const { nodes.get_geod_nodes( points ); } // normal list (for each point) in cart coords (for smooth // shading) - inline point_list get_point_normals() const { return nodes.get_normals(); } + inline void get_point_normals( std::vector<SGVec3f>& normals ) const { nodes.get_normals( normals ); } // Debug void set_debug( std::string path, std::vector<std::string> area_defs, std::vector<std::string> shape_defs ); diff --git a/src/BuildTiles/Main/tgconstruct_cleanup.cxx b/src/BuildTiles/Main/tgconstruct_cleanup.cxx index 00703ac3..15a6bb8f 100644 --- a/src/BuildTiles/Main/tgconstruct_cleanup.cxx +++ b/src/BuildTiles/Main/tgconstruct_cleanup.cxx @@ -199,7 +199,7 @@ void TGConstruct::AverageEdgeElevations( void ) elevation = elevation / num_elevations; /* Find this node, and update it's elevation */ - int idx = nodes.find( faces.node ); + int idx = nodes.find( faces.node.toSGGeod() ); if (idx != -1) { TGNode node = nodes.get_node( idx ); diff --git a/src/BuildTiles/Main/tgconstruct_clip.cxx b/src/BuildTiles/Main/tgconstruct_clip.cxx index 08f61892..2fd6df4b 100644 --- a/src/BuildTiles/Main/tgconstruct_clip.cxx +++ b/src/BuildTiles/Main/tgconstruct_clip.cxx @@ -54,19 +54,19 @@ bool TGConstruct::ClipLandclassPolys( void ) { p = Point3D(min.x(), min.y(), -9999.0); safety_base.add_node( 0, p ); - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); p = Point3D(max.x(), min.y(), -9999.0); safety_base.add_node( 0, p ); - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); p = Point3D(max.x(), max.y(), -9999.0); safety_base.add_node( 0, p ); - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); p = Point3D(min.x(), max.y(), -9999.0); safety_base.add_node( 0, p ); - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); // set up land mask, we clip most things to this since it is our // best representation of land vs. ocean. If we have other less @@ -274,7 +274,7 @@ bool TGConstruct::ClipLandclassPolys( void ) { for (int con=0; con < poly.contours(); con++) { for (int node = 0; node < poly.contour_size( con ); node++) { // ensure we have all nodes... - nodes.unique_add( poly.get_pt( con, node ) ); + nodes.unique_add( poly.get_pt( con, node ).toSGGeod() ); } } } diff --git a/src/BuildTiles/Main/tgconstruct_elevation.cxx b/src/BuildTiles/Main/tgconstruct_elevation.cxx index 81b456bd..4968677d 100644 --- a/src/BuildTiles/Main/tgconstruct_elevation.cxx +++ b/src/BuildTiles/Main/tgconstruct_elevation.cxx @@ -52,12 +52,12 @@ void TGConstruct::LoadElevationArray( bool add_nodes ) { if ( add_nodes ) { point_list corner_list = array.get_corner_list(); for (unsigned int i=0; i<corner_list.size(); i++) { - nodes.unique_add(corner_list[i]); + nodes.unique_add( corner_list[i].toSGGeod() ); } point_list fit_list = array.get_fitted_list(); for (unsigned int i=0; i<fit_list.size(); i++) { - nodes.unique_add(fit_list[i]); + nodes.unique_add( fit_list[i].toSGGeod() ); } } } @@ -72,13 +72,13 @@ static void calc_gc_course_dist( const Point3D& start, const Point3D& dest, // calculate spherical distance between two points (lon, lat specified // in degrees, result returned in meters) -static double distanceSphere( const Point3D p1, const Point3D p2 ) { - Point3D r1( p1.x() * SG_DEGREES_TO_RADIANS, - p1.y() * SG_DEGREES_TO_RADIANS, - p1.z() ); - Point3D r2( p2.x() * SG_DEGREES_TO_RADIANS, - p2.y() * SG_DEGREES_TO_RADIANS, - p2.z() ); +static double distanceSphere( const SGGeoc& p1, const SGGeod& p2 ) { + Point3D r1( p1.getLongitudeRad(), + p1.getLatitudeRad(), + p1.getRadiusM() ); + Point3D r2( p2.getLongitudeRad(), + p2.getLatitudeRad(), + p2.getElevationM() ); double course, dist_m; calc_gc_course_dist( r1, r2, &course, &dist_m ); @@ -94,22 +94,22 @@ void TGConstruct::CalcElevations( void ) TGPolyNodes tri_nodes; double e1, e2, e3, min; int n1, n2, n3; - point_list raw_nodes; - Point3D p; + std::vector<SGGeod> raw_nodes; + SGGeoc p; SG_LOG(SG_GENERAL, SG_ALERT, "fixing node heights"); for (int i = 0; i < (int)nodes.size(); ++i) { TGNode node = nodes.get_node( i ); - Point3D pos = node.GetPosition(); + SGGeod pos = node.GetPosition(); if ( !node.GetFixedPosition() ) { // set elevation as interpolated point from DEM data. - nodes.SetElevation( i, array.altitude_from_grid(pos.x() * 3600.0, pos.y() * 3600.0) ); + nodes.SetElevation( i, array.altitude_from_grid(pos.getLongitudeDeg() * 3600.0, pos.getLatitudeDeg() * 3600.0) ); } } - raw_nodes = nodes.get_geod_nodes(); + nodes.get_geod_nodes(raw_nodes); // now flatten some stuff for (unsigned int area = 0; area < TG_MAX_AREA_TYPES; area++) { if ( is_lake_area( (AreaType)area ) ) { @@ -126,11 +126,11 @@ void TGConstruct::CalcElevations( void ) } n1 = tri_nodes.get_pt( tri, 0 ); - e1 = nodes.get_node(n1).GetPosition().z(); + e1 = nodes.get_node(n1).GetPosition().getElevationM(); n2 = tri_nodes.get_pt( tri, 1 ); - e2 = nodes.get_node(n2).GetPosition().z(); + e2 = nodes.get_node(n2).GetPosition().getElevationM(); n3 = tri_nodes.get_pt( tri, 2 ); - e3 = nodes.get_node(n3).GetPosition().z(); + e3 = nodes.get_node(n3).GetPosition().getElevationM(); min = e1; if ( e2 < min ) { min = e2; } @@ -159,17 +159,17 @@ void TGConstruct::CalcElevations( void ) n1 = tri_nodes.get_pt( tri, 0 ); - e1 = nodes.get_node(n1).GetPosition().z(); + e1 = nodes.get_node(n1).GetPosition().getElevationM(); n2 = tri_nodes.get_pt( tri, 1 ); - e2 = nodes.get_node(n2).GetPosition().z(); + e2 = nodes.get_node(n2).GetPosition().getElevationM(); n3 = tri_nodes.get_pt( tri, 2 ); - e3 = nodes.get_node(n3).GetPosition().z(); + e3 = nodes.get_node(n3).GetPosition().getElevationM(); min = e1; - p = raw_nodes[n1]; + p = SGGeoc::fromGeod( raw_nodes[n1] ); - if ( e2 < min ) { min = e2; p = raw_nodes[n2]; } - if ( e3 < min ) { min = e3; p = raw_nodes[n3]; } + if ( e2 < min ) { min = e2; p = SGGeoc::fromGeod( raw_nodes[n2] ); } + if ( e3 < min ) { min = e3; p = SGGeoc::fromGeod( raw_nodes[n3] ); } double d1 = distanceSphere( p, raw_nodes[n1] ); double d2 = distanceSphere( p, raw_nodes[n2] ); @@ -202,17 +202,17 @@ void TGConstruct::CalcElevations( void ) n1 = tri_nodes.get_pt( tri, 0 ); - e1 = nodes.get_node(n1).GetPosition().z(); + e1 = nodes.get_node(n1).GetPosition().getElevationM(); n2 = tri_nodes.get_pt( tri, 1 ); - e2 = nodes.get_node(n2).GetPosition().z(); + e2 = nodes.get_node(n2).GetPosition().getElevationM(); n3 = tri_nodes.get_pt( tri, 2 ); - e3 = nodes.get_node(n3).GetPosition().z(); + e3 = nodes.get_node(n3).GetPosition().getElevationM(); min = e1; - p = raw_nodes[n1]; + p = SGGeoc::fromGeod( raw_nodes[n1] ); - if ( e2 < min ) { min = e2; p = raw_nodes[n2]; } - if ( e3 < min ) { min = e3; p = raw_nodes[n3]; } + if ( e2 < min ) { min = e2; p = SGGeoc::fromGeod( raw_nodes[n2] ); } + if ( e3 < min ) { min = e3; p = SGGeoc::fromGeod( raw_nodes[n3] ); } double d1 = distanceSphere( p, raw_nodes[n1] ); double d2 = distanceSphere( p, raw_nodes[n2] ); diff --git a/src/BuildTiles/Main/tgconstruct_lookup.cxx b/src/BuildTiles/Main/tgconstruct_lookup.cxx index 17efb617..4cf6b493 100644 --- a/src/BuildTiles/Main/tgconstruct_lookup.cxx +++ b/src/BuildTiles/Main/tgconstruct_lookup.cxx @@ -46,7 +46,7 @@ void TGConstruct::LookupNodesPerVertex( void ) for (int tri=0; tri < tris.contours(); tri++) { for (int vertex = 0; vertex < tris.contour_size(tri); vertex++) { - idx = nodes.find( tris.get_pt( tri, vertex ) ); + idx = nodes.find( tris.get_pt( tri, vertex ).toSGGeod() ); if (idx >= 0) { tri_nodes.add_node( tri, idx ); } else { @@ -73,7 +73,7 @@ void TGConstruct::LookupFacesPerNode( void ) for (int tri=0; tri < tris.contours(); tri++) { for (int sub = 0; sub < tris.contour_size(tri); sub++) { - int n = nodes.find( tris.get_pt( tri, sub ) ); + int n = nodes.find( tris.get_pt( tri, sub ).toSGGeod() ); nodes.AddFace( n, area, shape, segment, tri ); } } diff --git a/src/BuildTiles/Main/tgconstruct_math.cxx b/src/BuildTiles/Main/tgconstruct_math.cxx index c3773068..481c4907 100644 --- a/src/BuildTiles/Main/tgconstruct_math.cxx +++ b/src/BuildTiles/Main/tgconstruct_math.cxx @@ -32,14 +32,14 @@ //using std::string; double TGConstruct::calc_tri_area( int_list& triangle_nodes ) { - Point3D p1 = nodes.get_node( triangle_nodes[0] ).GetPosition(); - Point3D p2 = nodes.get_node( triangle_nodes[1] ).GetPosition(); - Point3D p3 = nodes.get_node( triangle_nodes[2] ).GetPosition(); + SGGeod p1 = nodes.get_node( triangle_nodes[0] ).GetPosition(); + SGGeod p2 = nodes.get_node( triangle_nodes[1] ).GetPosition(); + SGGeod p3 = nodes.get_node( triangle_nodes[2] ).GetPosition(); return triangle_area( p1, p2, p3 ); } -SGVec3d TGConstruct::calc_normal( double area, Point3D p1, Point3D p2, Point3D p3 ) { +SGVec3d TGConstruct::calc_normal( double area, const SGVec3d& p1, const SGVec3d& p2, const SGVec3d& p3 ) { SGVec3d v1, v2, normal; // do some sanity checking. With the introduction of landuse @@ -76,7 +76,7 @@ SGVec3d TGConstruct::calc_normal( double area, Point3D p1, Point3D p2, Point3D p return normal; } -void TGConstruct::calc_normals( point_list& wgs84_nodes, TGSuperPoly& sp ) { +void TGConstruct::calc_normals( std::vector<SGVec3d>& wgs84_nodes, TGSuperPoly& sp ) { // for each face in the superpoly, calculate a face normal SGVec3d normal; TGPolyNodes tri_nodes = sp.get_tri_idxs(); @@ -91,12 +91,11 @@ void TGConstruct::calc_normals( point_list& wgs84_nodes, TGSuperPoly& sp ) { for (int i=0; i<tri_nodes.contours(); i++) { face_nodes = tri_nodes.get_contour(i); - Point3D p1 = wgs84_nodes[ face_nodes[0] ]; - Point3D p2 = wgs84_nodes[ face_nodes[1] ]; - Point3D p3 = wgs84_nodes[ face_nodes[2] ]; + SGVec3d p1 = wgs84_nodes[ face_nodes[0] ]; + SGVec3d p2 = wgs84_nodes[ face_nodes[1] ]; + SGVec3d p3 = wgs84_nodes[ face_nodes[2] ]; area = calc_tri_area( face_nodes ); - normal = calc_normal( area, p1, p2, p3 ); face_normals.push_back( Point3D::fromSGVec3( normal ) ); @@ -110,7 +109,8 @@ void TGConstruct::calc_normals( point_list& wgs84_nodes, TGSuperPoly& sp ) { void TGConstruct::CalcFaceNormals( void ) { // traverse the superpols, and calc normals for each tri within - point_list wgs84_nodes = nodes.get_wgs84_nodes_as_Point3d(); + std::vector<SGVec3d> wgs84_nodes; + nodes.get_wgs84_nodes( wgs84_nodes ); for (unsigned int area = 0; area < TG_MAX_AREA_TYPES; area++) { for (unsigned int shape = 0; shape < polys_clipped.area_size(area); shape++ ) { @@ -127,10 +127,12 @@ void TGConstruct::CalcPointNormals( void ) // traverse triangle structure building the face normal table SG_LOG(SG_GENERAL, SG_ALERT, "Calculating point normals: 0%"); - Point3D normal; + SGVec3d normal; double face_area; - point_list wgs84_nodes = nodes.get_wgs84_nodes_as_Point3d(); + std::vector<SGVec3d> wgs84_nodes; + nodes.get_wgs84_nodes( wgs84_nodes ); + unsigned int one_percent = nodes.size() / 100; unsigned int cur_percent = 1; @@ -140,7 +142,7 @@ void TGConstruct::CalcPointNormals( void ) TGNeighborFaces* neighbor_faces = NULL; double total_area = 0.0; - Point3D average( 0.0 ); + SGVec3d average( 0.0, 0.0, 0.0 ); if ( i == one_percent ) { SG_LOG(SG_GENERAL, SG_ALERT, "Calculating point normals: " << cur_percent << "%" ); @@ -156,7 +158,7 @@ void TGConstruct::CalcPointNormals( void ) unsigned int tri = faces[j].tri; int_list face_nodes; - normal = polys_clipped.get_face_normal( at, shape, segment, tri ); + normal = polys_clipped.get_face_normal( at, shape, segment, tri ).toSGVec3d(); face_nodes = polys_clipped.get_tri_idxs( at, shape, segment ).get_contour( tri ) ; face_area = polys_clipped.get_face_area( at, shape, segment, tri ); @@ -166,11 +168,11 @@ void TGConstruct::CalcPointNormals( void ) } // if this node exists in the shared edge db, add the faces from the neighbooring tile - neighbor_faces = FindNeighborFaces( node.GetPosition() ); + neighbor_faces = FindNeighborFaces( Point3D::fromSGGeod( node.GetPosition() ) ); if ( neighbor_faces ) { int num_faces = neighbor_faces->face_areas.size(); for ( int j = 0; j < num_faces; j++ ) { - normal = neighbor_faces->face_normals[j]; + normal = neighbor_faces->face_normals[j].toSGVec3d(); face_area = neighbor_faces->face_areas[j]; normal *= face_area; @@ -180,6 +182,7 @@ void TGConstruct::CalcPointNormals( void ) } average /= total_area; - nodes.SetNormal( i, average ); + SGVec3f n = SGVec3f( average.x(), average.y(), average.z() ); + nodes.SetNormal( i, n ); } } diff --git a/src/BuildTiles/Main/tgconstruct_output.cxx b/src/BuildTiles/Main/tgconstruct_output.cxx index a87ce58f..4b27acf2 100644 --- a/src/BuildTiles/Main/tgconstruct_output.cxx +++ b/src/BuildTiles/Main/tgconstruct_output.cxx @@ -155,7 +155,7 @@ void TGConstruct::WriteBtgFile( void ) tri_v.push_back( index ); // add the node's normal - index = normals.unique_add( nodes.GetNormal( index ) ); + index = normals.unique_add( Point3D::fromSGVec3( nodes.GetNormal( index ) ) ); tri_n.push_back( index ); Point3D tc = tri_txs.get_pt( k, l ); @@ -174,7 +174,8 @@ void TGConstruct::WriteBtgFile( void ) } } - std::vector< SGVec3d > wgs84_nodes = nodes.get_wgs84_nodes_as_SGVec3d(); + std::vector< SGVec3d > wgs84_nodes; + nodes.get_wgs84_nodes( wgs84_nodes ); SGVec3d gbs_center = SGVec3d::fromGeod( bucket.get_center() ); double dist_squared, radius_squared = 0; for (int i = 0; i < (int)wgs84_nodes.size(); ++i) diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index d1a98dbc..9dbf4cac 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -177,9 +177,9 @@ bool TGConstruct::load_poly(const string& path) { poly.add_node( i, p ); if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); + nodes.unique_add_fixed_elevation( p.toSGGeod() ); } else { - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); } for ( j = 1; j < count - 1; ++j ) { @@ -194,9 +194,9 @@ bool TGConstruct::load_poly(const string& path) { p.snap( gSnap ); poly.add_node( i, p ); if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); + nodes.unique_add_fixed_elevation( p.toSGGeod() ); } else { - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); } } @@ -217,9 +217,9 @@ bool TGConstruct::load_poly(const string& path) { p.snap( gSnap ); poly.add_node( i, p ); if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); + nodes.unique_add_fixed_elevation( p.toSGGeod() ); } else { - nodes.unique_add( p ); + nodes.unique_add( p.toSGGeod() ); } } } diff --git a/src/BuildTiles/Main/tgconstruct_shared.cxx b/src/BuildTiles/Main/tgconstruct_shared.cxx index dfa2182a..ee591545 100644 --- a/src/BuildTiles/Main/tgconstruct_shared.cxx +++ b/src/BuildTiles/Main/tgconstruct_shared.cxx @@ -43,7 +43,7 @@ void TGConstruct::SaveSharedEdgeData( int stage ) switch( stage ) { case 1: { - point_list north, south, east, west; + std::vector<SGGeod> north, south, east, west; int nCount; nodes.get_geod_edge( bucket, north, south, east, west ); @@ -64,28 +64,28 @@ void TGConstruct::SaveSharedEdgeData( int stage ) nCount = north.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, north[i] ); + sgWriteGeod( fp, north[i] ); } // south nCount = south.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, south[i] ); + sgWriteGeod( fp, south[i] ); } // east nCount = east.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, east[i] ); + sgWriteGeod( fp, east[i] ); } // west nCount = west.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, west[i] ); + sgWriteGeod( fp, west[i] ); } gzclose(fp); @@ -113,7 +113,7 @@ void TGConstruct::SaveSharedEdgeData( int stage ) void TGConstruct::WriteNeighborFaces( gzFile& fp, Point3D pt ) { // find all neighboors of this point - int n = nodes.find( pt ); + int n = nodes.find( pt.toSGGeod() ); TGNode node = nodes.get_node( n ); TGFaceList faces = node.GetFaces(); @@ -130,19 +130,19 @@ void TGConstruct::WriteNeighborFaces( gzFile& fp, Point3D pt ) int_list face_nodes = polys_clipped.get_tri_idxs( at, shape, segment ).get_contour( tri ) ; { - Point3D p1 = nodes.get_node( face_nodes[0] ).GetPosition(); - Point3D p2 = nodes.get_node( face_nodes[1] ).GetPosition(); - Point3D p3 = nodes.get_node( face_nodes[2] ).GetPosition(); + SGGeod p1 = nodes.get_node( face_nodes[0] ).GetPosition(); + SGGeod p2 = nodes.get_node( face_nodes[1] ).GetPosition(); + SGGeod p3 = nodes.get_node( face_nodes[2] ).GetPosition(); - Point3D wgs_p1 = nodes.get_node( face_nodes[0] ).GetWgs84AsPoint3D(); - Point3D wgs_p2 = nodes.get_node( face_nodes[1] ).GetWgs84AsPoint3D(); - Point3D wgs_p3 = nodes.get_node( face_nodes[2] ).GetWgs84AsPoint3D(); + SGVec3d wgs_p1 = nodes.get_node( face_nodes[0] ).GetWgs84(); + SGVec3d wgs_p2 = nodes.get_node( face_nodes[1] ).GetWgs84(); + SGVec3d wgs_p3 = nodes.get_node( face_nodes[2] ).GetWgs84(); double face_area = triangle_area( p1, p2, p3 ); - Point3D face_normal = Point3D::fromSGVec3( calc_normal( face_area, wgs_p1, wgs_p2, wgs_p3 ) ); + SGVec3d face_normal = calc_normal( face_area, wgs_p1, wgs_p2, wgs_p3 ); sgWriteDouble( fp, face_area ); - sgWritePoint3D( fp, face_normal ); + sgWritedVec3( fp, face_normal ); } } } @@ -193,10 +193,10 @@ void TGConstruct::ReadNeighborFaces( gzFile& fp ) pFaces = AddNeighborFaces( node ); // new face - let's add our elevation first - int idx = nodes.find( node ); + int idx = nodes.find( node.toSGGeod() ); if (idx >= 0) { TGNode local = nodes.get_node( idx ); - pFaces->elevations.push_back( local.GetPosition().z() ); + pFaces->elevations.push_back( local.GetPosition().getElevationM() ); } } @@ -222,7 +222,7 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) { string dir; string file; - point_list north, south, east, west; + std::vector<SGGeod> north, south, east, west; int nCount; nodes.get_geod_edge( bucket, north, south, east, west ); @@ -233,7 +233,6 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) sgp.append( "dummy" ); sgp.create_dir( 0755 ); - // north edge file = dir + "/" + bucket.gen_index_str() + "_north_edge"; gzFile fp; @@ -247,8 +246,8 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { // write the 3d point - sgWritePoint3D( fp, north[i] ); - WriteNeighborFaces( fp, north[i] ); + sgWriteGeod( fp, north[i] ); + WriteNeighborFaces( fp, Point3D::fromSGGeod( north[i] ) ); } gzclose(fp); @@ -263,8 +262,8 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) nCount = south.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, south[i] ); - WriteNeighborFaces( fp, south[i] ); + sgWriteGeod( fp, south[i] ); + WriteNeighborFaces( fp, Point3D::fromSGGeod( south[i] ) ); } gzclose(fp); @@ -279,8 +278,8 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) nCount = east.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, east[i] ); - WriteNeighborFaces( fp, east[i] ); + sgWriteGeod( fp, east[i] ); + WriteNeighborFaces( fp, Point3D::fromSGGeod( east[i] ) ); } gzclose(fp); @@ -295,8 +294,8 @@ void TGConstruct::SaveSharedEdgeDataStage2( void ) nCount = west.size(); sgWriteInt( fp, nCount ); for (int i=0; i<nCount; i++) { - sgWritePoint3D( fp, west[i] ); - WriteNeighborFaces( fp, west[i] ); + sgWriteGeod( fp, west[i] ); + WriteNeighborFaces( fp, Point3D::fromSGGeod( west[i] ) ); } gzclose(fp); } @@ -497,28 +496,28 @@ void TGConstruct::LoadSharedEdgeData( int stage ) LoadNeighboorEdgeDataStage1( nb, north, south, east, west ); // Add southern nodes from northern tile for (unsigned int i=0; i<south.size(); i++) { - nodes.unique_add( south[i] ); + nodes.unique_add( south[i].toSGGeod() ); } // Read South Tile and add its northern nodes sb = sgBucketOffset(clon, clat, 0, -1); LoadNeighboorEdgeDataStage1( sb, north, south, east, west ); for (unsigned int i=0; i<north.size(); i++) { - nodes.unique_add( north[i] ); + nodes.unique_add( north[i].toSGGeod() ); } // Read East Tile and add its western nodes eb = sgBucketOffset(clon, clat, 1, 0); LoadNeighboorEdgeDataStage1( eb, north, south, east, west ); for (unsigned int i=0; i<west.size(); i++) { - nodes.unique_add( west[i] ); + nodes.unique_add( west[i].toSGGeod() ); } // Read West Tile and add its eastern nodes wb = sgBucketOffset(clon, clat, -1, 0); LoadNeighboorEdgeDataStage1( wb, north, south, east, west ); for (unsigned int i=0; i<east.size(); i++) { - nodes.unique_add( east[i] ); + nodes.unique_add( east[i].toSGGeod() ); } } break; diff --git a/src/BuildTiles/Main/tgconstruct_tesselate.cxx b/src/BuildTiles/Main/tgconstruct_tesselate.cxx index 3218564c..676f0da5 100644 --- a/src/BuildTiles/Main/tgconstruct_tesselate.cxx +++ b/src/BuildTiles/Main/tgconstruct_tesselate.cxx @@ -38,7 +38,7 @@ void TGConstruct::TesselatePolys( void ) { // tesselate the polygons and prepair them for final output - point_list poly_extra; + std::vector<SGGeod> poly_extra; SGGeod min, max; for (unsigned int area = 0; area < TG_MAX_AREA_TYPES; area++) { @@ -53,7 +53,7 @@ void TGConstruct::TesselatePolys( void ) TGPolygon poly = polys_clipped.get_poly(area, shape, segment); poly.get_bounding_box(min, max); - poly_extra = nodes.get_geod_inside( Point3D::fromSGGeod(min), Point3D::fromSGGeod(max) ); + nodes.get_geod_inside( min, max, poly_extra ); SG_LOG( SG_CLIPPER, SG_INFO, "Tesselating " << get_area_name( (AreaType)area ) << "(" << area << "): " << shape+1 << "-" << segment << " of " << (int)polys_clipped.area_size(area) << @@ -69,7 +69,7 @@ void TGConstruct::TesselatePolys( void ) for (int k=0; k< tri.contours(); k++) { for (int l = 0; l < tri.contour_size(k); l++) { // ensure we have all nodes... - nodes.unique_add( tri.get_pt( k, l ) ); + nodes.unique_add( tri.get_pt( k, l ).toSGGeod() ); } } diff --git a/src/Lib/Geometry/poly_cgal.cxx b/src/Lib/Geometry/poly_cgal.cxx index 19633fe1..a108927a 100644 --- a/src/Lib/Geometry/poly_cgal.cxx +++ b/src/Lib/Geometry/poly_cgal.cxx @@ -103,7 +103,7 @@ void insert_polygon(CDT& cdt,const Polygon_2& polygon) } } -TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const point_list& extra_nodes, bool verbose ) { +TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const std::vector<SGGeod>& extra_nodes, bool verbose ) { TGPolygon result; CDT cdt; @@ -118,10 +118,10 @@ TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const point_list& std::vector<Point> points; points.reserve(extra_nodes.size()); for (unsigned int n = 0; n < extra_nodes.size(); n++) { - points.push_back( Point(extra_nodes[n].x(), extra_nodes[n].y()) ); + points.push_back( Point(extra_nodes[n].getLongitudeDeg(), extra_nodes[n].getLatitudeDeg()) ); } cdt.insert(points.begin(), points.end()); - + // then insert each polygon as a constraint into the triangulation for (int c = 0; c < p.contours(); c++) { point_list contour = p.get_contour( c ); @@ -160,6 +160,6 @@ TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const point_list& TGPolygon polygon_tesselate_alt_cgal( TGPolygon &p, bool verbose ) { - point_list pl; pl.clear(); + std::vector<SGGeod> pl; pl.clear(); return ( polygon_tesselate_alt_with_extra_cgal(p, pl, verbose) ); } diff --git a/src/Lib/Geometry/poly_extra.cxx b/src/Lib/Geometry/poly_extra.cxx index eb09bcdf..c6297fff 100644 --- a/src/Lib/Geometry/poly_extra.cxx +++ b/src/Lib/Geometry/poly_extra.cxx @@ -59,6 +59,32 @@ void add_intermediate_nodes( int contour, const Point3D& start, } } +void add_intermediate_nodes( int contour, const Point3D& start, + const Point3D& end, std::vector<SGGeod>&tmp_nodes, + TGPolygon *result, + const double bbEpsilon, + const double errEpsilon + ) +{ + Point3D new_pt; + + SG_LOG(SG_GENERAL, SG_BULK, " " << start << " <==> " << end ); + + bool found_extra = find_intermediate_node( start, end, tmp_nodes, &new_pt, bbEpsilon, errEpsilon ); + + if ( found_extra ) { + // recurse with two sub segments + // SG_LOG(SG_GENERAL, SG_DEBUG, "dividing " << p0 << " " << nodes[extra_index] + // << " " << p1); + add_intermediate_nodes( contour, start, new_pt, tmp_nodes, result, bbEpsilon, errEpsilon ); + + result->add_node( contour, new_pt ); + SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt); + + add_intermediate_nodes( contour, new_pt, end, tmp_nodes, result, bbEpsilon, errEpsilon ); + } +} + // Search each segment for additional vertex points that may have been // created elsewhere that lie on the segment and split it there to @@ -113,12 +139,12 @@ TGPolygon add_tgnodes_to_poly( const TGPolygon& poly, TGPolygon result; result.erase(); SGGeod min, max; Point3D p0, p1; - point_list poly_points; + std::vector<SGGeod> poly_points; poly.get_bounding_box(min, max); SG_LOG(SG_GENERAL, SG_DEBUG, "add_tgnodes_to_poly : min " << min << " max " << max ); - poly_points = nodes->get_geod_inside( Point3D::fromSGGeod(min), Point3D::fromSGGeod(max) ); + nodes->get_geod_inside( min, max, poly_points ); for ( int i = 0; i < poly.contours(); ++i ) { SG_LOG(SG_GENERAL, SG_DEBUG, "contour = " << i); diff --git a/src/Lib/Geometry/poly_extra.hxx b/src/Lib/Geometry/poly_extra.hxx index 580151de..ba6dbec7 100644 --- a/src/Lib/Geometry/poly_extra.hxx +++ b/src/Lib/Geometry/poly_extra.hxx @@ -44,6 +44,14 @@ void add_intermediate_nodes( int contour, const Point3D& start, ); +// TEMP - converting tgnodes to SGGeod/SGVec3d +void add_intermediate_nodes( int contour, const Point3D& start, + const Point3D& end, const std::vector<SGGeod>& tmp_nodes, + TGPolygon *result, + const double bbEpsilon = SG_EPSILON*10, + const double errEpsilon = SG_EPSILON*4 + ); + // Search each segment for additional vertex points that may have been // created elsewhere that lie on the segment and split it there to // avoid "T" intersections. diff --git a/src/Lib/Geometry/poly_support.cxx b/src/Lib/Geometry/poly_support.cxx index 31a3cfa8..72c0e0f6 100644 --- a/src/Lib/Geometry/poly_support.cxx +++ b/src/Lib/Geometry/poly_support.cxx @@ -730,6 +730,126 @@ bool find_intermediate_node( const Point3D& start, const Point3D& end, return found_node; } +// temp +bool find_intermediate_node( const Point3D& start, const Point3D& end, + const std::vector<SGGeod>& nodes, Point3D *result, + double bbEpsilon, double errEpsilon + ) +{ + bool found_node = false; + double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min; + + Point3D p0 = start; + Point3D p1 = end; + + // cout << " find_intermediate_nodes() " << p0 << " <=> " << p1 << endl; + + double xdist = fabs(p0.x() - p1.x()); + double ydist = fabs(p0.y() - p1.y()); + // cout << "xdist = " << xdist << " ydist = " << ydist << endl; + x_err_min = xdist + 1.0; + y_err_min = ydist + 1.0; + + if ( xdist > ydist ) { + // cout << "use y = mx + b" << endl; + + // sort these in a sensible order + Point3D p_min, p_max; + if ( p0.x() < p1.x() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m = (p_min.y() - p_max.y()) / (p_min.x() - p_max.x()); + b = p_max.y() - m * p_max.x(); + + // cout << "m = " << m << " b = " << b << endl; + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + // cout << i << endl; + Point3D current = Point3D::fromSGGeod( nodes[i] ); + + if ( (current.x() > (p_min.x() + (bbEpsilon))) + && (current.x() < (p_max.x() - (bbEpsilon))) ) { + + // printf( "found a potential candidate %.7f %.7f %.7f\n", + // current.x(), current.y(), current.z() ); + + y_err = fabs(current.y() - (m * current.x() + b)); + // cout << "y_err = " << y_err << endl; + + if ( y_err < errEpsilon ) { + // cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; + // cout << p_min << " < " << current << " < " + // << p_max << endl; + found_node = true; + if ( y_err < y_err_min ) { + *result = current; + y_err_min = y_err; + } + } + } + } + } else { + // cout << "use x = m1 * y + b1" << endl; + + // sort these in a sensible order + Point3D p_min, p_max; + if ( p0.y() < p1.y() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m1 = (p_min.x() - p_max.x()) / (p_min.y() - p_max.y()); + b1 = p_max.x() - m1 * p_max.y(); + + // cout << " m1 = " << m1 << " b1 = " << b1 << endl; + // printf( " m = %.8f b = %.8f\n", 1/m1, -b1/m1); + + // cout << " should = 0 = " + // << fabs(p_min.x() - (m1 * p_min.y() + b1)) << endl; + // cout << " should = 0 = " + // << fabs(p_max.x() - (m1 * p_max.y() + b1)) << endl; + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + Point3D current = Point3D::fromSGGeod( nodes[i] ); + + if ( (current.y() > (p_min.y() + (bbEpsilon))) + && (current.y() < (p_max.y() - (bbEpsilon))) ) { + + // printf( "found a potential candidate %.7f %.7f %.7f\n", + // current.x(), current.y(), current.z() ); + + x_err = fabs(current.x() - (m1 * current.y() + b1)); + // cout << "x_err = " << x_err << endl; + + // if ( temp ) { + // cout << " (" << counter << ") x_err = " << x_err << endl; + // } + + if ( x_err < errEpsilon ) { + // cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; + // cout << p_min << " < " << current << " < " + // << p_max << endl; + found_node = true; + if ( x_err < x_err_min ) { + *result = current; + x_err_min = x_err; + } + } + } + } + } + + return found_node; +} + // Attempt to reduce degeneracies where a subsequent point of a // polygon lies *on* a previous line segment. diff --git a/src/Lib/Geometry/poly_support.hxx b/src/Lib/Geometry/poly_support.hxx index f0e66ea8..cacb551b 100644 --- a/src/Lib/Geometry/poly_support.hxx +++ b/src/Lib/Geometry/poly_support.hxx @@ -41,11 +41,11 @@ // Calculate the area of a triangle -inline double triangle_area( const Point3D& p1, const Point3D& p2, const Point3D& p3 ) +inline double triangle_area( const SGGeod& p1, const SGGeod& p2, const SGGeod& p3 ) { - return fabs(0.5 * ( p1.x() * p2.y() - p2.x() * p1.y() + - p2.x() * p3.y() - p3.x() * p2.y() + - p3.x() * p1.y() - p1.x() * p3.y() )); + return fabs(0.5 * ( p1.getLongitudeDeg() * p2.getLatitudeDeg() - p2.getLongitudeDeg() * p1.getLatitudeDeg() + + p2.getLongitudeDeg() * p3.getLatitudeDeg() - p3.getLongitudeDeg() * p2.getLatitudeDeg() + + p3.getLongitudeDeg() * p1.getLatitudeDeg() - p1.getLongitudeDeg() * p3.getLatitudeDeg() )); } @@ -53,7 +53,7 @@ inline double triangle_area( const Point3D& p1, const Point3D& p2, const Point3D // or splitting edges and without regard for holes. Returns a polygon // with one contour per tesselated triangle. TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, - const point_list &extra_nodes, bool verbose ); + const std::vector<SGGeod>& extra_nodes, bool verbose ); TGPolygon polygon_tesselate_alt_cgal( TGPolygon &p, bool verbose ); @@ -93,6 +93,13 @@ bool find_intermediate_node( const Point3D& start, const Point3D& end, const double errEpsilon = SG_EPSILON*4 ); +// TEMP +bool find_intermediate_node( const Point3D& start, const Point3D& end, + const std::vector<SGGeod>& nodes, Point3D *result, + const double bbEpsilon = SG_EPSILON*10, + const double errEpsilon = SG_EPSILON*4 + ); + // remove any degenerate contours TGPolygon remove_bad_contours( const TGPolygon &poly ); diff --git a/src/Lib/Geometry/tg_nodes.cxx b/src/Lib/Geometry/tg_nodes.cxx index ad1f4371..0de198b4 100644 --- a/src/Lib/Geometry/tg_nodes.cxx +++ b/src/Lib/Geometry/tg_nodes.cxx @@ -6,22 +6,22 @@ // compare node's positions (x, then y) int compare_position(const TGNode& n1, const TGNode& n2) { - Point3D pos1 = n1.GetPosition(); - Point3D pos2 = n2.GetPosition(); + SGGeod pos1 = n1.GetPosition(); + SGGeod pos2 = n2.GetPosition(); - if ( pos1.x() == pos2.x() ) { - return ( pos1.y() < pos2.y() ); + if ( pos1.getLongitudeDeg() == pos2.getLongitudeDeg() ) { + return ( pos1.getLatitudeDeg() < pos2.getLatitudeDeg() ); } else { - return ( pos1.x() < pos2.x() ); + return ( pos1.getLongitudeDeg() < pos2.getLongitudeDeg() ); } } int fuzzy_compare_xposition(const TGNode& n1, const TGNode& n2) { - Point3D pos1 = n1.GetPosition(); - Point3D pos2 = n2.GetPosition(); + SGGeod pos1 = n1.GetPosition(); + SGGeod pos2 = n2.GetPosition(); - if ( fabs(pos1.x() - pos2.x()) < FG_PROXIMITY_EPSILON ) { + if ( fabs(pos1.getLongitudeDeg() - pos2.getLongitudeDeg()) < FG_PROXIMITY_EPSILON ) { /* if x coords are within vacinity, then pos1 < pos2 */ return 1; } else { @@ -36,7 +36,7 @@ void TGNodes::SortNodes( void ) } // Find the index of the specified point (compair to the same // tolerance as unique_add(). Returns -1 if not found. -int TGNodes::sorted_find( const Point3D& p ) const { +int TGNodes::sorted_find( const SGGeod& p ) const { TGNode node( p ); const_node_list_iterator lb, ub; @@ -56,9 +56,9 @@ int TGNodes::sorted_find( const Point3D& p ) const { return -1; } -int TGNodes::linear_find( const Point3D& p ) const { +int TGNodes::linear_find( const SGGeod& p ) const { const_node_list_iterator current, last; - Point3D pos; + SGGeod pos; int counter = 0; // see if point already exists @@ -76,44 +76,7 @@ int TGNodes::linear_find( const Point3D& p ) const { return -1; } -#if 0 -void TGNodes::sorted_unique_add( const Point3D& p ) { -# if 0 - TGNode node( p ); - node.SetFixedPosition( false ); - - node_list_iterator idx = std::lower_bound( tg_node_list.begin(), tg_node_list.end(), node, compare_position ); - tg_node_list.insert( idx, node ); -#else - TGNode node( p ); - node.SetFixedPosition( false ); - - node_list_iterator lb, ub, idx; - - // first, find the range to search - ub=lower_bound( tg_node_list.begin(), tg_node_list.end(), node, fuzzy_compare_xposition ); - lb=upper_bound( tg_node_list.begin(), tg_node_list.end(), node, fuzzy_compare_xposition ); - - // then do a normal linear search in the range - if ( lb != tg_node_list.end() ) { - for ( ; lb != ub; ++lb ) { - if ( close_enough_2d(p, lb->GetPosition()) ) { -// return std::distance( tg_node_list.begin(), lb ); - return; - } - } - } - - // we didn't find an existing node - insert new one in correct order - idx = std::lower_bound( lb, ub, node, compare_position ); - tg_node_list.insert( idx, node ); - -// return std::distance( tg_node_list.begin(), idx ); -#endif -} -#endif - -void TGNodes::linear_unique_add( const Point3D& p ) { +void TGNodes::linear_unique_add( const SGGeod& p ) { node_list_iterator current, last; // see if point already exists @@ -134,56 +97,7 @@ void TGNodes::linear_unique_add( const Point3D& p ) { tg_node_list.push_back( node ); } -#if 0 -void TGNodes::sorted_unique_add_fixed_elevation( const Point3D& p ) { -# if 0 - TGNode node( p ); - node.SetFixedPosition(true); - - node_list_iterator idx = std::lower_bound( tg_node_list.begin(), tg_node_list.end(), node, compare_position ); - - if ( idx != tg_node_list.end() ) { - if ( close_enough_2d( p, idx->GetPosition() ) ) { - SG_LOG(SG_GENERAL, SG_ALERT, "AddFixedLocation: node " << p << " exists at " << std::distance(tg_node_list.begin(), idx) ); - idx->SetPosition( p ); - idx->SetFixedPosition( true ); - } else { - tg_node_list.insert( idx, node ); -} - } -#else - TGNode node( p ); - node.SetFixedPosition(true); - - node_list_iterator lb, ub, idx; - - // first, find the range to search - ub=lower_bound( tg_node_list.begin(), tg_node_list.end(), node, fuzzy_compare_xposition ); - lb=upper_bound( tg_node_list.begin(), tg_node_list.end(), node, fuzzy_compare_xposition ); - - // then do a normal linear search in the range - if ( lb != tg_node_list.end() ) { - for ( ; lb != ub; ++lb ) { - if ( close_enough_2d(p, lb->GetPosition()) ) { - lb->SetPosition( p ); - lb->SetFixedPosition( true ); - -// return std::distance( tg_node_list.begin(), lb ); - return; - } - } - } - - // we didn't find an existing node - insert new one in correct order - idx = std::lower_bound( lb, ub, node, compare_position ); - tg_node_list.insert( idx, node ); - -// return std::distance( tg_node_list.begin(), idx ); -#endif -} -#endif - -void TGNodes::linear_unique_add_fixed_elevation( const Point3D& p ) { +void TGNodes::linear_unique_add_fixed_elevation( const SGGeod& p ) { node_list_iterator current, last; // see if point already exists @@ -207,51 +121,60 @@ void TGNodes::linear_unique_add_fixed_elevation( const Point3D& p ) { // add to list tg_node_list.push_back( node ); - } -point_list TGNodes::get_geod_nodes( void ) const { - point_list points; +void TGNodes::get_geod_nodes( std::vector<SGGeod>& points ) const { const_node_list_iterator current, last; // see if point already exists current = tg_node_list.begin(); last = tg_node_list.end(); + points.clear(); for ( ; current != last; ++current ) { points.push_back( (*current).GetPosition() ); } - - return points; } // TODO: if the list is sorted, we should be able to get the range from x=min to x=max, // then within that range, add each point where y is within ymin, max // still linear search, but should be less points -point_list TGNodes::get_geod_inside( Point3D min, Point3D max ) const { - point_list points; +const double fgPoint3_Epsilon = 0.000001; + +static bool IsWithin( const SGGeod pt, double xmin, double xmax, double ymin, double ymax ) +{ + return ( (xmin <= pt.getLongitudeDeg()) && (ymin <= pt.getLatitudeDeg()) && + (xmax >= pt.getLongitudeDeg()) && (ymax >= pt.getLatitudeDeg()) ); +} + +static bool IsAlmostWithin( const SGGeod pt, const SGGeod& min, const SGGeod& max ) +{ + // make sure we take epsilon into account + return ( IsWithin(pt, + min.getLongitudeDeg() - fgPoint3_Epsilon, + max.getLongitudeDeg() + fgPoint3_Epsilon, + min.getLatitudeDeg() - fgPoint3_Epsilon, + max.getLatitudeDeg() + fgPoint3_Epsilon ) ); +} + +void TGNodes::get_geod_inside( const SGGeod& min, const SGGeod& max, std::vector<SGGeod>& points ) const { const_node_list_iterator current, last; // see if point already exists current = tg_node_list.begin(); last = tg_node_list.end(); + points.clear(); for ( ; current != last; ++current ) { - Point3D pt = (*current).GetPosition(); + SGGeod pt = (*current).GetPosition(); - if ( pt.IsAlmostWithin( min, max ) ) { + if ( IsAlmostWithin( pt, min, max ) ) { points.push_back( pt ); - } else { - if ( (pt < max) && (pt > min) ) { - SG_LOG(SG_GENERAL, SG_ALERT, "pt " << pt << " fails IsAlmostWithin, but sholdn't have: min " << min << " max " << max ); - } } } - - return points; } -void TGNodes::get_geod_edge( SGBucket b, point_list& north, point_list& south, point_list& east, point_list& west ) const { +void TGNodes::get_geod_edge( const SGBucket& b, std::vector<SGGeod>& north, std::vector<SGGeod>& south, std::vector<SGGeod>& east, std::vector<SGGeod>& west ) const { const_node_list_iterator current, last; double north_compare = b.get_center_lat() + 0.5 * b.get_height(); double south_compare = b.get_center_lat() - 0.5 * b.get_height(); @@ -262,91 +185,60 @@ void TGNodes::get_geod_edge( SGBucket b, point_list& north, point_list& south, p current = tg_node_list.begin(); last = tg_node_list.end(); + north.clear(); + south.clear(); + east.clear(); + west.clear(); + for ( ; current != last; ++current ) { - Point3D pt = (*current).GetPosition(); + SGGeod pt = (*current).GetPosition(); // may save the same point twice - so we get all the corners - if ( fabs(pt.y() - north_compare) < SG_EPSILON) { + if ( fabs(pt.getLatitudeDeg() - north_compare) < SG_EPSILON) { north.push_back( pt ); } - if ( fabs(pt.y() - south_compare) < SG_EPSILON) { + if ( fabs(pt.getLatitudeDeg() - south_compare) < SG_EPSILON) { south.push_back( pt ); } - if ( fabs(pt.x() - east_compare) < SG_EPSILON) { + if ( fabs(pt.getLongitudeDeg() - east_compare) < SG_EPSILON) { east.push_back( pt ); } - if ( fabs(pt.x() - west_compare) < SG_EPSILON) { + if ( fabs(pt.getLongitudeDeg() - west_compare) < SG_EPSILON) { west.push_back( pt ); } } } -std::vector< SGVec3d > TGNodes::get_wgs84_nodes_as_SGVec3d( void ) const { +void TGNodes::get_wgs84_nodes( std::vector<SGVec3d>& points ) const { const_node_list_iterator current, last; - std::vector< SGVec3d > points; current = tg_node_list.begin(); last = tg_node_list.end(); + points.clear(); for ( ; current != last; ++current ) { - points.push_back( (*current).GetWgs84AsSGVec3d() ); + points.push_back( (*current).GetWgs84() ); } - - return points; } - -point_list TGNodes::get_wgs84_nodes_as_Point3d( void ) const { - const_node_list_iterator current, last; - point_list points; - - current = tg_node_list.begin(); - last = tg_node_list.end(); - - for ( ; current != last; ++current ) { - points.push_back( (*current).GetWgs84AsPoint3D() ); - } - - return points; -} - -node_list TGNodes::get_fixed_elevation_nodes( void ) const { - node_list fixed_elev; +void TGNodes::get_normals( std::vector<SGVec3f>& normals ) const { const_node_list_iterator current, last; // see if point already exists current = tg_node_list.begin(); last = tg_node_list.end(); + normals.clear(); for ( ; current != last; ++current ) { - if ( (*current).GetFixedPosition() ) { - fixed_elev.push_back( (*current) ); - } + normals.push_back( (*current).GetNormal() ); } - - return fixed_elev; -} - -point_list TGNodes::get_normals( void ) const { - point_list points; - const_node_list_iterator current, last; - - // see if point already exists - current = tg_node_list.begin(); - last = tg_node_list.end(); - - for ( ; current != last; ++current ) { - points.push_back( (*current).GetNormal() ); - } - - return points; } void TGNodes::Dump( void ) { for (unsigned int i=0; i<tg_node_list.size(); i++) { TGNode node = tg_node_list[ i ]; std::string fixed; - + if ( node.GetFixedPosition() ) { fixed = " z is fixed elevation "; } else { @@ -363,43 +255,15 @@ void TGNodes::Dump( void ) { } } -// input from stream -std::istream& operator >> ( std::istream& in, TGNode& n ) -{ - int i, nCount; - - // Load a tgnode - in >> n.position; - n.CalcWgs84(); - - in >> n.normal; - in >> n.fixed_position; - in >> n.fixed_normal; - - in >> nCount; - for (i=0; i<nCount; i++) { - TGFaceLookup face; - - in >> face.area; - in >> face.shape; - in >> face.seg; - in >> face.tri; - - n.faces.push_back( face ); - } - - return in; -} - void TGNode::LoadFromGzFile(gzFile& fp) { int i, nCount; // Load a tgnode - sgReadPoint3D( fp, position); + sgReadGeod( fp, position); CalcWgs84(); - sgReadPoint3D( fp, normal ); + sgReadVec3( fp, normal ); sgReadInt( fp, (int*)&fixed_position ); sgReadInt( fp, (int*)&fixed_normal ); @@ -442,8 +306,8 @@ void TGNode::SaveToGzFile(gzFile& fp) int i, nCount; // Save a tgnode - sgWritePoint3D( fp, position ); - sgWritePoint3D( fp, normal ); + sgWriteGeod( fp, position ); + sgWriteVec3( fp, normal ); sgWriteInt( fp, (int)fixed_position ); sgWriteInt( fp, (int)fixed_normal ); @@ -457,26 +321,6 @@ void TGNode::SaveToGzFile(gzFile& fp) } } -// input from stream -std::istream& operator >> ( std::istream& in, TGNodes& ns ) -{ - int i, nCount; - - // Load sorted flag - in >> ns.sorted; - // Load all tgnodes - in >> nCount; - - for (i=0; i<nCount; i++) { - TGNode node; - in >> node; - - ns.tg_node_list.push_back( node ); - } - - return in; -} - void TGNodes::LoadFromGzFile(gzFile& fp) { int i, nCount; @@ -526,22 +370,4 @@ void TGNodes::SaveToGzFile(gzFile& fp) for (i=0; i<nCount; i++) { tg_node_list[i].SaveToGzFile( fp ); } -} - -#if 0 -bool TGNodes::LookupFixedElevation( Point3D p, double* z ) -{ - int index = find( p ); - bool found = false; - - if (index >= 0) { - TGNode node = tg_node_list[index]; - if ( node.GetFixedPosition() ) { - *z = tg_node_list[index].GetPosition().z(); - found = true; - } - } - - return found; -} -#endif +} \ No newline at end of file diff --git a/src/Lib/Geometry/tg_nodes.hxx b/src/Lib/Geometry/tg_nodes.hxx index 18e67534..7867cbaf 100644 --- a/src/Lib/Geometry/tg_nodes.hxx +++ b/src/Lib/Geometry/tg_nodes.hxx @@ -5,18 +5,16 @@ # include <config.h> #endif -#ifndef __cplusplus +#ifndef __cplusplus # error This library requires C++ -#endif +#endif #include <cstdlib> #include <simgear/compiler.h> #include <simgear/bucket/newbucket.hxx> -#include <simgear/math/sg_types.hxx> +//#include <simgear/math/sg_types.hxx> #include <simgear/io/lowlevel.hxx> -#include <Geometry/point3d.hxx> - #define FG_PROXIMITY_EPSILON 0.000001 #define FG_COURSE_EPSILON 0.0001 @@ -36,12 +34,11 @@ public: TGNode() { // constructor for serialization only } - - TGNode( Point3D p ) { + + TGNode( SGGeod p ) { position = p; - normal = Point3D(); CalcWgs84(); - + fixed_position = false; // no matter what - don't move x, y, or z (likely a hole around an airport generated ny genapts) fixed_normal = false; // no matter what - don't modify the normal - likely on a normal generated on a shared edge @@ -57,8 +54,7 @@ public: inline void CalcWgs84() { - SGGeod geod = SGGeod::fromDegM( position.x(), position.y(), position.z() ); - wgs84 = SGVec3d::fromGeod(geod); + wgs84 = SGVec3d::fromGeod(position); } inline void AddFace( unsigned int area, unsigned int shape, unsigned int segment, unsigned int tri ) @@ -76,10 +72,9 @@ public: inline bool GetFixedPosition( void ) const { return fixed_position; } inline void SetFixedNormal( bool fix ) { fixed_normal = fix; } inline bool GetFixedNormal( void ) const { return fixed_normal; } - inline SGVec3d GetWgs84AsSGVec3d( void ) const { return wgs84; } - inline Point3D GetWgs84AsPoint3D( void ) const { return Point3D::fromSGVec3( wgs84 ); } + inline SGVec3d GetWgs84( void ) const { return wgs84; } - inline void SetPosition( const Point3D& p ) + inline void SetPosition( const SGGeod& p ) { if (!fixed_position) { position = p; @@ -90,24 +85,23 @@ public: inline void SetElevation( double z ) { if (!fixed_position) { - position.setelev( z ); + position.setElevationM( z ); CalcWgs84(); } } - inline Point3D GetPosition( void ) const { return position; } - inline void SetNormal( const Point3D& n ) { normal = n; } - inline Point3D GetNormal( void ) const { return normal; } + inline SGGeod GetPosition( void ) const { return position; } + inline void SetNormal( const SGVec3f& n ) { normal = n; } + inline SGVec3f GetNormal( void ) const { return normal; } void SaveToGzFile( gzFile& fp ); void LoadFromGzFile( gzFile& fp ); // Friends for serialization - friend std::istream& operator>> ( std::istream&, TGNode& ); friend std::ostream& operator<< ( std::ostream&, const TGNode& ); private: - Point3D position; - Point3D normal; + SGGeod position; + SGVec3f normal; SGVec3d wgs84; bool fixed_position; @@ -119,9 +113,7 @@ typedef std::vector < TGNode > node_list; typedef node_list::iterator node_list_iterator; typedef node_list::const_iterator const_node_list_iterator; - /* This class handles ALL of the nodes in a tile : 3d nodes in elevation data, 2d nodes generated from landclass, etc) */ - class TGNodes { public: @@ -139,7 +131,7 @@ public: // Add a point to the point list if it doesn't already exist. // Returns the index (starting at zero) of the point in the list. - void unique_add( const Point3D& p ) { + void unique_add( const SGGeod& p ) { if ( !sorted ) { linear_unique_add( p ); } else { @@ -151,7 +143,7 @@ public: // Add a point to the point list if it doesn't already exist // (checking all three dimensions.) Returns the index (starting // at zero) of the point in the list. - void unique_add_fixed_elevation( const Point3D& p ) { + void unique_add_fixed_elevation( const SGGeod& p ) { if ( !sorted ) { linear_unique_add_fixed_elevation( p ); } else { @@ -162,7 +154,7 @@ public: // Find the index of the specified point (compair to the same // tolerance as unique_add(). Returns -1 if not found. - int find( const Point3D& p ) const { + int find( const SGGeod& p ) const { if ( sorted ) { return sorted_find( p ); } else { @@ -170,33 +162,30 @@ public: } } - node_list get_fixed_elevation_nodes( void ) const; - void SortNodes( void ); void SetElevation( int idx, double z ) { tg_node_list[idx].SetElevation( z ); } - Point3D GetNormal( int idx ) const { return tg_node_list[idx].GetNormal(); } - void SetNormal( int idx, Point3D n ) { tg_node_list[idx].SetNormal( n ); } + SGVec3f GetNormal( int idx ) const { return tg_node_list[idx].GetNormal(); } + void SetNormal( int idx, SGVec3f n ) { tg_node_list[idx].SetNormal( n ); } // return the master node list inline node_list& get_node_list() { return tg_node_list; } inline const node_list& get_node_list() const { return tg_node_list; } // return a point list of geodetic nodes - point_list get_geod_nodes() const; + void get_geod_nodes( std::vector<SGGeod>& points ) const; // Find all the nodes within a bounding box - point_list get_geod_inside( Point3D min, Point3D max ) const; + void get_geod_inside( const SGGeod& min, const SGGeod& max, std::vector<SGGeod>& points ) const; // Find a;; the nodes on the tile edges - void get_geod_edge( SGBucket b, point_list& north, point_list& south, point_list& east, point_list& west ) const; + void get_geod_edge( const SGBucket& b, std::vector<SGGeod>& north, std::vector<SGGeod>& south, std::vector<SGGeod>& east, std::vector<SGGeod>& west ) const; // return a point list of wgs84 nodes - std::vector< SGVec3d > get_wgs84_nodes_as_SGVec3d() const; - point_list get_wgs84_nodes_as_Point3d() const; + void get_wgs84_nodes( std::vector<SGVec3d>& points ) const; // return a point list of normals - point_list get_normals() const; + void get_normals( std::vector<SGVec3f>& normals ) const; // return the ith point inline TGNode get_node( int i ) const { return tg_node_list[i]; } @@ -214,67 +203,64 @@ public: void LoadFromGzFile( gzFile& fp ); // Friends for serialization - friend std::istream& operator>> ( std::istream&, TGNodes& ); friend std::ostream& operator<< ( std::ostream&, const TGNodes& ); private: - void linear_unique_add( const Point3D& p ); - void linear_unique_add_fixed_elevation( const Point3D& p ); + void linear_unique_add( const SGGeod& p ); + void linear_unique_add_fixed_elevation( const SGGeod& p ); - int sorted_find( const Point3D& p ) const; - int linear_find( const Point3D& p ) const; + int sorted_find( const SGGeod& p ) const; + int linear_find( const SGGeod& p ) const; node_list tg_node_list; bool sorted; // return true of the two points are "close enough" as defined by // FG_PROXIMITY_EPSILON - bool close_enough_2d( const Point3D& p1, const Point3D& p2 ) const; + bool close_enough_2d( const SGGeod& p1, const SGGeod& p2 ) const; // return true of the two points are "close enough" as defined by // FG_PROXIMITY_EPSILON - bool close_enough_3d( const Point3D& p1, const Point3D& p2 ) const; + bool close_enough_3d( const SGGeod& p1, const SGGeod& p2 ) const; // return true of the two points are "close enough" as defined by // FG_COURSE_EPSILON - bool course_close_enough( const Point3D& p1, const Point3D& p2 ); + bool course_close_enough( const SGGeod& p1, const SGGeod& p2 ); }; // return true of the two points are "close enough" as defined by // FG_PROXIMITY_EPSILON checking just x and y dimensions -inline bool TGNodes::close_enough_2d( const Point3D& p1, const Point3D& p2 ) +inline bool TGNodes::close_enough_2d( const SGGeod& p1, const SGGeod& p2 ) const { - if ( ( fabs(p1.x() - p2.x()) < FG_PROXIMITY_EPSILON ) && - ( fabs(p1.y() - p2.y()) < FG_PROXIMITY_EPSILON ) ) { + if ( ( fabs(p1.getLongitudeDeg() - p2.getLongitudeDeg()) < FG_PROXIMITY_EPSILON ) && + ( fabs(p1.getLatitudeDeg() - p2.getLatitudeDeg()) < FG_PROXIMITY_EPSILON ) ) { return true; } else { return false; } } - // return true of the two points are "close enough" as defined by // FG_PROXIMITY_EPSILON check all three dimensions -inline bool TGNodes::close_enough_3d( const Point3D& p1, const Point3D& p2 ) +inline bool TGNodes::close_enough_3d( const SGGeod& p1, const SGGeod& p2 ) const { - if ( ( fabs(p1.x() - p2.x()) < FG_PROXIMITY_EPSILON ) && - ( fabs(p1.y() - p2.y()) < FG_PROXIMITY_EPSILON ) && - ( fabs(p1.z() - p2.z()) < FG_PROXIMITY_EPSILON ) ) { + if ( ( fabs(p1.getLongitudeDeg() - p2.getLongitudeDeg()) < FG_PROXIMITY_EPSILON ) && + ( fabs(p1.getLatitudeDeg() - p2.getLatitudeDeg()) < FG_PROXIMITY_EPSILON ) && + ( fabs(p1.getElevationM() - p2.getElevationM()) < FG_PROXIMITY_EPSILON ) ) { return true; } else { return false; } } - // return true of the two points are "close enough" as defined by // FG_COURSE_EPSILON -inline bool TGNodes::course_close_enough( const Point3D& p1, const Point3D& p2 ) +inline bool TGNodes::course_close_enough( const SGGeod& p1, const SGGeod& p2 ) { - if ( ( fabs(p1.x() - p2.x()) < FG_COURSE_EPSILON ) && - ( fabs(p1.y() - p2.y()) < FG_COURSE_EPSILON ) ) { + if ( ( fabs(p1.getLongitudeDeg() - p2.getLongitudeDeg()) < FG_COURSE_EPSILON ) && + ( fabs(p1.getLatitudeDeg() - p2.getLatitudeDeg()) < FG_COURSE_EPSILON ) ) { return true; } else { return false;