Merge pull request #8242 from prusa3d/pm_occlusion_guided_seamstm_find_osx_10.12_bug

Implementation of occlusion guided seam placement - balance between visibility and angle sharpness.
Improvements in seam drawing and alignment, via cubic splines.
Spline/Polynomial fitting functions.
Refactoring of Point structure - angle functions.
Triangular mesh subdivision function.
Various small additions/refactoring.
! Macro which replaces std::deque with boost version on Windows builds.

31 files changed, 2551 insertions, 1430 deletions

diff --git a/doc/seam_placement/corner_penalty_function.png b/doc/seam_placement/corner_penalty_function.png
new file mode 100644
index 000000000..bf445d671
--- /dev/null
+++ b/doc/seam_placement/corner_penalty_function.png
diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp
index 217166f8c..3ea18de8d 100644
--- a/src/libslic3r/AABBTreeIndirect.hpp
+++ b/src/libslic3r/AABBTreeIndirect.hpp
@@ -709,10 +709,15 @@ inline bool intersect_ray_all_hits(
         origin, dir, VectorType(dir.cwiseInverse()),
         eps }
 	};
-	if (! tree.empty()) {
+	if (tree.empty()) {
+		hits.clear();
+	} else {
+		// Reusing the output memory if there is some memory already pre-allocated.
+        ray_intersector.hits = std::move(hits);
+        ray_intersector.hits.clear();
         ray_intersector.hits.reserve(8);
 		detail::intersect_ray_recursive_all_hits(ray_intersector, 0);
-		std::swap(hits, ray_intersector.hits);
+		hits = std::move(ray_intersector.hits);
 	    std::sort(hits.begin(), hits.end(), [](const auto &l, const auto &r) { return l.t < r.t; });
 	}
 	return ! hits.empty();
@@ -759,8 +764,8 @@ inline bool is_any_triangle_in_radius(
         const TreeType 						&tree,
         // Point to which the closest point on the indexed triangle set is searched for.
         const VectorType					&point,
-        // Maximum distance in which triangle is search for
-        typename VectorType::Scalar &max_distance)
+        //Square of maximum distance in which triangle is searched for
+        typename VectorType::Scalar &max_distance_squared)
 {
     using Scalar = typename VectorType::Scalar;
     auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
@@ -774,7 +779,7 @@ inline bool is_any_triangle_in_radius(
 		return false;
 	}
 
-	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance, hit_idx, hit_point);
+	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
 
     return hit_point.allFinite();
 }
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 838af44fb..a9ce8dd1f 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -138,10 +138,12 @@ set(SLIC3R_SOURCES
     GCodeWriter.hpp
     Geometry.cpp
     Geometry.hpp
+    Geometry/Bicubic.hpp
     Geometry/Circle.cpp
     Geometry/Circle.hpp
     Geometry/ConvexHull.cpp
     Geometry/ConvexHull.hpp
+    Geometry/Curves.hpp
     Geometry/MedialAxis.cpp
     Geometry/MedialAxis.hpp
     Geometry/Voronoi.hpp
@@ -225,6 +227,8 @@ set(SLIC3R_SOURCES
     SlicesToTriangleMesh.cpp
     SlicingAdaptive.cpp
     SlicingAdaptive.hpp
+    Subdivide.cpp
+    Subdivide.hpp
     SupportMaterial.cpp
     SupportMaterial.hpp
     Surface.cpp
diff --git a/src/libslic3r/Color.hpp b/src/libslic3r/Color.hpp
index 183705c4a..8044a0318 100644
--- a/src/libslic3r/Color.hpp
+++ b/src/libslic3r/Color.hpp
@@ -4,6 +4,8 @@
 #include <array>
 #include <algorithm>
 
+#include "Point.hpp"
+
 namespace Slic3r {
 
 class ColorRGB
@@ -133,41 +135,57 @@ public:
 	static const ColorRGBA Z()           { return { 0.0f, 0.0f, 0.75f, 1.0f }; }
 };
 
-extern ColorRGB operator * (float value, const ColorRGB& other);
-extern ColorRGBA operator * (float value, const ColorRGBA& other);
+ColorRGB operator * (float value, const ColorRGB& other);
+ColorRGBA operator * (float value, const ColorRGBA& other);
+
+ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t);
+ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t);
 
-extern ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t);
-extern ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t);
+ColorRGB complementary(const ColorRGB& color);
+ColorRGBA complementary(const ColorRGBA& color);
 
-extern ColorRGB complementary(const ColorRGB& color);
-extern ColorRGBA complementary(const ColorRGBA& color);
+ColorRGB saturate(const ColorRGB& color, float factor);
+ColorRGBA saturate(const ColorRGBA& color, float factor);
 
-extern ColorRGB saturate(const ColorRGB& color, float factor);
-extern ColorRGBA saturate(const ColorRGBA& color, float factor);
+ColorRGB opposite(const ColorRGB& color);
+ColorRGB opposite(const ColorRGB& a, const ColorRGB& b);
 
-extern ColorRGB opposite(const ColorRGB& color);
-extern ColorRGB opposite(const ColorRGB& a, const ColorRGB& b);
+bool can_decode_color(const std::string& color);
 
-extern bool can_decode_color(const std::string& color);
+bool decode_color(const std::string& color_in, ColorRGB& color_out);
+bool decode_color(const std::string& color_in, ColorRGBA& color_out);
 
-extern bool decode_color(const std::string& color_in, ColorRGB& color_out);
-extern bool decode_color(const std::string& color_in, ColorRGBA& color_out);
+bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGB>& colors_out);
+bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGBA>& colors_out);
 
-extern bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGB>& colors_out);
-extern bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGBA>& colors_out);
+std::string encode_color(const ColorRGB& color);
+std::string encode_color(const ColorRGBA& color);
 
-extern std::string encode_color(const ColorRGB& color);
-extern std::string encode_color(const ColorRGBA& color);
+ColorRGB  to_rgb(const ColorRGBA& other_rgba);
+ColorRGBA to_rgba(const ColorRGB& other_rgb);
+ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha);
 
-extern ColorRGB  to_rgb(const ColorRGBA& other_rgba);
-extern ColorRGBA to_rgba(const ColorRGB& other_rgb);
-extern ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha);
+// Color mapping of a value into RGB false colors.
+inline Vec3f value_to_rgbf(float minimum, float maximum, float value) 
+{
+    float ratio = 2.0f * (value - minimum) / (maximum - minimum);
+    float b = std::max(0.0f, (1.0f - ratio));
+    float r = std::max(0.0f, (ratio - 1.0f));
+    float g = 1.0f - b - r;
+    return Vec3f { r, g, b };
+}
+
+// Color mapping of a value into RGB false colors.
+inline Vec3i value_to_rgbi(float minimum, float maximum, float value)
+{
+    return (value_to_rgbf(minimum, maximum, value) * 255).cast<int>();
+}
 
-extern ColorRGBA picking_decode(unsigned int id);
-extern unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b);
+ColorRGBA picking_decode(unsigned int id);
+unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b);
 // Produce an alpha channel checksum for the red green blue components. The alpha channel may then be used to verify, whether the rgb components
 // were not interpolated by alpha blending or multi sampling.
-extern unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue);
+unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue);
 
 } // namespace Slic3r
 
diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp
index 3d0ca5069..1a6ee4b80 100644
--- a/src/libslic3r/GCode.cpp
+++ b/src/libslic3r/GCode.cpp
@@ -2312,7 +2312,6 @@ GCode::LayerResult GCode::process_layer(
     } // for objects
 
     // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders.
-    std::vector<std::unique_ptr<EdgeGrid::Grid>> lower_layer_edge_grids(layers.size());
     for (unsigned int extruder_id : layer_tools.extruders)
     {
         gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ?
@@ -2337,7 +2336,7 @@ GCode::LayerResult GCode::process_layer(
                     path.mm3_per_mm = mm3_per_mm;
                 }
                 //FIXME using the support_material_speed of the 1st object printed.
-                gcode += this->extrude_loop(loop, "skirt", m_config.support_material_speed.value);
+                gcode += this->extrude_loop(loop, "skirt"sv, m_config.support_material_speed.value);
             }
             m_avoid_crossing_perimeters.use_external_mp(false);
             // Allow a straight travel move to the first object point if this is the first layer (but don't in next layers).
@@ -2350,7 +2349,7 @@ GCode::LayerResult GCode::process_layer(
             this->set_origin(0., 0.);
             m_avoid_crossing_perimeters.use_external_mp();
             for (const ExtrusionEntity *ee : print.brim().entities) {
-                gcode += this->extrude_entity(*ee, "brim", m_config.support_material_speed.value);
+                gcode += this->extrude_entity(*ee, "brim"sv, m_config.support_material_speed.value);
             }
             m_brim_done = true;
             m_avoid_crossing_perimeters.use_external_mp(false);
@@ -2406,9 +2405,9 @@ GCode::LayerResult GCode::process_layer(
                     //FIXME the following code prints regions in the order they are defined, the path is not optimized in any way.
                     if (print.config().infill_first) {
                         gcode += this->extrude_infill(print, by_region_specific, false);
-                        gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]);
+                        gcode += this->extrude_perimeters(print, by_region_specific);
                     } else {
-                        gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]);
+                        gcode += this->extrude_perimeters(print, by_region_specific);
                         gcode += this->extrude_infill(print,by_region_specific, false);
                     }
                     // ironing
@@ -2543,51 +2542,28 @@ std::string GCode::change_layer(coordf_t print_z)
     return gcode;
 }
 
-
-
-static std::unique_ptr<EdgeGrid::Grid> calculate_layer_edge_grid(const Layer& layer)
-{
-    auto out = make_unique<EdgeGrid::Grid>();
-
-    // Create the distance field for a layer below.
-    const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5);
-    out->create(layer.lslices, distance_field_resolution);
-    out->calculate_sdf();
-#if 0
-        {
-            static int iRun = 0;
-            BoundingBox bbox = (*lower_layer_edge_grid)->bbox();
-            bbox.min(0) -= scale_(5.f);
-            bbox.min(1) -= scale_(5.f);
-            bbox.max(0) += scale_(5.f);
-            bbox.max(1) += scale_(5.f);
-            EdgeGrid::save_png(*(*lower_layer_edge_grid), bbox, scale_(0.1f), debug_out_path("GCode_extrude_loop_edge_grid-%d.png", iRun++));
-        }
-#endif
-    return out;
+static const auto comment_perimeter = "perimeter"sv;
+// Comparing string_view pointer & length for speed.
+static inline bool comment_is_perimeter(const std::string_view comment) {
+    return comment.data() == comment_perimeter.data() && comment.size() == comment_perimeter.size();
 }
 
-
-std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid)
+std::string GCode::extrude_loop(ExtrusionLoop loop, const std::string_view description, double speed)
 {
     // get a copy; don't modify the orientation of the original loop object otherwise
     // next copies (if any) would not detect the correct orientation
 
-    if (m_layer->lower_layer && lower_layer_edge_grid != nullptr && ! *lower_layer_edge_grid)
-        *lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
-
     // extrude all loops ccw
     bool was_clockwise = loop.make_counter_clockwise();
 
     // find the point of the loop that is closest to the current extruder position
     // or randomize if requested
     Point last_pos = this->last_pos();
-    if (m_config.spiral_vase) {
+    if (! m_config.spiral_vase && comment_is_perimeter(description)) {
+        assert(m_layer != nullptr);
+        m_seam_placer.place_seam(m_layer, loop, m_config.external_perimeters_first, this->last_pos());
+    } else
         loop.split_at(last_pos, false);
-    }
-    else
-        m_seam_placer.place_seam(loop, this->last_pos(), m_config.external_perimeters_first,
-                                 EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr);
 
     // clip the path to avoid the extruder to get exactly on the first point of the loop;
     // if polyline was shorter than the clipping distance we'd get a null polyline, so
@@ -2607,11 +2583,9 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
 
     // extrude along the path
     std::string gcode;
-    for (ExtrusionPaths::iterator path = paths.begin(); path != paths.end(); ++path) {
-//    description += ExtrusionLoop::role_to_string(loop.loop_role());
-//    description += ExtrusionEntity::role_to_string(path->role);
-        path->simplify(m_scaled_resolution);
-        gcode += this->_extrude(*path, description, speed);
+    for (ExtrusionPath &path : paths) {
+        path.simplify(m_scaled_resolution);
+        gcode += this->_extrude(path, description, speed);
     }
 
     // reset acceleration
@@ -2626,18 +2600,19 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
         // the side depends on the original winding order of the polygon (left for contours, right for holes)
         //FIXME improve the algorithm in case the loop is tiny.
         //FIXME improve the algorithm in case the loop is split into segments with a low number of points (see the Point b query).
-        Point a = paths.front().polyline.points[1];  // second point
-        Point b = *(paths.back().polyline.points.end()-3);       // second to last point
+        // Angle from the 2nd point to the last point.
+        double angle_inside = angle(paths.front().polyline.points[1]        - paths.front().first_point(),
+                                    *(paths.back().polyline.points.end()-3) - paths.front().first_point());
+        assert(angle_inside >= -M_PI && angle_inside <= M_PI);
+        // 3rd of this angle will be taken, thus make the angle monotonic before interpolation.
         if (was_clockwise) {
-            // swap points
-            Point c = a; a = b; b = c;
+            if (angle_inside > 0)
+                angle_inside -= 2.0 * M_PI;
+        } else {
+            if (angle_inside < 0)
+                angle_inside += 2.0 * M_PI;
         }
 
-        double angle = paths.front().first_point().ccw_angle(a, b) / 3;
-
-        // turn left if contour, turn right if hole
-        if (was_clockwise) angle *= -1;
-
         // create the destination point along the first segment and rotate it
         // we make sure we don't exceed the segment length because we don't know
         // the rotation of the second segment so we might cross the object boundary
@@ -2649,7 +2624,8 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
         // Shift by no more than a nozzle diameter.
         //FIXME Hiding the seams will not work nicely for very densely discretized contours!
         Point  pt = ((nd * nd >= l2) ? p2 : (p1 + v * (nd / sqrt(l2)))).cast<coord_t>();
-        pt.rotate(angle, paths.front().polyline.points.front());
+        // Rotate pt inside around the seam point.
+        pt.rotate(angle_inside / 3., paths.front().polyline.points.front());
         // generate the travel move
         gcode += m_writer.travel_to_xy(this->point_to_gcode(pt), "move inwards before travel");
     }
@@ -2657,13 +2633,11 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
     return gcode;
 }
 
-std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string description, double speed)
+std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed)
 {
     // extrude along the path
     std::string gcode;
     for (ExtrusionPath path : multipath.paths) {
-//    description += ExtrusionLoop::role_to_string(loop.loop_role());
-//    description += ExtrusionEntity::role_to_string(path->role);
         path.simplify(m_scaled_resolution);
         gcode += this->_extrude(path, description, speed);
     }
@@ -2676,22 +2650,21 @@ std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string
     return gcode;
 }
 
-std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid)
+std::string GCode::extrude_entity(const ExtrusionEntity &entity, const std::string_view description, double speed)
 {
     if (const ExtrusionPath* path = dynamic_cast<const ExtrusionPath*>(&entity))
         return this->extrude_path(*path, description, speed);
     else if (const ExtrusionMultiPath* multipath = dynamic_cast<const ExtrusionMultiPath*>(&entity))
         return this->extrude_multi_path(*multipath, description, speed);
     else if (const ExtrusionLoop* loop = dynamic_cast<const ExtrusionLoop*>(&entity))
-        return this->extrude_loop(*loop, description, speed, lower_layer_edge_grid);
+        return this->extrude_loop(*loop, description, speed);
     else
         throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()");
     return "";
 }
 
-std::string GCode::extrude_path(ExtrusionPath path, std::string description, double speed)
+std::string GCode::extrude_path(ExtrusionPath path, std::string_view description, double speed)
 {
-//    description += ExtrusionEntity::role_to_string(path.role());
     path.simplify(m_scaled_resolution);
     std::string gcode = this->_extrude(path, description, speed);
     if (m_wipe.enable) {
@@ -2704,24 +2677,15 @@ std::string GCode::extrude_path(ExtrusionPath path, std::string description, dou
 }
 
 // Extrude perimeters: Decide where to put seams (hide or align seams).
-std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid)
+std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region)
 {
     std::string gcode;
     for (const ObjectByExtruder::Island::Region &region : by_region)
         if (! region.perimeters.empty()) {
             m_config.apply(print.get_print_region(&region - &by_region.front()).config());
 
-            // plan_perimeters tries to place seams, it needs to have the lower_layer_edge_grid calculated already.
-            if (m_layer->lower_layer && ! lower_layer_edge_grid)
-                lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
-
-            m_seam_placer.plan_perimeters(std::vector<const ExtrusionEntity*>(region.perimeters.begin(), region.perimeters.end()),
-                *m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter),
-                (m_layer == NULL ? nullptr : m_layer->object()),
-                (lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr));
-
             for (const ExtrusionEntity* ee : region.perimeters)
-                gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid);
+                gcode += this->extrude_entity(*ee, comment_perimeter, -1.);
         }
     return gcode;
 }
