diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx
index 80997615..9bab7044 100644
--- a/src/Airports/GenAirports850/airport.cxx
+++ b/src/Airports/GenAirports850/airport.cxx
@@ -8,6 +8,8 @@
 #include <string>
 #include <ctime>
 
+#include <stdio.h>
+
 #include <simgear/compiler.h>
 #include <simgear/structure/exception.hxx>
 #include <simgear/debug/logstream.hxx>
@@ -18,22 +20,15 @@
 #include <simgear/misc/texcoord.hxx>
 
 #include <Polygon/polygon.hxx>
-#include <Polygon/texparams.hxx>
-#include <Polygon/superpoly.hxx>
-#include <Polygon/texparams.hxx>
-#include <Polygon/chop.hxx>
-
-#include <Geometry/poly_support.hxx>
-#include <Geometry/poly_extra.hxx>
+#include <Polygon/tg_unique_geod.hxx>
+#include <Polygon/tg_unique_vec3f.hxx>
+#include <Polygon/tg_unique_vec2f.hxx>
 
 #include <Output/output.hxx>
 
 #include "global.hxx"
 #include "elevations.hxx"
 
-
-#include <stdio.h>
-
 string SGLOG_GREEN  = "\033[0;32m";
 string SGLOG_NORMAL = "\033[0m";
 
@@ -57,7 +52,7 @@ Airport::Airport( int c, char* def)
     // we need to tokenize airports, since we can't scanf two strings next to each other...
     numParams = 0;
     bool done = false;
-    
+
     while (!done)
     {
         // trim leading whitespace
@@ -232,18 +227,29 @@ bool Airport::isDebugFeature( int feat )
 // TODO : Add somewhere
 // Determine node elevations of a point_list based on the provided
 // TGAptSurface.  Offset is added to the final elevation
-static point_list calc_elevations( TGAptSurface &surf,
-                                   const point_list& geod_nodes,
-                                   double offset )
+static std::vector<SGGeod> calc_elevations( TGAptSurface &surf, const std::vector<SGGeod>& geod_nodes, double offset )
+{
+    std::vector<SGGeod> result = geod_nodes;
+    for ( unsigned int i = 0; i < result.size(); ++i ) {
+        double elev = surf.query( result[i] );
+        result[i].setElevationM( elev + offset );
+    }
+
+    return result;
+}
+
+#if 1 // unused
+static point_list calc_elevations( TGAptSurface &surf, const point_list& geod_nodes, double offset )
 {
     point_list result = geod_nodes;
     for ( unsigned int i = 0; i < result.size(); ++i ) {
-        double elev = surf.query( SGGeod::fromDeg( result[i].lon(), result[i].lat()) );
+        double elev = surf.query( SGGeod::fromDeg( result[i].x(), result[i].y() ) );
         result[i].setelev( elev + offset );
     }
 
     return result;
 }
+#endif
 
 static tgContour calc_elevations( TGAptSurface &surf,
                                   const tgContour& geod_nodes,
@@ -308,7 +314,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     // runway lights
     tglightcontour_list rwy_lights;
 
-    Point3D p;
+    //Point3D p;
 
     bool make_shapefiles = false;
 
@@ -600,8 +606,11 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
 
     // add segments to polygons to remove any possible "T"
     // intersections
-    TGTriNodes tmp_pvmt_nodes;
-    TGTriNodes tmp_feat_nodes;
+    //TGTriNodes tmp_pvmt_nodes;
+    //TGTriNodes tmp_feat_nodes;
+
+    UniqueSGGeodSet tmp_pvmt_nodes;
+    UniqueSGGeodSet tmp_feat_nodes;
 
     // build temporary node list from runways...
     SG_LOG(SG_GENERAL, SG_INFO, "Build Node List " );
@@ -612,7 +621,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         {
     	    for ( unsigned int j = 0; j < poly.ContourSize( i ); ++j )
             {
-                tmp_pvmt_nodes.unique_add( Point3D::fromSGGeod( poly.GetNode(i, j) ) );
+                tmp_pvmt_nodes.add( poly.GetNode(i, j) );
     	    }
     	}
     }
@@ -625,7 +634,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         {
     	    for ( unsigned int j = 0; j < poly.ContourSize( i ); ++j )
             {
-        		tmp_pvmt_nodes.unique_add( Point3D::fromSGGeod( poly.GetNode(i, j) ) );
+        		tmp_pvmt_nodes.add( poly.GetNode(i, j) );
     	    }
     	}
     }
@@ -638,7 +647,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         {
     	    for ( unsigned int j = 0; j < poly.ContourSize( i ); ++j )
             {
-        		tmp_feat_nodes.unique_add( Point3D::fromSGGeod( poly.GetNode(i, j) ) );
+        		tmp_feat_nodes.add( poly.GetNode(i, j) );
     	    }
     	}
     }
@@ -648,7 +657,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     {
     	for ( unsigned int j = 0; j < base_poly.ContourSize( i ); ++j )
         {
-    	    tmp_pvmt_nodes.unique_add( Point3D::fromSGGeod( base_poly.GetNode(i, j) ) );
+    	    tmp_pvmt_nodes.add( base_poly.GetNode(i, j) );
     	}
     }
 
@@ -658,7 +667,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     {
         for ( unsigned int j = 0; j < divided_base.ContourSize( i ); ++j )
         {
-            tmp_pvmt_nodes.unique_add( Point3D::fromSGGeod( divided_base.GetNode(i, j) ) );
+            tmp_pvmt_nodes.add( divided_base.GetNode(i, j) );
         }
     }
 
@@ -669,7 +678,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     for ( unsigned int k = 0; k < rwy_polys.size(); ++k ) 
     {
         tgPolygon poly = rwy_polys[k];
-        poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes );
+        poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes.get_list() );
         SG_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes());
         rwy_polys[k] = poly;
     }
@@ -678,7 +687,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     for ( unsigned int k = 0; k < pvmt_polys.size(); ++k ) 
     {
         tgPolygon poly = pvmt_polys[k];
-        poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes );
+        poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes.get_list() );
         SG_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes());
         pvmt_polys[k] = poly;
     }
@@ -687,7 +696,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     for ( unsigned int k = 0; k < line_polys.size(); ++k ) 
     {
         tgPolygon poly = line_polys[k];
-        poly = tgPolygon::AddColinearNodes( poly, tmp_feat_nodes );
+        poly = tgPolygon::AddColinearNodes( poly, tmp_feat_nodes.get_list() );
         SG_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes());
         line_polys[k] = poly;
     }
@@ -715,7 +724,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     log_time = time(0);
     SG_LOG( SG_GENERAL, SG_ALERT, "Finished cleaning polys for " << icao << " at " << my_ctime(log_time) );
 
-    base_poly = tgPolygon::AddColinearNodes( base_poly, tmp_pvmt_nodes );
+    base_poly = tgPolygon::AddColinearNodes( base_poly, tmp_pvmt_nodes.get_list() );
     base_poly = tgPolygon::Snap( base_poly, gSnap );
 
     // Finally find slivers in base
@@ -815,11 +824,10 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
 
     // traverse the tri list and create ordered node and texture
     // coordinate lists
-
-    TGTriNodes nodes, normals, texcoords;
-    nodes.clear();
-    normals.clear();
-    texcoords.clear();
+    // start with just nodes
+    UniqueSGGeodSet  nodes;
+    UniqueSGVec3fSet normals;
+    UniqueSGVec2fSet texcoords;
 
     group_list pts_v; pts_v.clear();
     group_list pts_n; pts_n.clear();
@@ -841,9 +849,8 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     int_list tri_tc, strip_tc;
 
     // calculate "the" normal for this airport
-    SGVec3d vnt = SGVec3d::fromGeod( base_poly.GetNode(0, 0) );
+    SGVec3f vnt = SGVec3f::fromGeod( base_poly.GetNode(0, 0) );
     vnt = normalize(vnt);
-    Point3D vn = Point3D::fromSGVec3(vnt);
 
     SG_LOG(SG_GENERAL, SG_INFO, "Adding runway nodes and normals");
     for ( unsigned int k = 0; k < rwy_polys.size(); ++k ) 
@@ -860,17 +867,14 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
             // let's go out on a limb, and assume triangles have 3 points..
             for ( int j = 0; j < 3; ++j )
             {
-                p = Point3D::fromSGGeod( rwy_polys[k].GetTriNode( i, j ) );
-                index = nodes.unique_add( p );
-                SG_LOG(SG_GENERAL, SG_DEBUG, "added rwy point " << p << " at " << index );
+                index = nodes.add( rwy_polys[k].GetTriNode( i, j ) );
                 tri_v.push_back( index );
 
         		// use 'the' normal
-        		index = normals.unique_add( vn );
+        		index = normals.add( vnt );
         		tri_n.push_back( index );
 
-        		Point3D tc = Point3D::fromSGVec2( rwy_polys[k].GetTriTexCoord( i, j ) );
-        		index = texcoords.unique_add( tc );
+        		index = texcoords.add( rwy_polys[k].GetTriTexCoord(i,j) );
         		tri_tc.push_back( index );
     	    }
     	    tris_v.push_back( tri_v );
@@ -894,17 +898,14 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     	    tri_tc.clear();
     	    for ( int j = 0; j < 3; ++j )
             {
-                p = Point3D::fromSGGeod( pvmt_polys[k].GetTriNode( i, j ) );
-                index = nodes.unique_add( p );
-                SG_LOG(SG_GENERAL, SG_DEBUG, "added pvmt point " << p << " at " << index );
+                index = nodes.add( pvmt_polys[k].GetTriNode( i, j ) );
                 tri_v.push_back( index );
 
                 // use 'the' normal
-                index = normals.unique_add( vn );
+                index = normals.add( vnt );
                 tri_n.push_back( index );
 
-                Point3D tc = Point3D::fromSGVec2( pvmt_polys[k].GetTriTexCoord( i, j ) );
-                index = texcoords.unique_add( tc );
+                index = texcoords.add( pvmt_polys[k].GetTriTexCoord(i,j) );
                 tri_tc.push_back( index );
     	    }
     	    tris_v.push_back( tri_v );
