diff --git a/src/db/db/dbEdgeProcessor.cc b/src/db/db/dbEdgeProcessor.cc
index e4745add0..287ee022e 100644
--- a/src/db/db/dbEdgeProcessor.cc
+++ b/src/db/db/dbEdgeProcessor.cc
@@ -44,6 +44,8 @@
 namespace db
 {
 
+const double fill_factor = 1.5;
+
 // -------------------------------------------------------------------------------
 //  Some utilities ..
 
@@ -1095,14 +1097,14 @@ EdgeProcessor::clear ()
 }
 
 static void
-add_hparallel_cutpoints (WorkEdge &e1, WorkEdge &e2, std::vector <CutPoints> &cutpoints)
+add_hparallel_cutpoints (WorkEdge &e1, WorkEdge &e2, const db::Box &cell, std::vector <CutPoints> &cutpoints)
 {
   db::Coord e1_xmin = std::min (e1.x1 (), e1.x2 ());
   db::Coord e1_xmax = std::max (e1.x1 (), e1.x2 ());
-  if (e2.x1 () > e1_xmin && e2.x1 () < e1_xmax) {
+  if (e2.x1 () > e1_xmin && e2.x1 () < e1_xmax && cell.contains (e2.p1 ())) {
     e1.make_cutpoints (cutpoints)->add (e2.p1 (), &cutpoints, false);
   }
-  if (e2.x2 () > e1_xmin && e2.x2 () < e1_xmax) {
+  if (e2.x2 () > e1_xmin && e2.x2 () < e1_xmax && cell.contains (e2.p2 ())) {
     e1.make_cutpoints (cutpoints)->add (e2.p2 (), &cutpoints, false);
   }
 }
@@ -1125,14 +1127,28 @@ get_intersections_per_band_90 (std::vector <CutPoints> &cutpoints, std::vector <
   std::vector <WorkEdge>::iterator f = current;
   for (std::vector <WorkEdge>::iterator c = current; c != future; ) {
 
-    while (f != future && edge_xmin (*f) <= x) {
-      ++f;
-    }
+    size_t n = 0;
+    db::Coord xx = x;
 
-    db::Coord xx = std::numeric_limits <db::Coord>::max ();
-    if (f != future) {
-      xx = edge_xmin (*f);
-    }
+    //  fetch as many cells as to fill in roughly 50% more
+    //  (this is an empirical performance improvement factor)
+    do {
+
+      while (f != future && edge_xmin (*f) <= xx) {
+        ++f;
+      }
+
+      if (f != future) {
+        xx = edge_xmin (*f);
+      } else {
+        xx = std::numeric_limits <db::Coord>::max ();
+      }
+
+      if (n == 0) {
+        n = std::distance (c, f);
+      }
+
+    } while (f != future && std::distance (c, f) < long (n * fill_factor));
 
 #ifdef DEBUG_EDGE_PROCESSOR
     printf ("edges %d..%d:", x, xx); 
@@ -1165,15 +1181,15 @@ get_intersections_per_band_90 (std::vector <CutPoints> &cutpoints, std::vector <
 
                 //  parallel horizontal edges: produce the end points of each other edge as cutpoints
                 if (c1->p1 ().y () == c2->p1 ().y ()) {
-                  add_hparallel_cutpoints (*c1, *c2, cutpoints);
-                  add_hparallel_cutpoints (*c2, *c1, cutpoints);
+                  add_hparallel_cutpoints (*c1, *c2, cell, cutpoints);
+                  add_hparallel_cutpoints (*c2, *c1, cell, cutpoints);
                 }
 
               } else if (c1->p1 () != c2->p1 () && c1->p2 () != c2->p1 () &&
                          c1->p1 () != c2->p2 () && c1->p2 () != c2->p2 ()) {
 
                 std::pair <bool, db::Point> cp = c1->intersect_point (*c2);
-                if (cp.first) {
+                if (cp.first && cell.contains (cp.second)) {
 
                   //  add a cut point to c1 and c2 (c2 only if necessary)
                   c1->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, true);
@@ -1197,7 +1213,7 @@ get_intersections_per_band_90 (std::vector <CutPoints> &cutpoints, std::vector <
                            c1->p1 () != c2->p2 () && c1->p2 () != c2->p2 ()) {
 
               std::pair <bool, db::Point> cp = c1->intersect_point (*c2);
-              if (cp.first) {
+              if (cp.first && cell.contains (cp.second)) {
                 
                 //  add a cut point to c1 and c2
                 c2->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, true);
@@ -1336,10 +1352,9 @@ public:
 static void 
 get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector <WorkEdge>::iterator current, std::vector <WorkEdge>::iterator future, db::Coord y, db::Coord yy, bool with_h)
 {
-  std::vector <WorkEdge *> p1_weak; // holds weak interactions of edge endpoints with other edges
-  std::vector <WorkEdge *> ip_weak;
   double dy = y - 0.5;
   double dyy = yy + 0.5;
+  std::vector <std::pair<const WorkEdge *, WorkEdge *> > p1_weak;   // holds weak interactions of edge endpoints with other edges
 
   std::sort (current, future, edge_xmin_at_yinterval_double_compare<db::Coord> (dy, dyy));
 
@@ -1356,14 +1371,28 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
   std::vector <WorkEdge>::iterator f = current;
   for (std::vector <WorkEdge>::iterator c = current; c != future; ) {
 
-    while (f != future && edge_xmin_at_yinterval_double (*f, dy, dyy) <= x) {
-      ++f;
-    }
+    size_t n = 0;
+    db::Coord xx = x;
 
-    db::Coord xx = std::numeric_limits <db::Coord>::max ();
-    if (f != future) {
-      xx = edge_xmin_at_yinterval_double (*f, dy, dyy);
-    }
+    //  fetch as many cells as to fill in roughly 50% more
+    //  (this is an empirical performance improvement factor)
+    do {
+
+      while (f != future && edge_xmin_at_yinterval_double (*f, dy, dyy) <= xx) {
+        ++f;
+      }
+
+      if (f != future) {
+        xx = edge_xmin_at_yinterval_double (*f, dy, dyy);
+      } else {
+        xx = std::numeric_limits <db::Coord>::max ();
+      }
+
+      if (n == 0) {
+        n = std::distance (c, f);
+      }
+
+    } while (f != future && std::distance (c, f) < long (n * fill_factor));
 
 #ifdef DEBUG_EDGE_PROCESSOR
     printf ("edges %d..%d:", x, xx); 
@@ -1377,9 +1406,10 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
 
       db::Box cell (x, y, xx, yy);
 
-      for (std::vector <WorkEdge>::iterator c1 = c; c1 != f; ++c1) {
+      std::set<db::Point> weak_points;    // holds points that need to go in all other edges
+      p1_weak.clear ();
 
-        p1_weak.clear (); 
+      for (std::vector <WorkEdge>::iterator c1 = c; c1 != f; ++c1) {
 
         bool c1p1_in_cell = cell.contains (c1->p1 ());
         bool c1p2_in_cell = cell.contains (c1->p2 ());
@@ -1398,46 +1428,17 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
 
                 //  parallel horizontal edges: produce the end points of each other edge as cutpoints
                 if (c1->p1 ().y () == c2->p1 ().y ()) {
-                  add_hparallel_cutpoints (*c1, *c2, cutpoints);
-                  add_hparallel_cutpoints (*c2, *c1, cutpoints);
+                  add_hparallel_cutpoints (*c1, *c2, cell, cutpoints);
+                  add_hparallel_cutpoints (*c2, *c1, cell, cutpoints);
                 }
 
               } else if (c1->p1 () != c2->p1 () && c1->p2 () != c2->p1 () &&
                          c1->p1 () != c2->p2 () && c1->p2 () != c2->p2 ()) {
 
                 std::pair <bool, db::Point> cp = safe_intersect_point (*c1, *c2);
-                if (cp.first) {
-
-                  bool on_edge1 = is_point_on_exact (*c1, cp.second);
-
-                  //  add a cut point to c1 and c2 (points not on the edge give strong attractors)
-                  c1->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, !on_edge1);
-                  if (with_h) {
-                    c2->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, false);
-                  }
-
-#ifdef DEBUG_EDGE_PROCESSOR
-                  if (on_edge1) {
-                    printf ("weak intersection point %s between %s and %s.\n", cp.second.to_string ().c_str (), c1->to_string ().c_str (), c2->to_string ().c_str ());
-                  } else {
-                    printf ("intersection point %s between %s and %s.\n", cp.second.to_string ().c_str (), c1->to_string ().c_str (), c2->to_string ().c_str ());
-                  }
-#endif
-
-                  //  The new cutpoint must be inserted into other edges as well.
-                  ip_weak.clear ();
-                  for (std::vector <WorkEdge>::iterator cc = c; cc != f; ++cc) {
-                    if ((with_h || cc->dy () != 0) && cc != c1 && cc != c2 && is_point_on_fuzzy (*cc, cp.second)) {
-                      ip_weak.push_back (&*cc);
-                    }
-                  }
-                  for (std::vector <WorkEdge *>::iterator icc = ip_weak.begin (); icc != ip_weak.end (); ++icc) {
-                    (*icc)->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, true);
-#ifdef DEBUG_EDGE_PROCESSOR
-                    printf ("intersection point %s gives cutpoint in %s.\n", cp.second.to_string ().c_str (), (*icc)->to_string ().c_str ());
-#endif
-                  }
-
+                if (cp.first && cell.contains (cp.second)) {
+                  //  Stash the cutpoint as it must be inserted into other edges as well.
+                  weak_points.insert (cp.second);
                 }
 
               }
@@ -1478,41 +1479,9 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
                            c1->p1 () != c2->p2 () && c1->p2 () != c2->p2 ()) {
 
               std::pair <bool, db::Point> cp = safe_intersect_point (*c1, *c2);
-              if (cp.first) {
-                
-                bool on_edge1 = true;
-                bool on_edge2 = is_point_on_exact (*c2, cp.second);
-
-                //  add a cut point to c1 and c2
-                if (with_h || c1->dy () != 0) {
-                  on_edge1 = is_point_on_exact (*c1, cp.second);
-                  c1->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, !on_edge1);
-                }
-
-                c2->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, !on_edge2);
-
-#ifdef DEBUG_EDGE_PROCESSOR
-                if (!on_edge1 || !on_edge2) {
-                  printf ("intersection point %s between %s and %s.\n", cp.second.to_string ().c_str (), c1->to_string ().c_str (), c2->to_string ().c_str ()); 
-                } else {
-                  printf ("weak intersection point %s between %s and %s.\n", cp.second.to_string ().c_str (), c1->to_string ().c_str (), c2->to_string ().c_str ()); 
-                }
-#endif
-
-                //  The new cutpoint must be inserted into other edges as well.
-                ip_weak.clear ();
-                for (std::vector <WorkEdge>::iterator cc = c; cc != f; ++cc) {
-                  if ((with_h || cc->dy () != 0) && cc != c1 && cc != c2 && is_point_on_fuzzy (*cc, cp.second)) {
-                    ip_weak.push_back (&*cc);
-                  }
-                }
-                for (std::vector <WorkEdge *>::iterator icc = ip_weak.begin (); icc != ip_weak.end (); ++icc) {
-                  (*icc)->make_cutpoints (cutpoints)->add (cp.second, &cutpoints, true);
-#ifdef DEBUG_EDGE_PROCESSOR
-                  printf ("intersection point %s gives cutpoint in %s.\n", cp.second.to_string ().c_str (), (*icc)->to_string ().c_str ());
-#endif
-                }
-
+              if (cp.first && cell.contains (cp.second)) {
+                //  Stash the cutpoint as it must be inserted into other edges as well.
+                weak_points.insert (cp.second);
               }
 
             } 
@@ -1532,7 +1501,7 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
 #endif
                 c2->make_cutpoints (cutpoints)->add (c1->p1 (), &cutpoints, true);
               } else {
-                p1_weak.push_back (&*c2);
+                p1_weak.push_back (std::make_pair (c1.operator-> (), c2.operator-> ()));
               }
             }
 
@@ -1540,39 +1509,71 @@ get_intersections_per_band_any (std::vector <CutPoints> &cutpoints, std::vector
 
         }
 
-        if (! p1_weak.empty ()) {
+      }
 
-          bool strong = false;
-          for (std::vector<WorkEdge *>::const_iterator cp = p1_weak.begin (); cp != p1_weak.end () && ! strong; ++cp) {
-            if ((*cp)->data > 0 && cutpoints [(*cp)->data - 1].strong_cutpoints) {
-              strong = true;
-            }
-          }
+      //  insert weak intersection points into all relevant edges - weak into edges
+      //  where the point is on and strong into edges where the point is on in a fuzzy way.
 
-          p1_weak.back ()->make_cutpoints (cutpoints);
-          size_t n = p1_weak.back ()->data - 1;
-          for (std::vector<WorkEdge *>::const_iterator cp = p1_weak.begin (); cp != p1_weak.end (); ++cp) {
+      for (auto wp = weak_points.begin (); wp != weak_points.end (); ++wp) {
 
-            (*cp)->make_cutpoints (cutpoints);
-            size_t nn = (*cp)->data - 1;
-            if (strong) {
-              cutpoints [nn].add (c1->p1 (), &cutpoints);
+        for (std::vector <WorkEdge>::iterator cc = c; cc != f; ++cc) {
+          if ((with_h || cc->dy () != 0) && is_point_on_fuzzy (*cc, *wp)) {
+            bool on_edge = is_point_on_exact (*cc, *wp);
+            cc->make_cutpoints (cutpoints)->add (*wp, &cutpoints, !on_edge);
 #ifdef DEBUG_EDGE_PROCESSOR
-              printf ("Insert strong attractor %s in %s.\n", c1->p1 ().to_string ().c_str (), (*cp)->to_string ().c_str ()); 
-#endif
+            if (!on_edge) {
+              printf ("intersection point %s gives strong cutpoint in %s.\n", wp->to_string ().c_str (), cc->to_string ().c_str ());
             } else {
-              cutpoints [nn].add_attractor (c1->p1 (), n);
-#ifdef DEBUG_EDGE_PROCESSOR
-              printf ("Insert weak attractor %s in %s.\n", c1->p1 ().to_string ().c_str (), (*cp)->to_string ().c_str ()); 
-#endif
+              printf ("intersection point %s gives weak cutpoint in %s.\n", wp->to_string ().c_str (), cc->to_string ().c_str ());
             }
-
-            n = nn;
-
+#endif
           }
+        }
+
+      }
+
+      //  go through the list of "p1 to other edges" and insert p1 either as cutpoint
+      //  (if there are other strong cutpoints already) or weak attractor.
+
+      auto p1w_from = p1_weak.begin ();
+      while (p1w_from != p1_weak.end ()) {
+
+        bool strong = false;
+        auto p1w_to = p1w_from;
+        while (p1w_to != p1_weak.end () && p1w_to->first == p1w_from->first) {
+          if (p1w_to->second->data > 0 && cutpoints [p1w_to->second->data - 1].strong_cutpoints) {
+            strong = true;
+          }
+          ++p1w_to;
+        }
+
+        db::Point p1 = p1w_from->first->p1 ();
+
+        p1w_to [-1].second->make_cutpoints (cutpoints);
+        size_t n = p1w_to [-1].second->data - 1;
+
+        for (auto cp = p1w_from; cp != p1w_to; ++cp) {
+
+          cp->second->make_cutpoints (cutpoints);
+          size_t nn = cp->second->data - 1;
+          if (strong) {
+            cutpoints [nn].add (p1, &cutpoints);
+#ifdef DEBUG_EDGE_PROCESSOR
+            printf ("Insert strong attractor %s in %s.\n", cp->first->p1 ().to_string ().c_str (), cp->second->to_string ().c_str ());
+#endif
+          } else {
+            cutpoints [nn].add_attractor (p1, n);
+#ifdef DEBUG_EDGE_PROCESSOR
+            printf ("Insert weak attractor %s in %s.\n", cp->first->p1 ().to_string ().c_str (), cp->second->to_string ().c_str ());
+#endif
+          }
+
+          n = nn;
 
         }
 
+        p1w_from = p1w_to;
+
       }
 
     }
@@ -2259,7 +2260,7 @@ EdgeProcessor::redo_or_process (const std::vector<std::pair<db::EdgeSink *, db::
         progress->set (size_t (double (todo_next - todo) * p) + todo);
       }
 
-      size_t n = std::distance (current, future);
+      size_t n = 0;
       db::Coord yy = y;
 
       //  Use as many scanlines as to fetch approx. 50% new edges into the scanline (this
@@ -2276,7 +2277,11 @@ EdgeProcessor::redo_or_process (const std::vector<std::pair<db::EdgeSink *, db::
           yy = std::numeric_limits <db::Coord>::max ();
         }
 
-      } while (future != mp_work_edges->end () && std::distance (current, future) < long (n + n / 2));
+        if (n == 0) {
+          n = std::distance (current, future);
+        }
+
+      } while (future != mp_work_edges->end () && std::distance (current, future) < long (n * fill_factor));
 
       bool is90 = true;
 
diff --git a/src/db/db/dbPolygon.h b/src/db/db/dbPolygon.h
index 9eceb0a46..0b3a1f603 100644
--- a/src/db/db/dbPolygon.h
+++ b/src/db/db/dbPolygon.h
@@ -1044,6 +1044,19 @@ private:
   }
 };
 
+} // namespace db
+
+namespace std
+{
+
+//  injecting a global std::swap for polygons into the
+//  std namespace
+template <class C>
+void swap (db::polygon_contour<C> &a, db::polygon_contour<C> &b)
+{
+  a.swap (b);
+}
+
 }
 
 namespace db
