thin_walls : medial axis improvementsmerging

It's now an intensive post-processing of the raw voronoi diagram.
It picks each crossing and try to merge branch where it makes sense, updating the width.
+ new tests (it fail the all medial axis segments of a semicircumference have the same orientation, but it's intended)
+ filter for too small thin walls
+ edge case for too many thin/thick chunks

note: the algo do not know "the good direction". If the thing is more wide that long, it will extrude side-way.

8 files changed, 675 insertions, 123 deletions

diff --git a/t/thin.t b/t/thin.t
index 9147236ee..d091117a2 100644
--- a/t/thin.t
+++ b/t/thin.t
@@ -1,4 +1,4 @@
-use Test::More tests => 23;
+use Test::More tests => 28;
 use strict;
 use warnings;
 
@@ -110,6 +110,28 @@ if (0) {
 
 {
     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],
+    ));
+    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 $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],
+    ));
+    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
+	
+}
+
+{
+    my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
         [100, 100],
         [120, 100],
         [112, 200],
diff --git a/xs/src/libslic3r/ExPolygon.cpp b/xs/src/libslic3r/ExPolygon.cpp
index 0677a9a9a..9a4036ba3 100644
--- a/xs/src/libslic3r/ExPolygon.cpp
+++ b/xs/src/libslic3r/ExPolygon.cpp
@@ -206,8 +206,101 @@ void ExPolygon::simplify(double tolerance, ExPolygons* expolygons) const
     append(*expolygons, this->simplify(tolerance));
 }
 
+/// remove point that are at SCALED_EPSILON * 2 distance.
+void remove_point_too_near(ThickPolyline* to_reduce) {
+        const int32_t smallest = SCALED_EPSILON * 2;
+        uint32_t id = 1;
+        while (id < to_reduce->points.size() - 2) {
+            uint32_t newdist = 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);
+            } 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 (uint32_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 
+        uint32_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);
+            Point new_point;
+            new_point.x = (coord_t)((double)(to_modify->points[idx_other - 1].x) * (1 - percent_dist));
+            new_point.x += (coord_t)((double)(to_modify->points[idx_other].x) * (percent_dist));
+            new_point.y = (coord_t)((double)(to_modify->points[idx_other - 1].y) * (1 - percent_dist));
+            new_point.y += (coord_t)((double)(to_modify->points[idx_other].y) * (percent_dist));
+            to_modify->width.insert(to_modify->width.begin() + idx_other, new_width);
+            to_modify->points.insert(to_modify->points.begin() + idx_other, new_point);
+        }
+    }
+}
+
+/// 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) {
+    double nearestDist = point.distance_to(contour.contour.points.front());
+    Point nearest = contour.contour.points.front();
+    uint32_t id_nearest = 0;
+    double nearDist = nearestDist;
+    Point near = nearest;
+    uint32_t id_near=0;
+    for (uint32_t id_point = 1; id_point < contour.contour.points.size(); ++id_point) {
+        if (nearestDist > point.distance_to(contour.contour.points[id_point])) {
+            nearestDist = point.distance_to(contour.contour.points[id_point]);
+            near = nearest;
+            nearest = contour.contour.points[id_point];
+            id_near = id_nearest;
+            id_nearest = id_point;
+        }
+    }
+    double angle = 0;
+    Point point_before = id_nearest == 0 ? contour.contour.points.back() : contour.contour.points[id_nearest - 1];
+    Point point_after = id_nearest == contour.contour.points.size()-1 ? contour.contour.points.front() : contour.contour.points[id_nearest + 1];
+    //compute angle
+    angle = min(nearest.ccw_angle(point_before, point_after), nearest.ccw_angle(point_after, point_before));
+    //compute the diff from 90°
+    angle = abs(angle - PI / 2);
+    if (near != nearest && max(nearestDist, nearDist) + SCALED_EPSILON < nearest.distance_to(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 = min(nearest.ccw_angle(point_before, point_after), nearest.ccw_angle(point_after, point_before));
+        angle2 = abs(angle - PI / 2);
+        angle = (angle + angle2) / 2;
+    }
+
+    return 1-(angle/(PI/2));
+}
+
 void
-ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_width, ThickPolylines* polylines) const
+ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_width, ThickPolylines* polylines, double height) const
 {
     // init helper object
     Slic3r::Geometry::MedialAxis ma(max_width, min_width, this);
@@ -217,12 +310,16 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
     ThickPolylines pp;
     ma.build(&pp);
     
-    /*
-    SVG svg("medial_axis.svg");
-    svg.draw(*this);
-    svg.draw(pp);
-    svg.Close();
-    */
+
+    //{
+    //    stringstream stri;
+    //    stri << "medial_axis" << id << ".svg";
+    //    SVG svg(stri.str());
+    //    svg.draw(bounds);
+    //    svg.draw(*this);
+    //    svg.draw(pp);
+    //    svg.Close();
+    //}
     
     /* Find the maximum width returned; we're going to use this for validating and 
        filtering the output segments. */
@@ -230,51 +327,152 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
     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()));
 
-    /*  Aligned fusion: Fusion the bits at the end of lines by "increasing thikness"
-    *   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.
-    */
+    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) { return a.length() < b.length(); });
+
+    // 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.
+    //
+    int id_f = 0;
     bool changes = true;
+
+    
     while (changes) {
         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;
             int best_idx = 0;
-
+            double dot_poly_branch = 0;
+            double dot_candidate_branch = 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())) {
+                } 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())) {
+                } else if (polyline.first_point().coincides_with(other.first_point())) {
+                } else if (polyline.last_point().coincides_with(other.first_point())) {
                     polyline.reverse();
                 } else {
                     continue;
                 }
-
-                //only consider the other if the next point is near us
+                //std::cout << " try : " << i << ":" << j << " : " << 
+                //    (polyline.points.size() < 2 && other.points.size() < 2) <<
+                //    (!polyline.endpoints.second || !other.endpoints.second) <<
+                //    ((polyline.points.back().distance_to(other.points.back())
+                //    + (polyline.width.back() + other.width.back()) / 4)
+                //    > max_width*1.05) <<
+                //    (abs(polyline.length() - other.length()) > max_width / 2) << "\n";
+
+                //// mergeable tests
                 if (polyline.points.size() < 2 && other.points.size() < 2) continue;
                 if (!polyline.endpoints.second || !other.endpoints.second) continue;
