diff --git a/lib/Slic3r/GUI/PresetEditor.pm b/lib/Slic3r/GUI/PresetEditor.pm
index 447be9228b..d7d885bcae 100644
--- a/lib/Slic3r/GUI/PresetEditor.pm
+++ b/lib/Slic3r/GUI/PresetEditor.pm
@@ -439,7 +439,8 @@ sub options {
         adaptive_slicing adaptive_slicing_quality match_horizontal_surfaces
         perimeters spiral_vase
         top_solid_layers min_shell_thickness min_top_bottom_shell_thickness bottom_solid_layers
-        extra_perimeters avoid_crossing_perimeters thin_walls overhangs
+        extra_perimeters avoid_crossing_perimeters 
+        thin_walls thin_walls_min_width thin_walls_overlap overhangs
         seam_position external_perimeters_first
         fill_density fill_pattern top_infill_pattern bottom_infill_pattern fill_gaps
         infill_every_layers infill_only_where_needed
@@ -540,7 +541,13 @@ sub build {
             my $optgroup = $page->new_optgroup('Quality (slower slicing)');
             $optgroup->append_single_option_line('extra_perimeters');
             $optgroup->append_single_option_line('avoid_crossing_perimeters');
-            $optgroup->append_single_option_line('thin_walls');
+            my $line = Slic3r::GUI::OptionsGroup::Line->new(
+                label => 'Detect thin walls',
+            );
+            $line->append_option($optgroup->get_option('thin_walls'));
+            $line->append_option($optgroup->get_option('thin_walls_min_width'));
+            $line->append_option($optgroup->get_option('thin_walls_overlap'));
+            $optgroup->append_line($line);
             $optgroup->append_single_option_line('overhangs');
         }
         {
@@ -888,8 +895,9 @@ sub _update {
     my $have_perimeters = ($config->perimeters > 0) || ($config->min_shell_thickness > 0);
     if (any { /$opt_key/ } qw(all_keys perimeters)) {
         $self->get_field($_)->toggle($have_perimeters)
-            for qw(extra_perimeters thin_walls overhangs seam_position external_perimeters_first
-                    external_perimeter_extrusion_width
+            for qw(extra_perimeters thin_walls thin_walls_min_width thin_walls_overlap
+                    overhangs seam_position 
+                    external_perimeters_first external_perimeter_extrusion_width
                     perimeter_speed small_perimeter_speed external_perimeter_speed);
     }
 
@@ -975,6 +983,10 @@ sub _update {
     $self->get_field($_)->toggle($have_support_material && $have_support_interface)
         for qw(support_material_interface_spacing support_material_interface_extruder
             support_material_interface_speed);
+        
+    # thin walls settigns only when thins walls is activated
+    $self->get_field($_)->toggle($config->thin_walls)
+        for qw(thin_walls_min_width thin_walls_overlap);
     
     $self->get_field('perimeter_extrusion_width')->toggle($have_perimeters || $have_skirt || $have_brim);
     $self->get_field('support_material_extruder')->toggle($have_support_material || $have_skirt);
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d288999c1b..efccfb289a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -204,6 +204,7 @@ add_library(libslic3r STATIC
     ${LIBDIR}/libslic3r/LayerHeightSpline.cpp
     ${LIBDIR}/libslic3r/Line.cpp
     ${LIBDIR}/libslic3r/Log.cpp
+    ${LIBDIR}/libslic3r/MedialAxis.cpp
     ${LIBDIR}/libslic3r/Model.cpp
     ${LIBDIR}/libslic3r/MotionPlanner.cpp
     ${LIBDIR}/libslic3r/MultiPoint.cpp
@@ -316,6 +317,7 @@ set(SLIC3R_TEST_SOURCES
     ${TESTDIR}/libslic3r/test_printobject.cpp
     ${TESTDIR}/libslic3r/test_skirt_brim.cpp
     ${TESTDIR}/libslic3r/test_test_data.cpp
+    ${TESTDIR}/libslic3r/test_thin.cpp
     ${TESTDIR}/libslic3r/test_transformationmatrix.cpp
     ${TESTDIR}/libslic3r/test_trianglemesh.cpp
     ${TESTDIR}/libslic3r/test_extrusion_entity.cpp
diff --git a/src/GUI/Dialogs/PresetEditor.hpp b/src/GUI/Dialogs/PresetEditor.hpp
index 71c3eabda3..7287d7fc51 100644
--- a/src/GUI/Dialogs/PresetEditor.hpp
+++ b/src/GUI/Dialogs/PresetEditor.hpp
@@ -143,7 +143,8 @@ class PrintEditor : public PresetEditor {
             "adaptive_slicing"s, "adaptive_slicing_quality"s, "match_horizontal_surfaces"s,
             "perimeters"s, "spiral_vase"s,
             "top_solid_layers"s, "bottom_solid_layers"s,
-            "extra_perimeters"s, "avoid_crossing_perimeters"s, "thin_walls"s, "overhangs"s,
+            "extra_perimeters"s, "avoid_crossing_perimeters"s, 
+            "thin_walls"s, "thin_walls_min_width"s, "thin_walls_overlap"s, "overhangs"s,
             "seam_position"s, "external_perimeters_first"s,
             "fill_density"s, "fill_pattern"s, "top_infill_pattern"s, "bottom_infill_pattern"s, "fill_gaps"s,
             "infill_every_layers"s, "infill_only_where_needed"s,
diff --git a/src/test/libslic3r/test_fill.cpp b/src/test/libslic3r/test_fill.cpp
index a1ce98b4c7..0fe4751d63 100644
--- a/src/test/libslic3r/test_fill.cpp
+++ b/src/test/libslic3r/test_fill.cpp
@@ -421,6 +421,7 @@ bool test_if_solid_surface_filled(const ExPolygon& expolygon, double flow_spacin
     Flow flow(flow_spacing, 0.4, flow_spacing);
 
     filler->min_spacing = flow.spacing();
+    filler->density = density;
 
     Polylines paths {filler->fill_surface(surface)};
 
diff --git a/src/test/libslic3r/test_flow.cpp b/src/test/libslic3r/test_flow.cpp
index 1db8cd5fdf..8db6f462bf 100644
--- a/src/test/libslic3r/test_flow.cpp
+++ b/src/test/libslic3r/test_flow.cpp
@@ -93,10 +93,10 @@ SCENARIO(" Bridge flow specifics.", "[!mayfail]") {
 SCENARIO("Flow: Flow math for non-bridges", "[!mayfail]") {
     GIVEN("Nozzle Diameter of 0.4, a desired width of 1mm and layer height of 0.5") {
         auto width {ConfigOptionFloatOrPercent(1.0, false)};
-        float spacing {0.4};
-        float nozzle_diameter {0.4};
-        float bridge_flow {1.0};
-        float layer_height {0.5};
+        float spacing {0.4f};
+        float nozzle_diameter {0.4f};
+        float bridge_flow {1.0f};
+        float layer_height {0.5f};
 
         // Spacing for non-bridges is has some overlap
         THEN("External perimeter flow has spacing fixed to 1.1*nozzle_diameter") {
@@ -117,10 +117,10 @@ SCENARIO("Flow: Flow math for non-bridges", "[!mayfail]") {
     }
     /// Check the min/max
     GIVEN("Nozzle Diameter of 0.25") {
-        float spacing {0.4};
-        float nozzle_diameter {0.25};
-        float bridge_flow {0.0};
-        float layer_height {0.5};
+        float spacing {0.4f};
+        float nozzle_diameter {0.25f};
+        float bridge_flow {0.0f};
+        float layer_height {0.5f};
         WHEN("layer height is set to 0.2") {
             layer_height = 0.15f;
             THEN("Max width is set.") {
diff --git a/src/test/libslic3r/test_thin.cpp b/src/test/libslic3r/test_thin.cpp
new file mode 100644
index 0000000000..ab32db5ff4
--- /dev/null
+++ b/src/test/libslic3r/test_thin.cpp
@@ -0,0 +1,470 @@
+
+//#define CATCH_CONFIG_DISABLE
+
+#include <catch.hpp>
+#include "../test_data.hpp"
+#include "../../libslic3r/ClipperUtils.hpp"
+#include "../../libslic3r/MedialAxis.hpp"
+#include "../../libslic3r/SVG.hpp"
+
+using namespace Slic3r;
+using namespace Slic3r::Geometry;
+
+
+SCENARIO("thin walls: ")
+{
+
+
+    GIVEN("Square")
+    {
+        Points test_set;
+        test_set.reserve(4);
+        Points square {Point::new_scale(100, 100),
+			Point::new_scale(200, 100),
+			Point::new_scale(200, 200),
+			Point::new_scale(100, 200)};
+        Slic3r::Polygon hole_in_square{ Points{
+            Point::new_scale(140, 140),
+            Point::new_scale(140, 160),
+            Point::new_scale(160, 160),
+            Point::new_scale(160, 140) } };
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ square };
+        expolygon.holes = Slic3r::Polygons{ hole_in_square };
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(40), scale_(0.5), &res);
+
+            THEN("medial axis of a square shape is a single path"){
+                REQUIRE(res.size() == 1);
+            }
+            THEN("polyline forms a closed loop"){
+                REQUIRE(res[0].first_point().coincides_with(res[0].last_point()) == true);
+            }
+            THEN("medial axis loop has reasonable length"){
+                REQUIRE(res[0].length() > expolygon.holes[0].length());
+                REQUIRE(res[0].length() < expolygon.contour.length());
+            }
+        }
+    }
+
+    GIVEN("narrow rectangle"){
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(100, 100),
+            Point::new_scale(120, 100),
+            Point::new_scale(120, 200),
+            Point::new_scale(100, 200) } };
+        Polylines res;
+        expolygon.medial_axis(scale_(20), scale_(0.5), &res);
+
+        ExPolygon expolygon2;
+        expolygon2.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(100, 100),
+            Point::new_scale(120, 100),
+            Point::new_scale(120, 200),
+            Point::new_scale(105, 200), // extra point in the short side
+            Point::new_scale(100, 200) } };
+        Polylines res2;
+        expolygon.medial_axis(scale_(20), scale_(0.5), &res2);
+        WHEN("creating the medial axis"){
+
+            THEN("medial axis of a narrow rectangle is a single line"){
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has reasonable length") {
+                    REQUIRE(res[0].length() >= scale_(200 - 100 - (120 - 100)) - SCALED_EPSILON);
+                }
+            }
+            THEN("medial axis of a narrow rectangle with an extra vertex is still a single line"){
+                REQUIRE(res2.size() == 1);
+                THEN("medial axis of a narrow rectangle with an extra vertex has reasonable length") {
+                    REQUIRE(res2[0].length() >= scale_(200 - 100 - (120 - 100)) - SCALED_EPSILON);
+                }
+                THEN("extra vertices don't influence medial axis") {
+                    REQUIRE(res2[0].length() - res[0].length() < SCALED_EPSILON);
+                }
+            }
+        }
+    }
+
+    //TODO: compare with mainline slic3r
+    GIVEN("semicicumference") {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point{ 1185881, 829367 }, Point{ 1421988, 1578184 }, Point{ 1722442, 2303558 }, Point{ 2084981, 2999998 }, Point{ 2506843, 3662186 }, Point{ 2984809, 4285086 }, Point{ 3515250, 4863959 }, Point{ 4094122, 5394400 }, Point{ 4717018, 5872368 }, 
+            Point{ 5379210, 6294226 }, Point{ 6075653, 6656769 }, Point{ 6801033, 6957229 }, Point{ 7549842, 7193328 }, Point{ 8316383, 7363266 }, Point{ 9094809, 7465751 }, Point{ 9879211, 7500000 }, Point{ 10663611, 7465750 }, Point{ 11442038, 7363265 }, 
+            Point{ 12208580, 7193327 }, Point{ 12957389, 6957228 }, Point{ 13682769, 6656768 }, Point{ 14379209, 6294227 }, Point{ 15041405, 5872366 }, Point{ 15664297, 5394401 }, Point{ 16243171, 4863960 }, Point{ 16758641, 4301424 }, Point{ 17251579, 3662185 }, 
+            Point{ 17673439, 3000000 }, Point{ 18035980, 2303556 }, Point{ 18336441, 1578177 }, Point{ 18572539, 829368 }, Point{ 18750748, 0 }, Point{ 19758422, 0 }, Point{ 19727293, 236479 }, Point{ 19538467, 1088188 }, Point{ 19276136, 1920196 }, 
+            Point{ 18942292, 2726179 }, Point{ 18539460, 3499999 }, Point{ 18070731, 4235755 }, Point{ 17539650, 4927877 }, Point{ 16950279, 5571067 }, Point{ 16307090, 6160437 }, Point{ 15614974, 6691519 }, Point{ 14879209, 7160248 }, Point{ 14105392, 7563079 }, 
+            Point{ 13299407, 7896927 }, Point{ 12467399, 8159255 }, Point{ 11615691, 8348082 }, Point{ 10750769, 8461952 }, Point{ 9879211, 8500000 }, Point{ 9007652, 8461952 }, Point{ 8142729, 8348082 }, Point{ 7291022, 8159255 }, Point{ 6459015, 7896927 },
+            Point{ 5653029, 7563079 }, Point{ 4879210, 7160247 }, Point{ 4143447, 6691519 }, Point{ 3451331, 6160437 }, Point{ 2808141, 5571066 }, Point{ 2218773, 4927878 }, Point{ 1687689, 4235755 }, Point{ 1218962, 3499999 }, Point{ 827499, 2748020 }, 
+            Point{ 482284, 1920196 }, Point{ 219954, 1088186 }, Point{ 31126, 236479 }, Point{ 0, 0 }, Point{ 1005754, 0 }
+        } };
+
+        WHEN("creating the medial axis") {
+            Polylines res;
+            expolygon.medial_axis(scale_(1.324888), scale_(0.25), &res);
+
+            THEN("medial axis of a semicircumference is a single line") {
+                REQUIRE(res.size() == 1);
+            }
+            THEN("all medial axis segments of a semicircumference have the same orientation (but the 2 end points)") {
+                Lines lines = res[0].lines();
+                double min_angle = 1, max_angle = -1;
+                //std::cout << "first angle=" << lines[0].ccw(lines[1].b) << "\n";
+                for (int idx = 1; idx < lines.size() - 1; idx++) {
+                    double angle = lines[idx - 1].ccw(lines[idx].b);
+                    if (std::abs(angle) - EPSILON < 0) angle = 0;
+                    //if (angle < 0) std::cout << unscale_(lines[idx - 1].a.x()) << ":" << unscale_(lines[idx - 1].a.y()) << " -> " << unscale_(lines[idx - 1].b.x()) << ":" << unscale_(lines[idx - 1].b.y()) << " -> " << unscale_(lines[idx].b.x()) << ":" << unscale_(lines[idx].b.y()) << "\n";
+                    std::cout << "angle=" << 180 * lines[idx].a.ccw_angle(lines[idx - 1].a, lines[idx].b) / PI << "\n";
+                    min_angle = std::min(min_angle, angle);
+                    max_angle = std::max(max_angle, angle);
+                }
+                //std::cout << "last angle=" << lines[lines.size() - 2].ccw(lines[lines.size() - 1].b) << "\n";
+                // check whether turns are all CCW or all CW
+                bool allccw = (min_angle <= 0 && max_angle <= 0);
+                bool allcw = (min_angle >= 0 && max_angle >= 0);
+                bool allsame_orientation = allccw || allcw;
+                REQUIRE(allsame_orientation);
+            }
+        }
+    }
+    
+
+    GIVEN("round with large and very small distance between points"){
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(15.181601,-2.389639), Point::new_scale(15.112616,-1.320034), Point::new_scale(14.024491,-0.644338), Point::new_scale(13.978982,-0.624495), Point::new_scale(9.993299,0.855584), Point::new_scale(9.941970,0.871195), Point::new_scale(5.796743,1.872643),
+            Point::new_scale(5.743826,1.882168), Point::new_scale(1.509170,2.386464), Point::new_scale(1.455460,2.389639), Point::new_scale(-2.809359,2.389639), Point::new_scale(-2.862805,2.386464), Point::new_scale(-7.097726,1.882168), Point::new_scale(-7.150378,1.872643), Point::new_scale(-11.286344,0.873576),
+            Point::new_scale(-11.335028,0.858759), Point::new_scale(-14.348632,-0.237938), Point::new_scale(-14.360538,-0.242436), Point::new_scale(-15.181601,-0.737570), Point::new_scale(-15.171309,-2.388509)
+        } };
+        expolygon.holes.push_back(Slic3r::Polygon{ Points{
+            Point::new_scale( -11.023311,-1.034226 ), Point::new_scale( -6.920984,-0.042941 ), Point::new_scale( -2.768613,0.463207 ), Point::new_scale( 1.414714,0.463207 ), Point::new_scale( 5.567085,-0.042941 ), Point::new_scale( 9.627910,-1.047563 )
+            } });
+
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(2.5), scale_(0.5), &res);
+
+           THEN("medial axis of it is two line"){
+                REQUIRE(res.size() == 2);
+            }
+        }
+    }
+
+    GIVEN("french cross")
+    {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(4.3, 4), Point::new_scale(4.3, 0), Point::new_scale(4, 0), Point::new_scale(4, 4), Point::new_scale(0, 4), Point::new_scale(0, 4.5), Point::new_scale(4, 4.5), Point::new_scale(4, 10), Point::new_scale(4.3, 10), Point::new_scale(4.3, 4.5),
+            Point::new_scale(6, 4.5), Point::new_scale(6, 10), Point::new_scale(6.2, 10), Point::new_scale(6.2, 4.5), Point::new_scale(10, 4.5), Point::new_scale(10, 4), Point::new_scale(6.2, 4), Point::new_scale(6.2, 0), Point::new_scale(6, 0), Point::new_scale(6, 4),
+        } };
+        expolygon.contour.make_counter_clockwise();
+
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(0.55), scale_(0.25), &res);
+
+            THEN("medial axis of a (bit too narrow) french cross is two lines"){
+                REQUIRE(res.size() == 2);
+            }
+            THEN("medial axis has reasonable length"){
+                REQUIRE(res[0].length() >= scale_(9.9) - SCALED_EPSILON);
+                REQUIRE(res[1].length() >= scale_(9.9) - SCALED_EPSILON);
+            }
+
+            THEN("medial axis of a (bit too narrow) french cross is two lines has only strait lines (first line)"){
+                Lines lines = res[0].lines();
+                double min_angle = 1, max_angle = -1;
+                for (int idx = 1; idx < lines.size(); idx++){
+                    double angle = lines[idx - 1].ccw(lines[idx].b);
+                    min_angle = std::min(min_angle, angle);
+                    max_angle = std::max(max_angle, angle);
+                }
+                REQUIRE(min_angle == max_angle);
+                REQUIRE(min_angle == 0);
+            }
+            THEN("medial axis of a (bit too narrow) french cross is two lines has only strait lines (second line)"){
+                Lines lines = res[1].lines();
+                double min_angle = 1, max_angle = -1;
+                for (int idx = 1; idx < lines.size(); idx++){
+                    double angle = lines[idx - 1].ccw(lines[idx].b);
+                    min_angle = std::min(min_angle, angle);
+                    max_angle = std::max(max_angle, angle);
+                }
+                REQUIRE(min_angle == max_angle);
+                REQUIRE(min_angle == 0);
+            }
+        }
+
+    }
+
+
+    //TODO: compare with mainline slic3r
+    //GIVEN("tooth")
+    //{
+    //    ExPolygon expolygon;
+    //    expolygon.contour = Slic3r::Polygon{ Points{
+    //        Point::new_scale(0.86526705, 1.4509841), Point::new_scale(0.57696039, 1.8637021), 
+    //        Point::new_scale(0.4502297, 2.5569978), Point::new_scale(0.45626199, 3.2965596), 
+    //        Point::new_scale(1.1218851, 3.3049455), Point::new_scale(0.96681072, 2.8243202), 
+    //        Point::new_scale(0.86328971, 2.2056997), Point::new_scale(0.85367905, 1.7790778)
+    //    } };
+    //    expolygon.contour.make_counter_clockwise();
+    //    WHEN("creating the medial axis"){
+    //        Polylines res;
+    //        expolygon.medial_axis(scale_(1), scale_(0.25), &res);
+    //        THEN("medial axis of a tooth is two lines"){
+    //            REQUIRE(res.size() == 2);
+    //            THEN("medial axis has reasonable length") {
+    //                REQUIRE(res[0].length() >= scale_(1.4) - SCALED_EPSILON);
+    //                REQUIRE(res[1].length() >= scale_(1.4) - SCALED_EPSILON);
+    //                // TODO: check if min width is < 0.3 and max width is > 0.6 (min($res->[0]->width.front, $res->[0]->width.back) # problem: can't have access to width
+    //                //TODO: now i have access! correct it!
+    //            }
+    //        }
+    //    }
+    //} 
+
+    GIVEN("Anchor & Tapers")
+    {
+        ExPolygon tooth;
+        tooth.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(0,0), Point::new_scale(10,0), Point::new_scale(10,1.2), Point::new_scale(0,1.2)
+        } };
+        tooth.contour.make_counter_clockwise();
+        ExPolygon base_part;
+        base_part.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(0,-3), Point::new_scale(0,3), Point::new_scale(-2,3), Point::new_scale(-2,-3)
+        } };
+        base_part.contour.make_counter_clockwise();
+        //expolygon.contour = Slic3r::Polygon{ Points{
+        //    //Point::new_scale(0, 13), Point::new_scale(-1, 13), Point::new_scale(-1, 0), Point::new_scale(0.0,0.0),
+        //    Point::new_scale(0,0.2), Point::new_scale(3,0.2), Point::new_scale(3,0.4), Point::new_scale(0,0.4),
+        //    Point::new_scale(0,1), Point::new_scale(3,1), Point::new_scale(3,1.3), Point::new_scale(0,1.3),
+        //    Point::new_scale(0,2), Point::new_scale(3,2), Point::new_scale(3,2.4), Point::new_scale(0,2.4),
+        //    Point::new_scale(0,3), Point::new_scale(3,3), Point::new_scale(3,3.5), Point::new_scale(0,3.5),
+        //    Point::new_scale(0,4), Point::new_scale(3,4), Point::new_scale(3,4.6), Point::new_scale(0,4.6),
+        //    Point::new_scale(0,5), Point::new_scale(3,5), Point::new_scale(3,5.7), Point::new_scale(0,5.7),
+        //    Point::new_scale(0,6), Point::new_scale(3,6), Point::new_scale(3,6.8), Point::new_scale(0,6.8),
+        //    Point::new_scale(0,7.5), Point::new_scale(3,7.5), Point::new_scale(3,8.4), Point::new_scale(0,8.4),
+        //    Point::new_scale(0,9), Point::new_scale(3,9), Point::new_scale(3,10), Point::new_scale(0,10),
+        //    Point::new_scale(0,11), Point::new_scale(3,11), Point::new_scale(3,12.2), Point::new_scale(0,12.2),
+        //} };
+        WHEN("1 nozzle, 0.2 layer height") {
+            const coord_t nozzle_diam = scale_(1);
+            ExPolygon anchor = union_ex(ExPolygons{ tooth }, intersection_ex(ExPolygons{ base_part }, offset_ex(tooth, nozzle_diam / 2)), true)[0];
+            ThickPolylines res;
+            //expolygon.medial_axis(scale_(1), scale_(0.5), &res);
+            Slic3r::MedialAxis ma(tooth, nozzle_diam * 2, nozzle_diam/3, scale_(0.2));
+            ma.use_bounds(anchor)
+                .use_min_real_width(nozzle_diam)
+                .use_tapers(0.25*nozzle_diam);
+            ma.build(res);
+            THEN("medial axis of a simple line is one line") {
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has the length of the line + the length of the anchor") {
+                    std::cout << res[0].length() << "\n";
+                    REQUIRE(std::abs(res[0].length() - scale_(10.5)) < SCALED_EPSILON);
+                }
+                THEN("medial axis has the line width as max width") {
+                    double max_width = 0;
+                    for (coordf_t width : res[0].width) max_width = std::max(max_width, width);
+                    REQUIRE(std::abs(max_width - scale_(1.2)) < SCALED_EPSILON);
+                }
+                //compute the length of the tapers
+                THEN("medial axis has good tapers length") {
+                    int l1 = 0;
+                    for (size_t idx = 0; idx < res[0].width.size() - 1 && res[0].width[idx] - nozzle_diam < SCALED_EPSILON; ++idx)
+                        l1 += res[0].lines()[idx].length();
+                    int l2 = 0;
+                    for (size_t idx = res[0].width.size() - 1; idx > 0 && res[0].width[idx] - nozzle_diam < SCALED_EPSILON; --idx)
+                        l2 += res[0].lines()[idx - 1].length();
+                    REQUIRE(std::abs(l1 - l2) < SCALED_EPSILON);
+                    REQUIRE(std::abs(l1 - scale_(0.25 - 0.1)) < SCALED_EPSILON);
+                }
+            }
+        }
+
+        WHEN("1.2 nozzle, 0.6 layer height") {
+            const coord_t nozzle_diam = scale_(1.2);
+            ExPolygon anchor = union_ex(ExPolygons{ tooth }, intersection_ex(ExPolygons{ base_part }, offset_ex(tooth, nozzle_diam / 4)), true)[0];
+            ThickPolylines res;
+            //expolygon.medial_axis(scale_(1), scale_(0.5), &res);
+            Slic3r::MedialAxis ma(tooth, nozzle_diam * 2, nozzle_diam/3, scale_(0.6));
+            ma.use_bounds(anchor)
+                .use_min_real_width(nozzle_diam)
+                .use_tapers(1.0*nozzle_diam);
+            ma.build(res);
+            THEN("medial axis of a simple line is one line") {
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has the length of the line + the length of the anchor") {
+                    //0.3 because it's offseted by nozzle_diam / 4
+                    REQUIRE(std::abs(res[0].length() - scale_(10.3)) < SCALED_EPSILON);
+                }
+                THEN("medial axis can'ty have a line width below Flow::new_from_spacing(nozzle_diam).width") {
+                    double max_width = 0;
+                    for (coordf_t width : res[0].width) max_width = std::max(max_width, width);
+                    double min_width = Flow::new_from_spacing(float(unscale_(nozzle_diam)), float(unscale_(nozzle_diam)), 0.6f, false).scaled_width();
+                    REQUIRE(std::abs(max_width - min_width) < SCALED_EPSILON);
+                    REQUIRE(std::abs(max_width - nozzle_diam)  > SCALED_EPSILON);
+                }
+                //compute the length of the tapers
+                THEN("medial axis has a 45� taper and a shorter one") {
+                    int l1 = 0;
+                    for (size_t idx = 0; idx < res[0].width.size() - 1 && res[0].width[idx] - scale_(1.2) < SCALED_EPSILON; ++idx)
+                        l1 += res[0].lines()[idx].length();
+                    int l2 = 0;
+                    for (size_t idx = res[0].width.size() - 1; idx > 0 && res[0].width[idx] - scale_(1.2) < SCALED_EPSILON; --idx)
+                        l2 += res[0].lines()[idx - 1].length();
+                    //here the taper is limited by the 0-width spacing
+                    double min_width = Flow::new_from_spacing(float(unscale_(nozzle_diam)), float(unscale_(nozzle_diam)), 0.6f, false).scaled_width();
+                    REQUIRE(std::abs(l1 - l2) < SCALED_EPSILON);
+                    REQUIRE(l1 < scale_(0.6));
+                    REQUIRE(l1  > scale_(0.4));
+                }
+            }
+        }
+
+    }
+
+    GIVEN("1� rotated tooths")
+    {
+
+    }
+
+    GIVEN("narrow trapezoid")
+    {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(100, 100),
+            Point::new_scale(120, 100),
+            Point::new_scale(112, 200),
+            Point::new_scale(108, 200)
+        } };
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(20), scale_(0.5), &res);
+            THEN("medial axis of a narrow trapezoid is a single line"){
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has reasonable length") {
+                    REQUIRE(res[0].length() >= scale_(200 - 100 - (120 - 100)) - SCALED_EPSILON);
+                }
+            }
+        }
+    }
+
+    GIVEN("L shape")
+    {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(100, 100),
+            Point::new_scale(120, 100),
+            Point::new_scale(120, 180),
+            Point::new_scale(200, 180),
+            Point::new_scale(200, 200),
+            Point::new_scale(100, 200)
+        } };
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(20), scale_(0.5), &res);
+            THEN("medial axis of a L shape is a single line"){
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has reasonable length") {
+                    // 20 is the thickness of the expolygon, which is subtracted from the ends
+                    REQUIRE(res[0].length() + 20 > scale_(80 * 2) - SCALED_EPSILON);
+                    REQUIRE(res[0].length() + 20 < scale_(100 * 2) + SCALED_EPSILON);
+                }
+            }
+        }
+    }
+
+    GIVEN("shape"){
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point{ -203064906, -51459966 }, Point{ -219312231, -51459966 }, Point{ -219335477, -51459962 }, Point{ -219376095, -51459962 }, Point{ -219412047, -51459966 },
+            Point{ -219572094, -51459966 }, Point{ -219624814, -51459962 }, Point{ -219642183, -51459962 }, Point{ -219656665, -51459966 }, Point{ -220815482, -51459966 },
+            Point{ -220815482, -37738966 }, Point{ -221117540, -37738966 }, Point{ -221117540, -51762024 }, Point{ -203064906, -51762024 },
+        } };
+        WHEN("creating the medial axis"){
+            Polylines polylines;
+            expolygon.medial_axis(819998, 102499.75, &polylines);
+            double perimeter_len = expolygon.contour.split_at_first_point().length();
+            THEN("medial axis has reasonable length"){
+                double polyline_length = 0;
+                for (Slic3r::Polyline &poly : polylines) polyline_length += poly.length();
+                REQUIRE(polyline_length > perimeter_len * 3. / 8. - SCALED_EPSILON);
+            }
+        }
+    }
+
+    GIVEN("narrow triangle")
+    {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point::new_scale(50, 100),
+            Point::new_scale(1000, 102),
+            Point::new_scale(50, 104)
+        } };
+
+        WHEN("creating the medial axis"){
+            Polylines res;
+            expolygon.medial_axis(scale_(4), scale_(0.5), &res);
+            THEN("medial axis of a narrow triangle is a single line"){
+                REQUIRE(res.size() == 1);
+                THEN("medial axis has reasonable length") {
+                    REQUIRE(res[0].length() > scale_(200 - 100 - (120 - 100)) - SCALED_EPSILON);
+                }
+            }
+        }
+    }
+
+    GIVEN("GH #2474")
+    {
+        ExPolygon expolygon;
+        expolygon.contour = Slic3r::Polygon{ Points{
+            Point{91294454, 31032190},
+            Point{11294481, 31032190},
+            Point{11294481, 29967810},
+            Point{44969182, 29967810},
+            Point{89909960, 29967808},
+            Point{91294454, 29967808}
+        } };
+        WHEN("creating the medial axis") {
+            Polylines res;
+            expolygon.medial_axis(1871238, 500000, &res);
+            THEN("medial axis is a single polyline") {
+                REQUIRE(res.size() == 1);
+                Slic3r::Polyline polyline = res[0];
+                THEN("medial axis is horizontal and is centered") {
+                    double sum = 0;
+                    for (Line &l : polyline.lines()) sum += std::abs(l.b.y() - l.a.y());
+                    coord_t expected_y = expolygon.contour.bounding_box().center().y();
+                    REQUIRE((sum / polyline.size()) - expected_y < SCALED_EPSILON);
+                }
+
+                // order polyline from left to right
+                if (polyline.first_point().x() > polyline.last_point().x()) polyline.reverse();
+
+                THEN("expected x_min & x_max") {
+                    BoundingBox polyline_bb = polyline.bounding_box();
+                    REQUIRE(polyline.first_point().x() == polyline_bb.min.x());
+                    REQUIRE(polyline.last_point().x() == polyline_bb.max.x());
+                }
+
+                THEN("medial axis is not self-overlapping") {
+                    //TODO
+                    //REQUIRE(polyline.first_point().x() == polyline_bb.x_min());
+                    std::vector<coord_t> all_x;
+                    for (Point &p : polyline.points) all_x.push_back(p.x());
+                    std::vector<coord_t> sorted_x{ all_x };
+                    std::sort(sorted_x.begin(), sorted_x.end());
+                    for (size_t i = 0; i < all_x.size(); i++) {
+                        REQUIRE(all_x[i] == sorted_x[i]);
+                    }
+                }
+            }
+        }
+    }
+
+}
diff --git a/t/adaptive_width.t b/t/adaptive_width.t
index 7a0baa7527..14ffd3dd4c 100644
--- a/t/adaptive_width.t
+++ b/t/adaptive_width.t
@@ -32,6 +32,7 @@ sub scale_points (@) { map [scale $_->[X], scale $_->[Y]], @_ }
             expolygon       => $expolygon,
         ));
         my $config          = Slic3r::Config::Full->new;