@@ -2047,7 +2060,7 @@ public:
    *  @param remove_reflected True, if reflecting spikes shall be removed on compression
    */
   template <class I> 
-  void insert_hole (I start, I end, bool compress = default_compression<C> (), bool remove_reflected = false) 
+  void insert_hole (I start, I end, bool compress = default_compression<C> (), bool remove_reflected = false)
   {
     insert_hole (start, end, db::unit_trans<C> (), compress, remove_reflected);
   }
@@ -2073,21 +2086,17 @@ public:
     //  add the hole
     contour_type &h = add_hole ();
     h.assign (start, end, op, true, compress, true /*normalize*/, remove_reflected);
+  }
 
-    //  and keep the list sorted by swapping the 
-    //  elements and move the last one to the right 
-    //  position
-    //  KLUDGE: this is probably a performance bottleneck from applications
-    //  with polygons with may holes ..
-    if (holes () > 1) {
-      typename contour_list_type::iterator ins_pos = std::lower_bound (m_ctrs.begin () + 1, m_ctrs.end () - 1, h);
-      typename contour_list_type::iterator p = m_ctrs.end () - 1;
-      if (ins_pos != p) {
-        while (p != ins_pos) {
-          p->swap (p [-1]);
-          --p;
-        }
-      }
+  /**
+   *  @brief Sort the holes
+   *
+   *  Sorting the holes makes certain algorithms more effective.
+   */
+  void sort_holes ()
+  {
+    if (! m_ctrs.empty ()) {
+      std::sort (m_ctrs.begin () + 1, m_ctrs.end ());
     }
   }
 
@@ -2884,6 +2893,16 @@ public:
     tl_assert (false);
   }
 