@@ -928,17 +929,14 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     	    tri_tc.clear();
     	    for ( int j = 0; j < 3; ++j )
             {
-                p = Point3D::fromSGGeod( line_polys[k].GetTriNode( i, j ) );
-                index = nodes.unique_add( p );
-                SG_LOG(SG_GENERAL, SG_DEBUG, "added line point " << p << " at " << index );
+                index = nodes.add( line_polys[k].GetTriNode( i, j ) );
                 tri_v.push_back( index );
 
                 // use 'the' normal
-                index = normals.unique_add( vn );
+                index = normals.add( vnt );
                 tri_n.push_back( index );
 
-                Point3D tc = Point3D::fromSGVec2( line_polys[k].GetTriTexCoord( i, j ) );
-                index = texcoords.unique_add( tc );
+                index = texcoords.add( line_polys[k].GetTriTexCoord(i,j) );
                 tri_tc.push_back( index );
     	    }
     	    tris_v.push_back( tri_v );
@@ -960,33 +958,25 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     	tri_tc.clear();
     	for ( int j = 0; j < 3; ++j )
         {
-    	    p = Point3D::fromSGGeod( base_poly.GetTriNode( i, j ) );
-    	    index = nodes.unique_add( p );
-            SG_LOG(SG_GENERAL, SG_DEBUG, "added base point " << p << " at " << index );
+    	    index = nodes.add( base_poly.GetTriNode( i, j ) );
     	    tri_v.push_back( index );
 
-    	    index = normals.unique_add( vn );
+    	    index = normals.add( vnt );
     	    tri_n.push_back( index);
     	}
     	tris_v.push_back( tri_v );
     	tris_n.push_back( tri_n );
     	tri_materials.push_back( "Grass" );
 
-    	std::vector < SGGeod > geodNodes;
-    	for ( unsigned int j = 0; j < nodes.get_node_list().size(); j++ ) 
-        {
-    	    Point3D node = nodes.get_node_list()[j];
-    	    geodNodes.push_back( SGGeod::fromDegM( node.x(), node.y(), node.z() ) );
-    	}
-	    base_txs.clear();
+    	std::vector < SGGeod > geodNodes = nodes.get_list();
+
+        base_txs.clear();
     	base_txs = sgCalcTexCoords( b, geodNodes, tri_v );
 
     	base_tc.clear();
     	for ( unsigned int j = 0; j < base_txs.size(); ++j ) 
         {
-    	    SGVec2f tc = base_txs[j];
-    	    // SG_LOG(SG_GENERAL, SG_DEBUG, "base_tc = " << tc);
-    	    index = texcoords.simple_add( Point3D( tc.x(), tc.y(), 0 ) );
+    	    index = texcoords.add(  base_txs[j] );
     	    base_tc.push_back( index );
     	}
     	tris_tc.push_back( base_tc );
@@ -999,7 +989,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     {
     	for ( unsigned int j = 0; j < divided_base.ContourSize( i ); ++j )
         {
-    	    index = nodes.unique_add( Point3D::fromSGGeod( divided_base.GetNode(i, j) ) );
+    	    index = nodes.add( divided_base.GetNode(i, j) );
             SG_LOG(SG_GENERAL, SG_DEBUG, "added base point " << divided_base.GetNode(i, j) << " at " << index );
     	}
     }
@@ -1013,7 +1003,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
 
     SG_LOG(SG_GENERAL, SG_DEBUG, " calc average elevation");
     {
-        point_list dbg = nodes.get_node_list();
+        std::vector < SGGeod > dbg = nodes.get_list();
 
         // dump the node list
         SG_LOG(SG_GENERAL, SG_DEBUG, " node list size is " << dbg.size() );
@@ -1023,7 +1013,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         }
     }
 
-    double average = tgAverageElevation( root, elev_src, nodes.get_node_list() );
+    double average = tgAverageElevationGeod( root, elev_src, nodes.get_list() );
 
     // Now build the fitted airport surface ...
 
@@ -1032,28 +1022,28 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
 
     SGGeod min_deg = SGGeod::fromDeg(9999.0, 9999.0);
     SGGeod max_deg = SGGeod::fromDeg(-9999.0, -9999.0);
-    for ( unsigned int j = 0; j < nodes.get_node_list().size(); ++j ) 
+    for ( unsigned int j = 0; j < nodes.get_list().size(); ++j )
     {
-        Point3D p = nodes.get_node_list()[j];
-        if ( p.lon() < min_deg.getLongitudeDeg() )
+        SGGeod p = nodes.get_list()[j];
+        if ( p.getLongitudeDeg() < min_deg.getLongitudeDeg() )
         {
-            SG_LOG(SG_GENERAL, SG_DEBUG, "new min lon from node " << j << " is " << p.lon() );
-            min_deg.setLongitudeDeg( p.lon() );
+            SG_LOG(SG_GENERAL, SG_DEBUG, "new min lon from node " << j << " is " << p.getLongitudeDeg() );
+            min_deg.setLongitudeDeg( p.getLongitudeDeg() );
         }
-        if ( p.lon() > max_deg.getLongitudeDeg() )
+        if ( p.getLongitudeDeg() > max_deg.getLongitudeDeg() )
         {
-            SG_LOG(SG_GENERAL, SG_DEBUG, "new max lon from node " << j << " is " << p.lon() );
-            max_deg.setLongitudeDeg( p.lon() );
+            SG_LOG(SG_GENERAL, SG_DEBUG, "new max lon from node " << j << " is " << p.getLongitudeDeg() );
+            max_deg.setLongitudeDeg( p.getLongitudeDeg() );
         }
-        if ( p.lat() < min_deg.getLatitudeDeg() )
+        if ( p.getLatitudeDeg() < min_deg.getLatitudeDeg() )
         {
-            SG_LOG(SG_GENERAL, SG_DEBUG, "new min lat from node " << j << " is " << p.lat() );
-            min_deg.setLatitudeDeg( p.lat() );
+            SG_LOG(SG_GENERAL, SG_DEBUG, "new min lat from node " << j << " is " << p.getLatitudeDeg() );
+            min_deg.setLatitudeDeg( p.getLatitudeDeg() );
         }
-        if ( p.lat() > max_deg.getLatitudeDeg() )
+        if ( p.getLatitudeDeg() > max_deg.getLatitudeDeg() )
         {
-            SG_LOG(SG_GENERAL, SG_DEBUG, "new max lat from node " << j << " is " << p.lat() );
-            max_deg.setLatitudeDeg( p.lat() );
+            SG_LOG(SG_GENERAL, SG_DEBUG, "new max lat from node " << j << " is " << p.getLatitudeDeg() );
+            max_deg.setLatitudeDeg( p.getLatitudeDeg() );
         }
     }
 
@@ -1104,28 +1094,18 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     SG_LOG(SG_GENERAL, SG_DEBUG, " min: " << min_deg );
     SG_LOG(SG_GENERAL, SG_DEBUG, " max: " << max_deg );
     SG_LOG(SG_GENERAL, SG_DEBUG, " average: " << average );
-    
+
     tg::Rectangle aptBounds(min_deg, max_deg);
 
     TGAptSurface apt_surf( root, elev_src, aptBounds, average );
     SG_LOG(SG_GENERAL, SG_DEBUG, "Airport surface created");
 
-    // calculate node elevations
-    SG_LOG(SG_GENERAL, SG_DEBUG, "Computing airport node elevations");
-    point_list geod_nodes = calc_elevations( apt_surf,
-                                             nodes.get_node_list(),
-                                             0.0 );
-    divided_base = calc_elevations( apt_surf, divided_base, 0.0 );
-
-    SG_LOG(SG_GENERAL, SG_DEBUG, "Done with base calc_elevations()");
-
-
     // add base skirt (to hide potential cracks)
     // this has to happen after we've calculated the node elevations
     // but before we convert to wgs84 coordinates
 
+#if 0 // do we still need a skirt?
     int uindex, lindex;
