diff --git a/src/Airports/GenAirports850/linearfeature.cxx b/src/Airports/GenAirports850/linearfeature.cxx
index b8b0f7e4..bf6ceb3e 100644
--- a/src/Airports/GenAirports850/linearfeature.cxx
+++ b/src/Airports/GenAirports850/linearfeature.cxx
@@ -515,17 +515,20 @@ int LinearFeature::Finish( bool closed, unsigned int idx )
                 // first point on the mark - offset heading is 90deg
                 cur_outer = OffsetPointFirst( points.GetNode(j), points.GetNode(j+1), offset-width/2.0f );
                 cur_inner = OffsetPointFirst( points.GetNode(j), points.GetNode(j+1), offset+width/2.0f );
+                // OffsetPointsFirst( points.GetNode(j), points.GetNode(j+1), offset, width, cur_inner, cur_outer );
             }
             else if (j == marks[i]->end_idx)
             {
                 // last point on the mark - offset heading is 90deg
                 cur_outer = OffsetPointLast( points.GetNode(j-1), points.GetNode(j), offset-width/2.0f );
                 cur_inner = OffsetPointLast( points.GetNode(j-1), points.GetNode(j), offset+width/2.0f );
+                // OffsetPointsLast( points.GetNode(j-1), points.GetNode(j), offset, width, cur_inner, cur_outer );
             }
             else
             {
                 cur_outer = OffsetPointMiddle( points.GetNode(j-1), points.GetNode(j), points.GetNode(j+1), offset-width/2.0f );
                 cur_inner = OffsetPointMiddle( points.GetNode(j-1), points.GetNode(j), points.GetNode(j+1), offset+width/2.0f );
+                // OffsetPointsMiddle( points.GetNode(j-1), points.GetNode(j), points.GetNode(j+1), offset, width, cur_inner, cur_outer );
             }
 
             if ( markStarted )
diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt
index 7afacdd3..2e462776 100644
--- a/src/BuildTiles/Main/CMakeLists.txt
+++ b/src/BuildTiles/Main/CMakeLists.txt
@@ -37,19 +37,5 @@ target_link_libraries(tg-construct
 
 install(TARGETS tg-construct RUNTIME DESTINATION bin)
 
-
-
-add_executable(cliptst
-    cliptst.cxx)
-
-target_link_libraries(cliptst
-    terragear
-    ${GDAL_LIBRARY}
-    ${SIMGEAR_CORE_LIBRARIES}
-    ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
-)
-
-install(TARGETS cliptst  RUNTIME DESTINATION bin)
-
 INSTALL(FILES usgsmap.txt DESTINATION ${PKGDATADIR} )
 INSTALL(FILES default_priorities.txt DESTINATION ${PKGDATADIR} )
diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx
index 469f18a3..d72afc9f 100644
--- a/src/BuildTiles/Main/main.cxx
+++ b/src/BuildTiles/Main/main.cxx
@@ -166,8 +166,9 @@ int main(int argc, char **argv) {
         exit( -1 );
     }
 
-    // three identical work queues
-    SGLockedQueue<SGBucket> wq[3];
+    // tile work queue
+    std::vector<SGBucket> bucketList;
+    SGLockedQueue<SGBucket> wq;
 
     // First generate the workqueue of buckets to construct
     if (tile_id == -1) {
@@ -178,39 +179,27 @@ int main(int argc, char **argv) {
         SGBucket b_max( max );
 
         if ( b_min == b_max ) {
-            for (unsigned int q=0; q<3; q++) {
-                wq[q].push( b_min );
-            }
+			bucketList.push_back( b_min );
         } else {
-            int dx, dy;
-            int i, j;
-
-            sgBucketDiff(b_min, b_max, &dx, &dy);
-
             SG_LOG(SG_GENERAL, SG_ALERT, "  construction area spans tile boundaries");
-            SG_LOG(SG_GENERAL, SG_ALERT, "  dx = " << dx << "  dy = " << dy);
-
-            for ( j = 0; j <= dy; j++ ) {
-                for ( i = 0; i <= dx; i++ ) {
-                    for (unsigned int q=0; q<3; q++) {
-                        wq[q].push( sgBucketOffset(min.getLongitudeDeg(), min.getLatitudeDeg(), i, j) );
-                    }
-                }
-            }
+            sgGetBuckets( min, max, bucketList );
         }
     } else {
         // construct the specified tile
         SG_LOG(SG_GENERAL, SG_ALERT, "Building tile " << tile_id);
-        for (unsigned int q=0; q<3; q++) {
-            wq[q].push( SGBucket( tile_id ) );
-        }
+        bucketList.push_back( SGBucket( tile_id ) );
     }
 
+	/* fill the workqueue */
+	for (unsigned int i=0; i<bucketList.size(); i++) {
+	    wq.push( bucketList[i] );
+	}
+
     // now create the worker threads for stage 1
     std::vector<TGConstruct *> constructs;
 
     for (int i=0; i<num_threads; i++) {
-        TGConstruct* construct = new TGConstruct( areas, 1, wq[0] );
+        TGConstruct* construct = new TGConstruct( areas, 1, wq );
         //construct->set_cover( cover );
         construct->set_paths( work_dir, share_dir, output_dir, load_dirs );
         construct->set_options( ignoreLandmass, nudge );
@@ -223,7 +212,7 @@ int main(int argc, char **argv) {
         constructs[i]->start();
     }
     // wait for workqueue to empty
-    while( wq[0].size() ) {
+    while( wq.size() ) {
         tgSleep( 5 );
     }
     // wait for all threads to complete
@@ -237,8 +226,13 @@ int main(int argc, char **argv) {
     }
     constructs.clear();
 
+	/* fill the workqueue */
+	for (unsigned int i=0; i<bucketList.size(); i++) {
+	    wq.push( bucketList[i] );
+	}
+
     for (int i=0; i<num_threads; i++) {
-        TGConstruct* construct = new TGConstruct( areas, 2, wq[1] );
+        TGConstruct* construct = new TGConstruct( areas, 2, wq );
         //construct->set_cover( cover );
         construct->set_paths( work_dir, share_dir, output_dir, load_dirs );
         construct->set_options( ignoreLandmass, nudge );
@@ -251,7 +245,7 @@ int main(int argc, char **argv) {
         constructs[i]->start();
     }
     // wait for workqueue to empty
-    while( wq[1].size() ) {
+    while( wq.size() ) {
         tgSleep( 5 );
     }
     // wait for all threads to complete
@@ -264,8 +258,13 @@ int main(int argc, char **argv) {
     }
     constructs.clear();
 
+	/* fill the workqueue */
+	for (unsigned int i=0; i<bucketList.size(); i++) {
+	    wq.push( bucketList[i] );
+	}
+
     for (int i=0; i<num_threads; i++) {
-        TGConstruct* construct = new TGConstruct( areas, 3, wq[2] );
+        TGConstruct* construct = new TGConstruct( areas, 3, wq );
         //construct->set_cover( cover );
         construct->set_paths( work_dir, share_dir, output_dir, load_dirs );
         construct->set_options( ignoreLandmass, nudge );
@@ -278,7 +277,7 @@ int main(int argc, char **argv) {
         constructs[i]->start();
     }
     // wait for workqueue to empty
-    while( wq[2].size() ) {
+    while( wq.size() ) {
         tgSleep( 5 );
     }
     // wait for all threads to complete
diff --git a/src/BuildTiles/Main/tgconstruct_clip.cxx b/src/BuildTiles/Main/tgconstruct_clip.cxx
index 9d1753a8..6fc15d92 100644
--- a/src/BuildTiles/Main/tgconstruct_clip.cxx
+++ b/src/BuildTiles/Main/tgconstruct_clip.cxx
@@ -198,7 +198,10 @@ bool TGConstruct::ClipLandclassPolys( void ) {
 
 #if FIND_SLIVERS
     // Now, merge any slivers with clipped polys
-    merge_slivers(polys_clipped, slivers);
+    // merge_slivers(polys_clipped, slivers);
+    for ( unsigned int i = 0; i < area_defs.size(); i++ ) {
+        tgPolygon::MergeSlivers( polys_clipped.get_polys(i), slivers );
+    }
 #endif
 
     slivers.clear();
@@ -224,7 +227,9 @@ bool TGConstruct::ClipLandclassPolys( void ) {
                 }
             }
 
-            merge_slivers(polys_clipped, slivers);
+            for ( unsigned int i = 0; i < area_defs.size(); i++ ) {
+                tgPolygon::MergeSlivers( polys_clipped.get_polys(i), slivers );
+            }
         }
 #endif
 
@@ -253,4 +258,4 @@ bool TGConstruct::ClipLandclassPolys( void ) {
     }
 
     return true;
-}
\ No newline at end of file
+}
diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx
index cd75233a..89e180f7 100644
--- a/src/BuildTiles/Main/tgconstruct_poly.cxx
+++ b/src/BuildTiles/Main/tgconstruct_poly.cxx
@@ -29,6 +29,8 @@
 #include <simgear/misc/sg_dir.hxx>
 #include <simgear/debug/logstream.hxx>
 
+#include <terragear/tg_shapefile.hxx>
+
 #include "tgconstruct.hxx"
 
 using std::string;
@@ -103,6 +105,13 @@ int TGConstruct::LoadLandclassPolys( void ) {
                             }
                         }
                     }
+
+                    if (IsDebugShape( poly.GetId() )) {
+                        char layer[32];
+                        sprintf(layer, "loaded_%d", poly.GetId() );
+
+                        tgShapefile::FromPolygon( poly, ds_name, layer, material.c_str() );
+                    }
                 }
 
                 gzclose( fp );
diff --git a/src/BuildTiles/Main/tgconstruct_tesselate.cxx b/src/BuildTiles/Main/tgconstruct_tesselate.cxx
index 4fb0a9b5..b3006eaa 100644
--- a/src/BuildTiles/Main/tgconstruct_tesselate.cxx
+++ b/src/BuildTiles/Main/tgconstruct_tesselate.cxx
@@ -42,7 +42,9 @@ void TGConstruct::TesselatePolys( void )
             tgPolygon poly = polys_clipped.get_poly(area, p );
 
             if ( IsDebugShape( poly.GetId() ) ) {
-                tgShapefile::FromPolygon( poly, ds_name, "preteselate", "" );
+                char pn[32];
+                sprintf(pn, "%d", poly.GetId() );
+                tgShapefile::FromPolygon( poly, ds_name, "preteselate", pn );
             }
 
             tgRectangle rect = poly.GetBoundingBox();
diff --git a/src/Lib/terragear/clipper.cpp b/src/Lib/terragear/clipper.cpp
index 15e8b783..e3e61786 100644
--- a/src/Lib/terragear/clipper.cpp
+++ b/src/Lib/terragear/clipper.cpp
@@ -1,10 +1,10 @@
 /*******************************************************************************
 *                                                                              *
 * Author    :  Angus Johnson                                                   *
-* Version   :  4.9.1                                                           *
-* Date      :  9 October 2012                                                  *
+* Version   :  5.1.2                                                           *
+* Date      :  25 February 2013                                                *
 * Website   :  http://www.angusj.com                                           *
-* Copyright :  Angus Johnson 2010-2012                                         *
+* Copyright :  Angus Johnson 2010-2013                                         *
 *                                                                              *
 * License:                                                                     *
 * Use, modification & distribution is subject to Boost Software License Ver 1. *
@@ -26,7 +26,7 @@
 * Paper no. DETC2005-85513 pp. 565-575                                         *
 * ASME 2005 International Design Engineering Technical Conferences             *
 * and Computers and Information in Engineering Conference (IDETC/CIE2005)      *
-* September 24�28, 2005 , Long Beach, California, USA                          *
+* September 24-28, 2005 , Long Beach, California, USA                          *
 * http://www.me.berkeley.edu/~mcmains/pubs/DAC05OffsetPolygon.pdf              *
 *                                                                              *
 *******************************************************************************/
@@ -63,7 +63,89 @@ inline long64 Abs(long64 val)
 {
   return val < 0 ? -val : val;
 }
+
 //------------------------------------------------------------------------------
+// PolyTree methods ...
+//------------------------------------------------------------------------------
+
+void PolyTree::Clear()
+{
+    for (PolyNodes::size_type i = 0; i < AllNodes.size(); ++i)
+      delete AllNodes[i];
+    AllNodes.resize(0); 
+    Childs.resize(0);
+}
+//------------------------------------------------------------------------------
+
+PolyNode* PolyTree::GetFirst() const
+{
+  if (Childs.size() > 0)
+      return Childs[0];
+  else
+      return 0;
+}
+//------------------------------------------------------------------------------
+
+int PolyTree::Total() const
+{
+  return AllNodes.size();
+}
+
+//------------------------------------------------------------------------------
+// PolyNode methods ...
+//------------------------------------------------------------------------------
+
+PolyNode::PolyNode(): Childs(), Parent(0), Index(0)
+{
+}
+//------------------------------------------------------------------------------
+
+int PolyNode::ChildCount() const
+{
+  return Childs.size();
+}
+//------------------------------------------------------------------------------
+
+void PolyNode::AddChild(PolyNode& child)
+{
+  unsigned cnt = Childs.size();
+  Childs.push_back(&child);
+  child.Parent = this;
+  child.Index = cnt;
+}
+//------------------------------------------------------------------------------
+
+PolyNode* PolyNode::GetNext() const
+{ 
+  if (Childs.size() > 0) 
+      return Childs[0]; 
+  else
+      return GetNextSiblingUp();    
+}  
+//------------------------------------------------------------------------------
+
+PolyNode* PolyNode::GetNextSiblingUp() const
+{ 
+  if (!Parent) //protects against PolyTree.GetNextSiblingUp()
+      return 0;
+  else if (Index == Parent->Childs.size() - 1)
+      return Parent->GetNextSiblingUp();
+  else
+      return Parent->Childs[Index + 1];
+}  
+//------------------------------------------------------------------------------
+
+bool PolyNode::IsHole() const
+{ 
+  bool result = true;
+  PolyNode* node = Parent;
+  while (node)
+  {
+      result = !result;
+      node = node->Parent;
+  }
+  return result;
+}  
 
 //------------------------------------------------------------------------------
 // Int128 class (enables safe math on signed 64bit integers)
@@ -77,18 +159,24 @@ class Int128
 {
   public:
 
+    ulong64 lo;
+    long64 hi;
+
     Int128(long64 _lo = 0)
     {
-      lo = _lo;
-      if (lo < 0) hi = -1; else hi = 0;
+      lo = (ulong64)_lo;   
+      if (_lo < 0)  hi = -1; else hi = 0; 
     }
 
-    Int128(const Int128 &val): hi(val.hi), lo(val.lo){}
 
+    Int128(const Int128 &val): lo(val.lo), hi(val.hi){}
+
+    Int128(const long64& _hi, const ulong64& _lo): lo(_lo), hi(_hi){}
+    
     long64 operator = (const long64 &val)
     {
-      lo = val;
-      if (lo < 0) hi = -1; else hi = 0;
+      lo = (ulong64)val;
+      if (val < 0) hi = -1; else hi = 0;
       return val;
     }
 
@@ -124,7 +212,7 @@ class Int128
     {
       hi += rhs.hi;
       lo += rhs.lo;
-      if (ulong64(lo) < ulong64(rhs.lo)) hi++;
+      if (lo < rhs.lo) hi++;
       return *this;
     }
 
@@ -137,25 +225,10 @@ class Int128
 
     Int128& operator -= (const Int128 &rhs)
     {
-      Int128 tmp(rhs);
-      Negate(tmp);
-      *this += tmp;
+      *this += -rhs;
       return *this;
     }
 
-    //Int128 operator -() const
-    //{
-    //  Int128 result(*this);
-    //  if (result.lo == 0) {
-    //    if (result.hi != 0) result.hi = -1;
-    //  }
-    //  else {
-    //    result.lo = -result.lo;
-    //    result.hi = ~result.hi;
-    //  }
-    //  return result;
-    //}
-
     Int128 operator - (const Int128 &rhs) const
     {
       Int128 result(*this);
@@ -163,150 +236,115 @@ class Int128
       return result;
     }
 
-    Int128 operator * (const Int128 &rhs) const
+    Int128 operator-() const //unary negation
     {
-      if ( !(hi == 0 || hi == -1) || !(rhs.hi == 0 || rhs.hi == -1))
-        throw "Int128 operator*: overflow error";
-      bool negate = (hi < 0) != (rhs.hi < 0);
-
-      Int128 tmp(*this);
-      if (tmp.hi < 0) Negate(tmp);
-      ulong64 int1Hi = ulong64(tmp.lo) >> 32;
-      ulong64 int1Lo = ulong64(tmp.lo & 0xFFFFFFFF);
-
-      tmp = rhs;
-      if (tmp.hi < 0) Negate(tmp);
-      ulong64 int2Hi = ulong64(tmp.lo) >> 32;
-      ulong64 int2Lo = ulong64(tmp.lo & 0xFFFFFFFF);
-
-      //nb: see comments in clipper.pas
-      ulong64 a = int1Hi * int2Hi;
-      ulong64 b = int1Lo * int2Lo;
-      ulong64 c = int1Hi * int2Lo + int1Lo * int2Hi;
-
-      tmp.hi = long64(a + (c >> 32));
-      tmp.lo = long64(c << 32);
-      tmp.lo += long64(b);
-      if (ulong64(tmp.lo) < b) tmp.hi++;
-      if (negate) Negate(tmp);
-      return tmp;
+      if (lo == 0)
+        return Int128(-hi,0);
+      else 
+        return Int128(~hi,~lo +1);
     }
 
     Int128 operator/ (const Int128 &rhs) const
     {
       if (rhs.lo == 0 && rhs.hi == 0)
         throw "Int128 operator/: divide by zero";
-      bool negate = (rhs.hi < 0) != (hi < 0);
-      Int128 result(*this), denom(rhs);
-      if (result.hi < 0) Negate(result);
-      if (denom.hi < 0)  Negate(denom);
-      if (denom > result) return Int128(0); //result is only a fraction of 1
-      Negate(denom);
 
-      Int128 p(0);
-      for (int i = 0; i < 128; ++i)
+      bool negate = (rhs.hi < 0) != (hi < 0);
+      Int128 dividend = *this;
+      Int128 divisor = rhs;
+      if (dividend.hi < 0) dividend = -dividend;
+      if (divisor.hi < 0) divisor = -divisor;
+
+      if (divisor < dividend)
       {
-        p.hi = p.hi << 1;
-        if (p.lo < 0) p.hi++;
-        p.lo = long64(p.lo) << 1;
-        if (result.hi < 0) p.lo++;
-        result.hi = result.hi << 1;
-        if (result.lo < 0) result.hi++;
-        result.lo = long64(result.lo) << 1;
-        Int128 p2(p);
-        p += denom;
-        if (p.hi < 0) p = p2;
-        else result.lo++;
+          Int128 result = Int128(0);
+          Int128 cntr = Int128(1);
+          while (divisor.hi >= 0 && !(divisor > dividend))
+          {
+              divisor.hi <<= 1;
+              if ((long64)divisor.lo < 0) divisor.hi++;
+              divisor.lo <<= 1;
+
+              cntr.hi <<= 1;
+              if ((long64)cntr.lo < 0) cntr.hi++;
+              cntr.lo <<= 1;
+          }
+          divisor.lo >>= 1;
+          if ((divisor.hi & 1) == 1)
+              divisor.lo |= 0x8000000000000000LL; 
+          divisor.hi = (ulong64)divisor.hi >> 1;
+
+          cntr.lo >>= 1;
+          if ((cntr.hi & 1) == 1)
+              cntr.lo |= 0x8000000000000000LL; 
+          cntr.hi >>= 1;
+
+          while (cntr.hi != 0 || cntr.lo != 0)
+          {
+              if (!(dividend < divisor))
+              {
+                  dividend -= divisor;
+                  result.hi |= cntr.hi;
+                  result.lo |= cntr.lo;
+              }
+              divisor.lo >>= 1;
+              if ((divisor.hi & 1) == 1)
+                  divisor.lo |= 0x8000000000000000LL; 
+              divisor.hi >>= 1;
+
+              cntr.lo >>= 1;
+              if ((cntr.hi & 1) == 1)
+                  cntr.lo |= 0x8000000000000000LL; 
+              cntr.hi >>= 1;
+          }
+          if (negate) result = -result;
+          return result;
       }
-      if (negate) Negate(result);
-      return result;
+      else if (rhs.hi == this->hi && rhs.lo == this->lo)
+          return Int128(1);
+      else
+          return Int128(0);
     }
 
     double AsDouble() const
     {
       const double shift64 = 18446744073709551616.0; //2^64
-      const double bit64 = 9223372036854775808.0;
       if (hi < 0)
       {
-        Int128 tmp(*this);
-        Negate(tmp);
-        if (tmp.lo < 0)
-          return (double)tmp.lo - bit64 - tmp.hi * shift64;
-        else
-          return -(double)tmp.lo - tmp.hi * shift64;
+        if (lo == 0) return (double)hi * shift64;
+        else return -(double)(~lo + ~hi * shift64);
       }
-      else if (lo < 0)
-        return -(double)lo + bit64 + hi * shift64;
       else
-        return (double)lo + (double)hi * shift64;
+        return (double)(lo + hi * shift64);
     }
-
-    //for bug testing ...
-    //std::string AsString() const
-    //{
-    //  std::string result;
-    //  unsigned char r = 0;
-    //  Int128 tmp(0), val(*this);
-    //  if (hi < 0) Negate(val);
-    //  result.resize(50);
-    //  std::string::size_type i = result.size() -1;
-    //  while (val.hi != 0 || val.lo != 0)
-    //  {
-    //    Div10(val, tmp, r);
-    //    result[i--] = char('0' + r);
-    //    val = tmp;
-    //  }
-    //  if (hi < 0) result[i--] = '-';
-    //  result.erase(0,i+1);
-    //  if (result.size() == 0) result = "0";
-    //  return result;
-    //}
-
-private:
-    long64 hi;
-    long64 lo;
-
-    static void Negate(Int128 &val)
-    {
-      if (val.lo == 0) {
-        if (val.hi != 0) val.hi = -val.hi;;
-      }
-      else {
-        val.lo = -val.lo;
-        val.hi = ~val.hi;
-      }
-    }
-
-    //debugging only ...
-    //void Div10(const Int128 val, Int128& result, unsigned char & remainder) const
-    //{
-    //  remainder = 0;
-    //  result = 0;
-    //  for (int i = 63; i >= 0; --i)
-    //  {
-    //    if ((val.hi & ((long64)1 << i)) != 0)
-    //      remainder = char((remainder * 2) + 1); else
-    //      remainder *= char(2);
-    //    if (remainder >= 10)
-    //    {
-    //      result.hi += ((long64)1 << i);
-    //      remainder -= char(10);
-    //    }
-    //  }
-    //  for (int i = 63; i >= 0; --i)
-    //  {
-    //    if ((val.lo & ((long64)1 << i)) != 0)
-    //      remainder = char((remainder * 2) + 1); else
-    //      remainder *= char(2);
-    //    if (remainder >= 10)
-    //    {
-    //      result.lo += ((long64)1 << i);
-    //      remainder -= char(10);
-    //    }
-    //  }
-    //}
 };
 
+Int128 Int128Mul (long64 lhs, long64 rhs)
+{
+  bool negate = (lhs < 0) != (rhs < 0);
+
+  if (lhs < 0) lhs = -lhs;
+  ulong64 int1Hi = ulong64(lhs) >> 32;
+  ulong64 int1Lo = ulong64(lhs & 0xFFFFFFFF);
+
+  if (rhs < 0) rhs = -rhs;
+  ulong64 int2Hi = ulong64(rhs) >> 32;
+  ulong64 int2Lo = ulong64(rhs & 0xFFFFFFFF);
+
+  //nb: see comments in clipper.pas
+  ulong64 a = int1Hi * int2Hi;
+  ulong64 b = int1Lo * int2Lo;
+  ulong64 c = int1Hi * int2Lo + int1Lo * int2Hi;
+
+  Int128 tmp;
+  tmp.hi = long64(a + (c >> 32));
+  tmp.lo = long64(c << 32);
+  tmp.lo += long64(b);
+  if (tmp.lo < b) tmp.hi++;
+  if (negate) tmp = -tmp;
+  return tmp;
+}
+
 //------------------------------------------------------------------------------
 //------------------------------------------------------------------------------
 
@@ -326,39 +364,7 @@ bool FullRangeNeeded(const Polygon &pts)
   
 bool Orientation(const Polygon &poly)
 {
-  int highI = (int)poly.size() -1;
-  if (highI < 2) return false;
-
-  int j = 0, jplus, jminus;
-  for (int i = 0; i <= highI; ++i)
-  {
-    if (poly[i].Y < poly[j].Y) continue;
-    if ((poly[i].Y > poly[j].Y || poly[i].X < poly[j].X)) j = i;
-  };
-  if (j == highI) jplus = 0;
-  else jplus = j +1;
-  if (j == 0) jminus = highI;
-  else jminus = j -1;
-
-  IntPoint vec1, vec2;
-  //get cross product of vectors of the edges adjacent to highest point ...
-  vec1.X = poly[j].X - poly[jminus].X;
-  vec1.Y = poly[j].Y - poly[jminus].Y;
-  vec2.X = poly[jplus].X - poly[j].X;
-  vec2.Y = poly[jplus].Y - poly[j].Y;
-
-  if (Abs(vec1.X) > loRange || Abs(vec1.Y) > loRange ||
-    Abs(vec2.X) > loRange || Abs(vec2.Y) > loRange)
-  {
-    if (Abs(vec1.X) > hiRange || Abs(vec1.Y) > hiRange ||
-      Abs(vec2.X) > hiRange || Abs(vec2.Y) > hiRange)
-        throw "Coordinate exceeds range bounds.";
-    Int128 cross = Int128(vec1.X) * Int128(vec2.Y) -
-      Int128(vec2.X) * Int128(vec1.Y);
-    return cross >= 0;
-  }
-  else
-    return (vec1.X * vec2.Y - vec2.X * vec1.Y) >= 0;
+    return Area(poly) >= 0;
 }
 //------------------------------------------------------------------------------
 
@@ -368,44 +374,6 @@ inline bool PointsEqual( const IntPoint &pt1, const IntPoint &pt2)
 }
 //------------------------------------------------------------------------------
 
-bool Orientation(OutRec *outRec, bool UseFullInt64Range)
-{
-  //first make sure bottomPt is correctly assigned ...
-  OutPt *opBottom = outRec->pts, *op = outRec->pts->next;
-  while (op != outRec->pts)
-  {
-    if (op->pt.Y >= opBottom->pt.Y)
-    {
-      if (op->pt.Y > opBottom->pt.Y || op->pt.X < opBottom->pt.X)
-      opBottom = op;
-    }
-    op = op->next;
-  }
-  outRec->bottomPt = opBottom;
-  opBottom->idx = outRec->idx;
-
-  op = opBottom;
-  //find vertices either side of bottomPt (skipping duplicate points) ....
-  OutPt *opPrev = op->prev;
-  OutPt *opNext = op->next;
-  while (op != opPrev && PointsEqual(op->pt, opPrev->pt))
-    opPrev = opPrev->prev;
-  while (op != opNext && PointsEqual(op->pt, opNext->pt))
-    opNext = opNext->next;
-
-  IntPoint ip1, ip2;
-  ip1.X = op->pt.X - opPrev->pt.X;
-  ip1.Y = op->pt.Y - opPrev->pt.Y;
-  ip2.X = opNext->pt.X - op->pt.X;
-  ip2.Y = opNext->pt.Y - op->pt.Y;
-
-  if (UseFullInt64Range)
-    return Int128(ip1.X) * Int128(ip2.Y) - Int128(ip2.X) * Int128(ip1.Y) >= 0;
-  else
-    return (ip1.X * ip2.Y - ip2.X * ip1.Y) >= 0;
-}
-//------------------------------------------------------------------------------
-
 double Area(const Polygon &poly)
 {
   int highI = (int)poly.size() -1;
@@ -413,20 +381,18 @@ double Area(const Polygon &poly)
 
   if (FullRangeNeeded(poly)) {
     Int128 a;
-    a = (Int128(poly[highI].X) * Int128(poly[0].Y)) -
-      Int128(poly[0].X) * Int128(poly[highI].Y);
-    for (int i = 0; i < highI; ++i)
-      a += Int128(poly[i].X) * Int128(poly[i+1].Y) -
-        Int128(poly[i+1].X) * Int128(poly[i].Y);
+    a = Int128Mul(poly[highI].X + poly[0].X, poly[0].Y - poly[highI].Y);
+    for (int i = 1; i <= highI; ++i)
+      a += Int128Mul(poly[i - 1].X + poly[i].X, poly[i].Y - poly[i -1].Y);
     return a.AsDouble() / 2;
   }
   else
   {
     double a;
-    a = (double)poly[highI].X * poly[0].Y - (double)poly[0].X * poly[highI].Y;
-    for (int i = 0; i < highI; ++i)
-      a += (double)poly[i].X * poly[i+1].Y - (double)poly[i+1].X * poly[i].Y;
-    return a/2;
+    a = ((double)poly[highI].X + poly[0].X) * ((double)poly[0].Y - poly[highI].Y);
+    for (int i = 1; i <= highI; ++i)
+      a += ((double)poly[i - 1].X + poly[i].X) * ((double)poly[i].Y - poly[i - 1].Y);
+    return a / 2;
   }
 }
 //------------------------------------------------------------------------------
@@ -434,11 +400,11 @@ double Area(const Polygon &poly)
 double Area(const OutRec &outRec, bool UseFullInt64Range)
 {
   OutPt *op = outRec.pts;
+  if (!op) return 0;
   if (UseFullInt64Range) {
     Int128 a(0);
     do {
-      a += (Int128(op->prev->pt.X) * Int128(op->pt.Y)) -
-        Int128(op->pt.X) * Int128(op->prev->pt.Y);
+      a += Int128Mul(op->pt.X + op->prev->pt.X, op->prev->pt.Y - op->pt.Y);
       op = op->next;
     } while (op != outRec.pts);
     return a.AsDouble() / 2;
@@ -447,10 +413,10 @@ double Area(const OutRec &outRec, bool UseFullInt64Range)
   {
     double a = 0;
     do {
-      a += (op->prev->pt.X * op->pt.Y) - (op->pt.X * op->prev->pt.Y);
+      a = a + (op->pt.X + op->prev->pt.X) * (op->prev->pt.Y - op->pt.Y);
       op = op->next;
     } while (op != outRec.pts);
-    return a/2;
+    return a / 2;
   }
 }
 //------------------------------------------------------------------------------
@@ -477,8 +443,9 @@ bool PointInPolygon(const IntPoint &pt, OutPt *pp, bool UseFullInt64Range)
     {
       if ((((pp2->pt.Y <= pt.Y) && (pt.Y < pp2->prev->pt.Y)) ||
           ((pp2->prev->pt.Y <= pt.Y) && (pt.Y < pp2->pt.Y))) &&
-          Int128(pt.X - pp2->pt.X) < (Int128(pp2->prev->pt.X - pp2->pt.X) *
-          Int128(pt.Y - pp2->pt.Y)) / Int128(pp2->prev->pt.Y - pp2->pt.Y))
+          Int128(pt.X - pp2->pt.X) < 
+          Int128Mul(pp2->prev->pt.X - pp2->pt.X, pt.Y - pp2->pt.Y) / 
+          Int128(pp2->prev->pt.Y - pp2->pt.Y))
             result = !result;
       pp2 = pp2->next;
     }
@@ -503,8 +470,7 @@ bool PointInPolygon(const IntPoint &pt, OutPt *pp, bool UseFullInt64Range)
 bool SlopesEqual(TEdge &e1, TEdge &e2, bool UseFullInt64Range)
 {
   if (UseFullInt64Range)
-    return Int128(e1.deltaY) * Int128(e2.deltaX) ==
-      Int128(e1.deltaX) * Int128(e2.deltaY);
+    return Int128Mul(e1.deltaY, e2.deltaX) == Int128Mul(e1.deltaX, e2.deltaY);
   else return e1.deltaY * e2.deltaX == e1.deltaX * e2.deltaY;
 }
 //------------------------------------------------------------------------------
@@ -513,8 +479,7 @@ bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
   const IntPoint pt3, bool UseFullInt64Range)
 {
   if (UseFullInt64Range)
-    return Int128(pt1.Y-pt2.Y) * Int128(pt2.X-pt3.X) ==
-      Int128(pt1.X-pt2.X) * Int128(pt2.Y-pt3.Y);
+    return Int128Mul(pt1.Y-pt2.Y, pt2.X-pt3.X) == Int128Mul(pt1.X-pt2.X, pt2.Y-pt3.Y);
   else return (pt1.Y-pt2.Y)*(pt2.X-pt3.X) == (pt1.X-pt2.X)*(pt2.Y-pt3.Y);
 }
 //------------------------------------------------------------------------------