+  /**
+   *  @brief A dummy implementation of "sort_holes" provided for template instantiation
+   *
+   *  Asserts, if begin called.
+   */
+  void sort_holes ()
+  {
+    tl_assert (false);
+  }
+
   /**
    *  @brief A dummy implementation of "hole" provided for template instantiation
    *
@@ -3530,7 +3549,7 @@ inline void mem_stat (MemStatistics *stat, MemStatistics::purpose_t purpose, int
 namespace std
 {
 
-//  injecting a global std::swap for polygons into the 
+//  injecting a global std::swap for polygons into the
 //  std namespace
 template <class C>
 void swap (db::polygon<C> &a, db::polygon<C> &b)
diff --git a/src/db/db/dbPolygonGenerators.cc b/src/db/db/dbPolygonGenerators.cc
index 8849e3c1f..e14934415 100644
--- a/src/db/db/dbPolygonGenerators.cc
+++ b/src/db/db/dbPolygonGenerators.cc
@@ -865,6 +865,8 @@ PolygonGenerator::produce_poly (const PGPolyContour &c)
 
       }
 
+      m_poly.sort_holes ();
+
     }
 
     mp_psink->put (m_poly);
diff --git a/src/db/db/dbPolygonTools.cc b/src/db/db/dbPolygonTools.cc
index db18b69d3..5107b59bb 100644
--- a/src/db/db/dbPolygonTools.cc
+++ b/src/db/db/dbPolygonTools.cc
@@ -456,6 +456,8 @@ static bool _cut_polygon_internal (const PolygonType &input, const Edge &line, c
         }
       }
 
+      hull->sort_holes ();
+
       right_of_line->put (*hull);
 
     }
@@ -859,6 +861,8 @@ smooth (const db::Polygon &polygon, db::Coord d, bool keep_hv)
       }
     }
 
+    new_poly.sort_holes ();
+
   }
 
   return new_poly;
@@ -1351,6 +1355,8 @@ do_extract_rad (const db::polygon<C> &polygon, double &rinner, double &router, u
 
     }
 
+    new_polygon->sort_holes ();
+
   } else {
 
     if (! do_extract_rad_from_contour (polygon.begin_hull (), polygon.end_hull (), rinner, router, n, (std::vector<db::point<C> > *) 0, false)) {
@@ -1542,6 +1548,8 @@ do_compute_rounded (const db::polygon<C> &polygon, double rinner, double router,
     new_poly.insert_hole (new_pts.begin (), new_pts.end (), false /*don't compress*/);
   }
 
