genapts fixes

- remove splinter remaval / merging
- remove snap after adding colinear nodes
- added zfighting detector - intersects all generated triangles against one another
- this fixes ELLX, and to a large extent all 'large' airport self z-fighting.
- very small z-fighting still exists at 72 airports.  This is down from over 100 before the fix.

src/Airports/GenAirports850/airport.cxx

@@ -14,6 +14,7 @@
 
 #include <terragear/tg_polygon.hxx>
 #include <terragear/tg_chopper.hxx>
+#include <terragear/tg_shapefile.hxx>
 #include <terragear/tg_unique_geod.hxx>
 #include <terragear/tg_unique_vec3f.hxx>
 #include <terragear/tg_unique_vec2f.hxx>
@@ -233,6 +234,8 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
     tglightcontour_list rwy_lights;
 
     bool make_shapefiles = false;
+    char debug_root[32];
+    sprintf(debug_root, "./airport_dbg/%s/", icao.c_str() );
 
     // parse main airport information
     double apt_lon = 0.0, apt_lat = 0.0;
@@ -311,7 +314,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
         for ( unsigned int i=0; i<runways.size(); i++ )
         {
             TG_LOG(SG_GENERAL, SG_DEBUG, "Build Runway " << i + 1 << " of " << runways.size());
-            slivers.clear();
+            //slivers.clear();
 
             if ( isDebugRunway(i) ) {
                 sprintf( shapefile_name, "runway_%d", i );
@@ -330,7 +333,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
             }
 
             // Now try to merge any slivers we found
-            tgPolygon::MergeSlivers( rwy_polys, slivers );
+            // tgPolygon::MergeSlivers( rwy_polys, slivers );
         }
 
         log_time = time(0);
@@ -354,7 +357,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
         for ( unsigned int i=0; i<helipads.size(); i++ )
         {
             TG_LOG(SG_GENERAL, SG_DEBUG, "Build helipad " << i + 1 << " of " << helipads.size());
-            slivers.clear();
+            //slivers.clear();
 
             if (boundary.size())
             {
@@ -366,7 +369,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
             }
 
             // Now try to merge any slivers we found
-            tgPolygon::MergeSlivers( rwy_polys, slivers );
+            //tgPolygon::MergeSlivers( rwy_polys, slivers );
         }
     }
 
@@ -377,7 +380,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
         for ( unsigned int i=0; i<pavements.size(); i++ )
         {
             TG_LOG(SG_GENERAL, SG_DEBUG, "Build Pavement " << i + 1 << " of " << pavements.size() << " : " << pavements[i]->GetDescription());
-            slivers.clear();
+            //slivers.clear();
 
             if ( isDebugPavement(i) ) {
                 sprintf( shapefile_name, "pvmnt_%d", i );
@@ -396,8 +399,8 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
             }
 
             // Now try to merge any slivers we found
-            tgPolygon::MergeSlivers( rwy_polys, slivers );
+            //tgPolygon::MergeSlivers( rwy_polys, slivers );
-            tgPolygon::MergeSlivers( pvmt_polys, slivers );
+            //tgPolygon::MergeSlivers( pvmt_polys, slivers );
         }
 
         log_time = time(0);
@@ -411,7 +414,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
         for ( unsigned int i=0; i<taxiways.size(); i++ )
         {
             TG_LOG(SG_GENERAL, SG_DEBUG, "Build Taxiway " << i + 1 << " of " << taxiways.size());
-            slivers.clear();
+            //slivers.clear();
 
             if ( isDebugTaxiway(i) ) {
                 sprintf( shapefile_name, "taxiway_%d", i );
@@ -430,8 +433,8 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
             }
 
             // Now try to merge any slivers we found
-            tgPolygon::MergeSlivers( rwy_polys, slivers );
+            //tgPolygon::MergeSlivers( rwy_polys, slivers );
-            tgPolygon::MergeSlivers( pvmt_polys, slivers );
+            //tgPolygon::MergeSlivers( pvmt_polys, slivers );
         }
     }
 