@@ -523,8 +488,7 @@ bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
   const IntPoint pt3, const IntPoint pt4, bool UseFullInt64Range)
 {
   if (UseFullInt64Range)
-    return Int128(pt1.Y-pt2.Y) * Int128(pt3.X-pt4.X) ==
-      Int128(pt1.X-pt2.X) * Int128(pt3.Y-pt4.Y);
+    return Int128Mul(pt1.Y-pt2.Y, pt3.X-pt4.X) == Int128Mul(pt1.X-pt2.X, pt3.Y-pt4.Y);
   else return (pt1.Y-pt2.Y)*(pt3.X-pt4.X) == (pt1.X-pt2.X)*(pt3.Y-pt4.Y);
 }
 //------------------------------------------------------------------------------
@@ -532,7 +496,7 @@ bool SlopesEqual(const IntPoint pt1, const IntPoint pt2,
 double GetDx(const IntPoint pt1, const IntPoint pt2)
 {
   return (pt1.Y == pt2.Y) ?
-    HORIZONTAL : (double)(pt2.X - pt1.X) / (double)(pt2.Y - pt1.Y);
+    HORIZONTAL : (double)(pt2.X - pt1.X) / (pt2.Y - pt1.Y);
 }
 //---------------------------------------------------------------------------
 
@@ -542,7 +506,7 @@ void SetDx(TEdge &e)
   e.deltaY = (e.ytop - e.ybot);
 
   if (e.deltaY == 0) e.dx = HORIZONTAL;
-  else e.dx = (double)(e.deltaX) / (double)(e.deltaY);
+  else e.dx = (double)(e.deltaX) / e.deltaY;
 }
 //---------------------------------------------------------------------------
 
@@ -575,20 +539,6 @@ long64 TopX(TEdge &edge, const long64 currentY)
 }
 //------------------------------------------------------------------------------
 
-long64 TopX(const IntPoint pt1, const IntPoint pt2, const long64 currentY)
-{
-  //preconditions: pt1.Y <> pt2.Y and pt1.Y > pt2.Y
-  if (currentY >= pt1.Y) return pt1.X;
-  else if (currentY == pt2.Y) return pt2.X;
-  else if (pt1.X == pt2.X) return pt1.X;
-  else
-  {
-    double q = (double)(pt1.X-pt2.X)/(double)(pt1.Y-pt2.Y);
-    return Round(pt1.X + (currentY - pt1.Y) *q);
-  }
-}
-//------------------------------------------------------------------------------
-
 bool IntersectPoint(TEdge &edge1, TEdge &edge2,
   IntPoint &ip, bool UseFullInt64Range)
 {
@@ -602,8 +552,8 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
       ip.Y = edge2.ybot;
     } else
     {
-      b2 = edge2.ybot - (edge2.xbot/edge2.dx);
-      ip.Y = Round(ip.X/edge2.dx + b2);
+      b2 = edge2.ybot - (edge2.xbot / edge2.dx);
+      ip.Y = Round(ip.X / edge2.dx + b2);
     }
   }
   else if (NEAR_ZERO(edge2.dx))
@@ -614,16 +564,19 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
       ip.Y = edge1.ybot;
     } else
     {
-      b1 = edge1.ybot - (edge1.xbot/edge1.dx);
-      ip.Y = Round(ip.X/edge1.dx + b1);
+      b1 = edge1.ybot - (edge1.xbot / edge1.dx);
+      ip.Y = Round(ip.X / edge1.dx + b1);
     }
-  } else
+  } else 
   {
     b1 = edge1.xbot - edge1.ybot * edge1.dx;
     b2 = edge2.xbot - edge2.ybot * edge2.dx;
-    b2 = (b2-b1)/(edge1.dx - edge2.dx);
-    ip.Y = Round(b2);
-    ip.X = Round(edge1.dx * b2 + b1);
+    double q = (b2-b1) / (edge1.dx - edge2.dx);
+    ip.Y = Round(q);
+    if (std::fabs(edge1.dx) < std::fabs(edge2.dx))
+      ip.X = Round(edge1.dx * q + b1);
+    else 
+      ip.X = Round(edge2.dx * q + b2);
   }
 
   if (ip.Y < edge1.ytop || ip.Y < edge2.ytop) 
@@ -645,16 +598,17 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
 }
 //------------------------------------------------------------------------------
 
-void ReversePolyPtLinks(OutPt &pp)
+void ReversePolyPtLinks(OutPt *pp)
 {
+  if (!pp) return;
   OutPt *pp1, *pp2;
-  pp1 = &pp;
+  pp1 = pp;
   do {
   pp2 = pp1->next;
   pp1->next = pp1->prev;
   pp1->prev = pp2;
   pp1 = pp2;
-  } while( pp1 != &pp );
+  } while( pp1 != pp );
 }
 //------------------------------------------------------------------------------
 