-                if (polyline.points.back().distance_to(other.points.back()) > max_width) continue;
-                if (polyline.points.size() != other.points.size()) 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 / 2) 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);
+                }
+                double coeffSizeOtherJ = 1;
+                if (other.width.back() == 0) {
+                    coeffSizeOtherJ = 0.1+0.9*get_coeff_from_angle_countour(other.points.back(), *this);
+                }
+                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).
                 Pointf v_poly(polyline.lines().front().vector().x, polyline.lines().front().vector().y);
                 v_poly.scale(1 / std::sqrt(v_poly.x*v_poly.x + v_poly.y*v_poly.y));
                 Pointf v_other(other.lines().front().vector().x, 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 = v_poly.x*v_other.x + v_poly.y*v_other.y;
+
+                // 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
+                bool find_main_branch = false;
+                int biggest_main_branch_id = 0;
+                int biggest_main_branch_length = 0;
+                for (size_t k = 0; k < pp.size(); ++k) {
+                    //std::cout << "try to find main : " << k << " ? " << i << " " << j << " ";
+                    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 = 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 = 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;
+                    }
+                }
+                if (!find_main_branch && biggest_main_branch_length == 0) {
+                    // nothing -> it's impossible!
+                    dot_poly_branch = 0.707;
+                    dot_candidate_branch = 0.707;
+                    //std::cout << "no main branch... impossible!!\n";
+                } else if (!find_main_branch && 
+                    (pp[biggest_main_branch_id].length() < polyline.length() || pp[biggest_main_branch_id].length() < other.length()) ){
+                    //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)
+                    Pointf v_poly(polyline.lines().front().vector().x, polyline.lines().front().vector().y);
+                    v_poly.scale(1 / std::sqrt(v_poly.x*v_poly.x + v_poly.y*v_poly.y));
+                    Pointf v_candid(other.lines().front().vector().x, other.lines().front().vector().y);
+                    v_candid.scale(1 / std::sqrt(v_candid.x*v_candid.x + v_candid.y*v_candid.y));
+                    Pointf v_branch(-pp[biggest_main_branch_id].lines().front().vector().x, -pp[biggest_main_branch_id].lines().front().vector().y);
+                    v_branch.scale(1 / std::sqrt(v_branch.x*v_branch.x + v_branch.y*v_branch.y));
+                    dot_poly_branch = v_poly.x*v_branch.x + v_poly.y*v_branch.y;
+                    dot_candidate_branch = v_candid.x*v_branch.x + v_candid.y*v_branch.y;
+                    if (dot_poly_branch < 0) dot_poly_branch = 0;
+                    if (dot_candidate_branch < 0) dot_candidate_branch = 0;
+                }
+                //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)
+                if (dot_poly_branch < 0.1 || dot_candidate_branch < 0.1 ||
+                    (polyline.length()>other.length() ? polyline.length() / other.length() : other.length() / polyline.length()) > 4) {
+                    continue;
+                }
                 if (other_dot > best_dot) {
                     best_candidate = &other;
                     best_idx = j;
@@ -282,20 +480,48 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
                 }
             }
             if (best_candidate != nullptr) {
-
-                //TODO: witch if polyline.size > best_candidate->size
-                //doesn't matter rright now because a if in the selection process prevent this.
-
+                // 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 = (get_coeff_from_angle_countour(polyline.points.back(), *this));
+                const double coeff_angle_candi = (get_coeff_from_angle_countour(best_candidate->points.back(), *this));
+
+                //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() / max(polyline.length(), best_candidate->length());
+                double weight_candi = 2 - best_candidate->length() / 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
                 unsigned int idx_point = 1;
-                while (idx_point < polyline.points.size() && polyline.points[idx_point].distance_to(best_candidate->points[idx_point]) < max_width) {
+                while (idx_point < min(polyline.points.size(), best_candidate->points.size())) {
                     //fusion
-                    polyline.points[idx_point].x += best_candidate->points[idx_point].x;
-                    polyline.points[idx_point].x /= 2;
-                    polyline.points[idx_point].y += best_candidate->points[idx_point].y;
-                    polyline.points[idx_point].y /= 2;
-                    polyline.width[idx_point] += best_candidate->width[idx_point];
+                    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 / max(dot_poly_branch, dot_candidate_branch);
+                    value_from_current_width += 0.5*best_candidate->width[idx_point] * dot_candidate_branch / 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(min(dot_poly_branch, dot_candidate_branch) / max(dot_poly_branch, dot_candidate_branch));
+                    polyline.width[idx_point] = value_from_current_width + value_from_dist;
+                    //failsafe
+                    if (polyline.width[idx_point] > max_width) polyline.width[idx_point] = max_width;
+
                     ++idx_point;
                 }
                 if (idx_point < best_candidate->points.size()) {
@@ -323,21 +549,22 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
 
                 //remove points that are the same or too close each other, ie simplify
                 for (unsigned int idx_point = 1; idx_point < polyline.points.size(); ++idx_point) {
-                    //distance of 1 is on the sclaed coordinates, so it correspond to SCALE_FACTOR, so it's very small
-                    if (polyline.points[idx_point - 1].distance_to(polyline.points[idx_point]) < 1) {
+                    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.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 (unsigned int idx_point = 0; idx_point < polyline.points.size(); ++idx_point) {
-                    //distance of 1 is on the sclaed coordinates, so it correspond to SCALE_FACTOR, so it's very small
                     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;
                     }
                 }
@@ -350,31 +577,90 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
 
                 pp.erase(pp.begin() + best_idx);
                 changes = true;
+                break;
             }
         }
+        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(); });
+        }
     }
 