-
     for ( unsigned int i = 0; i < divided_base.Contours(); ++i )
     {
         strip_v.clear();
@@ -1133,8 +1113,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         strip_tc.clear();
 
         // prime the pump ...
-        p = Point3D::fromSGGeod( divided_base.GetNode( i, 0 ) );
-        uindex = nodes.find( p );
+        uindex = nodes.find( divided_base.GetNode( i, 0 ) );
         if ( uindex >= 0 )
         {
             Point3D lower = geod_nodes[uindex] - Point3D(0, 0, 20);
@@ -1222,89 +1201,59 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         base_tc.clear();
         for ( unsigned int j = 0; j < base_txs.size(); ++j )
         {
-            SGVec2f tc = base_txs[j];
-            // SG_LOG(SG_GENERAL, SG_DEBUG, "base_tc = " << tc);
-            index = texcoords.simple_add( Point3D( tc.x(), tc.y(), 0 ) );
+            index = texcoords.simple_add( base_txs[j] );
             base_tc.push_back( index );
         }
         strips_tc.push_back( base_tc );
     }
+#endif
 
     // add light points
-    superpoly_list tmp_light_list; 
-    tmp_light_list.clear();
-    typedef std::map < string, double, std::less<string> > elev_map_type;
-    typedef elev_map_type::const_iterator const_elev_map_iterator;
-    elev_map_type elevation_map;
-
-    SG_LOG(SG_GENERAL, SG_DEBUG, "Computing runway/approach lighting elevations");
-
     // pass one, calculate raw elevations from Array
-    for ( unsigned int i = 0; i < rwy_lights.size(); ++i ) 
-    {
-        TGTriNodes light_nodes;
-        light_nodes.clear();
-        point_list lights_v = rwy_lights[i].TempGetPosListAsPoint3D();
-        for ( unsigned int j = 0; j < lights_v.size(); ++j ) 
-        {
-            p = lights_v[j];
-            index = light_nodes.simple_add( p );
+    for ( unsigned int i = 0; i < rwy_lights.size(); ++i ) {
+        for ( unsigned int j = 0; j < rwy_lights[i].ContourSize(); j++ ) {
+            double light_elevation = calc_elevation( apt_surf, rwy_lights[i].GetNode(j), 0.0 );
+            rwy_lights[i].SetElevation(j, light_elevation);
         }
 
-        // calculate light node elevations
-        point_list geod_light_nodes = calc_elevations( apt_surf, light_nodes.get_node_list(), 0.0 );
-        TGPolygon p;
-        p.add_contour( geod_light_nodes, 0 );
-        TGSuperPoly s;
-        s.set_poly( p );
-        tmp_light_list.push_back( s );
-
-        string flag = rwy_lights[i].GetFlag();
-        if ( flag != (string)"" ) 
-        {
-            double max = -9999;
-            const_elev_map_iterator it = elevation_map.find( flag );
-            if ( it != elevation_map.end() ) 
-            {
-                max = elevation_map[flag];
-            }
-            for ( unsigned int j = 0; j < geod_light_nodes.size(); ++j ) 
-            {
-                if ( geod_light_nodes[j].z() > max ) 
-                {
-                    max = geod_light_nodes[j].z();
-                }
-            }
-            elevation_map[flag] = max;
-            SG_LOG( SG_GENERAL, SG_DEBUG, flag << " max = " << max );
-        }
+        // TODO : It doesn't look like this does anything... got back in history and make sure...
+        //string flag = rwy_lights[i].GetFlag();
+        //if ( flag != (string)"" ) 
+        //{
+        //    SG_LOG(SG_GENERAL, SG_INFO, "    flag " << flag);
+        //    double max = -9999;
+        //    const_elev_map_iterator it = elevation_map.find( flag );
+        //    if ( it != elevation_map.end() ) 
+        //    {
+        //       max = elevation_map[flag];
+        //    }
+        //    for ( unsigned int j = 0; j < geod_light_nodes.size(); ++j ) 
+        //    {
+        //        if ( geod_light_nodes[j].z() > max ) 
+        //        {
+        //            max = geod_light_nodes[j].z();
+        //        }
+        //    }
+        //    elevation_map[flag] = max;
+        //    SG_LOG( SG_GENERAL, SG_INFO, "      " << flag << " max = " << max );
+        //}
     }
 
-    SG_LOG(SG_GENERAL, SG_DEBUG, "Done with lighting calc_elevations() num light polys is " << rwy_lights.size() );
+    SG_LOG(SG_GENERAL, SG_INFO, "Done with lighting calc_elevations() num light polys is " << rwy_lights.size() );
 
     // pass two, for each light group check if we need to lift (based
     // on flag) and do so, then output next structures.
-    for ( unsigned int i = 0; i < rwy_lights.size(); ++i ) 
+    // for ( unsigned int i = 0; i < rwy_lights.size(); ++i )
+    for ( unsigned int i = 0; i < rwy_lights.size(); ++i )
     {
-        // tmp_light_list is a parallel structure to rwy_lights
-        point_list geod_light_nodes = tmp_light_list[i].get_poly().get_contour(0);
-        
-        // this is a little round about, but what we want to calculate the
-        // light node elevations as ground + an offset so we do them
-        // seperately, then we add them back into nodes to get the index
-        // out, but also add them to geod_nodes to maintain consistancy
-        // between these two lists.
-        point_list light_normals = rwy_lights[i].TempGetNormalListAsPoint3D();
         pt_v.clear();
         pt_n.clear();
-        for ( unsigned int j = 0; j < geod_light_nodes.size(); ++j ) 
+        for ( unsigned int j = 0; j < rwy_lights[i].ContourSize(); ++j )
         {
-            p = geod_light_nodes[j];
-            index = nodes.simple_add( p );
+            index = nodes.add( rwy_lights[i].GetPosition(j) );
             pt_v.push_back( index );
-            geod_nodes.push_back( p );
 
-            index = normals.unique_add( light_normals[j] );
+            index = normals.add( rwy_lights[i].GetNormal(j) );
             pt_n.push_back( index );
         }
         pts_v.push_back( pt_v );
@@ -1312,22 +1261,25 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
         pt_materials.push_back( rwy_lights[i].GetType() );
     }
 
+    // calculate node elevations
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Computing airport node elevations");
+
+    std::vector<SGGeod> geod_nodes = calc_elevations( apt_surf, nodes.get_list(), 0.0 );
+    divided_base = calc_elevations( apt_surf, divided_base, 0.0 );
+
     // calculate wgs84 mapping of nodes
-    std::vector< SGVec3d > wgs84_nodes;
+    std::vector<SGVec3d> wgs84_nodes;
     for ( unsigned int i = 0; i < geod_nodes.size(); ++i ) 
     {
-        SGGeod geod = SGGeod::fromDegM( geod_nodes[i].x(), geod_nodes[i].y(), geod_nodes[i].z() );
-    	SG_LOG(SG_GENERAL, SG_DEBUG, "geod pt = " << geod_nodes[i] );
-        SGVec3d cart = SGVec3d::fromGeod(geod);
-        SG_LOG(SG_GENERAL, SG_DEBUG, "  cart pt = " << cart );
-    	wgs84_nodes.push_back( cart );
+        wgs84_nodes.push_back( SGVec3d::fromGeod(geod_nodes[i]) );
     }
+
     SGSphered d;
     for ( unsigned int i = 0; i < wgs84_nodes.size(); ++i ) 
     {
         d.expandBy(wgs84_nodes[ i ]);
     }
-    
+
     SGVec3d gbs_center = d.getCenter();
     double gbs_radius = d.getRadius();
     SG_LOG(SG_GENERAL, SG_DEBUG, "gbs center = " << gbs_center);
@@ -1343,27 +1295,12 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src )
     string objpath = root + "/AirportObj";
     string name = icao + ".btg";
 
-    std::vector< SGVec3f > normals_3f;
-    for ( unsigned int i = 0; i < normals.get_node_list().size(); i++ ) 
-    {
-        Point3D node = normals.get_node_list()[i];
-        normals_3f.push_back( node.toSGVec3f() );
-    }
-
-    std::vector< SGVec2f > texcoords_2f;
-    for ( unsigned int i = 0; i < texcoords.get_node_list().size(); i++ ) 
-    {
-        Point3D node = texcoords.get_node_list()[i];
-        texcoords_2f.push_back( node.toSGVec2f() );
-    }
-
     SGBinObject obj;
-
     obj.set_gbs_center( gbs_center );
     obj.set_gbs_radius( gbs_radius );
     obj.set_wgs84_nodes( wgs84_nodes );
-    obj.set_normals( normals_3f );
-    obj.set_texcoords( texcoords_2f );
+    obj.set_normals( normals.get_list() );
+    obj.set_texcoords( texcoords.get_list() );
     obj.set_pts_v( pts_v );
     obj.set_pts_n( pts_n );
     obj.set_pt_materials( pt_materials );
diff --git a/src/Airports/GenAirports850/elevations.cxx b/src/Airports/GenAirports850/elevations.cxx
index cefa957a..2935de39 100644
--- a/src/Airports/GenAirports850/elevations.cxx
+++ b/src/Airports/GenAirports850/elevations.cxx
@@ -129,6 +129,98 @@ double tgAverageElevation( const std::string &root, const string_list elev_src,
     return average;
 }
 
+double tgAverageElevationGeod( const std::string &root, const string_list elev_src,
+                               const std::vector<SGGeod> points_source )
+{
+    bool done = false;
+    unsigned int i;
+    TGArray array;
+
+    // make a copy so our routine is non-destructive.
+    std::vector<SGGeod> points = points_source;
+
+    // just bail if no work to do
+    if ( points.size() == 0 ) {
+        return 0.0;
+    }
+
+    // set all elevations to -9999
+    for ( i = 0; i < points.size(); ++i ) {
+        points[i].setElevationM( -9999.0 );
+    }
+
+    while ( !done ) {
+        // find first node with -9999 elevation
+        SGGeod first = SGGeod();
+        bool found_one = false;
+        for ( i = 0; i < points.size(); ++i ) {
+            if ( points[i].getElevationM() < -9000.0 && !found_one ) {
+                first = points[i];
+                SG_LOG( SG_GENERAL, SG_DEBUG, "found first = " << first );
+
+                found_one = true;
+            }
+        }
+
+        if ( found_one ) {
+            SGBucket b( first );
+            std::string base = b.gen_base_path();
+
+            // try the various elevation sources
+            i = 0;
+            bool found_file = false;
+            while ( !found_file && i < elev_src.size() ) {
+                std::string array_path = root + "/" + elev_src[i] + "/" + base + "/" + b.gen_index_str();
+
+                if ( array.open(array_path) ) {
+                    found_file = true;
+                    SG_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = " << array_path );
+                }
+                i++;
+            }
+
+            // this will fill in a zero structure if no array data
+            // found/opened
+            array.parse( b );
+
+            // this will do a hasty job of removing voids by inserting
+            // data from the nearest neighbor (sort of)
+            array.remove_voids();
+
+            // update all the non-updated elevations that are inside
+            // this array file
+            double elev;
+            done = true;
+            for ( i = 0; i < points.size(); ++i ) {
+                if ( points[i].getElevationM() < -9000.0 ) {
+                    done = false;
+                    elev = array.altitude_from_grid( points[i].getLongitudeDeg() * 3600.0,
+                                                     points[i].getLatitudeDeg() * 3600.0 );
+                    if ( elev > -9000 ) {
+                        points[i].setElevationM( elev );
+                    }
+                }
+            }
+
+            array.close();
+        } else {
+            done = true;
+        }
+    }
+
+    // now find the average height of the queried points
+    double total = 0.0;
+    int count = 0;
+    for ( i = 0; i < points.size(); ++i ) {
+        total += points[i].getElevationM();
+        count++;
+    }
+    double average = total / (double) count;
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of point list = " << average);
+
+    return average;
+}
+
 
 // lookup node elevations for each point in the specified simple
 // matrix.  Returns average of all points.
diff --git a/src/Airports/GenAirports850/elevations.hxx b/src/Airports/GenAirports850/elevations.hxx
index 3226f7cb..31aae2f0 100644
--- a/src/Airports/GenAirports850/elevations.hxx
+++ b/src/Airports/GenAirports850/elevations.hxx
@@ -29,6 +29,10 @@
 double tgAverageElevation( const std::string &root, const string_list elev_src,
                            const point_list points_source );
 
+// same - with a geod list
+double tgAverageElevationGeod( const std::string &root, const string_list elev_src,
+                           const std::vector<SGGeod> points_source );
+
 // lookup node elevations for each point in the specified nurbs++
 // matrix.
 void tgCalcElevations( const std::string &root, const string_list elev_src, SimpleMatrix &Pts, double average );
diff --git a/src/Airports/GenAirports850/helipad.cxx b/src/Airports/GenAirports850/helipad.cxx
index df89a3db..77404114 100644
--- a/src/Airports/GenAirports850/helipad.cxx
+++ b/src/Airports/GenAirports850/helipad.cxx
@@ -44,7 +44,7 @@ tglightcontour_list Helipad::gen_helipad_lights(double maxsize) {
     tglightcontour_list result;
 
     // Vector calculation
-    SGVec3d vec = normalize(SGVec3d::fromGeod(GetLoc()));
+    SGVec3f vec = normalize(SGVec3f::fromGeod(GetLoc()));
 
     // Create yellow edge lights, 5m spacing
     tgContour area = gen_helipad_area_w_extend(0.0, 0.0);
diff --git a/src/Airports/GenAirports850/lights.cxx b/src/Airports/GenAirports850/lights.cxx
index be7f797a..fc9edbf2 100644
--- a/src/Airports/GenAirports850/lights.cxx
+++ b/src/Airports/GenAirports850/lights.cxx
@@ -34,7 +34,7 @@ using std::string;
 
 // calculate the runway light direction vector.  We take both runway
 // ends to get the direction of the runway.
-SGVec3d Runway::gen_runway_light_vector( float angle, bool recip ) {
+SGVec3f Runway::gen_runway_light_vector( float angle, bool recip ) {
     SGVec3f cart1, cart2;
     if ( !recip ) {
         cart1 = normalize(SGVec3f::fromGeod(GetStart()));
@@ -49,7 +49,7 @@ SGVec3d Runway::gen_runway_light_vector( float angle, bool recip ) {
     SGQuatf rotation = SGQuatf::fromAngleAxisDeg(angle, horizontal);
     SGVec3f light_vec = rotation.transform(runway_vec);
 
-    return (SGVec3d)light_vec;
+    return light_vec;
 }
 
 // generate runway edge lighting
@@ -68,7 +68,7 @@ tglightcontour_list Runway::gen_runway_edge_lights( bool recip )
     double step = dist / divs;
     SGGeod pt1, pt2;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     if ( recip ) {
         length_hdg = SGMiscd::normalizePeriodic(0, 360, rwy.heading + 180.0);
@@ -165,8 +165,8 @@ tglightcontour_list Runway::gen_runway_threshold_lights( const int kind, bool re
         ref2 = GetStart();
     }
 
-    SGVec3d normal1 = gen_runway_light_vector( 3.0, recip );
-    SGVec3d normal2 = gen_runway_light_vector( 3.0, !recip );
+    SGVec3f normal1 = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal2 = gen_runway_light_vector( 3.0, !recip );
 
     double left_hdg = SGMiscd::normalizePeriodic(0, 360, length_hdg - 90.0);
     int divs = (int)(rwy.width + 4) / 3.0;
@@ -244,7 +244,7 @@ tglightcontour_list Runway::gen_runway_center_line_lights( bool recip )
     tgLightContour r_lights;
     tglightcontour_list result;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     SGGeod pt1, pt2;
     double length_hdg;
@@ -301,7 +301,7 @@ tgLightContour Runway::gen_touchdown_zone_lights( bool recip )
 {
     tgLightContour lights;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     // determine the start point.
     SGGeod ref;
@@ -360,11 +360,11 @@ tgLightContour Runway::gen_reil( const int kind, bool recip )
 {
     tgLightContour lights;
     string flag = rwy.rwnum[get_thresh0(recip)];
-    SGVec3d normal;
+    SGVec3f normal;
 
     if (kind == 1) {
         // omnidirectional lights
-        normal = normalize(SGVec3d::fromGeod(GetStart()));
+        normal = normalize(SGVec3f::fromGeod(GetStart()));
     } else {
         // unidirectional lights
         normal = gen_runway_light_vector( 10, recip );
@@ -405,7 +405,7 @@ tglightcontour_list Runway::gen_calvert( const string &kind, bool recip )
     int i, j;
     string flag;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     // Generate long center bar of lights
     // determine the start point.
@@ -583,7 +583,7 @@ tglightcontour_list Runway::gen_alsf( const string &kind, bool recip )
     int i;
     string flag;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     // Generate long center bar of lights
     // determine the start point.
@@ -848,12 +848,12 @@ tgLightContour Runway::gen_odals( const int kind, bool recip )
 
     int i;
     string material;
-    SGVec3d normal;
+    SGVec3f normal;
 
     if (kind == 0) {
         // ODALS lighting is omni-directional, but we generate a normal as
         // a placeholder to keep everything happy.
-        normal = normalize(SGVec3d::fromGeod(GetStart()));
+        normal = normalize(SGVec3f::fromGeod(GetStart()));
         material = "RWY_ODALS_LIGHTS";
     } else {
         // RAIL lights are directional
@@ -906,7 +906,7 @@ tglightcontour_list Runway::gen_ssalx( const string& kind, bool recip )
     int i;
     string flag;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
     SGGeod ref_save, pt;
 
     // Generate long center bar of lights (every 200')
@@ -1023,7 +1023,7 @@ tglightcontour_list Runway::gen_malsx( const string& kind, bool recip )
     int i;
     string flag;
 
-    SGVec3d normal = gen_runway_light_vector( 3.0, recip );
+    SGVec3f normal = gen_runway_light_vector( 3.0, recip );
 
     // Generate long center bar of lights (every 60m)
     // determine the start point.
diff --git a/src/Airports/GenAirports850/linearfeature.cxx b/src/Airports/GenAirports850/linearfeature.cxx
index 9e4588de..050c1008 100644
--- a/src/Airports/GenAirports850/linearfeature.cxx
+++ b/src/Airports/GenAirports850/linearfeature.cxx
@@ -750,17 +750,17 @@ int LinearFeature::Finish( bool closed, unsigned int idx )
                         tmp = SGGeodesy::direct( prev_outer, heading, cur_light_dist );
                     }
 
-                    SGVec3d vec;
+                    SGVec3f vec;
                     if ( !directional_light )
                     {
                         // calculate the omnidirectional normal
-                        vec = normalize(SGVec3d::fromGeod(tmp));
+                        vec = normalize(SGVec3f::fromGeod(tmp));
                     } else
                     {
                         // calculate the directional normal
                         double heading_vec = SGMiscd::normalizePeriodic( 0, 360, heading + 90.0 );
-                        SGVec3d cart1 = SGVec3d::fromGeod(tmp);
-                        SGVec3d cart2 = SGVec3d::fromGeod( SGGeodesy::direct( tmp, heading_vec, 10 ) );
+                        SGVec3f cart1 = SGVec3f::fromGeod(tmp);
+                        SGVec3f cart2 = SGVec3f::fromGeod( SGGeodesy::direct( tmp, heading_vec, 10 ) );
                         vec = normalize(cart2 - cart1);
                     }
 
diff --git a/src/Airports/GenAirports850/main.cxx b/src/Airports/GenAirports850/main.cxx
index a6a98b22..d1b07db0 100644
--- a/src/Airports/GenAirports850/main.cxx
+++ b/src/Airports/GenAirports850/main.cxx
@@ -292,6 +292,8 @@ int main(int argc, char **argv)
     	}
     }
 
+    string airportareadir=work_dir+"/AirportArea";
+
     // check for output redirect
     if ( redirect_port >= 0 ) {
 
@@ -307,25 +309,37 @@ int main(int argc, char **argv)
     } else {
         // this is the main program -
         SG_LOG(SG_GENERAL, SG_INFO, "Launch command was " << argv[0] );
+        SG_LOG(SG_GENERAL, SG_INFO, "Run genapts with " << num_threads << " threads" );
+        SG_LOG(SG_GENERAL, SG_INFO, "Input file = " << input_file);
+        SG_LOG(SG_GENERAL, SG_INFO, "Terrain sources = ");
+        for ( unsigned int i = 0; i < elev_src.size(); ++i ) {
+            SG_LOG(SG_GENERAL, SG_INFO, "  " << work_dir << "/" << elev_src[i] );
+        }
+        SG_LOG(SG_GENERAL, SG_INFO, "Work directory = " << work_dir);
+        SG_LOG(SG_GENERAL, SG_INFO, "Nudge = " << nudge);
+        SG_LOG(SG_GENERAL, SG_INFO, "Longitude = " << min.getLongitudeDeg() << ':' << max.getLongitudeDeg());
+        SG_LOG(SG_GENERAL, SG_INFO, "Latitude = " << min.getLatitudeDeg() << ':' << max.getLatitudeDeg());
+
+        if (!max.isValid() || !min.isValid())
+        {
+            SG_LOG(SG_GENERAL, SG_ALERT, "Bad longitude or latitude");
+            exit(1);
+        }
+
+        // make work directory
+        SG_LOG(SG_GENERAL, SG_INFO, "Creating AirportArea directory");
+
+        SGPath sgp( airportareadir );
+        sgp.append( "dummy" );
+        sgp.create_dir( 0755 );
     }
 
-    SG_LOG(SG_GENERAL, SG_INFO, "Run genapts with " << num_threads << " threads" );
-    SG_LOG(SG_GENERAL, SG_INFO, "Input file = " << input_file);
-    SG_LOG(SG_GENERAL, SG_INFO, "Terrain sources = ");
-    for ( unsigned int i = 0; i < elev_src.size(); ++i ) 
-    {
-        SG_LOG(SG_GENERAL, SG_INFO, "  " << work_dir << "/" << elev_src[i] );
-    }
-    SG_LOG(SG_GENERAL, SG_INFO, "Work directory = " << work_dir);
-    SG_LOG(SG_GENERAL, SG_INFO, "Nudge = " << nudge);
-    SG_LOG(SG_GENERAL, SG_INFO, "Longitude = " << min.getLongitudeDeg() << ':' << max.getLongitudeDeg());
-    SG_LOG(SG_GENERAL, SG_INFO, "Latitude = " << min.getLatitudeDeg() << ':' << max.getLatitudeDeg());
+    string command = argv[0];
+    string lastaptfile = work_dir+"/last_apt";
 
-    if (!max.isValid() || !min.isValid())
-    {
-        SG_LOG(SG_GENERAL, SG_ALERT, "Bad longitude or latitude");
-    	exit(1);
-    }
+    // initialize persistant polygon counter
+    string counter_file = airportareadir+"/poly_counter";
+    tgPolygon_index_init( counter_file );
 
     tg::Rectangle boundingBox(min, max);
     boundingBox.sanify();
@@ -343,21 +357,6 @@ int main(int argc, char **argv)
     	exit(-1);
     }
 
-    // make work directory
-    SG_LOG(SG_GENERAL, SG_INFO, "Creating AirportArea directory");
-
-    string airportareadir=work_dir+"/AirportArea";
-    SGPath sgp( airportareadir );
-    sgp.append( "dummy" );
-    sgp.create_dir( 0755 );
-
-    string command = argv[0];
-    string lastaptfile = work_dir+"/last_apt";
-
-    // initialize persistant polygon counter
-    string counter_file = airportareadir+"/poly_counter";
-    tgPolygon_index_init( counter_file );
-
     // Initialize shapefile support (for debugging)
     tgShapefileInit();
 
diff --git a/src/Airports/GenAirports850/object.cxx b/src/Airports/GenAirports850/object.cxx
index b1e90c22..ba59e52e 100644
--- a/src/Airports/GenAirports850/object.cxx
+++ b/src/Airports/GenAirports850/object.cxx
@@ -20,9 +20,9 @@ void LightingObj::BuildBtg( tglightcontour_list& lights )
     // Calculate the light normal with the help of a second point
     // in the object heading direction.
     // SimGear takes care of the glideslope angle calculation from the normal
-    SGVec3d cart1 = SGVec3d::fromGeod(SGGeodesy::direct( ref, heading, 10));
-    SGVec3d cart2 = SGVec3d::fromGeod(ref);
-    SGVec3d normal = normalize(cart2 - cart1);
+    SGVec3f cart1 = SGVec3f::fromGeod(SGGeodesy::direct( ref, heading, 10));
+    SGVec3f cart2 = SGVec3f::fromGeod(ref);
+    SGVec3f normal = normalize(cart2 - cart1);
 
     // We know our normal, now create the lights
     SGGeod pt1;
diff --git a/src/Airports/GenAirports850/runway.hxx b/src/Airports/GenAirports850/runway.hxx
index 83f7641e..2690fa8c 100644
--- a/src/Airports/GenAirports850/runway.hxx
+++ b/src/Airports/GenAirports850/runway.hxx
@@ -164,7 +164,7 @@ private:
         return (rwy.threshold[get_thresh0(recip)] > 60.0) ? true : false;
     }
 
-    SGVec3d             gen_runway_light_vector( float angle, bool recip );
+    SGVec3f             gen_runway_light_vector( float angle, bool recip );
     tglightcontour_list gen_runway_edge_lights( bool recip );
     tglightcontour_list gen_runway_threshold_lights( const int kind, bool recip );
     tglightcontour_list gen_runway_center_line_lights( bool recip );
diff --git a/src/Airports/GenAirports850/scheduler.cxx b/src/Airports/GenAirports850/scheduler.cxx
index e13ef540..e11e984b 100644
--- a/src/Airports/GenAirports850/scheduler.cxx
+++ b/src/Airports/GenAirports850/scheduler.cxx
@@ -810,7 +810,7 @@ long Scheduler::FindAirport( string icao )
 
 void Scheduler::RetryAirport( AirportInfo* pai )
 {
-    retryList.push_back( *pai );
+    // retryList.push_back( *pai );
 }
 
 bool Scheduler::AddAirports( long start_pos, tg::Rectangle* boundingBox )
diff --git a/src/Lib/Polygon/polygon.cxx b/src/Lib/Polygon/polygon.cxx
index 0fce1acc..339e9248 100644
--- a/src/Lib/Polygon/polygon.cxx
+++ b/src/Lib/Polygon/polygon.cxx
@@ -1176,6 +1176,87 @@ bool FindIntermediateNode( const SGGeod& start, const SGGeod& end,
     return found_node;
 }
 
+bool FindIntermediateNode( const SGGeod& start, const SGGeod& end,
+                           const std::vector<SGGeod>& nodes, SGGeod& result,
+                           double bbEpsilon, double errEpsilon )
+{
+    bool found_node = false;
+    double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min;
+
+    SGGeod p0 = start;
+    SGGeod p1 = end;
+
+    double xdist = fabs(p0.getLongitudeDeg() - p1.getLongitudeDeg());
+    double ydist = fabs(p0.getLatitudeDeg()  - p1.getLatitudeDeg());
+
+    x_err_min = xdist + 1.0;
+    y_err_min = ydist + 1.0;
+
+    if ( xdist > ydist ) {
+        // sort these in a sensible order
+        SGGeod p_min, p_max;
+        if ( p0.getLongitudeDeg() < p1.getLongitudeDeg() ) {
+            p_min = p0;
+            p_max = p1;
+        } else {
+            p_min = p1;
+            p_max = p0;
+        }
+
+        m = (p_min.getLatitudeDeg() - p_max.getLatitudeDeg()) / (p_min.getLongitudeDeg() - p_max.getLongitudeDeg());
+        b = p_max.getLatitudeDeg() - m * p_max.getLongitudeDeg();
+
+        for ( int i = 0; i < (int)nodes.size(); ++i ) {
+            // cout << i << endl;
+            SGGeod current = nodes[i];
+
+            if ( (current.getLongitudeDeg() > (p_min.getLongitudeDeg() + (bbEpsilon))) && (current.getLongitudeDeg() < (p_max.getLongitudeDeg() - (bbEpsilon))) ) {
+                y_err = fabs(current.getLatitudeDeg() - (m * current.getLongitudeDeg() + b));
+
+                if ( y_err < errEpsilon ) {
+                    found_node = true;
+                    if ( y_err < y_err_min ) {
+                        result = current;
+                        y_err_min = y_err;
+                    }
+                }
+            }
+        }
+    } else {
+        // sort these in a sensible order
+        SGGeod p_min, p_max;
+        if ( p0.getLatitudeDeg() < p1.getLatitudeDeg() ) {
+            p_min = p0;
+            p_max = p1;
+        } else {
+            p_min = p1;
+            p_max = p0;
+        }
+
+        m1 = (p_min.getLongitudeDeg() - p_max.getLongitudeDeg()) / (p_min.getLatitudeDeg() - p_max.getLatitudeDeg());
+        b1 = p_max.getLongitudeDeg() - m1 * p_max.getLatitudeDeg();
+
+        for ( int i = 0; i < (int)nodes.size(); ++i ) {
+            SGGeod current = nodes[i];
+
+            if ( (current.getLatitudeDeg() > (p_min.getLatitudeDeg() + (bbEpsilon))) && (current.getLatitudeDeg() < (p_max.getLatitudeDeg() - (bbEpsilon))) ) {
+
+                x_err = fabs(current.getLongitudeDeg() - (m1 * current.getLatitudeDeg() + b1));
+
+                if ( x_err < errEpsilon ) {
+                    found_node = true;
+                    if ( x_err < x_err_min ) {
+                        result = current;
+                        x_err_min = x_err;
+                    }
+                }
+            }
+        }
+    }
+
+    return found_node;
+}
+
 void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, point_list& tmp_nodes, tgContour& result, double bbEpsilon, double errEpsilon )
 {
     SGGeod new_pt;
@@ -1194,6 +1275,24 @@ void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, point_list& tmp_n
     }
 }
 
+void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, std::vector<SGGeod>& nodes, tgContour& result, double bbEpsilon, double errEpsilon )
+{
+    SGGeod new_pt;
+
+    SG_LOG(SG_GENERAL, SG_BULK, "   " << p0 << " <==> " << p1 );
+
+    bool found_extra = FindIntermediateNode( p0, p1, nodes, new_pt, bbEpsilon, errEpsilon );
+
+    if ( found_extra ) {
+        AddIntermediateNodes( p0, new_pt, nodes, result, bbEpsilon, errEpsilon  );
+
+        result.AddNode( new_pt );
+        SG_LOG(SG_GENERAL, SG_BULK, "    adding = " << new_pt);
+
+        AddIntermediateNodes( new_pt, p1, nodes, result, bbEpsilon, errEpsilon  );
+    }
+}
+
 #define CLIPPER_FIXEDPT (10000000000000000)
 #define CLIPPER_FIXED1M (            90090)
 
@@ -1821,6 +1920,37 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, TGTriNodes node
     return result;
 }
 
+tgContour tgContour::AddColinearNodes( const tgContour& subject, std::vector<SGGeod>& nodes )
+{
+    SGGeod p0, p1;
+    tgContour result;
+
+    for ( unsigned int n = 0; n < subject.GetSize()-1; n++ ) {
+        p0 = subject.GetNode( n );
+        p1 = subject.GetNode( n+1 );
+
+        // add start of segment
+        result.AddNode( p0 );
+
+        // add intermediate points
+        AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
+    }
+
+    p0 = subject.GetNode( subject.GetSize() - 1 );
+    p1 = subject.GetNode( 0 );
+
+    // add start of segment
+    result.AddNode( p0 );
+
+    // add intermediate points
+    AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
+
+    // maintain original hole flag setting
+    result.SetHole( subject.GetHole() );
+
+    return result;
+}
+
 std::ostream& operator<< ( std::ostream& output, const tgContour& subject )
 {
     // Save the data
@@ -2379,6 +2509,20 @@ tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, TGTriNodes& nod
     return result;
 }
 
+tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, std::vector<SGGeod>& nodes )
+{
+    tgPolygon result;
+
+    result.SetMaterial( subject.GetMaterial() );
+    result.SetTexParams( subject.GetTexParams() );
+
+    for ( unsigned int c = 0; c < subject.Contours(); c++ ) {
+        result.AddContour( tgContour::AddColinearNodes( subject.GetContour(c), nodes ) );
+    }
+
+    return result;
+}
+
 void tgPolygon::InheritElevations( const tgPolygon& source )
 {
     TGTriNodes nodes;
@@ -2441,9 +2585,9 @@ void tgPolygon::InheritElevations( const tgPolygon& source )
 void tgPolygon::Texture( void )
 {
     SGGeod  p;
-    SGVec2d t;
+    SGVec2f t;
     double  x, y;
-    double  tx, ty;
+    float   tx, ty;
 
     SG_LOG(SG_GENERAL, SG_DEBUG, "Texture Poly with material " << material << " method " << tp.method << " tpref " << tp.ref << " heading " << tp.heading );
 
@@ -2491,29 +2635,29 @@ void tgPolygon::Texture( void )
                     //
                     // 4. Map x, y point into texture coordinates
                     //
-                    double tmp;
+                    float tmp;
 
-                    tmp = x / tp.width;
-                    tx = tmp * (tp.maxu - tp.minu) + tp.minu;
+                    tmp = (float)x / (float)tp.width;
+                    tx = tmp * (float)(tp.maxu - tp.minu) + (float)tp.minu;
                     SG_LOG(SG_GENERAL, SG_DEBUG, "  (" << tx << ")");
 
                     // clip u?
                     if ( (tp.method == TG_TEX_BY_TPS_CLIPU) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) {
-                        if ( tx < tp.min_clipu ) { tx = tp.min_clipu; }
-                        if ( tx > tp.max_clipu ) { tx = tp.max_clipu; }
+                        if ( tx < (float)tp.min_clipu ) { tx = (float)tp.min_clipu; }
+                        if ( tx > (float)tp.max_clipu ) { tx = (float)tp.max_clipu; }
                     }
 
-                    tmp = y / tp.length;
-                    ty = tmp * (tp.maxv - tp.minv) + tp.minv;
+                    tmp = (float)y / (float)tp.length;
+                    ty = tmp * (float)(tp.maxv - tp.minv) + (float)tp.minv;
                     SG_LOG(SG_GENERAL, SG_DEBUG, "  (" << ty << ")");
 
                     // clip v?
                     if ( (tp.method == TG_TEX_BY_TPS_CLIPV) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) {
-                        if ( ty < tp.min_clipv ) { ty = tp.min_clipv; }
-                        if ( ty > tp.max_clipv ) { ty = tp.max_clipv; }
+                        if ( ty < (float)tp.min_clipv ) { ty = (float)tp.min_clipv; }
+                        if ( ty > (float)tp.max_clipv ) { ty = (float)tp.max_clipv; }
                     }
 
-                    t = SGVec2d( tx, ty );
+                    t = SGVec2f( tx, ty );
                     SG_LOG(SG_GENERAL, SG_DEBUG, "  (" << tx << ", " << ty << ")");
 
                     triangles[i].SetTexCoord( j, t );
diff --git a/src/Lib/Polygon/polygon.hxx b/src/Lib/Polygon/polygon.hxx
index 92a29258..e00e7560 100644
--- a/src/Lib/Polygon/polygon.hxx
+++ b/src/Lib/Polygon/polygon.hxx
@@ -379,6 +379,7 @@ public:
     static void      AddToAccumulator( const tgContour& subject );
 
     static tgContour AddColinearNodes( const tgContour& subject, TGTriNodes nodes );
+    static tgContour AddColinearNodes( const tgContour& subject, std::vector<SGGeod>& nodes );
 
     // conversions
     static ClipperLib::Polygon ToClipper( const tgContour& subject );
@@ -406,7 +407,7 @@ public:
         node_list.push_back( p1 );
         node_list.push_back( p2 );
 
-        tc_list.resize(   3, SGVec2d(0.0, 0.0) );
+        tc_list.resize(   3, SGVec2f(0.0, 0.0) );
         norm_list.resize( 3, SGVec3d(0.0, 0.0, 0.0) );
         idx_list.resize(  3, -1 );
     }
@@ -414,10 +415,10 @@ public:
     SGGeod GetNode( unsigned int i ) const {
         return node_list[i];
     }
-    SGVec2d GetTexCoord( unsigned int i ) const {
+    SGVec2f GetTexCoord( unsigned int i ) const {
         return tc_list[i];
     }
-    void SetTexCoord( unsigned int i, const SGVec2d tc ) {
+    void SetTexCoord( unsigned int i, const SGVec2f tc ) {
         tc_list[i] = tc;
     }
 
@@ -426,7 +427,7 @@ public:
 
 private:
     std::vector<SGGeod>  node_list;
-    std::vector<SGVec2d> tc_list;
+    std::vector<SGVec2f> tc_list;
     std::vector<SGVec3d> norm_list;
     std::vector<int>     idx_list;
 
@@ -512,7 +513,7 @@ public:
     void SetNode( unsigned int c, unsigned int i, const SGGeod& n ) {
         contours[c].SetNode( i, n );
     }
-    
+
     void AddNode( unsigned int c, const SGGeod& n ) {
         // Make sure we have contour c.  If we don't add it
         while( contours.size() < c+1 ) {
@@ -538,7 +539,7 @@ public:
     SGGeod GetTriNode( unsigned int c, unsigned int i ) const {
         return triangles[c].GetNode( i );
     }
-    SGVec2d GetTriTexCoord( unsigned int c, unsigned int i ) const {
+    SGVec2f GetTriTexCoord( unsigned int c, unsigned int i ) const {
         return triangles[c].GetTexCoord( i );
     }
 
@@ -622,6 +623,7 @@ public:
     static void Tesselate( const tgPolygon& subject );
 
     static tgPolygon AddColinearNodes( const tgPolygon& subject, TGTriNodes& nodes );
+    static tgPolygon AddColinearNodes( const tgPolygon& subject, std::vector<SGGeod>& nodes );
 
     static void RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers );
     static void MergeSlivers( tgpolygon_list& subjects, tgcontour_list& slivers );
@@ -645,7 +647,7 @@ class tgLight
 {
 public:
     SGGeod      pos;
-    SGVec3d     norm;
+    SGVec3f     norm;
 };
 
 typedef std::vector <tgLight>  tglight_list;
@@ -659,17 +661,29 @@ public:
         return lights.size();
     }
 
-    void AddLight( SGGeod p, SGVec3d n ) {
+    void AddLight( SGGeod p, SGVec3f n ) {
         tgLight light;
         light.pos  = p;
         light.norm = n;
         lights.push_back(light);
     }
 
+    void SetElevation( unsigned int i, double elev ) {
+        lights[i].pos.setElevationM( elev );
+    }
+
     SGGeod GetNode( unsigned int i ) const {
         return lights[i].pos;
     }
 
+    SGGeod GetPosition( unsigned int i ) const {
+        return lights[i].pos;
+    }
+
+    SGVec3f GetNormal( unsigned int i ) const {
+        return lights[i].norm;
+    }
+
     std::string GetFlag( void ) const {
         return flag;
     }
@@ -695,6 +709,16 @@ public:
         return p3dlist;
     }
 
+    std::vector<SGGeod> GetPositionList( void ) {
+        std::vector<SGGeod> positions;
+
+        for (unsigned int i=0; i<lights.size(); i++) {
+            positions.push_back( lights[i].pos );
+        }
+
+        return positions;
+    }
+
     point_list TempGetNormalListAsPoint3D( void ) {
         point_list p3dlist;
 
@@ -704,6 +728,17 @@ public:
 
         return p3dlist;
     }
+
+    std::vector<SGVec3f> GetNormalList( void ) {
+        std::vector<SGVec3f> normals;
+
+        for (unsigned int i=0; i<lights.size(); i++) {
+            normals.push_back( lights[i].norm );
+        }
+
+        return normals;
+    }
+
     // END TEMP TEMP TEMP
 
     // Friend for output
diff --git a/src/Lib/Polygon/tg_unique_geod.hxx b/src/Lib/Polygon/tg_unique_geod.hxx
new file mode 100644
index 00000000..9ec0327e
--- /dev/null
+++ b/src/Lib/Polygon/tg_unique_geod.hxx
@@ -0,0 +1,128 @@
+#include <boost/unordered_set.hpp>
+
+// Implement Unique SGGeod list
+
+// We use a hash to store the indices.  All iterators returned from find are
+// constant, because modifying a value will modify the hash - rendering the
+// set invalid.
+// Selection of the bucket is tightly linked to the key equality function.
+// If any value is considered equal to another, than the hash function MUST
+// compute the same hash for each value.
+// So close_enough_2d will remain our testo of equality.  It simply rounds to
+// 6 decimal places.  Our hash function will round to 5, so at most 10 points in the
+// same bucket.
+// We could experiment with making it so 1 point per bucket, but I'm nervous...
+
+#define PROXIMITY_MULTIPLIER (100000)
+#define PROXIMITY_EPSILON    ((double) 1 / (double)( PROXIMITY_MULTIPLIER * 10 ) )
+
+class SGGeodIndex {
+public:
+    SGGeodIndex( SGGeod g ) {
+        geod = g;
+
+        std::size_t FNV_prime;
+        std::size_t offset_basis;
+
+        switch( sizeof( std::size_t ) ) {
+            case 4: // 32 bit system
+            default:
+                FNV_prime =    16777619ULL;
+                offset_basis = 2166136261ULL;
+                break;
+
+            case 8: // 64 bit system
+                FNV_prime =    1099511628211ULL;
+                offset_basis = 14695981039346656037ULL;
+                break;
+        }
+
+        hash = (std::size_t)offset_basis;
+
+        /* only hash lon, lat - we want to detect dups in 2d only */
+        unsigned long long raw_pt[2];
+        raw_pt[0] = (unsigned long long)( round(geod.getLongitudeDeg() * PROXIMITY_MULTIPLIER) );
+        raw_pt[1] = (unsigned long long)( round(geod.getLatitudeDeg() * PROXIMITY_MULTIPLIER) );
+
+        unsigned char* it = (unsigned char*)raw_pt;
+        for ( unsigned i=0; i<sizeof( raw_pt ); i++ ) {
+            hash = hash ^ *it++;
+            hash = hash * FNV_prime;
+        }
+
+        // SG_LOG(SG_GENERAL, SG_ALERT, " GetHash x: " << raw_pt[0] << " y: " << raw_pt[1] << " = " << hash );
+    }
+
+    inline void         SetOrderedIndex( unsigned int i ) { ordered_index = i; }
+    inline unsigned int GetOrderedIndex( void ) const     { return ordered_index; }
+    inline std::size_t  GetHash( void ) const             { return hash; }
+
+    friend bool operator == (const SGGeodIndex& a, const SGGeodIndex& b);
+    friend std::ostream& operator<< ( std::ostream&, const SGGeodIndex& );
+
+private:
+    SGGeod       geod;
+    std::size_t  hash;
+    unsigned int ordered_index;
+};
+
+inline bool operator == (const SGGeodIndex& a, const SGGeodIndex& b) {
+    return (( fabs(a.geod.getLongitudeDeg() - b.geod.getLongitudeDeg()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.geod.getLatitudeDeg()  - b.geod.getLatitudeDeg()) < PROXIMITY_EPSILON ));
+}
+
+struct SGGeodIndexHash : std::unary_function<SGGeodIndex, std::size_t>
+{
+    std::size_t operator()(SGGeodIndex const& gi) const {
+        return gi.GetHash();
+    }
+};
+
+typedef boost::unordered_set<SGGeodIndex, SGGeodIndexHash> unique_geod_set;
+typedef unique_geod_set::iterator unique_geod_set_iterator;
+typedef unique_geod_set::const_iterator const_unique_geod_set_iterator;
+
+class UniqueSGGeodSet {
+public:
+    UniqueSGGeodSet() {}
+
+    unsigned int add( const SGGeod& g );
+    int          find( const SGGeod& g ) const;
+    std::vector<SGGeod>& get_list( void ) { return geod_list; }
+    void get_node_list();
+
+private:
+    unique_geod_set         index_list;
+    std::vector<SGGeod>     geod_list;
+};
+
+unsigned int UniqueSGGeodSet::add( const SGGeod& g ) {
+    unique_geod_set_iterator it;
+    SGGeodIndex lookup( g );
+
+    it = index_list.find( lookup );
+    if ( it == index_list.end() ) {
+        lookup.SetOrderedIndex( geod_list.size() );
+        index_list.insert( lookup );
+
+        geod_list.push_back(g);
+    } else {
+        lookup = *it;
+    }
+
+    return lookup.GetOrderedIndex();
+}
+
+
+int UniqueSGGeodSet::find( const SGGeod& g ) const {
+    unique_geod_set_iterator it;
+    SGGeodIndex lookup( g );
+    int index = -1;
+
+    it = index_list.find( lookup );
+    if ( it != index_list.end() ) {
+        index = it->GetOrderedIndex();
+    }
+
+    return index;
+}
\ No newline at end of file
diff --git a/src/Lib/Polygon/tg_unique_vec2f.hxx b/src/Lib/Polygon/tg_unique_vec2f.hxx
new file mode 100644
index 00000000..4159a1d1
--- /dev/null
+++ b/src/Lib/Polygon/tg_unique_vec2f.hxx
@@ -0,0 +1,127 @@
+#include <boost/unordered_set.hpp>
+
+// Implement Unique SGVec2f list
+
+// We use a hash to store the indices.  All iterators returned from find are
+// constant, because modifying a value will modify the hash - rendering the
+// set invalid.
+// Selection of the bucket is tightly linked to the key equality function.
+// If any value is considered equal to another, than the hash function MUST
+// compute the same hash for each value.
+// So close_enough_2d will remain our testo of equality.  It simply rounds to
+// 6 decimal places.  Our hash function will round to 5, so at most 10 points in the
+// same bucket.
+// We could experiment with making it so 1 point per bucket, but I'm nervous...
+
+#define PROXIMITY_MULTIPLIER (100000)
+#define PROXIMITY_EPSILON    ((double) 1 / (double)( PROXIMITY_MULTIPLIER * 10 ) )
+
+class SGVec2fIndex {
+public:
+    SGVec2fIndex( SGVec2f v ) {
+        vec = v;
+
+        std::size_t FNV_prime;
+        std::size_t offset_basis;
+
+        switch( sizeof( std::size_t ) ) {
+            case 4: // 32 bit system
+            default:
+                FNV_prime =    16777619ULL;
+                offset_basis = 2166136261ULL;
+                break;
+
+            case 8: // 64 bit system
+                FNV_prime =    1099511628211ULL;
+                offset_basis = 14695981039346656037ULL;
+                break;
+        }
+
+        hash = (std::size_t)offset_basis;
+
+        /* only hash lon, lat - we want to detect dups in 2d only */
+        unsigned long long raw_pt[2];
+        raw_pt[0] = (unsigned long long)( round(vec.x() * PROXIMITY_MULTIPLIER) );
+        raw_pt[1] = (unsigned long long)( round(vec.y() * PROXIMITY_MULTIPLIER) );
+
+        unsigned char* it = (unsigned char*)raw_pt;
+        for ( unsigned i=0; i<sizeof( raw_pt ); i++ ) {
+            hash = hash ^ *it++;
+            hash = hash * FNV_prime;
+        }
+
+        // SG_LOG(SG_GENERAL, SG_ALERT, " GetHash x: " << raw_pt[0] << " y: " << raw_pt[1] << " = " << hash );
+    }
+
+    inline void         SetOrderedIndex( unsigned int i ) { ordered_index = i; }
+    inline unsigned int GetOrderedIndex( void ) const     { return ordered_index; }
+    inline std::size_t  GetHash( void ) const             { return hash; }
+
+    friend bool operator == (const SGVec2fIndex& a, const SGVec2fIndex& b);
+    friend std::ostream& operator<< ( std::ostream&, const SGVec2fIndex& );
+
+private:
+    SGVec2f      vec;
+    std::size_t  hash;
+    unsigned int ordered_index;
+};
+
+inline bool operator == (const SGVec2fIndex& a, const SGVec2fIndex& b) {
+    return (( fabs(a.vec.x() - b.vec.x()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.vec.y() - b.vec.y()) < PROXIMITY_EPSILON ));
+}
+
+struct SGVec2fIndexHash : std::unary_function<SGVec2fIndex, std::size_t>
+{
+    std::size_t operator()(SGVec2fIndex const& vi) const {
+        return vi.GetHash();
+    }
+};
+
+typedef boost::unordered_set<SGVec2fIndex, SGVec2fIndexHash> unique_vec2f_set;
+typedef unique_vec2f_set::iterator unique_vec2f_set_iterator;
+typedef unique_vec2f_set::const_iterator const_unique_vec2f_set_iterator;
+
+class UniqueSGVec2fSet {
+public:
+    UniqueSGVec2fSet() {}
+
+    unsigned int add( const SGVec2f& v );
+    int          find( const SGVec2f& v ) const;
+    std::vector<SGVec2f>& get_list( void ) { return vector_list; }
+
+private:
+    unique_vec2f_set        index_list;
+    std::vector<SGVec2f>    vector_list;
+};
+
+unsigned int UniqueSGVec2fSet::add( const SGVec2f& v ) {
+    unique_vec2f_set_iterator it;
+    SGVec2fIndex lookup( v );
+
+    it = index_list.find( lookup );
+    if ( it == index_list.end() ) {
+        lookup.SetOrderedIndex( vector_list.size() );
+        index_list.insert( lookup );
+
+        vector_list.push_back(v);
+    } else {
+        lookup = *it;
+    }
+
+    return lookup.GetOrderedIndex();
+}
+
+
+int UniqueSGVec2fSet::find( const SGVec2f& v ) const {
+    unique_vec2f_set_iterator it;
+    SGVec2fIndex lookup( v );
+    int index = -1;
+
+    it = index_list.find( lookup );
+    if ( it != index_list.end() ) {
+        index = it->GetOrderedIndex();
+    }
+
+    return index;
+}
\ No newline at end of file
diff --git a/src/Lib/Polygon/tg_unique_vec3d.hxx b/src/Lib/Polygon/tg_unique_vec3d.hxx
new file mode 100644
index 00000000..a94e4084
--- /dev/null
+++ b/src/Lib/Polygon/tg_unique_vec3d.hxx
@@ -0,0 +1,129 @@
+#include <boost/unordered_set.hpp>
+
+// Implement Unique SGVec3d list
+
+// We use a hash to store the indices.  All iterators returned from find are
+// constant, because modifying a value will modify the hash - rendering the
+// set invalid.
+// Selection of the bucket is tightly linked to the key equality function.
+// If any value is considered equal to another, than the hash function MUST
+// compute the same hash for each value.
+// So close_enough_2d will remain our testo of equality.  It simply rounds to
+// 6 decimal places.  Our hash function will round to 5, so at most 10 points in the
+// same bucket.
+// We could experiment with making it so 1 point per bucket, but I'm nervous...
+
+#define PROXIMITY_MULTIPLIER (100000)
+#define PROXIMITY_EPSILON    ((double) 1 / (double)( PROXIMITY_MULTIPLIER * 10 ) )
+
+class SGVec3dIndex {
+public:
+    SGVec3dIndex( SGVec3d v ) {
+        vec = v;
+
+        std::size_t FNV_prime;
+        std::size_t offset_basis;
+
+        switch( sizeof( std::size_t ) ) {
+            case 4: // 32 bit system
+            default:
+                FNV_prime =    16777619ULL;
+                offset_basis = 2166136261ULL;
+                break;
+
+            case 8: // 64 bit system
+                FNV_prime =    1099511628211ULL;
+                offset_basis = 14695981039346656037ULL;
+                break;
+        }
+
+        hash = (std::size_t)offset_basis;
+
+        /* only hash lon, lat - we want to detect dups in 2d only */
+        unsigned long long raw_pt[3];
+        raw_pt[0] = (unsigned long long)( round(vec.x() * PROXIMITY_MULTIPLIER) );
+        raw_pt[1] = (unsigned long long)( round(vec.y() * PROXIMITY_MULTIPLIER) );
+        raw_pt[2] = (unsigned long long)( round(vec.z() * PROXIMITY_MULTIPLIER) );
+
+        unsigned char* it = (unsigned char*)raw_pt;
+        for ( unsigned i=0; i<sizeof( raw_pt ); i++ ) {
+            hash = hash ^ *it++;
+            hash = hash * FNV_prime;
+        }
+
+        // SG_LOG(SG_GENERAL, SG_ALERT, " GetHash x: " << raw_pt[0] << " y: " << raw_pt[1] << " = " << hash );
+    }
+
+    inline void         SetOrderedIndex( unsigned int i ) { ordered_index = i; }
+    inline unsigned int GetOrderedIndex( void ) const     { return ordered_index; }
+    inline std::size_t  GetHash( void ) const             { return hash; }
+
+    friend bool operator == (const SGVec3dIndex& a, const SGVec3dIndex& b);
+    friend std::ostream& operator<< ( std::ostream&, const SGVec3dIndex& );
+
+private:
+    SGVec3d      vec;
+    std::size_t  hash;
+    unsigned int ordered_index;
+};
+
+inline bool operator == (const SGVec3dIndex& a, const SGVec3dIndex& b) {
+    return (( fabs(a.vec.x() - b.vec.x()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.vec.y() - b.vec.y()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.vec.z() - b.vec.z()) < PROXIMITY_EPSILON ));
+}
+
+struct SGVec3dIndexHash : std::unary_function<SGVec3dIndex, std::size_t>
+{
+    std::size_t operator()(SGVec3dIndex const& vi) const {
+        return vi.GetHash();
+    }
+};
+
+typedef boost::unordered_set<SGVec3dIndex, SGVec3dIndexHash> unique_vec3d_set;
+typedef unique_vec3d_set::iterator unique_vec3d_set_iterator;
+typedef unique_vec3d_set::const_iterator const_unique_vec3d_set_iterator;
+
+class UniqueSGVec3dSet {
+public:
+    UniqueSGVec3dSet() {}
+
+    unsigned int add( const SGVec3d& v );
+    int          find( const SGVec3d& v ) const;
+    std::vector<SGVec3d>& get_list( void ) { return vector_list; }
+
+private:
+    unique_vec3d_set        index_list;
+    std::vector<SGVec3d>    vector_list;
+};
+
+unsigned int UniqueSGVec3dSet::add( const SGVec3d& v ) {
+    unique_vec3d_set_iterator it;
+    SGVec3dIndex lookup( v );
+
+    it = index_list.find( lookup );
+    if ( it == index_list.end() ) {
+        lookup.SetOrderedIndex( vector_list.size() );
+        index_list.insert( lookup );
+
+        vector_list.push_back(v);
+    } else {
+        lookup = *it;
+    }
+
+    return lookup.GetOrderedIndex();
+}
+
+
+int UniqueSGVec3dSet::find( const SGVec3d& v ) const {
+    unique_vec3d_set_iterator it;
+    SGVec3dIndex lookup( v );
+    int index = -1;
+
+    it = index_list.find( lookup );
+    if ( it != index_list.end() ) {
+        index = it->GetOrderedIndex();
+    }
+
+    return index;
+}
\ No newline at end of file
diff --git a/src/Lib/Polygon/tg_unique_vec3f.hxx b/src/Lib/Polygon/tg_unique_vec3f.hxx
new file mode 100644
index 00000000..66d2c0b7
--- /dev/null
+++ b/src/Lib/Polygon/tg_unique_vec3f.hxx
@@ -0,0 +1,129 @@
+#include <boost/unordered_set.hpp>
+
+// Implement Unique SGVec3f list
+
+// We use a hash to store the indices.  All iterators returned from find are
+// constant, because modifying a value will modify the hash - rendering the
+// set invalid.
+// Selection of the bucket is tightly linked to the key equality function.
+// If any value is considered equal to another, than the hash function MUST
+// compute the same hash for each value.
+// So close_enough_2d will remain our testo of equality.  It simply rounds to
+// 6 decimal places.  Our hash function will round to 5, so at most 10 points in the
+// same bucket.
+// We could experiment with making it so 1 point per bucket, but I'm nervous...
+
+#define PROXIMITY_MULTIPLIER (100000)
+#define PROXIMITY_EPSILON    ((double) 1 / (double)( PROXIMITY_MULTIPLIER * 10 ) )
+
+class SGVec3fIndex {
+public:
+    SGVec3fIndex( SGVec3f v ) {
+        vec = v;
+
+        std::size_t FNV_prime;
+        std::size_t offset_basis;
+
+        switch( sizeof( std::size_t ) ) {
+            case 4: // 32 bit system
+            default:
+                FNV_prime =    16777619ULL;
+                offset_basis = 2166136261ULL;
+                break;
+
+            case 8: // 64 bit system
+                FNV_prime =    1099511628211ULL;
+                offset_basis = 14695981039346656037ULL;
+                break;
+        }
+
+        hash = (std::size_t)offset_basis;
+
+        /* only hash lon, lat - we want to detect dups in 2d only */
+        unsigned long long raw_pt[3];
+        raw_pt[0] = (unsigned long long)( round(vec.x() * PROXIMITY_MULTIPLIER) );
+        raw_pt[1] = (unsigned long long)( round(vec.y() * PROXIMITY_MULTIPLIER) );
+        raw_pt[2] = (unsigned long long)( round(vec.z() * PROXIMITY_MULTIPLIER) );
+
+        unsigned char* it = (unsigned char*)raw_pt;
+        for ( unsigned i=0; i<sizeof( raw_pt ); i++ ) {
+            hash = hash ^ *it++;
+            hash = hash * FNV_prime;
+        }
+
+        // SG_LOG(SG_GENERAL, SG_ALERT, " GetHash x: " << raw_pt[0] << " y: " << raw_pt[1] << " = " << hash );
+    }
+
+    inline void         SetOrderedIndex( unsigned int i ) { ordered_index = i; }
+    inline unsigned int GetOrderedIndex( void ) const     { return ordered_index; }
+    inline std::size_t  GetHash( void ) const             { return hash; }
+
+    friend bool operator == (const SGVec3fIndex& a, const SGVec3fIndex& b);
+    friend std::ostream& operator<< ( std::ostream&, const SGVec3fIndex& );
+
+private:
+    SGVec3f      vec;
+    std::size_t  hash;
+    unsigned int ordered_index;
+};
+
+inline bool operator == (const SGVec3fIndex& a, const SGVec3fIndex& b) {
+    return (( fabs(a.vec.x() - b.vec.x()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.vec.y() - b.vec.y()) < PROXIMITY_EPSILON ) &&
+            ( fabs(a.vec.z() - b.vec.z()) < PROXIMITY_EPSILON ));
+}
+
+struct SGVec3fIndexHash : std::unary_function<SGVec3fIndex, std::size_t>
+{
+    std::size_t operator()(SGVec3fIndex const& vi) const {
+        return vi.GetHash();
+    }
+};
+
+typedef boost::unordered_set<SGVec3fIndex, SGVec3fIndexHash> unique_vec3f_set;
+typedef unique_vec3f_set::iterator unique_vec3f_set_iterator;
+typedef unique_vec3f_set::const_iterator const_unique_vec3f_set_iterator;
+
+class UniqueSGVec3fSet {
+public:
+    UniqueSGVec3fSet() {}
+
+    unsigned int add( const SGVec3f& v );
+    int          find( const SGVec3f& v ) const;
+    std::vector<SGVec3f>& get_list( void ) { return vector_list; }
+
+private:
+    unique_vec3f_set        index_list;
+    std::vector<SGVec3f>    vector_list;
+};
+
+unsigned int UniqueSGVec3fSet::add( const SGVec3f& v ) {
+    unique_vec3f_set_iterator it;
+    SGVec3fIndex lookup( v );
+
+    it = index_list.find( lookup );
+    if ( it == index_list.end() ) {
+        lookup.SetOrderedIndex( vector_list.size() );
+        index_list.insert( lookup );
+
+        vector_list.push_back(v);
+    } else {
+        lookup = *it;
+    }
+
+    return lookup.GetOrderedIndex();
+}
+
+
+int UniqueSGVec3fSet::find( const SGVec3f& v ) const {
+    unique_vec3f_set_iterator it;
+    SGVec3fIndex lookup( v );
+    int index = -1;
+
+    it = index_list.find( lookup );
+    if ( it != index_list.end() ) {
+        index = it->GetOrderedIndex();
+    }
+
+    return index;
+}
\ No newline at end of file