+  new_poly.sort_holes ();
+
   return new_poly;
 }
 
@@ -2882,6 +2890,8 @@ snapped_polygon (const db::Polygon &poly, db::Coord gx, db::Coord gy, std::vecto
 
   }
 
+  pnew.sort_holes ();
+
   return pnew;
 }
 
@@ -2921,6 +2931,8 @@ scaled_and_snapped_polygon (const db::Polygon &poly, db::Coord gx, db::Coord mx,
 
   }
 
+  pnew.sort_holes ();
+
   return pnew;
 }
 
diff --git a/src/db/db/gsiDeclDbPolygon.cc b/src/db/db/gsiDeclDbPolygon.cc
index 4b3d1e7bd..eaa36e98a 100644
--- a/src/db/db/gsiDeclDbPolygon.cc
+++ b/src/db/db/gsiDeclDbPolygon.cc
@@ -926,6 +926,11 @@ struct polygon_defs
     return c->insert_hole (&pts[0], &pts[0] + sizeof (pts) / sizeof (pts[0]));
   }
 
+  static void sort_holes (C *c)
+  {
+    c->sort_holes ();
+  }
+
   static C *from_string (const char *s)
   {
     tl::Extractor ex (s);
@@ -1249,6 +1254,14 @@ struct polygon_defs
       "@brief Iterates over the points that make up the nth hole\n"
       "The hole number must be less than the number of holes (see \\holes)"
     ) +