+    // remove too small extrusion at start & end of polylines
+    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() < 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]) * percent_can_keep > max_width / 2
+                    && polyline.points.front().distance_to(polyline.points[1])* (1 - percent_can_keep) > max_width / 2) {
+                    //Can split => move the first point and assign a new weight.
+                    //the update of endpoints wil be performed in concatThickPolylines
+                    polyline.points.front().x = polyline.points.front().x + 
+                        (coord_t)((polyline.points[1].x - polyline.points.front().x) * percent_can_keep);
+                    polyline.points.front().y = polyline.points.front().y + 
+                        (coord_t)((polyline.points[1].y - polyline.points.front().y) * percent_can_keep);
+                    polyline.width.front() = min_width;
+                    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() < 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]) * percent_can_keep > max_width / 2
+                    && polyline.points.back().distance_to(polyline.points[polyline.points.size() - 2]) * (1-percent_can_keep) > max_width / 2) {
+                    //Can split => move the first point and assign a new weight.
+                    //the update of endpoints wil be performed in concatThickPolylines
+                    polyline.points.back().x = polyline.points.back().x + 
+                        (coord_t)((polyline.points[polyline.points.size() - 2].x - polyline.points.back().x) * percent_can_keep);
+                    polyline.points.back().y = polyline.points.back().y + 
+                        (coord_t)((polyline.points[polyline.points.size() - 2].y - polyline.points.back().y) * percent_can_keep);
+                    polyline.width.back() = min_width;
+                    changes = true;
+                    break;
+                }
+            }
+            polyline.points.erase(polyline.points.end()-1);
+            polyline.width.erase(polyline.width.end() - 1);
+            changes = true;
+        }
+        if (polyline.points.size() < 2) {
+            //remove self if too small
+            pp.erase(pp.begin() + i);
+            --i;
+        }
+    }
+    if (changes) concatThickPolylines(pp);
 
-    /* Loop through all returned polylines in order to extend their endpoints to the 
-       expolygon boundaries */
-    bool removed = false;
+    // Loop through all returned polylines in order to extend their endpoints to the 
+    //   expolygon boundaries
     for (size_t i = 0; i < pp.size(); ++i) {
         ThickPolyline& polyline = pp[i];
 
         // 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 */
+        // 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
         Point new_front = polyline.points.front();
         Point new_back = polyline.points.back();
         if (polyline.endpoints.first && !bounds.has_boundary_point(new_front)) {
-            Line line(polyline.points.front(), polyline.points[1]);
+            Line line(polyline.points[1], polyline.points.front());
 
             // prevent the line from touching on the other side, otherwise intersection() might return that solution
-            if (polyline.points.size() == 2) line.b = line.midpoint();
+            if (polyline.points.size() == 2) line.a = line.midpoint();
 
-            line.extend_start(max_width);
-            (void)bounds.contour.intersection(line, &new_front);
+            line.extend_end(max_width);
+            (void)bounds.contour.first_intersection(line, &new_front);
         }
         if (polyline.endpoints.second && !bounds.has_boundary_point(new_back)) {
             Line line(
@@ -386,7 +672,7 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
             if (polyline.points.size() == 2) line.a = line.midpoint();
             line.extend_end(max_width);
 
-            (void)bounds.contour.intersection(line, &new_back);
+            (void)bounds.contour.first_intersection(line, &new_back);
         }
         polyline.points.front() = new_front;
         polyline.points.back() = new_back;
@@ -394,7 +680,7 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
     }
 
 
-
+    // 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 
@@ -405,6 +691,7 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
         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.
         */
+    changes = false;
     for (size_t i = 0; i < pp.size(); ++i) {
         ThickPolyline& polyline = pp[i];
         if (polyline.endpoints.first && polyline.endpoints.second) continue; // optimization
@@ -441,32 +728,130 @@ ExPolygon::medial_axis(const ExPolygon &bounds, double max_width, double min_wid
         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.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()*2 - 2);
-                
+            assert(polyline.width.size() == polyline.points.size());
+            changes = true;
             pp.erase(pp.begin() + best_idx);
         }
     }
+    if (changes) concatThickPolylines(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 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;
-            removed = true;
-            continue;
+        // 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;
+            }
         }
+    }
+
+    //remove too short polyline
+    changes = true;
+    while (changes) {
+        changes = false;
+        
+        double shortest_size = max_w * 2;
+        int32_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 >= 0 && shortest_idx < pp.size()) {
+            pp.erase(pp.begin() + shortest_idx);
+            changes = true;
+        }
+        if (changes) concatThickPolylines(pp);
+    }
+
+    //TODO: reduce the flow at the intersection ( + ) points ?
+
+    //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 (auto hole = bounds.holes.begin(); hole != bounds.holes.end(); ++hole) {
+        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 (ThickLine l : polyline.thicklines()) {
+                l.a_width *= reduce_by;
+                l.b_width *= reduce_by;
+            }
+        }
+    }
     polylines->insert(polylines->end(), pp.begin(), pp.end());
 }
 