@@ -724,7 +678,7 @@ bool GetOverlapSegment(IntPoint pt1a, IntPoint pt1b, IntPoint pt2a,
   IntPoint pt2b, IntPoint &pt1, IntPoint &pt2)
 {
   //precondition: segments are colinear.
-  if ( pt1a.Y == pt1b.Y || Abs((pt1a.X - pt1b.X)/(pt1a.Y - pt1b.Y)) > 1 )
+  if (Abs(pt1a.X - pt1b.X) > Abs(pt1a.Y - pt1b.Y))
   {
     if (pt1a.X > pt1b.X) SwapPoints(pt1a, pt1b);
     if (pt2a.X > pt2b.X) SwapPoints(pt2a, pt2b);
@@ -799,7 +753,8 @@ OutPt* GetBottomPt(OutPt *pp)
 }
 //------------------------------------------------------------------------------
 
-bool FindSegment(OutPt* &pp, IntPoint &pt1, IntPoint &pt2)
+bool FindSegment(OutPt* &pp, bool UseFullInt64Range, 
+  IntPoint &pt1, IntPoint &pt2)
 {
   //outPt1 & outPt2 => the overlap segment (if the function returns true)
   if (!pp) return false;
@@ -807,8 +762,8 @@ bool FindSegment(OutPt* &pp, IntPoint &pt1, IntPoint &pt2)
   IntPoint pt1a = pt1, pt2a = pt2;
   do
   {
-    if (SlopesEqual(pt1a, pt2a, pp->pt, pp->prev->pt, true) &&
-      SlopesEqual(pt1a, pt2a, pp->pt, true) &&
+    if (SlopesEqual(pt1a, pt2a, pp->pt, pp->prev->pt, UseFullInt64Range) &&
+      SlopesEqual(pt1a, pt2a, pp->pt, UseFullInt64Range) &&
       GetOverlapSegment(pt1a, pt2a, pp->pt, pp->prev->pt, pt1, pt2))
         return true;
     pp = pp->next;
@@ -870,6 +825,7 @@ bool ClipperBase::AddPolygon( const Polygon &pg, PolyType polyType)
 {
   int len = (int)pg.size();
   if (len < 3) return false;
+
   Polygon p(len);
   p[0] = pg[0];
   int j = 0;
@@ -1209,89 +1165,46 @@ bool Clipper::Execute(ClipType clipType, Polygons &solution,
   m_SubjFillType = subjFillType;
   m_ClipFillType = clipFillType;
   m_ClipType = clipType;
-  bool succeeded = ExecuteInternal(false);
+  m_UsingPolyTree = false;
+  bool succeeded = ExecuteInternal();
   if (succeeded) BuildResult(solution);
   m_ExecuteLocked = false;
   return succeeded;
 }
 //------------------------------------------------------------------------------
 
-bool Clipper::Execute(ClipType clipType, ExPolygons &solution,
+bool Clipper::Execute(ClipType clipType, PolyTree& polytree,
     PolyFillType subjFillType, PolyFillType clipFillType)
 {
   if( m_ExecuteLocked ) return false;
   m_ExecuteLocked = true;
-  solution.resize(0);
   m_SubjFillType = subjFillType;
   m_ClipFillType = clipFillType;
   m_ClipType = clipType;
-  bool succeeded = ExecuteInternal(true);
-  if (succeeded) BuildResultEx(solution);
+  m_UsingPolyTree = true;
+  bool succeeded = ExecuteInternal();
+  if (succeeded) BuildResult2(polytree);
   m_ExecuteLocked = false;
   return succeeded;
 }
 //------------------------------------------------------------------------------
 
-bool PolySort(OutRec *or1, OutRec *or2)
+void Clipper::FixHoleLinkage(OutRec &outRec)
 {
-  if (or1 == or2) return false;
-  if (!or1->pts || !or2->pts)
-  {
-    if (or1->pts != or2->pts)
-    {
-      return or1->pts ? true : false;
-    }
-    else return false;
-  }
-  int i1, i2;
-  if (or1->isHole)
-    i1 = or1->FirstLeft->idx; else
-    i1 = or1->idx;
-  if (or2->isHole)
-    i2 = or2->FirstLeft->idx; else
-    i2 = or2->idx;
-  int result = i1 - i2;
-  if (result == 0 && (or1->isHole != or2->isHole))
-  {
-    return or1->isHole ? false : true;
-  }
-  else return result < 0;
+  //skip OutRecs that (a) contain outermost polygons or
+  //(b) already have the correct owner/child linkage ...
+  if (!outRec.FirstLeft ||                
+      (outRec.isHole != outRec.FirstLeft->isHole &&
+      outRec.FirstLeft->pts)) return;
+
+  OutRec* orfl = outRec.FirstLeft;
+  while (orfl && ((orfl->isHole == outRec.isHole) || !orfl->pts))
+      orfl = orfl->FirstLeft;
+  outRec.FirstLeft = orfl;
 }
 //------------------------------------------------------------------------------
 
-OutRec* FindAppendLinkEnd(OutRec *outRec)
-{
-  while (outRec->AppendLink) outRec = outRec->AppendLink;
-  return outRec;
-}
-//------------------------------------------------------------------------------
-
-void Clipper::FixHoleLinkage(OutRec *outRec)
-{
-  OutRec *tmp;
-  if (outRec->bottomPt)
-    tmp = m_PolyOuts[outRec->bottomPt->idx]->FirstLeft;
-  else
-    tmp = outRec->FirstLeft;
-  if (outRec == tmp) throw clipperException("HoleLinkage error");
-
-  if (tmp)
-  {
-    if (tmp->AppendLink) tmp = FindAppendLinkEnd(tmp);
-    if (tmp == outRec) tmp = 0;
-    else if (tmp->isHole)
-    {
-      FixHoleLinkage(tmp);
-      tmp = tmp->FirstLeft;
-    }
-  }
-  outRec->FirstLeft = tmp;
-  if (!tmp) outRec->isHole = false;
-  outRec->AppendLink = 0;
-}
-//------------------------------------------------------------------------------
-
-bool Clipper::ExecuteInternal(bool fixHoleLinkages)
+bool Clipper::ExecuteInternal()
 {
   bool succeeded;
   try {
@@ -1322,20 +1235,12 @@ bool Clipper::ExecuteInternal(bool fixHoleLinkages)
       if (!outRec->pts) continue;
       FixupOutPolygon(*outRec);
       if (!outRec->pts) continue;
-      if (outRec->isHole && fixHoleLinkages) FixHoleLinkage(outRec);
 
-      if (outRec->bottomPt == outRec->bottomFlag &&
-        (Orientation(outRec, m_UseFullRange) != (Area(*outRec, m_UseFullRange) > 0)))
-          DisposeBottomPt(*outRec);
-
-      if (outRec->isHole ==
-        (m_ReverseOutput ^ Orientation(outRec, m_UseFullRange)))
-          ReversePolyPtLinks(*outRec->pts);
+      if ((outRec->isHole ^ m_ReverseOutput) == (Area(*outRec, m_UseFullRange) > 0))
+        ReversePolyPtLinks(outRec->pts);
     }
 
-    JoinCommonEdges(fixHoleLinkages);
-    if (fixHoleLinkages)
-      std::sort(m_PolyOuts.begin(), m_PolyOuts.end(), PolySort);
+    if (m_Joins.size() > 0) JoinCommonEdges();
   }
 
   ClearJoins();
@@ -1624,14 +1529,10 @@ void Clipper::CopyAELToSEL()
 {
   TEdge* e = m_ActiveEdges;
   m_SortedEdges = e;
-  if (!m_ActiveEdges) return;
-  m_SortedEdges->prevInSEL = 0;
-  e = e->nextInAEL;
   while ( e )
   {
     e->prevInSEL = e->prevInAEL;
-    e->prevInSEL->nextInSEL = e;
-    e->nextInSEL = 0;
+    e->nextInSEL = e->nextInAEL;
     e = e->nextInAEL;
   }
 }
@@ -1679,7 +1580,7 @@ void Clipper::ClearHorzJoins()
 }
 //------------------------------------------------------------------------------
 
-void Clipper::InsertLocalMinimaIntoAEL( const long64 botY)
+void Clipper::InsertLocalMinimaIntoAEL(const long64 botY)
 {
   while(  m_CurrentLM  && ( m_CurrentLM->Y == botY ) )
   {
@@ -1716,21 +1617,18 @@ void Clipper::InsertLocalMinimaIntoAEL( const long64 botY)
       AddLocalMinPoly( lb, rb, IntPoint(lb->xcurr, m_CurrentLM->Y) );
 
     //if any output polygons share an edge, they'll need joining later ...
-    if (rb->outIdx >= 0)
+    if (rb->outIdx >= 0 && NEAR_EQUAL(rb->dx, HORIZONTAL))
     {
-      if (NEAR_EQUAL(rb->dx, HORIZONTAL))
+      for (HorzJoinList::size_type i = 0; i < m_HorizJoins.size(); ++i)
       {
-        for (HorzJoinList::size_type i = 0; i < m_HorizJoins.size(); ++i)
-        {
-          IntPoint pt, pt2; //returned by GetOverlapSegment() but unused here.
-          HorzJoinRec* hj = m_HorizJoins[i];
-          //if horizontals rb and hj.edge overlap, flag for joining later ...
-          if (GetOverlapSegment(IntPoint(hj->edge->xbot, hj->edge->ybot),
-            IntPoint(hj->edge->xtop, hj->edge->ytop),
-            IntPoint(rb->xbot, rb->ybot),
-            IntPoint(rb->xtop, rb->ytop), pt, pt2))
-              AddJoin(hj->edge, rb, hj->savedIdx);
-        }
+        IntPoint pt, pt2; //returned by GetOverlapSegment() but unused here.
+        HorzJoinRec* hj = m_HorizJoins[i];
+        //if horizontals rb and hj.edge overlap, flag for joining later ...
+        if (GetOverlapSegment(IntPoint(hj->edge->xbot, hj->edge->ybot),
+          IntPoint(hj->edge->xtop, hj->edge->ytop),
+          IntPoint(rb->xbot, rb->ybot),
+          IntPoint(rb->xtop, rb->ytop), pt, pt2))
+            AddJoin(hj->edge, rb, hj->savedIdx);
       }
     }
 
@@ -1783,7 +1681,7 @@ void Clipper::DeleteFromSEL(TEdge *e)
 //------------------------------------------------------------------------------
 
 void Clipper::IntersectEdges(TEdge *e1, TEdge *e2,
-     const IntPoint &pt, IntersectProtects protects)
+     const IntPoint &pt, const IntersectProtects protects)
 {
   //e1 will be to the left of e2 BELOW the intersection. Therefore e1 is before
   //e2 in AEL except when e1 is being inserted at the intersection point ...
@@ -1983,9 +1881,12 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2)
   OutRec *outRec2 = m_PolyOuts[e2->outIdx];
 
   OutRec *holeStateRec;
-  if (Param1RightOfParam2(outRec1, outRec2)) holeStateRec = outRec2;
-  else if (Param1RightOfParam2(outRec2, outRec1)) holeStateRec = outRec1;
-  else holeStateRec = GetLowermostRec(outRec1, outRec2);
+  if (Param1RightOfParam2(outRec1, outRec2)) 
+    holeStateRec = outRec2;
+  else if (Param1RightOfParam2(outRec2, outRec1)) 
+    holeStateRec = outRec1;
+  else 
+    holeStateRec = GetLowermostRec(outRec1, outRec2);
 
   OutPt* p1_lft = outRec1->pts;
   OutPt* p1_rt = p1_lft->prev;
@@ -1999,7 +1900,7 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2)
     if(  e2->side == esLeft )
     {
       //z y x a b c
-      ReversePolyPtLinks(*p2_lft);
+      ReversePolyPtLinks(p2_lft);
       p2_lft->next = p1_lft;
       p1_lft->prev = p2_lft;
       p1_rt->next = p2_rt;
@@ -2020,7 +1921,7 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2)
     if(  e2->side == esRight )
     {
       //a b c z y x
-      ReversePolyPtLinks( *p2_lft );
+      ReversePolyPtLinks(p2_lft);
       p1_rt->next = p2_rt;
       p2_rt->prev = p1_rt;
       p2_lft->next = p1_lft;
@@ -2046,7 +1947,9 @@ void Clipper::AppendPolygon(TEdge *e1, TEdge *e2)
   }
   outRec2->pts = 0;
   outRec2->bottomPt = 0;
-  outRec2->AppendLink = outRec1;
+
+  outRec2->FirstLeft = outRec1;
+
   int OKIdx = e1->outIdx;
   int ObsoleteIdx = e2->outIdx;
 
@@ -2085,29 +1988,13 @@ OutRec* Clipper::CreateOutRec()
   OutRec* result = new OutRec;
   result->isHole = false;
   result->FirstLeft = 0;
-  result->AppendLink = 0;
   result->pts = 0;
   result->bottomPt = 0;
-  result->sides = esNeither;
-  result->bottomFlag = 0;
-
+  result->polyNode = 0;
   return result;
 }
 //------------------------------------------------------------------------------
 
-void Clipper::DisposeBottomPt(OutRec &outRec)
-{
-  OutPt* next = outRec.bottomPt->next;
-  OutPt* prev = outRec.bottomPt->prev;
-  if (outRec.pts == outRec.bottomPt) outRec.pts = next;
-  delete outRec.bottomPt;
-  next->prev = prev;
-  prev->next = next;
-  outRec.bottomPt = next;
-  FixupOutPolygon(outRec);
-}
-//------------------------------------------------------------------------------
-
 void Clipper::AddOutPt(TEdge *e, const IntPoint &pt)
 {
   bool ToFront = (e->side == esLeft);
@@ -2132,53 +2019,6 @@ void Clipper::AddOutPt(TEdge *e, const IntPoint &pt)
     if ((ToFront && PointsEqual(pt, op->pt)) ||
       (!ToFront && PointsEqual(pt, op->prev->pt))) return;
 
-    if ((e->side | outRec->sides) != outRec->sides)
-    {
-      //check for 'rounding' artefacts ...
-      if (outRec->sides == esNeither && pt.Y == op->pt.Y)
-      {
-        if (ToFront)
-        {
-          if (pt.X == op->pt.X +1) return;    //ie wrong side of bottomPt
-        }
-        else if (pt.X == op->pt.X -1) return; //ie wrong side of bottomPt
-      }
-      outRec->sides = (EdgeSide)(outRec->sides | e->side);
-      if (outRec->sides == esBoth)
-      {
-        //A vertex from each side has now been added.
-        //Vertices of one side of an output polygon are quite commonly close to
-        //or even 'touching' edges of the other side of the output polygon.
-        //Very occasionally vertices from one side can 'cross' an edge on the
-        //the other side. The distance 'crossed' is always less that a unit
-        //and is purely an artefact of coordinate rounding. Nevertheless, this
-        //results in very tiny self-intersections. Because of the way
-        //orientation is calculated, even tiny self-intersections can cause
-        //the Orientation function to return the wrong result. Therefore, it's
-        //important to ensure that any self-intersections close to BottomPt are
-        //detected and removed before orientation is assigned.
-
-        OutPt *opBot, *op2;
-        if (ToFront)
-        {
-          opBot = outRec->pts;
-          op2 = opBot->next; //op2 == right side
-          if (opBot->pt.Y != op2->pt.Y && opBot->pt.Y != pt.Y &&
-            ((opBot->pt.X - pt.X)/(opBot->pt.Y - pt.Y) <
-            (opBot->pt.X - op2->pt.X)/(opBot->pt.Y - op2->pt.Y)))
-               outRec->bottomFlag = opBot;
-        } else
-        {
-          opBot = outRec->pts->prev;
-          op2 = opBot->prev; //op2 == left side
-          if (opBot->pt.Y != op2->pt.Y && opBot->pt.Y != pt.Y &&
-            ((opBot->pt.X - pt.X)/(opBot->pt.Y - pt.Y) >
-            (opBot->pt.X - op2->pt.X)/(opBot->pt.Y - op2->pt.Y)))
-               outRec->bottomFlag = opBot;
-        }
-      }
-    }
-
     OutPt* op2 = new OutPt;
     op2->pt = pt;
     op2->idx = outRec->idx;
@@ -2247,9 +2087,6 @@ TEdge *GetMaximaPair(TEdge *e)
 
 void Clipper::SwapPositionsInAEL(TEdge *edge1, TEdge *edge2)
 {
-  if(  !edge1->nextInAEL &&  !edge1->prevInAEL ) return;
-  if(  !edge2->nextInAEL &&  !edge2->prevInAEL ) return;
-
   if(  edge1->nextInAEL == edge2 )
   {
     TEdge* next = edge2->nextInAEL;
@@ -2474,7 +2311,7 @@ bool Clipper::ProcessIntersections(const long64 botY, const long64 topY)
   try {
     BuildIntersectList(botY, topY);
     if ( !m_IntersectNodes) return true;
-    if ( FixupIntersections() ) ProcessIntersectList();
+    if ( FixupIntersectionOrder() ) ProcessIntersectList();
     else return false;
   }
   catch(...) {
@@ -2503,15 +2340,11 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY)
 
   //prepare for sorting ...
   TEdge* e = m_ActiveEdges;
-  e->tmpX = TopX( *e, topY );
   m_SortedEdges = e;
-  m_SortedEdges->prevInSEL = 0;
-  e = e->nextInAEL;
   while( e )
   {
     e->prevInSEL = e->prevInAEL;
-    e->prevInSEL->nextInSEL = e;
-    e->nextInSEL = 0;
+    e->nextInSEL = e->nextInAEL;
     e->tmpX = TopX( *e, topY );
     e = e->nextInAEL;
   }
@@ -2548,7 +2381,7 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY)
 }
 //------------------------------------------------------------------------------
 
-bool ProcessParam1BeforeParam2(IntersectNode &node1, IntersectNode &node2)
+bool ProcessParam1BeforeParam2(const IntersectNode &node1, const IntersectNode &node2)
 {
   bool result;
   if (node1.pt.Y == node2.pt.Y)
@@ -2618,6 +2451,7 @@ void Clipper::DoMaxima(TEdge *e, long64 topY)
   {
     if (!eNext) throw clipperException("DoMaxima error");
     IntersectEdges( e, eNext, IntPoint(X, topY), ipBoth );
+    SwapPositionsInAEL(e, eNext);
     eNext = eNext->nextInAEL;
   }
   if( e->outIdx < 0 && eMaxPair->outIdx < 0 )
@@ -2778,7 +2612,7 @@ void Clipper::BuildResult(Polygons &polys)
       do
       {
         pg->push_back(p->pt);
-        p = p->next;
+        p = p->prev;
       } while (p != m_PolyOuts[i]->pts);
       //make sure each polygon has at least 3 vertices ...
       if (pg->size() < 3) pg->clear(); else k++;
@@ -2788,39 +2622,58 @@ void Clipper::BuildResult(Polygons &polys)
 }
 //------------------------------------------------------------------------------
 
-void Clipper::BuildResultEx(ExPolygons &polys)
+int PointCount(OutPt *pts)
 {
-  PolyOutList::size_type i = 0;
-  int k = 0;
-  polys.resize(0);
-  polys.reserve(m_PolyOuts.size());
-  while (i < m_PolyOuts.size() && m_PolyOuts[i]->pts)
-  {
-    ExPolygon epg;
-    OutPt* p = m_PolyOuts[i]->pts;
-    do {
-      epg.outer.push_back(p->pt);
-      p = p->next;
-    } while (p != m_PolyOuts[i]->pts);
-    i++;
-    //make sure polygons have at least 3 vertices ...
-    if (epg.outer.size() < 3) continue;
-    while (i < m_PolyOuts.size()
-      && m_PolyOuts[i]->pts && m_PolyOuts[i]->isHole)
+    if (!pts) return 0;
+    int result = 0;
+    OutPt* p = pts;
+    do
     {
-      Polygon pg;
-      p = m_PolyOuts[i]->pts;
-      do {
-        pg.push_back(p->pt);
+        result++;
         p = p->next;
-      } while (p != m_PolyOuts[i]->pts);
-      epg.holes.push_back(pg);
-      i++;
     }
-    polys.push_back(epg);
-    k++;
-  }
-  polys.resize(k);
+    while (p != pts);
+    return result;
+}
+//------------------------------------------------------------------------------
+
+void Clipper::BuildResult2(PolyTree& polytree)
+{
+    polytree.Clear();
+    polytree.AllNodes.reserve(m_PolyOuts.size());
+    //add each output polygon/contour to polytree ...
+    for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); i++)
+    {
+        OutRec* outRec = m_PolyOuts[i];
+        int cnt = PointCount(outRec->pts);
+        if (cnt < 3) continue;
+        FixHoleLinkage(*outRec);
+        PolyNode* pn = new PolyNode();
+        //nb: polytree takes ownership of all the PolyNodes
+        polytree.AllNodes.push_back(pn);
+        outRec->polyNode = pn;
+        pn->Parent = 0;
+        pn->Index = 0;
+        pn->Contour.reserve(cnt);
+        OutPt *op = outRec->pts;
+        for (int j = 0; j < cnt; j++)
+        {
+            pn->Contour.push_back(op->pt);
+            op = op->prev;
+        }
+    }
+
+    //fixup PolyNode links etc ...
+    polytree.Childs.reserve(m_PolyOuts.size());
+    for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); i++)
+    {
+        OutRec* outRec = m_PolyOuts[i];
+        if (!outRec->polyNode) continue;
+        if (outRec->FirstLeft) 
+          outRec->FirstLeft->polyNode->AddChild(*outRec->polyNode);
+        else
+          polytree.AddChild(*outRec->polyNode);
+    }
 }
 //------------------------------------------------------------------------------
 
@@ -2840,7 +2693,7 @@ void SwapIntersectNodes(IntersectNode &int1, IntersectNode &int2)
 }
 //------------------------------------------------------------------------------
 