@@ -2731,7 +2695,7 @@ std::string GCode::extrude_infill(const Print &print, const std::vector<ObjectBy
 {
     std::string 		 gcode;
     ExtrusionEntitiesPtr extrusions;
-    const char*          extrusion_name = ironing ? "ironing" : "infill";
+    const auto           extrusion_name = ironing ? "ironing"sv : "infill"sv;
     for (const ObjectByExtruder::Island::Region &region : by_region)
         if (! region.infills.empty()) {
             extrusions.clear();
@@ -2757,8 +2721,8 @@ std::string GCode::extrude_infill(const Print &print, const std::vector<ObjectBy
 
 std::string GCode::extrude_support(const ExtrusionEntityCollection &support_fills)
 {
-    static constexpr const char *support_label            = "support material";
-    static constexpr const char *support_interface_label  = "support material interface";
+    static constexpr const auto support_label            = "support material"sv;
+    static constexpr const auto support_interface_label  = "support material interface"sv;
 
     std::string gcode;
     if (! support_fills.entities.empty()) {
@@ -2767,7 +2731,7 @@ std::string GCode::extrude_support(const ExtrusionEntityCollection &support_fill
         for (const ExtrusionEntity *ee : support_fills.entities) {
             ExtrusionRole role = ee->role();
             assert(role == erSupportMaterial || role == erSupportMaterialInterface);
-            const char  *label = (role == erSupportMaterial) ? support_label : support_interface_label;
+            const auto   label = (role == erSupportMaterial) ? support_label : support_interface_label;
             const double speed = (role == erSupportMaterial) ? support_speed : support_interface_speed;
             const ExtrusionPath *path = dynamic_cast<const ExtrusionPath*>(ee);
             if (path)
@@ -2855,20 +2819,18 @@ void GCode::GCodeOutputStream::write_format(const char* format, ...)
     va_end(args);
 }
 
-std::string GCode::_extrude(const ExtrusionPath &path, std::string description, double speed)
+std::string GCode::_extrude(const ExtrusionPath &path, const std::string_view description, double speed)
 {
     std::string gcode;
-
-    if (is_bridge(path.role()))
-        description += " (bridge)";
+    const std::string_view description_bridge = is_bridge(path.role()) ? " (bridge)"sv : ""sv;
 
     // go to first point of extrusion path
     if (!m_last_pos_defined || m_last_pos != path.first_point()) {
-        gcode += this->travel_to(
-            path.first_point(),
-            path.role(),
-            "move to first " + description + " point"
-        );
+        std::string comment = "move to first ";
+        comment += description;
+        comment += description_bridge;
+        comment += " point";
+        gcode += this->travel_to(path.first_point(), path.role(), comment);
     }
 
     // compensate retraction
@@ -3006,7 +2968,11 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
     gcode += m_writer.set_speed(F, "", comment);
     double path_length = 0.;
     {
-        std::string comment = m_config.gcode_comments ? description : "";
+        std::string comment;
+        if (m_config.gcode_comments) {
+            comment = description;
+            comment += description_bridge;
+        }
         Vec2d prev = this->point_to_gcode_quantized(path.polyline.points.front());
         auto  it   = path.polyline.points.begin();
         auto  end  = path.polyline.points.end();
diff --git a/src/libslic3r/GCode.hpp b/src/libslic3r/GCode.hpp
index aa4682e34..5afe3d915 100644
--- a/src/libslic3r/GCode.hpp
+++ b/src/libslic3r/GCode.hpp
@@ -274,10 +274,10 @@ private:
     void            set_extruders(const std::vector<unsigned int> &extruder_ids);
     std::string     preamble();
     std::string     change_layer(coordf_t print_z);
-    std::string     extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr);
-    std::string     extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr);
-    std::string     extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.);
-    std::string     extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.);
+    std::string     extrude_entity(const ExtrusionEntity &entity, const std::string_view description, double speed = -1.);
+    std::string     extrude_loop(ExtrusionLoop loop, const std::string_view description, double speed = -1.);
+    std::string     extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed = -1.);
+    std::string     extrude_path(ExtrusionPath path, const std::string_view description, double speed = -1.);
 
     // Extruding multiple objects with soluble / non-soluble / combined supports
     // on a multi-material printer, trying to minimize tool switches.
@@ -342,7 +342,7 @@ private:
 		// For sequential print, the instance of the object to be printing has to be defined.
 		const size_t                     				 single_object_instance_idx);
 
-    std::string     extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid);
+    std::string     extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region);
     std::string     extrude_infill(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, bool ironing);
     std::string     extrude_support(const ExtrusionEntityCollection &support_fills);
 
@@ -428,7 +428,7 @@ private:
     // Processor
     GCodeProcessor m_processor;
 
-    std::string _extrude(const ExtrusionPath &path, std::string description = "", double speed = -1);
+    std::string _extrude(const ExtrusionPath &path, const std::string_view description, double speed = -1);
     void print_machine_envelope(GCodeOutputStream &file, Print &print);
     void _print_first_layer_bed_temperature(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait);
     void _print_first_layer_extruder_temperatures(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait);
diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp
index bcd238a72..6e0968145 100644
--- a/src/libslic3r/GCode/SeamPlacer.cpp
+++ b/src/libslic3r/GCode/SeamPlacer.cpp
@@ -1,1011 +1,1416 @@
 #include "SeamPlacer.hpp"
 
+#include "tbb/parallel_for.h"
+#include "tbb/blocked_range.h"
+#include "tbb/parallel_reduce.h"
+#include <boost/log/trivial.hpp>
+#include <random>
+#include <algorithm>
+#include <queue>
+
 #include "libslic3r/ExtrusionEntity.hpp"
 #include "libslic3r/Print.hpp"
 #include "libslic3r/BoundingBox.hpp"
+#include "libslic3r/Color.hpp"
 #include "libslic3r/EdgeGrid.hpp"
 #include "libslic3r/ClipperUtils.hpp"
-#include "libslic3r/SVG.hpp"
 #include "libslic3r/Layer.hpp"
+#include "libslic3r/QuadricEdgeCollapse.hpp"
+#include "libslic3r/Subdivide.hpp"
+
+#include "libslic3r/Geometry/Curves.hpp"
+
+#include "libslic3r/Utils.hpp"
+
+//#define DEBUG_FILES
+
+#ifdef DEBUG_FILES
+#include <boost/nowide/cstdio.hpp>
+#include <SVG.hpp>
+#endif
 
 namespace Slic3r {
 
-// This penalty is added to all points inside custom blockers (subtracted from pts inside enforcers).
-static constexpr float ENFORCER_BLOCKER_PENALTY = 100;
-
-// In case there are custom enforcers/blockers, the loop polygon shall always have
-// sides smaller than this (so it isn't limited to original resolution).
-static constexpr float MINIMAL_POLYGON_SIDE = scaled<float>(0.2f);
-
-// When spAligned is active and there is a support enforcer,
-// add this penalty to its center.
-static constexpr float ENFORCER_CENTER_PENALTY = -10.f;
-
-
-
-// This function was introduced in 2016 to assign penalties to overhangs.
-// LukasM thinks that it discriminated a bit too much, so especially external
-// seams were than placed in funny places (non-overhangs were preferred too much).
-// He implemented his own version (below) which applies fixed penalty for really big overlaps.
-// static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance)
-// {
-// // The extrudate is not fully supported by the lower layer. Fit a polynomial penalty curve.
-// // Solved by sympy package:
-// /*
-// from sympy import *
-// (x,a,b,c,d,r,z)=symbols('x a b c d r z')
-// p = a + b*x + c*x*x + d*x*x*x
-// p2 = p.subs(solve([p.subs(x, -r), p.diff(x).subs(x, -r), p.diff(x,x).subs(x, -r), p.subs(x, 0)-z], [a, b, c, d]))
-// from sympy.plotting import plot
-// plot(p2.subs(r,0.2).subs(z,1.), (x, -1, 3), adaptive=False, nb_of_points=400)
-// */
-//     if (overlap_distance < - nozzle_r) {
-//         // The extrudate is fully supported by the lower layer. This is the ideal case, therefore zero penalty.
-//         return 0.f;
-//     } else {
-//         float x  = overlap_distance / nozzle_r;
-//         float x2 = x * x;
-//         float x3 = x2 * x;
-//         return weight_zero * (1.f + 3.f * x + 3.f * x2 + x3);
-//     }
-// }
-static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance)
-{
-    return overlap_distance > nozzle_r ? weight_zero : 0.f;
+namespace SeamPlacerImpl {
+
+template<typename T> int sgn(T val) {
+    return int(T(0) < val) - int(val < T(0));
 }
 
+// base function: ((e^(((1)/(x^(2)+1)))-1)/(e-1))
+// checkout e.g. here: https://www.geogebra.org/calculator
+float gauss(float value, float mean_x_coord, float mean_value, float falloff_speed) {
+    float shifted = value - mean_x_coord;
+    float denominator = falloff_speed * shifted * shifted + 1.0f;
+    float exponent = 1.0f / denominator;
+    return mean_value * (std::exp(exponent) - 1.0f) / (std::exp(1.0f) - 1.0f);
+}
 
+/// Coordinate frame
+class Frame {
+public:
+    Frame() {
+        mX = Vec3f(1, 0, 0);
+        mY = Vec3f(0, 1, 0);
+        mZ = Vec3f(0, 0, 1);
+    }
 
-// Return a value in <0, 1> of a cubic B-spline kernel centered around zero.
-// The B-spline is re-scaled so it has value 1 at zero.
-static inline float bspline_kernel(float x)
-{
-    x = std::abs(x);
-    if (x < 1.f) {
-        return 1.f - (3.f / 2.f) * x * x + (3.f / 4.f) * x * x * x;
+    Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) :
+            mX(x), mY(y), mZ(z) {
     }
-    else if (x < 2.f) {
-        x -= 1.f;
-        float x2 = x * x;
-        float x3 = x2 * x;
-        return (1.f / 4.f) - (3.f / 4.f) * x + (3.f / 4.f) * x2 - (1.f / 4.f) * x3;
+
+    void set_from_z(const Vec3f &z) {
+        mZ = z.normalized();
+        Vec3f tmpZ = mZ;
+        Vec3f tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3f(0, 1, 0) : Vec3f(1, 0, 0);
+        mY = (tmpZ.cross(tmpX)).normalized();
+        mX = mY.cross(tmpZ);
     }
-    else
-        return 0;
-}
 
+    Vec3f to_world(const Vec3f &a) const {
+        return a.x() * mX + a.y() * mY + a.z() * mZ;
+    }
 
+    Vec3f to_local(const Vec3f &a) const {
+        return Vec3f(mX.dot(a), mY.dot(a), mZ.dot(a));
+    }
 
-static Points::const_iterator project_point_to_polygon_and_insert(Polygon &polygon, const Point &pt, double eps)
-{
-    assert(polygon.points.size() >= 2);
-    if (polygon.points.size() <= 1)
-    if (polygon.points.size() == 1)
-        return polygon.points.begin();
-
-    Point  pt_min;
-    double d_min = std::numeric_limits<double>::max();
-    size_t i_min = size_t(-1);
-
-    for (size_t i = 0; i < polygon.points.size(); ++ i) {
-        size_t j = i + 1;
-        if (j == polygon.points.size())
-            j = 0;
-        const Point &p1 = polygon.points[i];
-        const Point &p2 = polygon.points[j];
-        const Slic3r::Point v_seg = p2 - p1;
-        const Slic3r::Point v_pt  = pt - p1;
-        const int64_t l2_seg = int64_t(v_seg(0)) * int64_t(v_seg(0)) + int64_t(v_seg(1)) * int64_t(v_seg(1));
-        int64_t t_pt = int64_t(v_seg(0)) * int64_t(v_pt(0)) + int64_t(v_seg(1)) * int64_t(v_pt(1));
-        if (t_pt < 0) {
-            // Closest to p1.
-            double dabs = sqrt(int64_t(v_pt(0)) * int64_t(v_pt(0)) + int64_t(v_pt(1)) * int64_t(v_pt(1)));
-            if (dabs < d_min) {
-                d_min  = dabs;
-                i_min  = i;
-                pt_min = p1;
-            }
-        }
-        else if (t_pt > l2_seg) {
-            // Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the next step.
-            continue;
-        } else {
-            // Closest to the segment.
-            assert(t_pt >= 0 && t_pt <= l2_seg);
-            int64_t d_seg = int64_t(v_seg(1)) * int64_t(v_pt(0)) - int64_t(v_seg(0)) * int64_t(v_pt(1));
-            double d = double(d_seg) / sqrt(double(l2_seg));
-            double dabs = std::abs(d);
-            if (dabs < d_min) {
-                d_min  = dabs;
-                i_min  = i;
-                // Evaluate the foot point.
-                pt_min = p1;
-                double linv = double(d_seg) / double(l2_seg);
-                pt_min(0) = pt(0) - coord_t(floor(double(v_seg(1)) * linv + 0.5));
-                pt_min(1) = pt(1) + coord_t(floor(double(v_seg(0)) * linv + 0.5));
-                assert(Line(p1, p2).distance_to(pt_min) < scale_(1e-5));
-            }
-        }
+    const Vec3f& binormal() const {
+        return mX;
     }
 
-    assert(i_min != size_t(-1));
-    if ((pt_min - polygon.points[i_min]).cast<double>().norm() > eps) {
-        // Insert a new point on the segment i_min, i_min+1.
-        return polygon.points.insert(polygon.points.begin() + (i_min + 1), pt_min);
+    const Vec3f& tangent() const {
+        return mY;
     }
-    return polygon.points.begin() + i_min;
-}
 
+    const Vec3f& normal() const {
+        return mZ;
+    }
 
+private:
+    Vec3f mX, mY, mZ;
+};
 
-static std::vector<float> polygon_angles_at_vertices(const Polygon &polygon, const std::vector<float> &lengths, float min_arm_length)
-{
-    assert(polygon.points.size() + 1 == lengths.size());
-    if (min_arm_length > 0.25f * lengths.back())
-        min_arm_length = 0.25f * lengths.back();
+Vec3f sample_sphere_uniform(const Vec2f &samples) {
+    float term1 = 2.0f * float(PI) * samples.x();
+    float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
+    return {cos(term1) * term2, sin(term1) * term2,
+        1.0f - 2.0f * samples.y()};
+}
 
-    // Find the initial prev / next point span.
-    size_t idx_prev = polygon.points.size();
-    size_t idx_curr = 0;
-    size_t idx_next = 1;
-    while (idx_prev > idx_curr && lengths.back() - lengths[idx_prev] < min_arm_length)
-        -- idx_prev;
-    while (idx_next < idx_prev && lengths[idx_next] < min_arm_length)
-        ++ idx_next;
-
-    std::vector<float> angles(polygon.points.size(), 0.f);
-    for (; idx_curr < polygon.points.size(); ++ idx_curr) {
-        // Move idx_prev up until the distance between idx_prev and idx_curr is lower than min_arm_length.
-        if (idx_prev >= idx_curr) {
-            while (idx_prev < polygon.points.size() && lengths.back() - lengths[idx_prev] + lengths[idx_curr] > min_arm_length)
-                ++ idx_prev;
-            if (idx_prev == polygon.points.size())
-                idx_prev = 0;
-        }
-        while (idx_prev < idx_curr && lengths[idx_curr] - lengths[idx_prev] > min_arm_length)
-            ++ idx_prev;
-        // Move idx_prev one step back.
-        if (idx_prev == 0)
-            idx_prev = polygon.points.size() - 1;
-        else
-            -- idx_prev;
-        // Move idx_next up until the distance between idx_curr and idx_next is greater than min_arm_length.
-        if (idx_curr <= idx_next) {
-            while (idx_next < polygon.points.size() && lengths[idx_next] - lengths[idx_curr] < min_arm_length)
-                ++ idx_next;
-            if (idx_next == polygon.points.size())
-                idx_next = 0;
-        }
-        while (idx_next < idx_curr && lengths.back() - lengths[idx_curr] + lengths[idx_next] < min_arm_length)
-            ++ idx_next;
-        // Calculate angle between idx_prev, idx_curr, idx_next.
-        const Point &p0 = polygon.points[idx_prev];
-        const Point &p1 = polygon.points[idx_curr];
-        const Point &p2 = polygon.points[idx_next];
-        const Point  v1 = p1 - p0;
-        const Point  v2 = p2 - p1;
-        int64_t dot   = int64_t(v1(0))*int64_t(v2(0)) + int64_t(v1(1))*int64_t(v2(1));
-        int64_t cross = int64_t(v1(0))*int64_t(v2(1)) - int64_t(v1(1))*int64_t(v2(0));
-        float angle = float(atan2(double(cross), double(dot)));
-        angles[idx_curr] = angle;
-    }
+Vec3f sample_hemisphere_uniform(const Vec2f &samples) {
+    float term1 = 2.0f * float(PI) * samples.x();
+    float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
+    return {cos(term1) * term2, sin(term1) * term2,
+        abs(1.0f - 2.0f * samples.y())};
+}
 
-    return angles;
+Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) {
+    float term1 = 2.f * float(PI) * samples.x();
+    float term2 = pow(samples.y(), 1.f / (power + 1.f));
+    float term3 = sqrt(1.f - term2 * term2);
+
+    return Vec3f(cos(term1) * term3, sin(term1) * term3, term2);
 }
 
+std::vector<FaceVisibilityInfo> raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
+        const indexed_triangle_set &triangles, size_t negative_volumes_start_index) {
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start";
+
+    //prepare uniform samples of a hemisphere
+    float step_size = 1.0f / SeamPlacer::sqr_rays_per_triangle;
+    std::vector<Vec3f> precomputed_sample_directions(
+            SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle);
+    for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_triangle; ++x_idx) {
+        float sample_x = x_idx * step_size + step_size / 2.0;
+        for (size_t y_idx = 0; y_idx < SeamPlacer::sqr_rays_per_triangle; ++y_idx) {
+            size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_triangle + y_idx;
+            float sample_y = y_idx * step_size + step_size / 2.0;
+            precomputed_sample_directions[dir_index] = sample_hemisphere_uniform( { sample_x, sample_y });
+        }
+    }
 
+    bool model_contains_negative_parts = negative_volumes_start_index < triangles.indices.size();
+
+    std::vector<FaceVisibilityInfo> result(triangles.indices.size());
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, result.size()),
+            [&triangles, &precomputed_sample_directions, model_contains_negative_parts, negative_volumes_start_index,
+                    &raycasting_tree, &result](tbb::blocked_range<size_t> r) {
+                // Maintaining hits memory outside of the loop, so it does not have to be reallocated for each query.
+                std::vector<igl::Hit> hits;
+                for (size_t face_index = r.begin(); face_index < r.end(); ++face_index) {
+                    FaceVisibilityInfo &dest = result[face_index];
+                    dest.visibility = 1.0f;
+                    constexpr float decrease = 1.0f
+                            / (SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle);
+
+                    Vec3i face = triangles.indices[face_index];
+                    Vec3f A = triangles.vertices[face.x()];
+                    Vec3f B = triangles.vertices[face.y()];
+                    Vec3f C = triangles.vertices[face.z()];
+                    Vec3f center = (A + B + C) / 3.0f;
+                    Vec3f normal = ((B - A).cross(C - B)).normalized();
+                    // apply the local direction via Frame struct - the local_dir is with respect to +Z being forward
+                    Frame f;
+                    f.set_from_z(normal);
+
+                    for (const auto &dir : precomputed_sample_directions) {
+                        Vec3f final_ray_dir = (f.to_world(dir));
+                        if (!model_contains_negative_parts) {
+                            igl::Hit hitpoint;
+                            // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and
+                            // direction.
+                            Vec3d final_ray_dir_d = final_ray_dir.cast<double>();
+                            Vec3d ray_origin_d = (center + normal * 0.1).cast<double>(); // start above surface.
+                            bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices,
+                                    triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint);
+                            if (hit) {
+                                dest.visibility -= decrease;
+                            }
+                        } else { //TODO improve logic for order based boolean operations - consider order of volumes
+                            Vec3d ray_origin_d = (center + normal * 0.1).cast<double>(); // start above surface.
+                            if (face_index >= negative_volumes_start_index) { // if casting from negative volume face, invert direction, change start pos
+                                final_ray_dir = -1.0 * final_ray_dir;
+                                ray_origin_d = (center - normal * 0.1).cast<double>();
+                            }
+                            Vec3d final_ray_dir_d = final_ray_dir.cast<double>();
+                            bool some_hit = AABBTreeIndirect::intersect_ray_all_hits(triangles.vertices,
+                                    triangles.indices, raycasting_tree,
+                                    ray_origin_d, final_ray_dir_d, hits);
+                            if (some_hit) {
+                                int in_negative = 0;
+                                int in_positive = 0;
+                                // NOTE: iterating in reverse, from the last hit for one simple reason: We know the state of the ray at that point;
+                                //  It cannot be inside model, and it cannot be inside negative volume
+                                for (int hit_index = int(hits.size()) - 1; hit_index >= 0; --hit_index) {
+                                    Vec3f face_normal = its_face_normal(triangles, hits[hit_index].id);
+                                    if (hits[hit_index].id >= int(negative_volumes_start_index)) { //negative volume hit
+                                        in_negative += sgn(face_normal.dot(final_ray_dir)); // if volume face aligns with ray dir, we are leaving negative space
+                                        // which in reverse hit analysis means, that we are entering negative space :) and vice versa
+                                    } else {
+                                        in_positive += sgn(face_normal.dot(final_ray_dir));
+                                    }
+                                    if (in_positive > 0 && in_negative <= 0) {
+                                        dest.visibility -= decrease;
+                                        break;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            });
 
-void SeamPlacer::init(const Print& print)
-{
-    m_enforcers.clear();
-    m_blockers.clear();
-    m_seam_history.clear();
-    m_po_list.clear();
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end";
 
-    const std::vector<double>& nozzle_dmrs = print.config().nozzle_diameter.values;
-    float max_nozzle_dmr = *std::max_element(nozzle_dmrs.begin(), nozzle_dmrs.end());
+    return result;
+}
 
+std::vector<float> calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector<float> &lengths,
+        float min_arm_length) {
+    std::vector<float> result(polygon.size());
 
-    std::vector<ExPolygons> temp_enf;
-    std::vector<ExPolygons> temp_blk;
-    std::vector<Polygons>   temp_polygons;
+    if (polygon.size() == 1) {
+        result[0] = 0.0f;
+    }
 
-    for (const PrintObject* po : print.objects()) {
+    size_t idx_prev = 0;
+    size_t idx_curr = 0;
+    size_t idx_next = 0;
 
-        auto merge_and_offset = [po, &temp_polygons, max_nozzle_dmr](EnforcerBlockerType type, std::vector<ExPolygons>& out) {
-            // Offset the triangles out slightly.
-            auto offset_out = [](Polygon& input, float offset) -> ExPolygons {
-                ClipperLib::Paths out(1);
-                std::vector<float>  deltas(input.points.size(), offset);
-                input.make_counter_clockwise();
-                out.front() = mittered_offset_path_scaled(input.points, deltas, 3.);
-                return ClipperPaths_to_Slic3rExPolygons(out, true); // perform union
-            };
+    float distance_to_prev = 0;
+    float distance_to_next = 0;
 
+    //push idx_prev far enough back as initialization
+    while (distance_to_prev < min_arm_length) {
+        idx_prev = Slic3r::prev_idx_modulo(idx_prev, polygon.size());
+        distance_to_prev += lengths[idx_prev];
+    }
 
-            temp_polygons.clear();
-            po->project_and_append_custom_facets(true, type, temp_polygons);
-            out.clear();
-            out.reserve(temp_polygons.size());
-            float offset = scale_(max_nozzle_dmr + po->config().elefant_foot_compensation);
-            for (Polygons &src : temp_polygons) {
-                out.emplace_back(ExPolygons());
-                for (Polygon& plg : src) {
-                    ExPolygons offset_explg = offset_out(plg, offset);
-                    if (! offset_explg.empty())
-                        out.back().emplace_back(std::move(offset_explg.front()));
-                }
+    for (size_t _i = 0; _i < polygon.size(); ++_i) {
+        // pull idx_prev to current as much as possible, while respecting the min_arm_length
+        while (distance_to_prev - lengths[idx_prev] > min_arm_length) {
+            distance_to_prev -= lengths[idx_prev];
+            idx_prev = Slic3r::next_idx_modulo(idx_prev, polygon.size());
+        }
 
-                offset = scale_(max_nozzle_dmr);
-            }
-        };
-        merge_and_offset(EnforcerBlockerType::BLOCKER, temp_blk);
-        merge_and_offset(EnforcerBlockerType::ENFORCER, temp_enf);
-
-        // Remember this PrintObject and initialize a store of enforcers and blockers for it.
-        m_po_list.push_back(po);
-        size_t po_idx = m_po_list.size() - 1;
-        m_enforcers.emplace_back(std::vector<CustomTrianglesPerLayer>(temp_enf.size()));
-        m_blockers.emplace_back(std::vector<CustomTrianglesPerLayer>(temp_blk.size()));
-
-        // A helper class to store data to build the AABB tree from.
-        class CustomTriangleRef {
-        public:
-            CustomTriangleRef(size_t idx,
-                              Point&& centroid,
-                              BoundingBox&& bb)
-                : m_idx{idx}, m_centroid{centroid},
-                  m_bbox{AlignedBoxType(bb.min, bb.max)}
-            {}
-            size_t idx() const              { return m_idx;      }
-            const Point& centroid() const   { return m_centroid; }
-            const TreeType::BoundingBox& bbox() const { return m_bbox; }
-
-        private:
-            size_t m_idx;
-            Point m_centroid;
-            AlignedBoxType m_bbox;
-        };
+        //push idx_next forward as far as needed
+        while (distance_to_next < min_arm_length) {
+            distance_to_next += lengths[idx_next];
+            idx_next = Slic3r::next_idx_modulo(idx_next, polygon.size());
+        }
 
-        // A lambda to extract the ExPolygons and save them into the member AABB tree.
-        // Will be called for enforcers and blockers separately.
-        auto add_custom = [](std::vector<ExPolygons>& src, std::vector<CustomTrianglesPerLayer>& dest) {
-            // Go layer by layer, and append all the ExPolygons into the AABB tree.
-            size_t layer_idx = 0;
-            for (ExPolygons& expolys_on_layer : src) {
-                CustomTrianglesPerLayer& layer_data = dest[layer_idx];
-                std::vector<CustomTriangleRef> triangles_data;
-                layer_data.polys.reserve(expolys_on_layer.size());
-                triangles_data.reserve(expolys_on_layer.size());
-
-                for (ExPolygon& expoly : expolys_on_layer) {
-                    if (expoly.empty())
-                        continue;
-                    layer_data.polys.emplace_back(std::move(expoly));
-                    triangles_data.emplace_back(layer_data.polys.size() - 1,
-                                                layer_data.polys.back().centroid(),
-                                                layer_data.polys.back().bounding_box());
-                }
-                // All polygons are saved, build the AABB tree for them.
-                layer_data.tree.build(std::move(triangles_data));
-                ++layer_idx;
-            }
-        };
+        // Calculate angle between idx_prev, idx_curr, idx_next.
+        const Point &p0 = polygon.points[idx_prev];
+        const Point &p1 = polygon.points[idx_curr];
+        const Point &p2 = polygon.points[idx_next];
+        result[idx_curr] = float(angle(p1 - p0, p2 - p1));
 
-        add_custom(temp_enf, m_enforcers.at(po_idx));
-        add_custom(temp_blk, m_blockers.at(po_idx));
+        // increase idx_curr by one
+        float curr_distance = lengths[idx_curr];
+        idx_curr++;
+        distance_to_prev += curr_distance;
+        distance_to_next -= curr_distance;
     }
+
+    return result;
 }
 
+// structure to store global information about the model - occlusion hits, enforcers, blockers
+struct GlobalModelInfo {
+    indexed_triangle_set model;
+    AABBTreeIndirect::Tree<3, float> model_tree;
+    std::vector<FaceVisibilityInfo> visiblity_info;
+    indexed_triangle_set enforcers;
+    indexed_triangle_set blockers;
+    AABBTreeIndirect::Tree<3, float> enforcers_tree;
+    AABBTreeIndirect::Tree<3, float> blockers_tree;
+
+    bool is_enforced(const Vec3f &position, float radius) const {
+        if (enforcers.empty()) {
+            return false;
+        }
+        float radius_sqr = radius * radius;
+        return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices,
+                enforcers_tree, position, radius_sqr);
+    }
 
+    bool is_blocked(const Vec3f &position, float radius) const {
+        if (blockers.empty()) {
+            return false;
+        }
+        float radius_sqr = radius * radius;
+        return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices,
+                blockers_tree, position, radius_sqr);
+    }
 
-void SeamPlacer::plan_perimeters(const std::vector<const ExtrusionEntity*> perimeters,
-                            const Layer& layer, SeamPosition seam_position,
-                            Point last_pos, coordf_t nozzle_dmr, const PrintObject* po,
-                            const EdgeGrid::Grid* lower_layer_edge_grid)
-{
-    // When printing the perimeters, we want the seams on external and internal perimeters to match.
-    // We have a list of perimeters in the order to be printed. Each internal perimeter must inherit
-    // the seam from the previous external perimeter.
+    float calculate_point_visibility(const Vec3f &position) const {
+        size_t hit_idx;
+        Vec3f hit_point;
+        if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree,
+                position, hit_idx, hit_point) >= 0) {
+            return visiblity_info[hit_idx].visibility;
+        } else {
+            return 0.0f;
+        }
 
-    m_plan.clear();
-    m_plan_idx = 0;
+    }
 
-    if (perimeters.empty() || ! po)
-        return;
+#ifdef DEBUG_FILES
+    void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const {
+        indexed_triangle_set divided_mesh = obj_mesh;
+        Slic3r::CNumericLocalesSetter locales_setter;
 
-    m_plan.resize(perimeters.size());
+        FILE *fp = boost::nowide::fopen(file_name, "w");
+        if (fp == nullptr) {
+            BOOST_LOG_TRIVIAL(error)
+            << "stl_write_obj: Couldn't open " << file_name << " for writing";
+            return;
+        }
 
-    for (int i = 0; i < int(perimeters.size()); ++i) {
-        if (perimeters[i]->role() == erExternalPerimeter && perimeters[i]->is_loop()) {
-            last_pos = this->calculate_seam(
-                layer, seam_position, *dynamic_cast<const ExtrusionLoop*>(perimeters[i]), nozzle_dmr,
-                po, lower_layer_edge_grid, last_pos, false);
-            m_plan[i].external = true;
+        for (size_t i = 0; i < divided_mesh.vertices.size(); ++i) {
+            float visibility = calculate_point_visibility(divided_mesh.vertices[i]);
+            Vec3f color = value_to_rgbf(0.0f, 1.0f,
+                    visibility);
+            fprintf(fp, "v %f %f %f  %f %f %f\n",
+                    divided_mesh.vertices[i](0), divided_mesh.vertices[i](1), divided_mesh.vertices[i](2),
+                    color(0), color(1), color(2)
+                            );
         }
-        m_plan[i].seam_position = seam_position;
-        m_plan[i].layer = &layer;
-        m_plan[i].po = po;
-        m_plan[i].pt = last_pos;
+        for (size_t i = 0; i < divided_mesh.indices.size(); ++i)
+            fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1,
+                    divided_mesh.indices[i][2] + 1);
+        fclose(fp);
+
     }
+#endif
 }
-
-
-void SeamPlacer::place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
-                            const EdgeGrid::Grid* lower_layer_edge_grid)
-{
-    // const double seam_offset = nozzle_diameter;
-
-    Point seam = last_pos;
-    if (! m_plan.empty() && m_plan_idx < m_plan.size()) {
-        if (m_plan[m_plan_idx].external) {
-            seam = m_plan[m_plan_idx].pt;
-            // One more heuristics: if the seam is too far from current nozzle position,
-            // try to place it again. This can happen in cases where the external perimeter
-            // does not belong to the preceding ones and they are ordered so they end up
-            // far from each other.
-            if ((seam.cast<double>() - last_pos.cast<double>()).squaredNorm() > std::pow(scale_(5.*nozzle_diameter), 2.))
-                seam = this->calculate_seam(*m_plan[m_plan_idx].layer, m_plan[m_plan_idx].seam_position, loop, nozzle_diameter,
-                    m_plan[m_plan_idx].po, lower_layer_edge_grid, last_pos, false);
-
-            if (m_plan[m_plan_idx].seam_position == spAligned)
-                m_seam_history.add_seam(m_plan[m_plan_idx].po, m_plan[m_plan_idx].pt, loop.polygon().bounding_box());
-        }
-        else {
-            if (!external_first) {
-                // Internal perimeter printed before the external.
-                // First get list of external seams.
-                std::vector<size_t> ext_seams;
-                size_t external_cnt = 0;
-                for (size_t i = 0; i < m_plan.size(); ++i) {
-                    if (m_plan[i].external) {
-                        ext_seams.emplace_back(i);
-                        ++external_cnt;
+;
+
+//Extract perimeter polygons of the given layer
+Polygons extract_perimeter_polygons(const Layer *layer, const SeamPosition configured_seam_preference,
+        std::vector<const LayerRegion*> &corresponding_regions_out) {
+    Polygons polygons;
+    for (const LayerRegion *layer_region : layer->regions()) {
+        for (const ExtrusionEntity *ex_entity : layer_region->perimeters.entities) {
+            if (ex_entity->is_collection()) { //collection of inner, outer, and overhang perimeters
+                for (const ExtrusionEntity *perimeter : static_cast<const ExtrusionEntityCollection*>(ex_entity)->entities) {
+                    if (perimeter->role() == ExtrusionRole::erExternalPerimeter
+                            || (perimeter->role() == ExtrusionRole::erPerimeter
+                                    && configured_seam_preference == spRandom)) { //for random seam alignment, extract all perimeters
+                        Points p;
+                        perimeter->collect_points(p);
+                        polygons.emplace_back(std::move(p));
+                        corresponding_regions_out.push_back(layer_region);
                     }
                 }
-
-                if (!ext_seams.empty()) {
-                    // First find the line segment closest to an external seam:
-                    //int path_idx = 0;
-                    //int line_idx = 0;
-                    size_t ext_seam_idx = size_t(-1);
-                    double min_dist_sqr = std::numeric_limits<double>::max();
-                    std::vector<Lines> lines_vect;
-                    for (int i = 0; i < int(loop.paths.size()); ++i) {
-                        lines_vect.emplace_back(loop.paths[i].polyline.lines());
-                        const Lines& lines = lines_vect.back();
-                        for (int j = 0; j < int(lines.size()); ++j) {
-                            for (size_t k : ext_seams) {
-                                double d_sqr = lines[j].distance_to_squared(m_plan[k].pt);
-                                if (d_sqr < min_dist_sqr) {
-                                    //path_idx = i;
-                                    //line_idx = j;
-                                    ext_seam_idx = k;
-                                    min_dist_sqr = d_sqr;
-                                }
-                            }
-                        }
-                    }
-
-                    // Only accept seam that is reasonably close.
-                    if (ext_seam_idx != size_t(-1)) {
-                        // How many nozzle diameters is considered "close"?
-                        const double nozzle_d_limit = 2. * (1. + m_plan.size() / external_cnt);
-                        const double limit_dist_sqr = double(scale_(scale_((unscale(m_plan[ext_seam_idx].pt) - unscale(m_plan[m_plan_idx].pt)).squaredNorm() * std::pow(nozzle_d_limit * nozzle_diameter, 2.))));
-
-                        if (min_dist_sqr < limit_dist_sqr) {
-                            // Now find a projection of the external seam
-                            //const Lines& lines = lines_vect[path_idx];
-                            //Point closest = m_plan[ext_seam_idx].pt.projection_onto(lines[line_idx]);
-
-                            //                        This code does staggering of internal perimeters, turned off for now.
-                            // 
-                            //                       double dist = (closest.cast<double>() - lines[line_idx].b.cast<double>()).norm();
-                            //
-                            //                       // And walk along the perimeter until we make enough space for
-                            //                       // seams of all perimeters beforethe external one.
-                            //                       double offset = (ext_seam_idx - m_plan_idx) * scale_(seam_offset);
-                            //                       double last_offset = offset;
-                            //                       offset -= dist;
-                            //                       const Point* a = &closest;
-                            //                       const Point* b = &lines[line_idx].b;
-                            //                       while (++line_idx < int(lines.size()) && offset > 0.) {
-                            //                           last_offset = offset;
-                            //                           offset -= lines[line_idx].length();
-                            //                           a = &lines[line_idx].a;
-                            //                           b = &lines[line_idx].b;
-                            //                       }
-                            //
-                            //                        // We have walked far enough, too far maybe. Interpolate on the
-                            //                       // last segment to find the end precisely.
-                            //                        offset = std::min(0., offset); // In case that offset is still positive (we may have "wrapped around")
-                            //                        double ratio = last_offset / (last_offset - offset);
-                            //                        seam = (a->cast<double>() + ((b->cast<double>() - a->cast<double>()) * ratio)).cast<coord_t>();
-                            seam = m_plan[ext_seam_idx].pt;
-                        }
-                    }
+                if (polygons.empty()) {
+                    Points p;
+                    ex_entity->collect_points(p);
+                    polygons.emplace_back(std::move(p));
+                    corresponding_regions_out.push_back(layer_region);
                 }
+            } else {
+                Points p;
+                ex_entity->collect_points(p);
+                polygons.emplace_back(std::move(p));
+                corresponding_regions_out.push_back(layer_region);
             }
-            else {
-                // We should have a candidate ready from before. If not, use last_pos.
-                if (m_plan_idx > 0 && m_plan[m_plan_idx - 1].precalculated)
-                    seam = m_plan[m_plan_idx - 1].pt;
-            }
+        }
+    }
+
+    if (polygons.empty()) { // If there are no perimeter polygons for whatever reason (disabled perimeters .. ) insert dummy point
+        // it is easier than checking everywhere if the layer is not emtpy, no seam will be placed to this layer anyway
+        polygons.emplace_back(std::vector { Point { 0, 0 } });
+        corresponding_regions_out.push_back(nullptr);
+    }
 
-            // seam now contains a hot candidate for internal seam. Use it unless there is a sharp corner nearby.
-            // We will call the normal seam planning function, pretending that we are currently at the candidate point
-            // and set to spNearest. If the ideal seam it finds is close to current candidate, use it.
-            // This is to prevent having seams very close to corners, just because of external seam position.
-            seam = calculate_seam(*m_plan[m_plan_idx].layer, spNearest, loop, nozzle_diameter,
-                m_plan[m_plan_idx].po, lower_layer_edge_grid, seam, true);
-        }
-        m_plan[m_plan_idx].pt = seam;
-    }
-
-
-    // Split the loop at the point with a minium penalty.
-    if (!loop.split_at_vertex(seam))
-        // The point is not in the original loop. Insert it.
-        loop.split_at(seam, true);
-
-    if (external_first && m_plan_idx+1<m_plan.size() && ! m_plan[m_plan_idx+1].external) {
-//        This code does staggering of internal perimeters, turned off for now.
-//        Next perimeter should start near this one.
-//        const double dist_sqr = std::pow(double(scale_(seam_offset)), 2.);
-//        double running_sqr = 0.;
-//        double running_sqr_last = 0.;
-//        if (!loop.paths.empty() && loop.paths.back().polyline.points.size() > 1) {
-//            const ExtrusionPath& last = loop.paths.back();
-//            auto it = last.polyline.points.crbegin() + 1;
-//            for (; it != last.polyline.points.crend(); ++it) {
-//                running_sqr += (it->cast<double>() - (it - 1)->cast<double>()).squaredNorm();
-//                if (running_sqr > dist_sqr)
-//                    break;
-//                running_sqr_last = running_sqr;
-//            }
-//            if (running_sqr <= dist_sqr)
-//                it = last.polyline.points.crend() - 1;
-//            // Now interpolate.
-//            double ratio = (std::sqrt(dist_sqr) - std::sqrt(running_sqr_last)) / (std::sqrt(running_sqr) - std::sqrt(running_sqr_last));
-//            m_plan[m_plan_idx + 1].pt = ((it - 1)->cast<double>() + (it->cast<double>() - (it - 1)->cast<double>()) * std::min(ratio, 1.)).cast<coord_t>();
-//            m_plan[m_plan_idx + 1].precalculated = true;
-            m_plan[m_plan_idx + 1].pt = m_plan[m_plan_idx].pt;
-            m_plan[m_plan_idx + 1].precalculated = true;
-//        }
-    }
-
-    ++m_plan_idx;
+    return polygons;
 }
 
+// Insert SeamCandidates created from perimeter polygons in to the result vector.
+// Compute its type (Enfrocer,Blocker), angle, and position
+//each SeamCandidate also contains pointer to shared Perimeter structure representing the polygon
+// if Custom Seam modifiers are present, oversamples the polygon if necessary to better fit user intentions
+void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const LayerRegion *region,
+        const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) {
+    if (orig_polygon.size() == 0) {
+        return;
+    }
 
-// Returns "best" seam for a given perimeter.
-Point SeamPlacer::calculate_seam(const Layer& layer, const SeamPosition seam_position,
-               const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po,
-               const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest)
-{
-    Polygon polygon = loop.polygon();
+    Polygon polygon = orig_polygon;
     bool was_clockwise = polygon.make_counter_clockwise();
-    const coord_t  nozzle_r   = coord_t(scale_(0.5 * nozzle_dmr) + 0.5);
 
-    size_t po_idx = std::find(m_po_list.begin(), m_po_list.end(), po) - m_po_list.begin();
-
-    // Find current layer in respective PrintObject. Cache the result so the
-    // lookup is only done once per layer, not for each loop.
-    const Layer* layer_po = nullptr;
-    if (po == m_last_po && layer.print_z == m_last_print_z)
-        layer_po = m_last_layer_po;
-    else {
-        layer_po = po ? po->get_layer_at_printz(layer.print_z) : nullptr;
-        m_last_po = po;
-        m_last_print_z = layer.print_z;
-        m_last_layer_po = layer_po;
+    std::vector<float> lengths { };
+    for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) {
+        lengths.push_back(std::max((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm(), 0.01));
     }
-    if (! layer_po)
-        return last_pos;
-
-    // Index of this layer in the respective PrintObject.
-    size_t layer_idx = layer_po->id() - po->layers().front()->id(); // raft layers
+    lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.01));
 
-    assert(layer_idx < po->layer_count());
+    std::vector<float> local_angles = calculate_polygon_angles_at_vertices(polygon, lengths,
+            SeamPlacer::polygon_local_angles_arm_distance);
 
-    const bool custom_seam = loop.role() == erExternalPerimeter && this->is_custom_seam_on_layer(layer_idx, po_idx);
+    result.perimeters.push_back( { });
+    Perimeter &perimeter = result.perimeters.back();
 
-    if (custom_seam) {
-        // Seam enf/blockers can begin and end in between the original vertices.
-        // Let add extra points in between and update the leghths.
-        polygon.densify(MINIMAL_POLYGON_SIDE);
+    std::queue<Vec3f> orig_polygon_points { };
+    for (size_t index = 0; index < polygon.size(); ++index) {
+        Vec2f unscaled_p = unscale(polygon[index]).cast<float>();
+        orig_polygon_points.emplace(unscaled_p.x(), unscaled_p.y(), z_coord);
     }
+    Vec3f first = orig_polygon_points.front();
+    std::queue<Vec3f> oversampled_points { };
+    size_t orig_angle_index = 0;
+    perimeter.start_index = result.points.size();
+    perimeter.flow_width = region != nullptr ? region->flow(FlowRole::frExternalPerimeter).width() : 0.0f;
+    bool some_point_enforced = false;
+    while (!orig_polygon_points.empty() || !oversampled_points.empty()) {
+        EnforcedBlockedSeamPoint type = EnforcedBlockedSeamPoint::Neutral;
+        Vec3f position;
+        float local_ccw_angle = 0;
+        bool orig_point = false;
+        if (!oversampled_points.empty()) {
+            position = oversampled_points.front();
+            oversampled_points.pop();
+        } else {
+            position = orig_polygon_points.front();
+            orig_polygon_points.pop();
+            local_ccw_angle = was_clockwise ? -local_angles[orig_angle_index] : local_angles[orig_angle_index];
+            orig_angle_index++;
+            orig_point = true;
+        }
+
+        if (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) {
+            type = EnforcedBlockedSeamPoint::Enforced;
+            some_point_enforced = true;
+        }
 
-    if (seam_position != spRandom) {
-        // Retrieve the last start position for this object.
-        float last_pos_weight = 1.f;
+        if (global_model_info.is_blocked(position, SeamPlacer::enforcer_blocker_distance_tolerance)) {
+            type = EnforcedBlockedSeamPoint::Blocked;
+        }
 
-        if (seam_position == spAligned) {
-            // Seam is aligned to the seam at the preceding layer.
-            if (po != nullptr) {
-                std::optional<Point> pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, loop.polygon().bounding_box());
-                if (pos.has_value()) {
-                    last_pos = *pos;
-                    last_pos_weight = is_custom_enforcer_on_layer(layer_idx, po_idx) ? 0.f : 1.f;
+        if (orig_point) {
+            Vec3f pos_of_next = orig_polygon_points.empty() ? first : orig_polygon_points.front();
+            float distance_to_next = (position - pos_of_next).norm();
+            if (global_model_info.is_enforced(position, distance_to_next)) {
+                Vec3f vec_to_next = (pos_of_next - position).normalized();
+                float step_size = SeamPlacer::enforcer_oversampling_distance;
+                float step = step_size;
+                while (step < distance_to_next) {
+                    oversampled_points.push(position + vec_to_next * step);
+                    step += step_size;
                 }
             }
         }
-        else if (seam_position == spRear) {
-            // Object is centered around (0,0) in its current coordinate system.
-            last_pos.x() = 0;
-            last_pos.y() = coord_t(3. * po->bounding_box().radius());
-            last_pos_weight = 5.f;
-        } if (seam_position == spNearest) {
-            // last_pos already contains current nozzle position
-        }
-
-        // Insert a projection of last_pos into the polygon.
-        size_t last_pos_proj_idx;
-        {
-            auto it = project_point_to_polygon_and_insert(polygon, last_pos, 0.1 * nozzle_r);
-            last_pos_proj_idx = it - polygon.points.begin();
-        }
-
-        // Parametrize the polygon by its length.
-        std::vector<float> lengths = polygon.parameter_by_length();
-
-        // For each polygon point, store a penalty.
-        // First calculate the angles, store them as penalties. The angles are caluculated over a minimum arm length of nozzle_r.
-        std::vector<float> penalties = polygon_angles_at_vertices(polygon, lengths, float(nozzle_r));
-        // No penalty for reflex points, slight penalty for convex points, high penalty for flat surfaces.
-        const float penaltyConvexVertex = 1.f;
-        const float penaltyFlatSurface = 3.f;
-        const float penaltyOverhangHalf = 10.f;
-        // Penalty for visible seams.
-       for (size_t i = 0; i < polygon.points.size(); ++ i) {
-            float ccwAngle = penalties[i];
-            if (was_clockwise)
-                ccwAngle = - ccwAngle;
-            float penalty = 0;
-            if (ccwAngle <- float(0.6 * PI))
-                // Sharp reflex vertex. We love that, it hides the seam perfectly.
-                penalty = 0.f;
-            else if (ccwAngle > float(0.6 * PI))
-                // Seams on sharp convex vertices are more visible than on reflex vertices.
-                penalty = penaltyConvexVertex;
-            else if (ccwAngle < 0.f) {
-                // Interpolate penalty between maximum and zero.
-                penalty = penaltyFlatSurface * bspline_kernel(ccwAngle * float(PI * 2. / 3.));
-            } else {
-                assert(ccwAngle >= 0.f);
-                // Interpolate penalty between maximum and the penalty for a convex vertex.
-                penalty = penaltyConvexVertex + (penaltyFlatSurface - penaltyConvexVertex) * bspline_kernel(ccwAngle * float(PI * 2. / 3.));
-            }
-            // Give a negative penalty for points close to the last point or the prefered seam location.
-            float dist_to_last_pos_proj = (i < last_pos_proj_idx) ?
-                std::min(lengths[last_pos_proj_idx] - lengths[i], lengths.back() - lengths[last_pos_proj_idx] + lengths[i]) :
-                std::min(lengths[i] - lengths[last_pos_proj_idx], lengths.back() - lengths[i] + lengths[last_pos_proj_idx]);
-            float dist_max = 0.1f * lengths.back(); // 5.f * nozzle_dmr
-            penalty -= last_pos_weight * bspline_kernel(dist_to_last_pos_proj / dist_max);
-            penalties[i] = std::max(0.f, penalty);
-            if (prefer_nearest) {
-                // This hack limits the search around the nearest position projection.
-                penalties[i] += dist_to_last_pos_proj > 6.f * nozzle_r ? 100.f : 0.f;
+
+        result.points.emplace_back(position, perimeter, local_ccw_angle, type);
+    }
+
+    perimeter.end_index = result.points.size() - 1;
+
+    // We will find first patch of enforced points (patch: continuous section of enforced points) and select the middle
+    //      point, which will have priority during alignment
+    // If there are multiple enforced patches in the perimeter, others are ignored
+    if (some_point_enforced) {
+        size_t perimeter_size = perimeter.end_index - perimeter.start_index + 1;
+        const auto next_index = [&](size_t idx) {
+            return perimeter.start_index + Slic3r::next_idx_modulo(idx - perimeter.start_index, perimeter_size);
+        };
+
+        size_t first_enforced_idx = perimeter.start_index;
+        for (size_t _ = 0; _ < perimeter_size; ++_) {
+            if (result.points[first_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced &&
+                    result.points[next_index(first_enforced_idx)].type == EnforcedBlockedSeamPoint::Enforced) {
+                break;
             }
+            first_enforced_idx = next_index(first_enforced_idx);
         }
-       
-        // Penalty for overhangs.
-        if (lower_layer_edge_grid) {
-            // Use the edge grid distance field structure over the lower layer to calculate overhangs.
-            coord_t nozzle_r = coord_t(std::floor(scale_(0.5 * nozzle_dmr) + 0.5));
-            coord_t search_r = coord_t(std::floor(scale_(0.8 * nozzle_dmr) + 0.5));
-            for (size_t i = 0; i < polygon.points.size(); ++ i) {
-                const Point &p = polygon.points[i];
-                coordf_t dist;
-                // Signed distance is positive outside the object, negative inside the object.
-                // The point is considered at an overhang, if it is more than nozzle radius
-                // outside of the lower layer contour.
-                [[maybe_unused]] bool found = lower_layer_edge_grid->signed_distance(p, search_r, dist);
-                // If the approximate Signed Distance Field was initialized over lower_layer_edge_grid,
-                // then the signed distnace shall always be known.
-                assert(found);
-                penalties[i] += extrudate_overlap_penalty(float(nozzle_r), penaltyOverhangHalf, float(dist));
+        first_enforced_idx = next_index(first_enforced_idx);
+
+        // Gather also points with large angles (these are points from the original mesh, since oversampled points have zero angle)
+        // If there are any, the middle point will be picked from those (makes drawing over sharp corners easier)
+        std::vector<size_t> orig_large_angle_points_indices { };
+        std::vector<size_t> viable_points_indices { };
+        size_t last_enforced_idx = first_enforced_idx;
+        for (size_t _ = 0; _ < perimeter_size; ++_) {
+            if (result.points[last_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced) {
+                break;
             }
-        }
-        
-        // Custom seam. Huge (negative) constant penalty is applied inside
-        // blockers (enforcers) to rule out points that should not win.
-        if (custom_seam)
-            this->apply_custom_seam(polygon, po_idx, penalties, lengths, layer_idx, seam_position);
-
-        // Find a point with a minimum penalty.
-        size_t idx_min = std::min_element(penalties.begin(), penalties.end()) - penalties.begin();
-
-        if (seam_position != spAligned || ! is_custom_enforcer_on_layer(layer_idx, po_idx)) {
-            // Very likely the weight of idx_min is very close to the weight of last_pos_proj_idx.
-            // In that case use last_pos_proj_idx instead.
-            float penalty_aligned  = penalties[last_pos_proj_idx];
-            float penalty_min      = penalties[idx_min];
-            float penalty_diff_abs = std::abs(penalty_min - penalty_aligned);
-            float penalty_max      = std::max(std::abs(penalty_min), std::abs(penalty_aligned));
-            float penalty_diff_rel = (penalty_max == 0.f) ? 0.f : penalty_diff_abs / penalty_max;
-            // printf("Align seams, penalty aligned: %f, min: %f, diff abs: %f, diff rel: %f\n", penalty_aligned, penalty_min, penalty_diff_abs, penalty_diff_rel);
-            if (std::abs(penalty_diff_rel) < 0.05) {
-                // Penalty of the aligned point is very close to the minimum penalty.
-                // Align the seams as accurately as possible.
-                idx_min = last_pos_proj_idx;
+            viable_points_indices.push_back(last_enforced_idx);
+            if (abs(result.points[last_enforced_idx].local_ccw_angle) > 0.4 * PI) {
+                orig_large_angle_points_indices.push_back(last_enforced_idx);
             }
+            last_enforced_idx = next_index(last_enforced_idx);
         }
+        assert(viable_points_indices.size() > 0);
+        if (orig_large_angle_points_indices.empty()) {
+            size_t central_idx = viable_points_indices[viable_points_indices.size() / 2];
+            result.points[central_idx].central_enforcer = true;
+        } else {
+            size_t central_idx = orig_large_angle_points_indices.size() / 2;
+            result.points[orig_large_angle_points_indices[central_idx]].central_enforcer = true;
+        }
+    }
 
-
-        // Export the contour into a SVG file.
-        #if 0
-        {
-            static int iRun = 0;
-            SVG svg(debug_out_path("GCode_extrude_loop-%d.svg", iRun ++));
-            if (m_layer->lower_layer != NULL)
-                svg.draw(m_layer->lower_layer->slices);
-            for (size_t i = 0; i < loop.paths.size(); ++ i)
-                svg.draw(loop.paths[i].as_polyline(), "red");
-            Polylines polylines;
-            for (size_t i = 0; i < loop.paths.size(); ++ i)
-                polylines.push_back(loop.paths[i].as_polyline());
-            Slic3r::Polygons polygons;
-            coordf_t nozzle_dmr = EXTRUDER_CONFIG(nozzle_diameter);
-            coord_t delta = scale_(0.5*nozzle_dmr);
-            Slic3r::offset(polylines, &polygons, delta);
-//            for (size_t i = 0; i < polygons.size(); ++ i) svg.draw((Polyline)polygons[i], "blue");
-            svg.draw(last_pos, "green", 3);
-            svg.draw(polygon.points[idx_min], "yellow", 3);
-            svg.Close();
-        }
-        #endif
-        return polygon.points[idx_min];
-
-    } else
-        return this->get_random_seam(layer_idx, polygon, po_idx);
 }
 
+// Get index of previous and next perimeter point of the layer. Because SeamCandidates of all polygons of the given layer
+// are sequentially stored in the vector, each perimeter contains info about start and end index. These vales are used to
+// deduce index of previous and next neigbour in the corresponding perimeter.
+std::pair<size_t, size_t> find_previous_and_next_perimeter_point(const std::vector<SeamCandidate> &perimeter_points,
+        size_t point_index) {
+    const SeamCandidate &current = perimeter_points[point_index];
+    int prev = point_index - 1; //for majority of points, it is true that neighbours lie behind and in front of them in the vector
+    int next = point_index + 1;
+
+    if (point_index == current.perimeter.start_index) {
+        // if point_index is equal to start, it means that the previous neighbour is at the end
+        prev = current.perimeter.end_index;
+    }
 
-Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_idx,
-                                  bool* saw_custom) const
-{
-    // Parametrize the polygon by its length.
-    const std::vector<float> lengths = polygon.parameter_by_length();
-
-    // Which of the points are inside enforcers/blockers?
-    std::vector<size_t> enforcers_idxs;
-    std::vector<size_t> blockers_idxs;
-    this->get_enforcers_and_blockers(layer_idx, polygon, po_idx, enforcers_idxs, blockers_idxs);
-
-    bool has_enforcers = ! enforcers_idxs.empty();
-    bool has_blockers = ! blockers_idxs.empty();
-    if (saw_custom)
-        *saw_custom = has_enforcers || has_blockers;
-
-    assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end()));
-    assert(std::is_sorted(blockers_idxs.begin(), blockers_idxs.end()));
-    std::vector<float> edges;
-
-    // Lambda to calculate lengths of all edges of interest. Last parameter
-    // decides whether to measure edges inside or outside idxs.
-    // Negative number = not an edge of interest.
-    auto get_valid_length = [&lengths](const std::vector<size_t>& idxs,
-                                       std::vector<float>& edges,
-                                       bool measure_inside_edges) -> float
-    {
-        // First mark edges we are interested in by assigning a positive number.
-        edges.assign(lengths.size()-1, measure_inside_edges ? -1.f : 1.f);
-        for (size_t i=0; i<idxs.size(); ++i) {
-            size_t this_pt_idx = idxs[i];
-            // Two concurrent indices in the list -> the edge between them is the enforcer/blocker.
-            bool inside_edge = ((i != idxs.size()-1 && idxs[i+1] == this_pt_idx + 1)
-                             || (i == idxs.size()-1 && idxs.back() == lengths.size()-2 && idxs[0] == 0));
-            if (inside_edge)
-                edges[this_pt_idx] = measure_inside_edges ? 1.f : -1.f;
-        }
-        // Now measure them.
-        float running_total = 0.f;
-        for (size_t i=0; i<edges.size(); ++i) {
-            if (edges[i] > 0.f) {
-                edges[i] = lengths[i+1] - lengths[i];
-                running_total += edges[i];
-            }
-        }
-        return running_total;
-    };
+    if (point_index == current.perimeter.end_index) {
+        // if point_index is equal to end, than next neighbour is at the start
+        next = current.perimeter.start_index;
+    }
 
-    // Find all seam candidate edges and their lengths.
-    float valid_length = 0.f;
-    if (has_enforcers)
-        valid_length = get_valid_length(enforcers_idxs, edges, true);
-
-    if (! has_enforcers || valid_length == 0.f) {
-        // Second condition covers case with isolated enf points. Given how the painted
-        // triangles are projected, this should not happen. Stay on the safe side though.
-        if (has_blockers)
-            valid_length = get_valid_length(blockers_idxs, edges, false);
-        if (valid_length == 0.f) // No blockers or everything blocked - use the whole polygon.
-            valid_length = lengths.back();
-    }
-    assert(valid_length != 0.f);
-    // Now generate a random length and find the respective edge.
-    float rand_len = valid_length * (rand()/float(RAND_MAX));
-    size_t pt_idx = 0; // Index of the edge where to put the seam.
-    if (valid_length == lengths.back()) {
-        // Whole polygon is used for placing the seam.
-        auto it = std::lower_bound(lengths.begin(), lengths.end(), rand_len);
-        pt_idx = it == lengths.begin() ? 0 : (it-lengths.begin()-1); // this takes care of a corner case where rand() returns 0
-    } else {
-        float running = 0.f;
-        for (size_t i=0; i<edges.size(); ++i) {
-            running += edges[i] > 0.f ? edges[i] : 0.f;
-            if (running >= rand_len) {
-                pt_idx = i;
-                break;
+    assert(prev >= 0);
+    assert(next >= 0);
+    return {size_t(prev),size_t(next)};
+}
+
+// Computes all global model info - transforms object, performs raycasting
+void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) {
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: gather occlusion meshes: start";
+    auto obj_transform = po->trafo_centered();
+    indexed_triangle_set triangle_set;
+    indexed_triangle_set negative_volumes_set;
+    //add all parts
+    for (const ModelVolume *model_volume : po->model_object()->volumes) {
+        if (model_volume->type() == ModelVolumeType::MODEL_PART
+                || model_volume->type() == ModelVolumeType::NEGATIVE_VOLUME) {
+            auto model_transformation = model_volume->get_matrix();
+            indexed_triangle_set model_its = model_volume->mesh().its;
+            its_transform(model_its, model_transformation);
+            if (model_volume->type() == ModelVolumeType::MODEL_PART) {
+                its_merge(triangle_set, model_its);
+            } else {
+                its_merge(negative_volumes_set, model_its);
             }
         }
     }
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: gather occlusion meshes: end";
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: simplify occlusion meshes: start";
+
+    //simplify raycasting mesh
+    its_quadric_edge_collapse(triangle_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr, nullptr,
+            nullptr);
+    triangle_set = its_subdivide(triangle_set, SeamPlacer::raycasting_subdivision_target_length);
+
+    //simplify negative volumes
+    its_quadric_edge_collapse(negative_volumes_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr,
+            nullptr,
+            nullptr);
+    negative_volumes_set = its_subdivide(negative_volumes_set, SeamPlacer::raycasting_subdivision_target_length);
+
+    size_t negative_volumes_start_index = triangle_set.indices.size();
+    its_merge(triangle_set, negative_volumes_set);
+    its_transform(triangle_set, obj_transform);
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: simplify occlusion meshes: end";
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer:build AABB tree: start";
+    auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices,
+            triangle_set.indices);
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer:build AABB tree: end";
+    result.model = triangle_set;
+    result.model_tree = raycasting_tree;
+    result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set, negative_volumes_start_index);
+
+#ifdef DEBUG_FILES
+    auto filename = debug_out_path(("visiblity_of_" + std::to_string(po->id().id) + ".obj").c_str());
+    result.debug_export(triangle_set, filename.c_str());
+#endif
+}
 
-    if (! has_enforcers && ! has_blockers) {
-        // The polygon may be too coarse, calculate the point exactly.
-        assert(valid_length == lengths.back());
-        bool last_seg = pt_idx == polygon.points.size()-1;
-        size_t next_idx = last_seg ? 0 : pt_idx+1;
-        const Point& prev = polygon.points[pt_idx];
-        const Point& next = polygon.points[next_idx];
-        assert(next_idx == 0 || pt_idx+1 == next_idx);
-        coordf_t diff_x = next.x() - prev.x();
-        coordf_t diff_y = next.y() - prev.y();
-        coordf_t dist = lengths[last_seg ? pt_idx+1 : next_idx] - lengths[pt_idx];
-        return Point(prev.x() + (rand_len - lengths[pt_idx]) * (diff_x/dist),
-                     prev.y() + (rand_len - lengths[pt_idx]) * (diff_y/dist));
+void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) {
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: start";
 
-    } else {
-        // The polygon should be dense enough.
-        return polygon.points[pt_idx];
+    auto obj_transform = po->trafo_centered();
+
+    for (const ModelVolume *mv : po->model_object()->volumes) {
+        if (mv->is_seam_painted()) {
+            auto model_transformation = obj_transform * mv->get_matrix();
+
+            indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER);
+            its_transform(enforcers, model_transformation);
+            its_merge(result.enforcers, enforcers);
+
+            indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER);
+            its_transform(blockers, model_transformation);
+            its_merge(result.blockers, blockers);
+        }
     }
+
+    result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices,
+            result.enforcers.indices);
+    result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices,
+            result.blockers.indices);
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end";
 }
 
+struct SeamComparator {
+    SeamPosition setup;
 
+    SeamComparator(SeamPosition setup) :
+            setup(setup) {
+    }
 
+    // Standard comparator, must respect the requirements of comparators (e.g. give same result on same inputs) for sorting usage
+    // should return if a is better seamCandidate than b
+    bool is_first_better(const SeamCandidate &a, const SeamCandidate &b, const Vec2f &preffered_location = Vec2f { 0.0f,
+            0.0f }) const {
+        if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) {
+            return a.central_enforcer;
+        }
 
+        // Blockers/Enforcers discrimination, top priority
+        if (a.type != b.type) {
+            return a.type > b.type;
+        }
 
+        //avoid overhangs
+        if (a.overhang > 0.0f || b.overhang > 0.0f) {
+            return a.overhang < b.overhang;
+        }
 
+        // prefer hidden points (more than 1 mm inside)
+        if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) {
+            return true;
+        }
+        if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) {
+            return false;
+        }
 
+        if (setup == SeamPosition::spRear) {
+            return a.position.y() > b.position.y();
+        }
 
-void SeamPlacer::get_enforcers_and_blockers(size_t layer_id,
-                             const Polygon& polygon,
-                             size_t po_idx,
-                             std::vector<size_t>& enforcers_idxs,
-                             std::vector<size_t>& blockers_idxs) const
-{
-    enforcers_idxs.clear();
-    blockers_idxs.clear();
+        float distance_penalty_a = 1.0f;
+        float distance_penalty_b = 1.0f;
+        if (setup == spNearest) {
+            distance_penalty_a = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f);
+            distance_penalty_b = 1.1f - gauss((b.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f);
+        }
 
-    auto is_inside = [](const Point& pt,
-                        const CustomTrianglesPerLayer& custom_data) -> bool {
-        assert(! custom_data.polys.empty());
-        // Now ask the AABB tree which polygons we should check and check them.
-        std::vector<size_t> candidates;
-        AABBTreeIndirect::get_candidate_idxs(custom_data.tree, pt, candidates);
-        if (! candidates.empty())
-            for (size_t idx : candidates)
-                if (custom_data.polys[idx].contains(pt))
-                    return true;
-        return false;
-    };
+        //ranges:          [0 - 1]                                              (0 - 1.3]                               [0.1 - 1.1)
+        float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance)
+                * compute_angle_penalty(a.local_ccw_angle)
+                * distance_penalty_a;
+        float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance)
+                * compute_angle_penalty(b.local_ccw_angle)
+                * distance_penalty_b;
 
-    if (! m_enforcers[po_idx].empty()) {
-        const CustomTrianglesPerLayer& enforcers = m_enforcers[po_idx][layer_id];
-        if (! enforcers.polys.empty()) {
-            for (size_t i=0; i<polygon.points.size(); ++i) {
-                if (is_inside(polygon.points[i], enforcers))
-                    enforcers_idxs.emplace_back(i);
-            }
+        return penalty_a < penalty_b;
+    }
+
+    // Comparator used during alignment. If there is close potential aligned point, it is compared to the current
+    // seam point of the perimeter, to find out if the aligned point is not much worse than the current seam
+    // Also used by the random seam generator.
+    bool is_first_not_much_worse(const SeamCandidate &a, const SeamCandidate &b) const {
+        // Blockers/Enforcers discrimination, top priority
+        if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) {
+            // Prefer centers of enforcers.
+            return a.central_enforcer;
+        }
+
+        if (a.type == EnforcedBlockedSeamPoint::Enforced) {
+            return true;
+        }
+
+        if (a.type == EnforcedBlockedSeamPoint::Blocked) {
+            return false;
         }
+
+        if (a.type != b.type) {
+            return a.type > b.type;
+        }
+
+        //avoid overhangs
+        if (a.overhang > 0.0f || b.overhang > 0.0f) {
+            return a.overhang < b.overhang;
+        }
+
+        // prefer hidden points (more than 1 mm inside)
+        if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) {
+            return true;
+        }
+        if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) {
+            return false;
+        }
+
+        if (setup == SeamPosition::spRandom) {
+            return true;
+        }
+
+        if (setup == SeamPosition::spRear) {
+            return a.position.y() > b.position.y();
+        }
+
+        //ranges:          [0 - 1]                                          (0 - 1.3]                  ;
+        float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance)
+                * compute_angle_penalty(a.local_ccw_angle);
+        float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance)
+                * compute_angle_penalty(b.local_ccw_angle);
+
+        return penalty_a <= penalty_b || penalty_a - penalty_b < SeamPlacer::seam_align_score_tolerance;
     }
 
-    if (! m_blockers[po_idx].empty()) {
-        const CustomTrianglesPerLayer& blockers = m_blockers[po_idx][layer_id];
-        if (! blockers.polys.empty()) {
-            for (size_t i=0; i<polygon.points.size(); ++i) {
-                if (is_inside(polygon.points[i], blockers))
-                    blockers_idxs.emplace_back(i);
-            }
+    bool are_similar(const SeamCandidate &a, const SeamCandidate &b) const {
+        return is_first_not_much_worse(a, b) && is_first_not_much_worse(b, a);
+    }
+
+    float compute_angle_penalty(float ccw_angle) const {
+        // This function is used:
+        // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x)))
+        // looks scary, but it is gaussian combined with sigmoid,
+        // so that concave points have much smaller penalty over convex ones
+        // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png
+        return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) +
+                1.0f / (2 + std::exp(-ccw_angle)); // sigmoid, which heavily favourizes concave angles
+    }
+};
+
+#ifdef DEBUG_FILES
+void debug_export_points(const std::vector<PrintObjectSeamData::LayerSeams>  &layers,
+        const BoundingBox &bounding_box, std::string object_name, const SeamComparator &comparator) {
+    for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) {
+        std::string angles_file_name = debug_out_path(
+                (object_name + "_angles_" + std::to_string(layer_idx) + ".svg").c_str());
+        SVG angles_svg {angles_file_name, bounding_box};
+        float min_vis = 0;
+        float max_vis = min_vis;
+
+        float min_weight = std::numeric_limits<float>::min();
+        float max_weight = min_weight;
+
+        for (const SeamCandidate &point : layers[layer_idx].points) {
+            Vec3i color = value_to_rgbi(-PI, PI, point.local_ccw_angle);
+            std::string fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + ","
+            + std::to_string(color.z()) + ")";
+            angles_svg.draw(scaled(Vec2f(point.position.head<2>())), fill);
+            min_vis = std::min(min_vis, point.visibility);
+            max_vis = std::max(max_vis, point.visibility);
+
+            min_weight = std::min(min_weight, -comparator.compute_angle_penalty(point.local_ccw_angle));
+            max_weight = std::max(max_weight, -comparator.compute_angle_penalty(point.local_ccw_angle));
+
+        }
+
+        std::string visiblity_file_name = debug_out_path(
+                (object_name + "_visibility_" + std::to_string(layer_idx) + ".svg").c_str());
+        SVG visibility_svg {visiblity_file_name, bounding_box};
+        std::string weights_file_name = debug_out_path(
+                (object_name + "_weight_" + std::to_string(layer_idx) + ".svg").c_str());
+        SVG weight_svg {weights_file_name, bounding_box};
+        std::string overhangs_file_name = debug_out_path(
+                (object_name + "_overhang_" + std::to_string(layer_idx) + ".svg").c_str());
+        SVG overhangs_svg {overhangs_file_name, bounding_box};
+
+        for (const SeamCandidate &point : layers[layer_idx].points) {
+            Vec3i color = value_to_rgbi(min_vis, max_vis, point.visibility);
+            std::string visibility_fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + ","
+            + std::to_string(color.z()) + ")";
+            visibility_svg.draw(scaled(Vec2f(point.position.head<2>())), visibility_fill);
+
+            Vec3i weight_color = value_to_rgbi(min_weight, max_weight, -comparator.compute_angle_penalty(point.local_ccw_angle));
+            std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y())
+            + ","
+            + std::to_string(weight_color.z()) + ")";
+            weight_svg.draw(scaled(Vec2f(point.position.head<2>())), weight_fill);
+
+            Vec3i overhang_color = value_to_rgbi(-0.5, 0.5, std::clamp(point.overhang, -0.5f, 0.5f));
+            std::string overhang_fill = "rgb(" + std::to_string(overhang_color.x()) + ","
+            + std::to_string(overhang_color.y())
+            + ","
+            + std::to_string(overhang_color.z()) + ")";
+            overhangs_svg.draw(scaled(Vec2f(point.position.head<2>())), overhang_fill);
         }
     }
+}
+#endif
 
+// Pick best seam point based on the given comparator
+void pick_seam_point(std::vector<SeamCandidate> &perimeter_points, size_t start_index,
+        const SeamComparator &comparator) {
+    size_t end_index = perimeter_points[start_index].perimeter.end_index;
+
+    size_t seam_index = start_index;
+    for (size_t index = start_index; index <= end_index; ++index) {
+        if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index])) {
+            seam_index = index;
+        }
+    }
+    perimeter_points[start_index].perimeter.seam_index = seam_index;
 }
 
+size_t pick_nearest_seam_point_index(const std::vector<SeamCandidate> &perimeter_points, size_t start_index,
+        const Vec2f &preffered_location) {
+    size_t end_index = perimeter_points[start_index].perimeter.end_index;
+    SeamComparator comparator { spNearest };
 
-// Go through the polygon, identify points inside support enforcers and return
-// indices of points in the middle of each enforcer (measured along the contour).
-static std::vector<size_t> find_enforcer_centers(const Polygon& polygon,
-                                                 const std::vector<float>& lengths,
-                                                 const std::vector<size_t>& enforcers_idxs)
-{
-    std::vector<size_t> out;
-    assert(polygon.points.size()+1 == lengths.size());
-    assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end()));
-    if (polygon.size() < 2 || enforcers_idxs.empty())
-        return out;
-
-    auto get_center_idx = [&lengths](size_t start_idx, size_t end_idx) -> size_t {
-        assert(end_idx >= start_idx);
-        if (start_idx == end_idx)
-            return start_idx;
-        float t_c = lengths[start_idx] + 0.5f * (lengths[end_idx] - lengths[start_idx]);
-        auto it = std::lower_bound(lengths.begin() + start_idx, lengths.begin() + end_idx, t_c);
-        int ret = it - lengths.begin();
-        return ret;
-    };
+    size_t seam_index = start_index;
+    for (size_t index = start_index; index <= end_index; ++index) {
+        if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index], preffered_location)) {
+            seam_index = index;
+        }
+    }
+    return seam_index;
+}
 