@@ -445,12 +448,12 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
             if ( runways[i]->GetsShoulder() )
             {
-                slivers.clear();
+                //slivers.clear();
                 runways[i]->BuildShoulder( rwy_polys, slivers, pvmt_accum );
 
                 // Now try to merge any slivers we found
-                tgPolygon::MergeSlivers( rwy_polys, slivers );
+                //tgPolygon::MergeSlivers( rwy_polys, slivers );
-                tgPolygon::MergeSlivers( pvmt_polys, slivers );
+                //tgPolygon::MergeSlivers( pvmt_polys, slivers );
             }
         }
     }
@@ -465,12 +468,12 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
             if ( helipads[i]->GetsShoulder() )
             {
-                slivers.clear();
+                //slivers.clear();
                 helipads[i]->BuildShoulder( rwy_polys, slivers, pvmt_accum );
 
                 // Now try to merge any slivers we found
-                tgPolygon::MergeSlivers( rwy_polys, slivers );
+                //tgPolygon::MergeSlivers( rwy_polys, slivers );
-                tgPolygon::MergeSlivers( pvmt_polys, slivers );
+                //tgPolygon::MergeSlivers( pvmt_polys, slivers );
             }
         }
     }
@@ -628,14 +631,15 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
     TG_LOG( SG_GENERAL, SG_ALERT, "Finished cleaning polys for " << icao << " at " << DebugTimeToString(log_time) );
 
     base_poly = tgPolygon::AddColinearNodes( base_poly, tmp_pvmt_nodes );
-    base_poly = tgPolygon::Snap( base_poly, gSnap );
+    // base_poly = tgPolygon::Snap( base_poly, gSnap );
-
+    
     // Finally find slivers in base
-    slivers.clear();
+    //slivers.clear();
-    tgPolygon::RemoveSlivers( base_poly, slivers );
+    //tgPolygon::RemoveSlivers( base_poly, slivers );
-    tgPolygon::MergeSlivers( rwy_polys, slivers );
+    //tgPolygon::MergeSlivers( rwy_polys, slivers );
-    tgPolygon::MergeSlivers( pvmt_polys, slivers );
+    //tgPolygon::MergeSlivers( pvmt_polys, slivers );
-
+    
+#if 0
     // Then snap rwy and pavement to grid (was done right after adding intermediate nodes...)
     for ( unsigned int k = 0; k < rwy_polys.size(); ++k )
     {
@@ -653,10 +657,15 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
         poly = tgPolygon::RemoveBadContours( poly );
     	pvmt_polys[k] = poly;
     }
-
+#endif
+    
     cleanup_end.stamp();
     cleanup_time = cleanup_end - cleanup_start;
-
+    
+    /* before tessellating the base, make sure there are no
+       intersecting contours */
+    base_poly = tgPolygon::Simplify( base_poly );
+    
     triangulation_start.stamp();
 
     // tesselate the polygons and prepair them for final output
@@ -669,6 +678,14 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
             TG_LOG(SG_GENERAL, SG_DEBUG, "contours before " << rwy_polys[i].Contours() << " total points before = " << rwy_polys[i].TotalNodes());
             rwy_polys[i].Tesselate();
+            
+#if 0
+            // dump the triangles for each poly
+            char desc[128];
+            sprintf( desc, "poly_%06d", i );
+            tgShapefile::FromTriangles( rwy_polys[i], debug_root, "rwy_polys", desc );
+#endif
+
             TG_LOG(SG_GENERAL, SG_DEBUG, "triangles after = " << rwy_polys[i].Triangles());
             rwy_polys[i].Texture();
         }
@@ -684,6 +701,13 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
             TG_LOG(SG_GENERAL, SG_DEBUG, "contours before " << pvmt_polys[i].Contours() << " total points before = " << pvmt_polys[i].TotalNodes());
             pvmt_polys[i].Tesselate();
+
+#if 0            
+            char desc[128];
+            sprintf( desc, "poly_%06d", i );
+            tgShapefile::FromTriangles( pvmt_polys[i], debug_root, "pvmt_polys", desc );
+#endif
+
             TG_LOG(SG_GENERAL, SG_DEBUG, "triangles after = " << pvmt_polys[i].Triangles());
             pvmt_polys[i].Texture();
         }
@@ -699,17 +723,27 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
             TG_LOG(SG_GENERAL, SG_DEBUG, "contours before " << line_polys[i].Contours() << " total points before = " << line_polys[i].TotalNodes());
             line_polys[i].Tesselate();