-bool Clipper::FixupIntersections()
+bool Clipper::FixupIntersectionOrder()
 {
   if ( !m_IntersectNodes->next ) return true;
 
@@ -2855,8 +2708,8 @@ bool Clipper::FixupIntersections()
     else if (e1->nextInSEL == int1->edge2) e2 = e1->nextInSEL;
     else
     {
-      //The current intersection is out of order, so try and swap it with
-      //a subsequent intersection ...
+      //The current intersection (Int1) is out of order (since it doesn't
+      //contain adjacent edges), so swap it with a subsequent intersection ...
       while (int2)
       {
         if (int2->edge1->nextInSEL == int2->edge2 ||
@@ -2865,7 +2718,8 @@ bool Clipper::FixupIntersections()
       }
       if ( !int2 ) return false; //oops!!!
 
-      //found an intersect node that can be swapped ...
+      //found an intersect node (Int2) that does contain adjacent edges,
+      //so prepare to process it before Int1 ...
       SwapIntersectNodes(*int1, *int2);
       e1 = int1->edge1;
       e2 = int1->edge2;
@@ -2940,77 +2794,162 @@ void Clipper::DoBothEdges(TEdge *edge1, TEdge *edge2, const IntPoint &pt)
 }
 //----------------------------------------------------------------------
 
-void Clipper::JoinCommonEdges(bool fixHoleLinkages)
+bool Clipper::JoinPoints(const JoinRec *j, OutPt *&p1, OutPt *&p2)
+{
+  OutRec *outRec1 = m_PolyOuts[j->poly1Idx];
+  OutRec *outRec2 = m_PolyOuts[j->poly2Idx];
+  if (!outRec1 || !outRec2)  return false;  
+  OutPt *pp1a = outRec1->pts;
+  OutPt *pp2a = outRec2->pts;
+  IntPoint pt1 = j->pt2a, pt2 = j->pt2b;
+  IntPoint pt3 = j->pt1a, pt4 = j->pt1b;
+  if (!FindSegment(pp1a, m_UseFullRange, pt1, pt2)) return false;
+  if (outRec1 == outRec2)
+  {
+    //we're searching the same polygon for overlapping segments so
+    //segment 2 mustn't be the same as segment 1 ...
+    pp2a = pp1a->next;
+    if (!FindSegment(pp2a, m_UseFullRange, pt3, pt4) || (pp2a == pp1a)) 
+      return false;
+  }
+  else if (!FindSegment(pp2a, m_UseFullRange, pt3, pt4)) return false;
+
+  if (!GetOverlapSegment(pt1, pt2, pt3, pt4, pt1, pt2)) return false;
+
+  OutPt *p3, *p4, *prev = pp1a->prev;
+  //get p1 & p2 polypts - the overlap start & endpoints on poly1
+  if (PointsEqual(pp1a->pt, pt1)) p1 = pp1a;
+  else if (PointsEqual(prev->pt, pt1)) p1 = prev;
+  else p1 = InsertPolyPtBetween(pp1a, prev, pt1);
+
+  if (PointsEqual(pp1a->pt, pt2)) p2 = pp1a;
+  else if (PointsEqual(prev->pt, pt2)) p2 = prev;
+  else if ((p1 == pp1a) || (p1 == prev))
+    p2 = InsertPolyPtBetween(pp1a, prev, pt2);
+  else if (Pt3IsBetweenPt1AndPt2(pp1a->pt, p1->pt, pt2))
+    p2 = InsertPolyPtBetween(pp1a, p1, pt2); else
+    p2 = InsertPolyPtBetween(p1, prev, pt2);
+
+  //get p3 & p4 polypts - the overlap start & endpoints on poly2
+  prev = pp2a->prev;
+  if (PointsEqual(pp2a->pt, pt1)) p3 = pp2a;
+  else if (PointsEqual(prev->pt, pt1)) p3 = prev;
+  else p3 = InsertPolyPtBetween(pp2a, prev, pt1);
+
+  if (PointsEqual(pp2a->pt, pt2)) p4 = pp2a;
+  else if (PointsEqual(prev->pt, pt2)) p4 = prev;
+  else if ((p3 == pp2a) || (p3 == prev))
+    p4 = InsertPolyPtBetween(pp2a, prev, pt2);
+  else if (Pt3IsBetweenPt1AndPt2(pp2a->pt, p3->pt, pt2))
+    p4 = InsertPolyPtBetween(pp2a, p3, pt2); else
+    p4 = InsertPolyPtBetween(p3, prev, pt2);
+
+  //p1.pt == p3.pt and p2.pt == p4.pt so join p1 to p3 and p2 to p4 ...
+  if (p1->next == p2 && p3->prev == p4)
+  {
+    p1->next = p3;
+    p3->prev = p1;
+    p2->prev = p4;
+    p4->next = p2;
+    return true;
+  }
+  else if (p1->prev == p2 && p3->next == p4)
+  {
+    p1->prev = p3;
+    p3->next = p1;
+    p2->next = p4;
+    p4->prev = p2;
+    return true;
+  }
+  else
+    return false; //an orientation is probably wrong
+}
+//----------------------------------------------------------------------
+
+void Clipper::FixupJoinRecs(JoinRec *j, OutPt *pt, unsigned startIdx)
+{
+  for (JoinList::size_type k = startIdx; k < m_Joins.size(); k++)
+    {
+      JoinRec* j2 = m_Joins[k];
+      if (j2->poly1Idx == j->poly1Idx && PointIsVertex(j2->pt1a, pt))
+        j2->poly1Idx = j->poly2Idx;
+      if (j2->poly2Idx == j->poly1Idx && PointIsVertex(j2->pt2a, pt))
+        j2->poly2Idx = j->poly2Idx;
+    }
+}
+//----------------------------------------------------------------------
+
+bool Poly2ContainsPoly1(OutPt* outPt1, OutPt* outPt2, bool UseFullInt64Range)
+{
+  //find the first pt in outPt1 that isn't also a vertex of outPt2 ...
+  OutPt* outPt = outPt1;
+  do
+  {
+      if (!PointIsVertex(outPt->pt, outPt2)) break;
+      outPt = outPt->next;
+  }
+  while (outPt != outPt1);
+  bool result;
+  //sometimes a point on one polygon can be touching the other polygon 
+  //so to be totally confident outPt1 is inside outPt2 repeat ...
+  do
+  {
+    result = PointInPolygon(outPt->pt, outPt2, UseFullInt64Range);
+    outPt = outPt->next;
+  }
+  while (result && outPt != outPt1);
+  return result;
+}
+//----------------------------------------------------------------------
+
+void Clipper::FixupFirstLefts1(OutRec* OldOutRec, OutRec* NewOutRec)
+{ 
+  
+  for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i)
+  {
+    OutRec* outRec = m_PolyOuts[i];
+    if (outRec->pts && outRec->FirstLeft == OldOutRec) 
+    {
+      if (Poly2ContainsPoly1(outRec->pts, NewOutRec->pts, m_UseFullRange))
+        outRec->FirstLeft = NewOutRec;
+    }
+  }
+}
+//----------------------------------------------------------------------
+
+void Clipper::FixupFirstLefts2(OutRec* OldOutRec, OutRec* NewOutRec)
+{ 
+  for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i)
+  {
+    OutRec* outRec = m_PolyOuts[i];
+    if (outRec->FirstLeft == OldOutRec) outRec->FirstLeft = NewOutRec;
+  }
+}
+//----------------------------------------------------------------------
+
+void Clipper::JoinCommonEdges()
 {
   for (JoinList::size_type i = 0; i < m_Joins.size(); i++)
   {
     JoinRec* j = m_Joins[i];
+
     OutRec *outRec1 = m_PolyOuts[j->poly1Idx];
-    OutPt *pp1a = outRec1->pts;
     OutRec *outRec2 = m_PolyOuts[j->poly2Idx];
-    OutPt *pp2a = outRec2->pts;
-    IntPoint pt1 = j->pt2a, pt2 = j->pt2b;
-    IntPoint pt3 = j->pt1a, pt4 = j->pt1b;
-    if (!FindSegment(pp1a, pt1, pt2)) continue;
-    if (j->poly1Idx == j->poly2Idx)
-    {
-      //we're searching the same polygon for overlapping segments so
-      //segment 2 mustn't be the same as segment 1 ...
-      pp2a = pp1a->next;
-      if (!FindSegment(pp2a, pt3, pt4) || (pp2a == pp1a)) continue;
-    }
-    else if (!FindSegment(pp2a, pt3, pt4)) continue;
 
-    if (!GetOverlapSegment(pt1, pt2, pt3, pt4, pt1, pt2)) continue;
+    if (!outRec1->pts || !outRec2->pts) continue;
 
-    OutPt *p1, *p2, *p3, *p4;
-    OutPt *prev = pp1a->prev;
-    //get p1 & p2 polypts - the overlap start & endpoints on poly1
-    if (PointsEqual(pp1a->pt, pt1)) p1 = pp1a;
-    else if (PointsEqual(prev->pt, pt1)) p1 = prev;
-    else p1 = InsertPolyPtBetween(pp1a, prev, pt1);
+    //get the polygon fragment with the correct hole state (FirstLeft)
+    //before calling JoinPoints() ...
+    OutRec *holeStateRec;
+    if (outRec1 == outRec2) holeStateRec = outRec1;
+    else if (Param1RightOfParam2(outRec1, outRec2)) holeStateRec = outRec2;
+    else if (Param1RightOfParam2(outRec2, outRec1)) holeStateRec = outRec1;
+    else holeStateRec = GetLowermostRec(outRec1, outRec2);
 
-    if (PointsEqual(pp1a->pt, pt2)) p2 = pp1a;
-    else if (PointsEqual(prev->pt, pt2)) p2 = prev;
-    else if ((p1 == pp1a) || (p1 == prev))
-      p2 = InsertPolyPtBetween(pp1a, prev, pt2);
-    else if (Pt3IsBetweenPt1AndPt2(pp1a->pt, p1->pt, pt2))
-      p2 = InsertPolyPtBetween(pp1a, p1, pt2); else
-      p2 = InsertPolyPtBetween(p1, prev, pt2);
+    OutPt *p1, *p2;
+    if (!JoinPoints(j, p1, p2)) continue;
 
-    //get p3 & p4 polypts - the overlap start & endpoints on poly2
-    prev = pp2a->prev;
-    if (PointsEqual(pp2a->pt, pt1)) p3 = pp2a;
-    else if (PointsEqual(prev->pt, pt1)) p3 = prev;
-    else p3 = InsertPolyPtBetween(pp2a, prev, pt1);
-
-    if (PointsEqual(pp2a->pt, pt2)) p4 = pp2a;
-    else if (PointsEqual(prev->pt, pt2)) p4 = prev;
-    else if ((p3 == pp2a) || (p3 == prev))
-      p4 = InsertPolyPtBetween(pp2a, prev, pt2);
-    else if (Pt3IsBetweenPt1AndPt2(pp2a->pt, p3->pt, pt2))
-      p4 = InsertPolyPtBetween(pp2a, p3, pt2); else
-      p4 = InsertPolyPtBetween(p3, prev, pt2);
-
-    //p1.pt == p3.pt and p2.pt == p4.pt so join p1 to p3 and p2 to p4 ...
-    if (p1->next == p2 && p3->prev == p4)
-    {
-      p1->next = p3;
-      p3->prev = p1;
-      p2->prev = p4;
-      p4->next = p2;
-    }
-    else if (p1->prev == p2 && p3->next == p4)
-    {
-      p1->prev = p3;
-      p3->next = p1;
-      p2->next = p4;
-      p4->prev = p2;
-    }
-    else
-      continue; //an orientation is probably wrong
-
-    if (j->poly2Idx == j->poly1Idx)
+    if (outRec1 == outRec2)
     {
       //instead of joining two polygons, we've just created a new one by
       //splitting one polygon into two.
@@ -3025,111 +2964,75 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages)
       outRec2->bottomPt = outRec2->pts;
       outRec2->bottomPt->idx = outRec2->idx;
 
-      int K = 0;
-      if (PointInPolygon(outRec2->pts->pt, outRec1->pts, m_UseFullRange))
+      if (Poly2ContainsPoly1(outRec2->pts, outRec1->pts, m_UseFullRange))
       {
         //outRec2 is contained by outRec1 ...
-        K = 1;
         outRec2->isHole = !outRec1->isHole;
         outRec2->FirstLeft = outRec1;
-      } else if (PointInPolygon(outRec1->pts->pt, outRec2->pts, m_UseFullRange))
+
+        FixupJoinRecs(j, p2, i+1);
+
+        //fixup FirstLeft pointers that may need reassigning to OutRec1
+        if (m_UsingPolyTree) FixupFirstLefts2(outRec2, outRec1);
+
+        FixupOutPolygon(*outRec1); //nb: do this BEFORE testing orientation
+        FixupOutPolygon(*outRec2); //    but AFTER calling FixupJoinRecs()
+
+
+        if ((outRec2->isHole ^ m_ReverseOutput) == (Area(*outRec2, m_UseFullRange) > 0))
+          ReversePolyPtLinks(outRec2->pts);
+            
+      } else if (Poly2ContainsPoly1(outRec1->pts, outRec2->pts, m_UseFullRange))
       {
         //outRec1 is contained by outRec2 ...
-        K = 2;
         outRec2->isHole = outRec1->isHole;
         outRec1->isHole = !outRec2->isHole;
         outRec2->FirstLeft = outRec1->FirstLeft;
         outRec1->FirstLeft = outRec2;
-      } else
+
+        FixupJoinRecs(j, p2, i+1);
+
+        //fixup FirstLeft pointers that may need reassigning to OutRec1
+        if (m_UsingPolyTree) FixupFirstLefts2(outRec1, outRec2);
+
+        FixupOutPolygon(*outRec1); //nb: do this BEFORE testing orientation
+        FixupOutPolygon(*outRec2); //    but AFTER calling FixupJoinRecs()
+
+        if ((outRec1->isHole ^ m_ReverseOutput) == (Area(*outRec1, m_UseFullRange) > 0))
+          ReversePolyPtLinks(outRec1->pts);
+      } 
+      else
       {
+        //the 2 polygons are completely separate ...
         outRec2->isHole = outRec1->isHole;
         outRec2->FirstLeft = outRec1->FirstLeft;
-      }
 