@@ -474,7 +859,7 @@ void
 ExPolygon::medial_axis(double max_width, double min_width, Polylines* polylines) const
 {
     ThickPolylines tp;
-    this->medial_axis(*this, max_width, min_width, &tp);
+    this->medial_axis(*this, max_width, min_width, &tp, max_width/2.0);
     polylines->insert(polylines->end(), tp.begin(), tp.end());
 }
 
diff --git a/xs/src/libslic3r/ExPolygon.hpp b/xs/src/libslic3r/ExPolygon.hpp
index fe6e70c4a..0d4db06a4 100644
--- a/xs/src/libslic3r/ExPolygon.hpp
+++ b/xs/src/libslic3r/ExPolygon.hpp
@@ -53,7 +53,7 @@ public:
     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 medial_axis(const ExPolygon &bounds, double max_width, double min_width, ThickPolylines* polylines, double height) const;
     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/Fill/FillBase.cpp b/xs/src/libslic3r/Fill/FillBase.cpp
index 4e850c34f..f8647a0b0 100644
--- a/xs/src/libslic3r/Fill/FillBase.cpp
+++ b/xs/src/libslic3r/Fill/FillBase.cpp
@@ -174,6 +174,9 @@ void Fill::fill_surface_extrusion(const Surface *surface, const FillParams &para
         //    (poylineVolume) / extrudedVolume,
         //    this->no_overlap_expolygons.size());
         if (extrudedVolume != 0 && poylineVolume != 0) multFlow = poylineVolume / extrudedVolume;
+        //failsafe, it can happen
+        if (multFlow > 1.3) multFlow = 1.3;
+        if (multFlow < 0.8) multFlow = 0.8;
     }
 
     // Save into layer.
diff --git a/xs/src/libslic3r/Geometry.cpp b/xs/src/libslic3r/Geometry.cpp
index c3d3b85e1..a94bb62ff 100644
--- a/xs/src/libslic3r/Geometry.cpp
+++ b/xs/src/libslic3r/Geometry.cpp
@@ -907,7 +907,7 @@ MedialAxis::build(ThickPolylines* polylines)
             polyline.endpoints.first = rpolyline.endpoints.second;
         }
         