+    method_ext ("sort_holes", &sort_holes,
+      "@brief Brings the holes in a specific order\n"
+      "This function is normalize the hole order so the comparison of two "
+      "polygons does not depend on the order the holes were inserted. "
+      "Polygons generated by KLayout's alorithms have their holes sorted.\n"
+      "\n"
+      "This method has been introduced in version 0.28.8."
+    ) +
     method_ext ("size", &size_xy, gsi::arg ("dx"), gsi::arg ("dy"), gsi::arg ("mode"),
       "@brief Sizes the polygon (biasing)\n"
       "\n"
diff --git a/src/db/unit_tests/dbEdgeProcessorTests.cc b/src/db/unit_tests/dbEdgeProcessorTests.cc
index 238b854c9..ccadcd57e 100644
--- a/src/db/unit_tests/dbEdgeProcessorTests.cc
+++ b/src/db/unit_tests/dbEdgeProcessorTests.cc
@@ -2781,3 +2781,40 @@ TEST(135b)
   EXPECT_EQ (run_test135b (_this, db::Trans (db::Trans::m90)), "(-78,25;-33,34;-36,33;-37,33)");
   EXPECT_EQ (run_test135b (_this, db::Trans (db::Trans::m135)), "(-26,-78;-35,-33;-33,-36;-33,-37)");
 }