+        $config->set('thin_walls_overlap',0);
         my $loops           = Slic3r::ExtrusionPath::Collection->new;
         my $gap_fill        = Slic3r::ExtrusionPath::Collection->new;
         my $fill_surfaces   = Slic3r::Surface::Collection->new;
diff --git a/t/thin.t b/t/thin.t
index 2d256d2864..9bd4e9aa0f 100644
--- a/t/thin.t
+++ b/t/thin.t
@@ -1,4 +1,4 @@
-use Test::More tests => 23;
+use Test::More tests => 29;
 use strict;
 use warnings;
 
@@ -102,10 +102,44 @@ if (0) {
     is scalar(@$res), 1, 'medial axis of a semicircumference is a single line';
     
     # check whether turns are all CCW or all CW
-    my @lines = @{$res->[0]->lines};
+    my @all_lines = @{$res->[0]->lines};
+    # remove lines that are near the end.
+    my @lines = grep($_->a->y >= 1578184 || $_->b->y >= 1578184, @all_lines);
     my @angles = map { $lines[$_-1]->ccw($lines[$_]->b) } 1..$#lines;
     ok !!(none { $_ < 0 } @angles) || (none { $_ > 0 } @angles),
-        'all medial axis segments of a semicircumference have the same orientation';
+        'all medial axis segments of a semicircumference have the same orientation (but the 2 end points)';
+}
+
+{
+    my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
+        [4.3, 4], [4.3, 0], [4,0], [4,4], [0,4], [0,4.5], [4,4.5], [4,10], [4.3,10], [4.3, 4.5],
+        [6, 4.5], [6,10], [6.2,10], [6.2,4.5], [10,4.5], [10,4], [6.2,4], [6.2,0], [6, 0], [6, 4],
+    ));
+    $expolygon->contour->make_counter_clockwise();
+    my $res = $expolygon->medial_axis(scale 0.55, scale 0.25);
+    is scalar(@$res), 2, 'medial axis of a (bit too narrow) french cross is two lines';
+    ok unscale($res->[0]->length) >= (9.9) - epsilon, 'medial axis has reasonable length';
+    ok unscale($res->[1]->length) >= (9.9) - epsilon, 'medial axis has reasonable length';
+    
+    my @lines1 = @{$res->[0]->lines};
+    my @angles1 = map { $lines1[$_-1]->ccw($lines1[$_]->b) } 1..$#lines1;
+    my @lines2 = @{$res->[1]->lines};
+    my @angles2 = map { $lines2[$_-1]->ccw($lines2[$_]->b) } 1..$#lines2;
+    my @angles = (@angles1, @angles2);
+    ok !!(none { $_ != 0 } @angles),
+        'medial axis of a (bit too narrow) french cross is two lines has only strait lines';
+}
+
+{
+    my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
+        [0.86526705,1.4509841], [0.57696039,1.8637021], [0.4502297,2.5569978], [0.45626199,3.2965596], [1.1218851,3.3049455], [0.96681072,2.8243202], [0.86328971,2.2056997], [0.85367905,1.7790778],
+    ));
+    $expolygon->contour->make_counter_clockwise();
+    my $res = $expolygon->medial_axis(scale 1, scale 0.25);
+    is scalar(@$res), 1, 'medial axis of a (bit too narrow) french cross is two lines';
+    ok unscale($res->[0]->length) >= (1.4) - epsilon, 'medial axis has reasonable length';
+    # TODO: check if min width is < 0.3 and max width is > 0.6 (min($res->[0]->width.front, $res->[0]->width.back) # problem: can't have access to width
+    
 }
 
 {
diff --git a/xs/MANIFEST b/xs/MANIFEST
index a2b1940ab5..c144c822eb 100644
--- a/xs/MANIFEST
+++ b/xs/MANIFEST
@@ -121,6 +121,8 @@ src/libslic3r/libslic3r.h
 src/libslic3r/Line.cpp
 src/libslic3r/Line.hpp
 src/libslic3r/Log.hpp
+src/libslic3r/MedialAxis.cpp
+src/libslic3r/MedialAxis.hpp
 src/libslic3r/Model.cpp
 src/libslic3r/Model.hpp
 src/libslic3r/MotionPlanner.cpp
diff --git a/xs/src/libslic3r/ExPolygon.cpp b/xs/src/libslic3r/ExPolygon.cpp
index 4043b7dfe0..3767c6adad 100644
--- a/xs/src/libslic3r/ExPolygon.cpp
+++ b/xs/src/libslic3r/ExPolygon.cpp
@@ -1,5 +1,6 @@
 #include "BoundingBox.hpp"
 #include "ExPolygon.hpp"
+#include "MedialAxis.hpp"
 #include "Geometry.hpp"
 #include "Polygon.hpp"
 #include "Line.hpp"
@@ -192,6 +193,30 @@ ExPolygon::simplify(double tolerance, ExPolygons* expolygons) const
     expolygons->insert(expolygons->end(), ep.begin(), ep.end());
 }
 
+/// remove point that are at SCALED_EPSILON * 2 distance.
+//simplier than simplify
+void
+ExPolygon::remove_point_too_near(const coord_t tolerance) {
+    size_t id = 1;
+    while (id < this->contour.points.size() - 1) {
+        double newdist = min(this->contour.points[id].distance_to(this->contour.points[id - 1])
+            , this->contour.points[id].distance_to(this->contour.points[id + 1]));
+        if (newdist < (double)tolerance) {
+            this->contour.points.erase(this->contour.points.begin() + id);
+            newdist = this->contour.points[id].distance_to(this->contour.points[id - 1]);
+        }
+        //go to next one
+        //if you removed a point, it check if the next one isn't too near from the previous one.
+        // if not, it byepass it.
+        if (newdist > (double)tolerance) {
+            ++id;
+        }
+    }
+    if (this->contour.points.front().distance_to(this->contour.points.back()) < tolerance) {
+        this->contour.points.erase(this->contour.points.end() -1);
+    }
+}
+
 void
 ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_width, ThickPolylines* polylines) const
 {
@@ -399,42 +424,23 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
                 best_dot = other_dot;
             }
         }
-        if (best_candidate != nullptr) {
-
-            polyline.points.insert(polyline.points.end(), best_candidate->points.begin() + 1, best_candidate->points.end());
-            polyline.width.insert(polyline.width.end(), best_candidate->width.begin(), best_candidate->width.end());
-            polyline.endpoints.second = best_candidate->endpoints.second;
-            assert(polyline.width.size() == polyline.points.size()*2 - 2);
-                
-            pp.erase(pp.begin() + best_idx);
+        //go to next one
+        //if you removed a point, it check if the next one isn't too near from the previous one.
+        // if not, it byepass it.
+        if (newdist > (double)tolerance) {
+            ++id;
         }
     }
-    
-    for (size_t i = 0; i < pp.size(); ++i) {
-        ThickPolyline& polyline = pp[i];
-
-        //  remove too short polylines
-        // (we can't do this check before endpoints extension and clipping because we don't
-        // know how long will the endpoints be extended since it depends on polygon thickness
-        // which is variable - extension will be <= max_width/2 on each side)
-        if ((polyline.endpoints.first || polyline.endpoints.second)
-            && polyline.length() < max_w * 2) {
-            pp.erase(pp.begin() + i);
-            --i;
-            continue;
-        }
-
+    if (this->contour.points.front().distance_to(this->contour.points.back()) < tolerance) {
+        this->contour.points.erase(this->contour.points.end() -1);
     }
-
-    
-    polylines->insert(polylines->end(), pp.begin(), pp.end());
 }
 
 void
 ExPolygon::medial_axis(double max_width, double min_width, Polylines* polylines) const
 {
     ThickPolylines tp;
-    this->medial_axis(*this, max_width, min_width, &tp);
+    MedialAxis{*this, coord_t(max_width), coord_t(min_width), coord_t(max_width / 2.0)}.build(tp);
     polylines->insert(polylines->end(), tp.begin(), tp.end());
 }
 
diff --git a/xs/src/libslic3r/ExPolygon.hpp b/xs/src/libslic3r/ExPolygon.hpp
index 722c2392a6..2a2ae1a682 100644
--- a/xs/src/libslic3r/ExPolygon.hpp
+++ b/xs/src/libslic3r/ExPolygon.hpp
@@ -44,7 +44,8 @@ class ExPolygon
     Polygons simplify_p(double tolerance) const;
     ExPolygons simplify(double tolerance) const;
     void simplify(double tolerance, ExPolygons* expolygons) const;
-    void medial_axis(const ExPolygon &bounds, double max_width, double min_width, ThickPolylines* polylines) const;
+    void remove_point_too_near(const coord_t tolerance);
+    ///TODO: remove that function. It's kept only because of perl tests.
     void medial_axis(double max_width, double min_width, Polylines* polylines) const;
     void get_trapezoids(Polygons* polygons) const;
     void get_trapezoids(Polygons* polygons, double angle) const;