-      //now fixup any subsequent joins that match this polygon
-      for (JoinList::size_type k = i+1; k < m_Joins.size(); k++)
-      {
-        JoinRec* j2 = m_Joins[k];
-        if (j2->poly1Idx == j->poly1Idx && PointIsVertex(j2->pt1a, p2))
-          j2->poly1Idx = j->poly2Idx;
-        if (j2->poly2Idx == j->poly1Idx && PointIsVertex(j2->pt2a, p2))
-          j2->poly2Idx = j->poly2Idx;
-      }
+        FixupJoinRecs(j, p2, i+1);
 
-      FixupOutPolygon(*outRec1); //nb: do this BEFORE testing orientation
-      FixupOutPolygon(*outRec2); //    but AFTER calling PointIsVertex()
+        //fixup FirstLeft pointers that may need reassigning to OutRec2
+        if (m_UsingPolyTree) FixupFirstLefts1(outRec1, outRec2);
 
-      switch( K ) {
-        case 1: 
-          {
-            if (outRec2->isHole ==
-              (m_ReverseOutput ^ Orientation(outRec2, m_UseFullRange)))
-                ReversePolyPtLinks(*outRec2->pts);
-            break;
-          }
-        case 2: 
-          {
-            if (outRec1->isHole ==
-              (m_ReverseOutput ^ Orientation(outRec1, m_UseFullRange)))
-                ReversePolyPtLinks(*outRec1->pts);
-            //make sure any contained holes now link to the correct polygon ...
-            if (fixHoleLinkages && outRec1->isHole) 
-              for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k)
-              {
-                OutRec *orec = m_PolyOuts[k];
-                if (orec->isHole && orec->bottomPt && orec->FirstLeft == outRec1)
-                  orec->FirstLeft = outRec2;
-              }
-            break;
-          }
-        default: 
-          {
-            //make sure any contained holes now link to the correct polygon ...
-            if (fixHoleLinkages) 
-              for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k)
-              {
-                OutRec *orec = m_PolyOuts[k];
-                if (orec->isHole && orec->bottomPt && orec->FirstLeft == outRec1 &&
-                  !PointInPolygon(orec->bottomPt->pt, outRec1->pts, m_UseFullRange))
-                    orec->FirstLeft = outRec2;
-              }
-
-          }
+        FixupOutPolygon(*outRec1); //nb: do this BEFORE testing orientation
+        FixupOutPolygon(*outRec2); //    but AFTER calling FixupJoinRecs()
       }
      
-      if (Orientation(outRec1, m_UseFullRange) != (Area(*outRec1, m_UseFullRange) > 0))
-          DisposeBottomPt(*outRec1);
-      if (Orientation(outRec2, m_UseFullRange) != (Area(*outRec2, m_UseFullRange) > 0))
-          DisposeBottomPt(*outRec2);
-
     } else
     {
       //joined 2 polygons together ...
 
-      //make sure any holes contained by outRec2 now link to outRec1 ...
-      if (fixHoleLinkages) 
-        for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k)
-          if (m_PolyOuts[k]->isHole && m_PolyOuts[k]->bottomPt &&
-            m_PolyOuts[k]->FirstLeft == outRec2)
-              m_PolyOuts[k]->FirstLeft = outRec1;
-
-
-      //and cleanup redundant edges too ...
+      //cleanup redundant edges ...
       FixupOutPolygon(*outRec1);
 
-      if (outRec1->pts)
-      {
-        outRec1->isHole = !Orientation(outRec1, m_UseFullRange);
-        if (outRec1->isHole && !outRec1->FirstLeft)
-          outRec1->FirstLeft = outRec2->FirstLeft;
-      }
-
       //delete the obsolete pointer ...
       int OKIdx = outRec1->idx;
       int ObsoleteIdx = outRec2->idx;
       outRec2->pts = 0;
       outRec2->bottomPt = 0;
-      outRec2->AppendLink = outRec1;
+
+      outRec1->isHole = holeStateRec->isHole;
+      if (holeStateRec == outRec2) 
+        outRec1->FirstLeft = outRec2->FirstLeft;
+      outRec2->FirstLeft = outRec1;
 
       //now fixup any subsequent Joins that match this polygon
       for (JoinList::size_type k = i+1; k < m_Joins.size(); k++)
@@ -3138,6 +3041,9 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages)
         if (j2->poly1Idx == ObsoleteIdx) j2->poly1Idx = OKIdx;
         if (j2->poly2Idx == ObsoleteIdx) j2->poly2Idx = OKIdx;
       }
+
+      //fixup FirstLeft pointers that may need reassigning to OutRec1
+      if (m_UsingPolyTree) FixupFirstLefts2(outRec2, outRec1);
     }
   }
 }
@@ -3168,25 +3074,32 @@ struct DoublePoint
 //------------------------------------------------------------------------------
 
 Polygon BuildArc(const IntPoint &pt,
-  const double a1, const double a2, const double r)
+  const double a1, const double a2, const double r, double limit)
 {
-  long64 steps = std::max(6, int(std::sqrt(std::fabs(r)) * std::fabs(a2 - a1)));
-  if (steps > 0x100000) steps = 0x100000;
-  int n = (unsigned)steps;
-  Polygon result(n);
-  double da = (a2 - a1) / (n -1);
-  double a = a1;
-  for (int i = 0; i < n; ++i)
+  //see notes in clipper.pas regarding steps
+  double arcFrac = std::fabs(a2 - a1) / (2 * pi);
+  int steps = (int)(arcFrac * pi / std::acos(1 - limit / std::fabs(r)));
+  if (steps < 2) steps = 2;
+  else if (steps > (int)(222.0 * arcFrac)) steps = (int)(222.0 * arcFrac);
+
+  double x = std::cos(a1); 
+  double y = std::sin(a1);
+  double c = std::cos((a2 - a1) / steps);
+  double s = std::sin((a2 - a1) / steps);
+  Polygon result(steps +1);
+  for (int i = 0; i <= steps; ++i)
   {
-    result[i].X = pt.X + Round(std::cos(a)*r);
-    result[i].Y = pt.Y + Round(std::sin(a)*r);
-    a += da;
+      result[i].X = pt.X + Round(x * r);
+      result[i].Y = pt.Y + Round(y * r);
+      double x2 = x;
+      x = x * c - s * y;  //cross product
+      y = x2 * s + y * c; //dot product
   }
   return result;
 }
 //------------------------------------------------------------------------------
 