+
+//  issue #1366
+TEST(136)
+{
+  db::Layout layout_1;
+  unsigned int l_l20000d0;
+
+  {
+    std::string fn (tl::testdata ());
+    fn += "/bool/";
+    fn += "issue_1366.oas";
+    tl::InputStream stream (fn);
+    db::Reader reader (stream);
+
+    db::LoadLayoutOptions options;
+    reader.read (layout_1, options);
+
+    l_l20000d0 = layout_1.get_layer (db::LayerProperties (20000, 0));
+  }
+
+  db::ShapeProcessor proc;
+
+  db::Layout lr;
+  lr.dbu (0.0001);
+  db::Cell *lr_top = &lr.cell (lr.add_cell ("TOP"));
+
+  unsigned int lr_l100d0 = lr.insert_layer (db::LayerProperties (100, 0));
+
+  proc.merge (layout_1, layout_1.cell (*layout_1.begin_top_down ()), l_l20000d0,
+              lr_top->shapes (lr_l100d0), false /*hierarchical*/, 0, true /*resolve holes*/, true /*min coherence*/);
+
+  std::string au_fn (tl::testdata ());
+  au_fn += "/bool/";
+  au_fn += "issue_1366_au.gds";
+
+  db::compare_layouts (_this, lr, au_fn);
+}
diff --git a/src/db/unit_tests/dbPolygonTests.cc b/src/db/unit_tests/dbPolygonTests.cc
index c76e9630f..000b50f02 100644
--- a/src/db/unit_tests/dbPolygonTests.cc
+++ b/src/db/unit_tests/dbPolygonTests.cc
@@ -128,6 +128,7 @@ TEST(1)
   db::Polygon pp;
   pp.insert_hole (c3.begin (), c3.end ());
   pp.insert_hole (c2.begin (), c2.end ());