+            
+#if 0
+            char desc[128];
+            sprintf( desc, "poly_%06d", i );
+            tgShapefile::FromTriangles( line_polys[i], debug_root, "line_polys", desc );
+#endif
+
             TG_LOG(SG_GENERAL, SG_DEBUG, "triangles after = " << line_polys[i].Triangles());
             line_polys[i].Texture();
         }
     }
-
+    
-    /* before tessellating the base, make sure there are no
-       intersecting contours */
-    base_poly = tgPolygon::Simplify( base_poly );
-
     TG_LOG(SG_GENERAL, SG_INFO, "Tesselating base poly : " << base_poly.Contours() << " contours " );
     base_poly.Tesselate();
+
+#if 0    
+    char desc[128];
+    sprintf( desc, "poly" );
+    tgShapefile::FromTriangles( base_poly, debug_root, "base_poly", desc );
+#endif
+
     TG_LOG(SG_GENERAL, SG_INFO, "Tesselating base poly - done : Triangles = " << base_poly.Triangles());
     // should we texture base here?
     base_poly.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() );
@@ -726,6 +760,15 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
 
     triangulation_time = triangulation_end - triangulation_start;
 
+#if 0    
+    // CheckZFighting( "post triangulation", debug_root, base_poly, rwy_polys, pvmt_polys );
+    if ( CheckZFightingTriangles( "post triangulation tris", debug_root, base_poly, rwy_polys, pvmt_polys ) ) {
+        char cmd[128];
+        sprintf( cmd, "echo %s >> zfight.txt\n", icao.c_str() );
+        system ( cmd );
+    }
+#endif
+
     //
     // We should now have the runway polygons all generated with their
     // corresponding triangles and texture coordinates, and the
@@ -1166,3 +1209,142 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
     chopper.Add( apt_clearing, "Airport" );
     chopper.Save( false );
 }