-    int last_enforcer_start_idx = enforcers_idxs.front();
-    bool first_pt_in_list = enforcers_idxs.front() != 0;
-    bool last_pt_in_list = enforcers_idxs.back() == polygon.points.size() - 1;
-    bool wrap_around = last_pt_in_list && first_pt_in_list;
-
-    for (size_t i=0; i<enforcers_idxs.size(); ++i) {
-        if (i != enforcers_idxs.size() - 1) {
-            if (enforcers_idxs[i+1] != enforcers_idxs[i] + 1) {
-                // i is last point of current enforcer
-                out.push_back(get_center_idx(last_enforcer_start_idx, enforcers_idxs[i]));
-                last_enforcer_start_idx = enforcers_idxs[i+1];
-            }
+// picks random seam point uniformly, respecting enforcers blockers and overhang avoidance.
+void pick_random_seam_point(const std::vector<SeamCandidate> &perimeter_points, size_t start_index) {
+    SeamComparator comparator { spRandom };
+
+    // algorithm keeps a list of viable points and their lengths. If it finds a point
+    // that is much better than the viable_example_index (e.g. better type, no overhang; see is_first_not_much_worse)
+    // then it throws away stored lists and starts from start
+    // in the end, the list should contain points with same type (Enforced > Neutral > Blocked) and also only those which are not
+    // big overhang.
+    size_t viable_example_index = start_index;
+    size_t end_index = perimeter_points[start_index].perimeter.end_index;
+    struct Viable {
+        // Candidate seam point index.
+        size_t index;
+        float edge_length;
+        Vec3f edge;
+    };
+    std::vector<Viable> viables;
+
+    for (size_t index = start_index; index <= end_index; ++index) {
+        if (comparator.are_similar(perimeter_points[index], perimeter_points[viable_example_index])) {
+            // index ok, push info into viables
+            Vec3f edge_to_next { perimeter_points[index == end_index ? start_index : index + 1].position
+                    - perimeter_points[index].position };
+            float dist_to_next = edge_to_next.norm();
+            viables.push_back( { index, dist_to_next, edge_to_next });
+        } else if (comparator.is_first_not_much_worse(perimeter_points[viable_example_index],
+                perimeter_points[index])) {
+            // index is worse then viable_example_index, skip this point
         } else {
-            if (! wrap_around) {
-                // we can safely use the last enforcer point.
-                out.push_back(get_center_idx(last_enforcer_start_idx, enforcers_idxs[i]));
-            }
+            // index is better than viable example index, update example, clear gathered info, start again
+            // clear up all gathered info, start from scratch, update example index
+            viable_example_index = index;
+            viables.clear();
+
+            Vec3f edge_to_next = (perimeter_points[index == end_index ? start_index : index + 1].position
+                    - perimeter_points[index].position);
+            float dist_to_next = edge_to_next.norm();
+            viables.push_back( { index, dist_to_next, edge_to_next });
         }
     }
 
-    if (wrap_around) {
-        // Update first center already found.
-        if (out.empty()) {
-            // Probably an enforcer around the whole contour. Return nothing.
-            return out;
-        }
+    // now pick random point from the stored options
+    float len_sum = std::accumulate(viables.begin(), viables.end(), 0.0f, [](const float acc, const Viable &v) {
+        return acc + v.edge_length;
+    });
+    float picked_len = len_sum * (rand() / (float(RAND_MAX) + 1));
+
+    size_t point_idx = 0;
+    while (picked_len - viables[point_idx].edge_length > 0) {
+        picked_len = picked_len - viables[point_idx].edge_length;
+        point_idx++;
+    }
 
-        // find last point of the enforcer at the beginning:
-        size_t idx = 0;
-        while (enforcers_idxs[idx]+1 == enforcers_idxs[idx+1])
-            ++idx;
+    Perimeter &perimeter = perimeter_points[start_index].perimeter;
+    perimeter.seam_index = viables[point_idx].index;
+    perimeter.final_seam_position = perimeter_points[perimeter.seam_index].position
+            + viables[point_idx].edge.normalized() * picked_len;
+    perimeter.finalized = true;
+}
 
-        float t_s = lengths[last_enforcer_start_idx];
-        float t_e = lengths[idx];
-        float half_dist = 0.5f * (t_e + lengths.back() - t_s);
-        float t_c = (half_dist > t_e) ? t_s + half_dist : t_e - half_dist;
+struct EdgeGridWrapper {
+    explicit EdgeGridWrapper(ExPolygons ex_polys) :
+            ex_polys(ex_polys) {
 
-        auto it = std::lower_bound(lengths.begin(), lengths.end(), t_c);
-        out[0] = it - lengths.begin();
-        if (out[0] == lengths.size() - 1)
-            --out[0];
-        assert(out[0] < lengths.size() - 1);
+        grid.create(this->ex_polys, distance_field_resolution);
+        grid.calculate_sdf();
     }
-    return out;
+    const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5);
+    EdgeGrid::Grid grid;
+    ExPolygons ex_polys;
+}
+;
+
+EdgeGridWrapper compute_layer_merged_edge_grid(const Layer *layer) {
+    static const float eps = float(scale_(layer->object()->config().slice_closing_radius.value));
+    // merge with offset
+    ExPolygons merged = layer->merged(eps);
+    // ofsset back
+    ExPolygons layer_outline = offset_ex(merged, -eps);
+    return EdgeGridWrapper(layer_outline);
 }
 
+} // namespace SeamPlacerImpl
+
+// Parallel process and extract each perimeter polygon of the given print object.
+// Gather SeamCandidates of each layer into vector and build KDtree over them
+// Store results in the SeamPlacer variables m_seam_per_object
+void SeamPlacer::gather_seam_candidates(const PrintObject *po,
+        const SeamPlacerImpl::GlobalModelInfo &global_model_info, const SeamPosition configured_seam_preference) {
+    using namespace SeamPlacerImpl;
+
+    PrintObjectSeamData &seam_data = m_seam_per_object.emplace(po, PrintObjectSeamData { }).first->second;
+    seam_data.layers.resize(po->layer_count());
+
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, po->layers().size()),
+            [po, configured_seam_preference, &global_model_info, &seam_data](tbb::blocked_range<size_t> r) {
+                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
+                    PrintObjectSeamData::LayerSeams &layer_seams = seam_data.layers[layer_idx];
+                    const Layer *layer = po->get_layer(layer_idx);
+                    auto unscaled_z = layer->slice_z;
+                    std::vector<const LayerRegion*> regions;
+                    //NOTE corresponding region ptr may be null, if the layer has zero perimeters
+                    Polygons polygons = extract_perimeter_polygons(layer, configured_seam_preference, regions);
+                    for (size_t poly_index = 0; poly_index < polygons.size(); ++poly_index) {
+                        process_perimeter_polygon(polygons[poly_index], unscaled_z,
+                                regions[poly_index], global_model_info, layer_seams);
+                    }
+                    auto functor = SeamCandidateCoordinateFunctor { layer_seams.points };
+                    seam_data.layers[layer_idx].points_tree =
+                            std::make_unique<PrintObjectSeamData::SeamCandidatesTree>(functor,
+                                    layer_seams.points.size());
+                }
+            }
+    );
+}
 