+  pp.sort_holes ();
   pp.assign_hull (c1.begin (), c1.end ());
   EXPECT_EQ (pp.area (), 1000*100-2*380*80);
   EXPECT_EQ (pp.area2 (), 2*(1000*100-2*380*80));
diff --git a/testdata/bool/issue_1366.oas b/testdata/bool/issue_1366.oas
new file mode 100644
index 000000000..ee3f5741d
Binary files /dev/null and b/testdata/bool/issue_1366.oas differ
diff --git a/testdata/bool/issue_1366_au.gds b/testdata/bool/issue_1366_au.gds
new file mode 100644
index 000000000..fb136f5d7
Binary files /dev/null and b/testdata/bool/issue_1366_au.gds differ
diff --git a/testdata/bool/special2_au1.oas b/testdata/bool/special2_au1.oas
index c55543d9f..a5f502f42 100644
Binary files a/testdata/bool/special2_au1.oas and b/testdata/bool/special2_au1.oas differ
diff --git a/testdata/bool/special2_au1_tz.oas b/testdata/bool/special2_au1_tz.oas
index 86229ca34..3290873a0 100644
Binary files a/testdata/bool/special2_au1_tz.oas and b/testdata/bool/special2_au1_tz.oas differ
diff --git a/testdata/bool/special2_au2.oas b/testdata/bool/special2_au2.oas
index 9885b6e32..8b2e223d4 100644
Binary files a/testdata/bool/special2_au2.oas and b/testdata/bool/special2_au2.oas differ
diff --git a/testdata/bool/special2_au2_tz.oas b/testdata/bool/special2_au2_tz.oas
index 9ca99c96f..1e82c2993 100644
Binary files a/testdata/bool/special2_au2_tz.oas and b/testdata/bool/special2_au2_tz.oas differ
diff --git a/testdata/bool/special2_au3.oas b/testdata/bool/special2_au3.oas
index 7410b1c99..861bc779c 100644
Binary files a/testdata/bool/special2_au3.oas and b/testdata/bool/special2_au3.oas differ
diff --git a/testdata/bool/special2_au3_tz.oas b/testdata/bool/special2_au3_tz.oas
index b2d93d569..8af391f06 100644
Binary files a/testdata/bool/special2_au3_tz.oas and b/testdata/bool/special2_au3_tz.oas differ
diff --git a/testdata/bool/special2_au4.oas b/testdata/bool/special2_au4.oas
index 0e2d0b2d6..e8580a4f1 100644
Binary files a/testdata/bool/special2_au4.oas and b/testdata/bool/special2_au4.oas differ
diff --git a/testdata/bool/special2_au4_tz.oas b/testdata/bool/special2_au4_tz.oas
index 381af30ae..8f595465e 100644
Binary files a/testdata/bool/special2_au4_tz.oas and b/testdata/bool/special2_au4_tz.oas differ
diff --git a/testdata/bool/special2_au5.oas b/testdata/bool/special2_au5.oas
index 045eac670..50fa6a864 100644
Binary files a/testdata/bool/special2_au5.oas and b/testdata/bool/special2_au5.oas differ
diff --git a/testdata/bool/special2_au5_tz.oas b/testdata/bool/special2_au5_tz.oas
index 6081a28b4..a374b17aa 100644
Binary files a/testdata/bool/special2_au5_tz.oas and b/testdata/bool/special2_au5_tz.oas differ
diff --git a/testdata/python/dbPolygonTest.py b/testdata/python/dbPolygonTest.py
index 5ed0fc95a..feab569e6 100644
--- a/testdata/python/dbPolygonTest.py
+++ b/testdata/python/dbPolygonTest.py
@@ -122,6 +122,8 @@ class DBPolygonTests(unittest.TestCase):
     b.assign_hole(1, [ pya.DPoint( 15, 25 ), pya.DPoint( 25, 25 ), pya.DPoint( 25, 65 ) ])
     self.assertEqual( str(b), "(0,1;1,5;1,1/10,20;20,20;20,60)" )
     b.insert_hole( [ pya.DPoint( 1, 2 ), pya.DPoint( 2, 2 ), pya.DPoint( 2, 6 ) ] )