-        assert(polyline.width.size() == polyline.points.size()*2 - 2);
+        assert(polyline.width.size() == polyline.points.size());
         
         // prevent loop endpoints from being extended
         if (polyline.first_point().coincides_with(polyline.last_point())) {
@@ -968,8 +968,8 @@ MedialAxis::process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* pol
             
             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;
@@ -1049,34 +1049,39 @@ MedialAxis::validate_edge(const VD::edge_type* edge)
         ? 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;
-    
+    //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
     if (w0 > this->max_width && w1 > this->max_width)
         return false;
     
diff --git a/xs/src/libslic3r/PerimeterGenerator.cpp b/xs/src/libslic3r/PerimeterGenerator.cpp
index 31de8a8d3..0c2b0655a 100644
--- a/xs/src/libslic3r/PerimeterGenerator.cpp
+++ b/xs/src/libslic3r/PerimeterGenerator.cpp
@@ -127,7 +127,7 @@ void PerimeterGenerator::process()
                                 //and we want at least 1 perimeter of overlap
                                 ExPolygons bridge = unsupported_filtered;
                                 unsupported_filtered = intersection_ex(offset_ex(unsupported_filtered, (float)(perimeter_spacing)), last);
-                                // remove from the bridge & support the small inmperfections of the union
+                                // remove from the bridge & support the small imperfections of the union
                                 ExPolygons bridge_and_support = offset2_ex(union_ex(bridge, support, true), perimeter_spacing/2, -perimeter_spacing/2);
                                 // make him flush with perimeter area
                                 unsupported_filtered = intersection_ex(offset_ex(unsupported_filtered, (float)(perimeter_spacing / 2)), bridge_and_support);
@@ -208,37 +208,66 @@ void PerimeterGenerator::process()
                 //this variable stored the nexyt onion
                 ExPolygons next_onion;
                 if (i == 0) {
-                    // the minimum thickness of a single loop is:
-                    // ext_width/2 + ext_spacing/2 + spacing/2 + width/2
-                    next_onion = this->config->thin_walls ?
-                        offset2_ex(
-                            last,
-                            -(float)(ext_perimeter_width / 2 + ext_min_spacing / 2 - 1),
-                            +(float)(ext_min_spacing / 2 - 1)) :
-                            offset_ex(last, -(float)(ext_perimeter_width / 2));
+                    // compute next onion, without taking care of thin_walls : destroy too thin areas.
+                    if (!this->config->thin_walls)
+                        next_onion = offset_ex(last, -(float)(ext_perimeter_width / 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
+
+                        next_onion = offset2_ex(
+                            last,
+                            -(float)(ext_perimeter_width / 2 + ext_min_spacing / 2 - 1),
+                            +(float)(ext_min_spacing / 2 - 1));
+
+                        // detect edge case where a curve can be split in multiple small chunks.
+                        ExPolygons no_thin_onion = offset_ex(last, -(float)(ext_perimeter_width / 2));
+                        if (no_thin_onion.size()>0 && next_onion.size() > 3 * no_thin_onion.size()) {
+                            //use a sightly smaller spacing to try to drastically improve the split
+                            ExPolygons next_onion_secondTry = offset2_ex(
+                                last,
+                                -(float)(ext_perimeter_width / 2 + ext_min_spacing / 2.5 - 1),
+                                +(float)(ext_min_spacing / 2.5 - 1));
+                            if (abs(((int32_t)next_onion.size()) - ((int32_t)no_thin_onion.size())) > 
+                                2*abs(((int32_t)next_onion_secondTry.size()) - ((int32_t)no_thin_onion.size()))) {
+                                next_onion = next_onion_secondTry;
+                            }
+                        }
+
                         // 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 = (coord_t)scale_(this->ext_perimeter_flow.nozzle_diameter / 3);
                         
-                        Polygons no_thin_zone = offset(next_onion, (float)(ext_perimeter_width / 2));
-                        ExPolygons expp = offset2_ex(
-                            // medial axis requires non-overlapping geometry
-                            diff_ex(to_polygons(last),
-                                    no_thin_zone,
-                                    true),
-                                    (float)(-min_width / 2), (float)(min_width / 2));
+                        ExPolygons no_thin_zone = offset_ex(next_onion, (float)(ext_perimeter_width / 2));
+                        // medial axis requires non-overlapping geometry
+                        ExPolygons thin_zones = diff_ex(last, no_thin_zone, true);
+                        //don't use offset2_ex, because we don't want to merge the zones that have been separated.
+                        ExPolygons expp = offset_ex(thin_zones, (float)(-min_width / 2));
+                        //we push the bits removed and put them into what we will use as our anchor
+                        if (expp.size() > 0) {
+                            no_thin_zone = diff_ex(last, offset_ex(expp, (float)(min_width / 2)), true);
+                        }
                         // compute a bit of overlap to anchor thin walls inside the print.
-                        ExPolygons anchor = intersection_ex(to_polygons(offset_ex(expp, (float)(ext_perimeter_width / 2))), no_thin_zone, true);
                         for (ExPolygon &ex : expp) {
-                            ExPolygons bounds = union_ex(ExPolygons() = { ex }, anchor, true);
+                            //growing back the polygon
+                            //a vary little bit of overlap can be created here with other thin polygon, but it's more useful than worisome.
+                            ExPolygons ex_bigger = offset_ex(ex, (float)(min_width / 2));
+                            if (ex_bigger.size() != 1) continue; // impossible error, growing a single polygon can't create multiple or 0.
+                            ExPolygons anchor = intersection_ex(offset_ex(ex, (float)(min_width / 2) + 
+                                (float)(ext_perimeter_width / 2), jtSquare), no_thin_zone, true);
+                            ExPolygons bounds = union_ex(ex_bigger, 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_perimeter_width + ext_perimeter_spacing2, min_width, &thin_walls);
-                                    continue;
+                                if (!intersection_ex(ex_bigger[0], bound).empty()) {
+                                    //be sure it's not too small to extrude reliably
+                                    if (ex_bigger[0].area() > min_width*(ext_perimeter_width + ext_perimeter_spacing2)) {
+                                        // the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop
+                                        ex_bigger[0].medial_axis(bound, ext_perimeter_width + ext_perimeter_spacing2, min_width, 
+                                            &thin_walls, this->layer_height);
+                                    }
+                                    break;
                                 }
                             }
                         }
@@ -390,13 +419,17 @@ void PerimeterGenerator::process()
             double min = 0.2 * perimeter_width * (1 - INSET_OVERLAP_TOLERANCE);
             double max = 2. * perimeter_spacing;
             ExPolygons gaps_ex = diff_ex(
-                //FIXME offset2 would be enough and cheaper.
                 offset2_ex(gaps, -min/2, +min/2),
                 offset2_ex(gaps, -max/2, +max/2),
                 true);
             ThickPolylines polylines;
-            for (const ExPolygon &ex : gaps_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) {
+                    ex.medial_axis(ex, max, min, &polylines, this->layer_height);
+                }
+            }
             if (!polylines.empty()) {
                 ExtrusionEntityCollection gap_fill = this->_variable_width(polylines, 
                     erGapFill, this->solid_infill_flow);
diff --git a/xs/src/libslic3r/Polyline.cpp b/xs/src/libslic3r/Polyline.cpp
index 3432506c6..f015959ff 100644
--- a/xs/src/libslic3r/Polyline.cpp
+++ b/xs/src/libslic3r/Polyline.cpp
@@ -6,6 +6,7 @@
 #include "Polygon.hpp"
 #include <iostream>
 #include <utility>
+#include <algorithm>
 
 namespace Slic3r {
 
@@ -262,8 +263,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);
         }
     }
@@ -292,4 +293,96 @@ Lines3 Polyline3::lines() const
     return lines;
 }
 