-DoublePoint GetUnitNormal( const IntPoint &pt1, const IntPoint &pt2)
+DoublePoint GetUnitNormal(const IntPoint &pt1, const IntPoint &pt2)
 {
   if(pt2.X == pt1.X && pt2.Y == pt1.Y) 
     return DoublePoint(0, 0);
@@ -3216,7 +3129,7 @@ private:
 public:
 
 PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
-  double delta, JoinType jointype, double MiterLimit)
+  double delta, JoinType jointype, double limit, bool autoFix)
 {
     //nb precondition - out_polys != ptsin_polys
     if (NEAR_ZERO(delta))
@@ -3228,22 +3141,68 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
     this->m_p = in_polys;
     this->m_delta = delta;
     this->m_jointype = jointype;
-    if (MiterLimit <= 1) MiterLimit = 1;
-    m_RMin = 2/(MiterLimit*MiterLimit);
+
+    //ChecksInput - fixes polygon orientation if necessary and removes 
+    //duplicate vertices. Can be set false when you're sure that polygon
+    //orientation is correct and that there are no duplicate vertices.
+    if (autoFix) 
+    {
+      size_t Len = m_p.size(), botI = 0;
+      while (botI < Len && m_p[botI].size() == 0) botI++;
+      if (botI == Len) return;
+      
+      //botPt: used to find the lowermost (in inverted Y-axis) & leftmost point
+      //This point (on m_p[botI]) must be on an outer polygon ring and if 
+      //its orientation is false (counterclockwise) then assume all polygons 
+      //need reversing ...
+      IntPoint botPt = m_p[botI][0];      
+      for (size_t i = botI; i < Len; ++i)
+      {
+        if (m_p[i].size() < 3) continue;
+        if (UpdateBotPt(m_p[i][0], botPt)) botI = i;
+        Polygon::iterator it = m_p[i].begin() +1;
+        while (it != m_p[i].end())
+        {
+          if (PointsEqual(*it, *(it -1)))
+            it = m_p[i].erase(it);
+          else 
+          {
+            if (UpdateBotPt(*it, botPt)) botI = i;
+            ++it;
+          }
+        }
+      }
+      if (!Orientation(m_p[botI]))
+        ReversePolygons(m_p);
+    }
+
+    switch (jointype)
+    {
+        case jtRound:  //limit defaults to 0.125
+            if (limit <= 0) limit = 0.125;
+            else if (limit > std::fabs(delta)) limit = std::fabs(delta);
+            break;  
+        case jtMiter:  //limit defaults to twice delta's width ...
+            if (limit < 2) limit = 2; 
+            break;       
+        default:       //otherwise limit is unused
+            limit = 1;   
+    }
+    m_RMin = 2.0/(limit*limit);
  
     double deltaSq = delta*delta;
     out_polys.clear();
-    out_polys.resize(in_polys.size());
-    for (m_i = 0; m_i < in_polys.size(); m_i++)
+    out_polys.resize(m_p.size());
+    for (m_i = 0; m_i < m_p.size(); m_i++)
     {
         m_curr_poly = &out_polys[m_i];
-        size_t len = in_polys[m_i].size();
+        size_t len = m_p[m_i].size();
         if (len > 1 && m_p[m_i][0].X == m_p[m_i][len - 1].X &&
             m_p[m_i][0].Y == m_p[m_i][len-1].Y) len--;
 
         //when 'shrinking' polygons - to minimize artefacts
         //strip those polygons that have an area < pi * delta^2 ...
-        double a1 = Area(in_polys[m_i]);
+        double a1 = Area(m_p[m_i]);
         if (delta < 0) { if (a1 > 0 && a1 < deltaSq *pi) len = 0; }
         else if (a1 < 0 && -a1 < deltaSq *pi) len = 0; //holes have neg. area
 
@@ -3252,7 +3211,7 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
         else if (len == 1)
         {
             Polygon arc;
-            arc = BuildArc(in_polys[m_i][len-1], 0, 2 * pi, delta);
+            arc = BuildArc(m_p[m_i][len-1], 0, 2 * pi, delta, limit);
             out_polys[m_i] = arc;
             continue;
         }
@@ -3260,9 +3219,9 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
         //build normals ...
         normals.clear();
         normals.resize(len);
-        normals[len-1] = GetUnitNormal(in_polys[m_i][len-1], in_polys[m_i][0]);
+        normals[len-1] = GetUnitNormal(m_p[m_i][len-1], m_p[m_i][0]);
         for (m_j = 0; m_j < len -1; ++m_j)
-            normals[m_j] = GetUnitNormal(in_polys[m_i][m_j], in_polys[m_i][m_j+1]);
+            normals[m_j] = GetUnitNormal(m_p[m_i][m_j], m_p[m_i][m_j+1]);
         
         m_k = len -1;
         for (m_j = 0; m_j < len; ++m_j)
@@ -3273,11 +3232,11 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
             {
               m_R = 1 + (normals[m_j].X*normals[m_k].X + 
                 normals[m_j].Y*normals[m_k].Y);
-              if (m_R >= m_RMin) DoMiter(); else DoSquare(MiterLimit);
+              if (m_R >= m_RMin) DoMiter(); else DoSquare(limit);
               break;
             }
-            case jtSquare: DoSquare(); break;
-            case jtRound: DoRound(); break;
+            case jtSquare: DoSquare(1.0); break;
+            case jtRound: DoRound(limit); break;
           }
         m_k = m_j;
         }
@@ -3323,7 +3282,7 @@ void AddPoint(const IntPoint& pt)
 }
 //------------------------------------------------------------------------------
 
-void DoSquare(double mul = 1.0)
+void DoSquare(double mul)
 {
     IntPoint pt1 = IntPoint((long64)Round(m_p[m_i][m_j].X + normals[m_k].X * m_delta),
         (long64)Round(m_p[m_i][m_j].Y + normals[m_k].Y * m_delta));
@@ -3335,7 +3294,7 @@ void DoSquare(double mul = 1.0)
       double a2 = std::atan2(-normals[m_j].Y, -normals[m_j].X);
       a1 = std::fabs(a2 - a1);
       if (a1 > pi) a1 = pi * 2 - a1;
-      double dx = std::tan((pi - a1)/4) * std::fabs(m_delta * mul);
+      double dx = std::tan((pi - a1) / 4) * std::fabs(m_delta * mul);
       pt1 = IntPoint((long64)(pt1.X -normals[m_k].Y * dx),
         (long64)(pt1.Y + normals[m_k].X * dx));
       AddPoint(pt1);
@@ -3374,7 +3333,7 @@ void DoMiter()
 }
 //------------------------------------------------------------------------------
 
-void DoRound()
+void DoRound(double limit)
 {
     IntPoint pt1 = IntPoint((long64)Round(m_p[m_i][m_j].X + normals[m_k].X * m_delta),
         (long64)Round(m_p[m_i][m_j].Y + normals[m_k].Y * m_delta));
@@ -3390,7 +3349,7 @@ void DoRound()
         double a2 = std::atan2(normals[m_j].Y, normals[m_j].X);
         if (m_delta > 0 && a2 < a1) a2 += pi *2;
         else if (m_delta < 0 && a2 > a1) a2 -= pi *2;
-        Polygon arc = BuildArc(m_p[m_i][m_j], a1, a2, m_delta);
+        Polygon arc = BuildArc(m_p[m_i][m_j], a1, a2, m_delta, limit);
         for (Polygon::size_type m = 0; m < arc.size(); m++)
           AddPoint(arc[m]);
       }
@@ -3401,20 +3360,31 @@ void DoRound()
 }
 //--------------------------------------------------------------------------
 
+bool UpdateBotPt(const IntPoint &pt, IntPoint &botPt)
+{
+    if (pt.Y > botPt.Y || (pt.Y == botPt.Y && pt.X < botPt.X))
+    {
+        botPt = pt;
+        return true;
+    }
+    else return false;
+}
+//--------------------------------------------------------------------------
+
 }; //end PolyOffsetBuilder
 
 //------------------------------------------------------------------------------
 //------------------------------------------------------------------------------
 
 void OffsetPolygons(const Polygons &in_polys, Polygons &out_polys,
-  double delta, JoinType jointype, double MiterLimit)
+  double delta, JoinType jointype, double limit, bool autoFix)
 {
   if (&out_polys == &in_polys)
   {
     Polygons poly2(in_polys);
-    PolyOffsetBuilder(poly2, out_polys, delta, jointype, MiterLimit);
+    PolyOffsetBuilder(poly2, out_polys, delta, jointype, limit, autoFix);
   }
-  else PolyOffsetBuilder(in_polys, out_polys, delta, jointype, MiterLimit);
+  else PolyOffsetBuilder(in_polys, out_polys, delta, jointype, limit, autoFix);
 }
 //------------------------------------------------------------------------------
 
@@ -3440,6 +3410,62 @@ void SimplifyPolygons(Polygons &polys, PolyFillType fillType)
 }
 //------------------------------------------------------------------------------
 
+void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance)
+{
+  //delta = proximity in units/pixels below which vertices
+  //will be stripped. Default ~= sqrt(2) so when adjacent 
+  //vertices have both x & y coords within 1 unit, then 
+  //the second vertex will be stripped. 
+  int len = in_poly.size();
+  if (len < 3) 
+    out_poly.resize(0);
+  else
+    out_poly.resize(in_poly.size());
+
+  int d = (int)(distance * distance);
+  IntPoint p = in_poly[0];
+  int j = 1;
+  for (int i = 1; i < len; i++)
+  {
+      if ((in_poly[i].X - p.X) * (in_poly[i].X - p.X) +
+          (in_poly[i].Y - p.Y) * (in_poly[i].Y - p.Y) <= d)
+          continue;
+      out_poly[j] = in_poly[i];
+      p = in_poly[i];
+      j++;
+  }
+  p = in_poly[j - 1];
+  if ((in_poly[0].X - p.X) * (in_poly[0].X - p.X) +
+      (in_poly[0].Y - p.Y) * (in_poly[0].Y - p.Y) <= d)
+      j--;
+  if (j < len)
+    out_poly.resize(j);
+}
+//------------------------------------------------------------------------------
+void CleanPolygons(Polygons& in_polys, Polygons& out_polys, double distance)
+{
+  for (Polygons::size_type i = 0; i < in_polys.size(); ++i)
+    CleanPolygon(in_polys[i], out_polys[i], distance);
+}
+//------------------------------------------------------------------------------
+
+void AddPolyNodeToPolygons(PolyNode& polynode, Polygons& polygons)
+{
+  if (polynode.Contour.size() > 0)
+    polygons.push_back(polynode.Contour);
+  for (int i = 0; i < polynode.ChildCount(); ++i)
+    AddPolyNodeToPolygons(*polynode.Childs[i], polygons);
+}
+//------------------------------------------------------------------------------
+
+void PolyTreeToPolygons(PolyTree& polytree, Polygons& polygons)
+{
+  polygons.resize(0); 
+  polygons.reserve(polytree.Total());
+  AddPolyNodeToPolygons(polytree, polygons);
+}
+//------------------------------------------------------------------------------
+
 std::ostream& operator <<(std::ostream &s, IntPoint& p)
 {
   s << p.X << ' ' << p.Y << "\n";
diff --git a/src/Lib/terragear/clipper.hpp b/src/Lib/terragear/clipper.hpp
index a3427f74..1910ec7f 100644
--- a/src/Lib/terragear/clipper.hpp
+++ b/src/Lib/terragear/clipper.hpp
@@ -1,10 +1,10 @@
 /*******************************************************************************
 *                                                                              *
 * Author    :  Angus Johnson                                                   *
-* Version   :  4.9.1                                                           *
-* Date      :  9 October 2012                                                  *
+* Version   :  5.1.2                                                           *
+* Date      :  25 February 2013                                                *
 * Website   :  http://www.angusj.com                                           *
-* Copyright :  Angus Johnson 2010-2012                                         *
+* Copyright :  Angus Johnson 2010-2013                                         *
 *                                                                              *
 * License:                                                                     *
 * Use, modification & distribution is subject to Boost Software License Ver 1. *
@@ -26,7 +26,7 @@
 * Paper no. DETC2005-85513 pp. 565-575                                         *
 * ASME 2005 International Design Engineering Technical Conferences             *
 * and Computers and Information in Engineering Conference (IDETC/CIE2005)      *
-* September 24�28, 2005 , Long Beach, California, USA                          *
+* September 24-28, 2005 , Long Beach, California, USA                          *
 * http://www.me.berkeley.edu/~mcmains/pubs/DAC05OffsetPolygon.pdf              *
 *                                                                              *
 *******************************************************************************/
@@ -64,30 +64,64 @@ public:
 typedef std::vector< IntPoint > Polygon;
 typedef std::vector< Polygon > Polygons;
 
+
 std::ostream& operator <<(std::ostream &s, Polygon &p);
 std::ostream& operator <<(std::ostream &s, Polygons &p);
 
-struct ExPolygon {
-  Polygon  outer;
-  Polygons holes;
-};
-typedef std::vector< ExPolygon > ExPolygons;
+class PolyNode;
+typedef std::vector< PolyNode* > PolyNodes;
 
+class PolyNode 
+{ 
+public:
+    PolyNode();
+    Polygon Contour;
+    PolyNodes Childs;
+    PolyNode* Parent;
+    PolyNode* GetNext() const;
+    bool IsHole() const;
+    int ChildCount() const;
+private:
+    PolyNode* GetNextSiblingUp() const;
+    unsigned Index; //node index in Parent.Childs
+    void AddChild(PolyNode& child);
+    friend class Clipper; //to access Index
+};
+
+class PolyTree: public PolyNode
+{ 
+public:
+    ~PolyTree(){Clear();};
+    PolyNode* GetFirst() const;
+    void Clear();
+    int Total() const;
+private:
+    PolyNodes AllNodes;
+    friend class Clipper; //to access AllNodes
+};
+        
 enum JoinType { jtSquare, jtRound, jtMiter };
 
 bool Orientation(const Polygon &poly);
 double Area(const Polygon &poly);
+
 void OffsetPolygons(const Polygons &in_polys, Polygons &out_polys,
-  double delta, JoinType jointype = jtSquare, double MiterLimit = 2);
+  double delta, JoinType jointype = jtSquare, double limit = 0, bool autoFix = true);
+
 void SimplifyPolygon(const Polygon &in_poly, Polygons &out_polys, PolyFillType fillType = pftEvenOdd);
 void SimplifyPolygons(const Polygons &in_polys, Polygons &out_polys, PolyFillType fillType = pftEvenOdd);
 void SimplifyPolygons(Polygons &polys, PolyFillType fillType = pftEvenOdd);
 
+void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance = 1.415);
+void CleanPolygons(Polygons& in_polys, Polygons& out_polys, double distance = 1.415);
+
+void PolyTreeToPolygons(PolyTree& polytree, Polygons& polygons);
+
 void ReversePolygon(Polygon& p);
 void ReversePolygons(Polygons& p);
 
 //used internally ...
-enum EdgeSide { esNeither = 0, esLeft = 1, esRight = 2, esBoth = 3 };
+enum EdgeSide { esLeft = 1, esRight = 2};
 enum IntersectProtects { ipNone = 0, ipLeft = 1, ipRight = 2, ipBoth = 3 };
 
 struct TEdge {
@@ -140,12 +174,10 @@ struct OutPt; //forward declaration
 struct OutRec {
   int     idx;
   bool    isHole;
-  OutRec *FirstLeft;
-  OutRec *AppendLink;
+  OutRec *FirstLeft;  //see comments in clipper.pas
+  PolyNode *polyNode;
   OutPt  *pts;
   OutPt  *bottomPt;
-  OutPt  *bottomFlag;
-  EdgeSide sides;
 };
 
 struct OutPt {
@@ -206,19 +238,19 @@ public:
   Clipper();
   ~Clipper();
   bool Execute(ClipType clipType,
-  Polygons &solution,
-  PolyFillType subjFillType = pftEvenOdd,
-  PolyFillType clipFillType = pftEvenOdd);
+    Polygons &solution,
+    PolyFillType subjFillType = pftEvenOdd,
+    PolyFillType clipFillType = pftEvenOdd);
   bool Execute(ClipType clipType,
-  ExPolygons &solution,
-  PolyFillType subjFillType = pftEvenOdd,
-  PolyFillType clipFillType = pftEvenOdd);
+    PolyTree &polytree,
+    PolyFillType subjFillType = pftEvenOdd,
+    PolyFillType clipFillType = pftEvenOdd);
   void Clear();
   bool ReverseSolution() {return m_ReverseOutput;};
   void ReverseSolution(bool value) {m_ReverseOutput = value;};
 protected:
   void Reset();