+void SeamPlacer::calculate_candidates_visibility(const PrintObject *po,
+        const SeamPlacerImpl::GlobalModelInfo &global_model_info) {
+    using namespace SeamPlacerImpl;
+
+    std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
+            [&layers, &global_model_info](tbb::blocked_range<size_t> r) {
+                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
+                    for (auto &perimeter_point : layers[layer_idx].points) {
+                        perimeter_point.visibility = global_model_info.calculate_point_visibility(
+                                perimeter_point.position);
+                    }
+                }
+            });
+}
 
-void SeamPlacer::apply_custom_seam(const Polygon& polygon, size_t po_idx,
-                                   std::vector<float>& penalties,
-                                   const std::vector<float>& lengths,
-                                   int layer_id, SeamPosition seam_position) const
-{
-    if (! is_custom_seam_on_layer(layer_id, po_idx))
-        return;
+void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po) {
+    using namespace SeamPlacerImpl;
 
-    std::vector<size_t> enforcers_idxs;
-    std::vector<size_t> blockers_idxs;
-    this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs);
+    std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
+            [po, &layers](tbb::blocked_range<size_t> r) {
+                std::unique_ptr<EdgeGridWrapper> prev_layer_grid;
+                if (r.begin() > 0) { // previous layer exists
+                    prev_layer_grid = std::make_unique<EdgeGridWrapper>(
+                            compute_layer_merged_edge_grid(po->layers()[r.begin() - 1]));
+                }
 
-    for (size_t i : enforcers_idxs) {
-        assert(i < penalties.size());
-        penalties[i] -= float(ENFORCER_BLOCKER_PENALTY);
-    }
-    for (size_t i : blockers_idxs) {
-        assert(i < penalties.size());
-        penalties[i] += float(ENFORCER_BLOCKER_PENALTY);
+                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
+                    bool layer_has_multiple_loops =
+                            layers[layer_idx].points[0].perimeter.end_index
+                                    < layers[layer_idx].points.size() - 1;
+                    std::unique_ptr<EdgeGridWrapper> current_layer_grid = std::make_unique<EdgeGridWrapper>(
+                            compute_layer_merged_edge_grid(po->layers()[layer_idx]));
+
+                    for (SeamCandidate &perimeter_point : layers[layer_idx].points) {
+                        Point point = Point::new_scale(Vec2f { perimeter_point.position.head<2>() });
+                        if (prev_layer_grid.get() != nullptr) {
+                            coordf_t overhang_dist;
+                            prev_layer_grid->grid.signed_distance(point, scaled(perimeter_point.perimeter.flow_width),
+                                    overhang_dist);
+                            perimeter_point.overhang =
+                                    unscale<float>(overhang_dist) - perimeter_point.perimeter.flow_width;
+                        }
+
+                        if (layer_has_multiple_loops) { // search for embedded perimeter points (points hidden inside the print ,e.g. multimaterial join, best position for seam)
+                            coordf_t layer_embedded_distance;
+                            current_layer_grid->grid.signed_distance(point, scaled(1.0f),
+                                    layer_embedded_distance);
+                            perimeter_point.embedded_distance = unscale<float>(layer_embedded_distance);
+                        }
+                    }
+
+                    prev_layer_grid.swap(current_layer_grid);
+                }
+            }
+    );
+}
+
+// Estimates, if there is good seam point in the layer_idx which is close to last_point_pos
+// uses comparator.is_first_not_much_worse method to compare current seam with the closest point
+// (if current seam is too far away )
+// If the current chosen stream is close enough, it is stored in seam_string. returns true and updates last_point_pos
+// If the closest point is good enough to replace current chosen seam, it is stored in potential_string_seams, returns true and updates last_point_pos
+// Otherwise does nothing, returns false
+// Used by align_seam_points().
+bool SeamPlacer::find_next_seam_in_layer(
+        const std::vector<PrintObjectSeamData::LayerSeams> &layers,
+        std::pair<size_t, size_t> &last_point_indexes,
+        const size_t layer_idx, const float slice_z,
+        const SeamPlacerImpl::SeamComparator &comparator,
+        std::vector<std::pair<size_t, size_t>> &seam_string) const {
+    using namespace SeamPlacerImpl;
+
+    const SeamCandidate &last_point = layers[last_point_indexes.first].points[last_point_indexes.second];
+
+    Vec3f projected_position { last_point.position.x(), last_point.position.y(), slice_z };
+    std::vector<size_t> nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position,
+            SeamPlacer::seam_align_tolerable_dist);
+
+    if (nearby_points_indices.empty()) {
+        return false;
     }