diff --git a/xs/src/libslic3r/Geometry.cpp b/xs/src/libslic3r/Geometry.cpp
index b972e7e176..c71bc4bb6a 100644
--- a/xs/src/libslic3r/Geometry.cpp
+++ b/xs/src/libslic3r/Geometry.cpp
@@ -648,261 +648,5 @@ arrange(size_t total_parts, const Pointf &part_size, coordf_t dist, const Boundi
     return true;
 }
 
-void
-MedialAxis::build(ThickPolylines* polylines)
-{
-    construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd);
-    
-    /*
-    // DEBUG: dump all Voronoi edges
-    {
-        SVG svg("voronoi.svg");
-        svg.draw(*this->expolygon);
-        for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
-            if (edge->is_infinite()) continue;
-            
-            ThickPolyline polyline;
-            polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
-            polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
-            polyline.width.push_back(this->max_width);
-            polyline.width.push_back(this->max_width);
-            polylines->push_back(polyline);
-            
-            svg.draw(polyline);
-        }
-        svg.Close();
-        return;
-    }
-    */
-    
-    // collect valid edges (i.e. prune those not belonging to MAT)
-    // note: this keeps twins, so it inserts twice the number of the valid edges
-    this->valid_edges.clear();
-    {
-        std::set<const VD::edge_type*> seen_edges;
-        for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
-            // if we only process segments representing closed loops, none if the
-            // infinite edges (if any) would be part of our MAT anyway
-            if (edge->is_secondary() || edge->is_infinite()) continue;
-        
-            // don't re-validate twins
-            if (seen_edges.find(&*edge) != seen_edges.end()) continue;  // TODO: is this needed?
-            seen_edges.insert(&*edge);
-            seen_edges.insert(edge->twin());
-            
-            if (!this->validate_edge(&*edge)) continue;
-            this->valid_edges.insert(&*edge);
-            this->valid_edges.insert(edge->twin());
-        }
-    }
-    this->edges = this->valid_edges;
     
-    // iterate through the valid edges to build polylines
-    while (!this->edges.empty()) {
-        const VD::edge_type* edge = *this->edges.begin();
-        
-        // start a polyline
-        ThickPolyline polyline;
-        polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
-        polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
-        polyline.width.push_back(this->thickness[edge].first);
-        polyline.width.push_back(this->thickness[edge].second);
-        
-        // remove this edge and its twin from the available edges
-        (void)this->edges.erase(edge);
-        (void)this->edges.erase(edge->twin());
-        
-        // get next points
-        this->process_edge_neighbors(edge, &polyline);
-        
-        // get previous points
-        {
-            ThickPolyline rpolyline;
-            this->process_edge_neighbors(edge->twin(), &rpolyline);
-            polyline.points.insert(polyline.points.begin(), rpolyline.points.rbegin(), rpolyline.points.rend());
-            polyline.width.insert(polyline.width.begin(), rpolyline.width.rbegin(), rpolyline.width.rend());
-            polyline.endpoints.first = rpolyline.endpoints.second;
-        }
-        
-        assert(polyline.width.size() == polyline.points.size()*2 - 2);
-        
-        // prevent loop endpoints from being extended
-        if (polyline.first_point().coincides_with(polyline.last_point())) {
-            polyline.endpoints.first = false;
-            polyline.endpoints.second = false;
-        }
-        
-        // append polyline to result
-        polylines->push_back(polyline);
-    }
-}
-
-void
-MedialAxis::build(Polylines* polylines)
-{
-    ThickPolylines tp;
-    this->build(&tp);
-    polylines->insert(polylines->end(), tp.begin(), tp.end());
-}
-
-void
-MedialAxis::process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline)
-{
-    while (true) {
-        // Since rot_next() works on the edge starting point but we want
-        // to find neighbors on the ending point, we just swap edge with
-        // its twin.
-        const VD::edge_type* twin = edge->twin();
-    
-        // count neighbors for this edge
-        std::vector<const VD::edge_type*> neighbors;
-        for (const VD::edge_type* neighbor = twin->rot_next(); neighbor != twin;
-            neighbor = neighbor->rot_next()) {
-            if (this->valid_edges.count(neighbor) > 0) neighbors.push_back(neighbor);
-        }
-    
-        // if we have a single neighbor then we can continue recursively
-        if (neighbors.size() == 1) {
-            const VD::edge_type* neighbor = neighbors.front();
-            
-            // break if this is a closed loop
-            if (this->edges.count(neighbor) == 0) return;
-            
-            Point new_point(neighbor->vertex1()->x(), neighbor->vertex1()->y());
-            polyline->points.push_back(new_point);
-            polyline->width.push_back(this->thickness[neighbor].first);
-            polyline->width.push_back(this->thickness[neighbor].second);
-            (void)this->edges.erase(neighbor);
-            (void)this->edges.erase(neighbor->twin());
-            edge = neighbor;
-        } else if (neighbors.size() == 0) {
-            polyline->endpoints.second = true;
-            return;
-        } else {
-            // T-shaped or star-shaped joint
-            return;
-        }
-    }
-}
-
-bool
-MedialAxis::validate_edge(const VD::edge_type* edge)
-{
-    // prevent overflows and detect almost-infinite edges
-    if (std::abs(edge->vertex0()->x()) > (double)MAX_COORD
-     || std::abs(edge->vertex0()->y()) > (double)MAX_COORD
-     || std::abs(edge->vertex1()->x()) > (double)MAX_COORD
-     || std::abs(edge->vertex1()->y()) > (double)MAX_COORD)
-        return false;
-    
-    // construct the line representing this edge of the Voronoi diagram
-    const Line line(
-        Point( edge->vertex0()->x(), edge->vertex0()->y() ),
-        Point( edge->vertex1()->x(), edge->vertex1()->y() )
-    );
-    
-    // discard edge if it lies outside the supplied shape
-    // this could maybe be optimized (checking inclusion of the endpoints
-    // might give false positives as they might belong to the contour itself)
-    if (this->expolygon != NULL) {
-        if (line.a.coincides_with(line.b)) {
-            // in this case, contains(line) returns a false positive
-            if (!this->expolygon->contains(line.a)) return false;
-        } else {
-            if (!this->expolygon->contains(line)) return false;
-        }
-    }
-    
-    // retrieve the original line segments which generated the edge we're checking
-    const VD::cell_type* cell_l = edge->cell();
-    const VD::cell_type* cell_r = edge->twin()->cell();
-    const Line &segment_l = this->retrieve_segment(cell_l);
-    const Line &segment_r = this->retrieve_segment(cell_r);
-    
-    /*
-    SVG svg("edge.svg");
-    svg.draw(*this->expolygon);
-    svg.draw(line);
-    svg.draw(segment_l, "red");
-    svg.draw(segment_r, "blue");
-    svg.Close();
-    */
-    
-    /*  Calculate thickness of the cross-section at both the endpoints of this edge.
-        Our Voronoi edge is part of a CCW sequence going around its Voronoi cell 
-        located on the left side. (segment_l).
-        This edge's twin goes around segment_r. Thus, segment_r is 
-        oriented in the same direction as our main edge, and segment_l is oriented
-        in the same direction as our twin edge.
-        We used to only consider the (half-)distances to segment_r, and that works
-        whenever segment_l and segment_r are almost specular and facing. However, 
-        at curves they are staggered and they only face for a very little length
-        (our very short edge represents such visibility).
-        Both w0 and w1 can be calculated either towards cell_l or cell_r with equal
-        results by Voronoi definition.
-        When cell_l or cell_r don't refer to the segment but only to an endpoint, we
-        calculate the distance to that endpoint instead.  */
-    
-    coordf_t w0 = cell_r->contains_segment()
-        ? line.a.distance_to(segment_r)*2
-        : line.a.distance_to(this->retrieve_endpoint(cell_r))*2;
-    
-    coordf_t w1 = cell_l->contains_segment()
-        ? line.b.distance_to(segment_l)*2
-        : line.b.distance_to(this->retrieve_endpoint(cell_l))*2;
-    
-    if (cell_l->contains_segment() && cell_r->contains_segment()) {
-        // calculate the relative angle between the two boundary segments
-        double angle = fabs(segment_r.orientation() - segment_l.orientation());
-        if (angle > PI) angle = 2*PI - angle;
-        assert(angle >= 0 && angle <= PI);
-        
-        // fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
-        // we're interested only in segments close to the second case (facing segments)
-        // so we allow some tolerance.
-        // this filter ensures that we're dealing with a narrow/oriented area (longer than thick)
-        // we don't run it on edges not generated by two segments (thus generated by one segment
-        // and the endpoint of another segment), since their orientation would not be meaningful
-        if (PI - angle > PI/8) {
-            // angle is not narrow enough
-            
-            // only apply this filter to segments that are not too short otherwise their 
-            // angle could possibly be not meaningful
-            if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON || line.length() >= this->min_width)
-                return false;
-        }
-    } else {
-        if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON)
-            return false;
-    }
-    
-    if (w0 < this->min_width && w1 < this->min_width)
-        return false;
-    
-    if (w0 > this->max_width && w1 > this->max_width)
-        return false;
-    
-    this->thickness[edge]         = std::make_pair(w0, w1);
-    this->thickness[edge->twin()] = std::make_pair(w1, w0);
-    
-    return true;
-}
-
-const Line&
-MedialAxis::retrieve_segment(const VD::cell_type* cell) const
-{
-    return this->lines[cell->source_index()];
-}
-
-const Point&
-MedialAxis::retrieve_endpoint(const VD::cell_type* cell) const
-{
-    const Line& line = this->retrieve_segment(cell);
-    if (cell->source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) {
-        return line.a;
-    } else {
-        return line.b;
-    }
-}
-
 } }
diff --git a/xs/src/libslic3r/Geometry.hpp b/xs/src/libslic3r/Geometry.hpp
index 9d4579389c..8dce56a2fb 100644
--- a/xs/src/libslic3r/Geometry.hpp
+++ b/xs/src/libslic3r/Geometry.hpp
@@ -7,9 +7,6 @@
 #include "Polygon.hpp"
 #include "Polyline.hpp"
 
-#include "boost/polygon/voronoi.hpp"
-using boost::polygon::voronoi_builder;
-using boost::polygon::voronoi_diagram;
 
 namespace Slic3r { namespace Geometry {
 
@@ -45,28 +42,6 @@ bool arrange(
     // output
     Pointfs &positions);
 
-class MedialAxis {
-    public:
-    Lines lines;
-    const ExPolygon* expolygon;
-    double max_width;
-    double min_width;
-    MedialAxis(double _max_width, double _min_width, const ExPolygon* _expolygon = NULL)
-        : expolygon(_expolygon), max_width(_max_width), min_width(_min_width) {};
-    void build(ThickPolylines* polylines);
-    void build(Polylines* polylines);
-    
-    private:
-    typedef voronoi_diagram<double> VD;
-    VD vd;
-    std::set<const VD::edge_type*> edges, valid_edges;
-    std::map<const VD::edge_type*, std::pair<coordf_t,coordf_t> > thickness;
-    void process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline);
-    bool validate_edge(const VD::edge_type* edge);
-    const Line& retrieve_segment(const VD::cell_type* cell) const;
-    const Point& retrieve_endpoint(const VD::cell_type* cell) const;
-};
-
 } }
 
 #endif
diff --git a/xs/src/libslic3r/Layer.cpp b/xs/src/libslic3r/Layer.cpp
index 4b4c3d8533..9c48163378 100644
--- a/xs/src/libslic3r/Layer.cpp
+++ b/xs/src/libslic3r/Layer.cpp
@@ -189,6 +189,8 @@ Layer::make_perimeters()
                 && config.overhangs         == other_config.overhangs
                 && config.serialize("perimeter_extrusion_width").compare(other_config.serialize("perimeter_extrusion_width")) == 0
                 && config.thin_walls        == other_config.thin_walls
+                && config.thin_walls_min_width      == other_config.thin_walls_min_width
+                && config.thin_walls_overlap        == other_config.thin_walls_overlap
                 && config.external_perimeters_first == other_config.external_perimeters_first) {
                 layerms.push_back(other_layerm);
                 done.insert(it - this->regions.begin());
diff --git a/xs/src/libslic3r/LayerHeightSpline.cpp b/xs/src/libslic3r/LayerHeightSpline.cpp
index 0fadf35136..7fe114b2bb 100644
--- a/xs/src/libslic3r/LayerHeightSpline.cpp
+++ b/xs/src/libslic3r/LayerHeightSpline.cpp
@@ -84,7 +84,7 @@ bool LayerHeightSpline::updateLayerHeights(std::vector<coordf_t> heights)
         this->_layer_heights = heights;
         result = this->_updateBSpline();
     }else{
-        std::cerr << "Unable to update layer heights. You provided " << heights.size() << " layers, but " << this->_layers.size()-1 << " expected" << std::endl;
+        std::cout << "Unable to update layer heights. You provided " << heights.size() << " layers, but " << this->_layers.size()-1 << " expected" << std::endl;
     }
 
     this->_layers_updated = false;
@@ -193,7 +193,7 @@ bool LayerHeightSpline::_updateBSpline()
         result = true;
     } else {
         result = false;
-        std::cerr << "Spline setup failed." << std::endl;
+        std::cout << "Spline setup failed." << std::endl;
     }
 
     this->_is_valid = result;