+    self.assertEqual( str(b), "(0,1;1,5;1,1/10,20;20,20;20,60/1,2;2,2;2,6)" )
+    b.sort_holes()
     self.assertEqual( str(b), "(0,1;1,5;1,1/1,2;2,2;2,6/10,20;20,20;20,60)" )
     b.assign_hole(0, [ pya.DPoint( 15, 25 ), pya.DPoint( 25, 25 ), pya.DPoint( 25, 65 ) ])
     self.assertEqual( str(b), "(0,1;1,5;1,1/15,25;25,25;25,65/10,20;20,20;20,60)" )
@@ -234,6 +236,8 @@ class DBPolygonTests(unittest.TestCase):
     b.assign_hole(1, [ pya.Point( 15, 25 ), pya.Point( 25, 25 ), pya.Point( 25, 65 ) ])
     self.assertEqual( str(b), "(0,1;1,5;1,1/10,20;20,20;20,60)" )
     b.insert_hole( [ pya.Point( 1, 2 ), pya.Point( 2, 2 ), pya.Point( 2, 6 ) ] )
+    self.assertEqual( str(b), "(0,1;1,5;1,1/10,20;20,20;20,60/1,2;2,2;2,6)" )
+    b.sort_holes()
     self.assertEqual( str(b), "(0,1;1,5;1,1/1,2;2,2;2,6/10,20;20,20;20,60)" )
     b.assign_hole(0, [ pya.Point( 15, 25 ), pya.Point( 25, 25 ), pya.Point( 25, 65 ) ])
     self.assertEqual( str(b), "(0,1;1,5;1,1/15,25;25,25;25,65/10,20;20,20;20,60)" )
diff --git a/testdata/ruby/dbPolygonTest.rb b/testdata/ruby/dbPolygonTest.rb
index 97cdafc76..5d6911b73 100644
--- a/testdata/ruby/dbPolygonTest.rb
+++ b/testdata/ruby/dbPolygonTest.rb
@@ -128,6 +128,8 @@ class DBPolygon_TestClass < TestBase
     b.assign_hole(1, [ RBA::DPoint::new( 15, 25 ), RBA::DPoint::new( 25, 25 ), RBA::DPoint::new( 25, 65 ) ])
     assert_equal( b.to_s, "(0,1;1,5;1,1/10,20;20,20;20,60)" )
     b.insert_hole( [ RBA::DPoint::new( 1, 2 ), RBA::DPoint::new( 2, 2 ), RBA::DPoint::new( 2, 6 ) ] )
+    assert_equal( b.to_s, "(0,1;1,5;1,1/10,20;20,20;20,60/1,2;2,2;2,6)" )
+    b.sort_holes
     assert_equal( b.to_s, "(0,1;1,5;1,1/1,2;2,2;2,6/10,20;20,20;20,60)" )
     b.assign_hole(0, [ RBA::DPoint::new( 15, 25 ), RBA::DPoint::new( 25, 25 ), RBA::DPoint::new( 25, 65 ) ])
     assert_equal( b.to_s, "(0,1;1,5;1,1/15,25;25,25;25,65/10,20;20,20;20,60)" )
@@ -251,6 +253,8 @@ class DBPolygon_TestClass < TestBase
     b.assign_hole(1, [ RBA::Point::new( 15, 25 ), RBA::Point::new( 25, 25 ), RBA::Point::new( 25, 65 ) ])
     assert_equal( b.to_s, "(0,1;1,5;1,1/10,20;20,20;20,60)" )
     b.insert_hole( [ RBA::Point::new( 1, 2 ), RBA::Point::new( 2, 2 ), RBA::Point::new( 2, 6 ) ] )
+    assert_equal( b.to_s, "(0,1;1,5;1,1/10,20;20,20;20,60/1,2;2,2;2,6)" )
+    b.sort_holes
     assert_equal( b.to_s, "(0,1;1,5;1,1/1,2;2,2;2,6/10,20;20,20;20,60)" )
     b.assign_hole(0, [ RBA::Point::new( 15, 25 ), RBA::Point::new( 25, 25 ), RBA::Point::new( 25, 65 ) ])
     assert_equal( b.to_s, "(0,1;1,5;1,1/15,25;25,25;25,65/10,20;20,20;20,60)" )