-    if (seam_position == spAligned) {
-        std::vector<size_t> enf_centers = find_enforcer_centers(polygon, lengths, enforcers_idxs);
-        for (size_t idx : enf_centers) {
-            assert(idx < penalties.size());
-            penalties[idx] += ENFORCER_CENTER_PENALTY;
+
+    size_t best_nearby_point_index = nearby_points_indices[0];
+    size_t nearest_point_index = nearby_points_indices[0];
+
+    // Now find best nearby point, nearest point, and corresponding indices
+    for (const size_t &nearby_point_index : nearby_points_indices) {
+        const SeamCandidate &point = layers[layer_idx].points[nearby_point_index];
+        if (point.perimeter.finalized) {
+            continue; // skip over finalized perimeters, try to find some that is not finalized
+        }
+        if (comparator.is_first_better(point, layers[layer_idx].points[best_nearby_point_index],
+                projected_position.head<2>())
+                || layers[layer_idx].points[best_nearby_point_index].perimeter.finalized) {
+            best_nearby_point_index = nearby_point_index;
+        }
+        if ((point.position - projected_position).squaredNorm()
+                < (layers[layer_idx].points[nearest_point_index].position - projected_position).squaredNorm()
+                || layers[layer_idx].points[nearest_point_index].perimeter.finalized) {
+            nearest_point_index = nearby_point_index;
         }
     }
 
-#if 0
-    std::ostringstream os;
-    os << std::setw(3) << std::setfill('0') << layer_id;
-    int a = scale_(30.);
-    SVG svg("custom_seam" + os.str() + ".svg", BoundingBox(Point(-a, -a), Point(a, a)));
-    if (! m_enforcers[po_idx].empty())
-        svg.draw(m_enforcers[po_idx][layer_id].polys, "blue");
-    if (! m_blockers[po_idx].empty())
-        svg.draw(m_blockers[po_idx][layer_id].polys, "red");
+    const SeamCandidate &best_nearby_point = layers[layer_idx].points[best_nearby_point_index];
+    const SeamCandidate &nearest_point = layers[layer_idx].points[nearest_point_index];
 
-    if (! blockers_idxs.empty()) {
-        ;
+    if (nearest_point.perimeter.finalized) {
+        //all points are from already finalized perimeter, skip
+        return false;
     }
 
+    //from the nearest_point, deduce index of seam in the next layer
+    const SeamCandidate &next_layer_seam = layers[layer_idx].points[nearest_point.perimeter.seam_index];
 
-    size_t min_idx = std::min_element(penalties.begin(), penalties.end()) - penalties.begin();
+    // First try to pick central enforcer if any present
+    if (next_layer_seam.central_enforcer
+            && (next_layer_seam.position - projected_position).squaredNorm()
+                    < sqr(3 * SeamPlacer::seam_align_tolerable_dist)) {
+        last_point_indexes = std::pair<size_t, size_t> { layer_idx, nearest_point.perimeter.seam_index };
+        seam_string.push_back(last_point_indexes);
+        return true;
+    }
 
-    for (size_t i=0; i<polygon.points.size(); ++i) {
-        std::string fill;
-        coord_t size = 5e5;
-        if (min_idx == i)
-            fill = "yellow";
-        else
-            fill = (std::find(blockers_idxs.begin(), blockers_idxs.end(), i) != blockers_idxs.end() ? "green" : "black");
-        if (i != 0)
-            svg.draw(polygon.points[i], fill, size);
-        else
-            svg.draw(polygon.points[i], "red", 5e5);
+    // Next compare nearest and nearby point. If they are similar pick nearest, Otherwise expect curvy lines on smooth surfaces like chimney of benchy model
+    // We also compare it to the last point, to detect sharp changes in the scoring - that points to change in the model geometry and string should be ended.
+    if (comparator.are_similar(nearest_point, best_nearby_point)
+            && comparator.is_first_not_much_worse(nearest_point, next_layer_seam)
+            && comparator.are_similar(last_point, nearest_point)) {
+        last_point_indexes = std::pair<size_t, size_t> { layer_idx, nearest_point_index };
+        seam_string.push_back(last_point_indexes);
+        return true;
+    }
+    // If nearest point is not good enough, try it with the best nearby point.
+    if (comparator.is_first_not_much_worse(best_nearby_point, next_layer_seam)
+            && comparator.are_similar(last_point, nearest_point)) {
+        last_point_indexes = std::pair<size_t, size_t> { layer_idx, best_nearby_point_index };
+        seam_string.push_back(last_point_indexes);
+        return true;
     }
-#endif
 
+    return false;
 }
 
+// clusters already chosen seam points into strings across multiple layers, and then
+// aligns the strings via polynomial fit
+// Does not change the positions of the SeamCandidates themselves, instead stores
+// the new aligned position into the shared Perimeter structure of each perimeter
+// Note that this position does not necesarilly lay on the perimeter.
+void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) {
+    using namespace SeamPlacerImpl;
+
+    // Prepares Debug files for writing.
+#ifdef DEBUG_FILES
+    Slic3r::CNumericLocalesSetter locales_setter;
+    auto clusters_f = debug_out_path(("seam_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str());
+    FILE *clusters = boost::nowide::fopen(clusters_f.c_str(), "w");
+    if (clusters == nullptr) {
+        BOOST_LOG_TRIVIAL(error)
+        << "stl_write_obj: Couldn't open " << clusters_f << " for writing";
+        return;
+    }
+    auto aligned_f = debug_out_path(("aligned_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str());
+    FILE *aligns = boost::nowide::fopen(aligned_f.c_str(), "w");
+    if (aligns == nullptr) {
+        BOOST_LOG_TRIVIAL(error)
+        << "stl_write_obj: Couldn't open " << clusters_f << " for writing";
+        return;
+    }
+#endif
 
-
-std::optional<Point> SeamHistory::get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb)
-{
-    assert(layer_id >= m_layer_id || layer_id == 0);
-    if (layer_id != m_layer_id) {
-        // Get seam was called for different layer than last time.
-        if (layer_id == 0) // seq printing
-            m_data_this_layer.clear();
-        m_data_last_layer = m_data_this_layer;
-        m_data_this_layer.clear();
-        m_layer_id = layer_id;
+    //gather vector of all seams on the print_object - pair of layer_index and seam__index within that layer
+    const std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
+    std::vector<std::pair<size_t, size_t>> seams;
+    for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) {
+        const std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
+        size_t current_point_index = 0;
+        while (current_point_index < layer_perimeter_points.size()) {
+            seams.emplace_back(layer_idx, layer_perimeter_points[current_point_index].perimeter.seam_index);
+            current_point_index = layer_perimeter_points[current_point_index].perimeter.end_index + 1;
+        }
     }
 
-    std::optional<Point> out;
+    //sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice
+    std::sort(seams.begin(), seams.end(),
+            [&comparator, &layers](const std::pair<size_t, size_t> &left, const std::pair<size_t, size_t> &right) {
+                return comparator.is_first_better(layers[left.first].points[left.second],
+                        layers[right.first].points[right.second]);
+            }
+    );
+
+    //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment
+    // Keeping the vectors outside, so with a bit of luck they will not get reallocated after couple of for loop iterations.
+    std::vector<std::pair<size_t, size_t>> seam_string;
+    std::vector<Vec2f> observations;
+    std::vector<float> observation_points;
+    std::vector<float> weights;
+    for (const std::pair<size_t, size_t> &seam : seams) {
+        size_t layer_idx = seam.first;
+        size_t seam_index = seam.second;
+        const std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
+        if (layer_perimeter_points[seam_index].perimeter.finalized) {
+            // This perimeter is already aligned, skip seam
+            continue;
+        } else {
+
+            //initialize searching for seam string - cluster of nearby seams on previous and next layers
+            int skips = SeamPlacer::seam_align_tolerable_skips / 2;
+            int next_layer = layer_idx + 1;
+            std::pair<size_t, size_t> last_point_indexes = std::pair<size_t, size_t>(layer_idx, seam_index);
+
+            seam_string = { std::pair<size_t, size_t>(layer_idx, seam_index) };
+
+            //find seams or potential seams in forward direction; there is a budget of skips allowed
+            while (skips >= 0 && next_layer < int(layers.size())) {
+                if (find_next_seam_in_layer(layers, last_point_indexes, next_layer,
+                        float(po->get_layer(next_layer)->slice_z), comparator, seam_string)) {
+                    //String added, last_point_pos updated, nothing to be done
+                } else {
+                    // Layer skipped, reduce number of available skips
+                    skips--;
+                }
+                next_layer++;
+            }
 
-    auto seams_it = m_data_last_layer.find(po);
-    if (seams_it == m_data_last_layer.end())
-        return out;
+            //do additional check in back direction
+            next_layer = layer_idx - 1;
+            skips = SeamPlacer::seam_align_tolerable_skips / 2;
+            last_point_indexes = std::pair<size_t, size_t>(layer_idx, seam_index);
+            while (skips >= 0 && next_layer >= 0) {
+                if (find_next_seam_in_layer(layers, last_point_indexes, next_layer,
+                        float(po->get_layer(next_layer)->slice_z), comparator, seam_string)) {
+                    //String added, last_point_pos updated, nothing to be done
+                } else {
+                    // Layer skipped, reduce number of available skips
+                    skips--;
+                }
+                next_layer--;
+            }
 
-    const std::vector<SeamPoint>& seam_data_po = seams_it->second;
+            if (seam_string.size() < seam_align_minimum_string_seams) {
+                //string NOT long enough to be worth aligning, skip
+                continue;
+            }
 
-    // Find a bounding-box on the last layer that is close to one we see now.
-    double min_score = std::numeric_limits<double>::max();
-    for (const SeamPoint& sp : seam_data_po) {
-        const BoundingBox& bb = sp.m_island_bb;
+            // String is long enough, all string seams and potential string seams gathered, now do the alignment
+            //sort by layer index
+            std::sort(seam_string.begin(), seam_string.end(),
+                    [](const std::pair<size_t, size_t> &left, const std::pair<size_t, size_t> &right) {
+                        return left.first < right.first;
+                    });
+
+            // gather all positions of seams and their weights (weights are derived as negative penalty, they are made positive in next step)
+            observations.resize(seam_string.size());
+            observation_points.resize(seam_string.size());
+            weights.resize(seam_string.size());
+
+            //gather points positions and weights
+            // The algorithm uses only angle to compute penalty, to enforce snapping to sharp corners, if they are present
+            // after several experiments approach that gives best results is to snap the weight to one for sharp corners, and
+            //  leave it small for others. However, this can result in non-smooth line over area with a lot of unaligned sharp corners.
+            for (size_t index = 0; index < seam_string.size(); ++index) {
+                Vec3f pos = layers[seam_string[index].first].points[seam_string[index].second].position;
+                observations[index] = pos.head<2>();
+                observation_points[index] = pos.z();
+                weights[index] =
+                        (comparator.compute_angle_penalty(
+                                layers[seam_string[index].first].points[seam_string[index].second].local_ccw_angle)
+                                < comparator.compute_angle_penalty(0.4f * float(PI))) ? 1.0f : 0.1f;
+            }
 
-        if (! bb.overlap(island_bb)) {
-            // This bb does not even overlap. It is likely unrelated.
-            continue;
-        }
+            // Curve Fitting
+            size_t number_of_segments = std::max(size_t(1),
+                    size_t(observations.size() / SeamPlacer::seam_align_seams_per_segment));
+            auto curve = Geometry::fit_cubic_bspline(observations, observation_points, weights, number_of_segments);
+
+            // Do alignment - compute fitted point for each point in the string from its Z coord, and store the position into
+            // Perimeter structure of the point; also set flag aligned to true
+            for (size_t index = 0; index < seam_string.size(); ++index) {
+                const auto &pair = seam_string[index];
+                const float t = weights[index];
+                Vec3f current_pos = layers[pair.first].points[pair.second].position;
+                Vec2f fitted_pos = curve.get_fitted_value(current_pos.z());
+
+                //interpolate between current and fitted position, prefer current pos for large weights.
+                Vec3f final_position = t * current_pos + (1 - t) * to_3d(fitted_pos, current_pos.z());
+
+                Perimeter &perimeter = layers[pair.first].points[pair.second].perimeter;
+                perimeter.seam_index = pair.second;
+                perimeter.final_seam_position = final_position;
+                perimeter.finalized = true;
+            }
 
-        double score = std::pow(bb.min(0) - island_bb.min(0), 2.)
-                     + std::pow(bb.min(1) - island_bb.min(1), 2.)
-                     + std::pow(bb.max(0) - island_bb.max(0), 2.)
-                     + std::pow(bb.max(1) - island_bb.max(1), 2.);
+#ifdef DEBUG_FILES
+            auto randf = []() {
+                return float(rand()) / float(RAND_MAX);
+            };
+            Vec3f color { randf(), randf(), randf() };
+            for (size_t i = 0; i < seam_string.size(); ++i) {
+                auto orig_seam =  layers[seam_string[i].first].points[seam_string[i].second];
+                fprintf(clusters, "v %f %f %f %f %f %f \n", orig_seam.position[0],
+                        orig_seam.position[1],
+                        orig_seam.position[2], color[0], color[1],
+                        color[2]);
+            }
 
-        if (score < min_score) {
-            min_score = score;
-            out = sp.m_pos;
+            color = Vec3f { randf(), randf(), randf() };
+            for (size_t i = 0; i < seam_string.size(); ++i) {
+                const Perimeter &perimeter = layers[seam_string[i].first].points[seam_string[i].second].perimeter;
+                fprintf(aligns, "v %f %f %f %f %f %f \n", perimeter.final_seam_position[0],
+                        perimeter.final_seam_position[1],
+                        perimeter.final_seam_position[2], color[0], color[1],
+                        color[2]);
+            }
+#endif
         }
     }
 
-    return out;
+#ifdef DEBUG_FILES
+    fclose(clusters);
+    fclose(aligns);
+#endif
+
 }
 
+void SeamPlacer::init(const Print &print) {
+    using namespace SeamPlacerImpl;
+    m_seam_per_object.clear();
 
+    for (const PrintObject *po : print.objects()) {
 
-void SeamHistory::add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb)
-{
-    m_data_this_layer[po].push_back({pos, island_bb});;
-}
+        SeamPosition configured_seam_preference = po->config().seam_position.value;
+        SeamComparator comparator { configured_seam_preference };
+
+        GlobalModelInfo global_model_info { };
+        gather_enforcers_blockers(global_model_info, po);
+
+        if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) {
+            compute_global_occlusion(global_model_info, po);
+        }
+
+        BOOST_LOG_TRIVIAL(debug)
+        << "SeamPlacer: gather_seam_candidates: start";
+        gather_seam_candidates(po, global_model_info, configured_seam_preference);
+        BOOST_LOG_TRIVIAL(debug)
+        << "SeamPlacer: gather_seam_candidates: end";
+
+        if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) {
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: calculate_candidates_visibility : start";
+            calculate_candidates_visibility(po, global_model_info);
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: calculate_candidates_visibility : end";
+        }
 
+        BOOST_LOG_TRIVIAL(debug)
+        << "SeamPlacer: calculate_overhangs and layer embdedding : start";
+        calculate_overhangs_and_layer_embedding(po);
+        BOOST_LOG_TRIVIAL(debug)
+        << "SeamPlacer: calculate_overhangs and layer embdedding: end";
+
+        if (configured_seam_preference != spNearest) { // For spNearest, the seam is picked in the place_seam method with actual nozzle position information
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: pick_seam_point : start";
+            //pick seam point
+            std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
+                    [&layers, configured_seam_preference, comparator](tbb::blocked_range<size_t> r) {
+                        for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
+                            std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
+                            for (size_t current = 0; current < layer_perimeter_points.size();
+                                    current = layer_perimeter_points[current].perimeter.end_index + 1)
+                                if (configured_seam_preference == spRandom)
+                                    pick_random_seam_point(layer_perimeter_points, current);
+                                else
+                                    pick_seam_point(layer_perimeter_points, current, comparator);
+                        }
+                    });
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: pick_seam_point : end";
+        }
 
+        if (configured_seam_preference == spAligned) {
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: align_seam_points : start";
+            align_seam_points(po, comparator);
+            BOOST_LOG_TRIVIAL(debug)
+            << "SeamPlacer: align_seam_points : end";
+        }
 
-void SeamHistory::clear()
-{
-    m_layer_id = 0;
-    m_data_last_layer.clear();
-    m_data_this_layer.clear();
+#ifdef DEBUG_FILES
+        debug_export_points(layers, po->bounding_box(), std::to_string(po->id().id),
+                comparator);
+#endif
+    }
 }
 
+void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first,
+        const Point &last_pos) const {
+    using namespace SeamPlacerImpl;
+    const PrintObject *po = layer->object();
+    // Must not be called with supprot layer.
+    assert(dynamic_cast<const SupportLayer*>(layer) == nullptr);
+    // Object layer IDs are incremented by the number of raft layers.
+    assert(layer->id() >= po->slicing_parameters().raft_layers());
+    const size_t layer_index = layer->id() - po->slicing_parameters().raft_layers();
+    const double unscaled_z = layer->slice_z;
+
+    const PrintObjectSeamData::LayerSeams &layer_perimeters =
+            m_seam_per_object.find(layer->object())->second.layers[layer_index];
+
+    // Find the closest perimeter in the SeamPlacer to the first point of this loop.
+    size_t closest_perimeter_point_index;
+    {
+        const Point &fp = loop.first_point();
+        Vec2f unscaled_p = unscaled<float>(fp);
+        closest_perimeter_point_index = find_closest_point(*layer_perimeters.points_tree.get(),
+                to_3d(unscaled_p, float(unscaled_z)));
+    }
+
+    Vec3f seam_position;
+    size_t seam_index;
+    if (const Perimeter &perimeter = layer_perimeters.points[closest_perimeter_point_index].perimeter;
+    perimeter.finalized) {
+        seam_position = perimeter.final_seam_position;
+        seam_index = perimeter.seam_index;
+    } else {
+        seam_index =
+                po->config().seam_position == spNearest ?
+                        pick_nearest_seam_point_index(layer_perimeters.points, perimeter.start_index,
+                                unscaled<float>(last_pos)) :
+                        perimeter.seam_index;
+        seam_position = layer_perimeters.points[seam_index].position;
+    }
+
+    Point seam_point = Point::new_scale(seam_position.x(), seam_position.y());
+
+    if (const SeamCandidate &perimeter_point = layer_perimeters.points[seam_index];
+    (po->config().seam_position == spNearest || po->config().seam_position == spAligned) &&
+            loop.role() == ExtrusionRole::erPerimeter && //Hopefully internal perimeter
+            (seam_position - perimeter_point.position).squaredNorm() < 4.0f && // seam is on perimeter point
+            perimeter_point.local_ccw_angle < -EPSILON // In concave angles
+                    ) { // In this case, we are at internal perimeter, where the external perimeter has seam in concave angle. We want to align
+            // the internal seam into the concave corner, and not on the perpendicular projection on the closest edge (which is what the split_at function does)
+        size_t index_of_prev =
+                seam_index == perimeter_point.perimeter.start_index ?
+                                                                      perimeter_point.perimeter.end_index :
+                                                                      seam_index - 1;
+        size_t index_of_next =
+                seam_index == perimeter_point.perimeter.end_index ?
+                                                                    perimeter_point.perimeter.start_index :
+                                                                    seam_index + 1;
+
+        Vec2f dir_to_middle =
+                ((perimeter_point.position - layer_perimeters.points[index_of_prev].position).head<2>().normalized()
+                        + (perimeter_point.position - layer_perimeters.points[index_of_next].position).head<2>().normalized())
+                        * 0.5;
+
+        auto [_, projected_point] = loop.get_closest_path_and_point(seam_point, true);
+        //get closest projected point, determine depth of the seam point.
+        float depth = (float) unscale(Point(seam_point - projected_point)).norm();
+        float angle_factor = cos(-perimeter_point.local_ccw_angle / 2.0f); // There are some nice geometric identities in determination of the correct depth of new seam point.
+        //overshoot the target depth, in concave angles it will correctly snap to the corner; TODO: find out why such big overshoot is needed.
+        Vec2f final_pos = perimeter_point.position.head<2>() + (1.4142 * depth / angle_factor) * dir_to_middle;
+        seam_point = Point::new_scale(final_pos.x(), final_pos.y());
+    }
+
+    if (!loop.split_at_vertex(seam_point)) {
+        // The point is not in the original loop.
+        // Insert it.
+        loop.split_at(seam_point, true);
+    }
 
 }
+
+} // namespace Slic3r
diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp
index 0d72939b2..70881d558 100644
--- a/src/libslic3r/GCode/SeamPlacer.hpp
+++ b/src/libslic3r/GCode/SeamPlacer.hpp
@@ -3,12 +3,16 @@
 
 #include <optional>
 #include <vector>
+#include <memory>
+#include <atomic>
 
+#include "libslic3r/libslic3r.h"
 #include "libslic3r/ExtrusionEntity.hpp"
 #include "libslic3r/Polygon.hpp"
 #include "libslic3r/PrintConfig.hpp"
 #include "libslic3r/BoundingBox.hpp"
 #include "libslic3r/AABBTreeIndirect.hpp"
+#include "libslic3r/KDTreeIndirect.hpp"
 
 namespace Slic3r {
 
@@ -16,119 +20,148 @@ class PrintObject;
 class ExtrusionLoop;
 class Print;
 class Layer;
-namespace EdgeGrid { class Grid; }
 
+namespace EdgeGrid {
+class Grid;
+}
 
-class SeamHistory {
-public:
-    SeamHistory() { clear(); }
-    std::optional<Point> get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb);
-    void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb);
-    void clear();
+namespace SeamPlacerImpl {
 
-private:
-    struct SeamPoint {
-        Point m_pos;
-        BoundingBox m_island_bb;
-    };
+struct GlobalModelInfo;
+struct SeamComparator;
 
-    std::map<const PrintObject*, std::vector<SeamPoint>> m_data_last_layer;
-    std::map<const PrintObject*, std::vector<SeamPoint>> m_data_this_layer;
-    size_t m_layer_id;
+enum class EnforcedBlockedSeamPoint {
+    Blocked = 0,
+    Neutral = 1,
+    Enforced = 2,
 };
 
+// struct representing single perimeter loop
+struct Perimeter {
+    size_t start_index;
+    size_t end_index; //inclusive!
+    size_t seam_index;
+    float flow_width;
+
+    // During alignment, a final position may be stored here. In that case, finalized is set to true.
+    // Note that final seam position is not limited to points of the perimeter loop. In theory it can be any position
+    // Random position also uses this flexibility to set final seam point position
+    bool finalized = false;
+    Vec3f final_seam_position;
+};
 
+//Struct over which all processing of perimeters is done. For each perimeter point, its respective candidate is created,
+// then all the needed attributes are computed and finally, for each perimeter one point is chosen as seam.
+// This seam position can be then further aligned
+struct SeamCandidate {
+    SeamCandidate(const Vec3f &pos, Perimeter &perimeter,
+            float local_ccw_angle,
+            EnforcedBlockedSeamPoint type) :
+            position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), embedded_distance(0.0f), local_ccw_angle(
+                    local_ccw_angle), type(type), central_enforcer(false) {
+    }
+    const Vec3f position;
+    // pointer to Perimeter loop of this point. It is shared across all points of the loop
+    Perimeter &perimeter;
+    float visibility;
+    float overhang;
+    // distance inside the merged layer regions, for detecting perimeter points which are hidden indside the print (e.g. multimaterial join)
+    // Negative sign means inside the print, comes from EdgeGrid structure
+    float embedded_distance;
+    float local_ccw_angle;
+    EnforcedBlockedSeamPoint type;
+    bool central_enforcer; //marks this candidate as central point of enforced segment on the perimeter - important for alignment
+};
 
-class SeamPlacer {
-public:
-    void init(const Print& print);
-
-    // When perimeters are printed, first call this function with the respective
-    // external perimeter. SeamPlacer will find a location for its seam and remember it.
-    // Subsequent calls to get_seam will return this position.
-
-
-    void plan_perimeters(const std::vector<const ExtrusionEntity*> perimeters,
-        const Layer& layer, SeamPosition seam_position,
-        Point last_pos, coordf_t nozzle_dmr, const PrintObject* po,
-        const EdgeGrid::Grid* lower_layer_edge_grid);
-
-    void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
-                    const EdgeGrid::Grid* lower_layer_edge_grid);
-    
-
-    using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
-    using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
-
-private:
+struct FaceVisibilityInfo {
+    float visibility;
+};
 
-    // When given an external perimeter (!), returns the seam.
-    Point calculate_seam(const Layer& layer, const SeamPosition seam_position,
-        const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po,
-        const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest);
+struct SeamCandidateCoordinateFunctor {
+    SeamCandidateCoordinateFunctor(const std::vector<SeamCandidate> &seam_candidates) :
+            seam_candidates(seam_candidates) {
+    }
+    const std::vector<SeamCandidate> &seam_candidates;
+    float operator()(size_t index, size_t dim) const {
+        return seam_candidates[index].position[dim];
+    }
+};
+} // namespace SeamPlacerImpl
 
-    struct CustomTrianglesPerLayer {
-        Polygons polys;
-        TreeType tree;
-    };
+struct PrintObjectSeamData
+{
+    using SeamCandidatesTree = KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>;
 
-    // Just a cache to save some lookups.
-    const Layer* m_last_layer_po = nullptr;
-    coordf_t m_last_print_z = -1.;
-    const PrintObject* m_last_po = nullptr;
-
-    struct SeamPoint {
-        Point pt;
-        bool precalculated = false;
-        bool external = false;
-        const Layer* layer = nullptr;
-        SeamPosition seam_position;
-        const PrintObject* po = nullptr;
+    struct LayerSeams
+    {
+        Slic3r::deque<SeamPlacerImpl::Perimeter>    perimeters;
+        std::vector<SeamPlacerImpl::SeamCandidate>  points;
+        std::unique_ptr<SeamCandidatesTree>         points_tree;
     };
-    std::vector<SeamPoint> m_plan;
-    size_t m_plan_idx;
-
-    std::vector<std::vector<CustomTrianglesPerLayer>> m_enforcers;
-    std::vector<std::vector<CustomTrianglesPerLayer>> m_blockers;
-    std::vector<const PrintObject*> m_po_list;
-
-    //std::map<const PrintObject*, Point>  m_last_seam_position;
-    SeamHistory  m_seam_history;
-
-    // Get indices of points inside enforcers and blockers.
-    void get_enforcers_and_blockers(size_t layer_id,
-                                    const Polygon& polygon,
-                                    size_t po_id,
-                                    std::vector<size_t>& enforcers_idxs,
-                                    std::vector<size_t>& blockers_idxs) const;
-
-    // Apply penalties to points inside enforcers/blockers.
-    void apply_custom_seam(const Polygon& polygon, size_t po_id,
-                           std::vector<float>& penalties,
-                           const std::vector<float>& lengths,
-                           int layer_id, SeamPosition seam_position) const;
-
-    // Return random point of a polygon. The distribution will be uniform
-    // along the contour and account for enforcers and blockers.
-    Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id,
-                          bool* saw_custom = nullptr) const;
-
-    // Is there any enforcer/blocker on this layer?
-    bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const {
-        return is_custom_enforcer_on_layer(layer_id, po_idx)
-            || is_custom_blocker_on_layer(layer_id, po_idx);
+    // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter
+    std::vector<LayerSeams> layers;
+    // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD
+    // tree of all points of the given layer
+
+    void clear()
+    {
+        layers.clear();
     }
+};
 
-    bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty());
-    }
+class SeamPlacer {
+public:
+    static constexpr size_t raycasting_decimation_target_triangle_count = 10000;
+    static constexpr float raycasting_subdivision_target_length = 2.0f;
+    //square of number of rays per triangle
+    static constexpr size_t sqr_rays_per_triangle = 7;
+
+    // arm length used during angles computation
+    static constexpr float polygon_local_angles_arm_distance = 0.5f;
+
+    // increases angle importance at the cost of deacreasing visibility info importance. must be > 0
+    static constexpr float additional_angle_importance = 0.6f;
+
+    // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer
+    static constexpr float enforcer_blocker_distance_tolerance = 0.35f;
+    // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size
+    static constexpr float enforcer_oversampling_distance = 0.2f;
+
+    // When searching for seam clusters for alignment:
+    // following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer
+    static constexpr float seam_align_score_tolerance = 0.5f;
+    // seam_align_tolerable_dist - if next layer closes point is too far away, break string
+    static constexpr float seam_align_tolerable_dist = 1.0f;
+    // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped.
+    // this param limits the number of allowed skips
+    static constexpr size_t seam_align_tolerable_skips = 4;
+    // minimum number of seams needed in cluster to make alignment happen
+    static constexpr size_t seam_align_minimum_string_seams = 6;
+    // points covered by spline; determines number of splines for the given string
+    static constexpr size_t seam_align_seams_per_segment = 8;
+
+    //The following data structures hold all perimeter points for all PrintObject.
+    std::unordered_map<const PrintObject*, PrintObjectSeamData> m_seam_per_object;
+
+    void init(const Print &print);
+
+    void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const;
 
-    bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty());
-    }
+private:
+    void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info,
+            const SeamPosition configured_seam_preference);
+    void calculate_candidates_visibility(const PrintObject *po,
+            const SeamPlacerImpl::GlobalModelInfo &global_model_info);
+    void calculate_overhangs_and_layer_embedding(const PrintObject *po);
+    void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator);
+    bool find_next_seam_in_layer(
+            const std::vector<PrintObjectSeamData::LayerSeams> &layers,
+            std::pair<size_t, size_t> &last_point_indexes,
+            const size_t layer_idx, const float slice_z,
+            const SeamPlacerImpl::SeamComparator &comparator,
+            std::vector<std::pair<size_t, size_t>> &seam_string) const;
 };
 