diff --git a/xs/src/libslic3r/MedialAxis.cpp b/xs/src/libslic3r/MedialAxis.cpp
new file mode 100644
index 0000000000..2e18046279
--- /dev/null
+++ b/xs/src/libslic3r/MedialAxis.cpp
@@ -0,0 +1,1611 @@
+#include "MedialAxis.hpp"
+#include "BoundingBox.hpp"
+#include "ExPolygon.hpp"
+#include "Geometry.hpp"
+#include "Polygon.hpp"
+#include "Line.hpp"
+#include "ClipperUtils.hpp"
+#include "SVG.hpp"
+#include "polypartition.h"
+#include "poly2tri/poly2tri.h"
+#include <algorithm>
+#include <cassert>
+#include <list>
+
+namespace Slic3r {
+    
+void
+MedialAxis::build(Polylines &polylines)
+{
+    ThickPolylines tp;
+    this->build(tp);
+    polylines.insert(polylines.end(), tp.begin(), tp.end());
+}
+
+void
+MedialAxis::polyline_from_voronoi(const Lines& voronoi_edges, ThickPolylines* polylines)
+{
+    this->lines = voronoi_edges;
+    construct_voronoi(lines.begin(), lines.end(), &this->vd);
+
+    /*
+    // DEBUG: dump all Voronoi edges
+    {
+        for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
+            if (edge->is_infinite()) continue;
+            
+            ThickPolyline polyline;
+            polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
+            polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
+            polylines->push_back(polyline);
+        }
+        return;
+    }
+    */
+    
+    typedef const VD::edge_type   edge_t;
+    
+    // collect valid edges (i.e. prune those not belonging to MAT)
+    // note: this keeps twins, so it inserts twice the number of the valid edges
+    this->valid_edges.clear();
+    {
+        std::set<const VD::edge_type*> seen_edges;
+        for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
+            // if we only process segments representing closed loops, none if the
+            // infinite edges (if any) would be part of our MAT anyway
+            if (edge->is_secondary() || edge->is_infinite()) continue;
+        
+            // don't re-validate twins
+            if (seen_edges.find(&*edge) != seen_edges.end()) continue;  // TODO: is this needed?
+            seen_edges.insert(&*edge);
+            seen_edges.insert(edge->twin());
+            
+            if (!this->validate_edge(&*edge)) continue;
+            this->valid_edges.insert(&*edge);
+            this->valid_edges.insert(edge->twin());
+        }
+    }
+    this->edges = this->valid_edges;
+    
+    // iterate through the valid edges to build polylines
+    while (!this->edges.empty()) {
+        const edge_t* edge = *this->edges.begin();
+        
+        // start a polyline
+        ThickPolyline polyline;
+        polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
+        polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
+        polyline.width.push_back(this->thickness[edge].first);
+        polyline.width.push_back(this->thickness[edge].second);
+        
+        // remove this edge and its twin from the available edges
+        (void)this->edges.erase(edge);
+        (void)this->edges.erase(edge->twin());
+        
+        // get next points
+        this->process_edge_neighbors(edge, &polyline);
+        
+        // get previous points
+        {
+            ThickPolyline rpolyline;
+            this->process_edge_neighbors(edge->twin(), &rpolyline);
+            polyline.points.insert(polyline.points.begin(), rpolyline.points.rbegin(), rpolyline.points.rend());
+            polyline.width.insert(polyline.width.begin(), rpolyline.width.rbegin(), rpolyline.width.rend());
+            polyline.endpoints.first = rpolyline.endpoints.second;
+        }
+        
+        assert(polyline.width.size() == polyline.points.size());
+        // if loop, set endpoints to false
+        // prevent loop endpoints from being extended
+        if (polyline.first_point().coincides_with(polyline.last_point())) {
+            polyline.endpoints.first = false;
+            polyline.endpoints.second = false;
+        }
+        
+        // append polyline to result
+        polylines->push_back(polyline);
+    }
+
+    #ifdef SLIC3R_DEBUG
+    {
+        static int iRun = 0;
+        dump_voronoi_to_svg(this->lines, this->vd, polylines, debug_out_path("MedialAxis-%d.svg", iRun ++).c_str());
+        printf("Thick lines: ");
+        for (ThickPolylines::const_iterator it = polylines->begin(); it != polylines->end(); ++ it) {
+            ThickLines lines = it->thicklines();
+            for (ThickLines::const_iterator it2 = lines.begin(); it2 != lines.end(); ++ it2) {
+                printf("%f,%f ", it2->a_width, it2->b_width);
+            }
+        }
+        printf("\n");
+    }
+    #endif /* SLIC3R_DEBUG */
+}
+
+void
+MedialAxis::process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline)
+{
+    while (true) {
+        // Since rot_next() works on the edge starting point but we want
+        // to find neighbors on the ending point, we just swap edge with
+        // its twin.
+        const VD::edge_type* twin = edge->twin();
+    
+        // count neighbors for this edge
+        std::vector<const VD::edge_type*> neighbors;
+        for (const VD::edge_type* neighbor = twin->rot_next(); neighbor != twin;
+            neighbor = neighbor->rot_next()) {
+            if (this->valid_edges.count(neighbor) > 0) neighbors.push_back(neighbor);
+        }
+    
+        // if we have a single neighbor then we can continue recursively
+        if (neighbors.size() == 1) {
+            const VD::edge_type* neighbor = neighbors.front();
+            
+            // break if this is a closed loop
+            if (this->edges.count(neighbor) == 0) return;
+            
+            Point new_point(neighbor->vertex1()->x(), neighbor->vertex1()->y());
+            polyline->points.push_back(new_point);
+            polyline->width.push_back(this->thickness[neighbor].second);
+            
+            (void)this->edges.erase(neighbor);
+            (void)this->edges.erase(neighbor->twin());
+            edge = neighbor;
+        } else if (neighbors.size() == 0) {
+            polyline->endpoints.second = true;
+            return;
+        } else {
+            // T-shaped or star-shaped joint
+            return;
+        }
+    }
+}
+
+bool
+MedialAxis::validate_edge(const VD::edge_type* edge)
+{
+    // prevent overflows and detect almost-infinite edges
+    if (std::abs(edge->vertex0()->x()) > double(MAX_COORD) || 
+        std::abs(edge->vertex0()->y()) > double(MAX_COORD) || 
+        std::abs(edge->vertex1()->x()) > double(MAX_COORD) ||
+        std::abs(edge->vertex1()->y()) > double(MAX_COORD))
+        return false;
+
+    // construct the line representing this edge of the Voronoi diagram
+    const Line line(
+        Point( edge->vertex0()->x(), edge->vertex0()->y() ),
+        Point( edge->vertex1()->x(), edge->vertex1()->y() )
+    );
+    
+    // discard edge if it lies outside the supplied shape
+    // this could maybe be optimized (checking inclusion of the endpoints
+    // might give false positives as they might belong to the contour itself)
+    if (line.a.coincides_with(line.b)) {
+        // in this case, contains(line) returns a false positive
+        if (!this->expolygon.contains(line.a)) return false;
+    } else {
+        //test if  (!expolygon.contains(line))
+        Polyline line_a_b;
+        line_a_b.append(line.a);
+        line_a_b.append(line.b);
+        Polylines external_bits = diff_pl(Polylines{ line_a_b }, expolygon);
+        if (!external_bits.empty()){
+            //check if the bits that are not inside are under epsilon length
+            coordf_t max_length = 0;
+            for (Polyline &poly : external_bits){
+                max_length = std::max(max_length, poly.length());
+            }
+            if (max_length > SCALED_EPSILON)
+                return false;
+        }
+    }
+    
+    // retrieve the original line segments which generated the edge we're checking
+    const VD::cell_type* cell_l = edge->cell();
+    const VD::cell_type* cell_r = edge->twin()->cell();
+    const Line &segment_l = this->retrieve_segment(cell_l);
+    const Line &segment_r = this->retrieve_segment(cell_r);
+    
+    
+    //SVG svg("edge.svg");
+    //svg.draw(this->expolygon.expolygon);
+    //svg.draw(line);
+    //svg.draw(segment_l, "red");
+    //svg.draw(segment_r, "blue");
+    //svg.Close();
+    //
+    
+    /*  Calculate thickness of the cross-section at both the endpoints of this edge.
+        Our Voronoi edge is part of a CCW sequence going around its Voronoi cell 
+        located on the left side. (segment_l).
+        This edge's twin goes around segment_r. Thus, segment_r is 
+        oriented in the same direction as our main edge, and segment_l is oriented
+        in the same direction as our twin edge.
+        We used to only consider the (half-)distances to segment_r, and that works
+        whenever segment_l and segment_r are almost specular and facing. However, 
+        at curves they are staggered and they only face for a very little length
+        (our very short edge represents such visibility).
+        Both w0 and w1 can be calculated either towards cell_l or cell_r with equal
+        results by Voronoi definition.
+        When cell_l or cell_r don't refer to the segment but only to an endpoint, we
+        calculate the distance to that endpoint instead.  */
+    
+    coordf_t w0 = cell_r->contains_segment()
+        ? line.a.distance_to(segment_r)*2
+        : line.a.distance_to(this->retrieve_endpoint(cell_r))*2;
+    
+    coordf_t w1 = cell_l->contains_segment()
+        ? line.b.distance_to(segment_l)*2
+        : line.b.distance_to(this->retrieve_endpoint(cell_l))*2;
+    
+    //don't remove the line that goes to the intersection of the contour
+    // we use them to create nicer thin wall lines
+    //if (cell_l->contains_segment() && cell_r->contains_segment()) {
+    //    // calculate the relative angle between the two boundary segments
+    //    double angle = fabs(segment_r.orientation() - segment_l.orientation());
+    //    if (angle > PI) angle = 2*PI - angle;
+    //    assert(angle >= 0 && angle <= PI);
+    //    
+    //    // fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
+    //    // we're interested only in segments close to the second case (facing segments)
+    //    // so we allow some tolerance.
+    //    // this filter ensures that we're dealing with a narrow/oriented area (longer than thick)
+    //    // we don't run it on edges not generated by two segments (thus generated by one segment
+    //    // and the endpoint of another segment), since their orientation would not be meaningful
+    //    if (PI - angle > PI/8) {
+    //        // angle is not narrow enough
+    //        
+    //        // only apply this filter to segments that are not too short otherwise their 
+    //        // angle could possibly be not meaningful
+    //        if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON || line.length() >= this->min_width)
+    //            return false;
+    //    }
+    //} else {
+    //    if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON)
+    //        return false;
+    //}
+
+    // don't do that before we try to fusion them
+    //if (w0 < this->min_width && w1 < this->min_width)
+    //    return false;
+    //
+
+    //shouldn't occur if perimeter_generator is well made. *1.05 for a little wiggle room
+    if (w0 > this->max_width*1.05 && w1 > this->max_width*1.05)
+        return false;
+    
+    this->thickness[edge]         = std::make_pair(w0, w1);
+    this->thickness[edge->twin()] = std::make_pair(w1, w0);
+    
+    return true;
+}
+
+const Line&
+MedialAxis::retrieve_segment(const VD::cell_type* cell) const
+{
+    return lines[cell->source_index()];
+}
+
+const Point&
+MedialAxis::retrieve_endpoint(const VD::cell_type* cell) const
+{
+    const Line& line = this->retrieve_segment(cell);
+    if (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) {
+        return line.a;
+    } else {
+        return line.b;
+    }
+}
+
+
+/// remove point that are at SCALED_EPSILON * 2 distance.
+void
+remove_point_too_near(ThickPolyline* to_reduce)
+{
+    const coord_t smallest = (coord_t)SCALED_EPSILON * 2;
+    size_t id = 1;
+    while (id < to_reduce->points.size() - 1) {
+        coord_t newdist = (coord_t)std::min(to_reduce->points[id].distance_to(to_reduce->points[id - 1])
+            , to_reduce->points[id].distance_to(to_reduce->points[id + 1]));
+        if (newdist < smallest) {
+            to_reduce->points.erase(to_reduce->points.begin() + id);
+            to_reduce->width.erase(to_reduce->width.begin() + id);
+            newdist = (coord_t)to_reduce->points[id].distance_to(to_reduce->points[id - 1]);
+            //if you removed a point, it check if the next one isn't too near from the previous one.
+            // if not, it bypass it.
+            if (newdist > smallest) {
+                ++id;
+            }
+        }
+        //go to next one
+        else ++id;
+    }
+}
+
+/// add points  from pattern to to_modify at the same % of the length
+/// so not add if an other point is present at the correct position
+void
+add_point_same_percent(ThickPolyline* pattern, ThickPolyline* to_modify)
+{
+    const double to_modify_length = to_modify->length();
+    const double percent_epsilon = SCALED_EPSILON / to_modify_length;
+    const double pattern_length = pattern->length();
+
+    double percent_length = 0;
+    for (size_t idx_point = 1; idx_point < pattern->points.size() - 1; ++idx_point) {
+        percent_length += pattern->points[idx_point-1].distance_to(pattern->points[idx_point]) / pattern_length;
+        //find position 
+        size_t idx_other = 1;
+        double percent_length_other_before = 0;
+        double percent_length_other = 0;
+        while (idx_other < to_modify->points.size()) {
+            percent_length_other_before = percent_length_other;
+            percent_length_other += to_modify->points[idx_other-1].distance_to(to_modify->points[idx_other])
+                / to_modify_length;
+            if (percent_length_other > percent_length - percent_epsilon) {
+                //if higher (we have gone over it)
+                break;
+            }
+            ++idx_other;
+        }
+        if (percent_length_other > percent_length + percent_epsilon) {
+            //insert a new point before the position
+            double percent_dist = (percent_length - percent_length_other_before) / (percent_length_other - percent_length_other_before);
+            coordf_t new_width = to_modify->width[idx_other - 1] * (1 - percent_dist);
+            new_width += to_modify->width[idx_other] * (percent_dist);
+            to_modify->width.insert(to_modify->width.begin() + idx_other, new_width);
+            to_modify->points.insert(
+                to_modify->points.begin() + idx_other,
+                to_modify->points[idx_other - 1].interpolate(percent_dist, to_modify->points[idx_other]));
+        }
+    }
+}
+
+/// find the nearest angle in the contour (or 2 nearest if it's difficult to choose) 
+/// return 1 for an angle of 90° and 0 for an angle of 0° or 180°
+/// find the nearest angle in the contour (or 2 nearest if it's difficult to choose) 
+/// return 1 for an angle of 90° and 0 for an angle of 0° or 180°
+double
+get_coeff_from_angle_countour(Point &point, const ExPolygon &contour, coord_t min_dist_between_point) {
+    double nearest_dist = point.distance_to(contour.contour.points.front());
+    Point point_nearest = contour.contour.points.front();
+    size_t id_nearest = 0;
+    double near_dist = nearest_dist;
+    Point point_near = point_nearest;
+    size_t id_near = 0;
+    for (size_t id_point = 1; id_point < contour.contour.points.size(); ++id_point) {
+        if (nearest_dist > point.distance_to(contour.contour.points[id_point])) {
+            //update point_near
+            id_near = id_nearest;
+            point_near = point_nearest;
+            near_dist = nearest_dist;
+            //update nearest
+            nearest_dist = point.distance_to(contour.contour.points[id_point]);
+            point_nearest = contour.contour.points[id_point];
+            id_nearest = id_point;
+        }
+    }
+    double angle = 0;
+    size_t id_before = id_nearest == 0 ? contour.contour.points.size() - 1 : id_nearest - 1;
+    Point point_before = id_nearest == 0 ? contour.contour.points.back() : contour.contour.points[id_nearest - 1];
+    //Search one point far enough to be relevant
+    while (point_nearest.distance_to(point_before) < min_dist_between_point) {
+        point_before = id_before == 0 ? contour.contour.points.back() : contour.contour.points[id_before - 1];
+        id_before = id_before == 0 ? contour.contour.points.size() - 1 : id_before - 1;
+        //don't loop
+        if (id_before == id_nearest) {
+            id_before = id_nearest == 0 ? contour.contour.points.size() - 1 : id_nearest - 1;
+            point_before = id_nearest == 0 ? contour.contour.points.back() : contour.contour.points[id_nearest - 1];
+            break;
+        }
+    }
+    size_t id_after = id_nearest == contour.contour.points.size() - 1 ? 0 : id_nearest + 1;
+    Point point_after = id_nearest == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_nearest + 1];
+    //Search one point far enough to be relevant
+    while (point_nearest.distance_to(point_after) < min_dist_between_point) {
+        point_after = id_after == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_after + 1];
+        id_after = id_after == contour.contour.points.size() - 1 ? 0 : id_after + 1;
+        //don't loop
+        if (id_after == id_nearest) {
+            id_after = id_nearest == contour.contour.points.size() - 1 ? 0 : id_nearest + 1;
+            point_after = id_nearest == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_nearest + 1];
+            break;
+        }
+    }
+    //compute angle
+    angle = point_nearest.ccw_angle(point_before, point_after);
+    if (angle >= PI) angle = 2 * PI - angle;  // smaller angle
+    //compute the diff from 90°
+    angle = abs(angle - PI / 2);
+    if (point_near.coincides_with(point_nearest) && std::max(nearest_dist, near_dist) + SCALED_EPSILON < point_nearest.distance_to(point_near)) {
+        //not only nearest
+        Point point_before = id_near == 0 ? contour.contour.points.back() : contour.contour.points[id_near - 1];
+        Point point_after = id_near == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_near + 1];
+        double angle2 = std::min(point_nearest.ccw_angle(point_before, point_after), point_nearest.ccw_angle(point_after, point_before));
+        angle2 = abs(angle - PI / 2);
+        angle = (angle + angle2) / 2;
+    }
+
+    return 1 - (angle / (PI / 2));
+}
+
+double
+dot(Line l1, Line l2)
+{
+    Vectorf v_1 = normalize(Vectorf(l1.b.x - l1.a.x, l1.b.y - l1.a.y));
+    Vectorf v_2 = normalize(Vectorf(l2.b.x - l2.a.x, l2.b.y - l2.a.y));
+    return v_1.x*v_2.x + v_1.y*v_2.y;
+}
+
+void
+MedialAxis::fusion_curve(ThickPolylines &pp)
+{
+    //fusion Y with only 1 '0' value => the "0" branch "pull" the cross-point
+    bool changes = false;
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+        // only consider 2-point polyline with endpoint
+        //if (polyline.points.size() != 2) continue; // too restrictive.
+        if (polyline.endpoints.first) polyline.reverse();
+        else if (!polyline.endpoints.second) continue;
+        if (polyline.width.back() > EPSILON) continue;
+
+        //check my length is small
+        coord_t length = (coord_t)polyline.length();
+        if (length > max_width) continue;
+
+        size_t closest_point_idx = this->expolygon.contour.closest_point_index(polyline.points.back());
+
+        //check the 0-width point is on the contour.
+        if (closest_point_idx == (size_t)-1) continue;
+
+        size_t prev_idx = closest_point_idx == 0 ? this->expolygon.contour.points.size() - 1 : closest_point_idx - 1;
+        size_t next_idx = closest_point_idx == this->expolygon.contour.points.size() - 1 ? 0 : closest_point_idx + 1;
+        double mindot = 1;
+        mindot = std::min(mindot, abs(dot(Line(polyline.points[polyline.points.size() - 1], polyline.points[polyline.points.size() - 2]),
+            (Line(this->expolygon.contour.points[closest_point_idx], this->expolygon.contour.points[prev_idx])))));
+        mindot = std::min(mindot, abs(dot(Line(polyline.points[polyline.points.size() - 1], polyline.points[polyline.points.size() - 2]),
+            (Line(this->expolygon.contour.points[closest_point_idx], this->expolygon.contour.points[next_idx])))));
+
+        //compute angle
+        double coeff_contour_angle = this->expolygon.contour.points[closest_point_idx].ccw_angle(this->expolygon.contour.points[prev_idx], this->expolygon.contour.points[next_idx]);
+        if (coeff_contour_angle >= PI) coeff_contour_angle = 2 * PI - coeff_contour_angle;  // smaller angle
+        //compute the diff from 90°
+        coeff_contour_angle = abs(coeff_contour_angle - PI / 2);
+
+
+        // look if other end is a cross point with almost 90° angle
+        double sum_dot = 0;
+        double min_dot = 0;
+        // look if other end is a cross point with multiple other branch
+        std::vector<size_t> crosspoint;
+        for (size_t j = 0; j < pp.size(); ++j) {
+            if (j == i) continue;
+            ThickPolyline& other = pp[j];
+            if (polyline.first_point().coincides_with(other.last_point())) {
+                other.reverse();
+                crosspoint.push_back(j);
+                double dot_temp = dot(Line(polyline.points[0], polyline.points[1]), (Line(other.points[0], other.points[1])));
+                min_dot = std::min(min_dot, abs(dot_temp));
+                sum_dot += dot_temp;
+            } else if (polyline.first_point().coincides_with(other.first_point())) {
+                crosspoint.push_back(j);
+                double dot_temp = dot(Line(polyline.points[0], polyline.points[1]), (Line(other.points[0], other.points[1])));
+                min_dot = std::min(min_dot, abs(dot_temp));
+                sum_dot += dot_temp;
+            }
+        }
+        sum_dot = abs(sum_dot);
+
+        //only consider very shallow angle for contour
+        if (mindot > 0.15 &&
+            (1 - (coeff_contour_angle / (PI / 2))) > 0.2) continue;
+
+        //check if it's a line that we can pull
+        if (crosspoint.size() != 2) continue;
+        if (sum_dot > 0.2) continue;
+        if (min_dot > 0.5) continue;
+        //don't remove useful bits. TODO: use the mindot to know by how much to multiply (1 when 90°, 1.42 when 45+, 1 when 0°)
+        if (polyline.length() > polyline.width.front()*1.42) continue;
+
+        //don't pull, it distords the line if there are too many points.
+        //// pull it a bit, depends on my size, the dot?, and the coeff at my 0-end (~14% for a square, almost 0 for a gentle curve)
+        //coord_t length_pull = polyline.length();
+        //length_pull *= 0.144 * get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, min(min_width, polyline.length() / 2));
+
+        ////compute dir
+        //Vectorf pull_direction(polyline.points[1].x - polyline.points[0].x, polyline.points[1].y - polyline.points[0].y);
+        //pull_direction = normalize(pull_direction);
+        //pull_direction.x *= length_pull;
+        //pull_direction.y *= length_pull;
+
+        ////pull the points
+        //Point &p1 = pp[crosspoint[0]].points[0];
+        //p1.x = p1.x + (coord_t)pull_direction.x;
+        //p1.y = p1.y + (coord_t)pull_direction.y;
+
+        //Point &p2 = pp[crosspoint[1]].points[0];
+        //p2.x = p2.x + (coord_t)pull_direction.x;
+        //p2.y = p2.y + (coord_t)pull_direction.y;
+
+        //delete the now unused polyline
+        pp.erase(pp.begin() + i);
+        --i;
+        changes = true;
+    }
+    if (changes) {
+        concatThickPolylines(pp);
+        ///reorder, in case of change
+        std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { return a.length() < b.length(); });
+        //have to redo it to remove multi-branch bits.
+        fusion_curve(pp);
+    }
+}
+
+void
+MedialAxis::fusion_corners(ThickPolylines &pp)
+{
+
+    //fusion Y with only 1 '0' value => the "0" branch "pull" the cross-point
+    bool changes = false;
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+        // only consider polyline with 0-end
+        if (polyline.points.size() != 2) continue;
+        if (polyline.endpoints.first) polyline.reverse();
+        else if (!polyline.endpoints.second) continue;
+        if (polyline.width.back() > 0) continue;
+
+        //check my length is small
+        coord_t length = (coord_t)polyline.length();
+        if (length > max_width) continue;
+
+        // look if other end is a cross point with multiple other branch
+        std::vector<size_t> crosspoint;
+        for (size_t j = 0; j < pp.size(); ++j) {
+            if (j == i) continue;
+            ThickPolyline& other = pp[j];
+            if (polyline.first_point().coincides_with(other.last_point())) {
+                other.reverse();
+                crosspoint.push_back(j);
+            } else if (polyline.first_point().coincides_with(other.first_point())) {
+                crosspoint.push_back(j);
+            }
+        }
+        //check if it's a line that we can pull
+        if (crosspoint.size() != 2) continue;
+
+        // check if i am at the external side of a curve
+        double angle1 = polyline.points[0].ccw_angle(polyline.points[1], pp[crosspoint[0]].points[1]);
+        if (angle1 >= PI) angle1 = 2 * PI - angle1;  // smaller angle
+        double angle2 = polyline.points[0].ccw_angle(polyline.points[1], pp[crosspoint[1]].points[1]);
+        if (angle2 >= PI) angle2 = 2 * PI - angle2;  // smaller angle
+        if (angle1 + angle2 < PI) continue;
+
+        //check if is smaller or the other ones are not endpoits
+        if (pp[crosspoint[0]].endpoints.second && length > pp[crosspoint[0]].length()) continue;
+        if (pp[crosspoint[1]].endpoints.second && length > pp[crosspoint[1]].length()) continue;
+
+        //FIXME: also pull (a bit less) points that are near to this one.
+        // if true, pull it a bit, depends on my size, the dot?, and the coeff at my 0-end (~14% for a square, almost 0 for a gentle curve)
+        coord_t length_pull = (coord_t)polyline.length();
+        length_pull *= (coord_t)( 0.144 * get_coeff_from_angle_countour(
+            polyline.points.back(), 
+            this->expolygon, 
+            std::min(min_width, (coord_t)(polyline.length() / 2))));
+
+        //compute dir
+        Vectorf pull_direction(polyline.points[1].x - polyline.points[0].x, polyline.points[1].y - polyline.points[0].y);
+        pull_direction = normalize(pull_direction);
+        pull_direction.x *= length_pull;
+        pull_direction.y *= length_pull;
+
+        //pull the points
+        Point &p1 = pp[crosspoint[0]].points[0];
+        p1.x = p1.x + pull_direction.x;
+        p1.y = p1.y + pull_direction.y;
+
+        Point &p2 = pp[crosspoint[1]].points[0];
+        p2.x = p2.x + pull_direction.x;
+        p2.y = p2.y + pull_direction.y;
+
+        //delete the now unused polyline
+        pp.erase(pp.begin() + i);
+        --i;
+        changes = true;
+    }
+    if (changes) {
+        concatThickPolylines(pp);
+        ///reorder, in case of change
+        std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { return a.length() < b.length(); });
+    }
+}
+
+void
+MedialAxis::extends_line_both_side(ThickPolylines& pp) {
+    const ExPolygons anchors = offset2_ex(to_polygons(diff_ex(*this->bounds, this->expolygon)), -SCALED_RESOLUTION, SCALED_RESOLUTION);
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+        this->extends_line(polyline, anchors, this->min_width);
+        if (!polyline.points.empty()) {
+            polyline.reverse();
+            this->extends_line(polyline, anchors, this->min_width);
+        }
+        if (polyline.points.empty()) {
+            pp.erase(pp.begin() + i);
+            --i;
+        }
+    }
+}
+
+void
+MedialAxis::extends_line(ThickPolyline& polyline, const ExPolygons& anchors, const coord_t join_width)
+{
+    // extend initial and final segments of each polyline if they're actual endpoints
+    // We assign new endpoints to temporary variables because in case of a single-line
+    // polyline, after we extend the start point it will be caught by the intersection()
+    // call, so we keep the inner point until we perform the second intersection() as well
+    if (polyline.endpoints.second && !bounds->has_boundary_point(polyline.points.back())) {
+        size_t first_idx = polyline.points.size() - 2;
+        Line line(*(polyline.points.begin() + first_idx), polyline.points.back());
+        while (line.length() < SCALED_RESOLUTION && first_idx>0) {
+            first_idx--;
+            line.a = *(polyline.points.begin() + first_idx);
+        }
+        // prevent the line from touching on the other side, otherwise intersection() might return that solution
+        if (polyline.points.size() == 2 && this->expolygon.contains(line.midpoint())) line.a = line.midpoint();
+
+        line.extend_end(max_width);
+        Point new_back;
+        if (this->expolygon.contour.has_boundary_point(polyline.points.back())) {
+            new_back = polyline.points.back();
+        } else {
+            bool finded = this->expolygon.contour.first_intersection(line, &new_back);
+            //verify also for holes.
+            Point new_back_temp;
+            for (Polygon hole : this->expolygon.holes) {
+                if (hole.first_intersection(line, &new_back_temp)) {
+                    if (!finded || line.a.distance_to(new_back_temp) < line.a.distance_to(new_back)) {
+                        finded = true;
+                        new_back = new_back_temp;
+                    }
+                }
+            }
+            // safety check if no intersection
+            if (!finded){
+                if (!this->expolygon.contains(line.b)) {
+                    //it's not possible to print that
+                    polyline.points.clear();
+                    polyline.width.clear();
+                    return;
+                }
+                new_back = line.b;
+            }
+
+            polyline.points.push_back(new_back);
+            polyline.width.push_back(polyline.width.back());
+        }
+        Point new_bound;
+        bool finded = bounds->contour.first_intersection(line, &new_bound);
+        //verify also for holes.
+        Point new_bound_temp;
+        for (Polygon hole : bounds->holes) {
+            if (hole.first_intersection(line, &new_bound_temp)) {
+                if (!finded || line.a.distance_to(new_bound_temp) < line.a.distance_to(new_bound)) {
+                    finded = true;
+                    new_bound = new_bound_temp;
+                }
+            }
+        }
+        // safety check if no intersection
+        if (!finded) {
+            if (line.b.coincides_with_epsilon(polyline.points.back())) return;
+            //check if we don't over-shoot inside us
+            bool is_in_anchor = false;
+            for (const ExPolygon& a : anchors) {
+                if (a.contains(line.b)) {
+                    is_in_anchor = true;
+                    break;
+                }
+            }
+            if (!is_in_anchor) return;
+            new_bound = line.b;
+        }
+        // find anchor
+        Point best_anchor;
+        double shortest_dist = max_width;
+        for (const ExPolygon& a : anchors) {
+            Point p_maybe_inside = a.contour.centroid();
+            double test_dist = new_bound.distance_to(p_maybe_inside) + new_back.distance_to(p_maybe_inside);
+            //if (test_dist < max_width / 2 && (test_dist < shortest_dist || shortest_dist < 0)) {
+            double angle_test = new_back.ccw_angle(p_maybe_inside, line.a);
+            if (angle_test > PI) angle_test = 2 * PI - angle_test;
+            if (test_dist < max_width && test_dist<shortest_dist && abs(angle_test) > PI / 2) {
+                shortest_dist = test_dist;
+                best_anchor = p_maybe_inside;
+            }
+        }
+        if (best_anchor.x != 0 && best_anchor.y != 0) {
+            Point p_obj = best_anchor + new_bound;
+            p_obj.x /= 2;
+            p_obj.y /= 2;
+            Line l2 = Line(new_back, p_obj);
+            l2.extend_end(max_width);
+            (void)bounds->contour.first_intersection(l2, &new_bound);
+        }
+        if (new_bound.coincides_with_epsilon(new_back)) {
+            return;
+        }
+        polyline.points.push_back(new_bound);
+        //polyline.width.push_back(join_width);
+        //it thickens the line a bit too early, imo
+        polyline.width.push_back(polyline.width.back());
+    }
+}
+
+void
+MedialAxis::main_fusion(ThickPolylines& pp)
+{
+
+    bool changes = true;
+    std::map<Point, double> coeff_angle_cache;
+    while (changes) {
+        concatThickPolylines(pp);
+        //reoder pp by length (ascending) It's really important to do that to avoid building the line from the width insteand of the length
+        std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) {
+            bool ahas0 = a.width.front() == 0 || a.width.back() == 0;
+            bool bhas0 = b.width.front() == 0 || b.width.back() == 0;
+            if (ahas0 && !bhas0) return true;
+            if (!ahas0 && bhas0) return false;
+            return a.length() < b.length();
+        });
+        changes = false;
+        for (size_t i = 0; i < pp.size(); ++i) {
+            ThickPolyline& polyline = pp[i];
+
+            //simple check to see if i can be fusionned
+            if (!polyline.endpoints.first && !polyline.endpoints.second) continue;
+
+
+            ThickPolyline* best_candidate = nullptr;
+            float best_dot = -1;
+            size_t best_idx = 0;
+            double dot_poly_branch = 0;
+            double dot_candidate_branch = 0;
+
+            bool find_main_branch = false;
+            size_t biggest_main_branch_id = 0;
+            coord_t biggest_main_branch_length = 0;
+
+            // find another polyline starting here
+            for (size_t j = i + 1; j < pp.size(); ++j) {
+                ThickPolyline& other = pp[j];
+                if (polyline.last_point().coincides_with(other.last_point())) {
+                    polyline.reverse();
+                    other.reverse();
+                } else if (polyline.first_point().coincides_with(other.last_point())) {
+                    other.reverse();
+                } else if (polyline.first_point().coincides_with(other.first_point())) {
+                } else if (polyline.last_point().coincides_with(other.first_point())) {
+                    polyline.reverse();
+                } else {
+                    continue;
+                }
+                
+                //// mergeable tests
+                if (polyline.points.size() < 2 && other.points.size() < 2) continue;
+                if (!polyline.endpoints.second || !other.endpoints.second) continue;
+                // test if the new width will not be too big if a fusion occur
+                //note that this isn't the real calcul. It's just to avoid merging lines too far apart.
+                if (
+                    ((polyline.points.back().distance_to(other.points.back())
+                    + (polyline.width.back() + other.width.back()) / 4)
+                > max_width*1.05))
+                    continue;
+                // test if the lines are not too different in length.
+                if (abs(polyline.length() - other.length()) > max_width) continue;
+
+
+                //test if we don't merge with something too different and without any relevance.
+                double coeffSizePolyI = 1;
+                if (polyline.width.back() == 0) {
+                    coeffSizePolyI = 0.1 + 0.9*get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, std::min(min_width, (coord_t)(polyline.length() / 2)));
+                }
+                double coeffSizeOtherJ = 1;
+                if (other.width.back() == 0) {
+                    coeffSizeOtherJ = 0.1 + 0.9*get_coeff_from_angle_countour(other.points.back(), this->expolygon, std::min(min_width, (coord_t)(polyline.length() / 2)));
+                }
+                if (abs(polyline.length()*coeffSizePolyI - other.length()*coeffSizeOtherJ) > max_width / 2) continue;
+
+
+                //compute angle to see if it's better than previous ones (straighter = better).
+                //we need to add how strait we are from our main.
+                float test_dot = (float)(dot(polyline.lines().front(), other.lines().front()));
+
+                // Get the branch/line in wich we may merge, if possible
+                // with that, we can decide what is important, and how we can merge that.
+                // angle_poly - angle_candi =90� => one is useless
+                // both angle are equal => both are useful with same strength
+                // ex: Y => | both are useful to crete a nice line
+                // ex2: TTTTT => -----  these 90� useless lines should be discarded
+                find_main_branch = false;
+                biggest_main_branch_id = 0;
+                biggest_main_branch_length = 0;
+                for (size_t k = 0; k < pp.size(); ++k) {
+                    if (k == i || k == j) continue;
+                    ThickPolyline& main = pp[k];
+                    if (polyline.first_point().coincides_with(main.last_point())) {
+                        main.reverse();
+                        if (!main.endpoints.second)
+                            find_main_branch = true;
+                        else if (biggest_main_branch_length < main.length()) {
+                            biggest_main_branch_id = k;
+                            biggest_main_branch_length = (coord_t)main.length();
+                        }
+                    } else if (polyline.first_point().coincides_with(main.first_point())) {
+                        if (!main.endpoints.second)
+                            find_main_branch = true;
+                        else if (biggest_main_branch_length < main.length()) {
+                            biggest_main_branch_id = k;
+                            biggest_main_branch_length = (coord_t)main.length();
+                        }
+                    }
+                    if (find_main_branch) {
+                        //use this variable to store the good index and break to compute it
+                        biggest_main_branch_id = k;
+                        break;
+                    }
+                }
+                double dot_poly_branch_test = 0.707;
+                double dot_candidate_branch_test = 0.707;
+                if (!find_main_branch && biggest_main_branch_length == 0) {
+                    // nothing -> it's impossible!
+                    dot_poly_branch_test = 0.707;
+                    dot_candidate_branch_test = 0.707;
+                } else if (!find_main_branch && (
+                    (pp[biggest_main_branch_id].length() < polyline.length() && (polyline.width.back() != 0 || pp[biggest_main_branch_id].width.back() ==0)) 
+                    || (pp[biggest_main_branch_id].length() < other.length() && (other.width.back() != 0 || pp[biggest_main_branch_id].width.back() == 0)))) {
+                    //the main branch should have no endpoint or be bigger!
+                    //here, it have an endpoint, and is not the biggest -> bad!
+                    continue;
+                } else {
+                    //compute the dot (biggest_main_branch_id)
+                    dot_poly_branch_test = -dot(Line(polyline.points[0], polyline.points[1]), Line(pp[biggest_main_branch_id].points[0], pp[biggest_main_branch_id].points[1]));
+                    dot_candidate_branch_test = -dot(Line(other.points[0], other.points[1]), Line(pp[biggest_main_branch_id].points[0], pp[biggest_main_branch_id].points[1]));
+                    if (dot_poly_branch_test < 0) dot_poly_branch_test = 0;
+                    if (dot_candidate_branch_test < 0) dot_candidate_branch_test = 0;
+                    if (pp[biggest_main_branch_id].width.back()>0)
+                        test_dot += 2 * (float)dot_poly_branch;
+                }
+                //test if it's useful to merge or not
+                //ie, don't merge  'T' but ok for 'Y', merge only lines of not disproportionate different length (ratio max: 4) (or they are both with 0-width end)
+                if (dot_poly_branch_test < 0.1 || dot_candidate_branch_test < 0.1 ||
+                    (
+                        ((polyline.length()>other.length() ? polyline.length() / other.length() : other.length() / polyline.length()) > 4) 
+                        && !(polyline.width.back() == 0 && other.width.back()==0)
+                        )) {
+                    continue;
+                }
+                if (test_dot > best_dot) {
+                    best_candidate = &other;
+                    best_idx = j;
+                    best_dot = test_dot;
+                    dot_poly_branch = dot_poly_branch_test;
+                    dot_candidate_branch = dot_candidate_branch_test;
+                }
+            }
+            if (best_candidate != nullptr) {
+                // delete very near points
+                remove_point_too_near(&polyline);
+                remove_point_too_near(best_candidate);
+
+                // add point at the same pos than the other line to have a nicer fusion
+                add_point_same_percent(&polyline, best_candidate);
+                add_point_same_percent(best_candidate, &polyline);
+
+                //get the angle of the nearest points of the contour to see : _| (good) \_ (average) __(bad)
+                //sqrt because the result are nicer this way: don't over-penalize /_ angles
+                //TODO: try if we can achieve a better result if we use a different algo if the angle is <90�
+                const double coeff_angle_poly = (coeff_angle_cache.find(polyline.points.back()) != coeff_angle_cache.end())
+                    ? coeff_angle_cache[polyline.points.back()]
+                    : (get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, std::min(min_width, (coord_t)(polyline.length() / 2))));
+                const double coeff_angle_candi = (coeff_angle_cache.find(best_candidate->points.back()) != coeff_angle_cache.end())
+                    ? coeff_angle_cache[best_candidate->points.back()]
+                    : (get_coeff_from_angle_countour(best_candidate->points.back(), this->expolygon, std::min(min_width, (coord_t)(best_candidate->length() / 2))));
+
+                //this will encourage to follow the curve, a little, because it's shorter near the center
+                //without that, it tends to go to the outter rim.
+                double weight_poly = 2 - (polyline.length() / std::max(polyline.length(), best_candidate->length()));
+                double weight_candi = 2 - (best_candidate->length() / std::max(polyline.length(), best_candidate->length()));
+                weight_poly *= coeff_angle_poly;
+                weight_candi *= coeff_angle_candi;
+                const double coeff_poly = (dot_poly_branch * weight_poly) / (dot_poly_branch * weight_poly + dot_candidate_branch * weight_candi);
+                const double coeff_candi = 1.0 - coeff_poly;
+                //iterate the points
+                // as voronoi should create symetric thing, we can iterate synchonously
+                size_t idx_point = 1;
+                while (idx_point < std::min(polyline.points.size(), best_candidate->points.size())) {
+                    //fusion
+                    polyline.points[idx_point].x = polyline.points[idx_point].x * coeff_poly + best_candidate->points[idx_point].x * coeff_candi;
+                    polyline.points[idx_point].y = polyline.points[idx_point].y * coeff_poly + best_candidate->points[idx_point].y * coeff_candi;
+
+                    // The width decrease with distance from the centerline.
+                    // This formula is what works the best, even if it's not perfect (created empirically).  0->3% error on a gap fill on some tests.
+                    //If someone find  an other formula based on the properties of the voronoi algorithm used here, and it works better, please use it.
+                    //or maybe just use the distance to nearest edge in bounds->..
+                    double value_from_current_width = 0.5*polyline.width[idx_point] * dot_poly_branch / std::max(dot_poly_branch, dot_candidate_branch);
+                    value_from_current_width += 0.5*best_candidate->width[idx_point] * dot_candidate_branch / std::max(dot_poly_branch, dot_candidate_branch);
+                    double value_from_dist = 2 * polyline.points[idx_point].distance_to(best_candidate->points[idx_point]);
+                    value_from_dist *= sqrt(std::min(dot_poly_branch, dot_candidate_branch) / std::max(dot_poly_branch, dot_candidate_branch));
+                    polyline.width[idx_point] = value_from_current_width + value_from_dist;
+                    //failsafes
+                    if (polyline.width[idx_point] > max_width)
+                        polyline.width[idx_point] = max_width;
+                    //failsafe: try to not go out of the radius of the section, take the width of the merging point for that. (and with some offset)
+                    coord_t main_branch_width = pp[biggest_main_branch_id].width.front();
+                    coord_t main_branch_dist = pp[biggest_main_branch_id].points.front().distance_to(polyline.points[idx_point]);
+                    coord_t max_width_from_main = std::sqrt(main_branch_width*main_branch_width + main_branch_dist*main_branch_dist);
+                    if (find_main_branch && polyline.width[idx_point] > max_width_from_main)
+                        polyline.width[idx_point] = max_width_from_main;
+                    if (find_main_branch && polyline.width[idx_point] > pp[biggest_main_branch_id].width.front() * 1.1)
+                        polyline.width[idx_point] = pp[biggest_main_branch_id].width.front() * 1.1;
+
+                    ++idx_point;
+                }
+                if (idx_point < best_candidate->points.size()) {
+                    if (idx_point + 1 < best_candidate->points.size()) {
+                        //create a new polyline
+                        pp.emplace_back();
+                        pp.back().endpoints.first = true;
+                        pp.back().endpoints.second = best_candidate->endpoints.second;
+                        for (size_t idx_point_new_line = idx_point; idx_point_new_line < best_candidate->points.size(); ++idx_point_new_line) {
+                            pp.back().points.push_back(best_candidate->points[idx_point_new_line]);
+                            pp.back().width.push_back(best_candidate->width[idx_point_new_line]);
+                        }
+                    } else {
+                        //Add last point
+                        polyline.points.push_back(best_candidate->points[idx_point]);
+                        polyline.width.push_back(best_candidate->width[idx_point]);
+                        //select if an end occur
+                        polyline.endpoints.second &= best_candidate->endpoints.second;
+                    }
+
+                } else {
+                    //select if an end occur
+                    polyline.endpoints.second &= best_candidate->endpoints.second;
+                }
+
+                //remove points that are the same or too close each other, ie simplify
+                for (size_t idx_point = 1; idx_point < polyline.points.size(); ++idx_point) {
+                    if (polyline.points[idx_point - 1].distance_to(polyline.points[idx_point]) < SCALED_EPSILON) {
+                        if (idx_point < polyline.points.size() - 1) {
+                            polyline.points.erase(polyline.points.begin() + idx_point);
+                            polyline.width.erase(polyline.width.begin() + idx_point);
+                        } else {
+                            polyline.points.erase(polyline.points.begin() + idx_point - 1);
+                            polyline.width.erase(polyline.width.begin() + idx_point - 1);
+                        }
+                        --idx_point;
+                    }
+                }
+                //remove points that are outside of the geometry
+                for (size_t idx_point = 0; idx_point < polyline.points.size(); ++idx_point) {
+                    if (!bounds->contains_b(polyline.points[idx_point])) {
+                        polyline.points.erase(polyline.points.begin() + idx_point);
+                        polyline.width.erase(polyline.width.begin() + idx_point);
+                        --idx_point;
+                    }
+                }
+
+                if (polyline.points.size() < 2) {
+                    //remove self
+                    pp.erase(pp.begin() + i);
+                    --i;
+                    --best_idx;
+                } else {
+                    //update cache
+                    coeff_angle_cache[polyline.points.back()] = coeff_angle_poly * coeff_poly + coeff_angle_candi * coeff_candi;
+                }
+
+                pp.erase(pp.begin() + best_idx);
+                changes = true;
+                break;
+            }
+        }
+    }
+}
+
+void
+MedialAxis::remove_too_thin_extrusion(ThickPolylines& pp)
+{
+    // remove too thin extrusion at start & end of polylines
+    bool changes = false;
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+        // remove bits with too small extrusion
+        while (polyline.points.size() > 1 && polyline.width.front() < this->min_width && polyline.endpoints.first) {
+            //try to split if possible
+            if (polyline.width[1] > min_width) {
+                double percent_can_keep = (min_width - polyline.width[0]) / (polyline.width[1] - polyline.width[0]);
+                if (polyline.points.front().distance_to(polyline.points[1]) * (1 - percent_can_keep) > SCALED_RESOLUTION) {
+                    //Can split => move the first point and assign a new weight.
+                    //the update of endpoints wil be performed in concatThickPolylines
+                    polyline.points.front() = polyline.points.front().interpolate(percent_can_keep, polyline.points[1]);
+                    polyline.width.front() = min_width;
+                } else {
+                    /// almost 0-length, Remove
+                    polyline.points.erase(polyline.points.begin());
+                    polyline.width.erase(polyline.width.begin());
+                }
+                changes = true;
+                break;
+            }
+            polyline.points.erase(polyline.points.begin());
+            polyline.width.erase(polyline.width.begin());
+            changes = true;
+        }
+        while (polyline.points.size() > 1 && polyline.width.back() < this->min_width && polyline.endpoints.second) {
+            //try to split if possible
+            if (polyline.width[polyline.points.size() - 2] > min_width) {
+                double percent_can_keep = (min_width - polyline.width.back()) / (polyline.width[polyline.points.size() - 2] - polyline.width.back());
+                if (polyline.points.back().distance_to(polyline.points[polyline.points.size() - 2]) * (1 - percent_can_keep) > SCALED_RESOLUTION) {
+                    //Can split => move the first point and assign a new weight.
+                    //the update of endpoints wil be performed in concatThickPolylines
+                    polyline.points.back() = polyline.points.back().interpolate(percent_can_keep, polyline.points[polyline.points.size() - 2]);
+                    polyline.width.back() = min_width;
+                } else {
+                    /// almost 0-length, Remove
+                    polyline.points.erase(polyline.points.end() - 1);
+                    polyline.width.erase(polyline.width.end() - 1);
+                }
+                changes = true;
+                break;
+            }
+            polyline.points.erase(polyline.points.end() - 1);
+            polyline.width.erase(polyline.width.end() - 1);
+            changes = true;
+        }
+        //remove points and bits that comes from a "main line"
+        if (polyline.points.size() < 2 || (changes && polyline.length() < max_width && polyline.points.size() ==2)) {
+            //remove self if too small
+            pp.erase(pp.begin() + i);
+            --i;
+        }
+    }
+    if (changes) concatThickPolylines(pp);
+}
+
+void
+MedialAxis::concatenate_polylines_with_crossing(ThickPolylines& pp)
+{
+
+    // concatenate, but even where multiple thickpolyline join, to create nice long strait polylines
+    /*  If we removed any short polylines we now try to connect consecutive polylines
+    in order to allow loop detection. Note that this algorithm is greedier than
+    MedialAxis::process_edge_neighbors() as it will connect random pairs of
+    polylines even when more than two start from the same point. This has no
+    drawbacks since we optimize later using nearest-neighbor which would do the
+    same, but should we use a more sophisticated optimization algorithm we should
+    not connect polylines when more than two meet.
+    Optimisation of the old algorithm : now we select the most "strait line" choice
+    when we merge with an other line at a point with more than two meet.
+    */
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+        if (polyline.endpoints.first && polyline.endpoints.second) continue; // optimization
+
+        ThickPolyline* best_candidate = nullptr;
+        float best_dot = -1;
+        size_t best_idx = 0;
+
+        // find another polyline starting here
+        for (size_t j = 0; j < pp.size(); ++j) {
+            if (j == i) continue;
+            ThickPolyline& other = pp[j];
+            if (other.endpoints.first && other.endpoints.second) continue;
+            bool me_reverse = false;
+            bool other_reverse = false;
+            if (polyline.last_point().coincides_with(other.last_point())) {
+                other_reverse = true;
+            } else if (polyline.first_point().coincides_with(other.last_point())) {
+                me_reverse = true;
+                other_reverse = true;
+            } else if (polyline.first_point().coincides_with(other.first_point())) {
+                me_reverse = true;
+            } else if (!polyline.last_point().coincides_with(other.first_point())) {
+                continue;
+            }
+
+            Pointf v_poly(me_reverse ? polyline.lines().front().vector().x : polyline.lines().back().vector().x, 
+                me_reverse ? polyline.lines().front().vector().y : polyline.lines().back().vector().y);
+            v_poly.scale(1 / std::sqrt(v_poly.x*v_poly.x + v_poly.y*v_poly.y));
+            Pointf v_other(other_reverse ? other.lines().back().vector().x : other.lines().front().vector().x,
+                other_reverse ? other.lines().back().vector().y : other.lines().front().vector().y);
+            v_other.scale(1 / std::sqrt(v_other.x*v_other.x + v_other.y*v_other.y));
+            float other_dot = std::abs(float(v_poly.x*v_other.x + v_poly.y*v_other.y));
+            if (other_dot > best_dot) {
+                best_candidate = &other;
+                best_idx = j;
+                best_dot = other_dot;
+            }
+        }
+        if (best_candidate != nullptr && best_candidate->points.size() > 1) {
+            if (polyline.last_point().coincides_with(best_candidate->last_point())) {
+                best_candidate->reverse();
+            } else if (polyline.first_point().coincides_with(best_candidate->last_point())) {
+                polyline.reverse();
+                best_candidate->reverse();
+            } else if (polyline.first_point().coincides_with(best_candidate->first_point())) {
+                polyline.reverse();
+            }
+            //intersections may create over-extrusion because the included circle can be a bit larger. We have to make it short again if needed.
+            if (polyline.points.size() > 1 && best_candidate->points.size() > 1 
+                    && polyline.width.back() > polyline.width[polyline.width.size() - 2] 
+                    && polyline.width.back() > best_candidate->width[1]) {
+                polyline.width.back() = std::min(polyline.width[polyline.width.size() - 2], best_candidate->width[1]);
+            }
+            polyline.points.insert(polyline.points.end(), best_candidate->points.begin() + 1, best_candidate->points.end());
+            polyline.width.insert(polyline.width.end(), best_candidate->width.begin() + 1, best_candidate->width.end());
+            polyline.endpoints.second = best_candidate->endpoints.second;
+            assert(polyline.width.size() == polyline.points.size());
+            if (best_idx < i) i--;
+            pp.erase(pp.begin() + best_idx);
+        }
+    }
+}
+
+void
+MedialAxis::remove_too_thin_points(ThickPolylines& pp)
+{
+    //remove too thin polylines points (inside a polyline : split it)
+    for (size_t i = 0; i < pp.size(); ++i) {
+        ThickPolyline& polyline = pp[i];
+
+        // remove bits with too small extrusion
+        size_t idx_point = 0;
+        while (idx_point<polyline.points.size()) {
+            if (polyline.width[idx_point] < min_width) {
+                if (idx_point == 0) {
+                    //too thin at start
+                    polyline.points.erase(polyline.points.begin());
+                    polyline.width.erase(polyline.width.begin());
+                    idx_point = 0;
+                } else if (idx_point == 1) {
+                    //too thin at start
+                    polyline.points.erase(polyline.points.begin());
+                    polyline.width.erase(polyline.width.begin());
+                    polyline.points.erase(polyline.points.begin());
+                    polyline.width.erase(polyline.width.begin());
+                    idx_point = 0;
+                } else if (idx_point == polyline.points.size() - 2) {
+                    //too thin at (near) end
+                    polyline.points.erase(polyline.points.end() - 1);
+                    polyline.width.erase(polyline.width.end() - 1);
+                    polyline.points.erase(polyline.points.end() - 1);
+                    polyline.width.erase(polyline.width.end() - 1);
+                } else if (idx_point == polyline.points.size() - 1) {
+                    //too thin at end
+                    polyline.points.erase(polyline.points.end() - 1);
+                    polyline.width.erase(polyline.width.end() - 1);
+                } else {
+                    //too thin in middle : split
+                    pp.emplace_back();
+                    ThickPolyline &newone = pp.back();
+                    newone.points.insert(newone.points.begin(), polyline.points.begin() + idx_point + 1, polyline.points.end());
+                    newone.width.insert(newone.width.begin(), polyline.width.begin() + idx_point + 1, polyline.width.end());
+                    polyline.points.erase(polyline.points.begin() + idx_point, polyline.points.end());
+                    polyline.width.erase(polyline.width.begin() + idx_point, polyline.width.end());
+                }
+            } else idx_point++;
+
+            if (polyline.points.size() < 2) {
+                //remove self if too small
+                pp.erase(pp.begin() + i);
+                --i;
+                break;
+            }
+        }
+    }
+}
+
+void
+MedialAxis::remove_too_short_polylines(ThickPolylines& pp, const coord_t min_size)
+{
+    //remove too short polyline
+    bool changes = true;
+    while (changes) {
+        changes = false;
+
+        double shortest_size = min_size;
+        size_t shortest_idx = -1;
+        for (size_t i = 0; i < pp.size(); ++i) {
+            ThickPolyline& polyline = pp[i];
+            // Remove the shortest polylines : polyline that are shorter than wider
+            // (we can't do this check before endpoints extension and clipping because we don't
+            // know how long will the endpoints be extended since it depends on polygon thickness
+            // which is variable - extension will be <= max_width/2 on each side) 
+            if ((polyline.endpoints.first || polyline.endpoints.second)
+                && polyline.length() < max_width / 2) {
+                if (shortest_size > polyline.length()) {
+                    shortest_size = polyline.length();
+                    shortest_idx = i;
+                }
+
+            }
+        }
+        if (shortest_idx < pp.size()) {
+            pp.erase(pp.begin() + shortest_idx);
+            changes = true;
+        }
+        if (changes) concatThickPolylines(pp);
+    }
+}
+
+void
+MedialAxis::ensure_not_overextrude(ThickPolylines& pp)
+{
+    //ensure the volume extruded is correct for what we have been asked
+    // => don't over-extrude
+    double surface = 0;
+    double volume = 0;
+    for (ThickPolyline& polyline : pp) {
+        for (ThickLine &l : polyline.thicklines()) {
+            surface += l.length() * (l.a_width + l.b_width) / 2;
+            double width_mean = (l.a_width + l.b_width) / 2;
+            volume += height * (width_mean - height * (1. - 0.25 * PI)) * l.length();
+        }
+    }
+
+    // compute bounds volume
+    double boundsVolume = 0;
+    boundsVolume += height*bounds->area();
+    // add external "perimeter gap"
+    double perimeterRoundGap = bounds->contour.length() * height * (1 - 0.25*PI) * 0.5;
+    // add holes "perimeter gaps"
+    double holesGaps = 0;
+    for (const Polygon &hole : bounds->holes) {
+        holesGaps += hole.length() * height * (1 - 0.25*PI) * 0.5;
+    }
+    boundsVolume += perimeterRoundGap + holesGaps;
+
+    if (boundsVolume < volume) {
+        //reduce width
+        double reduce_by = boundsVolume / volume;
+        for (ThickPolyline& polyline : pp) {
+            for (coordf_t &width : polyline.width) {
+                width *= reduce_by;
+            }
+        }
+    }
+}
+
+void 
+MedialAxis::simplify_polygon_frontier()
+{
+
+    //simplify the boundary between us and the bounds->
+    //it will remove every point in the surface contour that aren't on the bounds contour
+    this->expolygon = this->surface;
+    this->expolygon.contour.remove_collinear_points();
+    for (Polygon &hole : this->expolygon.holes)
+        hole.remove_collinear_points();
+    if (&this->surface != this->bounds) {
+        bool need_intersect = false;
+        for (size_t i = 0; i < this->expolygon.contour.points.size(); i++) {
+            Point &p_check = this->expolygon.contour.points[i];
+            //if (!find) {
+            if (!bounds->has_boundary_point(p_check)) {
+                //check if we put it at a bound point instead of delete it
+                size_t prev_i = i == 0 ? this->expolygon.contour.points.size() - 1 : (i - 1);
+                size_t next_i = i == this->expolygon.contour.points.size() - 1 ? 0 : (i + 1);
+                const Point* closest = bounds->contour.closest_point(p_check);
+                if (closest != nullptr && closest->distance_to(p_check) + SCALED_EPSILON
+                    < std::min(p_check.distance_to(this->expolygon.contour.points[prev_i]), p_check.distance_to(this->expolygon.contour.points[next_i])) / 2) {
+                    p_check.x = closest->x;
+                    p_check.y = closest->y;
+                    need_intersect = true;
+                } else {
+                    this->expolygon.contour.points.erase(this->expolygon.contour.points.begin() + i);
+                    i--;
+                }
+            }
+        }
+        if (need_intersect) {
+            ExPolygons simplified_polygons = intersection_ex(this->expolygon, *bounds);
+            if (simplified_polygons.size() == 1) {
+                this->expolygon = simplified_polygons[0];
+            } else {
+                //can't simply that much, reuse the given one
+                this->expolygon = this->surface;
+                this->expolygon.contour.remove_collinear_points();
+                for (Polygon &hole : this->expolygon.holes)
+                    hole.remove_collinear_points();
+            }
+        }
+    }
+
+    if (!this->expolygon.contour.points.empty())
+        this->expolygon.remove_point_too_near((coord_t)SCALED_RESOLUTION);
+}
+
+/// Grow the extrusion to at least nozzle_diameter*1.05 (lowest safe extrusion width)
+/// Do not grow points inside the anchor.
+void
+MedialAxis::grow_to_nozzle_diameter(ThickPolylines& pp, const ExPolygons& anchors)
+{
+    //compute the min width
+    coord_t min_width = this->nozzle_diameter;
+    if (this->height > 0) min_width = Flow::new_from_spacing(
+        float(unscale(this->nozzle_diameter)),
+        float(unscale(this->nozzle_diameter)),
+        float(unscale(this->height)), false).scaled_width();
+    //ensure the width is not lower than min_width.
+    for (ThickPolyline& polyline : pp) {
+        for (int i = 0; i < polyline.points.size(); ++i) {
+            bool is_anchored = false;
+            for (const ExPolygon &poly : anchors) {
+                if (poly.contains(polyline.points[i])) {
+                    is_anchored = true;
+                    break;
+                }
+            }
+            if (!is_anchored && polyline.width[i] < min_width)
+                polyline.width[i] = min_width;
+        }
+    }
+}
+
+void
+MedialAxis::taper_ends(ThickPolylines& pp)
+{
+    // minimum size of the taper: be sure to extrude at least the "round edges" of the extrusion (0-spacing extrusion).
+    const coord_t min_size = std::max(this->nozzle_diameter * 0.1, this->height * (1. - 0.25 * PI));
+    const coordf_t length = std::min(this->taper_size, (this->nozzle_diameter - min_size) / 2);
+    if (length <= SCALED_RESOLUTION) return;
+    //ensure the width is not lower than min_size.
+    for (ThickPolyline& polyline : pp) {
+        if (polyline.length() < length * 2.2) continue;
+        if (polyline.endpoints.first) {
+            polyline.width[0] = min_size;
+            coord_t current_dist = min_size;
+            coord_t last_dist = min_size;
+            for (size_t i = 1; i<polyline.width.size(); ++i) {
+                current_dist += (coord_t) polyline.points[i - 1].distance_to(polyline.points[i]);
+                if (current_dist > length) {
+                    //create a new point if not near enough
+                    if (current_dist > length + SCALED_RESOLUTION) {
+                        coordf_t percent_dist = (length - last_dist) / (current_dist - last_dist);
+                        polyline.points.insert(polyline.points.begin() + i, polyline.points[i - 1].interpolate(percent_dist, polyline.points[i]));
+                        polyline.width.insert(polyline.width.begin() + i, polyline.width[i]);
+                    }
+                    break;
+                }
+                polyline.width[i] = std::max((coordf_t)min_size, min_size + (polyline.width[i] - min_size) * current_dist / length);
+                last_dist = current_dist;
+            }
+        }
+        if (polyline.endpoints.second) {
+            polyline.width[polyline.width.size() - 1] = min_size;
+            coord_t current_dist = min_size;
+            coord_t last_dist = min_size;
+            for (size_t i = polyline.width.size()-1; i > 0; --i) {
+                current_dist += (coord_t)polyline.points[i].distance_to(polyline.points[i - 1]);
+                if (current_dist > length) {
+                    //create new point if not near enough
+                    if (current_dist > length + SCALED_RESOLUTION) {
+                        coordf_t percent_dist = (length - last_dist) / (current_dist - last_dist);
+                        polyline.points.insert(polyline.points.begin() + i, polyline.points[i].interpolate(percent_dist, polyline.points[i - 1]));
+                        polyline.width.insert(polyline.width.begin() + i, polyline.width[i - 1]);
+                    }
+                    break;
+                }
+                polyline.width[i - 1] = std::max((coordf_t)min_size, min_size + (polyline.width[i - 1] - min_size) * current_dist / length);
+                last_dist = current_dist;
+            }
+        }
+    }
+}
+
+void
+MedialAxis::build(ThickPolylines &polylines_out)
+{
+
+    simplify_polygon_frontier();
+    //safety check
+    if (this->expolygon.area() < this->min_width * this->min_width) this->expolygon = this->surface;
+    if (this->expolygon.area() < this->min_width * this->min_width) return;
+
+
+    // compute the Voronoi diagram and extract medial axis polylines
+    ThickPolylines pp;
+    this->polyline_from_voronoi(this->expolygon.lines(), &pp);
+
+    //sanity check, as the voronoi can return (abeit very rarely) randomly high values.
+    for (ThickPolyline &tp : pp) {
+        for (int i = 0; i < tp.width.size(); i++) {
+            if (tp.width[i] > this->max_width) {
+                tp.width[i] = this->max_width;
+            }
+        }
+    }
+    
+    concatThickPolylines(pp);
+    
+    /* Find the maximum width returned; we're going to use this for validating and 
+       filtering the output segments. */
+    double max_w = 0;
+    for (ThickPolylines::const_iterator it = pp.begin(); it != pp.end(); ++it)
+        max_w = fmaxf(max_w, *std::max_element(it->width.begin(), it->width.end()));
+
+
+    fusion_curve(pp);
+
+    // Aligned fusion: Fusion the bits at the end of lines by "increasing thickness"
+    // For that, we have to find other lines,
+    // and with a next point no more distant than the max width.
+    // Then, we can merge the bit from the first point to the second by following the mean.
+    //
+    main_fusion(pp);
+
+    //fusion right-angle corners.
+    fusion_corners(pp);
+
+    // Loop through all returned polylines in order to extend their endpoints to the 
+    //   expolygon boundaries (if done here, it may be cut later if not thick enough)
+    if (stop_at_min_width) {
+        extends_line_both_side(pp);
+    }
+
+    //reduce extrusion when it's too thin to be printable
+    remove_too_thin_extrusion(pp);
+
+    remove_too_thin_points(pp);
+
+    // Loop through all returned polylines in order to extend their endpoints to the 
+    //   expolygon boundaries
+    if (!stop_at_min_width) {
+        extends_line_both_side(pp);
+    }
+
+    concatenate_polylines_with_crossing(pp);
+
+    remove_too_short_polylines(pp, max_w * 2);
+
+    //TODO: reduce the flow at the intersection ( + ) points ?
+    ensure_not_overextrude(pp);
+
+    if (nozzle_diameter != min_width) {
+        grow_to_nozzle_diameter(pp, diff_ex(*this->bounds, this->expolygon));
+    }
+    if(this->taper_size != 0){
+        taper_ends(pp);
+    }
+
+    polylines_out.insert(polylines_out.end(), pp.begin(), pp.end());
+
+}
+
+ExtrusionEntityCollection
+discretize_variable_width(const ThickPolylines &polylines, ExtrusionRole role, Flow flow)
+{
+    // this value determines granularity of adaptive width, as G-code does not allow
+    // variable extrusion within a single move; this value shall only affect the amount
+    // of segments, and any pruning shall be performed before we apply this tolerance
+    const double tolerance = 4*SCALED_RESOLUTION;//scale_(0.05);
+    
+    ExtrusionEntityCollection coll;
+    for (const ThickPolyline &p : polylines) {
+        ExtrusionPaths paths;
+        ExtrusionPath path(role);
+        ThickLines lines = p.thicklines();
+        
+        for (int i = 0; i < (int)lines.size(); ++i) {
+            ThickLine& line = lines[i];
+            
+            const coordf_t line_len = line.length();
+            if (line_len < SCALED_EPSILON) continue;
+            
+            assert(line.a_width >= 0);
+            assert(line.b_width >= 0);
+            double thickness_delta = fabs(line.a_width - line.b_width);
+            if (thickness_delta > tolerance && ceil(thickness_delta / tolerance) > 2) {
+                const size_t segments = 1 + std::min(16000.0, ceil(thickness_delta / tolerance));
+                Points pp;
+                std::vector<coordf_t> width;
+                {
+                    for (size_t j = 0; j < segments; ++j) {
+                        pp.push_back(line.a.interpolate(((double)j) / segments, line.b));
+                        double percent_width = ((double)j) / (segments-1);
+                        width.push_back(line.a_width*(1 - percent_width) + line.b_width*percent_width);
+                    }
+                    pp.push_back(line.b);
+
+                    assert(pp.size() == segments + 1);
+                    assert(width.size() == segments);
+                }
+                
+                // delete this line and insert new ones
+                lines.erase(lines.begin() + i);
+                for (size_t j = 0; j < segments; ++j) {
+                    ThickLine new_line(pp[j], pp[j + 1]);
+                    new_line.a_width = width[j];
+                    new_line.b_width = width[j];
+                    lines.insert(lines.begin() + i + j, new_line);
+                }
+                
+                --i;
+                continue;
+            } else if (thickness_delta > 0) {
+                //create a middle point
+                ThickLine new_line(line.a.interpolate(0.5, line.b), line.b);
+                new_line.a_width = line.b_width;
+                new_line.b_width = line.b_width;
+                line.b = new_line.a;
+                line.b_width = line.a_width;
+                lines.insert(lines.begin() + i + 1, new_line);
+
+                --i;
+                continue;
+            }
+            //gapfill : we want to be able to fill the voids (touching the perimeters), so the spacing is what we want.
+            //thinwall: we want the extrusion to not go out of the polygon, so the width is what we want.
+            //  but we can't extrude with a negative spacing, so we have to gradually fall back to spacing if the width is too small.
+
+            // default: extrude a thin wall that doesn't go outside of the specified width.
+            coordf_t wanted_width = unscale(line.a_width);
+            if (role == erGapFill) {
+                // Convert from spacing to extrusion width based on the extrusion model
+                // of a square extrusion ended with semi circles.
+                wanted_width = unscale(line.a_width) + flow.height * (1. - 0.25 * PI);
+            } else if (unscale(line.a_width) < 2 * flow.height * (1. - 0.25 * PI)) {
+                //width (too) small, be sure to not extrude with negative spacing.
+                //we began to fall back to spacing gradually even before the spacing go into the negative
+                //  to make extrusion1 < extrusion2 if width1 < width2 even if width2 is too small. 
+                wanted_width = unscale(line.a_width)*0.35 + 1.3 * flow.height * (1. - 0.25 * PI);
+            }
+
+            if (path.polyline.points.empty()) {
+                flow.width = wanted_width;
+                path.polyline.append(line.a);
+                path.polyline.append(line.b);
+                path.mm3_per_mm = flow.mm3_per_mm();
+                path.width = flow.width;
+                path.height = flow.height;
+            } else {
+                thickness_delta = scale_(fabs(flow.width - wanted_width));
+                if (thickness_delta <= tolerance / 2) {
+                    // the width difference between this line and the current flow width is 
+                    // within the accepted tolerance
+                    path.polyline.append(line.b);
+                } else {
+                    // we need to initialize a new line
+                    paths.emplace_back(std::move(path));
+                    path = ExtrusionPath(role);
+                    --i;
+                }
+            }
+        }
+        if (path.polyline.is_valid())
+            paths.emplace_back(std::move(path));
+        // Append paths to collection.
+        if (!paths.empty()) {
+            if (paths.front().first_point().coincides_with(paths.back().last_point())) {
+                coll.append(ExtrusionLoop(paths));
+            } else {
+                //TODO: add them to an unsortable collection, to be able to keep the order anchor->outside
+                // and to keep the order of extrusion (do not stop extruding at intersection to turn => bad)
+                // BUT this need the ironing code to be able to have multiple collection inside each other
+                // BUT this need the ironing code to be able to have the no_sort flag to be useful
+                //ExtrusionEntityCollection unsortable_coll(paths);
+                //unsortable_coll.no_sort = true;
+                //coll.append(unsortable_coll);
+                coll.append(paths);
+            }
+        }
+    }
+    
+    return coll;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/MedialAxis.hpp b/xs/src/libslic3r/MedialAxis.hpp
new file mode 100644
index 0000000000..36efb49cae
--- /dev/null
+++ b/xs/src/libslic3r/MedialAxis.hpp
@@ -0,0 +1,120 @@
+#ifndef slic3r_MedialAxis_hpp_
+#define slic3r_MedialAxis_hpp_
+
+#include "libslic3r.h"
+#include "ExPolygon.hpp"
+#include "Polyline.hpp"
+#include "Geometry.hpp"
+#include "ExtrusionEntityCollection.hpp"
+#include "Flow.hpp"
+#include <vector>
+
+#include "boost/polygon/voronoi.hpp"
+using boost::polygon::voronoi_builder;
+using boost::polygon::voronoi_diagram;
+
+namespace Slic3r {
+
+/// This class is used to create single-line extrusion pattern with variable width to cover a ExPolygon.
+/// The ends can enter a boundary area if neded, and can have a taper at each end.
+/// The constructor initialize the mandatory variable.
+/// you must use the setter to add the opptional settings before calling build().
+class MedialAxis {
+    public:
+        /// _expolygon: the polygon to fill
+        /// _max_width : maximum width of the extrusion. _expolygon shouldn't have a spot where a circle diameter is higher than that (or almost).
+        /// _min_width : minimum width of the extrusion, every spot where a circle diameter is lower than that will be ignored (unless it's the tip of the extrusion)
+        /// _height: height of the extrusion, used to compute the difference between width and spacing.
+        MedialAxis(const ExPolygon &_expolygon, const coord_t _max_width, const coord_t _min_width, const coord_t _height)
+            : surface(_expolygon), max_width(_max_width), min_width(_min_width), height(_height),
+            bounds(&_expolygon), nozzle_diameter(_min_width), taper_size(0), stop_at_min_width(true){};
+
+        /// create the polylines_out collection of variable-width polyline to extrude.
+        void build(ThickPolylines &polylines_out);
+        /// You shouldn't use this method as it doesn't give you the variable width. Can be useful for debugging.
+        void build(Polylines &polylines);
+
+        /// optional parameter: anchor area in which the extrusion should extends into. Default : expolygon (no bound)
+        MedialAxis& use_bounds(const ExPolygon & _bounds) { this->bounds = &_bounds; return *this; }
+        /// optional parameter: the real minimum width : it will grow the width of every extrusion that has a width lower than that. Default : min_width (same min)
+        MedialAxis& use_min_real_width(const coord_t nozzle_diameter) { this->nozzle_diameter = nozzle_diameter; return *this; }
+        /// optional parameter: create a taper of this length at each end (inside a bound or not). Default : 0 (no taper)
+        MedialAxis& use_tapers(const coord_t taper_size) { this->taper_size = taper_size; return *this; }
+        /// optional parameter: if true, the entension inside the bounds can be cut if the width is too small. Default : true
+        MedialAxis& set_stop_at_min_width(const bool stop_at_min_width) { this->stop_at_min_width = stop_at_min_width; return *this; }
+
+    private:
+        /// Cache value: lines is here only to avoid passing it in argument of many methods. Initialized in polyline_from_voronoi.
+        Lines lines;
+
+        /// input polygon to fill
+        const ExPolygon& surface;
+        /// the copied expolygon from surface, it's modified in build() to simplify it. It's then used to create the voronoi diagram.
+        ExPolygon expolygon;
+        const ExPolygon* bounds;
+        /// maximum width of the extrusion. _expolygon shouldn't have a spot where a circle diameter is higher than that (or almost).
+        const coord_t max_width;
+        /// minimum width of the extrusion, every spot where a circle diameter is lower than that will be ignored (unless it's the tip of the extrusion)
+        const coord_t min_width;
+        /// height of the extrusion, used to compute the diufference between width and spacing.
+        const coord_t height;
+        /// Used to compute the real minimum width we can extrude. if != min_width, it activate grow_to_nozzle_diameter(). 
+        coord_t nozzle_diameter;
+        /// if != , it activates taper_ends(). Can use nozzle_diameter.
+        coord_t taper_size;
+        //if true, remove_too_* can shorten the bits created by extends_line.
+        bool stop_at_min_width;
+
+        //voronoi stuff
+        class VD : public voronoi_diagram<double> {
+        public:
+            typedef double                                          coord_type;
+            typedef boost::polygon::point_data<coordinate_type>     point_type;
+            typedef boost::polygon::segment_data<coordinate_type>   segment_type;
+            typedef boost::polygon::rectangle_data<coordinate_type> rect_type;
+        };
+        VD vd;
+        std::set<const VD::edge_type*> edges, valid_edges;
+        std::map<const VD::edge_type*, std::pair<coordf_t, coordf_t> > thickness;
+        void process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline);
+        bool validate_edge(const VD::edge_type* edge);
+        const Line& retrieve_segment(const VD::cell_type* cell) const;
+        const Point& retrieve_endpoint(const VD::cell_type* cell) const;
+        void polyline_from_voronoi(const Lines& voronoi_edges, ThickPolylines* polylines_out);
+
+        // functions called by build:
+
+        /// create a simplied version of surface, store it in expolygon
+        void simplify_polygon_frontier();
+        /// fusion little polylines created (by voronoi) on the external side of a curve inside the main polyline.
+        void fusion_curve(ThickPolylines &pp);
+        /// fusion polylines created by voronoi, where needed.
+        void main_fusion(ThickPolylines& pp);
+        /// like fusion_curve but for sharp angles like a square corner.
+        void fusion_corners(ThickPolylines &pp);
+        /// extends the polylines inside bounds, use extends_line on both end
+        void extends_line_both_side(ThickPolylines& pp);
+        /// extends the polylines inside bounds (anchors)
+        void extends_line(ThickPolyline& polyline, const ExPolygons& anchors, const coord_t join_width);
+        /// remove too thin bits at start & end of polylines
+        void remove_too_thin_extrusion(ThickPolylines& pp);
+        /// instead of keeping polyline split at each corssing, we try to create long strait polylines that can cross each other.
+        void concatenate_polylines_with_crossing(ThickPolylines& pp);
+        /// remove bits around points that are too thin (can be inside the polyline)
+        void remove_too_thin_points(ThickPolylines& pp);
+        /// delete polylines that are too short
+        void remove_too_short_polylines(ThickPolylines& pp, const coord_t min_size);
+        /// be sure we didn't try to push more plastic than the volume defined by surface * height can receive. If overextruded, reduce all widths by the correct %.
+        void ensure_not_overextrude(ThickPolylines& pp);
+        /// if nozzle_diameter > min_width, grow bits that are < width(nozzle_diameter) to width(nozzle_diameter) (don't activate that for gapfill)
+        void grow_to_nozzle_diameter(ThickPolylines& pp, const ExPolygons& anchors);
+        /// taper the ends of polylines (don't activate that for gapfill)
+        void taper_ends(ThickPolylines& pp);
+};
+    
+    /// create a ExtrusionEntityCollection from ThickPolylines, discretizing the variable width into little sections (of 4*SCALED_RESOLUTION length) where needed.
+    ExtrusionEntityCollection discretize_variable_width(const ThickPolylines &polylines, ExtrusionRole role, Flow flow);
+}
+
+
+#endif
diff --git a/xs/src/libslic3r/MultiPoint.cpp b/xs/src/libslic3r/MultiPoint.cpp
index d8ee1d09d6..4625bcffd2 100644
--- a/xs/src/libslic3r/MultiPoint.cpp
+++ b/xs/src/libslic3r/MultiPoint.cpp
@@ -158,6 +158,30 @@ MultiPoint::intersection(const Line& line, Point* intersection) const
     return false;
 }
 