-  virtual bool ExecuteInternal(bool fixHoleLinkages);
+  virtual bool ExecuteInternal();
 private:
   PolyOutList       m_PolyOuts;
   JoinList          m_Joins;
@@ -227,11 +259,12 @@ private:
   Scanbeam         *m_Scanbeam;
   TEdge           *m_ActiveEdges;
   TEdge           *m_SortedEdges;
-  IntersectNode    *m_IntersectNodes;
-  bool              m_ExecuteLocked;
-  PolyFillType      m_ClipFillType;
-  PolyFillType      m_SubjFillType;
-  bool              m_ReverseOutput;
+  IntersectNode   *m_IntersectNodes;
+  bool             m_ExecuteLocked;
+  PolyFillType     m_ClipFillType;
+  PolyFillType     m_SubjFillType;
+  bool             m_ReverseOutput;
+  bool             m_UsingPolyTree; 
   void DisposeScanbeamList();
   void SetWindingCount(TEdge& edge);
   bool IsEvenOddFillType(const TEdge& edge) const;
@@ -259,10 +292,9 @@ private:
   void DoEdge2(TEdge *edge1, TEdge *edge2, const IntPoint &pt);
   void DoBothEdges(TEdge *edge1, TEdge *edge2, const IntPoint &pt);
   void IntersectEdges(TEdge *e1, TEdge *e2,
-    const IntPoint &pt, IntersectProtects protects);
+    const IntPoint &pt, const IntersectProtects protects);
   OutRec* CreateOutRec();
   void AddOutPt(TEdge *e, const IntPoint &pt);
-  void DisposeBottomPt(OutRec &outRec);
   void DisposeAllPolyPts();
   void DisposeOutRec(PolyOutList::size_type index);
   bool ProcessIntersections(const long64 botY, const long64 topY);
@@ -271,18 +303,22 @@ private:
   void ProcessIntersectList();
   void ProcessEdgesAtTopOfScanbeam(const long64 topY);
   void BuildResult(Polygons& polys);
-  void BuildResultEx(ExPolygons& polys);
+  void BuildResult2(PolyTree& polytree);
   void SetHoleState(TEdge *e, OutRec *OutRec);
   void DisposeIntersectNodes();
-  bool FixupIntersections();
+  bool FixupIntersectionOrder();
   void FixupOutPolygon(OutRec &outRec);
   bool IsHole(TEdge *e);
-  void FixHoleLinkage(OutRec *outRec);
+  void FixHoleLinkage(OutRec &outRec);
   void AddJoin(TEdge *e1, TEdge *e2, int e1OutIdx = -1, int e2OutIdx = -1);
   void ClearJoins();
   void AddHorzJoin(TEdge *e, int idx);
   void ClearHorzJoins();
-  void JoinCommonEdges(bool fixHoleLinkages);
+  bool JoinPoints(const JoinRec *j, OutPt *&p1, OutPt *&p2);
+  void FixupJoinRecs(JoinRec *j, OutPt *pt, unsigned startIdx);
+  void JoinCommonEdges();
+  void FixupFirstLefts1(OutRec* OldOutRec, OutRec* NewOutRec);
+  void FixupFirstLefts2(OutRec* OldOutRec, OutRec* NewOutRec);
 };
 
 //------------------------------------------------------------------------------
diff --git a/src/Lib/terragear/tg_accumulator.cxx b/src/Lib/terragear/tg_accumulator.cxx
index c8037516..11fbb0ba 100644
--- a/src/Lib/terragear/tg_accumulator.cxx
+++ b/src/Lib/terragear/tg_accumulator.cxx
@@ -61,8 +61,8 @@ void tgAccumulator::Add( const tgContour& subject )
 
 void tgAccumulator::ToShapefiles( const std::string& path, const std::string& layer_prefix, bool individual )
 {
-    char shapefile[16];
-    char layer[16];
+    char shapefile[32];
+    char layer[32];
 
     if ( individual ) {
         for (unsigned int i=0; i < accum.size(); i++) {
diff --git a/src/Lib/terragear/tg_misc.cxx b/src/Lib/terragear/tg_misc.cxx
index 566d8b30..cd6b8eaa 100644
--- a/src/Lib/terragear/tg_misc.cxx
+++ b/src/Lib/terragear/tg_misc.cxx
@@ -194,6 +194,7 @@ SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod&
     return OffsetPointMiddle( gPrev, gCur, gNext, offset_by, unused );
 }
 
+
 SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by )
 {
     double courseOffset;
@@ -232,6 +233,106 @@ SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by
     return pt;
 }
 
+
+void OffsetPointsMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, double width, int& turn_dir, SGGeod& inner, SGGeod& outer )
+{
+    double  courseCur, courseNext, courseAvg, theta;
+    SGVec3d dirCur, dirNext, dirAvg, cp;
+    double  courseOffset, distOffsetInner, distOffsetOuter;
+    SGGeod  pt;
+
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Find average angle for contour: prev (" << gPrev << "), "
+                                                                  "cur (" << gCur  << "), "
+                                                                 "next (" << gNext << ")" );
+
+    // first, find if the line turns left or right ar src
+    // for this, take the cross product of the vectors from prev to src, and src to next.
+    // if the cross product is negetive, we've turned to the left
+    // if the cross product is positive, we've turned to the right
+    courseCur = SGGeodesy::courseDeg( gCur, gPrev );
+    dirCur = SGVec3d( sin( courseCur*SGD_DEGREES_TO_RADIANS ), cos( courseCur*SGD_DEGREES_TO_RADIANS ), 0.0f );
+
+    courseNext = SGGeodesy::courseDeg( gCur, gNext );
+    dirNext = SGVec3d( sin( courseNext*SGD_DEGREES_TO_RADIANS ), cos( courseNext*SGD_DEGREES_TO_RADIANS ), 0.0f );
+
+    // Now find the average
+    dirAvg = normalize( dirCur + dirNext );
+    courseAvg = SGMiscd::rad2deg( atan( dirAvg.x()/dirAvg.y() ) );
+    if (courseAvg < 0) {
+        courseAvg += 180.0f;
+    }
+
+    // check the turn direction
+    cp    = cross( dirCur, dirNext );
+    theta = SGMiscd::rad2deg(CalculateTheta( dirCur, dirNext ) );
+
+    if ( (abs(theta - 180.0) < 0.1) || (abs(theta) < 0.1) || (isnan(theta)) ) {
+        // straight line blows up math - offset 90 degree and dist is as given
+        courseOffset = SGMiscd::normalizePeriodic(0, 360, courseNext-90.0);
+        distOffsetInner  = offset_by+width/2.0f;
+        distOffsetOuter  = offset_by-width/2.0f;
+    }  else  {
+        // calculate correct distance for the offset point
+        if (cp.z() < 0.0f) {
+            courseOffset = SGMiscd::normalizePeriodic(0, 360, courseAvg+180);
+            turn_dir = 0;
+        } else {
+            courseOffset = SGMiscd::normalizePeriodic(0, 360, courseAvg);
+            turn_dir = 1;
+        }
+        distOffsetInner = (offset_by+width/2.0f)/sin(SGMiscd::deg2rad(courseNext-courseOffset));
+        distOffsetOuter = (offset_by-width/2.0f)/sin(SGMiscd::deg2rad(courseNext-courseOffset));
+    }
+
+    // calculate the points from cur
+    inner = SGGeodesy::direct(gCur, courseOffset, distOffsetInner);
+    outer = SGGeodesy::direct(gCur, courseOffset, distOffsetOuter);
+}
+
+void OffsetPointsMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, double width, SGGeod& inner, SGGeod& outer )
+{
+    int unused;
+    return OffsetPointsMiddle( gPrev, gCur, gNext, offset_by, width, unused, inner, outer );
+}
+
+
+void OffsetPointsFirst( const SGGeod& cur, const SGGeod& next, double offset_by, double width, SGGeod& inner, SGGeod& outer )
+{
+    double courseOffset;
+
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at Start : cur (" << cur  << "), "
+                                                            "next (" << next << ")" );
+
+    // find the offset angle
+    courseOffset = SGGeodesy::courseDeg( cur, next ) - 90;
+    courseOffset = SGMiscd::normalizePeriodic(0, 360, courseOffset);
+
+    // calculate the point from cur
+    inner = SGGeodesy::direct( cur, courseOffset, offset_by+width/2.0f );
+    outer = SGGeodesy::direct( cur, courseOffset, offset_by-width/2.0f );
+}
+
+void OffsetPointsLast( const SGGeod& prev, const SGGeod& cur, double offset_by, double width, SGGeod& inner, SGGeod& outer )
+{
+    double courseOffset;
+
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at End   : prev (" << prev  << "), "
+                                                              "cur (" << cur << ")" );
+
+    // find the offset angle
+    courseOffset = SGGeodesy::courseDeg( prev, cur ) - 90;
+    courseOffset = SGMiscd::normalizePeriodic(0, 360, courseOffset);
+
+    // calculate the point from cur
+    inner = SGGeodesy::direct( cur, courseOffset, offset_by+width/2.0f );
+    outer = SGGeodesy::direct( cur, courseOffset, offset_by-width/2.0f );
+}
+
+
+
+
+
+
 SGGeod midpoint( const SGGeod& p0, const SGGeod& p1 )
 {
     return SGGeod::fromDegM( (p0.getLongitudeDeg() + p1.getLongitudeDeg()) / 2,
diff --git a/src/Lib/terragear/tg_polygon.hxx b/src/Lib/terragear/tg_polygon.hxx
index 1125f7f7..6a91f4df 100644
--- a/src/Lib/terragear/tg_polygon.hxx
+++ b/src/Lib/terragear/tg_polygon.hxx
@@ -71,6 +71,11 @@ SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod&
 SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by );
 SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by );
 
+void OffsetPointsMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, double width, int& turn_dir, SGGeod& inner, SGGeod& outer );
+void OffsetPointsMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, double width, SGGeod& inner, SGGeod& outer );
+void OffsetPointsFirst( const SGGeod& cur, const SGGeod& next, double offset_by, double width, SGGeod& inner, SGGeod& outer );
+void OffsetPointsLast( const SGGeod& prev, const SGGeod& cur, double offset_by, double width, SGGeod& inner, SGGeod& outer );
+
 // what abount this?
 
 // Save clipper to shapefile
diff --git a/src/Lib/terragear/tg_polygon_clean.cxx b/src/Lib/terragear/tg_polygon_clean.cxx
index 11785628..1855391b 100644
--- a/src/Lib/terragear/tg_polygon_clean.cxx
+++ b/src/Lib/terragear/tg_polygon_clean.cxx
@@ -174,7 +174,6 @@ tgPolygon tgPolygon::RemoveSpikes( const tgPolygon& subject )
 // Move slivers from in polygon to out polygon.
 void tgPolygon::RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers )
 {
-#if 0
     // traverse each contour of the polygon and attempt to identify
     // likely slivers
     SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygon::RemoveSlivers()");
@@ -226,7 +225,6 @@ void tgPolygon::RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers )
             }
         }
     }
-#endif
 }
 
 tgcontour_list tgPolygon::MergeSlivers( tgpolygon_list& polys, tgcontour_list& sliver_list ) {
diff --git a/src/Lib/terragear/tg_shapefile.cxx b/src/Lib/terragear/tg_shapefile.cxx
index 5b5d3dbc..6638db3c 100644
--- a/src/Lib/terragear/tg_shapefile.cxx
+++ b/src/Lib/terragear/tg_shapefile.cxx
@@ -128,6 +128,50 @@ void tgShapefile::FromClipper( const ClipperLib::Polygons& subject, const std::s
     }
 }
 
+void tgShapefile::FromContour( const tgContour& 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::OpenLayer( ds_id, layer.c_str() );
+    SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id);
+
+    OGRPolygon*    polygon = new OGRPolygon();
+
+    if (subject.GetSize() < 3) {
+        SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points");
+    } else {
+        // 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=new OGRLinearRing();
+        for (unsigned int pt = 0; pt < subject.GetSize(); pt++) {
+            OGRPoint *point=new OGRPoint();
+
+            point->setX( subject.GetNode(pt).getLongitudeDeg() );
+            point->setY( subject.GetNode(pt).getLatitudeDeg() );
+            point->setZ( 0.0 );
+            ring->addPoint(point);
+        }
+        ring->closeRings();
+
+        polygon->addRingDirectly(ring);
+
+        OGRFeature* feature = NULL;
+        feature = new OGRFeature( l_id->GetLayerDefn() );
+        feature->SetField("ID", description.c_str());
+        feature->SetGeometry(polygon);
+        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() );
diff --git a/src/Lib/terragear/tg_shapefile.hxx b/src/Lib/terragear/tg_shapefile.hxx
index f3917c50..8e059fab 100644
--- a/src/Lib/terragear/tg_shapefile.hxx
+++ b/src/Lib/terragear/tg_shapefile.hxx
@@ -6,6 +6,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 tgPolygon ToPolygon( const void* subject );