-
-}
+} // namespace Slic3r
 
 #endif // libslic3r_SeamPlacer_hpp_
diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp
index c825f8254..82ffbd8d1 100644
--- a/src/libslic3r/Geometry.hpp
+++ b/src/libslic3r/Geometry.hpp
@@ -36,9 +36,9 @@ enum Orientation
 static inline Orientation orient(const Point &a, const Point &b, const Point &c)
 {
     static_assert(sizeof(coord_t) * 2 == sizeof(int64_t), "orient works with 32 bit coordinates");
-    int64_t u = int64_t(b(0)) * int64_t(c(1)) - int64_t(b(1)) * int64_t(c(0));
-    int64_t v = int64_t(a(0)) * int64_t(c(1)) - int64_t(a(1)) * int64_t(c(0));
-    int64_t w = int64_t(a(0)) * int64_t(b(1)) - int64_t(a(1)) * int64_t(b(0));
+    int64_t u = int64_t(b.x()) * int64_t(c.y()) - int64_t(b.y()) * int64_t(c.x());
+    int64_t v = int64_t(a.x()) * int64_t(c.y()) - int64_t(a.y()) * int64_t(c.x());
+    int64_t w = int64_t(a.x()) * int64_t(b.y()) - int64_t(a.y()) * int64_t(b.x());
     int64_t d = u - v + w;
     return (d > 0) ? ORIENTATION_CCW : ((d == 0) ? ORIENTATION_COLINEAR : ORIENTATION_CW);
 }
diff --git a/src/libslic3r/Geometry/Bicubic.hpp b/src/libslic3r/Geometry/Bicubic.hpp
new file mode 100644
index 000000000..98c6f8bb2
--- /dev/null
+++ b/src/libslic3r/Geometry/Bicubic.hpp
@@ -0,0 +1,291 @@
+#ifndef BICUBIC_HPP
+#define BICUBIC_HPP
+
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+#include <Eigen/Dense>
+
+namespace Slic3r {
+
+namespace Geometry {
+
+namespace BicubicInternal {
+// Linear kernel, to be able to test cubic methods with hat kernels.
+template<typename T>
+struct LinearKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return T(0.);
+    }
+    static T a01() {
+        return T(0.);
+    }
+    static T a02() {
+        return T(0.);
+    }
+    static T a03() {
+        return T(0.);
+    }
+    static T a10() {
+        return T(1.);
+    }
+    static T a11() {
+        return T(-1.);
+    }
+    static T a12() {
+        return T(0.);
+    }
+    static T a13() {
+        return T(0.);
+    }
+    static T a20() {
+        return T(0.);
+    }
+    static T a21() {
+        return T(1.);
+    }
+    static T a22() {
+        return T(0.);
+    }
+    static T a23() {
+        return T(0.);
+    }
+    static T a30() {
+        return T(0.);
+    }
+    static T a31() {
+        return T(0.);
+    }
+    static T a32() {
+        return T(0.);
+    }
+    static T a33() {
+        return T(0.);
+    }
+};
+
+// Interpolation kernel aka Catmul-Rom aka Keyes kernel.
+template<typename T>
+struct CubicCatmulRomKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return 0;
+    }
+    static T a01() {
+        return T( -0.5);
+    }
+    static T a02() {
+        return T( 1.);
+    }
+    static T a03() {
+        return T( -0.5);
+    }
+    static T a10() {
+        return T( 1.);
+    }
+    static T a11() {
+        return 0;
+    }
+    static T a12() {
+        return T( -5. / 2.);
+    }
+    static T a13() {
+        return T( 3. / 2.);
+    }
+    static T a20() {
+        return 0;
+    }
+    static T a21() {
+        return T( 0.5);
+    }
+    static T a22() {
+        return T( 2.);
+    }
+    static T a23() {
+        return T( -3. / 2.);
+    }
+    static T a30() {
+        return 0;
+    }
+    static T a31() {
+        return 0;
+    }
+    static T a32() {
+        return T( -0.5);
+    }
+    static T a33() {
+        return T( 0.5);
+    }
+};
+
+// B-spline kernel
+template<typename T>
+struct CubicBSplineKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return T( 1. / 6.);
+    }
+    static T a01() {
+        return T( -3. / 6.);
+    }
+    static T a02() {
+        return T( 3. / 6.);
+    }
+    static T a03() {
+        return T( -1. / 6.);
+    }
+    static T a10() {
+        return T( 4. / 6.);
+    }
+    static T a11() {
+        return 0;
+    }
+    static T a12() {
+        return T( -6. / 6.);
+    }
+    static T a13() {
+        return T( 3. / 6.);
+    }
+    static T a20() {
+        return T( 1. / 6.);
+    }
+    static T a21() {
+        return T( 3. / 6.);
+    }
+    static T a22() {
+        return T( 3. / 6.);
+    }
+    static T a23() {
+        return T( -3. / 6.);
+    }
+    static T a30() {
+        return 0;
+    }
+    static T a31() {
+        return 0;
+    }
+    static T a32() {
+        return 0;
+    }
+    static T a33() {
+        return T( 1. / 6.);
+    }
+};
+
+template<class T>
+inline T clamp(T a, T lower, T upper)
+        {
+    return (a < lower) ? lower :
+           (a > upper) ? upper : a;
+}
+}
+
+template<typename Kernel>
+struct CubicKernelWrapper
+{
+    typedef typename Kernel::FloatType FloatType;
+
+    static constexpr size_t kernel_span = 4;
+
+    static FloatType kernel(FloatType x)
+            {
+        x = fabs(x);
+        if (x >= (FloatType) 2.)
+            return 0.0f;
+        if (x <= (FloatType) 1.) {
+            FloatType x2 = x * x;
+            FloatType x3 = x2 * x;
+            return Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3;
+        }
+        assert(x > (FloatType )1. && x < (FloatType )2.);
+        x -= (FloatType) 1.;
+        FloatType x2 = x * x;
+        FloatType x3 = x2 * x;
+        return Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3;
+    }
+
+    static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x)
+            {
+        const FloatType x2 = x * x;
+        const FloatType x3 = x * x * x;
+        return f0 * (Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3) +
+                f1 * (Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3) +
+                f2 * (Kernel::a20() + Kernel::a21() * x + Kernel::a22() * x2 + Kernel::a23() * x3) +
+                f3 * (Kernel::a30() + Kernel::a31() * x + Kernel::a32() * x2 + Kernel::a33() * x3);
+    }
+};
+
+// Linear splines
+template<typename NumberType>
+using LinearKernel = CubicKernelWrapper<BicubicInternal::LinearKernel<NumberType>>;
+
+// Catmul-Rom splines
+template<typename NumberType>
+using CubicCatmulRomKernel = CubicKernelWrapper<BicubicInternal::CubicCatmulRomKernel<NumberType>>;
+
+// Cubic B-splines
+template<typename NumberType>
+using CubicBSplineKernel = CubicKernelWrapper<BicubicInternal::CubicBSplineKernel<NumberType>>;
+
+template<typename KernelWrapper>
+static typename KernelWrapper::FloatType cubic_interpolate(const Eigen::ArrayBase<typename KernelWrapper::FloatType> &F,
+        const typename KernelWrapper::FloatType pt) {
+    typedef typename KernelWrapper::FloatType T;
+    const int w = int(F.size());
+    const int ix = (int) floor(pt);
+    const T s = pt - T( ix);
+
+    if (ix > 1 && ix + 2 < w) {
+        // Inside the fully interpolated region.
+        return KernelWrapper::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
+    }
+    // Transition region. Extend with a constant function.
+    auto f = [&F, w](T x) {
+        return F[BicubicInternal::clamp(x, 0, w - 1)];
+    };
+    return KernelWrapper::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s);
+}
+
+template<typename Kernel, typename Derived>
+static float bicubic_interpolate(const Eigen::MatrixBase<Derived> &F,
+        const Eigen::Matrix<typename Kernel::FloatType, 2, 1, Eigen::DontAlign> &pt) {
+    typedef typename Kernel::FloatType T;
+    const int w = F.cols();
+    const int h = F.rows();
+    const int ix = (int) floor(pt[0]);
+    const int iy = (int) floor(pt[1]);
+    const T s = pt[0] - T( ix);
+    const T t = pt[1] - T( iy);
+
+    if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) {
+        // Inside the fully interpolated region.
+        return Kernel::interpolate(
+                Kernel::interpolate(F(ix - 1, iy - 1), F(ix, iy - 1), F(ix + 1, iy - 1), F(ix + 2, iy - 1), s),
+                Kernel::interpolate(F(ix - 1, iy), F(ix, iy), F(ix + 1, iy), F(ix + 2, iy), s),
+                Kernel::interpolate(F(ix - 1, iy + 1), F(ix, iy + 1), F(ix + 1, iy + 1), F(ix + 2, iy + 1), s),
+                Kernel::interpolate(F(ix - 1, iy + 2), F(ix, iy + 2), F(ix + 1, iy + 2), F(ix + 2, iy + 2), s), t);
+    }
+    // Transition region. Extend with a constant function.
+    auto f = [&F, w, h](int x, int y) {
+        return F(BicubicInternal::clamp(x, 0, w - 1), BicubicInternal::clamp(y, 0, h - 1));
+    };
+    return Kernel::interpolate(
+            Kernel::interpolate(f(ix - 1, iy - 1), f(ix, iy - 1), f(ix + 1, iy - 1), f(ix + 2, iy - 1), s),
+            Kernel::interpolate(f(ix - 1, iy), f(ix, iy), f(ix + 1, iy), f(ix + 2, iy), s),
+            Kernel::interpolate(f(ix - 1, iy + 1), f(ix, iy + 1), f(ix + 1, iy + 1), f(ix + 2, iy + 1), s),
+            Kernel::interpolate(f(ix - 1, iy + 2), f(ix, iy + 2), f(ix + 1, iy + 2), f(ix + 2, iy + 2), s), t);
+}
+
+} //namespace Geometry
+
+} // namespace Slic3r
+
+#endif /* BICUBIC_HPP */
diff --git a/src/libslic3r/Geometry/ConvexHull.cpp b/src/libslic3r/Geometry/ConvexHull.cpp
index b1ff77f80..2e92535f2 100644
--- a/src/libslic3r/Geometry/ConvexHull.cpp
+++ b/src/libslic3r/Geometry/ConvexHull.cpp
@@ -1,6 +1,7 @@
 #include "libslic3r.h"
 #include "ConvexHull.hpp"
 #include "BoundingBox.hpp"
+#include "../Geometry.hpp"
 
 #include <boost/multiprecision/integer.hpp>
 
@@ -19,13 +20,13 @@ Polygon convex_hull(Points pts)
         hull.points.resize(2 * n);
         // Build lower hull
         for (int i = 0; i < n; ++ i) {
-            while (k >= 2 && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
+            while (k >= 2 && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
                 -- k;
             hull[k ++] = pts[i];
         }
         // Build upper hull
         for (int i = n-2, t = k+1; i >= 0; i--) {
-            while (k >= t && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
+            while (k >= t && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
                 -- k;
             hull[k ++] = pts[i];
         }
@@ -58,7 +59,7 @@ Pointf3s convex_hull(Pointf3s points)
                 Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
                 Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
 
-                if (p.ccw(k2, k1) <= 0)
+                if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
                     --k;
                 else
                     break;
@@ -76,7 +77,7 @@ Pointf3s convex_hull(Pointf3s points)
                 Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
                 Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
 
-                if (p.ccw(k2, k1) <= 0)
+                if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
                     --k;
                 else
                     break;
diff --git a/src/libslic3r/Geometry/Curves.hpp b/src/libslic3r/Geometry/Curves.hpp
new file mode 100644
index 000000000..6542cb706
--- /dev/null
+++ b/src/libslic3r/Geometry/Curves.hpp
@@ -0,0 +1,205 @@
+#ifndef SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
+#define SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
+
+#include "libslic3r/Point.hpp"
+#include "Bicubic.hpp"
+
+#include <iostream>
+
+//#define LSQR_DEBUG
+
+namespace Slic3r {
+namespace Geometry {
+
+template<int Dimension, typename NumberType>
+struct PolynomialCurve {
+    Eigen::MatrixXf coefficients;
+
+    Vec3f get_fitted_value(const NumberType value) const {
+        auto result = Vec<Dimension, NumberType>::Zero();
+        size_t order = this->coefficients.rows() - 1;
+        auto x = NumberType(1.);
+        for (size_t index = 0; index < order + 1; ++index, x *= value)
+            result += x * this->coefficients.col(index);
+        return result;
+    }
+};
+
+//https://towardsdatascience.com/least-square-polynomial-CURVES-using-c-eigen-package-c0673728bd01
+template<int Dimension, typename NumberType>
+PolynomialCurve<Dimension, NumberType> fit_polynomial(const std::vector<Vec<Dimension, NumberType>> &observations,
+        const std::vector<NumberType> &observation_points,
+        const std::vector<NumberType> &weights, size_t order) {
+    // check to make sure inputs are correct
+    size_t cols = order + 1;
+    assert(observation_points.size() >= cols);
+    assert(observation_points.size() == weights.size());
+    assert(observations.size() == weights.size());
+
+    Eigen::MatrixXf data_points(Dimension, observations.size());
+    Eigen::MatrixXf T(observations.size(), cols);
+    for (size_t i = 0; i < weights.size(); ++i) {
+        auto squared_weight = sqrt(weights[i]);
+        data_points.col(i) = observations[i] * squared_weight;
+        // Populate the matrix
+        auto x = squared_weight;
+        auto c = observation_points[i];
+        for (size_t j = 0; j < cols; ++j, x *= c)
+            T(i, j) = x;
+    }
+
+    const auto QR = T.householderQr();
+    Eigen::MatrixXf coefficients(Dimension, cols);
+    // Solve for linear least square fit
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+        coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
+    }
+
+    return {std::move(coefficients)};
+}
+
+template<size_t Dimension, typename NumberType, typename KernelType>
+struct PiecewiseFittedCurve {
+    using Kernel = KernelType;
+
+    Eigen::MatrixXf coefficients;
+    NumberType start;
+    NumberType segment_size;
+    size_t endpoints_level_of_freedom;
+
+    Vec<Dimension, NumberType> get_fitted_value(const NumberType &observation_point) const {
+        Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
+
+        //find corresponding segment index; expects kernels to be centered
+        int middle_right_segment_index = floor((observation_point - start) / segment_size);
+        //find index of first segment that is affected by the point i; this can be deduced from kernel_span
+        int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1;
+        for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
+                segment_index++) {
+            NumberType segment_start = start + segment_index * segment_size;
+            NumberType normalized_segment_distance = (segment_start - observation_point) / segment_size;
+
+            int parameter_index = segment_index + endpoints_level_of_freedom;
+            parameter_index = std::clamp(parameter_index, 0, int(coefficients.cols()) - 1);
+            result += Kernel::kernel(normalized_segment_distance) * coefficients.col(parameter_index);
+        }
+        return result;
+    }
+};
+
+// observations: data to be fitted by the curve
+// observation points: growing sequence of points where the observations were made.
+//      In other words, for function f(x) = y, observations are y0...yn, and observation points are x0...xn
+// weights: how important the observation is
+// segments_count: number of segments inside the valid length of the curve
+// endpoints_level_of_freedom: number of additional parameters at each end; reasonable values depend on the kernel span
+template<typename Kernel, int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, Kernel> fit_curve(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        const std::vector<NumberType> &observation_points,
+        const std::vector<NumberType> &weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom) {
+
+    // check to make sure inputs are correct
+    assert(segments_count > 0);
+    assert(observations.size() > 0);
+    assert(observation_points.size() == observations.size());
+    assert(observation_points.size() == weights.size());
+    assert(segments_count <= observations.size());
+
+    //prepare sqrt of weights, which will then be applied to both matrix T and observed data: https://en.wikipedia.org/wiki/Weighted_least_squares
+    std::vector<NumberType> sqrt_weights(weights.size());
+    for (size_t index = 0; index < weights.size(); ++index) {
+        assert(weights[index] > 0);
+        sqrt_weights[index] = sqrt(weights[index]);
+    }
+
+    // prepare result and compute metadata
+    PiecewiseFittedCurve<Dimension, NumberType, Kernel> result { };
+
+    NumberType valid_length = observation_points.back() - observation_points.front();
+    NumberType segment_size = valid_length / NumberType(segments_count);
+    result.start = observation_points.front();
+    result.segment_size = segment_size;
+    result.endpoints_level_of_freedom = endpoints_level_of_freedom;
+
+    // prepare observed data
+    // Eigen defaults to column major memory layout.
+    Eigen::MatrixXf data_points(Dimension, observations.size());
+    for (size_t index = 0; index < observations.size(); ++index) {
+        data_points.col(index) = observations[index] * sqrt_weights[index];
+    }
+    // parameters count is always increased by one to make the parametric space of the curve symmetric.
+    // without this fix, the end of the curve is less flexible than the beginning
+    size_t parameters_count = segments_count + 1 + 2 * endpoints_level_of_freedom;
+    //Create weight matrix T for each point and each segment;
+    Eigen::MatrixXf T(observation_points.size(), parameters_count);
+    T.setZero();
+    //Fill the weight matrix
+    for (size_t i = 0; i < observation_points.size(); ++i) {
+        NumberType observation_point = observation_points[i];
+        //find corresponding segment index; expects kernels to be centered
+        int middle_right_segment_index = floor((observation_point - result.start) / result.segment_size);
+        //find index of first segment that is affected by the point i; this can be deduced from kernel_span
+        int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1;
+        for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
+                segment_index++) {
+            NumberType segment_start = result.start + segment_index * result.segment_size;
+            NumberType normalized_segment_distance = (segment_start - observation_point) / result.segment_size;
+
+            int parameter_index = segment_index + endpoints_level_of_freedom;
+            parameter_index = std::clamp(parameter_index, 0, int(parameters_count) - 1);
+            T(i, parameter_index) += Kernel::kernel(normalized_segment_distance) * sqrt_weights[i];
+        }
+    }
+
+#ifdef LSQR_DEBUG
+    std::cout << "weight matrix: " << std::endl;
+    for (int obs = 0; obs < observation_points.size(); ++obs) {
+        std::cout << std::endl;
+        for (int segment = 0; segment < parameters_count; ++segment) {
+            std::cout << T(obs, segment) << "  ";
+        }
+    }
+    std::cout << std::endl;
+#endif
+
+    // Solve for linear least square fit
+    result.coefficients.resize(Dimension, parameters_count);
+    const auto QR = T.fullPivHouseholderQr();
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+        result.coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
+    }
+
+    return result;
+}
+
+template<int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, CubicBSplineKernel<NumberType>>
+fit_cubic_bspline(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        std::vector<NumberType> observation_points,
+        std::vector<NumberType> weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom = 0) {
+    return fit_curve<CubicBSplineKernel<NumberType>>(observations, observation_points, weights, segments_count,
+            endpoints_level_of_freedom);
+}
+
+template<int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, CubicCatmulRomKernel<NumberType>>
+fit_catmul_rom_spline(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        std::vector<NumberType> observation_points,
+        std::vector<NumberType> weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom = 0) {
+    return fit_curve<CubicCatmulRomKernel<NumberType>>(observations, observation_points, weights, segments_count,
+            endpoints_level_of_freedom);
+}
+
+}
+}
+
+#endif /* SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ */
diff --git a/src/libslic3r/KDTreeIndirect.hpp b/src/libslic3r/KDTreeIndirect.hpp
index 12e462569..36a157456 100644
--- a/src/libslic3r/KDTreeIndirect.hpp
+++ b/src/libslic3r/KDTreeIndirect.hpp
@@ -59,7 +59,7 @@ public:
 		CONTINUE_RIGHT  = 2,
 		STOP 			= 4,
 	};
-	template<typename CoordType> 
+	template<typename CoordType>
 	unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
 	{
 		CoordType dist = point_coord - this->coordinate(idx, dimension);
@@ -110,8 +110,8 @@ private:
 
 	// Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
 	// https://en.wikipedia.org/wiki/Quickselect
-	// Items left of the k'th item are lower than the k'th item in the "dimension", 
-	// items right of the k'th item are higher than the k'th item in the "dimension", 
+	// Items left of the k'th item are lower than the k'th item in the "dimension",
+	// items right of the k'th item are higher than the k'th item in the "dimension",
 	void partition_input(std::vector<size_t> &input, const size_t dimension, size_t left, size_t right, const size_t k) const
 	{
 		while (left < right) {
@@ -231,6 +231,53 @@ size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& poi
 	return find_closest_point(kdtree, point, [](size_t) { return true; });
 }
 
+// Find nearby points (spherical neighbourhood) using Euclidian metrics.
+template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+       const typename KDTreeIndirectType::CoordType& max_distance, FilterFn filter)
+        {
+    using CoordType = typename KDTreeIndirectType::CoordType;
+
+    struct Visitor {
+        const KDTreeIndirectType &kdtree;
+        const PointType center;
+        const CoordType max_distance_squared;
+        const FilterFn filter;
+        std::vector<size_t> result;
+
+        Visitor(const KDTreeIndirectType &kdtree, const PointType& center, const CoordType &max_distance,
+                FilterFn filter) :
+                kdtree(kdtree), center(center), max_distance_squared(max_distance*max_distance), filter(filter) {
+        }
+        unsigned int operator()(size_t idx, size_t dimension) {
+            if (this->filter(idx)) {
+                auto dist = CoordType(0);
+                for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
+                    CoordType d = center[i] - kdtree.coordinate(idx, i);
+                    dist += d * d;
+                }
+                if (dist < max_distance_squared) {
+                    result.push_back(idx);
+                }
+            }
+           return kdtree.descent_mask(center[dimension], max_distance_squared, idx, dimension);
+        }
+    } visitor(kdtree, center, max_distance, filter);
+
+    kdtree.visit(visitor);
+    return visitor.result;
+}
+
+template<typename KDTreeIndirectType, typename PointType>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+       const typename KDTreeIndirectType::CoordType& max_distance)
+        {
+    return find_nearby_points(kdtree, center, max_distance, [](size_t) {
+        return true;
+    });
+}
+
+
 } // namespace Slic3r
 
 #endif /* slic3r_KDTreeIndirect_hpp_ */
diff --git a/src/libslic3r/Line.hpp b/src/libslic3r/Line.hpp
index 751e59458..8631ee08b 100644
--- a/src/libslic3r/Line.hpp
+++ b/src/libslic3r/Line.hpp
@@ -113,7 +113,6 @@ public:
     Vector vector() const { return this->b - this->a; }
     Vector normal() const { return Vector((this->b(1) - this->a(1)), -(this->b(0) - this->a(0))); }
     bool   intersection(const Line& line, Point* intersection) const;
-    double ccw(const Point& point) const { return point.ccw(*this); }
     // Clip a line with a bounding box. Returns false if the line is completely outside of the bounding box.
 	bool   clip_with_bbox(const BoundingBox &bbox);
     // Extend the line from both sides by an offset.
diff --git a/src/libslic3r/Point.cpp b/src/libslic3r/Point.cpp
index b2427d46d..9f56f96d6 100644
--- a/src/libslic3r/Point.cpp
+++ b/src/libslic3r/Point.cpp
@@ -108,36 +108,6 @@ bool Point::nearest_point(const Points &points, Point* point) const
     return true;
 }
 
-/* Three points are a counter-clockwise turn if ccw > 0, clockwise if
- * ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
- * gives the signed area of the triangle formed by p1, p2 and this point.
- * In other words it is the 2D cross product of p1-p2 and p1-this, i.e.
- * z-component of their 3D cross product.
- * We return double because it must be big enough to hold 2*max(|coordinate|)^2
- */
-double Point::ccw(const Point &p1, const Point &p2) const
-{
-    static_assert(sizeof(coord_t) == 4, "Point::ccw() requires a 32 bit coord_t");
-    return cross2((p2 - p1).cast<int64_t>(), (*this - p1).cast<int64_t>());
-//    return cross2((p2 - p1).cast<double>(), (*this - p1).cast<double>());
-}
-
-double Point::ccw(const Line &line) const
-{
-    return this->ccw(line.a, line.b);
-}
-
-// returns the CCW angle between this-p1 and this-p2
-// i.e. this assumes a CCW rotation from p1 to p2 around this
-double Point::ccw_angle(const Point &p1, const Point &p2) const
-{
-    //FIXME this calculates an atan2 twice! Project one vector into the other!
-    double angle = atan2(p1.x() - (*this).x(), p1.y() - (*this).y())
-                 - atan2(p2.x() - (*this).x(), p2.y() - (*this).y());
-    // we only want to return only positive angles
-    return angle <= 0 ? angle + 2*PI : angle;
-}
-
 Point Point::projection_onto(const MultiPoint &poly) const
 {
     Point running_projection = poly.first_point();
diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp
index 6b430c2fe..84152ee9c 100644
--- a/src/libslic3r/Point.hpp
+++ b/src/libslic3r/Point.hpp
@@ -28,6 +28,9 @@ using Mat = Eigen::Matrix<T, N, M, Eigen::DontAlign, N, M>;
 
 template<int N, class T> using Vec = Mat<N, 1, T>;
 
+template<typename NumberType>
+using DynVec = Eigen::Matrix<NumberType, Eigen::Dynamic, 1>;
+
 // Eigen types, to replace the Slic3r's own types in the future.
 // Vector types with a fixed point coordinate base type.
 using Vec2crd = Eigen::Matrix<coord_t,  2, 1, Eigen::DontAlign>;
@@ -76,24 +79,35 @@ inline const auto &identity3d = identity<3, double>;
 
 inline bool operator<(const Vec2d &lhs, const Vec2d &rhs) { return lhs.x() < rhs.x() || (lhs.x() == rhs.x() && lhs.y() < rhs.y()); }
 
-template<int Options>
-int32_t cross2(const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v2) = delete;
-
-template<typename T, int Options>
-inline T cross2(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v2)
-{
-    return v1.x() * v2.y() - v1.y() * v2.x();
-}
-
+// Cross product of two 2D vectors.
+// None of the vectors may be of int32_t type as the result would overflow.
 template<typename Derived, typename Derived2>
 inline typename Derived::Scalar cross2(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2)
 {
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
+    static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
+    static_assert(! std::is_same<typename Derived::Scalar, int32_t>::value, "cross2(): Scalar type must not be int32_t, otherwise the cross product would overflow.");
     static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value, "cross2(): Scalar types of 1st and 2nd operand must be equal.");
     return v1.x() * v2.y() - v1.y() * v2.x();
 }
 
-template<typename T, int Options>
-inline Eigen::Matrix<T, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v) { return Eigen::Matrix<T, 2, 1, Eigen::DontAlign>(- v.y(), v.x()); }
+// 2D vector perpendicular to the argument.
+template<typename Derived>
+inline Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Derived> &v)
+{ 
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "perp(): parameter is not a 2D vector");
+    return { - v.y(), v.x() };
+}
+
+// Angle from v1 to v2, returning double atan2(y, x) normalized to <-PI, PI>.
+template<typename Derived, typename Derived2>
+inline double angle(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2) {
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "angle(): first parameter is not a 2D vector");
+    static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "angle(): second parameter is not a 2D vector");
+    auto v1d = v1.template cast<double>();
+    auto v2d = v2.template cast<double>();
+    return atan2(cross2(v1d, v2d), v1d.dot(v2d));
+}
 
 template<class T, int N, int Options>
 Eigen::Matrix<T, 2, 1, Eigen::DontAlign> to_2d(const Eigen::MatrixBase<Eigen::Matrix<T, N, 1, Options>> &ptN) { return { ptN.x(), ptN.y() }; }
@@ -165,9 +179,6 @@ public:
     int    nearest_point_index(const PointConstPtrs &points) const;
     int    nearest_point_index(const PointPtrs &points) const;
     bool   nearest_point(const Points &points, Point* point) const;
-    double ccw(const Point &p1, const Point &p2) const;
-    double ccw(const Line &line) const;
-    double ccw_angle(const Point &p1, const Point &p2) const;
     Point  projection_onto(const MultiPoint &poly) const;
     Point  projection_onto(const Line &line) const;
 };
diff --git a/src/libslic3r/Polygon.cpp b/src/libslic3r/Polygon.cpp
index c7bbe9706..09f1c393d 100644
--- a/src/libslic3r/Polygon.cpp
+++ b/src/libslic3r/Polygon.cpp
@@ -171,50 +171,56 @@ Point Polygon::centroid() const
     return Point(Vec2d(c / (3. * area_sum)));
 }
 
-// find all concave vertices (i.e. having an internal angle greater than the supplied angle)
-// (external = right side, thus we consider ccw orientation)
-Points Polygon::concave_points(double angle) const
-{
-    Points points;
-    angle = 2. * PI - angle + EPSILON;
-    
-    // check whether first point forms a concave angle
-    if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) <= angle)
-        points.push_back(this->points.front());
-    
-    // check whether points 1..(n-1) form concave angles
-    for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++ p)
-        if (p->ccw_angle(*(p-1), *(p+1)) <= angle)
-        	points.push_back(*p);
-    
-    // check whether last point forms a concave angle
-    if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) <= angle)
-        points.push_back(this->points.back());
+// Filter points from poly to the output with the help of FilterFn.
+// filter function receives two vectors:
+// v1: this_point - previous_point
+// v2: next_point - this_point
+// and returns true if the point is to be copied to the output.
+template<typename FilterFn>
+Points filter_points_by_vectors(const Points &poly, FilterFn filter)
+{
+    // Last point is the first point visited.
+    Point p1 = poly.back();
+    // Previous vector to p1.
+    Vec2d v1 = (p1 - *(poly.end() - 2)).cast<double>();
+
+    Points out;
+    for (Point p2 : poly) {
+        // p2 is next point to the currently visited point p1.
+        Vec2d v2 = (p2 - p1).cast<double>();
+        if (filter(v1, v2))
+            out.emplace_back(p2);
+        v1 = v2;
+        p1 = p2;
+    }
     
-    return points;
+    return out;
 }
 
-// find all convex vertices (i.e. having an internal angle smaller than the supplied angle)
-// (external = right side, thus we consider ccw orientation)
-Points Polygon::convex_points(double angle) const
+template<typename ConvexConcaveFilterFn>
+Points filter_convex_concave_points_by_angle_threshold(const Points &poly, double angle_threshold, ConvexConcaveFilterFn convex_concave_filter)
 {
-    Points points;
-    angle = 2*PI - angle - EPSILON;
-    
-    // check whether first point forms a convex angle
-    if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) >= angle)
-        points.push_back(this->points.front());
-    
-    // check whether points 1..(n-1) form convex angles
-    for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++p) {
-        if (p->ccw_angle(*(p-1), *(p+1)) >= angle) points.push_back(*p);
+    assert(angle_threshold >= 0.);
+    if (angle_threshold < EPSILON) {
+        double cos_angle  = cos(angle_threshold);
+        return filter_points_by_vectors(poly, [convex_concave_filter, cos_angle](const Vec2d &v1, const Vec2d &v2){
+            return convex_concave_filter(v1, v2) && v1.normalized().dot(v2.normalized()) < cos_angle;
+        });
+    } else {
+        return filter_points_by_vectors(poly, [convex_concave_filter](const Vec2d &v1, const Vec2d &v2){
+            return convex_concave_filter(v1, v2);
+        });
     }
-    
-    // check whether last point forms a convex angle
-    if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) >= angle)
-        points.push_back(this->points.back());
-    
-    return points;
+}
+
+Points Polygon::convex_points(double angle_threshold) const
+{
+    return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) > 0.; });
+}
+
+Points Polygon::concave_points(double angle_threshold) const
+{
+    return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) < 0.; });
 }
 
 // Projection of a point onto the polygon.
diff --git a/src/libslic3r/Polygon.hpp b/src/libslic3r/Polygon.hpp
index 089820565..77f3020be 100644
--- a/src/libslic3r/Polygon.hpp
+++ b/src/libslic3r/Polygon.hpp
@@ -65,8 +65,11 @@ public:
     void densify(float min_length, std::vector<float>* lengths = nullptr);
     void triangulate_convex(Polygons* polygons) const;
     Point centroid() const;
-    Points concave_points(double angle = PI) const;
-    Points convex_points(double angle = PI) const;
+    // Considering CCW orientation of this polygon, find all convex resp. concave points
+    // with the angle at the vertex larger than a threshold.
+    // Zero angle_threshold means to accept all convex resp. concave points.
+    Points convex_points(double angle_threshold = 0.) const;
+    Points concave_points(double angle_threshold = 0.) const;
     // Projection of a point onto the polygon.
     Point point_projection(const Point &point) const;
     std::vector<float> parameter_by_length() const;