+bool
+MultiPoint::first_intersection(const Line& line, Point* intersection) const
+{
+    bool   found = false;
+    double dmin  = 0.;
+    for (const Line &l : this->lines()) {
+        Point ip;
+        if (l.intersection(line, &ip)) {
+            if (! found) {
+                found = true;
+                dmin = ip.distance_to(line.a);
+                *intersection = ip;
+            } else {
+                double d = ip.distance_to(line.a);
+                if (d < dmin) {
+                    dmin = d;
+                    *intersection = ip;
+                }
+            }
+        }
+    }
+    return found;
+}
+
 std::string
 MultiPoint::dump_perl() const
 {
diff --git a/xs/src/libslic3r/MultiPoint.hpp b/xs/src/libslic3r/MultiPoint.hpp
index acb409b205..07abe553ac 100644
--- a/xs/src/libslic3r/MultiPoint.hpp
+++ b/xs/src/libslic3r/MultiPoint.hpp
@@ -31,6 +31,25 @@ class MultiPoint
 
     int find_point(const Point &point) const;
     bool has_boundary_point(const Point &point) const;
+    /// return the index of the closest point in this polygon in relation with "point"
+    /// \param point the point to compare with. 
+    /// \return the index of the closest point in the points vector, or max(size_t) if points is empty. 
+    size_t closest_point_index(const Point &point) const {
+        size_t idx = -1;
+        if (! this->points.empty()) {
+            idx = 0;
+            double dist_min = this->points.front().distance_to_sq(point);
+            for (size_t i = 1; i < this->points.size(); ++ i) {
+                double d = this->points[i].distance_to_sq(point);
+                if (d < dist_min) {
+                    dist_min = d;
+                    idx = i;
+                }
+            }
+        }
+        return idx;
+    }
+    const Point* closest_point(const Point &point) const { return this->points.empty() ? nullptr : &this->points[this->closest_point_index(point)]; }
     BoundingBox bounding_box() const;
     
     // Return true if there are exact duplicates.
@@ -43,6 +62,7 @@ class MultiPoint
     void append(const Points &points);
     void append(const Points::const_iterator &begin, const Points::const_iterator &end);
     bool intersection(const Line& line, Point* intersection) const;
+    bool first_intersection(const Line& line, Point* intersection) const;
     std::string dump_perl() const;
     
     static Points _douglas_peucker(const Points &points, const double tolerance);
diff --git a/xs/src/libslic3r/PerimeterGenerator.cpp b/xs/src/libslic3r/PerimeterGenerator.cpp
index e789a94dce..57e1fe4513 100644
--- a/xs/src/libslic3r/PerimeterGenerator.cpp
+++ b/xs/src/libslic3r/PerimeterGenerator.cpp
@@ -1,6 +1,7 @@
 #include "PerimeterGenerator.hpp"
 #include "ClipperUtils.hpp"
 #include "ExtrusionEntityCollection.hpp"
+#include "MedialAxis.hpp"
 #include <cmath>
 #include <cassert>
 
@@ -26,6 +27,9 @@ PerimeterGenerator::process()
     // solid infill
     coord_t ispacing            = this->solid_infill_flow.scaled_spacing();
     
+    //nozzle diameter
+    const double nozzle_diameter = this->print_config->nozzle_diameter.get_at(this->config->perimeter_extruder-1);
+    
     // Calculate the minimum required spacing between two adjacent traces.
     // This should be equal to the nominal flow spacing but we experiment
     // with some tolerance in order to avoid triggering medial axis when
@@ -47,8 +51,6 @@ PerimeterGenerator::process()
         // We consider overhang any part where the entire nozzle diameter is not supported by the
         // lower layer, so we take lower slices and offset them by half the nozzle diameter used 
         // in the current layer
-        double nozzle_diameter = this->print_config->nozzle_diameter.get_at(this->config->perimeter_extruder-1);
-        
         this->_lower_slices_p = offset(*this->lower_slices, scale_(+nozzle_diameter/2));
     }
     