+
+bool Airport::CheckZFightingTriangles( const char* prefix, const char* debug_root, const tgPolygon& base_poly, const tgpolygon_list& rwy_polys, const tgpolygon_list& pvmt_polys )
+{
+    char layer[128];
+    char desc[128];
+
+    sprintf( layer, "zfighting_%s", prefix );
+    double min_area_thresh = 1.0e-10;
+    
+    bool zfighting = false;
+    
+    if ( rwy_polys.size() ) {
+        for ( unsigned int i = 0; i < rwy_polys.size(); ++i ) {
+            tgPolygon   subject = rwy_polys[i];
+            for ( unsigned int j = 0; j < subject.Triangles(); ++j ) {            
+                tgTriangle  subTri = subject.GetTriangle( j );
+                tgRectangle subBB  = subTri.GetBoundingBox();
+    
+                for ( unsigned int k = 0; k < pvmt_polys.size(); ++k ) {
+                    tgPolygon   test   = pvmt_polys[k];
+                    for ( unsigned int l = 0; l < test.Triangles(); ++l ) {
+                        tgTriangle  testTri = test.GetTriangle( l );
+                        tgRectangle testBB  = testTri.GetBoundingBox();
+                
+                        if ( subBB.intersects( testBB ) ) {
+                            // find the intersection
+                            tgPolygon intersection = tgTriangle::Intersect( subTri, testTri );
+                    
+                            for ( unsigned int m = 0; m < intersection.Contours(); m++ ) {
+                                tgContour intContour = intersection.GetContour( m );
+                                if ( (intContour.GetSize() > 2) && (intContour.GetArea() > min_area_thresh) ) {                        
+                                    TG_LOG( SG_GENERAL, SG_ALERT, prefix << "Z-FIGHTING between runway poly " << i << " and pavement poly " << k << " contour has " << intContour.GetSize() << " nodes " << " area is " << intContour.GetArea() );
+                                    sprintf( desc, "rwy_%06d_pvmt_%06d", i, j );
+                                    tgShapefile::FromContour( intContour, debug_root, layer, desc );
+                                    zfighting = true;
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // intersect each pavement poly with each pavement, runway, and base poly
+    if ( pvmt_polys.size() ) {
+        for ( unsigned int i = 0; i < pvmt_polys.size(); ++i ) {
+            tgPolygon   subject = pvmt_polys[i];
+            for ( unsigned int j = 0; j < subject.Triangles(); ++j ) {            
+                tgTriangle  subTri = subject.GetTriangle( j );
+                tgRectangle subBB  = subTri.GetBoundingBox();
+    
+                for ( unsigned int k = 0; k < rwy_polys.size(); ++k ) {
+                    tgPolygon   test   = rwy_polys[k];
+                    for ( unsigned int l = 0; l < test.Triangles(); ++l ) {
+                        tgTriangle  testTri = test.GetTriangle( l );
+                        tgRectangle testBB  = testTri.GetBoundingBox();
+                
+                        if ( subBB.intersects( testBB ) ) {
+                            // find the intersection
+                            tgPolygon intersection = tgTriangle::Intersect( subTri, testTri );
+
+                            for ( unsigned int m = 0; m < intersection.Contours(); m++ ) {
+                                tgContour intContour = intersection.GetContour( m );
+                                if ( (intContour.GetSize() > 2) && (intContour.GetArea() > min_area_thresh) ) {                        
+                                    TG_LOG( SG_GENERAL, SG_ALERT, prefix << "Z-FIGHTING between pavement poly " << i << " and runway poly " << k << " contour has " << intContour.GetSize() << " nodes " << " area is " << intContour.GetArea() );
+                                    sprintf( desc, "pvmt_%06d_rwy_%06d", i, j );
+                                    tgShapefile::FromContour( intContour, debug_root, layer, desc );
+                                    zfighting = true;
+                                }
+                            }                            
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    // intersect base poly with each runway and pavement poly
+    tgPolygon   subject = base_poly;
+
+    for ( unsigned int i = 0; i < subject.Triangles(); ++i ) {
+        tgTriangle  subTri = subject.GetTriangle( i );
+        tgRectangle subBB  = subTri.GetBoundingBox();
+
+        for ( unsigned int j = 0; j < rwy_polys.size(); ++j ) {
+            tgPolygon   test   = rwy_polys[j];
+            for ( unsigned int k = 0; k < test.Triangles(); ++k ) {
+                tgTriangle  testTri = test.GetTriangle( k );
+                tgRectangle testBB  = testTri.GetBoundingBox();
+        
+                if ( subBB.intersects( testBB ) ) {
+                    // find the intersection
+                    tgPolygon intersection = tgTriangle::Intersect( subTri, testTri );
+            
+                    for ( unsigned int m = 0; m < intersection.Contours(); m++ ) {
+                        tgContour intContour = intersection.GetContour( m );
+                        if ( (intContour.GetSize() > 2) && (intContour.GetArea() > min_area_thresh) ) {
+                            TG_LOG( SG_GENERAL, SG_ALERT, prefix << "Z-FIGHTING between base poly and runway poly " << j << " contour has " << intContour.GetSize() << " nodes " << " area is " << intContour.GetArea() );
+                            sprintf( desc, "base_rwy_%06d", j );
+                            tgShapefile::FromContour( intContour, debug_root, layer, desc );
+                            zfighting = true;
+                        }
+                    }                    
+                }
+            }
+        }
+    }
+
+    for ( unsigned int i = 0; i < subject.Triangles(); ++i ) {
+        tgTriangle  subTri = subject.GetTriangle( i );
+        tgRectangle subBB  = subTri.GetBoundingBox();
+
+        for ( unsigned int j = 0; j < pvmt_polys.size(); ++j ) {
+            tgPolygon   test   = pvmt_polys[j];
+            for ( unsigned int k = 0; k < test.Triangles(); ++k ) {
+                tgTriangle  testTri = test.GetTriangle( k );
+                tgRectangle testBB  = testTri.GetBoundingBox();
+        
+                if ( subBB.intersects( testBB ) ) {
+                    // find the intersection
+                    tgPolygon intersection = tgTriangle::Intersect( subTri, testTri );
+            
+                    for ( unsigned int m = 0; m < intersection.Contours(); m++ ) {
+                        tgContour intContour = intersection.GetContour( m );
+                        if ( (intContour.GetSize() > 2) && (intContour.GetArea() > min_area_thresh) ) {                        
+                            TG_LOG( SG_GENERAL, SG_ALERT, prefix << "Z-FIGHTING between base poly and pavement poly " << j << " contour has " << intContour.GetSize() << " nodes " << " area is " << intContour.GetArea() );
+                            sprintf( desc, "base_pvmt_%06d", j );
+                            tgShapefile::FromContour( intContour, debug_root, layer, desc );
+                            zfighting = true;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    return zfighting;
+}

src/Airports/GenAirports850/airport.hxx

@@ -131,6 +131,8 @@ public:
     bool isDebugFeature ( int i );
 
 private:
+    bool CheckZFightingTriangles( const char* prefix, const char* debug_root, const tgPolygon& base_poly, const tgpolygon_list& rwy_polys, const tgpolygon_list& pvmt_polys );
+    
     int          code;               // airport, heliport or sea port
     int          altitude;           // in meters
     std::string  icao;               // airport code

src/Airports/GenAirports850/closedpoly.cxx

@@ -1,3 +1,4 @@
+#include <stdio.h>
 #include <stdlib.h>
 
 #include <simgear/debug/logstream.hxx>
@@ -457,21 +458,36 @@ int ClosedPoly::BuildBtg( tgpolygon_list& rwy_polys, tgcontour_list& slivers, tg
 
 int ClosedPoly::BuildBtg( tgpolygon_list& rwy_polys, tgcontour_list& slivers, tgAccumulator& accum, std::string& shapefile_name )
 {
+    char layer[128];
+    
     if ( is_pavement && pre_tess.Contours() )
     {
         if(  shapefile_name.size() ) {
-            tgShapefile::FromPolygon( pre_tess, "./airport_dbg", std::string("preclip"), shapefile_name );
+            sprintf( layer, "%s_preclip", shapefile_name.c_str() );
-            accum.ToShapefiles( "./airport_dbg", "accum", true );
+            tgShapefile::FromPolygon( pre_tess, "./airport_dbg", layer, std::string("preclip") );
+            
+            pre_tess.Tesselate();
+            sprintf( layer, "%s_preclip_tris", shapefile_name.c_str() );
+            tgShapefile::FromTriangles( pre_tess, "./airport_dbg", layer, std::string("preclip") );
+
+            //accum.ToShapefiles( "./airport_dbg", "accum", true );
         }
 
         tgPolygon clipped = accum.Diff( pre_tess );
         if ( clipped.Contours() ) {
             if(  shapefile_name.size() ) {
-                tgShapefile::FromPolygon( clipped, "./airport_dbg", std::string("postclip"), shapefile_name );
+                sprintf( layer, "%s_postclip", shapefile_name.c_str() );
+                tgShapefile::FromPolygon( clipped, "./airport_dbg", layer, std::string("postclip") );
             }
 
-            tgPolygon::RemoveSlivers( clipped, slivers );
+            // tgPolygon::RemoveSlivers( clipped, slivers );
 
+            if(  shapefile_name.size() ) {
+                clipped.Tesselate();
+                sprintf( layer, "%s_postclip_tris", shapefile_name.c_str() );
+                tgShapefile::FromTriangles( clipped, "./airport_dbg", layer, std::string("postclip") );                
+            }
+            
             clipped.SetMaterial( GetMaterial( surface_type ) );
             clipped.SetTexParams( clipped.GetNode(0,0), 5.0, 5.0, texture_heading );
             clipped.SetTexLimits( 0.0, 0.0, 1.0, 1.0 );

src/Airports/GenAirports850/helipad.cxx

@@ -146,7 +146,7 @@ void Helipad::BuildBtg( tgpolygon_list& rwy_polys,
 
     // Clip the new polygon against what ever has already been created.
     tgPolygon clipped = accum.Diff( helipad );
-    tgPolygon::RemoveSlivers( clipped, slivers );
+    // tgPolygon::RemoveSlivers( clipped, slivers );
 
     // Split long edges to create an object that can better flow with
     // the surface terrain
@@ -171,7 +171,7 @@ void Helipad::BuildBtg( tgpolygon_list& rwy_polys,
 
     tgPolygon outer_poly = tgContour::Diff( outer_area, heli_poly );
     clipped = accum.Diff( outer_poly );
-    tgPolygon::RemoveSlivers( clipped, slivers );
+    // tgPolygon::RemoveSlivers( clipped, slivers );
 
     // Split long edges to create an object that can better flow with
     // the surface terrain
@@ -249,7 +249,7 @@ void Helipad::BuildShoulder( tgpolygon_list& rwy_polys,
 
         // Clip the new polygon against what ever has already been created.
         tgPolygon clipped = accum.Diff( shoulder );
-        tgPolygon::RemoveSlivers( clipped, slivers );
+        // tgPolygon::RemoveSlivers( clipped, slivers );
 
         // Split long edges to create an object that can better flow with
         // the surface terrain

src/Airports/GenAirports850/rwy_gen.cxx

@@ -313,7 +313,7 @@ void Runway::gen_runway_section( const tgPolygon& runway,
         tgShapefile::FromPolygon( clipped, "./airport_dbg", std::string("postclip"), shapefile_name );
     }
 
-    tgPolygon::RemoveSlivers( clipped, slivers );
+    //tgPolygon::RemoveSlivers( clipped, slivers );
 
     // Split long edges to create an object that can better flow with
     // the surface terrain
@@ -451,7 +451,7 @@ void Runway::gen_runway_section( const tgPolygon& runway,
     // Clip the new polygon against what ever has already been created.
     tgPolygon clipped = accum.Diff( section );
 
-    tgPolygon::RemoveSlivers( clipped, slivers );
+    //tgPolygon::RemoveSlivers( clipped, slivers );
 
     // Split long edges to create an object that can better flow with
     // the surface terrain
@@ -875,7 +875,7 @@ void Runway::BuildShoulder( tgpolygon_list& rwy_polys,
 
         // Clip the new polygon against what ever has already been created.
         tgPolygon clipped = accum.Diff( shoulder );
-        tgPolygon::RemoveSlivers( clipped, slivers );
+        //tgPolygon::RemoveSlivers( clipped, slivers );
 
         // Split long edges to create an object that can better flow with
         // the surface terrain

src/BuildTiles/Main/tgconstruct.cxx

@@ -93,7 +93,7 @@ void TGConstruct::run()
         if ( debug_shapes.size() || debug_all ) {
             sprintf(ds_name, "%s/constructdbg_%s", debug_path.c_str(), bucket.gen_index_str().c_str() );
         } else {
-            strcpy( ds_name, "" );
+            sprintf(ds_name, "%s/constructdbg_%s", debug_path.c_str(), bucket.gen_index_str().c_str() );
         }
 
         if ( stage == 1 ) {

src/BuildTiles/Main/tgconstruct_clip.cxx

@@ -218,6 +218,8 @@ bool TGConstruct::ClipLandclassPolys( void ) {
     // finally, what ever is left over goes to ocean
     remains = accum.Diff( safety_base );
     
+    // tgShapefile::FromPolygon( remains, ds_name, "remains", "poly" );
+    
     if ( debug_all || debug_shapes.size() || debug_areas.size() ) {
         char layer[32];
 	    char name[32];

src/Lib/terragear/tg_misc.cxx

@@ -8,7 +8,7 @@
 #include "tg_polygon.hxx"
 #include "tg_misc.hxx"
 
-const double isEqual2D_Epsilon = 0.000001;
+const double isEqual2D_Epsilon = 0.000000001;
 
 #define CLIPPER_FIXEDPT           (1000000000000)
 #define CLIPPER_METERS_PER_DEGREE (111000)

src/Lib/terragear/tg_polygon.hxx

@@ -112,6 +112,29 @@ public:
         idx_list.resize(  3, -1 );
     }
 
+    tgRectangle GetBoundingBox( void ) const
+    {
+        SGGeod min, max;
+
+        double minx =  std::numeric_limits<double>::infinity();
+        double miny =  std::numeric_limits<double>::infinity();
+        double maxx = -std::numeric_limits<double>::infinity();
+        double maxy = -std::numeric_limits<double>::infinity();
+
+        for (unsigned int i = 0; i < node_list.size(); i++) {
+            SGGeod pt = GetNode(i);
+            if ( pt.getLongitudeDeg() < minx ) { minx = pt.getLongitudeDeg(); }
+            if ( pt.getLongitudeDeg() > maxx ) { maxx = pt.getLongitudeDeg(); }
+            if ( pt.getLatitudeDeg()  < miny ) { miny = pt.getLatitudeDeg(); }
+            if ( pt.getLatitudeDeg()  > maxy ) { maxy = pt.getLatitudeDeg(); }
+        }
+
+        min = SGGeod::fromDeg( minx, miny );
+        max = SGGeod::fromDeg( maxx, maxy );
+
+        return tgRectangle( min, max );        
+    }
+
     SGGeod const& GetNode( unsigned int i ) const {
         return node_list[i];
     }
@@ -155,6 +178,9 @@ public:
                             p3.getLongitudeDeg() * p1.getLatitudeDeg() - p1.getLongitudeDeg() * p3.getLatitudeDeg() ));
     }
 
+    static ClipperLib::Path ToClipper( const tgTriangle& subject );
+    static tgPolygon Intersect( const tgTriangle& subject, const tgTriangle& clip );
+    
     void SaveToGzFile( gzFile& fp ) const;
     void LoadFromGzFile( gzFile& fp );
 
@@ -282,6 +308,9 @@ public:
     void AddTriangle( const SGGeod& p1, const SGGeod p2, const SGGeod p3 ) {
         triangles.push_back( tgTriangle( p1, p2, p3 ) );
     }
+    tgTriangle GetTriangle( unsigned int t ) {
+        return triangles[t];
+    }
 
     SGGeod GetTriNode( unsigned int c, unsigned int i ) const {
         return triangles[c].GetNode( i );

src/Lib/terragear/tg_polygon_clip.cxx

@@ -4,6 +4,7 @@
 #include <simgear/debug/logstream.hxx>
 
 #include "tg_polygon.hxx"
+#include "tg_misc.hxx"
 
 static bool clipper_dump = false;
 void tgPolygon::SetClipperDump( bool dmp )
@@ -197,3 +198,52 @@ tgPolygon tgPolygon::FromClipper( const ClipperLib::Paths& subject )
 
     return result;
 }
+
+ClipperLib::Path tgTriangle::ToClipper( const tgTriangle& subject )
+{
+    ClipperLib::Path  contour;
+
+    for ( unsigned int i=0; i<3; i++)
+    {
+        SGGeod p = subject.GetNode( i );
+        contour.push_back( SGGeod_ToClipper(p) );
+    }
+
+    // boundaries need to be orientation: true
+    if ( !Orientation( contour ) ) {
+        //SG_LOG(SG_GENERAL, SG_INFO, "Building clipper contour - boundary contour needs to be reversed" );
+        ReversePath( contour );
+    }
+
+    return contour;
+}
+
+tgPolygon tgTriangle::Intersect( const tgTriangle& subject, const tgTriangle& clip )
+{
+    tgPolygon result;
+    UniqueSGGeodSet all_nodes;
+
+    /* before diff - gather all nodes */
+    for ( unsigned int i = 0; i < 3; ++i ) {
+        all_nodes.add( subject.GetNode(i) );
+    }
+
+    for ( unsigned int i = 0; i < 3; ++i ) {
+        all_nodes.add( clip.GetNode(i) );
+    }
+
+    ClipperLib::Path  clipper_subject = tgTriangle::ToClipper( subject );
+    ClipperLib::Path  clipper_clip    = tgTriangle::ToClipper( clip );
+    ClipperLib::Paths clipper_result;
+
+    ClipperLib::Clipper c;
+    c.Clear();
+    c.AddPath(clipper_subject, ClipperLib::ptSubject, true);
+    c.AddPath(clipper_clip, ClipperLib::ptClip, true);
+    c.Execute(ClipperLib::ctIntersection, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
+
+    result = tgPolygon::FromClipper( clipper_result );
+    result = tgPolygon::AddColinearNodes( result, all_nodes );
+
+    return result;
+}

src/Lib/terragear/tg_polygon_tesselate.cxx

@@ -105,6 +105,9 @@ void tgPolygon::Tesselate( const std::vector<SGGeod>& extra )
 
     SG_LOG( SG_GENERAL, SG_DEBUG, "Tess with extra" );
 
+    // clear any triangles from previous tesselation
+    triangles.clear();
+    
     // Bail right away if polygon is empty
     if ( contours.size() != 0 ) {
         // First, convert the extra points to cgal Points
@@ -162,6 +165,9 @@ void tgPolygon::Tesselate()
 
     SG_LOG( SG_GENERAL, SG_DEBUG, "Tess" );
 
+    // clear any triangles from previous tesselation
+    triangles.clear();
+
     // Bail right away if polygon is empty
     if ( contours.size() != 0 ) {
         // insert each polygon as a constraint into the triangulation

src/Lib/terragear/tg_shapefile.cxx

@@ -65,6 +65,34 @@ void* tgShapefile::OpenLayer( void* ds_id, const char* layer_name ) {
     return (void*)layer;
 }
 
+void* tgShapefile::OpenLineLayer( void* ds_id, const char* layer_name ) {
+    OGRDataSource* datasource = ( OGRDataSource * )ds_id;
+    OGRLayer* layer;
+
+    OGRSpatialReference srs;
+    srs.SetWellKnownGeogCS("WGS84");
+
+    layer = datasource->GetLayerByName( layer_name );
+
+    if ( !layer ) {
+        layer = datasource->CreateLayer( layer_name, &srs, wkbLineString25D, NULL );
+
+        OGRFieldDefn descriptionField( "ID", OFTString );
+        descriptionField.SetWidth( 128 );
+
+        if( layer->CreateField( &descriptionField ) != OGRERR_NONE ) {
+            SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'Description' failed" );
+        }
+    }
+
+    if ( !layer ) {
+        SG_LOG( SG_GENERAL, SG_ALERT, "Creation of layer '" << layer_name << "' failed" );
+        return NULL;
+    }
+
+    return (void*)layer;
+}
+
 void* tgShapefile::CloseDatasource( void* ds_id )
 {
     OGRDataSource* datasource = ( OGRDataSource * )ds_id;
@@ -172,6 +200,46 @@ void tgShapefile::FromContour( const tgContour& subject, const std::string& data
     ds_id = tgShapefile::CloseDatasource( ds_id );
 }
 
+void tgShapefile::FromTriangles( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description )
+{
+    void*          ds_id = tgShapefile::OpenDatasource( datasource.c_str() );
+    SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenDatasource returned " << (unsigned long)ds_id);
+
+    OGRLayer*      l_id  = (OGRLayer *)tgShapefile::OpenLineLayer( ds_id, layer.c_str() );
+    SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id);
+
+    SG_LOG(SG_GENERAL, SG_DEBUG, "subject has " << subject.Contours() << " contours ");
+
+    for ( unsigned int i = 0; i < subject.Triangles(); i++ ) {
+        // FIXME: Current we ignore the hole-flag and instead assume
+        //        that the first ring is not a hole and the rest
+        //        are holes
+        OGRLinearRing  ring;
+        for (unsigned int pt = 0; pt < 3; pt++) {
+            OGRPoint point;
+
+            point.setX( subject.GetTriNode(i, pt).getLongitudeDeg() );
+            point.setY( subject.GetTriNode(i, pt).getLatitudeDeg() );
+            point.setZ( 0.0 );
+            ring.addPoint(&point);
+        }
+        ring.closeRings();
+
+        OGRFeature* feature = NULL;
+        feature = new OGRFeature( l_id->GetLayerDefn() );
+        feature->SetField("ID", description.c_str());
+        feature->SetGeometry(&ring);
+        if( l_id->CreateFeature( feature ) != OGRERR_NONE )
+        {
+            SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile");
+        }
+        OGRFeature::DestroyFeature(feature);
+    }
+
+    // close after each write
+    ds_id = tgShapefile::CloseDatasource( ds_id );    
+}
+
 void tgShapefile::FromPolygon( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description )
 {
     void*          ds_id = tgShapefile::OpenDatasource( datasource.c_str() );

src/Lib/terragear/tg_shapefile.hxx

@@ -8,6 +8,7 @@ class tgShapefile
 public:
     static void  FromContour( const tgContour& subject, const std::string& datasource, const std::string& layer, const std::string& description );
     static void  FromPolygon( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description );
+    static void  FromTriangles( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description );
     static tgPolygon ToPolygon( const void* subject );
 
     static void  FromClipper( const ClipperLib::Paths& subject, const std::string& datasource, const std::string& layer, const std::string& description );
@@ -17,5 +18,6 @@ private:
 
     static void* OpenDatasource( const char* datasource_name );
     static void* OpenLayer( void* ds_id, const char* layer_name );
+    static void* OpenLineLayer( void* ds_id, const char* layer_name );
     static void* CloseDatasource( void* ds_id );
 };