diff --git a/src/libslic3r/SLA/bicubic.h b/src/libslic3r/SLA/bicubic.h
deleted file mode 100644
index 870d00dbd..000000000
--- a/src/libslic3r/SLA/bicubic.h
+++ /dev/null
@@ -1,186 +0,0 @@
-#ifndef BICUBIC_HPP
-#define BICUBIC_HPP
-
-#include <algorithm>
-#include <vector>
-#include <cmath>
-
-#include <Eigen/Dense>
-
-namespace Slic3r {
-
-namespace BicubicInternal {
-	// Linear kernel, to be able to test cubic methods with hat kernels.
-	template<typename T>
-	struct LinearKernel
-	{
-		typedef T	FloatType;
-
-		static T a00() { return T(0.); }
-		static T a01() { return T(0.); }
-		static T a02() { return T(0.); }
-		static T a03() { return T(0.); }
-		static T a10() { return T(1.); }
-		static T a11() { return T(-1.); }
-		static T a12() { return T(0.); }
-		static T a13() { return T(0.); }
-		static T a20() { return T(0.); }
-		static T a21() { return T(1.); }
-		static T a22() { return T(0.); }
-		static T a23() { return T(0.); }
-		static T a30() { return T(0.); }
-		static T a31() { return T(0.); }
-		static T a32() { return T(0.); }
-		static T a33() { return T(0.); }
-	};
-
-	// Interpolation kernel aka Catmul-Rom aka Keyes kernel.
-	template<typename T>
-	struct CubicCatmulRomKernel
-	{
-		typedef T	FloatType;
-
-		static T a00() { return     0;     }
-		static T a01() { return (T)-0.5;   }
-		static T a02() { return (T) 1.;    }
-		static T a03() { return (T)-0.5;   }
-		static T a10() { return (T) 1.;    }
-		static T a11() { return     0;     }
-		static T a12() { return (T)-5./2.; }
-		static T a13() { return (T) 3./2.; }
-		static T a20() { return     0;     }
-		static T a21() { return (T) 0.5;   }
-		static T a22() { return (T) 2.;    }
-		static T a23() { return (T)-3./2.; }
-		static T a30() { return     0;     }
-		static T a31() { return     0;     }
-		static T a32() { return (T)-0.5;   }
-		static T a33() { return (T) 0.5;   }
-	};
-
-	// B-spline kernel
-	template<typename T>
-	struct CubicBSplineKernel
-	{
-		typedef T	FloatType;
-
-		static T a00() { return (T)  1./6.; }
-		static T a01() { return (T) -3./6.; }
-		static T a02() { return (T)  3./6.; }
-		static T a03() { return (T) -1./6.; }
-		static T a10() { return (T)  4./6.; }
-		static T a11() { return      0;     }
-		static T a12() { return (T) -6./6.; }
-		static T a13() { return (T)  3./6.; }
-		static T a20() { return (T)  1./6.; }
-		static T a21() { return (T)  3./6.; }
-		static T a22() { return (T)  3./6.; }
-		static T a23() { return (T)- 3./6.; }
-		static T a30() { return      0;     }
-		static T a31() { return      0;     }
-		static T a32() { return      0;     }
-		static T a33() { return (T)  1./6.; }
-	};
-
-	template<class T>
-	inline T clamp(T a, T lower, T upper)
-	{
-		return (a < lower) ? lower : 
-	   		   (a > upper) ? upper : a;
-	}
-}
-
-template<typename KERNEL>
-struct CubicKernel
-{
-	typedef typename KERNEL					KernelInternal;
-	typedef typename KERNEL::FloatType		FloatType;
-
-	static FloatType kernel(FloatType x)
-	{
-		x = fabs(x);
-		if (x >= (FloatType)2.)
-			return 0.0f;
-		if (x <= (FloatType)1.) {
-			FloatType x2 = x * x;
-			FloatType x3 = x2 * x;
-			return KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3;
-		}
-		assert(x > (FloatType)1. && x < (FloatType)2.);
-		x -= (FloatType)1.;
-		FloatType x2 = x * x;
-		FloatType x3 = x2 * x;
-		return KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3;
-	}
-
-	static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x)
-	{
-		const FloatType  x2 = x*x;
-		const FloatType  x3 = x*x*x;
-		return f0*(KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3) +
-			   f1*(KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3) +
-			   f2*(KERNEL::a20() + KERNEL::a21() * x + KERNEL::a22() * x2 + KERNEL::a23() * x3) + 
-			   f3*(KERNEL::a30() + KERNEL::a31() * x + KERNEL::a32() * x2 + KERNEL::a33() * x3);
-	}
-};
-
-// Linear splines
-typedef CubicKernel<BicubicInternal::LinearKernel<float>>					LinearKernelf;
-typedef CubicKernel<BicubicInternal::LinearKernel<double>>					LinearKerneld;
-// Catmul-Rom splines
-typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>>			CubicCatmulRomKernelf;
-typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>>			CubicCatmulRomKerneld;
-typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>>			CubicInterpolationKernelf;
-typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>>			CubicInterpolationKerneld;
-// Cubic B-splines
-typedef CubicKernel<BicubicInternal::CubicBSplineKernel<float>>				CubicBSplineKernelf;
-typedef CubicKernel<BicubicInternal::CubicBSplineKernel<double>>			CubicBSplineKerneld;
-
-template<typename KERNEL, typename Derived>
-static float cubic_interpolate(const Eigen::ArrayBase<Derived> &F, const typename KERNEL::FloatType pt, const typename KERNEL::FloatType dx)
-{
-	typedef typename KERNEL::FloatType T;
-	const int w  = int(F.size());
-	const int ix = (int)floor(pt);
-	const T   s  = pt - (T)ix;
-
-	if (ix > 1 && ix + 2 < w) {
-		// Inside the fully interpolated region.
-		return KERNEL::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
-	}
-	// Transition region. Extend with a constant function.
-	auto f = [&F, w](x) { return F[BicubicInternal::clamp(x, 0, w - 1)]; }
-	return KERNEL::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s);
-}
-
-template<typename KERNEL, typename Derived>
-static float bicubic_interpolate(const Eigen::MatrixBase<Derived> &F, const Eigen::Matrix<typename KERNEL::FloatType, 2, 1, Eigen::DontAlign> &pt, const typename KERNEL::FloatType dx)
-{
-	typedef typename KERNEL::FloatType T;
-	const int w  = F.cols();
-	const int h  = F.rows();
-	const int ix = (int)floor(pt[0]);
-	const int iy = (int)floor(pt[1]);
-	const T   s  = pt[0] - (T)ix;
-	const T   t  = pt[1] - (T)iy;
-
-	if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) {
-		// Inside the fully interpolated region.
-		return KERNEL::interpolate(
-			KERNEL::interpolate(F(ix-1,iy-1),F(ix  ,iy-1),F(ix+1,iy-1),F(ix+2,iy-1),s),
-			KERNEL::interpolate(F(ix-1,iy  ),F(ix  ,iy  ),F(ix+1,iy  ),F(ix+2,iy  ),s),
-			KERNEL::interpolate(F(ix-1,iy+1),F(ix  ,iy+1),F(ix+1,iy+1),F(ix+2,iy+1),s),
-			KERNEL::interpolate(F(ix-1,iy+2),F(ix  ,iy+2),F(ix+1,iy+2),F(ix+2,iy+2),s),t);
-	}
-	// Transition region. Extend with a constant function.
-	auto f = [&f, w, h](int x, int y) { return F(BicubicInternal::clamp(x,0,w-1),BicubicInternal::clamp(y,0,h-1)); }
-	return KERNEL::interpolate(
-		KERNEL::interpolate(f(ix-1,iy-1),f(ix  ,iy-1),f(ix+1,iy-1),f(ix+2,iy-1),s),
-		KERNEL::interpolate(f(ix-1,iy  ),f(ix  ,iy  ),f(ix+1,iy  ),f(ix+2,iy  ),s),
-		KERNEL::interpolate(f(ix-1,iy+1),f(ix  ,iy+1),f(ix+1,iy+1),f(ix+2,iy+1),s),
-		KERNEL::interpolate(f(ix-1,iy+2),f(ix  ,iy+2),f(ix+1,iy+2),f(ix+2,iy+2),s),t);
-}
-
-} // namespace Slic3r
-
-#endif /* BICUBIC_HPP */
diff --git a/src/libslic3r/Subdivide.cpp b/src/libslic3r/Subdivide.cpp
new file mode 100644
index 000000000..31988a851
--- /dev/null
+++ b/src/libslic3r/Subdivide.cpp
@@ -0,0 +1,218 @@
+#include "Subdivide.hpp"
+#include "Point.hpp"
+
+namespace Slic3r{
+
+indexed_triangle_set its_subdivide(
+    const indexed_triangle_set &its, float max_length)
+{
+    // same order as key order in Edge Divides
+    struct VerticesSequence
+    {
+        size_t   start_index;
+        bool     positive_order;
+        VerticesSequence(size_t start_index, bool positive_order = true)
+            : start_index(start_index), positive_order(positive_order){}
+    };
+    //                         vertex index small, big      vertex index from key.first to key.second
+    using EdgeDivides = std::map<std::pair<size_t, size_t>, VerticesSequence>;
+    struct Edges
+    {
+        Vec3f data[3];
+        Vec3f lengths;
+        Edges(const Vec3crd &indices, const std::vector<Vec3f> &vertices)
+            : lengths(-1.f,-1.f,-1.f)
+        {
+            const Vec3f &v0 = vertices[indices[0]];
+            const Vec3f &v1 = vertices[indices[1]];
+            const Vec3f &v2 = vertices[indices[2]];
+            data[0] = v0 - v1;
+            data[1] = v1 - v2;
+            data[2] = v2 - v0;
+        }
+        float abs_sum(const Vec3f &v)
+        {
+            return abs(v[0]) + abs(v[1]) + abs(v[2]);
+        }
+        bool is_dividable(const float& max_length) {
+            Vec3f sum(abs_sum(data[0]), abs_sum(data[1]), abs_sum(data[2]));
+            Vec3i biggest_index = (sum[0] > sum[1]) ?
+                                      ((sum[0] > sum[2]) ?
+                                           ((sum[2] > sum[1]) ?
+                                                Vec3i(0, 2, 1) :
+                                                Vec3i(0, 1, 2)) :
+                                           Vec3i(2, 0, 1)) :
+                                      ((sum[1] > sum[2]) ?
+                                           ((sum[2] > sum[0]) ?
+                                                Vec3i(1, 2, 0) :
+                                                Vec3i(1, 0, 2)) :
+                                           Vec3i(2, 1, 0));
+            for (int i = 0; i < 3; i++) {
+                int index = biggest_index[i];
+                if (sum[index] <= max_length) return false;
+                lengths[index] = data[index].norm();
+                if (lengths[index] <= max_length) continue;
+
+                // calculate rest of lengths
+                for (int j = i + 1; j < 3; j++) {
+                    index     = biggest_index[j];
+                    lengths[index] = data[index].norm();
+                }
+                return true;
+            }
+            return false;
+        }
+    };
+    struct TriangleLengths
+    {
+        Vec3crd indices;
+        Vec3f l; // lengths
+        TriangleLengths(const Vec3crd &indices, const Vec3f &lengths)
+            : indices(indices), l(lengths)
+        {}
+
+        int get_divide_index(float max_length) {
+            if (l[0] > l[1] && l[0] > l[2]) {
+                if (l[0] > max_length) return 0;
+            } else if (l[1] > l[2]) {
+                if (l[1] > max_length) return 1;
+            } else {
+                if (l[2] > max_length) return 2;
+            }
+            return -1;
+        }
+
+        // divide triangle add new vertex to vertices
+        std::pair<TriangleLengths, TriangleLengths> divide(
+            int divide_index, float max_length,
+            std::vector<Vec3f> &vertices,
+            EdgeDivides &edge_divides)
+        {
+            // index to lengths and indices
+            size_t i0 = divide_index;
+            size_t i1 = (divide_index + 1) % 3;
+            size_t vi0   = indices[i0];
+            size_t vi1   = indices[i1];
+            std::pair<size_t, size_t> key(vi0, vi1);
+            bool key_swap = false;
+            if (key.first > key.second) {
+                std::swap(key.first, key.second);
+                key_swap = true;
+            }
+
+            float length = l[divide_index];
+            size_t count_edge_vertices  = static_cast<size_t>(floor(length / max_length));
+            float count_edge_segments = static_cast<float>(count_edge_vertices + 1);
+
+            auto it = edge_divides.find(key);
+            if (it == edge_divides.end()) {
+                // Create new vertices
+                VerticesSequence new_vs(vertices.size());
+                Vec3f vf = vertices[key.first]; // copy
+                const Vec3f &vs = vertices[key.second];
+                Vec3f dir = vs - vf;
+                for (size_t i = 1; i <= count_edge_vertices; ++i) {
+                    float ratio = i / count_edge_segments;
+                    vertices.push_back(vf + dir * ratio);
+                }
+                bool     success;
+                std::tie(it,success) = edge_divides.insert({key, new_vs});
+                assert(success);
+            }
+            const VerticesSequence &vs = it->second;
+
+            int index_offset = count_edge_vertices/2;
+            size_t i2 = (divide_index + 2) % 3;
+            if (count_edge_vertices % 2 == 0 && key_swap == (l[i1] < l[i2])) {
+                --index_offset;
+            }
+            int sign = (vs.positive_order) ? 1 : -1;
+            size_t new_index = vs.start_index + sign*index_offset;
+
+            size_t vi2   = indices[i2];
+            const Vec3f &v2 = vertices[vi2];
+            Vec3f        new_edge = v2 - vertices[new_index];
+            float        new_len  = new_edge.norm();
+
+            float ratio = (1 + index_offset) / count_edge_segments;
+            float len1 = l[i0] * ratio;
+            float len2 = l[i0] - len1;
+            if (key_swap) std::swap(len1, len2);
+
+            Vec3crd indices1(vi0, new_index, vi2);
+            Vec3f lengths1(len1, new_len, l[i2]);
+
+            Vec3crd indices2(new_index, vi1, vi2);
+            Vec3f lengths2(len2, l[i1], new_len);
+
+            // append key for divided edge when neccesary
+            if (index_offset > 0) {
+                std::pair<size_t, size_t> new_key(key.first, new_index);
+                bool new_key_swap = false;
+                if (new_key.first > new_key.second) {
+                    std::swap(new_key.first, new_key.second);
+                    new_key_swap = true;
+                }
+                if (edge_divides.find(new_key) == edge_divides.end()) {
+                    // insert new
+                    edge_divides.insert({new_key, (new_key_swap) ?
+                        VerticesSequence(new_index - sign, !vs.positive_order)
+                        : vs});
+                }
+            }
+
+            if (index_offset < int(count_edge_vertices)-1) {
+                std::pair<size_t, size_t> new_key(new_index, key.second);
+                bool new_key_swap = false;
+                if (new_key.first > new_key.second) {
+                    std::swap(new_key.first, new_key.second);
+                    new_key_swap = true;
+                }
+                // bad order
+                if (edge_divides.find(new_key) == edge_divides.end()) {
+                    edge_divides.insert({new_key, (new_key_swap) ?
+                        VerticesSequence(vs.start_index + sign*(count_edge_vertices-1), !vs.positive_order)
+                        : VerticesSequence(new_index + sign, vs.positive_order)});
+                }
+            }
+
+            return {TriangleLengths(indices1, lengths1),
+                    TriangleLengths(indices2, lengths2)};
+        }
+    };
+    indexed_triangle_set result;
+    result.indices.reserve(its.indices.size());
+    const std::vector<Vec3f> &vertices = its.vertices;
+    result.vertices = vertices; // copy
+    std::queue<TriangleLengths> tls;
+
+    EdgeDivides edge_divides;
+    for (const Vec3crd &indices : its.indices) {
+        Edges edges(indices, vertices);
+        // speed up only sum not sqrt is apply
+        if (!edges.is_dividable(max_length)) {
+             // small triangle
+            result.indices.push_back(indices);
+            continue;
+        }
+        TriangleLengths tl(indices, edges.lengths);
+        do {
+            int divide_index = tl.get_divide_index(max_length);
+            if (divide_index < 0) {
+                // no dividing
+                result.indices.push_back(tl.indices);
+                if (tls.empty()) break;
+                tl = tls.front(); // copy
+                tls.pop();
+            } else {
+                auto [tl1, tl2] = tl.divide(divide_index, max_length,
+                                            result.vertices, edge_divides);
+                tl = tl1;
+                tls.push(tl2);
+            }
+        } while (true);
+    }
+    return result;
+}
+
+}
diff --git a/src/libslic3r/Subdivide.hpp b/src/libslic3r/Subdivide.hpp
new file mode 100644
index 000000000..f97e4b3c5
--- /dev/null
+++ b/src/libslic3r/Subdivide.hpp
@@ -0,0 +1,12 @@
+#ifndef libslic3r_Subdivide_hpp_
+#define libslic3r_Subdivide_hpp_
+
+#include "TriangleMesh.hpp"
+
+namespace Slic3r {
+
+indexed_triangle_set its_subdivide(const indexed_triangle_set &its, float max_length);
+
+}
+
+#endif //libslic3r_Subdivide_hpp_
diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h
index 5e0fb67b3..2131a9233 100644
--- a/src/libslic3r/libslic3r.h
+++ b/src/libslic3r/libslic3r.h
@@ -23,6 +23,13 @@
 #include <cmath>
 #include <type_traits>
 
+#ifdef _WIN32
+// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
+// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740
+// Thus it is recommended to use boost::container::deque instead.
+#include <boost/container/deque.hpp>
+#endif // _WIN32
+
 #include "Technologies.hpp"
 #include "Semver.hpp"
 
@@ -73,6 +80,16 @@ namespace Slic3r {
 
 extern Semver SEMVER;
 
+// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
+template<class T, class Allocator = std::allocator<T>>
+using deque = 
+#ifdef _WIN32
+    // Use boost implementation, which allocates blocks of 512 bytes instead of blocks of 8 bytes.
+    boost::container::deque<T, Allocator>;
+#else // _WIN32
+    std::deque<T, Allocator>;
+#endif // _WIN32
+
 template<typename T, typename Q>
 inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); }
 
diff --git a/src/libslic3r/pchheader.hpp b/src/libslic3r/pchheader.hpp
index e6591f574..aeb79e3d8 100644
--- a/src/libslic3r/pchheader.hpp
+++ b/src/libslic3r/pchheader.hpp
@@ -64,6 +64,12 @@
 #include <boost/config.hpp>
 #include <boost/config/warning_disable.hpp>
 #include <boost/container/small_vector.hpp>
+#ifdef _WIN32
+// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
+// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740
+// Thus it is recommended to use boost::container::deque instead.
+#include <boost/container/deque.hpp>
+#endif // _WIN32
 #include <boost/date_time/local_time/local_time.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/filesystem.hpp>
diff --git a/t/geometry.t b/t/geometry.t
index 0f37c0aa3..bb72b22e1 100644
--- a/t/geometry.t
+++ b/t/geometry.t
@@ -2,7 +2,7 @@ use Test::More;
 use strict;
 use warnings;
 
-plan tests => 42;
+plan tests => 30;
 
 BEGIN {
     use FindBin;
@@ -189,63 +189,6 @@ my $polygons = [
     is +Slic3r::Point->new(10, 0)->distance_to_line($line), 0, 'distance_to';
 }
 
-#==========================================================
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [100,100],
-        [200,100],
-        [200,200],
-        [100,200],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in ccw square';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in ccw square';
-    
-    $square->make_clockwise;
-    is scalar(@{$square->concave_points(PI*4/3)}), 4, 'fuor concave vertices detected in cw square';
-    is scalar(@{$square->convex_points(PI*2/3)}), 0, 'no convex vertices detected in cw square';
-}
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [150,100],
-        [200,100],
-        [200,200],
-        [100,200],
-        [100,100],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
-}
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [200,200],
-        [100,200],
-        [100,100],
-        [150,100],
-        [200,100],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
-}
-
-{
-    my $triangle = Slic3r::Polygon->new(
-        [16000170,26257364], [714223,461012], [31286371,461008],
-    );
-    is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle';
-    is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle';
-}
-
-{
-    my $triangle = Slic3r::Polygon->new(
-        [16000170,26257364], [714223,461012], [20000000,461012], [31286371,461012],
-    );
-    is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle having collinear point';
-    is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle having collinear point';
-}
-
 {
     my $triangle = Slic3r::Polygon->new(
         [16000170,26257364], [714223,461012], [31286371,461008],
diff --git a/t/perimeters.t b/t/perimeters.t
index d3f96f122..adc2a6cec 100644
--- a/t/perimeters.t
+++ b/t/perimeters.t
@@ -223,7 +223,7 @@ use Slic3r::Test;
                         
                         if ($model eq 'cube_with_concave_hole') {
                             # check that loop starts at a concave vertex
-                            my $ccw_angle = $loop->first_point->ccw_angle(@$loop[-2,1]);
+                            my $ccw_angle = $loop->[-2]->ccw($loop->first_point, $loop->[1]);
                             my $convex = ($ccw_angle > PI);  # whether the angle on the *right* side is convex
                             $starts_on_convex_point = 1
                                 if ($convex && $is_contour) || (!$convex && $is_hole);
diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt
index 69470aad1..248245182 100644
--- a/tests/libslic3r/CMakeLists.txt
+++ b/tests/libslic3r/CMakeLists.txt
@@ -8,6 +8,7 @@ add_executable(${_TEST_NAME}_tests
 	test_clipper_utils.cpp
 	test_color.cpp
 	test_config.cpp
+	test_curve_fitting.cpp
 	test_elephant_foot_compensation.cpp
 	test_geometry.cpp
 	test_placeholder_parser.cpp
diff --git a/tests/libslic3r/test_curve_fitting.cpp b/tests/libslic3r/test_curve_fitting.cpp
new file mode 100644
index 000000000..faf7839c7
--- /dev/null
+++ b/tests/libslic3r/test_curve_fitting.cpp
@@ -0,0 +1,118 @@
+#include <catch2/catch.hpp>
+#include <test_utils.hpp>
+
+#include <libslic3r/Geometry/Curves.hpp>
+#include <libslic3r/Utils.hpp>
+#include <libslic3r/SVG.hpp>
+
+TEST_CASE("Curves: cubic b spline fit test", "[Curves]") {
+    using namespace Slic3r;
+    using namespace Slic3r::Geometry;
+
+    auto fx = [&](size_t index) {
+        return float(index) / 200.0f;
+    };
+
+    auto fy = [&](size_t index) {
+        return 1.0f;
+    };
+
+    std::vector<Vec<1, float>> observations { };
+    std::vector<float> observation_points { };
+    std::vector<float> weights { };
+    for (size_t index = 0; index < 200; ++index) {
+        observations.push_back(Vec<1, float> { fy(index) });
+        observation_points.push_back(fx(index));
+        weights.push_back(1);
+    }
+
+    Vec2f fmin { fx(0), fy(0) };
+    Vec2f fmax { fx(200), fy(200) };
+
+    auto bspline = fit_cubic_bspline(observations, observation_points, weights, 1);
+
+    Approx ap(1.0f);
+    ap.epsilon(0.1f);
+
+    for (int p = 0; p < 200; ++p) {
+        float fitted_val = bspline.get_fitted_value(fx(p))(0);
+        float expected = fy(p);
+
+        REQUIRE(fitted_val == ap(expected));
+
+    }
+}
+
+TEST_CASE("Curves: quadratic f cubic b spline fit test", "[Curves]") {
+    using namespace Slic3r;
+    using namespace Slic3r::Geometry;
+
+    auto fx = [&](size_t index) {
+        return float(index) / 100.0f;
+    };
+
+    auto fy = [&](size_t index) {
+        return (fx(index) - 1) * (fx(index) - 1);
+    };
+
+    std::vector<Vec<1, float>> observations { };
+    std::vector<float> observation_points { };
+    std::vector<float> weights { };
+    for (size_t index = 0; index < 200; ++index) {
+        observations.push_back(Vec<1, float> { fy(index) });
+        observation_points.push_back(fx(index));
+        weights.push_back(1);
+    }
+
+    Vec2f fmin { fx(0), fy(0) };
+    Vec2f fmax { fx(200), fy(200) };
+
+    auto bspline = fit_cubic_bspline(observations, observation_points, weights, 10);
+
+    for (int p = 0; p < 200; ++p) {
+        float fitted_val = bspline.get_fitted_value(fx(p))(0);
+        float expected = fy(p);
+
+        auto check = [](float a, float b) {
+            return abs(a - b) < 0.2f;
+        };
+        //Note: checking is problematic, splines will not perfectly align
+        REQUIRE(check(fitted_val, expected));
+
+    }
+}
+
+TEST_CASE("Curves: polynomial fit test", "[Curves]") {
+    using namespace Slic3r;
+    using namespace Slic3r::Geometry;
+
+    auto fx = [&](size_t index) {
+        return float(index) / 100.0f;
+    };
+
+    auto fy = [&](size_t index) {
+        return (fx(index) - 1) * (fx(index) - 1);
+    };
+
+    std::vector<Vec<1, float>> observations { };
+    std::vector<float> observation_points { };
+    std::vector<float> weights { };
+    for (size_t index = 0; index < 200; ++index) {
+        observations.push_back(Vec<1, float> { fy(index) });
+        observation_points.push_back(fx(index));
+        weights.push_back(1);
+    }
+
+    Vec2f fmin { fx(0), fy(0) };
+    Vec2f fmax { fx(200), fy(200) };
+
+    Approx ap(1.0f);
+    ap.epsilon(0.1f);
+
+    auto poly = fit_polynomial(observations, observation_points, weights, 2);
+
+    REQUIRE(poly.coefficients(0, 0) == ap(1));
+    REQUIRE(poly.coefficients(0, 1) == ap(-2));
+    REQUIRE(poly.coefficients(0, 2) == ap(1));
+}
+
diff --git a/tests/libslic3r/test_geometry.cpp b/tests/libslic3r/test_geometry.cpp
index 21f9a9663..34e625e6d 100644
--- a/tests/libslic3r/test_geometry.cpp
+++ b/tests/libslic3r/test_geometry.cpp
@@ -373,7 +373,43 @@ SCENARIO("Line distances", "[Geometry]"){
     }
 }
 
+SCENARIO("Calculating angles", "[Geometry]")
+{
+    GIVEN(("Vectors 30 degrees apart"))
+    {
+        std::vector<std::pair<Point, Point>> pts {
+            { {1000, 0}, { 866, 500 } },
+            { { 866, 500 }, { 500, 866 } },
+            { { 500, 866 }, { 0, 1000 } },
+            { { -500, 866 }, { -866, 500 } }
+        };
+
+        THEN("Angle detected is 30 degrees")
+        {
+            for (auto &p : pts)
+                REQUIRE(is_approx(angle(p.first, p.second), M_PI / 6.));
+        }
+    }
+
+    GIVEN(("Vectors 30 degrees apart"))
+    {
+        std::vector<std::pair<Point, Point>> pts {
+            { { 866, 500 }, {1000, 0} },
+            { { 500, 866 }, { 866, 500 } },
+            { { 0, 1000 }, { 500, 866 } },
+            { { -866, 500 }, { -500, 866 } }
+        };
+
+        THEN("Angle detected is -30 degrees")
+        {
+            for (auto &p : pts)
+                REQUIRE(is_approx(angle(p.first, p.second), - M_PI / 6.));
+        }
+    }
+}
+
 SCENARIO("Polygon convex/concave detection", "[Geometry]"){
+    static constexpr const double angle_threshold = M_PI / 3.;
     GIVEN(("A Square with dimension 100")){
         auto square = Slic3r::Polygon /*new_scale*/(std::vector<Point>({
             Point(100,100),
@@ -381,13 +417,13 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(200,200),
             Point(100,200)}));
         THEN("It has 4 convex points counterclockwise"){
-            REQUIRE(square.concave_points(PI*4/3).size() == 0);
-            REQUIRE(square.convex_points(PI*2/3).size() == 4);
+            REQUIRE(square.concave_points(angle_threshold).size() == 0);
+            REQUIRE(square.convex_points(angle_threshold).size() == 4);
         }
         THEN("It has 4 concave points clockwise"){
             square.make_clockwise();
-            REQUIRE(square.concave_points(PI*4/3).size() == 4);
-            REQUIRE(square.convex_points(PI*2/3).size() == 0);
+            REQUIRE(square.concave_points(angle_threshold).size() == 4);
+            REQUIRE(square.convex_points(angle_threshold).size() == 0);
         }
     }
     GIVEN("A Square with an extra colinearvertex"){
@@ -398,8 +434,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(100,200),
             Point(100,100)}));
         THEN("It has 4 convex points counterclockwise"){
-            REQUIRE(square.concave_points(PI*4/3).size() == 0);
-            REQUIRE(square.convex_points(PI*2/3).size() == 4);
+            REQUIRE(square.concave_points(angle_threshold).size() == 0);
+            REQUIRE(square.convex_points(angle_threshold).size() == 4);
         }
     }
     GIVEN("A Square with an extra collinear vertex in different order"){
@@ -410,8 +446,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(150,100),
             Point(200,100)}));
         THEN("It has 4 convex points counterclockwise"){
-            REQUIRE(square.concave_points(PI*4/3).size() == 0);
-            REQUIRE(square.convex_points(PI*2/3).size() == 4);
+            REQUIRE(square.concave_points(angle_threshold).size() == 0);
+            REQUIRE(square.convex_points(angle_threshold).size() == 4);
         }
     }
 
@@ -422,8 +458,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(31286371,461008)
         }));
         THEN("it has three convex vertices"){
-            REQUIRE(triangle.concave_points(PI*4/3).size() == 0);
-            REQUIRE(triangle.convex_points(PI*2/3).size() == 3);
+            REQUIRE(triangle.concave_points(angle_threshold).size() == 0);
+            REQUIRE(triangle.convex_points(angle_threshold).size() == 3);
         }
     }
 
@@ -435,8 +471,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(31286371,461012)
         }));
         THEN("it has three convex vertices"){
-            REQUIRE(triangle.concave_points(PI*4/3).size() == 0);
-            REQUIRE(triangle.convex_points(PI*2/3).size() == 3);
+            REQUIRE(triangle.concave_points(angle_threshold).size() == 0);
+            REQUIRE(triangle.convex_points(angle_threshold).size() == 3);
         }
     }
     GIVEN("A polygon with concave vertices with angles of specifically 4/3pi"){
@@ -453,8 +489,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
             Point(38092663,692699),Point(52100125,692699)
         }));
         THEN("the correct number of points are detected"){
-            REQUIRE(polygon.concave_points(PI*4/3).size() == 6);
-            REQUIRE(polygon.convex_points(PI*2/3).size() == 10);
+            REQUIRE(polygon.concave_points(angle_threshold).size() == 6);
+            REQUIRE(polygon.convex_points(angle_threshold).size() == 10);
         }
     }
 }
diff --git a/xs/xsp/Line.xsp b/xs/xsp/Line.xsp
index 777dc41fa..36181c3ba 100644
--- a/xs/xsp/Line.xsp
+++ b/xs/xsp/Line.xsp
@@ -41,7 +41,7 @@
     Clone<Point> normal();
     Clone<Point> vector();
     double ccw(Point* point)
-        %code{% RETVAL = THIS->ccw(*point); %};
+        %code{% RETVAL = cross2((THIS->a - *point).cast<double>(), (THIS->b - THIS->a).cast<double>()); %};
 %{
 
 Line*
diff --git a/xs/xsp/Point.xsp b/xs/xsp/Point.xsp
index beefc6249..2d70e0203 100644
--- a/xs/xsp/Point.xsp
+++ b/xs/xsp/Point.xsp
@@ -39,13 +39,7 @@
     double perp_distance_to_line(Line* line)
         %code{% RETVAL = line->perp_distance_to(*THIS); %};
     double ccw(Point* p1, Point* p2)
-        %code{% RETVAL = THIS->ccw(*p1, *p2); %};
-    double ccw_angle(Point* p1, Point* p2)
-        %code{% RETVAL = THIS->ccw_angle(*p1, *p2); %};
-    Point* projection_onto_polygon(Polygon* polygon)
-        %code{% RETVAL = new Point(THIS->projection_onto(*polygon)); %};
-    Point* projection_onto_polyline(Polyline* polyline)
-        %code{% RETVAL = new Point(THIS->projection_onto(*polyline)); %};
+        %code{% RETVAL = cross2((*p1 - *THIS).cast<double>(), (*p2 - *p1).cast<double>()); %};
     Point* projection_onto_line(Line* line)
         %code{% RETVAL = new Point(THIS->projection_onto(*line)); %};
     Point* negative()
diff --git a/xs/xsp/Polygon.xsp b/xs/xsp/Polygon.xsp
index a94425477..6b6d52524 100644
--- a/xs/xsp/Polygon.xsp
+++ b/xs/xsp/Polygon.xsp
@@ -39,8 +39,6 @@
         %code{% THIS->triangulate_convex(&RETVAL); %};
     Clone<Point> centroid();
     Clone<BoundingBox> bounding_box();
-    Points concave_points(double angle);
-    Points convex_points(double angle);
     Clone<Point> point_projection(Point* point)
         %code{% RETVAL = THIS->point_projection(*point); %};
     Clone<Point> intersection(Line* line)