@@ -85,44 +87,80 @@ PerimeterGenerator::process()
             for (int i = 0; i <= loop_number+1; ++i) {  // outer loop is 0
                 Polygons offsets;
                 if (i == 0) {
-                    // the minimum thickness of a single loop is:
-                    // ext_width/2 + ext_spacing/2 + spacing/2 + width/2
+                    // compute next onion, without taking care of thin_walls : destroy too thin areas.
+                    if (!this->config->thin_walls)
+                        offsets = offset(last, -(float)(ext_pwidth / 2));
+
+
+                    // look for thin walls
                     if (this->config->thin_walls) {
+                        // the minimum thickness of a single loop is:
+                        // ext_width/2 + ext_spacing/2 + spacing/2 + width/2
+                        //here, we shrink & grow by ext_min_spacing to remove areas where the current loop can't be extruded 
                         offsets = offset2(
                             last,
                             -(ext_pwidth/2 + ext_min_spacing/2 - 1),
-                            +(ext_min_spacing/2 - 1)
-                        );
-                    } else {
-                        offsets = offset(last, -ext_pwidth/2);
-                    }
-                    
-                    // look for thin walls
-                    if (this->config->thin_walls) {
-                        Polygons no_thin_zone = offset(offsets, +ext_pwidth/2);
-                        Polygons diffpp = diff(
-                            last,
-                            no_thin_zone,
-                            true  // medial axis requires non-overlapping geometry
-                        );
+                            +(ext_min_spacing/2 - 1));
+
+                        // detect edge case where a curve can be split in multiple small chunks.
+                        Polygons no_thin_onion = offset(last, -(float)(ext_pwidth / 2));
+                        float div = 2;
+                        while (no_thin_onion.size()>0 && offsets.size() > no_thin_onion.size() && no_thin_onion.size() + offsets.size() > 3) {
+                            div += 0.5;
+                            //use a sightly smaller offset2 spacing to try to improve the split, but with a little bit of over-extrusion
+                            Polygons next_onion_secondTry = offset2(
+                                last,
+                                -(float)(ext_pwidth / 2 + ext_min_spacing / div - 1),
+                                +(float)(ext_min_spacing / div - 1));
+                            if (offsets.size() >  next_onion_secondTry.size() * 1.1) {
+                                offsets = next_onion_secondTry;
+                            }
+                            if (div > 3) break;
+                        }
                         
                         // the following offset2 ensures almost nothing in @thin_walls is narrower than $min_width
                         // (actually, something larger than that still may exist due to mitering or other causes)
-                        coord_t min_width = scale_(this->ext_perimeter_flow.nozzle_diameter / 3);
-                        ExPolygons expp = offset2_ex(diffpp, -min_width/2, +min_width/2);
-						
-                         // compute a bit of overlap to anchor thin walls inside the print.
-                        ExPolygons anchor = intersection_ex(to_polygons(offset_ex(expp, (float)(ext_pwidth / 2))), no_thin_zone, true);
+                        coord_t min_width = scale_(this->config->thin_walls_min_width.get_abs_value(this->ext_perimeter_flow.nozzle_diameter));
                         
-                        // the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop
-                        for (ExPolygons::const_iterator ex = expp.begin(); ex != expp.end(); ++ex) {
-                            ExPolygons bounds = _clipper_ex(ClipperLib::ctUnion, (Polygons)*ex, to_polygons(anchor), true);
-							//search our bound
+                        Polygons no_thin_zone = offset(offsets, ext_pwidth/2, CLIPPER_OFFSET_SCALE, jtSquare, 3);
+                        // medial axis requires non-overlapping geometry
+                        Polygons thin_zones = diff(last, no_thin_zone, true);
+                        //don't use offset2_ex, because we don't want to merge the zones that have been separated.
+                        //a very little bit of overlap can be created here with other thin polygons, but it's more useful than worisome.
+                        ExPolygons half_thins = offset_ex(thin_zones, (float)(-min_width / 2));
+                        //simplify them
+                        for (ExPolygon &half_thin : half_thins) {
+                            half_thin.remove_point_too_near(SCALED_RESOLUTION);
+                        }
+                        //we push the bits removed and put them into what we will use as our anchor
+                        if (half_thins.size() > 0) {
+                            no_thin_zone = diff(last, to_polygons(offset_ex(half_thins, (float)(min_width / 2) - SCALED_EPSILON)), true);
+                        }
+                        // compute a bit of overlap to anchor thin walls inside the print.
+                        for (ExPolygon &half_thin : half_thins) {
+                            //growing back the polygon
+                            ExPolygons thin = offset_ex(half_thin, (float)(min_width / 2));
+                            assert(thin.size() <= 1);
+                            if (thin.empty()) continue;
+                            double overlap = (coord_t)scale_(this->config->thin_walls_overlap.get_abs_value(this->ext_perimeter_flow.nozzle_diameter));
+                            ExPolygons anchor = intersection_ex(
+                                to_polygons(offset_ex(half_thin, (float)(min_width / 2 + overlap), CLIPPER_OFFSET_SCALE, jtSquare, 3)), 
+                                no_thin_zone, true);
+                            ExPolygons bounds = _clipper_ex(ClipperLib::ctUnion, to_polygons(thin), to_polygons(anchor), true);
                             for (ExPolygon &bound : bounds) {
-                                if (!intersection_ex(*ex, bound).empty()) {
-                                    // the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop
-                                    ex->medial_axis(bound, ext_pwidth + ext_pspacing2, min_width, &thin_walls);
-                                    continue;
+                                if (!intersection_ex(thin[0], bound).empty()) {
+                                    //be sure it's not too small to extrude reliably
+                                    thin[0].remove_point_too_near(SCALED_RESOLUTION);
+                                    if (thin[0].area() > min_width*(ext_pwidth + ext_pspacing2)) {
+                                        bound.remove_point_too_near(SCALED_RESOLUTION);
+                                        // the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop (*1.1 because of circles approx.)
+                                        Slic3r::MedialAxis ma{ thin[0], (ext_pwidth + ext_pspacing2)*1.1, min_width, coord_t(this->layer_height) };
+                                        ma.use_bounds(bound)
+                                            .use_min_real_width((coord_t)scale_(this->ext_perimeter_flow.nozzle_diameter))
+                                            .use_tapers(overlap)
+                                            .build(thin_walls);
+                                    }
+                                    break;
                                 }
                             }
                         }
@@ -284,6 +322,8 @@ PerimeterGenerator::process()
             
             // collapse 
             double min = 0.2*pwidth * (1 - INSET_OVERLAP_TOLERANCE);
+            //be sure we don't gapfill where the perimeters are already touching each other (negative spacing).
+            min = std::max(min, double(Flow::new_from_spacing(EPSILON, nozzle_diameter, this->layer_height, false).scaled_width()));
             double max = 2*pspacing;
             ExPolygons gaps_ex = diff_ex(
                 offset2(gaps, -min/2, +min/2),
@@ -292,12 +332,15 @@ PerimeterGenerator::process()
             );
             
             ThickPolylines polylines;
-            for (ExPolygons::const_iterator ex = gaps_ex.begin(); ex != gaps_ex.end(); ++ex)
-                ex->medial_axis(*ex, max, min, &polylines);
-            
+            for (const ExPolygon &ex : gaps_ex) {
+                //remove too small gaps that are too hard to fill.
+                //ie one that are smaller than an extrusion with width of min and a length of max.
+                if (ex.area() > min*max) {
+                    MedialAxis{ex, coord_t(max),coord_t(min), coord_t(this->layer_height)}.build(polylines);
+                }
+            }
             if (!polylines.empty()) {
-                ExtrusionEntityCollection gap_fill = this->_variable_width(polylines, 
-                    erGapFill, this->solid_infill_flow);
+                ExtrusionEntityCollection gap_fill = discretize_variable_width(polylines, erGapFill, this->solid_infill_flow);
                 
                 this->gap_fill->append(gap_fill.entities);
             
@@ -421,8 +464,7 @@ PerimeterGenerator::_traverse_loops(const PerimeterGeneratorLoops &loops,
     
     // append thin walls to the nearest-neighbor search (only for first iteration)
     if (!thin_walls.empty()) {
-        ExtrusionEntityCollection tw = this->_variable_width
-            (thin_walls, erExternalPerimeter, this->ext_perimeter_flow);
+        ExtrusionEntityCollection tw = discretize_variable_width(thin_walls, erExternalPerimeter, this->ext_perimeter_flow);
         
         coll.append(tw.entities);
         thin_walls.clear();
@@ -463,110 +505,6 @@ PerimeterGenerator::_traverse_loops(const PerimeterGeneratorLoops &loops,
     return entities;
 }
 
-ExtrusionEntityCollection
-PerimeterGenerator::_variable_width(const ThickPolylines &polylines, ExtrusionRole role, Flow flow) const
-{
-    // this value determines granularity of adaptive width, as G-code does not allow
-    // variable extrusion within a single move; this value shall only affect the amount
-    // of segments, and any pruning shall be performed before we apply this tolerance
-    const double tolerance = scale_(0.05);
-    
-    ExtrusionEntityCollection coll;
-    for (ThickPolylines::const_iterator p = polylines.begin(); p != polylines.end(); ++p) {
-        ExtrusionPaths paths;
-        ExtrusionPath path(role);
-        ThickLines lines = p->thicklines();
-        
-        for (int i = 0; i < (int)lines.size(); ++i) {
-            const ThickLine& line = lines[i];
-            
-            const coordf_t line_len = line.length();
-            if (line_len < SCALED_EPSILON) continue;
-            
-            double thickness_delta = fabs(line.a_width - line.b_width);
-            if (thickness_delta > tolerance) {
-                const size_t segments = ceil(thickness_delta / tolerance);
-                const coordf_t seg_len = line_len / segments;
-                Points pp;
-                std::vector<coordf_t> width;
-                {
-                    pp.push_back(line.a);
-                    width.push_back(line.a_width);
-                    for (size_t j = 1; j < segments; ++j) {
-                        pp.push_back(line.point_at(j*seg_len));
-                        
-                        coordf_t w = line.a_width + (j*seg_len) * (line.b_width-line.a_width) / line_len;
-                        width.push_back(w);
-                        width.push_back(w);
-                    }
-                    pp.push_back(line.b);
-                    width.push_back(line.b_width);
-                    
-                    assert(pp.size() == segments + 1);
-                    assert(width.size() == segments*2);
-                }
-                
-                // delete this line and insert new ones
-                lines.erase(lines.begin() + i);
-                for (size_t j = 0; j < segments; ++j) {
-                    ThickLine new_line(pp[j], pp[j+1]);
-                    new_line.a_width = width[2*j];
-                    new_line.b_width = width[2*j+1];
-                    lines.insert(lines.begin() + i + j, new_line);
-                }
-                
-                --i;
-                continue;
-            }
-            
-            const double w = fmax(line.a_width, line.b_width);
-            
-            if (path.polyline.points.empty()) {
-                flow.width = unscale(w);
-                #ifdef SLIC3R_DEBUG
-                printf("  filling %f gap\n", flow.width);
-                #endif
-                
-                // make sure we don't include too thin segments which
-                // may cause even slightly negative mm3_per_mm because of floating point math
-                path.mm3_per_mm  = flow.mm3_per_mm();
-                if (path.mm3_per_mm < EPSILON) continue;
-                
-                path.width       = flow.width;
-                path.height      = flow.height;
-                path.polyline.append(line.a);
-                path.polyline.append(line.b);
-            } else {
-                thickness_delta = fabs(scale_(flow.width) - w);
-                if (thickness_delta <= tolerance) {
-                    // the width difference between this line and the current flow width is 
-                    // within the accepted tolerance
-                
-                    path.polyline.append(line.b);
-                } else {
-                    // we need to initialize a new line
-                    paths.push_back(path);
-                    path = ExtrusionPath(role);
-                    --i;
-                }
-            }
-        }
-        if (path.polyline.is_valid())
-            paths.push_back(path);
-        
-        // append paths to collection
-        if (!paths.empty()) {
-            if (paths.front().first_point().coincides_with(paths.back().last_point())) {
-                coll.append(ExtrusionLoop(paths));
-            } else {
-                coll.append(paths);
-            }
-        }
-    }
-    
-    return coll;
-}
-
 bool
 PerimeterGeneratorLoop::is_internal_contour() const
 {
diff --git a/xs/src/libslic3r/PerimeterGenerator.hpp b/xs/src/libslic3r/PerimeterGenerator.hpp
index 0e7fbd3e4b..4bd7802c3b 100644
--- a/xs/src/libslic3r/PerimeterGenerator.hpp
+++ b/xs/src/libslic3r/PerimeterGenerator.hpp
@@ -86,8 +86,6 @@ class PerimeterGenerator {
     
     ExtrusionEntityCollection _traverse_loops(const PerimeterGeneratorLoops &loops,
         ThickPolylines &thin_walls) const;
-    ExtrusionEntityCollection _variable_width
-        (const ThickPolylines &polylines, ExtrusionRole role, Flow flow) const;
 };
 
 }
diff --git a/xs/src/libslic3r/Point.cpp b/xs/src/libslic3r/Point.cpp
index bd5fcbda7f..7381f9e832 100644
--- a/xs/src/libslic3r/Point.cpp
+++ b/xs/src/libslic3r/Point.cpp
@@ -193,6 +193,14 @@ Point::distance_to(const Point &point) const
     return sqrt(dx*dx + dy*dy);
 }
 
+double
+Point::distance_to_sq(const Point &point) const
+{
+    double dx = ((double)point.x - this->x);
+    double dy = ((double)point.y - this->y);
+    return dx*dx + dy*dy;
+}
+
 /* distance to the closest point of line */
 double
 Point::distance_to(const Line &line) const
@@ -268,8 +276,8 @@ Point::projection_onto(const MultiPoint &poly) const
     for (Lines::const_iterator line = lines.begin(); line != lines.end(); ++line) {
         Point point_temp = this->projection_onto(*line);
         if (this->distance_to(point_temp) < running_min) {
-	        running_projection = point_temp;
-	        running_min = this->distance_to(running_projection);
+            running_projection = point_temp;
+            running_min = this->distance_to(running_projection);
         }
     }
     return running_projection;
@@ -303,6 +311,19 @@ Point::projection_onto(const Line &line) const
     }
 }
 
+/// This method create a new point on the line defined by this and p2.
+/// The new point is place at position defined by |p2-this| * percent, starting from this
+/// \param percent the proportion of the segment length to place the point
+/// \param p2 the second point, forming a segment with this
+/// \return a new point, == this if percent is 0 and == p2 if percent is 1
+Point Point::interpolate(const double percent, const Point &p2) const
+{
+    Point p_out;
+    p_out.x = this->x*(1 - percent) + p2.x*(percent);
+    p_out.y = this->y*(1 - percent) + p2.y*(percent);
+    return p_out;
+}
+
 Point
 Point::negative() const
 {
diff --git a/xs/src/libslic3r/Point.hpp b/xs/src/libslic3r/Point.hpp
index 45d6ba33a8..a02349f400 100644
--- a/xs/src/libslic3r/Point.hpp
+++ b/xs/src/libslic3r/Point.hpp
@@ -47,6 +47,7 @@ class Point
     static Point new_scale(Pointf p);
     bool operator==(const Point& rhs) const;
     bool operator!=(const Point& rhs) const { return !(*this == rhs); }
+    bool operator<(const Point& rhs) const { return this->x < rhs.x || (this->x == rhs.x && this->y < rhs.y); }
     std::string wkt() const;
     std::string dump_perl() const;
     void scale(double factor);
@@ -73,6 +74,7 @@ class Point
     bool nearest_point(const Points &points, Point* point) const;
     bool nearest_waypoint(const Points &points, const Point &dest, Point* point) const;
     double distance_to(const Point &point) const;
+    double distance_to_sq(const Point &point) const;
     double distance_to(const Line &line) const;
     double perp_distance_to(const Line &line) const;
     double ccw(const Point &p1, const Point &p2) const;
@@ -80,6 +82,7 @@ class Point
     double ccw_angle(const Point &p1, const Point &p2) const;
     Point projection_onto(const MultiPoint &poly) const;
     Point projection_onto(const Line &line) const;
+    Point interpolate(const double percent, const Point &p) const;
     Point negative() const;
     Vector vector_to(const Point &point) const;
     void align_to_grid(const Point &spacing, const Point &base = Point(0,0));
@@ -129,6 +132,8 @@ class Pointf
     bool operator==(const Pointf& rhs) const;
     bool coincides_with_epsilon(const Pointf& rhs) const;
     Pointf& operator/=(const double& scalar); 
+    bool operator!=(const Pointf &rhs) const { return ! (*this == rhs); }
+    bool operator< (const Pointf& rhs) const { return this->x < rhs.x || (this->x == rhs.x && this->y < rhs.y); }
 
     std::string wkt() const;
     std::string dump_perl() const;
@@ -143,8 +148,13 @@ class Pointf
 
 Pointf operator+(const Pointf& point1, const Pointf& point2);
 Pointf operator/(const Pointf& point1, const double& scalar);
-
-std::ostream& operator<<(std::ostream &stm, const Pointf3 &pointf3);
+inline Pointf operator*(double scalar, const Pointf& p) { return Pointf(scalar * p.x, scalar * p.y); }
+inline Pointf operator*(const Pointf& p, double scalar) { return Pointf(scalar * p.x, scalar * p.y); }
+inline Vectorf normalize(const Vectorf& v)
+{
+    coordf_t len = sqrt((v.x*v.x) + (v.y*v.y));
+    return (len != 0.0) ? 1.0 / len * v : Vectorf(0.0, 0.0);
+}
 
 class Pointf3 : public Pointf
 {
@@ -162,6 +172,8 @@ class Pointf3 : public Pointf
     Vectorf3 vector_to(const Pointf3 &point) const;
 };
 
+std::ostream& operator<<(std::ostream &stm, const Pointf3 &pointf3);
+
 template <class T>
 inline Points
 to_points(const std::vector<T> &items)
diff --git a/xs/src/libslic3r/Polyline.cpp b/xs/src/libslic3r/Polyline.cpp
index b80335073f..f0f33dc55c 100644
--- a/xs/src/libslic3r/Polyline.cpp
+++ b/xs/src/libslic3r/Polyline.cpp
@@ -6,6 +6,7 @@
 #include "ClipperUtils.hpp"
 #include <iostream>
 #include <utility>
+#include <algorithm>
 
 namespace Slic3r {
 
@@ -170,9 +171,9 @@ Polyline::split_at(const Point &point, Polyline* p1, Polyline* p2) const
     for (Lines::const_iterator line = lines.begin(); line != lines.end(); ++line) {
         Point p_tmp = point.projection_onto(*line);
         if (point.distance_to(p_tmp) < min) {
-	        p = p_tmp;
-	        min = point.distance_to(p);
-	        line_idx = line - lines.begin();
+            p = p_tmp;
+            min = point.distance_to(p);
+            line_idx = line - lines.begin();
         }
     }
     
@@ -235,8 +236,8 @@ ThickPolyline::thicklines() const
         lines.reserve(this->points.size() - 1);
         for (size_t i = 0; i < this->points.size()-1; ++i) {
             ThickLine line(this->points[i], this->points[i+1]);
-            line.a_width = this->width[2*i];
-            line.b_width = this->width[2*i+1];
+            line.a_width = this->width[i];
+            line.b_width = this->width[i + 1];
             lines.push_back(line);
         }
     }
@@ -251,4 +252,105 @@ ThickPolyline::reverse()
     std::swap(this->endpoints.first, this->endpoints.second);
 }
 
+void
+concatThickPolylines(ThickPolylines& pp) {
+    bool changes = true;
+    while (changes){
+        changes = false;
+        //concat polyline if only 2 polyline at a point
+        for (size_t i = 0; i < pp.size(); ++i) {
+            ThickPolyline *polyline = &pp[i];
+            if (polyline->first_point().coincides_with(polyline->last_point())) {
+                polyline->endpoints.first = false;
+                polyline->endpoints.second = false;
+                continue;
+            }
+
+            size_t id_candidate_first_point = -1;
+            size_t id_candidate_last_point = -1;
+            size_t nbCandidate_first_point = 0;
+            size_t nbCandidate_last_point = 0;
+            // find another polyline starting here
+            for (size_t j = 0; j < pp.size(); ++j) {
+                if (j == i) continue;
+                ThickPolyline *other = &pp[j];
+                if (other->first_point().coincides_with(other->last_point())) continue;
+                if (polyline->last_point().coincides_with(other->last_point())) {
+                    //other->reverse();
+                    id_candidate_last_point = j;
+                    nbCandidate_last_point++;
+                }
+                if (polyline->last_point().coincides_with(other->first_point())) {
+                    id_candidate_last_point = j;
+                    nbCandidate_last_point++;
+                }
+                if (polyline->first_point().coincides_with(other->last_point())) {
+                    id_candidate_first_point = j;
+                    nbCandidate_first_point++;
+                }
+                if (polyline->first_point().coincides_with(other->first_point())) {
+                    id_candidate_first_point = j;
+                    nbCandidate_first_point++;
+                    //other->reverse();
+                }
+            }
+            if (id_candidate_last_point == id_candidate_first_point && nbCandidate_first_point == 1 && nbCandidate_last_point == 1) {
+                if (polyline->first_point().coincides_with(pp[id_candidate_first_point].first_point())) pp[id_candidate_first_point].reverse();
+                // it's a trap! it's a  loop!
+                polyline->points.insert(polyline->points.end(), pp[id_candidate_first_point].points.begin() + 1, pp[id_candidate_first_point].points.end());
+                polyline->width.insert(polyline->width.end(), pp[id_candidate_first_point].width.begin() + 1, pp[id_candidate_first_point].width.end());
+                pp.erase(pp.begin() + id_candidate_first_point);
+                changes = true;
+                polyline->endpoints.first = false;
+                polyline->endpoints.second = false;
+            } else {
+
+                if (nbCandidate_first_point == 1) {
+                    if (polyline->first_point().coincides_with(pp[id_candidate_first_point].first_point())) pp[id_candidate_first_point].reverse();
+                    //concat at front
+                    polyline->width[0] = std::max(polyline->width.front(), pp[id_candidate_first_point].width.back());
+                    polyline->points.insert(polyline->points.begin(), pp[id_candidate_first_point].points.begin(), pp[id_candidate_first_point].points.end() - 1);
+                    polyline->width.insert(polyline->width.begin(), pp[id_candidate_first_point].width.begin(), pp[id_candidate_first_point].width.end() - 1);
+                    polyline->endpoints.first = pp[id_candidate_first_point].endpoints.first;
+                    pp.erase(pp.begin() + id_candidate_first_point);
+                    changes = true;
+                    if (id_candidate_first_point < i) {
+                        i--;
+                        polyline = &pp[i];
+                    }
+                    if (id_candidate_last_point > id_candidate_first_point) {
+                        id_candidate_last_point--;
+                    }
+                } else if (nbCandidate_first_point == 0) {
+                    //update endpoint
+                    polyline->endpoints.first = true;
+                }
+                if (nbCandidate_last_point == 1) {
+                    if (polyline->last_point().coincides_with(pp[id_candidate_last_point].last_point())) pp[id_candidate_last_point].reverse();
+                    //concat at back
+                    polyline->width[polyline->width.size() - 1] = std::max(polyline->width.back(), pp[id_candidate_last_point].width.front());
+                    polyline->points.insert(polyline->points.end(), pp[id_candidate_last_point].points.begin() + 1, pp[id_candidate_last_point].points.end());
+                    polyline->width.insert(polyline->width.end(), pp[id_candidate_last_point].width.begin() + 1, pp[id_candidate_last_point].width.end());
+                    polyline->endpoints.second = pp[id_candidate_last_point].endpoints.second;
+                    pp.erase(pp.begin() + id_candidate_last_point);
+                    changes = true;
+                    if (id_candidate_last_point < i) {
+                        i--;
+                        polyline = &pp[i];
+                    }
+                } else if (nbCandidate_last_point == 0) {
+                    //update endpoint
+                    polyline->endpoints.second = true;
+                }
+
+                if (polyline->last_point().coincides_with(polyline->first_point())) {
+                    //the concat has created a loop : update endpoints
+                    polyline->endpoints.first = false;
+                    polyline->endpoints.second = false;
+                }
+            }
+        }
+    }
+}
+
 }
diff --git a/xs/src/libslic3r/Polyline.hpp b/xs/src/libslic3r/Polyline.hpp
index 49e261bdf7..db4111db3b 100644
--- a/xs/src/libslic3r/Polyline.hpp
+++ b/xs/src/libslic3r/Polyline.hpp
@@ -36,15 +36,26 @@ class Polyline : public MultiPoint {
     Polygons grow(double delta, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtSquare, double miterLimit = 3.0) const;
 };
 
+
+/// ThickPolyline : a polyline with a width for each point
+/// This class has a vector of coordf_t, it must be the same size as points.
+/// it's used to store the size of the line at this point.
+/// Also, the endpoint let us know if the front() and back() of the polyline 
+/// join something or is a dead-end.
 class ThickPolyline : public Polyline {
     public:
+    /// width size must be == point size
     std::vector<coordf_t> width;
+    /// if true => it's an endpoint, if false it join an other ThickPolyline. first is at front(), second is at back()
     std::pair<bool,bool> endpoints;
     ThickPolyline() : endpoints(std::make_pair(false, false)) {};
     ThickLines thicklines() const;
     void reverse();
 };
 
+/// concatenate poylines if possible and refresh the endpoints
+void concatThickPolylines(ThickPolylines &polylines);
+
 inline Polylines
 to_polylines(const Polygons &polygons)
 {
diff --git a/xs/src/libslic3r/PrintConfig.cpp b/xs/src/libslic3r/PrintConfig.cpp
index 9b72137a23..a8e68b3623 100644
--- a/xs/src/libslic3r/PrintConfig.cpp
+++ b/xs/src/libslic3r/PrintConfig.cpp
@@ -821,16 +821,16 @@ PrintConfigDef::PrintConfigDef()
     def->default_value = new ConfigOptionInt(100);
 
     def = this->add("max_layer_height", coFloats);
-	def->label = "Max";
-	def->tooltip = "This is the highest printable layer height for this extruder and limits the resolution for adaptive slicing. Typical values are slightly smaller than nozzle_diameter.";
-	def->sidetext = "mm";
-	def->cli = "max-layer-height=f@";
-	def->min = 0;
-	{
-		ConfigOptionFloats* opt = new ConfigOptionFloats();
-		opt->values.push_back(0.3);
-		def->default_value = opt;
-	}
+    def->label = "Max";
+    def->tooltip = "This is the highest printable layer height for this extruder and limits the resolution for adaptive slicing. Typical values are slightly smaller than nozzle_diameter.";
+    def->sidetext = "mm";
+    def->cli = "max-layer-height=f@";
+    def->min = 0;
+    {
+        ConfigOptionFloats* opt = new ConfigOptionFloats();
+        opt->values.push_back(0.3);
+        def->default_value = opt;
+    }
 
     def = this->add("max_print_speed", coFloat);
     def->label = __TRANS("Max print speed");
@@ -1619,12 +1619,28 @@ PrintConfigDef::PrintConfigDef()
     }
     
     def = this->add("thin_walls", coBool);
-    def->label = __TRANS("Detect thin walls");
+    def->label = __TRANS("");
     def->category = __TRANS("Layers and Perimeters");
     def->tooltip = __TRANS("Detect single-width walls (parts where two extrusions don't fit and we need to collapse them into a single trace).");
     def->cli = "thin-walls!";
     def->default_value = new ConfigOptionBool(true);
 
+    def = this->add("thin_walls_min_width", coFloatOrPercent);
+    def->label = __TRANS("min width");
+    def->category = __TRANS("Layers and Perimeters");
+    def->tooltip = __TRANS("Minimum width for the extrusion to be extruded (widths lower than the nozzle diameter will be over-extruded at the nozzle diameter). Can be percent of the nozzle size.");
+    def->cli = "thin-walls-min-width=s";
+    def->min = 0;
+    def->default_value = new ConfigOptionFloatOrPercent(33, true);
+
+    def = this->add("thin_walls_overlap", coFloatOrPercent);
+    def->label = __TRANS("overlap");
+    def->category = __TRANS("Layers and Perimeters");
+    def->tooltip = __TRANS("Overlap between the thin walls and the perimeters. Can be a % of the external perimeter width (default 50%)");
+    def->cli = "thin-walls-overlap=s";
+    def->min = 0;
+    def->default_value = new ConfigOptionFloatOrPercent(50, true);
+
     def = this->add("threads", coInt);
     def->label = __TRANS("Threads");
     def->tooltip = __TRANS("Threads are used to parallelize long-running tasks. Optimal threads number is slightly above the number of available cores/processors.");
diff --git a/xs/src/libslic3r/PrintConfig.hpp b/xs/src/libslic3r/PrintConfig.hpp
index 9597814845..1376f900fb 100644
--- a/xs/src/libslic3r/PrintConfig.hpp
+++ b/xs/src/libslic3r/PrintConfig.hpp
@@ -269,6 +269,8 @@ class PrintRegionConfig : public virtual StaticPrintConfig
     ConfigOptionInt                 solid_infill_every_layers;
     ConfigOptionFloatOrPercent      solid_infill_speed;
     ConfigOptionBool                thin_walls;
+    ConfigOptionFloatOrPercent      thin_walls_min_width;
+    ConfigOptionFloatOrPercent      thin_walls_overlap;
     ConfigOptionFloatOrPercent      top_infill_extrusion_width;
     ConfigOptionEnum<InfillPattern> top_infill_pattern;
     ConfigOptionInt                 top_solid_layers;
@@ -311,6 +313,8 @@ class PrintRegionConfig : public virtual StaticPrintConfig
         OPT_PTR(solid_infill_every_layers);
         OPT_PTR(solid_infill_speed);
         OPT_PTR(thin_walls);
+        OPT_PTR(thin_walls_min_width);
+        OPT_PTR(thin_walls_overlap);
         OPT_PTR(top_infill_extrusion_width);
         OPT_PTR(top_infill_pattern);
         OPT_PTR(top_solid_infill_speed);
diff --git a/xs/src/libslic3r/PrintRegion.cpp b/xs/src/libslic3r/PrintRegion.cpp
index 6c5c6ed467..bda491383c 100644
--- a/xs/src/libslic3r/PrintRegion.cpp
+++ b/xs/src/libslic3r/PrintRegion.cpp
@@ -82,6 +82,8 @@ PrintRegion::invalidate_state_by_config(const PrintConfigBase &config)
             || opt_key == "first_layer_extrusion_width"
             || opt_key == "perimeter_extrusion_width"
             || opt_key == "thin_walls"
+            || opt_key == "thin_walls_min_width"
+            || opt_key == "thin_walls_overlap"
             || opt_key == "external_perimeters_first") {
             steps.insert(posPerimeters);
         } else if (opt_key == "first_layer_extrusion_width") {