+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];
+
+            int32_t id_candidate_first_point = -1;
+            int32_t id_candidate_last_point = -1;
+            int32_t nbCandidate_first_point = 0;
+            int32_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 (polyline->last_point().coincides_with(other->last_point())) {
+                    other->reverse();
+                    id_candidate_last_point = j;
+                    nbCandidate_last_point++;
+                } else if (polyline->first_point().coincides_with(other->last_point())) {
+                    id_candidate_first_point = j;
+                    nbCandidate_first_point++;
+                } else if (polyline->first_point().coincides_with(other->first_point())) {
+                    id_candidate_first_point = j;
+                    nbCandidate_first_point++;
+                    other->reverse();
+                } else if (polyline->last_point().coincides_with(other->first_point())) {
+                    id_candidate_last_point = j;
+                    nbCandidate_last_point++;
+                } else {
+                    continue;
+                }
+            }
+            if (id_candidate_last_point == id_candidate_first_point && nbCandidate_first_point == 1 && nbCandidate_last_point == 1) {
+                // it's a trap! it's a  loop!
+                if (pp[id_candidate_first_point].points.size() > 2) {
+                    polyline->points.insert(polyline->points.begin(), pp[id_candidate_first_point].points.begin() + 1, pp[id_candidate_first_point].points.end() - 1);
+                    polyline->width.insert(polyline->width.begin(), pp[id_candidate_first_point].width.begin() + 1, pp[id_candidate_first_point].width.end() - 1);
+                }
+                pp.erase(pp.begin() + id_candidate_first_point);
+                changes = true;
+                polyline->endpoints.first = false;
+                polyline->endpoints.second = false;
+                continue;
+            }
+
+            if (nbCandidate_first_point == 1) {
+                //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 && !polyline->endpoints.first && !polyline->first_point().coincides_with(polyline->last_point())) {
+                //update endpoint
+                polyline->endpoints.first = true;
+            }
+            if (nbCandidate_last_point == 1) {
+                //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 && !polyline->endpoints.second && !polyline->first_point().coincides_with(polyline->last_point())) {
+                //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 b64743d84..860c49a7a 100644
--- a/xs/src/libslic3r/Polyline.hpp
+++ b/xs/src/libslic3r/Polyline.hpp
@@ -128,15 +128,26 @@ inline void polylines_append(Polylines &dst, Polylines &&src)
 
 bool remove_degenerate(Polylines &polylines);
 
+
+/// ThickPolyline : a polyline with a width for each point
+/// This calss has a vector of coordf_t, it must be the same size than 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);
+
 class Polyline3 : public MultiPoint3
 {
 public: