Merge remote-tracking branch 'remote/master' into SoftFever

src/libslic3r/AABBTreeIndirect.hpp

@@ -446,6 +446,57 @@ namespace detail {
 		}
 	}
 
+    // Real-time collision detection, Ericson, Chapter 5
+    template<typename Vector>
+    static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c)
+    {
+        using Scalar = typename Vector::Scalar;
+        // Check if P in vertex region outside A
+        Vector ab = b - a;
+        Vector ac = c - a;
+        Vector ap = p - a;
+        Scalar d1 = ab.dot(ap);
+        Scalar d2 = ac.dot(ap);
+        if (d1 <= 0 && d2 <= 0)
+          return a;
+        // Check if P in vertex region outside B
+        Vector bp = p - b;
+        Scalar d3 = ab.dot(bp);
+        Scalar d4 = ac.dot(bp);
+        if (d3 >= 0 && d4 <= d3)
+          return b;
+        // Check if P in edge region of AB, if so return projection of P onto AB
+        Scalar vc = d1*d4 - d3*d2;
+        if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) {
+            Scalar v = d1 / (d1 - d3);
+            return a + v * ab;
+        }
+        // Check if P in vertex region outside C
+        Vector cp = p - c;
+        Scalar d5 = ab.dot(cp);
+        Scalar d6 = ac.dot(cp);
+        if (d6 >= 0 && d5 <= d6)
+          return c;
+        // Check if P in edge region of AC, if so return projection of P onto AC
+        Scalar vb = d5*d2 - d1*d6;
+        if (vb <= 0 && d2 >= 0 && d6 <= 0) {
+          Scalar w = d2 / (d2 - d6);
+          return a + w * ac;
+        }
+        // Check if P in edge region of BC, if so return projection of P onto BC
+        Scalar va = d3*d6 - d5*d4;
+        if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
+          Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
+          return b + w * (c - b);
+        }
+        // P inside face region. Compute Q through its barycentric coordinates (u,v,w)
+        Scalar denom = Scalar(1.0) / (va + vb + vc);
+        Scalar v = vb * denom;
+        Scalar w = vc * denom;
+        return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
+    };
+
+
 	// Nothing to do with COVID-19 social distancing.
 	template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
 	struct IndexedTriangleSetDistancer {
@@ -453,74 +504,36 @@ namespace detail {
 		using IndexedFaceType 	= AIndexedFaceType;
 		using TreeType			= ATreeType;
 		using VectorType 		= AVectorType;
+		using ScalarType 		= typename VectorType::Scalar;
 
 		const std::vector<VertexType> 		&vertices;
 		const std::vector<IndexedFaceType> 	&faces;
 		const TreeType 						&tree;
 
 		const VectorType					 origin;
+
+		inline VectorType closest_point_to_origin(size_t primitive_index,
+		        ScalarType& squared_distance){
+		    const auto &triangle = this->faces[primitive_index];
+		    VectorType closest_point = closest_point_to_triangle<VectorType>(origin,
+		            this->vertices[triangle(0)].template cast<ScalarType>(),
+		            this->vertices[triangle(1)].template cast<ScalarType>(),
+		            this->vertices[triangle(2)].template cast<ScalarType>());
+		    squared_distance = (origin - closest_point).squaredNorm();
+		    return closest_point;
+		}
 	};
 
-	// Real-time collision detection, Ericson, Chapter 5
-	template<typename Vector>
-	static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c)
-	{
-		using Scalar = typename Vector::Scalar;
-		// Check if P in vertex region outside A
-		Vector ab = b - a;
-		Vector ac = c - a;
-		Vector ap = p - a;
-		Scalar d1 = ab.dot(ap);
-		Scalar d2 = ac.dot(ap);
-		if (d1 <= 0 && d2 <= 0)
-		  return a;
-		// Check if P in vertex region outside B
-		Vector bp = p - b;
-		Scalar d3 = ab.dot(bp);
-		Scalar d4 = ac.dot(bp);
-		if (d3 >= 0 && d4 <= d3)
-		  return b;
-		// Check if P in edge region of AB, if so return projection of P onto AB
-		Scalar vc = d1*d4 - d3*d2;
-		if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) {
-		    Scalar v = d1 / (d1 - d3);
-		    return a + v * ab;
-		}
-		// Check if P in vertex region outside C
-		Vector cp = p - c;
-		Scalar d5 = ab.dot(cp);
-		Scalar d6 = ac.dot(cp);
-		if (d6 >= 0 && d5 <= d6)
-		  return c;
-		// Check if P in edge region of AC, if so return projection of P onto AC
-		Scalar vb = d5*d2 - d1*d6;
-		if (vb <= 0 && d2 >= 0 && d6 <= 0) {
-		  Scalar w = d2 / (d2 - d6);
-		  return a + w * ac;
-		}
-		// Check if P in edge region of BC, if so return projection of P onto BC
-		Scalar va = d3*d6 - d5*d4;
-		if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
-		  Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
-		  return b + w * (c - b);
-		}
-		// P inside face region. Compute Q through its barycentric coordinates (u,v,w)
-		Scalar denom = Scalar(1.0) / (va + vb + vc);
-		Scalar v = vb * denom;
-		Scalar w = vc * denom;
-		return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
-	};
-
-	template<typename IndexedTriangleSetDistancerType, typename Scalar>
-    static inline Scalar squared_distance_to_indexed_triangle_set_recursive(
-        IndexedTriangleSetDistancerType	&distancer,
+	template<typename IndexedPrimitivesDistancerType, typename Scalar>
+    static inline Scalar squared_distance_to_indexed_primitives_recursive(
+        IndexedPrimitivesDistancerType	&distancer,
 		size_t 							 node_idx,
 		Scalar 							 low_sqr_d,
   		Scalar 							 up_sqr_d,
 		size_t 							&i,
-  		Eigen::PlainObjectBase<typename IndexedTriangleSetDistancerType::VectorType> &c)
+  		Eigen::PlainObjectBase<typename IndexedPrimitivesDistancerType::VectorType> &c)
 	{
-		using Vector = typename IndexedTriangleSetDistancerType::VectorType;
+		using Vector = typename IndexedPrimitivesDistancerType::VectorType;
 
   		if (low_sqr_d > up_sqr_d)
 			return low_sqr_d;
@@ -538,13 +551,9 @@ namespace detail {
 		assert(node.is_valid());
   		if (node.is_leaf()) 
   		{
-            const auto &triangle = distancer.faces[node.idx];
-            Vector c_candidate = closest_point_to_triangle<Vector>(
-				distancer.origin, 
-                distancer.vertices[triangle(0)].template cast<Scalar>(),
-                distancer.vertices[triangle(1)].template cast<Scalar>(),
-                distancer.vertices[triangle(2)].template cast<Scalar>());
-            set_min((c_candidate - distancer.origin).squaredNorm(), node.idx, c_candidate);
+            Scalar sqr_dist;
+            Vector c_candidate = distancer.closest_point_to_origin(node.idx, sqr_dist);
+            set_min(sqr_dist, node.idx, c_candidate);
   		} 
   		else
   		{
@@ -561,7 +570,7 @@ namespace detail {
 			{
                 size_t	i_left;
                 Vector 	c_left = c;
-                Scalar	sqr_d_left = squared_distance_to_indexed_triangle_set_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left);
+                Scalar	sqr_d_left = squared_distance_to_indexed_primitives_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left);
 				set_min(sqr_d_left, i_left, c_left);
 				looked_left = true;
 			};
@@ -569,13 +578,13 @@ namespace detail {
 			{
                 size_t	i_right;
                 Vector	c_right = c;
-                Scalar	sqr_d_right = squared_distance_to_indexed_triangle_set_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right);
+                Scalar	sqr_d_right = squared_distance_to_indexed_primitives_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right);
 				set_min(sqr_d_right, i_right, c_right);
 				looked_right = true;
 			};
 
 			// must look left or right if in box
-            using BBoxScalar = typename IndexedTriangleSetDistancerType::TreeType::BoundingBox::Scalar;
+            using BBoxScalar = typename IndexedPrimitivesDistancerType::TreeType::BoundingBox::Scalar;
             if (node_left.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
 			  	look_left();
             if (node_right.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
@@ -709,10 +718,15 @@ inline bool intersect_ray_all_hits(
         origin, dir, VectorType(dir.cwiseInverse()),
         eps }
 	};
-	if (! tree.empty()) {
+	if (tree.empty()) {
+		hits.clear();
+	} else {
+		// Reusing the output memory if there is some memory already pre-allocated.
+        ray_intersector.hits = std::move(hits);
+        ray_intersector.hits.clear();
         ray_intersector.hits.reserve(8);
 		detail::intersect_ray_recursive_all_hits(ray_intersector, 0);
-		std::swap(hits, ray_intersector.hits);
+		hits = std::move(ray_intersector.hits);
 	    std::sort(hits.begin(), hits.end(), [](const auto &l, const auto &r) { return l.t < r.t; });
 	}
 	return ! hits.empty();
@@ -742,7 +756,7 @@ inline typename VectorType::Scalar squared_distance_to_indexed_triangle_set(
     auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
         { vertices, faces, tree, point };
     return tree.empty() ? Scalar(-1) : 
-    	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
+    	detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
 }
 
 // Decides if exists some triangle in defined radius on a 3D indexed triangle set using a pre-built AABBTreeIndirect::Tree.
@@ -759,22 +773,22 @@ inline bool is_any_triangle_in_radius(
         const TreeType 						&tree,
         // Point to which the closest point on the indexed triangle set is searched for.
         const VectorType					&point,
-        // Maximum distance in which triangle is search for
-        typename VectorType::Scalar &max_distance)
+        //Square of maximum distance in which triangle is searched for
+        typename VectorType::Scalar &max_distance_squared)
 {
     using Scalar = typename VectorType::Scalar;
     auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
             { vertices, faces, tree, point };
 
-	size_t hit_idx;
-	VectorType hit_point = VectorType::Ones() * (std::nan(""));
+    size_t hit_idx;
+    VectorType hit_point = VectorType::Ones() * (NaN<typename VectorType::Scalar>);
 
 	if(tree.empty())
 	{
 		return false;
 	}
 
-	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance, hit_idx, hit_point);
+	detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
 
     return hit_point.allFinite();
 }

src/libslic3r/AABBTreeLines.hpp

@@ -0,0 +1,112 @@
+#ifndef SRC_LIBSLIC3R_AABBTREELINES_HPP_
+#define SRC_LIBSLIC3R_AABBTREELINES_HPP_
+
+#include "libslic3r/Point.hpp"
+#include "libslic3r/EdgeGrid.hpp"
+#include "libslic3r/AABBTreeIndirect.hpp"
+#include "libslic3r/Line.hpp"
+
+namespace Slic3r {
+
+namespace AABBTreeLines {
+
+namespace detail {
+
+template<typename ALineType, typename ATreeType, typename AVectorType>
+struct IndexedLinesDistancer {
+    using LineType = ALineType;
+    using TreeType = ATreeType;
+    using VectorType = AVectorType;
+    using ScalarType = typename VectorType::Scalar;
+
+    const std::vector<LineType> &lines;
+    const TreeType &tree;
+
+    const VectorType origin;
+
+    inline VectorType closest_point_to_origin(size_t primitive_index,
+            ScalarType &squared_distance) {
+        VectorType nearest_point;
+        const LineType &line = lines[primitive_index];
+        squared_distance = line_alg::distance_to_squared(line, origin, &nearest_point);
+        return nearest_point;
+    }
+};
+
+}
+
+// Build a balanced AABB Tree over a vector of float lines, balancing the tree
+// on centroids of the lines.
+// Epsilon is applied to the bounding boxes of the AABB Tree to cope with numeric inaccuracies
+// during tree traversal.
+template<typename LineType>
+inline AABBTreeIndirect::Tree<2, typename LineType::Scalar> build_aabb_tree_over_indexed_lines(
+        const std::vector<LineType> &lines,
+        //FIXME do we want to apply an epsilon?
+        const float eps = 0)
+        {
+    using TreeType = AABBTreeIndirect::Tree<2, typename LineType::Scalar>;
+//    using              CoordType      = typename TreeType::CoordType;
+    using VectorType = typename TreeType::VectorType;
+    using BoundingBox = typename TreeType::BoundingBox;
+
+    struct InputType {
+        size_t idx() const {
+            return m_idx;
+        }
+        const BoundingBox& bbox() const {
+            return m_bbox;
+        }
+        const VectorType& centroid() const {
+            return m_centroid;
+        }
+
+        size_t m_idx;
+        BoundingBox m_bbox;
+        VectorType m_centroid;
+    };
+
+    std::vector<InputType> input;
+    input.reserve(lines.size());
+    const VectorType veps(eps, eps);
+    for (size_t i = 0; i < lines.size(); ++i) {
+        const LineType &line = lines[i];
+        InputType n;
+        n.m_idx = i;
+        n.m_centroid = (line.a + line.b) * 0.5;
+        n.m_bbox = BoundingBox(line.a, line.a);
+        n.m_bbox.extend(line.b);
+        n.m_bbox.min() -= veps;
+        n.m_bbox.max() += veps;
+        input.emplace_back(n);
+    }
+
+    TreeType out;
+    out.build(std::move(input));
+    return out;
+}
+
+// Finding a closest line, its closest point and squared distance to the closest point
+// Returns squared distance to the closest point or -1 if the input is empty.
+template<typename LineType, typename TreeType, typename VectorType>
+inline typename VectorType::Scalar squared_distance_to_indexed_lines(
+        const std::vector<LineType> &lines,
+        const TreeType &tree,
+        const VectorType &point,
+        size_t &hit_idx_out,
+        Eigen::PlainObjectBase<VectorType> &hit_point_out)
+        {
+    using Scalar = typename VectorType::Scalar;
+    auto distancer = detail::IndexedLinesDistancer<LineType, TreeType, VectorType>
+            { lines, tree, point };
+    return tree.empty() ?
+                          Scalar(-1) :
+                          AABBTreeIndirect::detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0),
+                                  std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
+}
+
+}
+
+}
+
+#endif /* SRC_LIBSLIC3R_AABBTREELINES_HPP_ */

src/libslic3r/AppConfig.cpp

@@ -184,10 +184,10 @@ void AppConfig::set_defaults()
 
 #ifdef _WIN32
 
-#ifdef SUPPORT_3D_CONNEXION
+//#ifdef SUPPORT_3D_CONNEXION
     if (get("use_legacy_3DConnexion").empty())
-        set_bool("use_legacy_3DConnexion", false);
-#endif
+        set_bool("use_legacy_3DConnexion", true);
+//#endif
 
 #ifdef SUPPORT_DARK_MODE
     if (get("dark_color_mode").empty())

src/libslic3r/Arrange.cpp

@@ -447,7 +447,7 @@ protected:
         return std::make_tuple(score, fullbb);
     }
     
-    std::function<double(const Item&)> get_objfn();
+    std::function<double(const Item&, const ItemGroup&)> get_objfn();
     
 public:
     AutoArranger(const TBin &                  bin,
@@ -508,7 +508,8 @@ public:
             bin_poly.contour.points.emplace_back(c0.x(), c1.y());
             return bin_poly;
         };
-
+        
+        // preload fixed items (and excluded regions) on plate
         m_pconf.on_preload = [this](const ItemGroup &items, PConfig &cfg) {
             if (items.empty()) return;
 
@@ -527,8 +528,12 @@ public:
                 }
             }
 
-            cfg.object_function = [this, bb, starting_point](const Item& item) {
-                return fixed_overfit_topright_sliding(objfunc(item, starting_point), bb);
+            cfg.object_function = [this, bb, starting_point](const Item& item, const ItemGroup& packed_items) {
+                bool packed_are_excluded_region = std::all_of(packed_items.begin(), packed_items.end(), [](Item& itm) { return itm.is_virt_object && !itm.is_wipe_tower; });
+                if(packed_are_excluded_region)
+                    return fixed_overfit_topright_sliding(objfunc(item, starting_point), bb);
+                else
+                    return fixed_overfit(objfunc(item, starting_point), bb);
             };
         };
 
@@ -597,11 +602,11 @@ public:
     }
 };
 
-template<> std::function<double(const Item&)> AutoArranger<Box>::get_objfn()
+template<> std::function<double(const Item&, const ItemGroup&)> AutoArranger<Box>::get_objfn()
 {
     auto origin_pack = m_pconf.starting_point == PConfig::Alignment::CENTER ? m_bin.center() : m_bin.minCorner();
 
-    return [this, origin_pack](const Item &itm) {
+    return [this, origin_pack](const Item &itm, const ItemGroup&) {
         auto result = objfunc(itm, origin_pack);
         
         double score = std::get<0>(result);
@@ -623,11 +628,11 @@ template<> std::function<double(const Item&)> AutoArranger<Box>::get_objfn()
     };
 }
 
-template<> std::function<double(const Item&)> AutoArranger<Circle>::get_objfn()
+template<> std::function<double(const Item&, const ItemGroup&)> AutoArranger<Circle>::get_objfn()
 {
     auto bb = sl::boundingBox(m_bin);
     auto origin_pack = m_pconf.starting_point == PConfig::Alignment::CENTER ? bb.center() : bb.minCorner();
-    return [this, origin_pack](const Item &item) {
+    return [this, origin_pack](const Item &item, const ItemGroup&) {
         
         auto result = objfunc(item, origin_pack);
         
@@ -653,11 +658,11 @@ template<> std::function<double(const Item&)> AutoArranger<Circle>::get_objfn()
 // Specialization for a generalized polygon.
 // Warning: this is much slower than with Box bed. Need further speedup.
 template<>
-std::function<double(const Item &)> AutoArranger<ExPolygon>::get_objfn()
+std::function<double(const Item &, const ItemGroup&)> AutoArranger<ExPolygon>::get_objfn()
 {
     auto bb = sl::boundingBox(m_bin);
     auto origin_pack = m_pconf.starting_point == PConfig::Alignment::CENTER ? bb.center() : bb.minCorner();
-    return [this, origin_pack](const Item &itm) {
+    return [this, origin_pack](const Item &itm, const ItemGroup&) {
         auto result = objfunc(itm, origin_pack);
 
         double score = std::get<0>(result);

src/libslic3r/Brim.cpp

@@ -270,7 +270,7 @@ static ExPolygons top_level_outer_brim_area(const Print                   &print
     return diff_ex(brim_area, no_brim_area);
 }
 
-// BBS: the brims of different objs will not overlapped with each other, and are stored by objs and by extruders 
+// BBS: the brims of different objs will not overlapped with each other, and are stored by objs and by extruders
 static ExPolygons top_level_outer_brim_area(const Print& print, const ConstPrintObjectPtrs& top_level_objects_with_brim,
     const float no_brim_offset, double& brim_width_max, std::map<ObjectID, double>& brim_width_map,
     std::map<ObjectID, ExPolygons>& brimAreaMap,
@@ -453,7 +453,7 @@ static ExPolygons inner_brim_area(const Print                   &print,
     return diff_ex(intersection_ex(to_polygons(std::move(brim_area)), holes), no_brim_area);
 }
 
-// BBS: the brims of different objs will not overlapped with each other, and are stored by objs and by extruders 
+// BBS: the brims of different objs will not overlapped with each other, and are stored by objs and by extruders
 static ExPolygons inner_brim_area(const Print& print, const ConstPrintObjectPtrs& top_level_objects_with_brim,
     const float no_brim_offset, std::map<ObjectID, ExPolygons>& brimAreaMap,
     std::map<ObjectID, ExPolygons>& supportBrimAreaMap,
@@ -875,7 +875,7 @@ static ExPolygons outer_inner_brim_area(const Print& print,
     for (const auto& objectWithExtruder : objPrintVec)
         brimToWrite.insert({ objectWithExtruder.first, {true,true} });
 
-    std::map<ObjectID, ExPolygons> objectIslandMap;
+    ExPolygons objectIslands;
 
     for (unsigned int extruderNo : printExtruders) {
         ++extruderNo;
@@ -906,7 +906,11 @@ static ExPolygons outer_inner_brim_area(const Print& print,
                     std::vector<ModelVolume*> groupVolumePtrs;
                     for (auto& volumeID : volumeGroup.volume_ids) {
                         ModelVolume* currentModelVolumePtr = nullptr;
-                        for (auto volumePtr : object->model_object()->volumes) {
+                        //BBS: support shared object logic
+                        const PrintObject* shared_object = object->get_shared_object();
+                        if (!shared_object)
+                            shared_object = object;
+                        for (auto volumePtr : shared_object->model_object()->volumes) {
                             if (volumePtr->id() == volumeID) {
                                 currentModelVolumePtr = volumePtr;
                                 break;
@@ -968,7 +972,7 @@ static ExPolygons outer_inner_brim_area(const Print& print,
                         append_and_translate(brim_area, brim_area_object, instance, print, brimAreaMap);
                     append_and_translate(no_brim_area, no_brim_area_object, instance);
                     append_and_translate(holes, holes_object, instance);
-                    append_and_translate(objectIslandMap[instance.print_object->id()], objectIsland, instance);
+                    append_and_translate(objectIslands, objectIsland, instance);
 
                 }
                 if (brimAreaMap.find(object->id()) != brimAreaMap.end())
@@ -1032,27 +1036,44 @@ static ExPolygons outer_inner_brim_area(const Print& print,
             }
         }
     }
-    for (const PrintObject* object : print.objects()) {
+    for (const PrintObject* object : print.objects()) 
         if (brimAreaMap.find(object->id()) != brimAreaMap.end()) {
             brimAreaMap[object->id()] = diff_ex(brimAreaMap[object->id()], no_brim_area);
 
-            // BBS: brim should be contacted to at least one object island
-            if (objectIslandMap.find(object->id()) != objectIslandMap.end() && !objectIslandMap[object->id()].empty()) {
-                auto tempArea = brimAreaMap[object->id()];
-                brimAreaMap[object->id()].clear();
-                // the error bound is set to 2x flow width
-                for (auto& ta : tempArea) {
-                    auto offsetedTa = offset_ex(ta, print.brim_flow().scaled_spacing() * 2, jtRound, SCALED_RESOLUTION);
-                    if (!intersection_ex(offsetedTa, objectIslandMap[object->id()]).empty())
-                        brimAreaMap[object->id()].push_back(ta);
-                }
-            }
-        }
-
         if (supportBrimAreaMap.find(object->id()) != supportBrimAreaMap.end())
             supportBrimAreaMap[object->id()] = diff_ex(supportBrimAreaMap[object->id()], no_brim_area);
     }
-    //brim_area = diff_ex(brim_area, no_brim_area);
+
+    brim_area.clear();
+    for (const PrintObject* object : print.objects()) {
+        // BBS: brim should be contacted to at least one object's island or brim area
+        if (brimAreaMap.find(object->id()) != brimAreaMap.end()) {
+            // find other objects' brim area
+            ExPolygons otherExPolys;
+            for (const PrintObject* otherObject : print.objects()) {
+                if ((otherObject->id() != object->id()) && (brimAreaMap.find(otherObject->id()) != brimAreaMap.end())) {
+                    expolygons_append(otherExPolys, brimAreaMap[otherObject->id()]);
+                }
+            }
+
+            auto tempArea = brimAreaMap[object->id()];
+            brimAreaMap[object->id()].clear();
+
+            for (int ia = 0; ia != tempArea.size(); ++ia) {
+                // find this object's other brim area
+                ExPolygons otherExPoly;
+                for (int iao = 0; iao != tempArea.size(); ++iao)
+                    if (iao != ia) otherExPoly.push_back(tempArea[iao]);
+
+                auto offsetedTa = offset_ex(tempArea[ia], print.brim_flow().scaled_spacing() * 2, jtRound, SCALED_RESOLUTION);
+                if (!intersection_ex(offsetedTa, objectIslands).empty() ||
+                    !intersection_ex(offsetedTa, otherExPoly).empty() ||
+                    !intersection_ex(offsetedTa, otherExPolys).empty())
+                    brimAreaMap[object->id()].push_back(tempArea[ia]);
+            }
+            expolygons_append(brim_area, brimAreaMap[object->id()]);
+        }
+    }
     return brim_area;
 }
 // Flip orientation of open polylines to minimize travel distance.
@@ -1066,7 +1087,7 @@ static void optimize_polylines_by_reversing(Polylines *polylines)
             double dist_to_start = (next.first_point() - prev.last_point()).cast<double>().norm();
             double dist_to_end   = (next.last_point() - prev.last_point()).cast<double>().norm();
 
-            if (dist_to_end < dist_to_start) 
+            if (dist_to_end < dist_to_start)
                 next.reverse();
         }
     }

src/libslic3r/CMakeLists.txt

@@ -142,10 +142,12 @@ set(lisbslic3r_sources
     GCodeWriter.hpp
     Geometry.cpp
     Geometry.hpp
+    Geometry/Bicubic.hpp
     Geometry/Circle.cpp
     Geometry/Circle.hpp
     Geometry/ConvexHull.cpp
     Geometry/ConvexHull.hpp
+    Geometry/Curves.hpp
     Geometry/MedialAxis.cpp
     Geometry/MedialAxis.hpp
     Geometry/Voronoi.hpp
@@ -176,6 +178,8 @@ set(lisbslic3r_sources
     CustomGCode.hpp
     Arrange.hpp
     Arrange.cpp
+    NormalUtils.cpp
+    NormalUtils.hpp
     Orient.hpp
     Orient.cpp
     MultiPoint.cpp
@@ -222,6 +226,8 @@ set(lisbslic3r_sources
     QuadricEdgeCollapse.cpp
     QuadricEdgeCollapse.hpp
     Semver.cpp
+    ShortEdgeCollapse.cpp
+    ShortEdgeCollapse.hpp
     ShortestPath.cpp
     ShortestPath.hpp
     SLAPrint.cpp
@@ -264,6 +270,8 @@ set(lisbslic3r_sources
     Thread.hpp
     TriangleSelector.cpp
     TriangleSelector.hpp
+    TriangleSetSampling.cpp
+    TriangleSetSampling.hpp
     MTUtils.hpp
     VariableWidth.cpp
     VariableWidth.hpp
@@ -315,7 +323,6 @@ set(lisbslic3r_sources
     SLA/Clustering.hpp
     SLA/Clustering.cpp
     SLA/ReprojectPointsOnMesh.hpp
-
     Arachne/BeadingStrategy/BeadingStrategy.hpp
     Arachne/BeadingStrategy/BeadingStrategy.cpp
     Arachne/BeadingStrategy/BeadingStrategyFactory.hpp
@@ -356,6 +363,8 @@ set(lisbslic3r_sources
     Arachne/SkeletalTrapezoidationJoint.hpp
     Arachne/WallToolPaths.hpp
     Arachne/WallToolPaths.cpp
+    Shape/TextShape.hpp
+    Shape/TextShape.cpp
 )
 
 if (APPLE)

src/libslic3r/ExtrusionEntity.cpp

@@ -158,11 +158,10 @@ double ExtrusionLoop::length() const
     return len;
 }
 
-bool ExtrusionLoop::split_at_vertex(const Point &point)
+bool ExtrusionLoop::split_at_vertex(const Point &point, const double scaled_epsilon)
 {
     for (ExtrusionPaths::iterator path = this->paths.begin(); path != this->paths.end(); ++path) {
-        int idx = path->polyline.find_point(point);
-        if (idx != -1) {
+        if (int idx = path->polyline.find_point(point, scaled_epsilon); idx != -1) {
             if (this->paths.size() == 1) {
                 // just change the order of points
                 Polyline p1, p2;
@@ -207,46 +206,57 @@ bool ExtrusionLoop::split_at_vertex(const Point &point)
     return false;
 }
 
-std::pair<size_t, Point> ExtrusionLoop::get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const
+ExtrusionLoop::ClosestPathPoint ExtrusionLoop::get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const
 {
     // Find the closest path and closest point belonging to that path. Avoid overhangs, if asked for.
-    size_t path_idx = 0;
-    Point  p;
-    {
-        double min = std::numeric_limits<double>::max();
-        Point  p_non_overhang;
-        size_t path_idx_non_overhang = 0;
-        double min_non_overhang = std::numeric_limits<double>::max();
-        for (const ExtrusionPath& path : this->paths) {
-            Point p_tmp = point.projection_onto(path.polyline);
-            double dist = (p_tmp - point).cast<double>().norm();
-            if (dist < min) {
-                p = p_tmp;
-                min = dist;
-                path_idx = &path - &this->paths.front();
-            }
-            if (prefer_non_overhang && !is_bridge(path.role()) && dist < min_non_overhang) {
-                p_non_overhang = p_tmp;
-                min_non_overhang = dist;
-                path_idx_non_overhang = &path - &this->paths.front();
-            }
+    ClosestPathPoint out{0, 0};
+    double           min2 = std::numeric_limits<double>::max();
+    ClosestPathPoint best_non_overhang{0, 0};
+    double           min2_non_overhang = std::numeric_limits<double>::max();
+    for (const ExtrusionPath &path : this->paths) {
+        std::pair<int, Point> foot_pt_ = foot_pt(path.polyline.points, point);
+        double                d2       = (foot_pt_.second - point).cast<double>().squaredNorm();
+        if (d2 < min2) {
+            out.foot_pt     = foot_pt_.second;
+            out.path_idx    = &path - &this->paths.front();
+            out.segment_idx = foot_pt_.first;
+            min2            = d2;
         }
-        if (prefer_non_overhang && min_non_overhang != std::numeric_limits<double>::max()) {
-            // Only apply the non-overhang point if there is one.
-            path_idx = path_idx_non_overhang;
-            p = p_non_overhang;
+        if (prefer_non_overhang && !is_bridge(path.role()) && d2 < min2_non_overhang) {
+            best_non_overhang.foot_pt     = foot_pt_.second;
+            best_non_overhang.path_idx    = &path - &this->paths.front();
+            best_non_overhang.segment_idx = foot_pt_.first;
+            min2_non_overhang             = d2;
         }
     }
-    return std::make_pair(path_idx, p);
+    if (prefer_non_overhang && min2_non_overhang != std::numeric_limits<double>::max())
+        // Only apply the non-overhang point if there is one.
+        out = best_non_overhang;
+    return out;
 }
 
 // Splitting an extrusion loop, possibly made of multiple segments, some of the segments may be bridging.
-void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang)
+void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon)
 {
     if (this->paths.empty())
         return;
     
-    auto [path_idx, p] = get_closest_path_and_point(point, prefer_non_overhang);
+    auto [path_idx, segment_idx, p] = get_closest_path_and_point(point, prefer_non_overhang);
+
+    // Snap p to start or end of segment_idx if closer than scaled_epsilon.
+    {
+        const Point *p1 = this->paths[path_idx].polyline.points.data() + segment_idx;
+        const Point *p2 = p1;
+        ++p2;
+        double       d2_1 = (point - *p1).cast<double>().squaredNorm();
+        double       d2_2 = (point - *p2).cast<double>().squaredNorm();
+        const double thr2 = scaled_epsilon * scaled_epsilon;
+        if (d2_1 < d2_2) {
+            if (d2_1 < thr2) p = *p1;
+        } else {
+            if (d2_2 < thr2) p = *p2;
+        }
+    }
     
     // now split path_idx in two parts
     const ExtrusionPath &path = this->paths[path_idx];

src/libslic3r/ExtrusionEntity.hpp

@@ -314,9 +314,15 @@ public:
     const Point& last_point() const override { assert(this->first_point() == this->paths.back().polyline.points.back()); return this->first_point(); }
     Polygon polygon() const;
     double length() const override;
-    bool split_at_vertex(const Point &point);
-    void split_at(const Point &point, bool prefer_non_overhang);
-    std::pair<size_t, Point> get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const;
+    bool split_at_vertex(const Point &point, const double scaled_epsilon = scaled<double>(0.001));
+    void split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon = scaled<double>(0.001));
+    struct ClosestPathPoint
+    {
+        size_t path_idx;
+        size_t segment_idx;
+        Point  foot_pt;
+    };
+    ClosestPathPoint         get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const;
     void clip_end(double distance, ExtrusionPaths* paths) const;
     // Test, whether the point is extruded by a bridging flow.
     // This used to be used to avoid placing seams on overhangs, but now the EdgeGrid is used instead.

src/libslic3r/Format/STEP.cpp

@@ -31,9 +31,6 @@
 const double STEP_TRANS_CHORD_ERROR = 0.005;
 const double STEP_TRANS_ANGLE_RES = 1;
 
-const int LOAD_STEP_STAGE_READ_FILE          = 0;
-const int LOAD_STEP_STAGE_GET_SOLID          = 1;
-const int LOAD_STEP_STAGE_GET_MESH           = 2;
 
 namespace Slic3r {
 
@@ -213,11 +210,11 @@ static void getNamedSolids(const TopLoc_Location& location, const std::string& p
     }
 }
 
-bool load_step(const char *path, Model *model, ImportStepProgressFn proFn, StepIsUtf8Fn isUtf8Fn)
+bool load_step(const char *path, Model *model, ImportStepProgressFn stepFn, StepIsUtf8Fn isUtf8Fn)
 {
     bool cb_cancel = false;
-    if (proFn) {
-        proFn(LOAD_STEP_STAGE_READ_FILE, 0, 1, cb_cancel);
+    if (stepFn) {
+        stepFn(LOAD_STEP_STAGE_READ_FILE, 0, 1, cb_cancel);
         if (cb_cancel)
             return false;
     }
@@ -245,9 +242,13 @@ bool load_step(const char *path, Model *model, ImportStepProgressFn proFn, StepI
 
     unsigned int id{1};
     Standard_Integer topShapeLength = topLevelShapes.Length() + 1;
+    auto stage_unit2 = topShapeLength / LOAD_STEP_STAGE_UNIT_NUM + 1;
+
     for (Standard_Integer iLabel = 1; iLabel < topShapeLength; ++iLabel) {
-        if (proFn) {
-            proFn(LOAD_STEP_STAGE_GET_SOLID, iLabel, topShapeLength, cb_cancel);
+        if (stepFn) {
+            if ((iLabel % stage_unit2) == 0) {
+                stepFn(LOAD_STEP_STAGE_GET_SOLID, iLabel, topShapeLength, cb_cancel);
+            }
             if (cb_cancel) {
                 shapeTool.reset(nullptr);
                 application->Close(document);
@@ -257,14 +258,94 @@ bool load_step(const char *path, Model *model, ImportStepProgressFn proFn, StepI
         getNamedSolids(TopLoc_Location{}, "", id, shapeTool, topLevelShapes.Value(iLabel), namedSolids);
     }
 
-    ModelObject* new_object = model->add_object();
-    const char *last_slash = strrchr(path, DIR_SEPARATOR);
+    std::vector<stl_file> stl;
+    stl.resize(namedSolids.size());
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, namedSolids.size()), [&](const tbb::blocked_range<size_t> &range) {
+        for (size_t i = range.begin(); i < range.end(); i++) {
+            BRepMesh_IncrementalMesh mesh(namedSolids[i].solid, STEP_TRANS_CHORD_ERROR, false, STEP_TRANS_ANGLE_RES, true);
+            // BBS: calculate total number of the nodes and triangles
+            int aNbNodes     = 0;
+            int aNbTriangles = 0;
+            for (TopExp_Explorer anExpSF(namedSolids[i].solid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
+                TopLoc_Location aLoc;
+                Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(anExpSF.Current()), aLoc);
+                if (!aTriangulation.IsNull()) {
+                    aNbNodes += aTriangulation->NbNodes();
+                    aNbTriangles += aTriangulation->NbTriangles();
+                }
+            }
+
+            if (aNbTriangles == 0)
+                // BBS: No triangulation on the shape.
+                continue;
+
+            stl[i].stats.type                = inmemory;
+            stl[i].stats.number_of_facets    = (uint32_t) aNbTriangles;
+            stl[i].stats.original_num_facets = stl[i].stats.number_of_facets;
+            stl_allocate(&stl[i]);
+
+            std::vector<Vec3f> points;
+            points.reserve(aNbNodes);
+            // BBS: count faces missing triangulation
+            Standard_Integer aNbFacesNoTri = 0;
+            // BBS: fill temporary triangulation
+            Standard_Integer aNodeOffset    = 0;
+            Standard_Integer aTriangleOffet = 0;
+            for (TopExp_Explorer anExpSF(namedSolids[i].solid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
+                const TopoDS_Shape &aFace = anExpSF.Current();
+                TopLoc_Location     aLoc;
+                Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(aFace), aLoc);
+                if (aTriangulation.IsNull()) {
+                    ++aNbFacesNoTri;
+                    continue;
+                }
+                // BBS: copy nodes
+                gp_Trsf aTrsf = aLoc.Transformation();
+                for (Standard_Integer aNodeIter = 1; aNodeIter <= aTriangulation->NbNodes(); ++aNodeIter) {
+                    gp_Pnt aPnt = aTriangulation->Node(aNodeIter);
+                    aPnt.Transform(aTrsf);
+                    points.emplace_back(std::move(Vec3f(aPnt.X(), aPnt.Y(), aPnt.Z())));
+                }
+                // BBS: copy triangles
+                const TopAbs_Orientation anOrientation = anExpSF.Current().Orientation();
+                Standard_Integer anId[3];
+                for (Standard_Integer aTriIter = 1; aTriIter <= aTriangulation->NbTriangles(); ++aTriIter) {
+                    Poly_Triangle aTri = aTriangulation->Triangle(aTriIter);
+
+                    aTri.Get(anId[0], anId[1], anId[2]);
+                    if (anOrientation == TopAbs_REVERSED)
+                        std::swap(anId[1], anId[2]);
+                    // BBS: save triangles facets
+                    stl_facet facet;
+                    facet.vertex[0] = points[anId[0] + aNodeOffset - 1].cast<float>();
+                    facet.vertex[1] = points[anId[1] + aNodeOffset - 1].cast<float>();
+                    facet.vertex[2] = points[anId[2] + aNodeOffset - 1].cast<float>();
+                    facet.extra[0]  = 0;
+                    facet.extra[1]  = 0;
+                    stl_normal normal;
+                    stl_calculate_normal(normal, &facet);
+                    stl_normalize_vector(normal);
+                    facet.normal                                      = normal;
+                    stl[i].facet_start[aTriangleOffet + aTriIter - 1] = facet;
+                }
+
+                aNodeOffset += aTriangulation->NbNodes();
+                aTriangleOffet += aTriangulation->NbTriangles();
+            }
+        }
+    });
+
+    ModelObject *new_object = model->add_object();
+    const char * last_slash = strrchr(path, DIR_SEPARATOR);
     new_object->name.assign((last_slash == nullptr) ? path : last_slash + 1);
     new_object->input_file = path;
 
-    for (size_t i = 0; i < namedSolids.size(); ++i) {
-        if (proFn) {
-            proFn(LOAD_STEP_STAGE_GET_MESH, i, namedSolids.size(), cb_cancel);
+    auto stage_unit3 = stl.size() / LOAD_STEP_STAGE_UNIT_NUM + 1;
+    for (size_t i = 0; i < stl.size(); i++) {
+        if (stepFn) {
+            if ((i % stage_unit3) == 0) {
+                stepFn(LOAD_STEP_STAGE_GET_MESH, i, stl.size(), cb_cancel);
+            }
             if (cb_cancel) {
                 model->delete_object(new_object);
                 shapeTool.reset(nullptr);
@@ -273,94 +354,13 @@ bool load_step(const char *path, Model *model, ImportStepProgressFn proFn, StepI
             }
         }
 
-        BRepMesh_IncrementalMesh mesh(namedSolids[i].solid, STEP_TRANS_CHORD_ERROR, false, STEP_TRANS_ANGLE_RES, true);
-        //BBS: calculate total number of the nodes and triangles
-        int aNbNodes = 0;
-        int aNbTriangles = 0;
-        for (TopExp_Explorer anExpSF(namedSolids[i].solid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
-            TopLoc_Location aLoc;
-            Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(anExpSF.Current()), aLoc);
-            if (!aTriangulation.IsNull()) {
-                aNbNodes += aTriangulation->NbNodes();
-                aNbTriangles += aTriangulation->NbTriangles();
-            }
-        }
-
-        if (aNbTriangles == 0) {
-            //BBS: No triangulation on the shape.
-            continue;
-        }
-
-        stl_file stl;
-        stl.stats.type = inmemory;
-        stl.stats.number_of_facets = (uint32_t)aNbTriangles;
-        stl.stats.original_num_facets = stl.stats.number_of_facets;
-        stl_allocate(&stl);
-
-        std::vector<Vec3f> points;
-        points.reserve(aNbNodes);
-        //BBS: count faces missing triangulation
-        Standard_Integer aNbFacesNoTri = 0;
-        //BBS: fill temporary triangulation
-        Standard_Integer aNodeOffset = 0;
-        Standard_Integer aTriangleOffet = 0;
-        for (TopExp_Explorer anExpSF(namedSolids[i].solid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
-            const TopoDS_Shape& aFace = anExpSF.Current();
-            TopLoc_Location aLoc;
-            Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(aFace), aLoc);
-            if (aTriangulation.IsNull()) {
-                ++aNbFacesNoTri;
-                continue;
-            }
-            //BBS: copy nodes
-            gp_Trsf aTrsf = aLoc.Transformation();
-            for (Standard_Integer aNodeIter = 1; aNodeIter <= aTriangulation->NbNodes(); ++aNodeIter) {
-                gp_Pnt aPnt = aTriangulation->Node(aNodeIter);
-                aPnt.Transform(aTrsf);
-                points.emplace_back(std::move(Vec3f(aPnt.X(), aPnt.Y(), aPnt.Z())));
-            }
-            //BBS: copy triangles
-            const TopAbs_Orientation anOrientation = anExpSF.Current().Orientation();
-            for (Standard_Integer aTriIter = 1; aTriIter <= aTriangulation->NbTriangles(); ++aTriIter) {
-                Poly_Triangle aTri = aTriangulation->Triangle(aTriIter);
-
-                Standard_Integer anId[3];
-                aTri.Get(anId[0], anId[1], anId[2]);
-                if (anOrientation == TopAbs_REVERSED) {
-                    //BBS: swap 1, 2.
-                    Standard_Integer aTmpIdx = anId[1];
-                    anId[1] = anId[2];
-                    anId[2] = aTmpIdx;
-                }
-                //BBS: Update nodes according to the offset.
-                anId[0] += aNodeOffset;
-                anId[1] += aNodeOffset;
-                anId[2] += aNodeOffset;
-                //BBS: save triangles facets
-                stl_facet facet;
-                facet.vertex[0] = points[anId[0] - 1].cast<float>();
-                facet.vertex[1] = points[anId[1] - 1].cast<float>();
-                facet.vertex[2] = points[anId[2] - 1].cast<float>();
-                facet.extra[0] = 0;
-                facet.extra[1] = 0;
-                stl_normal normal;
-                stl_calculate_normal(normal, &facet);
-                stl_normalize_vector(normal);
-                facet.normal = normal;
-                stl.facet_start[aTriangleOffet + aTriIter - 1] = facet;
-            }
-
-            aNodeOffset += aTriangulation->NbNodes();
-            aTriangleOffet += aTriangulation->NbTriangles();
-        }
-
         TriangleMesh triangle_mesh;
-        triangle_mesh.from_stl(stl);
-        ModelVolume* new_volume = new_object->add_volume(std::move(triangle_mesh));
-        new_volume->name = namedSolids[i].name;
+        triangle_mesh.from_stl(stl[i]);
+        ModelVolume *new_volume       = new_object->add_volume(std::move(triangle_mesh));
+        new_volume->name              = namedSolids[i].name;
         new_volume->source.input_file = path;
-        new_volume->source.object_idx = (int)model->objects.size() - 1;
-        new_volume->source.volume_idx = (int)new_object->volumes.size() - 1;
+        new_volume->source.object_idx = (int) model->objects.size() - 1;
+        new_volume->source.volume_idx = (int) new_object->volumes.size() - 1;
     }
 
     shapeTool.reset(nullptr);

src/libslic3r/Format/STEP.hpp

@@ -6,6 +6,13 @@ namespace Slic3r {
 class TriangleMesh;
 class ModelObject;
 
+// load step stage
+const int LOAD_STEP_STAGE_READ_FILE          = 0;
+const int LOAD_STEP_STAGE_GET_SOLID          = 1;
+const int LOAD_STEP_STAGE_GET_MESH           = 2;
+const int LOAD_STEP_STAGE_NUM                = 3;
+const int LOAD_STEP_STAGE_UNIT_NUM           = 5;
+
 typedef std::function<void(int load_stage, int current, int total, bool& cancel)> ImportStepProgressFn;
 typedef std::function<void(bool isUtf8)> StepIsUtf8Fn;
 

src/libslic3r/Format/STL.cpp

@@ -14,10 +14,10 @@
 
 namespace Slic3r {
 
-bool load_stl(const char *path, Model *model, const char *object_name_in)
+bool load_stl(const char *path, Model *model, const char *object_name_in, ImportstlProgressFn stlFn)
 {
     TriangleMesh mesh;
-    if (! mesh.ReadSTLFile(path)) {
+    if (! mesh.ReadSTLFile(path, true, stlFn)) {
 //    die "Failed to open $file\n" if !-e $path;
         return false;
     }

src/libslic3r/Format/STL.hpp

@@ -1,13 +1,16 @@
 #ifndef slic3r_Format_STL_hpp_
 #define slic3r_Format_STL_hpp_
 
+#include <admesh/stl.h>
+
 namespace Slic3r {
 
+class Model;
 class TriangleMesh;
 class ModelObject;
 
 // Load an STL file into a provided model.
-extern bool load_stl(const char *path, Model *model, const char *object_name = nullptr);
+extern bool load_stl(const char *path, Model *model, const char *object_name = nullptr, ImportstlProgressFn stlFn = nullptr);
 
 extern bool store_stl(const char *path, TriangleMesh *mesh, bool binary);
 extern bool store_stl(const char *path, ModelObject *model_object, bool binary);

src/libslic3r/Format/bbs_3mf.cpp

@@ -258,6 +258,8 @@ static constexpr const char* SOURCE_OFFSET_Z_KEY = "source_offset_z";
 static constexpr const char* SOURCE_IN_INCHES    = "source_in_inches";
 static constexpr const char* SOURCE_IN_METERS    = "source_in_meters";
 
+static constexpr const char* MESH_SHARED_KEY = "mesh_shared";
+
 static constexpr const char* MESH_STAT_EDGES_FIXED          = "edges_fixed";
 static constexpr const char* MESH_STAT_DEGENERATED_FACETS   = "degenerate_facets";
 static constexpr const char* MESH_STAT_FACETS_REMOVED       = "facets_removed";
@@ -702,6 +704,8 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
         std::string m_thumbnail_path;
         std::vector<std::string> m_sub_model_paths;
 
+        std::map<int, ModelVolume*> m_shared_meshes;
+
         //BBS: plater related structures
         bool m_is_bbl_3mf { false };
         bool m_parsing_slice_info { false };
@@ -847,8 +851,8 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
 
         bool _handle_start_relationship(const char** attributes, unsigned int num_attributes);
 
-        void _generate_current_object_list(std::vector<Id> &sub_objects, Id object_id, IdToCurrentObjectMap current_objects);
-        bool _generate_volumes_new(ModelObject& object, const std::vector<Id> &sub_objects, const ObjectMetadata::VolumeMetadataList& volumes, ConfigSubstitutionContext& config_substitutions);
+        void _generate_current_object_list(std::vector<Component> &sub_objects, Id object_id, IdToCurrentObjectMap current_objects);
+        bool _generate_volumes_new(ModelObject& object, const std::vector<Component> &sub_objects, const ObjectMetadata::VolumeMetadataList& volumes, ConfigSubstitutionContext& config_substitutions);
         bool _generate_volumes(ModelObject& object, const Geometry& geometry, const ObjectMetadata::VolumeMetadataList& volumes, ConfigSubstitutionContext& config_substitutions);
 
         // callbacks to parse the .model file
@@ -1281,7 +1285,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                         add_error("3rd 3mf, can not find object, id " + std::to_string(object.first.second));
                         return false;
                     }
-                    std::vector<Id> object_id_list;
+                    std::vector<Component> object_id_list;
                     _generate_current_object_list(object_id_list, object.first, m_current_objects);
 
                     ObjectMetadata::VolumeMetadataList volumes;
@@ -1289,7 +1293,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
 
                     for (int k = 0; k < object_id_list.size(); k++)
                     {
-                        Id object_id = object_id_list[k];
+                        Id object_id = object_id_list[k].object_id;
                         volumes.emplace_back(object_id.second);
                     }
 
@@ -1344,7 +1348,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                 model_object->sla_drain_holes = std::move(obj_drain_holes->second);
             }*/
 
-            std::vector<Id> object_id_list;
+            std::vector<Component> object_id_list;
             _generate_current_object_list(object_id_list, object.first, m_current_objects);
 
             ObjectMetadata::VolumeMetadataList volumes;
@@ -1375,7 +1379,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                 //volumes.emplace_back(0, (int)obj_geometry->second.triangles.size() - 1);
                 for (int k = 0; k < object_id_list.size(); k++)
                 {
-                    Id object_id = object_id_list[k];
+                    Id object_id = object_id_list[k].object_id;
                     volumes.emplace_back(object_id.second);
                 }
 
@@ -2757,31 +2761,31 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
             check_painting_version(m_mm_painting_version, MM_PAINTING_VERSION,
                 _(L("The selected 3MF contains multi-material painted object using a newer version of BambuStudio and is not compatible.")));*/
         } else if (m_curr_metadata_name == BBL_MODEL_ID_TAG) {
-            m_model_id = m_curr_characters;
+            m_model_id = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_MODEL_NAME_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found model name = " << m_curr_characters;
-            model_info.model_name = m_curr_characters;
+            model_info.model_name = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_DESIGNER_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found designer = " << m_curr_characters;
-            m_designer = m_curr_characters;
+            m_designer = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_DESIGNER_USER_ID_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found designer_user_id = " << m_curr_characters;
-            m_designer_user_id = m_curr_characters;
+            m_designer_user_id = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_DESIGNER_COVER_FILE_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found designer_cover = " << m_curr_characters;
-            model_info.cover_file = m_curr_characters;
+            model_info.cover_file = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_DESCRIPTION_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found description = " << m_curr_characters;
-            model_info.description = m_curr_characters;
+            model_info.description = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_LICENSE_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found license = " << m_curr_characters;
-            model_info.license = m_curr_characters;
+            model_info.license = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_COPYRIGHT_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found copyright = " << m_curr_characters;
-            model_info.copyright = m_curr_characters;
+            model_info.copyright = xml_unescape(m_curr_characters);
         } else if (m_curr_metadata_name == BBL_REGION_TAG) {
             BOOST_LOG_TRIVIAL(trace) << "design_info, load_3mf found region = " << m_curr_characters;
-            m_contry_code = m_curr_characters;
+            m_contry_code = xml_unescape(m_curr_characters);
         }
 
         return true;
@@ -3279,34 +3283,34 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
         return true;
     }
 
-    void _BBS_3MF_Importer::_generate_current_object_list(std::vector<Id> &sub_objects, Id object_id, IdToCurrentObjectMap current_objects)
+    void _BBS_3MF_Importer::_generate_current_object_list(std::vector<Component> &sub_objects, Id object_id, IdToCurrentObjectMap current_objects)
     {
-        std::list<Id> id_list;
-        id_list.push_back(object_id);
+        std::list<Component> id_list;
+        id_list.push_back({ object_id, Transform3d::Identity() });
 
         while (!id_list.empty())
         {
-            Id current_id = id_list.front();
+            Component current_id = id_list.front();
             id_list.pop_front();
-            IdToCurrentObjectMap::iterator current_object = current_objects.find(current_id);
+            IdToCurrentObjectMap::iterator current_object = current_objects.find(current_id.object_id);
             if (current_object != current_objects.end()) {
                 //found one
                 if (!current_object->second.components.empty()) {
                     for (const Component& comp: current_object->second.components)
                     {
-                        id_list.push_back(comp.object_id);
+                        id_list.push_back(comp);
                     }
                 }
                 else if (!(current_object->second.geometry.empty())) {
                     //CurrentObject* ptr = &(current_objects[current_id]);
                     //CurrentObject* ptr2 = &(current_object->second);
-                    sub_objects.push_back(current_object->first);
+                    sub_objects.push_back({ current_object->first, current_id.transform });
                 }
             }
         }
     }
 
-    bool _BBS_3MF_Importer::_generate_volumes_new(ModelObject& object, const std::vector<Id> &sub_objects, const ObjectMetadata::VolumeMetadataList& volumes, ConfigSubstitutionContext& config_substitutions)
+    bool _BBS_3MF_Importer::_generate_volumes_new(ModelObject& object, const std::vector<Component> &sub_objects, const ObjectMetadata::VolumeMetadataList& volumes, ConfigSubstitutionContext& config_substitutions)
     {
         if (!object.volumes.empty()) {
             add_error("object already built with parts");
@@ -3319,7 +3323,8 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
         for (unsigned int index = 0; index < sub_objects.size(); index++)
         {
             //find the volume metadata firstly
-            Id object_id = sub_objects[index];
+            Component sub_comp = sub_objects[index];
+            Id object_id = sub_comp.object_id;
             IdToCurrentObjectMap::iterator current_object = m_current_objects.find(object_id);
             if (current_object == m_current_objects.end()) {
                 add_error("sub_objects can not be found, id=" + std::to_string(object_id.second));
@@ -3338,70 +3343,114 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
 
             Transform3d volume_matrix_to_object = Transform3d::Identity();
             bool        has_transform 		    = false;
+            int         shared_mesh_id = -1;
             if (volume_data)
             {
+                int found_count = 0;
                 // extract the volume transformation from the volume's metadata, if present
                 for (const Metadata& metadata : volume_data->metadata) {
                     if (metadata.key == MATRIX_KEY) {
                         volume_matrix_to_object = Slic3r::Geometry::transform3d_from_string(metadata.value);
                         has_transform 			= ! volume_matrix_to_object.isApprox(Transform3d::Identity(), 1e-10);
-                        break;
+                        found_count++;
                     }
+                    else if (metadata.key == MESH_SHARED_KEY){
+                        //add the shared mesh logic
+                        shared_mesh_id = ::atoi(metadata.value.c_str());
+                        found_count++;
+                    }
+
+                    if (found_count >= 2)
+                        break;
                 }
             }
             else {
                 //create a volume_data
                 volume_data = &default_volume_data;
             }
-            // splits volume out of imported geometry
-            indexed_triangle_set its;
-            its.indices.assign(sub_object->geometry.triangles.begin(), sub_object->geometry.triangles.end());
-            const size_t triangles_count = its.indices.size();
+
+            ModelVolume* volume = nullptr;
+            ModelVolume *shared_volume = nullptr;
+            if (shared_mesh_id != -1) {
+                std::map<int, ModelVolume*>::iterator iter = m_shared_meshes.find(shared_mesh_id);
+                if (iter != m_shared_meshes.end()) {
+                    shared_volume = iter->second;
+                }
+            }
+
+            const size_t triangles_count = sub_object->geometry.triangles.size();
             if (triangles_count == 0) {
                 add_error("found no trianges in the object " + std::to_string(sub_object->id));
                 return false;
             }
-            for (const Vec3i& face : its.indices) {
-                for (const int tri_id : face) {
-                    if (tri_id < 0 || tri_id >= int(sub_object->geometry.vertices.size())) {
-                        add_error("invalid vertex id in object " + std::to_string(sub_object->id));
-                        return false;
+            if (!shared_volume){
+                // splits volume out of imported geometry
+                indexed_triangle_set its;
+                its.indices.assign(sub_object->geometry.triangles.begin(), sub_object->geometry.triangles.end());
+                //const size_t triangles_count = its.indices.size();
+                //if (triangles_count == 0) {
+                //    add_error("found no trianges in the object " + std::to_string(sub_object->id));
+                //    return false;
+                //}
+                for (const Vec3i& face : its.indices) {
+                    for (const int tri_id : face) {
+                        if (tri_id < 0 || tri_id >= int(sub_object->geometry.vertices.size())) {
+                            add_error("invalid vertex id in object " + std::to_string(sub_object->id));
+                            return false;
+                        }
                     }
                 }
-            }
 
-            its.vertices.assign(sub_object->geometry.vertices.begin(), sub_object->geometry.vertices.end());
+                its.vertices.assign(sub_object->geometry.vertices.begin(), sub_object->geometry.vertices.end());
 
-            // BBS
-            for (const std::string prop_str : sub_object->geometry.face_properties) {
-                FaceProperty face_prop;
-                face_prop.from_string(prop_str);
-                its.properties.push_back(face_prop);
-            }
-
-            TriangleMesh triangle_mesh(std::move(its), volume_data->mesh_stats);
-
-            if (m_version == 0) {
-                // if the 3mf was not produced by BambuStudio and there is only one instance,
-                // bake the transformation into the geometry to allow the reload from disk command
-                // to work properly
-                if (object.instances.size() == 1) {
-                    triangle_mesh.transform(object.instances.front()->get_transformation().get_matrix(), false);
-                    object.instances.front()->set_transformation(Slic3r::Geometry::Transformation());
-                    //FIXME do the mesh fixing?
+                // BBS
+                for (const std::string prop_str : sub_object->geometry.face_properties) {
+                    FaceProperty face_prop;
+                    face_prop.from_string(prop_str);
+                    its.properties.push_back(face_prop);
                 }
-            }
-            if (triangle_mesh.volume() < 0)
-                triangle_mesh.flip_triangles();
 
-			ModelVolume* volume = object.add_volume(std::move(triangle_mesh));
+                TriangleMesh triangle_mesh(std::move(its), volume_data->mesh_stats);
+
+                if (m_version == 0) {
+                    // if the 3mf was not produced by BambuStudio and there is only one instance,
+                    // bake the transformation into the geometry to allow the reload from disk command
+                    // to work properly
+                    if (object.instances.size() == 1) {
+                        triangle_mesh.transform(object.instances.front()->get_transformation().get_matrix(), false);
+                        object.instances.front()->set_transformation(Slic3r::Geometry::Transformation());
+                        //FIXME do the mesh fixing?
+                    }
+                }
+                if (triangle_mesh.volume() < 0)
+                    triangle_mesh.flip_triangles();
+
+                volume = object.add_volume(std::move(triangle_mesh));
+
+                m_shared_meshes[sub_object->id] = volume;
+            }
+            else {
+                //create volume to use shared mesh
+                volume = object.add_volume_with_shared_mesh(*shared_volume);
+            }
             // stores the volume matrix taken from the metadata, if present
             if (has_transform)
                 volume->source.transform = Slic3r::Geometry::Transformation(volume_matrix_to_object);
+
             volume->calculate_convex_hull();
 
+            //set transform from 3mf
+            Slic3r::Geometry::Transformation comp_transformatino(sub_comp.transform);
+            volume->set_transformation(volume->get_transformation() * comp_transformatino);
+            if (shared_volume) {
+                const TriangleMesh& trangle_mesh = volume->mesh();
+                Vec3d shift = trangle_mesh.get_init_shift();
+                if (!shift.isApprox(Vec3d::Zero()))
+                    volume->translate(shift);
+            }
+
             // recreate custom supports, seam and mmu segmentation from previously loaded attribute
-            if (m_load_config) {
+            {
                 volume->supported_facets.reserve(triangles_count);
                 volume->seam_facets.reserve(triangles_count);
                 volume->mmu_segmentation_facets.reserve(triangles_count);
@@ -3448,7 +3497,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                     volume->source.is_converted_from_inches = metadata.value == "1";
                 else if (metadata.key == SOURCE_IN_METERS)
                     volume->source.is_converted_from_meters = metadata.value == "1";
-                else if (metadata.key == MATRIX_KEY)
+                else if ((metadata.key == MATRIX_KEY) || (metadata.key == MESH_SHARED_KEY))
                     continue;
                 else
                     volume->config.set_deserialize(metadata.key, metadata.value, config_substitutions);
@@ -3814,18 +3863,38 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
 
         auto filename = boost::format("3D/Objects/%s_%d.model") % object.name % obj_id;
         std::string filepath = temp_path + "/" + filename.str();
-        if (!open_zip_writer(&archive, filepath)) {
+        std::string filepath_tmp = filepath + ".tmp";
+        boost::system::error_code ec;
+        boost::filesystem::remove(filepath_tmp, ec);
+        if (!open_zip_writer(&archive, filepath_tmp)) {
             add_error("Unable to open the file");
             BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << ":" << __LINE__ << boost::format(", Unable to open the file\n");
             return false;
         }
 
+        struct close_lock
+        {
+            mz_zip_archive & archive;
+            std::string const * filename;
+            void close() {
+                close_zip_writer(&archive);
+                filename = nullptr;
+            }
+            ~close_lock() {
+                if (filename) {
+                    close_zip_writer(&archive);
+                    boost::filesystem::remove(*filename);
+                }
+            }
+        } lock{archive, &filepath_tmp};
+
         IdToObjectDataMap objects_data;
         objects_data.insert({obj_id, {&object, obj_id}});
         _add_model_file_to_archive(filename.str(), archive, model, objects_data);
 
         mz_zip_writer_finalize_archive(&archive);
-        close_zip_writer(&archive);
+        lock.close();
+        boost::filesystem::rename(filepath_tmp, filepath, ec);
         return true;
     }
 
@@ -4363,6 +4432,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
             return false;
         }
 
+
         {
             std::stringstream stream;
             reset_stream(stream);
@@ -4413,27 +4483,29 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                 model_id = project->project_model_id;
                 region_code = project->project_country_code;
             }
+            if (!sub_model) {
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_MODEL_NAME_TAG          << "\">" << xml_escape(name)         << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_TAG            << "\">" << xml_escape(user_name)    << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_USER_ID_TAG    << "\">" << user_id                  << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_COVER_FILE_TAG << "\">" << xml_escape(design_cover) << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_DESCRIPTION_TAG         << "\">" << xml_escape(description)  << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_COPYRIGHT_TAG           << "\">" << xml_escape(copyright)    << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"" << BBL_LICENSE_TAG             << "\">" << xml_escape(license)      << "</" << METADATA_TAG << ">\n";
 
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_MODEL_NAME_TAG          << "\">" << xml_escape(name)         << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_TAG            << "\">" << xml_escape(user_name)    << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_USER_ID_TAG    << "\">" << user_id                  << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_DESIGNER_COVER_FILE_TAG << "\">" << xml_escape(design_cover) << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_DESCRIPTION_TAG         << "\">" << xml_escape(description)  << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_COPYRIGHT_TAG           << "\">" << xml_escape(copyright)    << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"" << BBL_LICENSE_TAG             << "\">" << xml_escape(license)      << "</" << METADATA_TAG << ">\n";
+                /* save model info */
+                if (!model_id.empty()) {
+                    stream << " <" << METADATA_TAG << " name=\"" << BBL_MODEL_ID_TAG << "\">" << model_id    << "</" << METADATA_TAG << ">\n";
+                    stream << " <" << METADATA_TAG << " name=\"" << BBL_REGION_TAG << "\">"   << region_code << "</" << METADATA_TAG << ">\n";
+                }
 
-            /* save model info */
-            if (!model_id.empty()) {
-                stream << " <" << METADATA_TAG << " name=\"" << BBL_MODEL_ID_TAG << "\">" << model_id    << "</" << METADATA_TAG << ">\n";
-                stream << " <" << METADATA_TAG << " name=\"" << BBL_REGION_TAG << "\">"   << region_code << "</" << METADATA_TAG << ">\n";
+                std::string date = Slic3r::Utils::utc_timestamp(Slic3r::Utils::get_current_time_utc());
+                // keep only the date part of the string
+                date = date.substr(0, 10);
+                stream << " <" << METADATA_TAG << " name=\"CreationDate\">" << date << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"ModificationDate\">" << date << "</" << METADATA_TAG << ">\n";
+                stream << " <" << METADATA_TAG << " name=\"Application\">" << SLIC3R_APP_KEY << "-" << SLIC3R_VERSION << "</" << METADATA_TAG << ">\n";
             }
 
-            std::string date = Slic3r::Utils::utc_timestamp(Slic3r::Utils::get_current_time_utc());
-            // keep only the date part of the string
-            date = date.substr(0, 10);
-            stream << " <" << METADATA_TAG << " name=\"CreationDate\">" << date << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"ModificationDate\">" << date << "</" << METADATA_TAG << ">\n";
-            stream << " <" << METADATA_TAG << " name=\"Application\">" << SLIC3R_APP_KEY << "-" << SLIC3R_VERSION << "</" << METADATA_TAG << ">\n";
             stream << " <" << RESOURCES_TAG << ">\n";
             std::string buf = stream.str();
             if (! buf.empty() && ! mz_zip_writer_add_staged_data(&context, buf.data(), buf.size())) {
@@ -4645,12 +4717,37 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                 if (m_from_backup_save) {
                     for (unsigned int index = 1; index <= object.volumes.size(); index ++) {
                         unsigned int ref_id = object_id | (index << 16);
-                        stream << "    <" << COMPONENT_TAG << " objectid=\"" << ref_id << "\"/>\n";
+                        stream << "    <" << COMPONENT_TAG << " objectid=\"" << ref_id; // << "\"/>\n";
+                        //add the transform of the volume
+                        ModelVolume* volume = object.volumes[index - 1];
+                        const Transform3d& transf = volume->get_matrix();
+                        stream << "\" " << TRANSFORM_ATTR << "=\"";
+                        for (unsigned c = 0; c < 4; ++c) {
+                            for (unsigned r = 0; r < 3; ++r) {
+                                stream << transf(r, c);
+                                if (r != 2 || c != 3)
+                                    stream << " ";
+                            }
+                        }
+                        stream << "\"/>\n";
                     }
                 }
                 else {
-                    for (unsigned int index = object_id; index < volume_start_id; index ++)
-                        stream << "    <" << COMPONENT_TAG << " objectid=\"" << index << "\"/>\n";
+                    for (unsigned int index = object_id; index < volume_start_id; index ++) {
+                        stream << "    <" << COMPONENT_TAG << " objectid=\"" << index; // << "\"/>\n";
+                        //add the transform of the volume
+                        ModelVolume* volume = object.volumes[index - object_id];
+                        const Transform3d& transf = volume->get_matrix();
+                        stream << "\" " << TRANSFORM_ATTR << "=\"";
+                        for (unsigned c = 0; c < 4; ++c) {
+                            for (unsigned r = 0; r < 3; ++r) {
+                                stream << transf(r, c);
+                                if (r != 2 || c != 3)
+                                    stream << " ";
+                            }
+                        }
+                        stream << "\"/>\n";
+                    }
                 }
             }
             else {
@@ -4800,7 +4897,10 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
             const Transform3d& matrix = volume->get_matrix();
 
             for (size_t i = 0; i < its.vertices.size(); ++i) {
-                Vec3f v = (matrix * its.vertices[i].cast<double>()).cast<float>();
+                //don't save the volume's matrix into vertex data
+                //add the shared mesh logic
+                //Vec3f v = (matrix * its.vertices[i].cast<double>()).cast<float>();
+                Vec3f v = its.vertices[i];
                 char* ptr = buf;
                 boost::spirit::karma::generate(ptr, boost::spirit::lit("     <") << VERTEX_TAG << " x=\"");
                 ptr = format_coordinate(v.x(), ptr);
@@ -5202,6 +5302,7 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
     bool _BBS_3MF_Exporter::_add_model_config_file_to_archive(mz_zip_archive& archive, const Model& model, PlateDataPtrs& plate_data_list, const IdToObjectDataMap &objects_data, int export_plate_idx, bool save_gcode)
     {
         std::stringstream stream;
+        std::map<const TriangleMesh*, int> shared_meshes;
         // Store mesh transformation in full precision, as the volumes are stored transformed and they need to be transformed back
         // when loaded as accurately as possible.
 		stream << std::setprecision(std::numeric_limits<double>::max_digits10);
@@ -5294,6 +5395,17 @@ void PlateData::parse_filament_info(GCodeProcessorResult *result)
                                 stream << "      <" << METADATA_TAG << " "<< KEY_ATTR << "=\"" << key << "\" " << VALUE_ATTR << "=\"" << volume->config.opt_serialize(key) << "\"/>\n";
                             }
 
+                            //add the shared mesh logic
+                            const TriangleMesh* current_mesh = volume->mesh_ptr();
+                            std::map<const TriangleMesh*,int>::iterator mesh_iter;
+                            mesh_iter = shared_meshes.find(current_mesh);
+                            if (mesh_iter != shared_meshes.end()) {
+                                stream << "      <" << METADATA_TAG << " "<< KEY_ATTR << "=\"" << MESH_SHARED_KEY << "\" " << VALUE_ATTR << "=\"" << mesh_iter->second << "\"/>\n";
+                            }
+                            else {
+                                shared_meshes[current_mesh] = it->second;
+                            }
+
                             // stores mesh's statistics
                             const RepairedMeshErrors& stats = volume->mesh().stats().repaired_errors;
                             stream << "      <" << MESH_STAT_TAG << " ";

src/libslic3r/Format/bbs_3mf.hpp

@@ -166,6 +166,7 @@ const int IMPORT_STAGE_CHECK_MODE_GCODE = 9;
 const int UPDATE_GCODE_RESULT           = 10;
 const int IMPORT_LOAD_CONFIG            = 11;
 const int IMPORT_LOAD_MODEL_OBJECTS     = 12;
+const int IMPORT_STAGE_MAX              = 13;
 
 //BBS export 3mf progress
 typedef std::function<void(int export_stage, int current, int total, bool& cancel)> Export3mfProgressFn;

src/libslic3r/GCode.cpp

@@ -554,6 +554,9 @@ bool GCode::gcode_label_objects = false;
     {
         std::string gcode;
 
+        assert(m_layer_idx >= 0);
+        if (m_layer_idx >= (int) m_tool_changes.size()) return gcode;
+
         // Calculate where the wipe tower layer will be printed. -1 means that print z will not change,
         // resulting in a wipe tower with sparse layers.
         double wipe_tower_z  = -1;
@@ -571,16 +574,12 @@ bool GCode::gcode_label_objects = false;
             m_is_first_print = false;
         }
 
-        assert(m_layer_idx >= 0);
         if (gcodegen.writer().need_toolchange(extruder_id) || finish_layer) {
-            if (m_layer_idx < (int)m_tool_changes.size()) {
-                if (!(size_t(m_tool_change_idx) < m_tool_changes[m_layer_idx].size()))
-                    throw Slic3r::RuntimeError("Wipe tower generation failed, possibly due to empty first layer.");
+            if (!(size_t(m_tool_change_idx) < m_tool_changes[m_layer_idx].size())) throw Slic3r::RuntimeError("Wipe tower generation failed, possibly due to empty first layer.");
 
-                if (!ignore_sparse) {
-                    gcode += append_tcr(gcodegen, m_tool_changes[m_layer_idx][m_tool_change_idx++], extruder_id, wipe_tower_z);
-                    m_last_wipe_tower_print_z = wipe_tower_z;
-                }
+            if (!ignore_sparse) {
+                gcode += append_tcr(gcodegen, m_tool_changes[m_layer_idx][m_tool_change_idx++], extruder_id, wipe_tower_z);
+                m_last_wipe_tower_print_z = wipe_tower_z;
             }
         }
 
@@ -667,6 +666,7 @@ std::vector<GCode::LayerToPrint> GCode::collect_layers_to_print(const PrintObjec
             --idx_tree_support_layer;
         }
 
+        layer_to_print.original_object = &object;
         layers_to_print.push_back(layer_to_print);
 
         bool has_extrusions = (layer_to_print.object_layer && layer_to_print.object_layer->has_extrusions())
@@ -1261,6 +1261,7 @@ enum BambuBedType {
     bbtCoolPlate = 1,
     bbtEngineeringPlate = 2,
     bbtHighTemperaturePlate = 3,
+    bbtTexturedPEIPlate         = 4,
 };
 
 static BambuBedType to_bambu_bed_type(BedType type)
@@ -1272,6 +1273,8 @@ static BambuBedType to_bambu_bed_type(BedType type)
         bambu_bed_type = bbtEngineeringPlate;
     else if (type == btPEI)
         bambu_bed_type = bbtHighTemperaturePlate;
+    else if (type == btPTE)
+        bambu_bed_type = bbtTexturedPEIPlate;
 
     return bambu_bed_type;
 }
@@ -1570,7 +1573,8 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
     print.throw_if_canceled();
 
     // Collect custom seam data from all objects.
-    m_seam_placer.init(print);
+    std::function<void(void)> throw_if_canceled_func = [&print]() { print.throw_if_canceled(); };
+    m_seam_placer.init(print, throw_if_canceled_func);
 
     // BBS: priming logic is removed, always set first extruer here.
     //if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming))
@@ -2153,7 +2157,9 @@ std::vector<GCode::InstanceToPrint> GCode::sort_print_object_instances(
         // Sequential print, single object is being printed.
         for (ObjectByExtruder &object_by_extruder : objects_by_extruder) {
             const size_t       layer_id     = &object_by_extruder - objects_by_extruder.data();
-            const PrintObject *print_object = layers[layer_id].object();
+            //BBS:add the support of shared print object
+            const PrintObject *print_object = layers[layer_id].original_object;
+            //const PrintObject *print_object = layers[layer_id].object();
             if (print_object)
                 out.emplace_back(object_by_extruder, layer_id, *print_object, single_object_instance_idx);
         }
@@ -2163,7 +2169,9 @@ std::vector<GCode::InstanceToPrint> GCode::sort_print_object_instances(
         sorted.reserve(objects_by_extruder.size());
         for (ObjectByExtruder &object_by_extruder : objects_by_extruder) {
             const size_t       layer_id     = &object_by_extruder - objects_by_extruder.data();
-            const PrintObject *print_object = layers[layer_id].object();
+            //BBS:add the support of shared print object
+            const PrintObject *print_object = layers[layer_id].original_object;
+            //const PrintObject *print_object = layers[layer_id].object();
             if (print_object)
                 sorted.emplace_back(print_object, &object_by_extruder);
         }
@@ -2173,6 +2181,10 @@ std::vector<GCode::InstanceToPrint> GCode::sort_print_object_instances(
             out.reserve(sorted.size());
             for (const PrintInstance *instance : *ordering) {
                 const PrintObject &print_object = *instance->print_object;
+                //BBS:add the support of shared print object
+                //const PrintObject* print_obj_ptr = &print_object;
+                //if (print_object.get_shared_object())
+                //    print_obj_ptr = print_object.get_shared_object();
                 std::pair<const PrintObject*, ObjectByExtruder*> key(&print_object, nullptr);
                 auto it = std::lower_bound(sorted.begin(), sorted.end(), key);
                 if (it != sorted.end() && it->first == &print_object)
@@ -2795,7 +2807,6 @@ GCode::LayerResult GCode::process_layer(
         m_wipe_tower->set_is_first_print(true);
 
     // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders.
-    std::vector<std::unique_ptr<EdgeGrid::Grid>> lower_layer_edge_grids(layers.size());
     for (unsigned int extruder_id : layer_tools.extruders)
     {
         gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ?
@@ -2910,12 +2921,7 @@ GCode::LayerResult GCode::process_layer(
                         m_layer = layers[instance_to_print.layer_id].tree_support_layer;
                     }
                     m_object_layer_over_raft = false;
-                    // BBS. Keep paths order
-#if 0
-                    gcode += this->extrude_support(
-                        // support_extrusion_role is erSupportMaterial, erSupportTransition, erSupportMaterialInterface or erMixed for all extrusion paths.
-                        instance_to_print.object_by_extruder.support->chained_path_from(m_last_pos, instance_to_print.object_by_extruder.support_extrusion_role));
-#else
+
                     //BBS: print supports' brims first
                     if (this->m_objSupportsWithBrim.find(instance_to_print.print_object.id()) != this->m_objSupportsWithBrim.end() && !print_wipe_extrusions) {
                         this->set_origin(0., 0.);
@@ -2945,12 +2951,10 @@ GCode::LayerResult GCode::process_layer(
                     ExtrusionRole support_extrusion_role = instance_to_print.object_by_extruder.support_extrusion_role;
                     bool is_overridden = support_extrusion_role == erSupportMaterialInterface ? support_intf_overridden : support_overridden;
                     if (is_overridden == (print_wipe_extrusions != 0))
-                        support_eec.entities = filter_by_extrusion_role(instance_to_print.object_by_extruder.support->entities, instance_to_print.object_by_extruder.support_extrusion_role);
+                        gcode += this->extrude_support(
+                            // support_extrusion_role is erSupportMaterial, erSupportTransition, erSupportMaterialInterface or erMixed for all extrusion paths.
+                            instance_to_print.object_by_extruder.support->chained_path_from(m_last_pos, support_extrusion_role));
 
-                    for (auto& ptr : support_eec.entities)
-                        ptr = ptr->clone();
-                    gcode += this->extrude_support(support_eec);
-#endif
                     m_layer = layer_to_print.layer();
                     m_object_layer_over_raft = object_layer_over_raft;
                 }
@@ -2984,9 +2988,9 @@ GCode::LayerResult GCode::process_layer(
                     //This behaviour is same with cura
                     if (is_infill_first && !first_layer) {
                         gcode += this->extrude_infill(print, by_region_specific, false);
-                        gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]);
+                        gcode += this->extrude_perimeters(print, by_region_specific);
                     } else {
-                        gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]);
+                        gcode += this->extrude_perimeters(print, by_region_specific);
                         gcode += this->extrude_infill(print,by_region_specific, false);
                     }
                     // ironing
@@ -3164,14 +3168,11 @@ static std::unique_ptr<EdgeGrid::Grid> calculate_layer_edge_grid(const Layer& la
 }
 
 
-std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid)
+std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed)
 {
     // get a copy; don't modify the orientation of the original loop object otherwise
     // next copies (if any) would not detect the correct orientation
 
-    if (m_layer->lower_layer && lower_layer_edge_grid != nullptr && ! *lower_layer_edge_grid)
-        *lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
-
     //BBS: extrude contour of wall ccw, hole of wall cw, except spiral mode
     bool was_clockwise = loop.is_clockwise();
     if (m_config.spiral_mode || !is_perimeter(loop.role()))
@@ -3181,17 +3182,13 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
     // find the point of the loop that is closest to the current extruder position
     // or randomize if requested
     Point last_pos = this->last_pos();
-    if (m_config.spiral_mode) {
+    if (!m_config.spiral_mode && description == "perimeter") {
+        assert(m_layer != nullptr);
+        bool is_outer_wall_first = m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill
+                                || m_config.wall_infill_order == WallInfillOrder::InfillOuterInner;
+        m_seam_placer.place_seam(m_layer, loop, is_outer_wall_first, this->last_pos());
+    } else
         loop.split_at(last_pos, false);
-    }
-    else {
-        //BBS
-        bool is_outer_wall_first =
-                m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill ||
-                m_config.wall_infill_order == WallInfillOrder::InfillOuterInner;
-        m_seam_placer.place_seam(loop, this->last_pos(), is_outer_wall_first,
-            EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr);
-    }
 
     // clip the path to avoid the extruder to get exactly on the first point of the loop;
     // if polyline was shorter than the clipping distance we'd get a null polyline, so
@@ -3301,14 +3298,14 @@ std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string
     return gcode;
 }
 
-std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid)
+std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed)
 {
     if (const ExtrusionPath* path = dynamic_cast<const ExtrusionPath*>(&entity))
         return this->extrude_path(*path, description, speed);
     else if (const ExtrusionMultiPath* multipath = dynamic_cast<const ExtrusionMultiPath*>(&entity))
         return this->extrude_multi_path(*multipath, description, speed);
     else if (const ExtrusionLoop* loop = dynamic_cast<const ExtrusionLoop*>(&entity))
-        return this->extrude_loop(*loop, description, speed, lower_layer_edge_grid);
+        return this->extrude_loop(*loop, description, speed);
     else
         throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()");
     return "";
@@ -3330,24 +3327,15 @@ std::string GCode::extrude_path(ExtrusionPath path, std::string description, dou
 }
 
 // Extrude perimeters: Decide where to put seams (hide or align seams).
-std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid)
+std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region)
 {
     std::string gcode;
     for (const ObjectByExtruder::Island::Region &region : by_region)
         if (! region.perimeters.empty()) {
             m_config.apply(print.get_print_region(&region - &by_region.front()).config());
 
-            // plan_perimeters tries to place seams, it needs to have the lower_layer_edge_grid calculated already.
-            if (m_layer->lower_layer && ! lower_layer_edge_grid)
-                lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
-
-            m_seam_placer.plan_perimeters(std::vector<const ExtrusionEntity*>(region.perimeters.begin(), region.perimeters.end()),
-                *m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter),
-                (m_layer == NULL ? nullptr : m_layer->object()),
-                (lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr));
-
             for (const ExtrusionEntity* ee : region.perimeters)
-                gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid);
+                gcode += this->extrude_entity(*ee, "perimeter", -1.);
         }
     return gcode;
 }
@@ -3817,6 +3805,11 @@ bool GCode::needs_retraction(const Polyline &travel, ExtrusionRole role)
         return false;
     }
 
+    //BBS: force to retract when leave from external perimeter for a long travel
+    //Better way is judging whether the travel move direction is same with last extrusion move.
+    if (is_perimeter(m_last_processor_extrusion_role) && m_last_processor_extrusion_role != erPerimeter)
+        return true;
+
     if (role == erSupportMaterial || role == erSupportTransition) {
         const SupportLayer* support_layer = dynamic_cast<const SupportLayer*>(m_layer);
         //FIXME support_layer->support_islands.contains should use some search structure!

src/libslic3r/GCode.hpp

@@ -41,11 +41,11 @@ class OozePrevention {
 public:
     bool enable;
     Points standby_points;
-    
+
     OozePrevention() : enable(false) {}
     std::string pre_toolchange(GCode &gcodegen);
     std::string post_toolchange(GCode &gcodegen);
-    
+
 private:
     int _get_temp(GCode &gcodegen);
 };
@@ -54,7 +54,7 @@ class Wipe {
 public:
     bool enable;
     Polyline path;
-    
+
     Wipe() : enable(false) {}
     bool has_path() const { return !this->path.points.empty(); }
     void reset_path() { this->path = Polyline(); }
@@ -136,15 +136,15 @@ public:
 };
 
 class GCode {
-public:        
+public:
     GCode() :
     	m_origin(Vec2d::Zero()),
-        m_enable_loop_clipping(true), 
-        m_enable_cooling_markers(false), 
+        m_enable_loop_clipping(true),
+        m_enable_cooling_markers(false),
         m_enable_extrusion_role_markers(false),
         m_last_processor_extrusion_role(erNone),
         m_layer_count(0),
-        m_layer_index(-1), 
+        m_layer_index(-1),
         m_layer(nullptr),
         m_object_layer_over_raft(false),
         //m_volumetric_speed(0),
@@ -201,10 +201,11 @@ public:
     // public, so that it could be accessed by free helper functions from GCode.cpp
     struct LayerToPrint
     {
-        LayerToPrint() : object_layer(nullptr), support_layer(nullptr), tree_support_layer(nullptr) {}
+        LayerToPrint() : object_layer(nullptr), support_layer(nullptr), tree_support_layer(nullptr), original_object(nullptr) {}
         const Layer* 		object_layer;
         const SupportLayer* support_layer;
         const TreeSupportLayer* tree_support_layer;
+        const PrintObject*  original_object; //BBS: used for shared object logic
         const Layer* 		layer()   const
         {
             if (object_layer != nullptr)
@@ -218,7 +219,11 @@ public:
 
             return nullptr;
         }
-        const PrintObject* 	object()  const { return (this->layer() != nullptr) ? this->layer()->object() : nullptr; }
+
+        const PrintObject* 	object()   const
+        {
+            return (this->layer() != nullptr) ? this->layer()->object() : nullptr;
+        }
         coordf_t            print_z() const
         {
             coordf_t sum_z = 0.;
@@ -249,7 +254,7 @@ private:
 
         bool is_open() const { return f; }
         bool is_error() const;
-        
+
         void flush();
         void close();
 
@@ -257,12 +262,12 @@ private:
         void write(const std::string& what) { this->write(what.c_str()); }
         void write(const char* what);
 
-        // Write a string into a file. 
+        // Write a string into a file.
         // Add a newline, if the string does not end with a newline already.
         // Used to export a custom G-code section processed by the PlaceholderParser.
         void writeln(const std::string& what);
 
-        // Formats and write into a file the given data. 
+        // Formats and write into a file the given data.
         void write_format(const char* format, ...);
 
     private:
@@ -325,8 +330,8 @@ private:
     std::string     preamble();
     // BBS
     std::string     change_layer(coordf_t print_z, bool lazy_raise = false);
-    std::string     extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr);
-    std::string     extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr);
+    std::string     extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1.);
+    std::string     extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1.);
     std::string     extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.);
     std::string     extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.);
 
@@ -375,7 +380,7 @@ private:
 		InstanceToPrint(ObjectByExtruder &object_by_extruder, size_t layer_id, const PrintObject &print_object, size_t instance_id) :
 			object_by_extruder(object_by_extruder), layer_id(layer_id), print_object(print_object), instance_id(instance_id) {}
 
-		// Repository 
+		// Repository
 		ObjectByExtruder		&object_by_extruder;
 		// Index into std::vector<LayerToPrint>, which contains Object and Support layers for the current print_z, collected for a single object, or for possibly multiple objects with multiple instances.
 		const size_t       		 layer_id;
@@ -393,7 +398,7 @@ private:
 		// For sequential print, the instance of the object to be printing has to be defined.
 		const size_t                     				 single_object_instance_idx);
 
-    std::string     extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid);
+    std::string     extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region);
     std::string     extrude_infill(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, bool ironing);
     std::string     extrude_support(const ExtrusionEntityCollection &support_fills);
 
@@ -500,14 +505,14 @@ private:
     bool                                object_layer_over_raft() const { return m_object_layer_over_raft; }
 
     friend ObjectByExtruder& object_by_extruder(
-        std::map<unsigned int, std::vector<ObjectByExtruder>> &by_extruder, 
-        unsigned int                                           extruder_id, 
-        size_t                                                 object_idx, 
+        std::map<unsigned int, std::vector<ObjectByExtruder>> &by_extruder,
+        unsigned int                                           extruder_id,
+        size_t                                                 object_idx,
         size_t                                                 num_objects);
     friend std::vector<ObjectByExtruder::Island>& object_islands_by_extruder(
-        std::map<unsigned int, std::vector<ObjectByExtruder>>  &by_extruder, 
-        unsigned int                                            extruder_id, 
-        size_t                                                  object_idx, 
+        std::map<unsigned int, std::vector<ObjectByExtruder>>  &by_extruder,
+        unsigned int                                            extruder_id,
+        size_t                                                  object_idx,
         size_t                                                  num_objects,
         size_t                                                  num_islands);
 

src/libslic3r/GCode/SeamPlacer.hpp

@@ -3,12 +3,16 @@
 
 #include <optional>
 #include <vector>
+#include <memory>
+#include <atomic>
 
+#include "libslic3r/libslic3r.h"
 #include "libslic3r/ExtrusionEntity.hpp"
 #include "libslic3r/Polygon.hpp"
 #include "libslic3r/PrintConfig.hpp"
 #include "libslic3r/BoundingBox.hpp"
 #include "libslic3r/AABBTreeIndirect.hpp"
+#include "libslic3r/KDTreeIndirect.hpp"
 
 namespace Slic3r {
 
@@ -16,119 +20,147 @@ class PrintObject;
 class ExtrusionLoop;
 class Print;
 class Layer;
-namespace EdgeGrid { class Grid; }
-
-
-class SeamHistory {
-public:
-    SeamHistory() { clear(); }
-    std::optional<Point> get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb);
-    void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb);
-    void clear();
-
-private:
-    struct SeamPoint {
-        Point m_pos;
-        BoundingBox m_island_bb;
-    };
-
-    std::map<const PrintObject*, std::vector<SeamPoint>> m_data_last_layer;
-    std::map<const PrintObject*, std::vector<SeamPoint>> m_data_this_layer;
-    size_t m_layer_id;
-};
-
-
-
-class SeamPlacer {
-public:
-    void init(const Print& print);
-
-    // When perimeters are printed, first call this function with the respective
-    // external perimeter. SeamPlacer will find a location for its seam and remember it.
-    // Subsequent calls to get_seam will return this position.
-
-
-    void plan_perimeters(const std::vector<const ExtrusionEntity*> perimeters,
-        const Layer& layer, SeamPosition seam_position,
-        Point last_pos, coordf_t nozzle_dmr, const PrintObject* po,
-        const EdgeGrid::Grid* lower_layer_edge_grid);
-
-    void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
-                    const EdgeGrid::Grid* lower_layer_edge_grid);
-    
-
-    using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
-    using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
-
-private:
-
-    // When given an external perimeter (!), returns the seam.
-    Point calculate_seam(const Layer& layer, const SeamPosition seam_position,
-        const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po,
-        const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos);
-
-    struct CustomTrianglesPerLayer {
-        Polygons polys;
-        TreeType tree;
-    };
-
-    // Just a cache to save some lookups.
-    const Layer* m_last_layer_po = nullptr;
-    coordf_t m_last_print_z = -1.;
-    const PrintObject* m_last_po = nullptr;
-
-    struct SeamPoint {
-        Point pt;
-        bool precalculated = false;
-        bool external = false;
-        const Layer* layer = nullptr;
-        SeamPosition seam_position;
-        const PrintObject* po = nullptr;
-    };
-    std::vector<SeamPoint> m_plan;
-    size_t m_plan_idx;
-
-    std::vector<std::vector<CustomTrianglesPerLayer>> m_enforcers;
-    std::vector<std::vector<CustomTrianglesPerLayer>> m_blockers;
-    std::vector<const PrintObject*> m_po_list;
-
-    //std::map<const PrintObject*, Point>  m_last_seam_position;
-    SeamHistory  m_seam_history;
-
-    // Get indices of points inside enforcers and blockers.
-    void get_enforcers_and_blockers(size_t layer_id,
-                                    const Polygon& polygon,
-                                    size_t po_id,
-                                    std::vector<size_t>& enforcers_idxs,
-                                    std::vector<size_t>& blockers_idxs) const;
-
-    // Apply penalties to points inside enforcers/blockers.
-    void apply_custom_seam(const Polygon& polygon, size_t po_id,
-                           std::vector<float>& penalties,
-                           const std::vector<float>& lengths,
-                           int layer_id, SeamPosition seam_position) const;
-
-    // Return random point of a polygon. The distribution will be uniform
-    // along the contour and account for enforcers and blockers.
-    Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id,
-                          bool* saw_custom = nullptr) const;
-
-    // Is there any enforcer/blocker on this layer?
-    bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const {
-        return is_custom_enforcer_on_layer(layer_id, po_idx)
-            || is_custom_blocker_on_layer(layer_id, po_idx);
-    }
-
-    bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty());
-    }
-
-    bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty());
-    }
-};
-
 
+namespace EdgeGrid {
+class Grid;
 }
 
+namespace SeamPlacerImpl {
+
+// ************  FOR BACKPORT COMPATIBILITY ONLY ***************
+// Angle from v1 to v2, returning double atan2(y, x) normalized to <-PI, PI>.
+template<typename Derived, typename Derived2> inline double angle(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2)
+{
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "angle(): first parameter is not a 2D vector");
+    static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "angle(): second parameter is not a 2D vector");
+    auto v1d = v1.template cast<double>();
+    auto v2d = v2.template cast<double>();
+    return atan2(cross2(v1d, v2d), v1d.dot(v2d));
+}
+// ***************************
+
+struct GlobalModelInfo;
+struct SeamComparator;
+
+enum class EnforcedBlockedSeamPoint {
+    Blocked  = 0,
+    Neutral  = 1,
+    Enforced = 2,
+};
+
+// struct representing single perimeter loop
+struct Perimeter
+{
+    size_t start_index{};
+    size_t end_index{}; // inclusive!
+    size_t seam_index{};
+    float  flow_width{};
+
+    // During alignment, a final position may be stored here. In that case, finalized is set to true.
+    // Note that final seam position is not limited to points of the perimeter loop. In theory it can be any position
+    // Random position also uses this flexibility to set final seam point position
+    bool  finalized           = false;
+    Vec3f final_seam_position = Vec3f::Zero();
+};
+
+// Struct over which all processing of perimeters is done. For each perimeter point, its respective candidate is created,
+// then all the needed attributes are computed and finally, for each perimeter one point is chosen as seam.
+// This seam position can be then further aligned
+struct SeamCandidate
+{
+    SeamCandidate(const Vec3f &pos, Perimeter &perimeter, float local_ccw_angle, EnforcedBlockedSeamPoint type)
+        : position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), embedded_distance(0.0f), local_ccw_angle(local_ccw_angle), type(type), central_enforcer(false)
+    {}
+    const Vec3f position;
+    // pointer to Perimeter loop of this point. It is shared across all points of the loop
+    Perimeter &perimeter;
+    float      visibility;
+    float      overhang;
+    // distance inside the merged layer regions, for detecting perimeter points which are hidden indside the print (e.g. multimaterial join)
+    // Negative sign means inside the print, comes from EdgeGrid structure
+    float                    embedded_distance;
+    float                    local_ccw_angle;
+    EnforcedBlockedSeamPoint type;
+    bool                     central_enforcer; // marks this candidate as central point of enforced segment on the perimeter - important for alignment
+};
+
+struct SeamCandidateCoordinateFunctor
+{
+    SeamCandidateCoordinateFunctor(const std::vector<SeamCandidate> &seam_candidates) : seam_candidates(seam_candidates) {}
+    const std::vector<SeamCandidate> &seam_candidates;
+    float                             operator()(size_t index, size_t dim) const { return seam_candidates[index].position[dim]; }
+};
+} // namespace SeamPlacerImpl
+
+struct PrintObjectSeamData
+{
+    using SeamCandidatesTree = KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>;
+
+    struct LayerSeams
+    {
+        Slic3r::deque<SeamPlacerImpl::Perimeter>   perimeters;
+        std::vector<SeamPlacerImpl::SeamCandidate> points;
+        std::unique_ptr<SeamCandidatesTree>        points_tree;
+    };
+    // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter
+    std::vector<LayerSeams> layers;
+    // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD
+    // tree of all points of the given layer
+
+    void clear() { layers.clear(); }
+};
+
+class SeamPlacer
+{
+public:
+    // Number of samples generated on the mesh. There are sqr_rays_per_sample_point*sqr_rays_per_sample_point rays casted from each samples
+    static constexpr size_t raycasting_visibility_samples_count = 30000;
+    // square of number of rays per sample point
+    static constexpr size_t sqr_rays_per_sample_point = 5;
+
+    // snapping angle - angles larger than this value will be snapped to during seam painting
+    static constexpr float sharp_angle_snapping_threshold = 55.0f * float(PI) / 180.0f;
+    // overhang angle for seam placement that still yields good results, in degrees, measured from vertical direction
+    //BBS
+    static constexpr float overhang_angle_threshold = 45.0f * float(PI) / 180.0f;
+
+    // determines angle importance compared to visibility ( neutral value is 1.0f. )
+    static constexpr float angle_importance_aligned = 0.6f;
+    static constexpr float angle_importance_nearest = 1.0f; // use much higher angle importance for nearest mode, to combat the visibility info noise
+
+    // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size
+    static constexpr float enforcer_oversampling_distance = 0.2f;
+
+    // When searching for seam clusters for alignment:
+    // following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer
+    static constexpr float seam_align_score_tolerance = 0.3f;
+    // seam_align_tolerable_dist_factor - how far to search for seam from current position, final dist is seam_align_tolerable_dist_factor * flow_width
+    static constexpr float seam_align_tolerable_dist_factor = 4.0f;
+    // minimum number of seams needed in cluster to make alignment happen
+    static constexpr size_t seam_align_minimum_string_seams = 6;
+    // millimeters covered by spline; determines number of splines for the given string
+    static constexpr size_t seam_align_mm_per_segment = 4.0f;
+
+    // The following data structures hold all perimeter points for all PrintObject.
+    std::unordered_map<const PrintObject *, PrintObjectSeamData> m_seam_per_object;
+
+    void init(const Print &print, std::function<void(void)> throw_if_canceled_func);
+
+    void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const;
+
+private:
+    void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info, const SeamPosition configured_seam_preference);
+    void calculate_candidates_visibility(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info);
+    void calculate_overhangs_and_layer_embedding(const PrintObject *po);
+    void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator);
+    std::vector<std::pair<size_t, size_t>> find_seam_string(const PrintObject *po, std::pair<size_t, size_t> start_seam, const SeamPlacerImpl::SeamComparator &comparator) const;
+    std::optional<std::pair<size_t, size_t>> find_next_seam_in_layer(const std::vector<PrintObjectSeamData::LayerSeams> &layers,
+                                                                     const Vec3f &                                       projected_position,
+                                                                     const size_t                                        layer_idx,
+                                                                     const float                                         max_distance,
+                                                                     const SeamPlacerImpl::SeamComparator &              comparator) const;
+};
+
+} // namespace Slic3r
+
 #endif // libslic3r_SeamPlacer_hpp_

src/libslic3r/GCode/ToolOrdering.cpp

@@ -474,7 +474,8 @@ void ToolOrdering::reorder_extruders(std::vector<unsigned int> tool_order_layer0
     if (m_layer_tools.empty())
         return;
 
-    assert(!tool_order_layer0.empty());
+    if (tool_order_layer0.empty())
+        return;
 
     // Reorder the extruders of first layer
     {

src/libslic3r/Geometry/Bicubic.hpp

@@ -0,0 +1,291 @@
+#ifndef BICUBIC_HPP
+#define BICUBIC_HPP
+
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+#include <Eigen/Dense>
+
+namespace Slic3r {
+
+namespace Geometry {
+
+namespace BicubicInternal {
+// Linear kernel, to be able to test cubic methods with hat kernels.
+template<typename T>
+struct LinearKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return T(0.);
+    }
+    static T a01() {
+        return T(0.);
+    }
+    static T a02() {
+        return T(0.);
+    }
+    static T a03() {
+        return T(0.);
+    }
+    static T a10() {
+        return T(1.);
+    }
+    static T a11() {
+        return T(-1.);
+    }
+    static T a12() {
+        return T(0.);
+    }
+    static T a13() {
+        return T(0.);
+    }
+    static T a20() {
+        return T(0.);
+    }
+    static T a21() {
+        return T(1.);
+    }
+    static T a22() {
+        return T(0.);
+    }
+    static T a23() {
+        return T(0.);
+    }
+    static T a30() {
+        return T(0.);
+    }
+    static T a31() {
+        return T(0.);
+    }
+    static T a32() {
+        return T(0.);
+    }
+    static T a33() {
+        return T(0.);
+    }
+};
+
+// Interpolation kernel aka Catmul-Rom aka Keyes kernel.
+template<typename T>
+struct CubicCatmulRomKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return 0;
+    }
+    static T a01() {
+        return T( -0.5);
+    }
+    static T a02() {
+        return T( 1.);
+    }
+    static T a03() {
+        return T( -0.5);
+    }
+    static T a10() {
+        return T( 1.);
+    }
+    static T a11() {
+        return 0;
+    }
+    static T a12() {
+        return T( -5. / 2.);
+    }
+    static T a13() {
+        return T( 3. / 2.);
+    }
+    static T a20() {
+        return 0;
+    }
+    static T a21() {
+        return T( 0.5);
+    }
+    static T a22() {
+        return T( 2.);
+    }
+    static T a23() {
+        return T( -3. / 2.);
+    }
+    static T a30() {
+        return 0;
+    }
+    static T a31() {
+        return 0;
+    }
+    static T a32() {
+        return T( -0.5);
+    }
+    static T a33() {
+        return T( 0.5);
+    }
+};
+
+// B-spline kernel
+template<typename T>
+struct CubicBSplineKernel
+{
+    typedef T FloatType;
+
+    static T a00() {
+        return T( 1. / 6.);
+    }
+    static T a01() {
+        return T( -3. / 6.);
+    }
+    static T a02() {
+        return T( 3. / 6.);
+    }
+    static T a03() {
+        return T( -1. / 6.);
+    }
+    static T a10() {
+        return T( 4. / 6.);
+    }
+    static T a11() {
+        return 0;
+    }
+    static T a12() {
+        return T( -6. / 6.);
+    }
+    static T a13() {
+        return T( 3. / 6.);
+    }
+    static T a20() {
+        return T( 1. / 6.);
+    }
+    static T a21() {
+        return T( 3. / 6.);
+    }
+    static T a22() {
+        return T( 3. / 6.);
+    }
+    static T a23() {
+        return T( -3. / 6.);
+    }
+    static T a30() {
+        return 0;
+    }
+    static T a31() {
+        return 0;
+    }
+    static T a32() {
+        return 0;
+    }
+    static T a33() {
+        return T( 1. / 6.);
+    }
+};
+
+template<class T>
+inline T clamp(T a, T lower, T upper)
+        {
+    return (a < lower) ? lower :
+           (a > upper) ? upper : a;
+}
+}
+
+template<typename Kernel>
+struct CubicKernelWrapper
+{
+    typedef typename Kernel::FloatType FloatType;
+
+    static constexpr size_t kernel_span = 4;
+
+    static FloatType kernel(FloatType x)
+            {
+        x = fabs(x);
+        if (x >= (FloatType) 2.)
+            return 0.0f;
+        if (x <= (FloatType) 1.) {
+            FloatType x2 = x * x;
+            FloatType x3 = x2 * x;
+            return Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3;
+        }
+        assert(x > (FloatType )1. && x < (FloatType )2.);
+        x -= (FloatType) 1.;
+        FloatType x2 = x * x;
+        FloatType x3 = x2 * x;
+        return Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3;
+    }
+
+    static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x)
+            {
+        const FloatType x2 = x * x;
+        const FloatType x3 = x * x * x;
+        return f0 * (Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3) +
+                f1 * (Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3) +
+                f2 * (Kernel::a20() + Kernel::a21() * x + Kernel::a22() * x2 + Kernel::a23() * x3) +
+                f3 * (Kernel::a30() + Kernel::a31() * x + Kernel::a32() * x2 + Kernel::a33() * x3);
+    }
+};
+
+// Linear splines
+template<typename NumberType>
+using LinearKernel = CubicKernelWrapper<BicubicInternal::LinearKernel<NumberType>>;
+
+// Catmul-Rom splines
+template<typename NumberType>
+using CubicCatmulRomKernel = CubicKernelWrapper<BicubicInternal::CubicCatmulRomKernel<NumberType>>;
+
+// Cubic B-splines
+template<typename NumberType>
+using CubicBSplineKernel = CubicKernelWrapper<BicubicInternal::CubicBSplineKernel<NumberType>>;
+
+template<typename KernelWrapper>
+static typename KernelWrapper::FloatType cubic_interpolate(const Eigen::ArrayBase<typename KernelWrapper::FloatType> &F,
+        const typename KernelWrapper::FloatType pt) {
+    typedef typename KernelWrapper::FloatType T;
+    const int w = int(F.size());
+    const int ix = (int) floor(pt);
+    const T s = pt - T( ix);
+
+    if (ix > 1 && ix + 2 < w) {
+        // Inside the fully interpolated region.
+        return KernelWrapper::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
+    }
+    // Transition region. Extend with a constant function.
+    auto f = [&F, w](T x) {
+        return F[BicubicInternal::clamp(x, 0, w - 1)];
+    };
+    return KernelWrapper::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s);
+}
+
+template<typename Kernel, typename Derived>
+static float bicubic_interpolate(const Eigen::MatrixBase<Derived> &F,
+        const Eigen::Matrix<typename Kernel::FloatType, 2, 1, Eigen::DontAlign> &pt) {
+    typedef typename Kernel::FloatType T;
+    const int w = F.cols();
+    const int h = F.rows();
+    const int ix = (int) floor(pt[0]);
+    const int iy = (int) floor(pt[1]);
+    const T s = pt[0] - T( ix);
+    const T t = pt[1] - T( iy);
+
+    if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) {
+        // Inside the fully interpolated region.
+        return Kernel::interpolate(
+                Kernel::interpolate(F(ix - 1, iy - 1), F(ix, iy - 1), F(ix + 1, iy - 1), F(ix + 2, iy - 1), s),
+                Kernel::interpolate(F(ix - 1, iy), F(ix, iy), F(ix + 1, iy), F(ix + 2, iy), s),
+                Kernel::interpolate(F(ix - 1, iy + 1), F(ix, iy + 1), F(ix + 1, iy + 1), F(ix + 2, iy + 1), s),
+                Kernel::interpolate(F(ix - 1, iy + 2), F(ix, iy + 2), F(ix + 1, iy + 2), F(ix + 2, iy + 2), s), t);
+    }
+    // Transition region. Extend with a constant function.
+    auto f = [&F, w, h](int x, int y) {
+        return F(BicubicInternal::clamp(x, 0, w - 1), BicubicInternal::clamp(y, 0, h - 1));
+    };
+    return Kernel::interpolate(
+            Kernel::interpolate(f(ix - 1, iy - 1), f(ix, iy - 1), f(ix + 1, iy - 1), f(ix + 2, iy - 1), s),
+            Kernel::interpolate(f(ix - 1, iy), f(ix, iy), f(ix + 1, iy), f(ix + 2, iy), s),
+            Kernel::interpolate(f(ix - 1, iy + 1), f(ix, iy + 1), f(ix + 1, iy + 1), f(ix + 2, iy + 1), s),
+            Kernel::interpolate(f(ix - 1, iy + 2), f(ix, iy + 2), f(ix + 1, iy + 2), f(ix + 2, iy + 2), s), t);
+}
+
+} //namespace Geometry
+
+} // namespace Slic3r
+
+#endif /* BICUBIC_HPP */

src/libslic3r/Geometry/Curves.hpp

@@ -0,0 +1,218 @@
+#ifndef SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
+#define SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
+
+#include "libslic3r/Point.hpp"
+#include "Bicubic.hpp"
+
+#include <iostream>
+
+//#define LSQR_DEBUG
+
+namespace Slic3r {
+namespace Geometry {
+
+template<int Dimension, typename NumberType>
+struct PolynomialCurve {
+    Eigen::MatrixXf coefficients;
+
+    Vec<Dimension, NumberType> get_fitted_value(const NumberType& value) const {
+        Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
+        size_t order = this->coefficients.rows() - 1;
+        auto x = NumberType(1.);
+        for (size_t index = 0; index < order + 1; ++index, x *= value)
+            result += x * this->coefficients.col(index);
+        return result;
+    }
+};
+
+//https://towardsdatascience.com/least-square-polynomial-CURVES-using-c-eigen-package-c0673728bd01
+template<int Dimension, typename NumberType>
+PolynomialCurve<Dimension, NumberType> fit_polynomial(const std::vector<Vec<Dimension, NumberType>> &observations,
+        const std::vector<NumberType> &observation_points,
+        const std::vector<NumberType> &weights, size_t order) {
+    // check to make sure inputs are correct
+    size_t cols = order + 1;
+    assert(observation_points.size() >= cols);
+    assert(observation_points.size() == weights.size());
+    assert(observations.size() == weights.size());
+
+    Eigen::MatrixXf data_points(Dimension, observations.size());
+    Eigen::MatrixXf T(observations.size(), cols);
+    for (size_t i = 0; i < weights.size(); ++i) {
+        auto squared_weight = sqrt(weights[i]);
+        data_points.col(i) = observations[i] * squared_weight;
+        // Populate the matrix
+        auto x = squared_weight;
+        auto c = observation_points[i];
+        for (size_t j = 0; j < cols; ++j, x *= c)
+            T(i, j) = x;
+    }
+
+    const auto QR = T.householderQr();
+    Eigen::MatrixXf coefficients(Dimension, cols);
+    // Solve for linear least square fit
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+        coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
+    }
+
+    return {std::move(coefficients)};
+}
+
+template<size_t Dimension, typename NumberType, typename KernelType>
+struct PiecewiseFittedCurve {
+    using Kernel = KernelType;
+
+    Eigen::MatrixXf coefficients;
+    NumberType start;
+    NumberType segment_size;
+    size_t endpoints_level_of_freedom;
+
+    Vec<Dimension, NumberType> get_fitted_value(const NumberType &observation_point) const {
+        Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
+
+        //find corresponding segment index; expects kernels to be centered
+        int middle_right_segment_index = floor((observation_point - start) / segment_size);
+        //find index of first segment that is affected by the point i; this can be deduced from kernel_span
+        int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1;
+        for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
+                segment_index++) {
+            NumberType segment_start = start + segment_index * segment_size;
+            NumberType normalized_segment_distance = (segment_start - observation_point) / segment_size;
+
+            int parameter_index = segment_index + endpoints_level_of_freedom;
+            parameter_index = std::clamp(parameter_index, 0, int(coefficients.cols()) - 1);
+            result += Kernel::kernel(normalized_segment_distance) * coefficients.col(parameter_index);
+        }
+        return result;
+    }
+};
+
+// observations: data to be fitted by the curve
+// observation points: growing sequence of points where the observations were made.
+//      In other words, for function f(x) = y, observations are y0...yn, and observation points are x0...xn
+// weights: how important the observation is
+// segments_count: number of segments inside the valid length of the curve
+// endpoints_level_of_freedom: number of additional parameters at each end; reasonable values depend on the kernel span
+template<typename Kernel, int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, Kernel> fit_curve(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        const std::vector<NumberType> &observation_points,
+        const std::vector<NumberType> &weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom) {
+
+    // check to make sure inputs are correct
+    assert(segments_count > 0);
+    assert(observations.size() > 0);
+    assert(observation_points.size() == observations.size());
+    assert(observation_points.size() == weights.size());
+    assert(segments_count <= observations.size());
+
+    //prepare sqrt of weights, which will then be applied to both matrix T and observed data: https://en.wikipedia.org/wiki/Weighted_least_squares
+    std::vector<NumberType> sqrt_weights(weights.size());
+    for (size_t index = 0; index < weights.size(); ++index) {
+        assert(weights[index] > 0);
+        sqrt_weights[index] = sqrt(weights[index]);
+    }
+
+    // prepare result and compute metadata
+    PiecewiseFittedCurve<Dimension, NumberType, Kernel> result { };
+
+    NumberType valid_length = observation_points.back() - observation_points.front();
+    NumberType segment_size = valid_length / NumberType(segments_count);
+    result.start = observation_points.front();
+    result.segment_size = segment_size;
+    result.endpoints_level_of_freedom = endpoints_level_of_freedom;
+
+    // prepare observed data
+    // Eigen defaults to column major memory layout.
+    Eigen::MatrixXf data_points(Dimension, observations.size());
+    for (size_t index = 0; index < observations.size(); ++index) {
+        data_points.col(index) = observations[index] * sqrt_weights[index];
+    }
+    // parameters count is always increased by one to make the parametric space of the curve symmetric.
+    // without this fix, the end of the curve is less flexible than the beginning
+    size_t parameters_count = segments_count + 1 + 2 * endpoints_level_of_freedom;
+    //Create weight matrix T for each point and each segment;
+    Eigen::MatrixXf T(observation_points.size(), parameters_count);
+    T.setZero();
+    //Fill the weight matrix
+    for (size_t i = 0; i < observation_points.size(); ++i) {
+        NumberType observation_point = observation_points[i];
+        //find corresponding segment index; expects kernels to be centered
+        int middle_right_segment_index = floor((observation_point - result.start) / result.segment_size);
+        //find index of first segment that is affected by the point i; this can be deduced from kernel_span
+        int start_segment_idx = middle_right_segment_index - int(Kernel::kernel_span / 2) + 1;
+        for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
+                segment_index++) {
+            NumberType segment_start = result.start + segment_index * result.segment_size;
+            NumberType normalized_segment_distance = (segment_start - observation_point) / result.segment_size;
+
+            int parameter_index = segment_index + endpoints_level_of_freedom;
+            parameter_index = std::clamp(parameter_index, 0, int(parameters_count) - 1);
+            T(i, parameter_index) += Kernel::kernel(normalized_segment_distance) * sqrt_weights[i];
+        }
+    }
+
+#ifdef LSQR_DEBUG
+    std::cout << "weight matrix: " << std::endl;
+    for (int obs = 0; obs < observation_points.size(); ++obs) {
+        std::cout << std::endl;
+        for (int segment = 0; segment < parameters_count; ++segment) {
+            std::cout << T(obs, segment) << "  ";
+        }
+    }
+    std::cout << std::endl;
+#endif
+
+    // Solve for linear least square fit
+    result.coefficients.resize(Dimension, parameters_count);
+    const auto QR = T.fullPivHouseholderQr();
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+        result.coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
+    }
+
+    return result;
+}
+
+
+template<int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, LinearKernel<NumberType>>
+fit_linear_spline(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        std::vector<NumberType> observation_points,
+        std::vector<NumberType> weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom = 0) {
+    return fit_curve<LinearKernel<NumberType>>(observations, observation_points, weights, segments_count,
+            endpoints_level_of_freedom);
+}
+
+template<int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, CubicBSplineKernel<NumberType>>
+fit_cubic_bspline(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        std::vector<NumberType> observation_points,
+        std::vector<NumberType> weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom = 0) {
+    return fit_curve<CubicBSplineKernel<NumberType>>(observations, observation_points, weights, segments_count,
+            endpoints_level_of_freedom);
+}
+
+template<int Dimension, typename NumberType>
+PiecewiseFittedCurve<Dimension, NumberType, CubicCatmulRomKernel<NumberType>>
+fit_catmul_rom_spline(
+        const std::vector<Vec<Dimension, NumberType>> &observations,
+        std::vector<NumberType> observation_points,
+        std::vector<NumberType> weights,
+        size_t segments_count,
+        size_t endpoints_level_of_freedom = 0) {
+    return fit_curve<CubicCatmulRomKernel<NumberType>>(observations, observation_points, weights, segments_count,
+            endpoints_level_of_freedom);
+}
+
+}
+}
+
+#endif /* SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ */

src/libslic3r/KDTreeIndirect.hpp

@@ -11,224 +11,362 @@
 
 namespace Slic3r {
 
+enum class VisitorReturnMask : unsigned int {
+    CONTINUE_LEFT  = 1,
+    CONTINUE_RIGHT = 2,
+    STOP           = 4,
+};
+
 // KD tree for N-dimensional closest point search.
 template<size_t ANumDimensions, typename ACoordType, typename ACoordinateFn>
 class KDTreeIndirect
 {
 public:
-	static constexpr size_t NumDimensions = ANumDimensions;
-	using					CoordinateFn  = ACoordinateFn;
-	using					CoordType     = ACoordType;
+    static constexpr size_t NumDimensions = ANumDimensions;
+    using					CoordinateFn  = ACoordinateFn;
+    using					CoordType     = ACoordType;
     // Following could be static constexpr size_t, but that would not link in C++11
     enum : size_t {
         npos = size_t(-1)
     };
 
-	KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {}
-	KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t>   indices) : coordinate(coordinate) { this->build(std::move(indices)); }
-	KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> &&indices) : coordinate(coordinate) { this->build(std::move(indices)); }
-	KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); }
-	KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {}
-	KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; }
-	void clear() { m_nodes.clear(); }
+    KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {}
+    KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> indices) : coordinate(coordinate) { this->build(indices); }
+    KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); }
+    KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {}
+    KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; }
+    void clear() { m_nodes.clear(); }
 
-	void build(size_t num_indices)
-	{
-		std::vector<size_t> indices;
-		indices.reserve(num_indices);
-		for (size_t i = 0; i < num_indices; ++ i)
-			indices.emplace_back(i);
-		this->build(std::move(indices));
-	}
+    void build(size_t num_indices)
+    {
+        std::vector<size_t> indices;
+        indices.reserve(num_indices);
+        for (size_t i = 0; i < num_indices; ++ i)
+            indices.emplace_back(i);
+        this->build(indices);
+    }
 
-	void build(std::vector<size_t> &&indices)
-	{
-		if (indices.empty())
-			clear();
-		else {
-			// Allocate enough memory for a full binary tree.
-			m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos);
-			build_recursive(indices, 0, 0, 0, indices.size() - 1);
-		}
-		indices.clear();
-	}
+    void build(std::vector<size_t> &indices)
+    {
+        if (indices.empty())
+            clear();
+        else {
+            // Allocate enough memory for a full binary tree.
+            m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos);
+            build_recursive(indices, 0, 0, 0, indices.size() - 1);
+        }
+        indices.clear();
+    }
 
-	enum class VisitorReturnMask : unsigned int
-	{
-		CONTINUE_LEFT   = 1,
-		CONTINUE_RIGHT  = 2,
-		STOP 			= 4,
-	};
-	template<typename CoordType> 
-	unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
-	{
-		CoordType dist = point_coord - this->coordinate(idx, dimension);
-		return (dist * dist < search_radius + CoordType(EPSILON)) ?
-			// The plane intersects a hypersphere centered at point_coord of search_radius.
-			((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) :
-			// The plane does not intersect the hypersphere.
-			(dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT);
-	}
+    template<typename CoordType>
+    unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
+    {
+        CoordType dist = point_coord - this->coordinate(idx, dimension);
+        return (dist * dist < search_radius + CoordType(EPSILON)) ?
+                                                                    // The plane intersects a hypersphere centered at point_coord of search_radius.
+                   ((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) :
+                   // The plane does not intersect the hypersphere.
+                   (dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT);
+    }
 
-	// Visitor is supposed to return a bit mask of VisitorReturnMask.
-	template<typename Visitor>
-	void visit(Visitor &visitor) const
-	{
+       // Visitor is supposed to return a bit mask of VisitorReturnMask.
+    template<typename Visitor>
+    void visit(Visitor &visitor) const
+    {
         visit_recursive(0, 0, visitor);
-	}
+    }
 
-	CoordinateFn coordinate;
+    CoordinateFn coordinate;
 
 private:
-	// Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension.
-	void build_recursive(std::vector<size_t> &input, size_t node, const size_t dimension, const size_t left, const size_t right)
-	{
-		if (left > right)
-			return;
+    // Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension.
+    void build_recursive(std::vector<size_t> &input, size_t node, const size_t dimension, const size_t left, const size_t right)
+    {
+        if (left > right)
+            return;
 
-		assert(node < m_nodes.size());
+        assert(node < m_nodes.size());
 
-		if (left == right) {
-			// Insert a node into the balanced tree.
-			m_nodes[node] = input[left];
-			return;
-		}
+        if (left == right) {
+            // Insert a node into the balanced tree.
+            m_nodes[node] = input[left];
+            return;
+        }
 
-		// Partition the input to left / right pieces of the same length to produce a balanced tree.
-		size_t center = (left + right) / 2;
-		partition_input(input, dimension, left, right, center);
-		// Insert a node into the tree.
-		m_nodes[node] = input[center];
-		// Build up the left / right subtrees.
-		size_t next_dimension = dimension;
-		if (++ next_dimension == NumDimensions)
-			next_dimension = 0;
-		if (center > left)
-			build_recursive(input, node * 2 + 1, next_dimension, left, center - 1);
-		build_recursive(input, node * 2 + 2, next_dimension, center + 1, right);
-	}
+           // Partition the input to left / right pieces of the same length to produce a balanced tree.
+        size_t center = (left + right) / 2;
+        partition_input(input, dimension, left, right, center);
+        // Insert a node into the tree.
+        m_nodes[node] = input[center];
+        // Build up the left / right subtrees.
+        size_t next_dimension = dimension;
+        if (++ next_dimension == NumDimensions)
+            next_dimension = 0;
+        if (center > left)
+            build_recursive(input, node * 2 + 1, next_dimension, left, center - 1);
+        build_recursive(input, node * 2 + 2, next_dimension, center + 1, right);
+    }
 
-	// Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
-	// https://en.wikipedia.org/wiki/Quickselect
-	// Items left of the k'th item are lower than the k'th item in the "dimension", 
-	// items right of the k'th item are higher than the k'th item in the "dimension", 
-	void partition_input(std::vector<size_t> &input, const size_t dimension, size_t left, size_t right, const size_t k) const
-	{
-		while (left < right) {
-			size_t center = (left + right) / 2;
-			CoordType pivot;
-			{
-				// Bubble sort the input[left], input[center], input[right], so that a median of the three values
-				// will end up in input[center].
-				CoordType left_value   = this->coordinate(input[left],   dimension);
-				CoordType center_value = this->coordinate(input[center], dimension);
-				CoordType right_value  = this->coordinate(input[right],  dimension);
-				if (left_value > center_value) {
-					std::swap(input[left], input[center]);
-					std::swap(left_value,  center_value);
-				}
-				if (left_value > right_value) {
-					std::swap(input[left], input[right]);
-					right_value = left_value;
-				}
-				if (center_value > right_value) {
-					std::swap(input[center], input[right]);
-					center_value = right_value;
-				}
-				pivot = center_value;
-			}
-			if (right <= left + 2)
-				// The <left, right> interval is already sorted.
-				break;
-			size_t i = left;
-			size_t j = right - 1;
-			std::swap(input[center], input[j]);
-			// Partition the set based on the pivot.
-			for (;;) {
-				// Skip left points that are already at correct positions.
-				// Search will certainly stop at position (right - 1), which stores the pivot.
-				while (this->coordinate(input[++ i], dimension) < pivot) ;
-				// Skip right points that are already at correct positions.
-				while (this->coordinate(input[-- j], dimension) > pivot && i < j) ;
-				if (i >= j)
-					break;
-				std::swap(input[i], input[j]);
-			}
-			// Restore pivot to the center of the sequence.
-			std::swap(input[i], input[right - 1]);
-			// Which side the kth element is in?
-			if (k < i)
-				right = i - 1;
-			else if (k == i)
-				// Sequence is partitioned, kth element is at its place.
-				break;
-			else
-				left = i + 1;
-		}
-	}
+       // Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
+       // https://en.wikipedia.org/wiki/Quickselect
+       // Items left of the k'th item are lower than the k'th item in the "dimension",
+       // items right of the k'th item are higher than the k'th item in the "dimension",
+    void partition_input(std::vector<size_t> &input, const size_t dimension, size_t left, size_t right, const size_t k) const
+    {
+        while (left < right) {
+            size_t center = (left + right) / 2;
+            CoordType pivot;
+            {
+                // Bubble sort the input[left], input[center], input[right], so that a median of the three values
+                // will end up in input[center].
+                CoordType left_value   = this->coordinate(input[left],   dimension);
+                CoordType center_value = this->coordinate(input[center], dimension);
+                CoordType right_value  = this->coordinate(input[right],  dimension);
+                if (left_value > center_value) {
+                    std::swap(input[left], input[center]);
+                    std::swap(left_value,  center_value);
+                }
+                if (left_value > right_value) {
+                    std::swap(input[left], input[right]);
+                    right_value = left_value;
+                }
+                if (center_value > right_value) {
+                    std::swap(input[center], input[right]);
+                    center_value = right_value;
+                }
+                pivot = center_value;
+            }
+            if (right <= left + 2)
+                // The <left, right> interval is already sorted.
+                break;
+            size_t i = left;
+            size_t j = right - 1;
+            std::swap(input[center], input[j]);
+            // Partition the set based on the pivot.
+            for (;;) {
+                // Skip left points that are already at correct positions.
+                // Search will certainly stop at position (right - 1), which stores the pivot.
+                while (this->coordinate(input[++ i], dimension) < pivot) ;
+                // Skip right points that are already at correct positions.
+                while (this->coordinate(input[-- j], dimension) > pivot && i < j) ;
+                if (i >= j)
+                    break;
+                std::swap(input[i], input[j]);
+            }
+            // Restore pivot to the center of the sequence.
+            std::swap(input[i], input[right - 1]);
+            // Which side the kth element is in?
+            if (k < i)
+                right = i - 1;
+            else if (k == i)
+                // Sequence is partitioned, kth element is at its place.
+                break;
+            else
+                left = i + 1;
+        }
+    }
 
-	template<typename Visitor>
-	void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const
-	{
-		assert(! m_nodes.empty());
-		if (node >= m_nodes.size() || m_nodes[node] == npos)
-			return;
+    template<typename Visitor>
+    void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const
+    {
+        assert(! m_nodes.empty());
+        if (node >= m_nodes.size() || m_nodes[node] == npos)
+            return;
 
-		// Left / right child node index.
-		size_t left  = node * 2 + 1;
-		size_t right = left + 1;
-		unsigned int mask = visitor(m_nodes[node], dimension);
-		if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) {
-			size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension;
-			if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT)
-				visit_recursive(left,  next_dimension, visitor);
-			if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT)
-				visit_recursive(right, next_dimension, visitor);
-		}
-	}
+           // Left / right child node index.
+        size_t left  = node * 2 + 1;
+        size_t right = left + 1;
+        unsigned int mask = visitor(m_nodes[node], dimension);
+        if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) {
+            size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension;
+            if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT)
+                visit_recursive(left,  next_dimension, visitor);
+            if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT)
+                visit_recursive(right, next_dimension, visitor);
+        }
+    }
 
-	std::vector<size_t> m_nodes;
+    std::vector<size_t> m_nodes;
 };
 
 // Find a closest point using Euclidian metrics.
 // Returns npos if not found.
-template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
-size_t find_closest_point(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter)
+template<size_t K,
+         typename PointType,
+         typename FilterFn,
+         size_t D,
+         typename CoordT,
+         typename CoordFn>
+std::array<size_t, K> find_closest_points(
+    const KDTreeIndirect<D, CoordT, CoordFn> &kdtree,
+    const PointType                          &point,
+    FilterFn                                  filter)
 {
-	using CoordType = typename KDTreeIndirectType::CoordType;
+    using Tree = KDTreeIndirect<D, CoordT, CoordFn>;
 
-	struct Visitor {
-		const KDTreeIndirectType   &kdtree;
-		const PointType    		   &point;
-		const FilterFn				filter;
-		size_t 						min_idx  = KDTreeIndirectType::npos;
-		CoordType					min_dist = std::numeric_limits<CoordType>::max();
+    struct Visitor
+    {
+        const Tree      &kdtree;
+        const PointType &point;
+        const FilterFn   filter;
 
-		Visitor(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) : kdtree(kdtree), point(point), filter(filter) {}
-		unsigned int operator()(size_t idx, size_t dimension) {
-			if (this->filter(idx)) {
-				auto dist = CoordType(0);
-				for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++ i) {
-					CoordType d = point[i] - kdtree.coordinate(idx, i);
-					dist += d * d;
-				}
-				if (dist < min_dist) {
-					min_dist = dist;
-					min_idx  = idx;
-				}
-			}
-			return kdtree.descent_mask(point[dimension], min_dist, idx, dimension);
-		}
-	} visitor(kdtree, point, filter);
+        std::array<std::pair<size_t, CoordT>, K> results;
 
-	kdtree.visit(visitor);
-	return visitor.min_idx;
+        Visitor(const Tree &kdtree, const PointType &point, FilterFn filter)
+            : kdtree(kdtree), point(point), filter(filter)
+        {
+            results.fill(std::make_pair(Tree::npos,
+                                        std::numeric_limits<CoordT>::max()));
+        }
+        unsigned int operator()(size_t idx, size_t dimension)
+        {
+            if (this->filter(idx)) {
+                auto dist = CoordT(0);
+                for (size_t i = 0; i < D; ++i) {
+                    CoordT d = point[i] - kdtree.coordinate(idx, i);
+                    dist += d * d;
+                }
+
+                auto res = std::make_pair(idx, dist);
+                auto it  = std::lower_bound(results.begin(), results.end(),
+                                            res, [](auto &r1, auto &r2) {
+                                               return r1.second < r2.second;
+                                            });
+
+                if (it != results.end()) {
+                    std::rotate(it, std::prev(results.end()), results.end());
+                    *it = res;
+                }
+            }
+            return kdtree.descent_mask(point[dimension],
+                                       results.front().second, idx,
+                                       dimension);
+        }
+    } visitor(kdtree, point, filter);
+
+    kdtree.visit(visitor);
+    std::array<size_t, K> ret;
+    for (size_t i = 0; i < K; i++) ret[i] = visitor.results[i].first;
+
+    return ret;
+}
+
+template<size_t K, typename PointType, size_t D, typename CoordT, typename CoordFn>
+std::array<size_t, K> find_closest_points(
+    const KDTreeIndirect<D, CoordT, CoordFn> &kdtree, const PointType &point)
+{
+    return find_closest_points<K>(kdtree, point, [](size_t) { return true; });
+}
+
+template<typename PointType,
+         typename FilterFn,
+         size_t D,
+         typename CoordT,
+         typename CoordFn>
+size_t find_closest_point(const KDTreeIndirect<D, CoordT, CoordFn> &kdtree,
+                          const PointType                          &point,
+                          FilterFn                                  filter)
+{
+    return find_closest_points<1>(kdtree, point, filter)[0];
 }
 
 template<typename KDTreeIndirectType, typename PointType>
 size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& point)
 {
-	return find_closest_point(kdtree, point, [](size_t) { return true; });
+    return find_closest_point(kdtree, point, [](size_t) { return true; });
+}
+
+// Find nearby points (spherical neighbourhood) using Euclidian metrics.
+template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+                                       const typename KDTreeIndirectType::CoordType& max_distance, FilterFn filter)
+{
+    using CoordType = typename KDTreeIndirectType::CoordType;
+
+    struct Visitor {
+        const KDTreeIndirectType &kdtree;
+        const PointType center;
+        const CoordType max_distance_squared;
+        const FilterFn filter;
+        std::vector<size_t> result;
+
+        Visitor(const KDTreeIndirectType &kdtree, const PointType& center, const CoordType &max_distance,
+                FilterFn filter) :
+            kdtree(kdtree), center(center), max_distance_squared(max_distance*max_distance), filter(filter) {
+        }
+        unsigned int operator()(size_t idx, size_t dimension) {
+            if (this->filter(idx)) {
+                auto dist = CoordType(0);
+                for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
+                    CoordType d = center[i] - kdtree.coordinate(idx, i);
+                    dist += d * d;
+                }
+                if (dist < max_distance_squared) {
+                    result.push_back(idx);
+                }
+            }
+            return kdtree.descent_mask(center[dimension], max_distance_squared, idx, dimension);
+        }
+    } visitor(kdtree, center, max_distance, filter);
+
+    kdtree.visit(visitor);
+    return visitor.result;
+}
+
+template<typename KDTreeIndirectType, typename PointType>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+                                       const typename KDTreeIndirectType::CoordType& max_distance)
+{
+    return find_nearby_points(kdtree, center, max_distance, [](size_t) {
+        return true;
+    });
+}
+
+// Find nearby points (spherical neighbourhood) using Euclidian metrics.
+template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree,
+                                       const PointType          &bb_min,
+                                       const PointType          &bb_max,
+                                       FilterFn                  filter)
+{
+    struct Visitor {
+        const KDTreeIndirectType &kdtree;
+        const PointType &bb_min, &bb_max;
+        const FilterFn filter;
+        std::vector<size_t> result;
+
+        Visitor(const KDTreeIndirectType &kdtree, const PointType& bbmin, const PointType& bbmax,
+                FilterFn filter) :
+            kdtree(kdtree), bb_min{bbmin}, bb_max{bbmax}, filter(filter) {
+        }
+        unsigned int operator()(size_t idx, size_t dimension) {
+            unsigned int ret =
+                static_cast<unsigned int>(VisitorReturnMask::CONTINUE_LEFT) |
+                static_cast<unsigned int>(VisitorReturnMask::CONTINUE_RIGHT);
+
+            if (this->filter(idx)) {
+                PointType p;
+                bool contains = true;
+                for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
+                    p(i) = kdtree.coordinate(idx, i);
+                    contains = contains && bb_min(i) <= p(i) && p(i) <= bb_max(i);
+                }
+
+                if (p(dimension) < bb_min(dimension))
+                    ret = static_cast<unsigned int>(VisitorReturnMask::CONTINUE_RIGHT);
+                if (p(dimension) > bb_max(dimension))
+                    ret = static_cast<unsigned int>(VisitorReturnMask::CONTINUE_LEFT);
+
+                if (contains)
+                    result.emplace_back(idx);
+            }
+
+            return ret;
+        }
+    } visitor(kdtree, bb_min, bb_max, filter);
+
+    kdtree.visit(visitor);
+    return visitor.result;
 }
 
 } // namespace Slic3r

src/libslic3r/Model.cpp

@@ -137,7 +137,7 @@ Model::~Model()
 // Loading model from a file, it may be a simple geometry file as STL or OBJ, however it may be a project file as well.
 Model Model::read_from_file(const std::string& input_file, DynamicPrintConfig* config, ConfigSubstitutionContext* config_substitutions,
                             LoadStrategy options, PlateDataPtrs* plate_data, std::vector<Preset*>* project_presets, bool *is_xxx, Semver* file_version, Import3mfProgressFn proFn,
-                            ImportStepProgressFn stepFn, StepIsUtf8Fn stepIsUtf8Fn, BBLProject* project)
+                            ImportstlProgressFn stlFn, ImportStepProgressFn stepFn, StepIsUtf8Fn stepIsUtf8Fn, BBLProject* project)
 {
     Model model;
 
@@ -163,7 +163,7 @@ Model Model::read_from_file(const std::string& input_file, DynamicPrintConfig* c
         boost::algorithm::iends_with(input_file, ".step"))
         result = load_step(input_file.c_str(), &model, stepFn, stepIsUtf8Fn);
     else if (boost::algorithm::iends_with(input_file, ".stl"))
-        result = load_stl(input_file.c_str(), &model);
+        result = load_stl(input_file.c_str(), &model, nullptr, stlFn);
     else if (boost::algorithm::iends_with(input_file, ".obj"))
         result = load_obj(input_file.c_str(), &model);
     //BBS: remove the old .amf.xml files
@@ -185,10 +185,10 @@ Model Model::read_from_file(const std::string& input_file, DynamicPrintConfig* c
 
     if (model.objects.empty())
         throw Slic3r::RuntimeError("The supplied file couldn't be read because it's empty");
-    
+
     for (ModelObject *o : model.objects)
         o->input_file = input_file;
-    
+
     if (options & LoadStrategy::AddDefaultInstances)
         model.add_default_instances();
 
@@ -259,14 +259,14 @@ Model Model::read_from_archive(const std::string& input_file, DynamicPrintConfig
 
     //BBS
     //CustomGCode::update_custom_gcode_per_print_z_from_config(model.custom_gcode_per_print_z, config);
-    
+
     BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << ":" << __LINE__ << boost::format("import 3mf IMPORT_STAGE_UPDATE_GCODE\n");
     if (proFn) {
         proFn(IMPORT_STAGE_UPDATE_GCODE, 0, 1, cb_cancel);
         if (cb_cancel)
             throw Slic3r::RuntimeError("Canceled");
     }
-    
+
     CustomGCode::check_mode_for_custom_gcode_per_print_z(model.custom_gcode_per_print_z);
 
     BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << ":" << __LINE__ << boost::format("import 3mf IMPORT_STAGE_CHECK_MODE_GCODE\n");
@@ -416,7 +416,7 @@ void Model::collect_reusable_objects(std::vector<ObjectBase*>& objects)
                        std::mem_fn(&ObjectBase::id));
         model_object->volumes.clear();
     }
-    // we never own these objects 
+    // we never own these objects
     this->objects.clear();
 }
 
@@ -591,7 +591,7 @@ void Model::convert_multipart_object(unsigned int max_extruders)
     assert(this->objects.size() >= 2);
     if (this->objects.size() < 2)
         return;
-    
+
     ModelObject* object = new ModelObject(this);
     object->input_file = this->objects.front()->input_file;
     object->name = boost::filesystem::path(this->objects.front()->input_file).stem().string();
@@ -601,7 +601,7 @@ void Model::convert_multipart_object(unsigned int max_extruders)
 
 	for (const ModelObject* o : this->objects)
     	for (const ModelVolume* v : o->volumes) {
-            // If there are more than one object, put all volumes together 
+            // If there are more than one object, put all volumes together
             // Each object may contain any number of volumes and instances
             // The volumes transformations are relative to the object containing them...
             Geometry::Transformation trafo_volume = v->get_transformation();
@@ -620,7 +620,7 @@ void Model::convert_multipart_object(unsigned int max_extruders)
             } else {
                 for (const ModelInstance* i : o->instances)
                     // ...so, transform everything to a common reference system (world)
-                	copy_volume(object->add_volume(*v))->set_transformation(i->get_transformation() * trafo_volume);                    
+                    copy_volume(object->add_volume(*v))->set_transformation(i->get_transformation() * trafo_volume);
             }
         }
 
@@ -1010,6 +1010,20 @@ ModelVolume* ModelObject::add_volume(const ModelVolume &other, TriangleMesh &&me
     return v;
 }
 
+ModelVolume* ModelObject::add_volume_with_shared_mesh(const ModelVolume &other, ModelVolumeType type /*= ModelVolumeType::INVALID*/)
+{
+    ModelVolume* v = new ModelVolume(this, other.m_mesh);
+    if (type != ModelVolumeType::INVALID && v->type() != type)
+        v->set_type(type);
+    this->volumes.push_back(v);
+	// The volume should already be centered at this point of time when copying shared pointers of the triangle mesh and convex hull.
+//	v->center_geometry_after_creation();
+//    this->invalidate_bounding_box();
+    // BBS: backup
+    Slic3r::save_object_mesh(*this);
+    return v;
+}
+
 void ModelObject::delete_volume(size_t idx)
 {
     ModelVolumePtrs::iterator i = this->volumes.begin() + idx;
@@ -1324,7 +1338,7 @@ Polygon ModelObject::convex_hull_2d(const Transform3d& trafo_instance) const
 
 void ModelObject::center_around_origin(bool include_modifiers)
 {
-    // calculate the displacements needed to 
+    // calculate the displacements needed to
     // center this object around the origin
     const BoundingBoxf3 bb = include_modifiers ? full_raw_mesh_bounding_box() : raw_mesh_bounding_box();
 
@@ -1476,8 +1490,8 @@ void ModelObject::convert_units(ModelObjectPtrs& new_objects, ConversionType con
 
             // Perform conversion only if the target "imperial" state is different from the current one.
             // This check supports conversion of "mixed" set of volumes, each with different "imperial" state.
-            if (//vol->source.is_converted_from_inches != from_imperial && 
-                (volume_idxs.empty() || 
+            if (//vol->source.is_converted_from_inches != from_imperial &&
+                (volume_idxs.empty() ||
                  std::find(volume_idxs.begin(), volume_idxs.end(), vol_idx) != volume_idxs.end())) {
                 vol->scale_geometry_after_creation(koef);
                 vol->set_offset(Vec3d(koef, koef, koef).cwiseProduct(volume->get_offset()));
@@ -1538,7 +1552,7 @@ ModelObjectPtrs ModelObject::cut(size_t instance, std::array<Vec3d, 4> plane_poi
     bool keep_lower = attributes.has(ModelObjectCutAttribute::KeepLower);
     bool cut_to_parts = attributes.has(ModelObjectCutAttribute::CutToParts);
     ModelObject* upper = keep_upper ? ModelObject::new_clone(*this) : nullptr;
-    ModelObject* lower = cut_to_parts ? upper : (keep_lower ? ModelObject::new_clone(*this) : nullptr);
+    ModelObject* lower = (cut_to_parts&&upper!=nullptr) ? upper : (keep_lower ? ModelObject::new_clone(*this) : nullptr);
 
     if (attributes.has(ModelObjectCutAttribute::KeepUpper)) {
         upper->set_model(nullptr);
@@ -1598,7 +1612,7 @@ ModelObjectPtrs ModelObject::cut(size_t instance, std::array<Vec3d, 4> plane_poi
             if (attributes.has(ModelObjectCutAttribute::KeepLower))
                 lower->add_volume(*volume);
         }
-        else if (! volume->mesh().empty()) {            
+        else if (! volume->mesh().empty()) {
             // Transform the mesh by the combined transformation matrix.
             // Flip the triangles in case the composite transformation is left handed.
 			TriangleMesh mesh(volume->mesh());
@@ -1744,7 +1758,7 @@ ModelObjectPtrs ModelObject::segment(size_t instance, unsigned int max_extruders
             // Modifiers are not cut, but we still need to add the instance transformation
             // to the modifier volume transformation to preserve their shape properly.
             volume->set_transformation(Geometry::Transformation(instance_matrix * volume_matrix));
-            upper->add_volume(*volume); 
+            upper->add_volume(*volume);
         }
         else if (!volume->mesh().empty()) {
             // Transform the mesh by the combined transformation matrix.
@@ -2185,7 +2199,7 @@ void ModelObject::print_info() const
     using namespace std;
     cout << fixed;
     boost::nowide::cout << "[" << boost::filesystem::path(this->input_file).filename().string() << "]" << endl;
-    
+
     TriangleMesh mesh = this->raw_mesh();
     BoundingBoxf3 bb = mesh.bounding_box();
     Vec3d size = bb.size();
@@ -2203,7 +2217,7 @@ void ModelObject::print_info() const
     cout << "manifold = "   << (mesh.stats().manifold() ? "yes" : "no") << endl;
     if (! mesh.stats().manifold())
         cout << "open_edges = " << mesh.stats().open_edges << endl;
-    
+
     if (mesh.stats().repaired()) {
         const RepairedMeshErrors& stats = mesh.stats().repaired_errors;
         if (stats.degenerate_facets > 0)
@@ -2286,7 +2300,7 @@ void ModelVolume::set_material_id(t_model_material_id material_id)
 }
 
 ModelMaterial* ModelVolume::material() const
-{ 
+{
     return this->object->get_model()->get_material(m_material_id);
 }
 
@@ -2362,8 +2376,10 @@ void ModelVolume::center_geometry_after_creation(bool update_source_offset)
     Vec3d shift = this->mesh().bounding_box().center();
     if (!shift.isApprox(Vec3d::Zero()))
     {
-    	if (m_mesh)
-        	const_cast<TriangleMesh*>(m_mesh.get())->translate(-(float)shift(0), -(float)shift(1), -(float)shift(2));
+        if (m_mesh) {
+            const_cast<TriangleMesh*>(m_mesh.get())->translate(-(float)shift(0), -(float)shift(1), -(float)shift(2));
+            const_cast<TriangleMesh*>(m_mesh.get())->set_init_shift(shift);
+        }
         if (m_convex_hull)
 			const_cast<TriangleMesh*>(m_convex_hull.get())->translate(-(float)shift(0), -(float)shift(1), -(float)shift(2));
         translate(shift);
@@ -2803,7 +2819,7 @@ double Model::getThermalLength(const std::vector<ModelVolume*> modelVolumePtrs)
     }
     return thermalLength;
 }
-// max printing speed, difference in bed temperature and envirument temperature and bed adhension coefficients are considered 
+// max printing speed, difference in bed temperature and envirument temperature and bed adhension coefficients are considered
 double ModelInstance::get_auto_brim_width(double deltaT, double adhension) const
 {
     BoundingBoxf3 raw_bbox = object->raw_mesh_bounding_box();
@@ -2896,7 +2912,7 @@ double ModelInstance::get_auto_brim_width() const
 void ModelInstance::get_arrange_polygon(void* ap) const
 {
 //    static const double SIMPLIFY_TOLERANCE_MM = 0.1;
-    
+
     Vec3d rotation = get_rotation();
     rotation.z()   = 0.;
     Transform3d trafo_instance =
@@ -2909,7 +2925,7 @@ void ModelInstance::get_arrange_polygon(void* ap) const
 //        pp = p.simplify(scaled<double>(SIMPLIFY_TOLERANCE_MM));
 //        if (!pp.empty()) p = pp.front();
 //    }
-   
+
     arrangement::ArrangePolygon& ret = *(arrangement::ArrangePolygon*)ap;
     ret.poly.contour = std::move(p);
     ret.translation  = Vec2crd{scaled(get_offset(X)), scaled(get_offset(Y))};
@@ -3039,6 +3055,12 @@ void FacetsAnnotation::set_triangle_from_string(int triangle_id, const std::stri
     }
 }
 
+bool FacetsAnnotation::equals(const FacetsAnnotation &other) const
+{
+    const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>>& data = other.get_data();
+    return (m_data == data);
+}
+
 // Test whether the two models contain the same number of ModelObjects with the same set of IDs
 // ordered in the same order. In that case it is not necessary to kill the background processing.
 bool model_object_list_equal(const Model &model_old, const Model &model_new)
@@ -3140,22 +3162,22 @@ bool model_property_changed(const ModelObject &model_object_old, const ModelObje
 
 bool model_custom_supports_data_changed(const ModelObject& mo, const ModelObject& mo_new)
 {
-    return model_property_changed(mo, mo_new, 
-        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; }, 
+    return model_property_changed(mo, mo_new,
+        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; },
         [](const ModelVolume &mv_old, const ModelVolume &mv_new){ return mv_old.supported_facets.timestamp_matches(mv_new.supported_facets); });
 }
 
 bool model_custom_seam_data_changed(const ModelObject& mo, const ModelObject& mo_new)
 {
-    return model_property_changed(mo, mo_new, 
-        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; }, 
+    return model_property_changed(mo, mo_new,
+        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; },
         [](const ModelVolume &mv_old, const ModelVolume &mv_new){ return mv_old.seam_facets.timestamp_matches(mv_new.seam_facets); });
 }
 
 bool model_mmu_segmentation_data_changed(const ModelObject& mo, const ModelObject& mo_new)
 {
-    return model_property_changed(mo, mo_new, 
-        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; }, 
+    return model_property_changed(mo, mo_new,
+        [](const ModelVolumeType t) { return t == ModelVolumeType::MODEL_PART; },
         [](const ModelVolume &mv_old, const ModelVolume &mv_new){ return mv_old.mmu_segmentation_facets.timestamp_matches(mv_new.mmu_segmentation_facets); });
 }
 
@@ -3189,7 +3211,7 @@ bool model_has_advanced_features(const Model &model)
 void check_model_ids_validity(const Model &model)
 {
     std::set<ObjectID> ids;
-    auto check = [&ids](ObjectID id) { 
+    auto check = [&ids](ObjectID id) {
         assert(id.valid());
         assert(ids.find(id) == ids.end());
         ids.insert(id);

src/libslic3r/Model.hpp

@@ -18,6 +18,8 @@
 #include "Format/bbs_3mf.hpp"
 //BBS: add step
 #include "Format/STEP.hpp"
+//BBS: add stl
+#include "Format/STL.hpp"
 
 #include <map>
 #include <memory>
@@ -306,6 +308,7 @@ public:
     ModelVolume*            add_volume(TriangleMesh &&mesh, ModelVolumeType type = ModelVolumeType::MODEL_PART);
     ModelVolume*            add_volume(const ModelVolume &volume, ModelVolumeType type = ModelVolumeType::INVALID);
     ModelVolume*            add_volume(const ModelVolume &volume, TriangleMesh &&mesh);
+    ModelVolume*            add_volume_with_shared_mesh(const ModelVolume &other, ModelVolumeType type = ModelVolumeType::MODEL_PART);
     void                    delete_volume(size_t idx);
     void                    clear_volumes();
     void                    sort_volumes(bool full_sort);
@@ -615,6 +618,7 @@ public:
     void set_triangle_from_string(int triangle_id, const std::string& str);
     // After deserializing the last triangle, shrink data to fit.
     void shrink_to_fit() { m_data.first.shrink_to_fit(); m_data.second.shrink_to_fit(); }
+    bool equals(const FacetsAnnotation &other) const;
 
 private:
     // Constructors to be only called by derived classes.
@@ -674,6 +678,7 @@ public:
 
     // The triangular model.
     const TriangleMesh& mesh() const { return *m_mesh.get(); }
+    const TriangleMesh* mesh_ptr() const { return m_mesh.get(); }
     void                set_mesh(const TriangleMesh &mesh) { m_mesh = std::make_shared<const TriangleMesh>(mesh); }
     void                set_mesh(TriangleMesh &&mesh) { m_mesh = std::make_shared<const TriangleMesh>(std::move(mesh)); }
     void                set_mesh(const indexed_triangle_set &mesh) { m_mesh = std::make_shared<const TriangleMesh>(mesh); }
@@ -857,6 +862,18 @@ private:
         if (mesh.facets_count() > 1)
             calculate_convex_hull();
     }
+    ModelVolume(ModelObject *object, const std::shared_ptr<const TriangleMesh> &mesh, ModelVolumeType type = ModelVolumeType::MODEL_PART) : m_mesh(mesh), m_type(type), object(object)
+    {
+		assert(this->id().valid());
+        assert(this->config.id().valid());
+        assert(this->supported_facets.id().valid());
+        assert(this->seam_facets.id().valid());
+        assert(this->mmu_segmentation_facets.id().valid());
+        assert(this->id() != this->config.id());
+        assert(this->id() != this->supported_facets.id());
+        assert(this->id() != this->seam_facets.id());
+        assert(this->id() != this->mmu_segmentation_facets.id());
+    }
     ModelVolume(ModelObject *object, TriangleMesh &&mesh, TriangleMesh &&convex_hull, ModelVolumeType type = ModelVolumeType::MODEL_PART) :
 		m_mesh(new TriangleMesh(std::move(mesh))), m_convex_hull(new TriangleMesh(std::move(convex_hull))), m_type(type), object(object) {
 		assert(this->id().valid());
@@ -1282,7 +1299,7 @@ public:
         DynamicPrintConfig* config = nullptr, ConfigSubstitutionContext* config_substitutions = nullptr,
         LoadStrategy options = LoadStrategy::AddDefaultInstances, PlateDataPtrs* plate_data = nullptr,
         std::vector<Preset*>* project_presets = nullptr, bool* is_xxx = nullptr, Semver* file_version = nullptr, Import3mfProgressFn proFn = nullptr,
-        ImportStepProgressFn stepFn = nullptr, StepIsUtf8Fn stepIsUtf8Fn = nullptr, BBLProject* project = nullptr);
+        ImportstlProgressFn stlFn = nullptr, ImportStepProgressFn stepFn = nullptr, StepIsUtf8Fn stepIsUtf8Fn = nullptr, BBLProject* project = nullptr);
     // BBS
     static double findMaxSpeed(const ModelObject* object);
     static double getThermalLength(const ModelVolume* modelVolumePtr);

src/libslic3r/MultiPoint.cpp

@@ -63,6 +63,23 @@ int MultiPoint::find_point(const Point &point) const
     return -1;  // not found
 }
 
+int MultiPoint::find_point(const Point &point, double scaled_epsilon) const
+{
+    if (scaled_epsilon == 0) return this->find_point(point);
+
+    auto dist2_min = std::numeric_limits<double>::max();
+    auto eps2      = scaled_epsilon * scaled_epsilon;
+    int  idx_min   = -1;
+    for (const Point &pt : this->points) {
+        double d2 = (pt - point).cast<double>().squaredNorm();
+        if (d2 < dist2_min) {
+            idx_min   = int(&pt - &this->points.front());
+            dist2_min = d2;
+        }
+    }
+    return dist2_min < eps2 ? idx_min : -1;
+}
+
 bool MultiPoint::has_boundary_point(const Point &point) const
 {
     double dist = (point.projection_onto(*this) - point).cast<double>().norm();

src/libslic3r/MultiPoint.hpp

@@ -43,7 +43,12 @@ public:
     double length() const;
     bool   is_valid() const { return this->points.size() >= 2; }
 
+    // Return index of a polygon point exactly equal to point.
+    // Return -1 if no such point exists.
     int  find_point(const Point &point) const;
+    // Return index of the closest point to point closer than scaled_epsilon.
+    // Return -1 if no such point exists.
+    int  find_point(const Point &point, const double scaled_epsilon) const;
     bool has_boundary_point(const Point &point) const;
     int  closest_point_index(const Point &point) const {
         int idx = -1;

src/libslic3r/NormalUtils.cpp

@@ -0,0 +1,142 @@
+#include "NormalUtils.hpp"
+
+using namespace Slic3r;
+
+Vec3f NormalUtils::create_triangle_normal(
+    const stl_triangle_vertex_indices &indices,
+    const std::vector<stl_vertex> &    vertices)
+{
+    const stl_vertex &v0        = vertices[indices[0]];
+    const stl_vertex &v1        = vertices[indices[1]];
+    const stl_vertex &v2        = vertices[indices[2]];
+    Vec3f             direction = (v1 - v0).cross(v2 - v0);
+    direction.normalize();
+    return direction;
+}
+
+std::vector<Vec3f> NormalUtils::create_triangle_normals(
+    const indexed_triangle_set &its)
+{
+    std::vector<Vec3f> normals;
+    normals.reserve(its.indices.size());
+    for (const Vec3crd &index : its.indices) {
+        normals.push_back(create_triangle_normal(index, its.vertices));
+    }
+    return normals;
+}
+
+NormalUtils::Normals NormalUtils::create_normals_average_neighbor(
+    const indexed_triangle_set &its)
+{
+    size_t             count_vertices = its.vertices.size();
+    std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
+    std::vector<unsigned int> count(count_vertices, 0);
+    for (const Vec3crd &indice : its.indices) {
+        Vec3f normal = create_triangle_normal(indice, its.vertices);
+        for (int i = 0; i < 3; ++i) {
+            normals[indice[i]] += normal;
+            ++count[indice[i]];
+        }
+    }
+    // normalize to size 1
+    for (auto &normal : normals) {
+        size_t index = &normal - &normals.front();
+        normal /= static_cast<float>(count[index]);
+    }
+    return normals;
+}
+
+// calc triangle angle of vertex defined by index to triangle indices
+float NormalUtils::indice_angle(int                            i,
+                                const Vec3crd &                indice,
+                                const std::vector<stl_vertex> &vertices)
+{
+    int i1 = (i == 0) ? 2 : (i - 1);
+    int i2 = (i == 2) ? 0 : (i + 1);
+
+    Vec3f v1 = vertices[i1] - vertices[i];
+    Vec3f v2 = vertices[i2] - vertices[i];
+
+    v1.normalize();
+    v2.normalize();
+
+    float w = v1.dot(v2);
+    if (w > 1.f)
+        w = 1.f;
+    else if (w < -1.f)
+        w = -1.f;
+    return acos(w);
+}
+
+NormalUtils::Normals NormalUtils::create_normals_angle_weighted(
+    const indexed_triangle_set &its)
+{
+    size_t             count_vertices = its.vertices.size();
+    std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
+    std::vector<float> count(count_vertices, 0.f);
+    for (const Vec3crd &indice : its.indices) {
+        Vec3f normal = create_triangle_normal(indice, its.vertices);
+        Vec3f angles(indice_angle(0, indice, its.vertices),
+                     indice_angle(1, indice, its.vertices), 0.f);
+        angles[2] = (M_PI - angles[0] - angles[1]);
+        for (int i = 0; i < 3; ++i) {
+            const float &weight = angles[i];
+            normals[indice[i]] += normal * weight;
+            count[indice[i]] += weight;
+        }
+    }
+    // normalize to size 1
+    for (auto &normal : normals) {
+        size_t index = &normal - &normals.front();
+        normal /= count[index];
+    }
+    return normals;
+}
+
+NormalUtils::Normals NormalUtils::create_normals_nelson_weighted(
+    const indexed_triangle_set &its)
+{
+    size_t             count_vertices = its.vertices.size();
+    std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
+    std::vector<float> count(count_vertices, 0.f);
+    const std::vector<stl_vertex> &vertices = its.vertices;
+    for (const Vec3crd &indice : its.indices) {
+        Vec3f normal = create_triangle_normal(indice, vertices);
+
+        const stl_vertex &v0 = vertices[indice[0]];
+        const stl_vertex &v1 = vertices[indice[1]];
+        const stl_vertex &v2 = vertices[indice[2]];
+
+        float e0 = (v0 - v1).norm();
+        float e1 = (v1 - v2).norm();
+        float e2 = (v2 - v0).norm();
+
+        Vec3f coefs(e0 * e2, e0 * e1, e1 * e2);
+        for (int i = 0; i < 3; ++i) {
+            const float &weight = coefs[i];
+            normals[indice[i]] += normal * weight;
+            count[indice[i]] += weight;
+        }
+    }
+    // normalize to size 1
+    for (auto &normal : normals) {
+        size_t index = &normal - &normals.front();
+        normal /= count[index];
+    }
+    return normals;
+}
+
+// calculate normals by averaging normals of neghbor triangles
+std::vector<Vec3f> NormalUtils::create_normals(
+    const indexed_triangle_set &its, VertexNormalType type)
+{
+    switch (type) {
+    case VertexNormalType::AverageNeighbor:
+        return create_normals_average_neighbor(its);
+    case VertexNormalType::AngleWeighted:
+        return create_normals_angle_weighted(its);
+    case VertexNormalType::NelsonMaxWeighted:
+    default:
+        return create_normals_nelson_weighted(its);
+    }
+}

src/libslic3r/NormalUtils.hpp

@@ -0,0 +1,69 @@
+#ifndef slic3r_NormalUtils_hpp_
+#define slic3r_NormalUtils_hpp_
+
+#include "Point.hpp"
+#include "Model.hpp"
+
+namespace Slic3r {
+
+/// <summary>
+/// Collection of static function
+/// to create normals
+/// </summary>
+class NormalUtils
+{
+public:
+    using Normal = Vec3f;
+    using Normals = std::vector<Normal>;
+    NormalUtils() = delete; // only static functions
+
+    enum class VertexNormalType {
+        AverageNeighbor,
+        AngleWeighted,
+        NelsonMaxWeighted
+    };
+
+    /// <summary>
+    /// Create normal for triangle defined by indices from vertices
+    /// </summary>
+    /// <param name="indices">index into vertices</param>
+    /// <param name="vertices">vector of vertices</param>
+    /// <returns>normal to triangle(normalized to size 1)</returns>
+    static Normal create_triangle_normal(
+        const stl_triangle_vertex_indices &indices,
+        const std::vector<stl_vertex> &    vertices);
+
+    /// <summary>
+    /// Create normals for each vertices
+    /// </summary>
+    /// <param name="its">indices and vertices</param>
+    /// <returns>Vector of normals</returns>
+    static Normals create_triangle_normals(const indexed_triangle_set &its);
+
+    /// <summary>
+    /// Create normals for each vertex by averaging neighbor triangles normal
+    /// </summary>
+    /// <param name="its">Triangle indices and vertices</param>
+    /// <param name="type">Type of calculation normals</param>
+    /// <returns>Normal for each vertex</returns>
+    static Normals create_normals(
+        const indexed_triangle_set &its,
+        VertexNormalType type = VertexNormalType::NelsonMaxWeighted);
+    static Normals create_normals_average_neighbor(const indexed_triangle_set &its);
+    static Normals create_normals_angle_weighted(const indexed_triangle_set &its);
+    static Normals create_normals_nelson_weighted(const indexed_triangle_set &its);
+
+    /// <summary>
+    /// Calculate angle of trinagle side.
+    /// </summary>
+    /// <param name="i">index to indices, define angle point</param>
+    /// <param name="indice">address to vertices</param>
+    /// <param name="vertices">vertices data</param>
+    /// <returns>Angle [in radian]</returns>
+    static float indice_angle(int                            i,
+                              const Vec3crd &                indice,
+                              const std::vector<stl_vertex> &vertices);
+};
+
+} // namespace Slic3r
+#endif // slic3r_NormalUtils_hpp_

src/libslic3r/Polyline.cpp

@@ -517,6 +517,28 @@ bool remove_degenerate(Polylines &polylines)
     return modified;
 }
 
+std::pair<int, Point> foot_pt(const Points &polyline, const Point &pt)
+{
+    if (polyline.size() < 2) return std::make_pair(-1, Point(0, 0));
+
+    auto  d2_min = std::numeric_limits<double>::max();
+    Point foot_pt_min;
+    Point prev    = polyline.front();
+    auto  it      = polyline.begin();
+    auto  it_proj = polyline.begin();
+    for (++it; it != polyline.end(); ++it) {
+        Point  foot_pt = pt.projection_onto(Line(prev, *it));
+        double d2      = (foot_pt - pt).cast<double>().squaredNorm();
+        if (d2 < d2_min) {
+            d2_min      = d2;
+            foot_pt_min = foot_pt;
+            it_proj     = it;
+        }
+        prev = *it;
+    }
+    return std::make_pair(int(it_proj - polyline.begin()) - 1, foot_pt_min);
+}
+
 ThickLines ThickPolyline::thicklines() const
 {
     ThickLines lines;

src/libslic3r/Polyline.hpp

@@ -222,6 +222,9 @@ const Point& leftmost_point(const Polylines &polylines);
 
 bool remove_degenerate(Polylines &polylines);
 
+// Returns index of a segment of a polyline and foot point of pt on polyline.
+std::pair<int, Point> foot_pt(const Points &polyline, const Point &pt);
+
 class ThickPolyline : public Polyline {
 public:
     ThickPolyline() : endpoints(std::make_pair(false, false)) {}

src/libslic3r/Preset.cpp

@@ -691,7 +691,7 @@ static std::vector<std::string> s_Preset_filament_options {
     "filament_flow_ratio", "filament_density", "filament_cost", "filament_minimal_purge_on_wipe_tower",
     "chamber_temperature", "nozzle_temperature", "nozzle_temperature_initial_layer",
     // BBS
-    "cool_plate_temp", "eng_plate_temp", "hot_plate_temp", "cool_plate_temp_initial_layer", "eng_plate_temp_initial_layer", "hot_plate_temp_initial_layer",
+    "cool_plate_temp", "eng_plate_temp", "hot_plate_temp", "textured_plate_temp", "cool_plate_temp_initial_layer", "eng_plate_temp_initial_layer", "hot_plate_temp_initial_layer","textured_plate_temp_initial_layer",
     // "bed_type",
     //BBS:temperature_vitrification
     "temperature_vitrification", "reduce_fan_stop_start_freq", "slow_down_for_layer_cooling", "fan_min_speed",
@@ -1400,7 +1400,7 @@ bool PresetCollection::load_user_preset(std::string name, std::map<std::string,
 
     //filament_id
     std::string cloud_filament_id;
-    if ((m_type == Preset::TYPE_FILAMENT) && preset_values.find(BBL_JSON_KEY_FILAMENT_ID) == preset_values.end()) {
+    if ((m_type == Preset::TYPE_FILAMENT) && preset_values.find(BBL_JSON_KEY_FILAMENT_ID) != preset_values.end()) {
         cloud_filament_id = preset_values[BBL_JSON_KEY_FILAMENT_ID];
     }
 
@@ -1830,9 +1830,12 @@ std::pair<Preset*, bool> PresetCollection::load_external_preset(
 
 Preset& PresetCollection::load_preset(const std::string &path, const std::string &name, DynamicPrintConfig &&config, bool select)
 {
+    lock();
     auto it = this->find_preset_internal(name);
     if (it == m_presets.end() || it->name != name) {
         // The preset was not found. Create a new preset.
+        if (m_presets.begin() + m_idx_selected >= it)
+            ++m_idx_selected;
         it = m_presets.emplace(it, Preset(m_type, name, false));
     }
     Preset &preset = *it;
@@ -1842,6 +1845,7 @@ Preset& PresetCollection::load_preset(const std::string &path, const std::string
     preset.is_dirty = false;
     if (select)
         this->select_preset_by_name(name, true);
+    unlock();
     //BBS: add config related logs
     BOOST_LOG_TRIVIAL(debug) << __FUNCTION__ << boost::format(", preset type %1%, name %2%, path %3%, is_system %4%, is_default %5% is_visible %6%")%Preset::get_type_string(m_type) %preset.name %preset.file %preset.is_system %preset.is_default %preset.is_visible;
     return preset;

src/libslic3r/PresetBundle.cpp

@@ -592,6 +592,102 @@ PresetsConfigSubstitutions PresetBundle::load_user_presets(AppConfig &config, st
     return substitutions;
 }
 
+PresetsConfigSubstitutions PresetBundle::import_presets(std::vector<std::string> &              files,
+                                                        std::function<int(std::string const &)> override_confirm,
+                                                        ForwardCompatibilitySubstitutionRule    rule)
+{
+    PresetsConfigSubstitutions substitutions;
+    int overwrite = 0;
+    std::vector<std::string>   result;
+    for (auto &file : files) {
+        if (Slic3r::is_json_file(file)) {
+            try {
+                DynamicPrintConfig config;
+                // BBS: change to json format
+                // ConfigSubstitutions config_substitutions = config.load_from_ini(preset.file, substitution_rule);
+                std::map<std::string, std::string> key_values;
+                std::string                        reason;
+                ConfigSubstitutions                config_substitutions = config.load_from_json(file, rule, key_values, reason);
+                std::string name = key_values[BBL_JSON_KEY_NAME];
+                std::string version_str = key_values[BBL_JSON_KEY_VERSION];
+                boost::optional<Semver> version = Semver::parse(version_str);
+                if (!version) continue;
+                Semver app_version = *(Semver::parse(SLIC3R_VERSION));
+                if (version->maj() != app_version.maj()) {
+                    BOOST_LOG_TRIVIAL(warning) << "Preset incompatibla, not loading: " << name;
+                    continue;
+                }
+
+                PresetCollection * collection = nullptr;
+                if (config.has("printer_settings_id"))
+                    collection = &printers;
+                else if (config.has("print_settings_id"))
+                    collection = &prints;
+                else if (config.has("filament_settings_id"))
+                    collection = &filaments;
+                if (collection == nullptr) {
+                    BOOST_LOG_TRIVIAL(warning) << "Preset type is unknown, not loading: " << name;
+                    continue;
+                }
+                if (auto p = collection->find_preset(name, false)) {
+                    if (p->is_default || p->is_system) {
+                        BOOST_LOG_TRIVIAL(warning) << "Preset already present and is system preset, not loading: " << name;
+                        continue;
+                    }
+                    overwrite = override_confirm(name);
+                    if (overwrite == 0 || overwrite == 2) {
+                        BOOST_LOG_TRIVIAL(warning) << "Preset already present, not loading: " << name;
+                        continue;
+                    }
+                }
+
+                DynamicPrintConfig new_config;
+                Preset *      inherit_preset  = nullptr;
+                ConfigOption *inherits_config = config.option(BBL_JSON_KEY_INHERITS);
+                if (inherits_config) {
+                    ConfigOptionString *option_str     = dynamic_cast<ConfigOptionString *>(inherits_config);
+                    std::string         inherits_value = option_str->value;
+                    inherit_preset = collection->find_preset(inherits_value, false, true);
+                }
+                if (inherit_preset) {
+                    new_config = inherit_preset->config;
+                } else {
+                    // Find a default preset for the config. The PrintPresetCollection provides different default preset based on the "printer_technology" field.
+                    // new_config = default_preset.config;
+                    // we should skip this preset here
+                    BOOST_LOG_TRIVIAL(warning) << __FUNCTION__ << boost::format(", can not find inherit preset for user preset %1%, just skip") % name;
+                    continue;
+                }
+                new_config.apply(std::move(config));
+                
+                Preset &preset     = collection->load_preset(collection->path_from_name(name), name, std::move(new_config), false);
+                preset.is_external = true;
+                preset.version     = *version;
+                if (inherit_preset)
+                    preset.base_id     = inherit_preset->setting_id;
+                Preset::normalize(preset.config);
+                // Report configuration fields, which are misplaced into a wrong group.
+                const Preset &default_preset = collection->default_preset_for(new_config);
+                std::string incorrect_keys = Preset::remove_invalid_keys(preset.config, default_preset.config);
+                if (!incorrect_keys.empty())
+                    BOOST_LOG_TRIVIAL(error) << "Error in a preset file: The preset \"" << preset.file << "\" contains the following incorrect keys: " << incorrect_keys
+                                                << ", which were removed";
+                if (!config_substitutions.empty())
+                    substitutions.push_back({name, collection->type(), PresetConfigSubstitutions::Source::UserFile, file, std::move(config_substitutions)});
+
+                preset.save(inherit_preset ? &inherit_preset->config : nullptr);
+                result.push_back(file);
+            } catch (const std::ifstream::failure &err) {
+                BOOST_LOG_TRIVIAL(error) << boost::format("The config cannot be loaded: %1%. Reason: %2%") % file % err.what();
+            } catch (const std::runtime_error &err) {
+                BOOST_LOG_TRIVIAL(error) << boost::format("Failed importing config file: %1%. Reason: %2%") % file % err.what();
+            }
+        }
+    }
+    files = result;
+    return substitutions;
+}
+
 //BBS save user preset to user_id preset folder
 void PresetBundle::save_user_presets(AppConfig& config, std::vector<std::string>& need_to_delete_list)
 {
@@ -626,6 +722,17 @@ void PresetBundle::update_user_presets_directory(const std::string preset_folder
     BOOST_LOG_TRIVIAL(debug) << __FUNCTION__ << boost::format(" finished");
 }
 
+void PresetBundle::remove_user_presets_directory(const std::string preset_folder)
+{
+    const std::string dir_user_presets = data_dir() + "/" + PRESET_USER_DIR + "/" + preset_folder;
+
+    BOOST_LOG_TRIVIAL(debug) << __FUNCTION__ << boost::format(" enter, delete directory : %1%") % dir_user_presets;
+    fs::path folder(dir_user_presets);
+    if (fs::exists(folder)) {
+        fs::remove_all(folder);
+    }
+}
+
 void PresetBundle::update_system_preset_setting_ids(std::map<std::string, std::map<std::string, std::string>>& system_presets)
 {
     for (auto iterator: system_presets)
@@ -3215,6 +3322,36 @@ void PresetBundle::update_compatible(PresetSelectCompatibleType select_other_pri
     }
 }*/
 
+std::vector<std::string> PresetBundle::export_current_configs(const std::string &                     path,
+                                                              std::function<int(std::string const &)> override_confirm,
+                                                              bool                                    include_modify,
+                                                              bool                                    export_system_settings)
+{
+    const Preset &print_preset    = include_modify ? prints.get_edited_preset() : prints.get_selected_preset();
+    const Preset &printer_preset  = include_modify ? printers.get_edited_preset() : printers.get_selected_preset();
+    std::set<Preset const *> presets { &print_preset, &printer_preset };
+    for (auto &f : filament_presets) {
+        auto filament_preset = filaments.find_preset(f, include_modify);
+        if (filament_preset) presets.insert(filament_preset);
+    }
+
+    int overwrite = 0;
+    std::vector<std::string> result;
+    for (auto preset : presets) {
+        if ((preset->is_system  && !export_system_settings) || preset->is_default)
+            continue;
+        std::string file = path + "/" + preset->name + ".json";
+        if (boost::filesystem::exists(file) && overwrite < 2) {
+            overwrite = override_confirm(preset->name);
+            if (overwrite == 0 || overwrite == 2)
+                continue;
+        }
+        preset->config.save_to_json(file, preset->name, "", preset->version.to_string());
+        result.push_back(file);
+    }
+    return result;
+}
+
 // Set the filament preset name. As the name could come from the UI selection box,
 // an optional "(modified)" suffix will be removed from the filament name.
 void PresetBundle::set_filament_preset(size_t idx, const std::string &name)

src/libslic3r/PresetBundle.hpp

@@ -46,9 +46,11 @@ public:
 
     // BBS Load user presets
     PresetsConfigSubstitutions load_user_presets(AppConfig &config, std::map<std::string, std::map<std::string, std::string>>& my_presets, ForwardCompatibilitySubstitutionRule rule);
+    PresetsConfigSubstitutions import_presets(std::vector<std::string> &files, std::function<int(std::string const &)> override_confirm, ForwardCompatibilitySubstitutionRule rule);
     void save_user_presets(AppConfig& config, std::vector<std::string>& need_to_delete_list);
     void remove_users_preset(AppConfig &config);
-    void  update_user_presets_directory(const std::string preset_folder);
+    void update_user_presets_directory(const std::string preset_folder);
+    void remove_user_presets_directory(const std::string preset_folder);
     void update_system_preset_setting_ids(std::map<std::string, std::map<std::string, std::string>>& system_presets);
 
     //BBS: add API to get previous machine
@@ -164,6 +166,8 @@ public:
     //void export_current_configbundle(const std::string &path);
     //BBS: add a function to export system presets for cloud-slicer
     //void export_system_configs(const std::string &path);
+    std::vector<std::string> export_current_configs(const std::string &path, std::function<int(std::string const &)> override_confirm, 
+        bool include_modify, bool export_system_settings = false);
 
     // Enable / disable the "- default -" preset.
     void                        set_default_suppressed(bool default_suppressed);

src/libslic3r/Print.cpp

@@ -96,6 +96,7 @@ bool Print::invalidate_state_by_config_options(const ConfigOptionResolver & /* n
         "cool_plate_temp_initial_layer",
         "eng_plate_temp_initial_layer",
         "hot_plate_temp_initial_layer",
+        "textured_plate_temp_initial_layer",
         "gcode_add_line_number",
         "layer_change_gcode",
         "fan_min_speed",
@@ -175,6 +176,7 @@ bool Print::invalidate_state_by_config_options(const ConfigOptionResolver & /* n
             || opt_key == "cool_plate_temp"
             || opt_key == "eng_plate_temp"
             || opt_key == "hot_plate_temp"
+            || opt_key == "textured_plate_temp"
             || opt_key == "enable_prime_tower"
             || opt_key == "prime_tower_width"
             || opt_key == "prime_tower_brim_width"
@@ -1132,6 +1134,46 @@ void Print::auto_assign_extruders(ModelObject* model_object) const
     }
 }
 
+void  PrintObject::set_shared_object(PrintObject *object)
+{
+    m_shared_object = object;
+    BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": this=%1%, found shared object from %2%")%this%m_shared_object;
+}
+
+void  PrintObject::clear_shared_object()
+{
+    if (m_shared_object) {
+        BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": this=%1%, clear previous shared object data %2%")%this %m_shared_object;
+        m_layers.clear();
+        m_support_layers.clear();
+        m_tree_support_layers.clear();
+
+        m_shared_object = nullptr;
+
+        invalidate_all_steps_without_cancel();
+    }
+}
+
+void  PrintObject::copy_layers_from_shared_object()
+{
+    if (m_shared_object) {
+        m_layers.clear();
+        m_support_layers.clear();
+        m_tree_support_layers.clear();
+
+        firstLayerObjSliceByVolume.clear();
+        firstLayerObjSliceByGroups.clear();
+
+        BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": this=%1%, copied layers from object %2%")%this%m_shared_object;
+        m_layers = m_shared_object->layers();
+        m_support_layers = m_shared_object->support_layers();
+        m_tree_support_layers = m_shared_object->tree_support_layers();
+
+        firstLayerObjSliceByVolume = m_shared_object->firstLayerObjSlice();
+        firstLayerObjSliceByGroups = m_shared_object->firstLayerObjGroups();
+    }
+}
+
 // BBS
 BoundingBox PrintObject::get_first_layer_bbox(float& a, float& layer_height, std::string& name)
 {
@@ -1181,15 +1223,115 @@ void Print::process()
 {
     name_tbb_thread_pool_threads_set_locale();
 
+    //compute the PrintObject with the same geometries
+    BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": this=%1%, enter")%this;
+    for (PrintObject *obj : m_objects)
+        obj->clear_shared_object();
+
+    //add the print_object share check logic
+    auto is_print_object_the_same = [this](const PrintObject* object1, const PrintObject* object2) -> bool{
+        if (object1->trafo().matrix() != object2->trafo().matrix())
+            return false;
+        const ModelObject* model_obj1 = object1->model_object();
+        const ModelObject* model_obj2 = object2->model_object();
+        if (model_obj1->volumes.size() != model_obj2->volumes.size())
+            return false;
+        bool has_extruder1 = model_obj1->config.has("extruder");
+        bool has_extruder2 = model_obj2->config.has("extruder");
+        if ((has_extruder1 != has_extruder2)
+            || (has_extruder1 && model_obj1->config.extruder() != model_obj2->config.extruder()))
+            return false;
+        for (int index = 0; index < model_obj1->volumes.size(); index++) {
+            const ModelVolume &model_volume1 = *model_obj1->volumes[index];
+            const ModelVolume &model_volume2 = *model_obj2->volumes[index];
+            if (model_volume1.type() != model_volume2.type())
+                return false;
+            if (model_volume1.mesh_ptr() != model_volume2.mesh_ptr())
+                return false;
+            has_extruder1 = model_volume1.config.has("extruder");
+            has_extruder2 = model_volume2.config.has("extruder");
+            if ((has_extruder1 != has_extruder2)
+                || (has_extruder1 && model_volume1.config.extruder() != model_volume2.config.extruder()))
+                return false;
+            if (!model_volume1.supported_facets.equals(model_volume2.supported_facets))
+                return false;
+            if (!model_volume1.seam_facets.equals(model_volume2.seam_facets))
+                return false;
+            if (!model_volume1.mmu_segmentation_facets.equals(model_volume2.mmu_segmentation_facets))
+                return false;
+            if (model_volume1.config.get() != model_volume2.config.get())
+                return false;
+        }
+        //if (!object1->config().equals(object2->config()))
+        //    return false;
+        if (model_obj1->config.get() != model_obj2->config.get())
+            return false;
+        return true;
+    };
+    int object_count = m_objects.size();
+    std::set<PrintObject*> need_slicing_objects;
+    for (int index = 0; index < object_count; index++)
+    {
+        PrintObject *obj =  m_objects[index];
+        for (PrintObject *slicing_obj : need_slicing_objects)
+        {
+            if (is_print_object_the_same(obj, slicing_obj)) {
+                obj->set_shared_object(slicing_obj);
+                break;
+            }
+        }
+        if (!obj->get_shared_object())
+            need_slicing_objects.insert(obj);
+    }
+
+    BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": total object counts %1% in current print, need to slice %2%")%m_objects.size()%need_slicing_objects.size();
     BOOST_LOG_TRIVIAL(info) << "Starting the slicing process." << log_memory_info();
+    for (PrintObject *obj : m_objects) {
+        if (need_slicing_objects.count(obj) != 0) {
+            obj->make_perimeters();
+        }
+        else {
+            if (obj->set_started(posSlice))
+                obj->set_done(posSlice);
+            if (obj->set_started(posPerimeters))
+                obj->set_done(posPerimeters);
+        }
+    }
+    for (PrintObject *obj : m_objects) {
+        if (need_slicing_objects.count(obj) != 0) {
+            obj->infill();
+        }
+        else {
+            if (obj->set_started(posPrepareInfill))
+                obj->set_done(posPrepareInfill);
+            if (obj->set_started(posInfill))
+                obj->set_done(posInfill);
+        }
+    }
+    for (PrintObject *obj : m_objects) {
+        if (need_slicing_objects.count(obj) != 0) {
+            obj->ironing();
+        }
+        else {
+            if (obj->set_started(posIroning))
+                obj->set_done(posIroning);
+        }
+    }
+    for (PrintObject *obj : m_objects) {
+        if (need_slicing_objects.count(obj) != 0) {
+            obj->generate_support_material();
+        }
+        else {
+            if (obj->set_started(posSupportMaterial))
+                obj->set_done(posSupportMaterial);
+        }
+    }
+
     for (PrintObject *obj : m_objects)
-        obj->make_perimeters();
-    for (PrintObject *obj : m_objects)
-        obj->infill();
-    for (PrintObject *obj : m_objects)
-        obj->ironing();
-    for (PrintObject *obj : m_objects)
-        obj->generate_support_material();
+    {
+        if (need_slicing_objects.count(obj) == 0)
+            obj->copy_layers_from_shared_object();
+    }
     if (this->set_started(psWipeTower)) {
         m_wipe_tower_data.clear();
         m_tool_ordering.clear();
@@ -1292,8 +1434,17 @@ void Print::process()
         this->set_done(psSkirtBrim);
     }
     //BBS
-    for (PrintObject *obj : m_objects)
-        obj->simplify_extrusion_path();
+    for (PrintObject *obj : m_objects) {
+        if (need_slicing_objects.count(obj) != 0) {
+            obj->simplify_extrusion_path();
+        }
+        else {
+            if (obj->set_started(posSimplifyPath))
+                obj->set_done(posSimplifyPath);
+            if (obj->set_started(posSimplifySupportPath))
+                obj->set_done(posSimplifySupportPath);
+        }
+    }
 
     BOOST_LOG_TRIVIAL(info) << "Slicing process finished." << log_memory_info();
 }
@@ -1634,7 +1785,7 @@ void Print::_make_wipe_tower()
         for (LayerTools& layer_tools : layer_tools_array) {
             layer_tools.has_wipe_tower = true;
             if (layer_tools.wipe_tower_partitions == 0) {
-                layer_tools.wipe_tower_partitions = 1; 
+                layer_tools.wipe_tower_partitions = 1;
             }
         }
     }

src/libslic3r/Print.hpp

@@ -49,6 +49,12 @@ struct groupedVolumeSlices
     ExPolygons              slices;
 };
 
+enum SupportNecessaryType {
+    NoNeedSupp=0,
+    SharpTail,
+    LargeOverhang,
+};
+
 namespace FillAdaptive {
     struct Octree;
     struct OctreeDeleter;
@@ -409,6 +415,11 @@ public:
     //BBS
     BoundingBox get_first_layer_bbox(float& area, float& layer_height, std::string& name);
 
+    PrintObject* get_shared_object() const { return m_shared_object; }
+    void         set_shared_object(PrintObject *object);
+    void         clear_shared_object();
+    void         copy_layers_from_shared_object();
+
     // BBS: Boundingbox of the first layer
     BoundingBox                 firstLayerObjectBrimBoundingBox;
 private:
@@ -457,7 +468,7 @@ private:
     std::pair<FillAdaptive::OctreePtr, FillAdaptive::OctreePtr> prepare_adaptive_infill_data();
 
     // BBS
-    bool is_support_necessary();
+    SupportNecessaryType is_support_necessary();
 
     // XYZ in scaled coordinates
     Vec3crd									m_size;
@@ -489,6 +500,8 @@ private:
     // BBS: per object skirt
     ExtrusionEntityCollection               m_skirt;
 
+    PrintObject*                            m_shared_object{ nullptr };
+
  public:
     //BBS: When printing multi-material objects, this settings will make slicer to clip the overlapping object parts one by the other.
     //(2nd part will be clipped by the 1st, 3rd part will be clipped by the 1st and 2nd etc).

src/libslic3r/PrintBase.hpp

@@ -629,6 +629,8 @@ protected:
         { return m_state.invalidate_multiple(il.begin(), il.end(), PrintObjectBase::cancel_callback(m_print)); }
     bool            invalidate_all_steps()
         { return m_state.invalidate_all(PrintObjectBase::cancel_callback(m_print)); }
+    bool            invalidate_all_steps_without_cancel()
+        { return m_state.invalidate_all([](){}); }
 
     bool            is_step_started_unguarded(PrintObjectStepEnum step) const { return m_state.is_started_unguarded(step); }
     bool            is_step_done_unguarded(PrintObjectStepEnum step) const { return m_state.is_done_unguarded(step); }

src/libslic3r/PrintConfig.cpp

@@ -230,7 +230,8 @@ static const t_config_enum_values s_keys_map_OverhangFanThreshold = {
 static const t_config_enum_values s_keys_map_BedType = {
     { "Cool Plate",         btPC },
     { "Engineering Plate",  btEP  },
-    { "High Temp Plate",    btPEI  }
+    { "High Temp Plate",    btPEI  },
+    { "Textured PEI Plate",      btPTE }
 };
 CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(BedType)
 
@@ -462,6 +463,16 @@ void PrintConfigDef::init_fff_params()
     def->max = 120;
     def->set_default_value(new ConfigOptionInts{ 45 });
 
+    def             = this->add("textured_plate_temp", coInts);
+    def->label      = L("Other layers");
+    def->tooltip    = L("Bed temperature for layers except the initial one. "
+                     "Value 0 means the filament does not support to print on the Textured PEI Plate");
+    def->sidetext   = L("°C");
+    def->full_label = L("Bed temperature");
+    def->min        = 0;
+    def->max        = 120;
+    def->set_default_value(new ConfigOptionInts{45});
+
     def = this->add("cool_plate_temp_initial_layer", coInts);
     def->label = L("Initial layer");
     def->full_label = L("Initial layer bed temperature");
@@ -492,6 +503,15 @@ void PrintConfigDef::init_fff_params()
     def->max = 120;
     def->set_default_value(new ConfigOptionInts{ 45 });
 
+    def             = this->add("textured_plate_temp_initial_layer", coInts);
+    def->label      = L("Initial layer");
+    def->full_label = L("Initial layer bed temperature");
+    def->tooltip    = L("Bed temperature of the initial layer. "
+                     "Value 0 means the filament does not support to print on the Textured PEI Plate");
+    def->sidetext   = L("°C");
+    def->max        = 0;
+    def->max        = 120;
+    def->set_default_value(new ConfigOptionInts{45});
 
     def = this->add("curr_bed_type", coEnums);
     def->label = L("Bed type");
@@ -501,9 +521,11 @@ void PrintConfigDef::init_fff_params()
     def->enum_values.emplace_back("Cool Plate");
     def->enum_values.emplace_back("Engineering Plate");
     def->enum_values.emplace_back("High Temp Plate");
+    def->enum_values.emplace_back("Textured PEI Plate");
     def->enum_labels.emplace_back(L("Cool Plate"));
     def->enum_labels.emplace_back(L("Engineering Plate"));
     def->enum_labels.emplace_back(L("High Temp Plate"));
+    def->enum_labels.emplace_back(L("Textured PEI Plate"));
     def->set_default_value(new ConfigOptionEnum<BedType>(btPC));
 
     def = this->add("before_layer_change_gcode", coString);
@@ -1091,7 +1113,7 @@ void PrintConfigDef::init_fff_params()
     def->tooltip = L("Filament price. For statistics only");
     def->sidetext = L("money/kg");
     def->min = 0;
-    def->mode = comDevelop;
+    def->mode = comAdvanced;
     def->set_default_value(new ConfigOptionFloats { 0. });
 
     def = this->add("filament_settings_id", coStrings);
@@ -2000,9 +2022,11 @@ void PrintConfigDef::init_fff_params()
     def->enum_values.push_back("nearest");
     def->enum_values.push_back("aligned");
     def->enum_values.push_back("back");
+    def->enum_values.push_back("random");
     def->enum_labels.push_back(L("Nearest"));
     def->enum_labels.push_back(L("Aligned"));
     def->enum_labels.push_back(L("Back"));
+    def->enum_labels.push_back(L("Random"));
     def->mode = comSimple;
     def->set_default_value(new ConfigOptionEnum<SeamPosition>(spAligned));
 
@@ -2108,10 +2132,12 @@ void PrintConfigDef::init_fff_params()
 
     def = this->add("timelapse_no_toolhead", coBool);
     def->label = L("Timelapse");
-    def->tooltip = L("Record timelapse video of printing without showing toolhead. In this mode "
-                    "the toolhead docks near the excess chute at each layer change, and then "
-                    "a snapshot is taken with the chamber camera. When printing finishes a timelapse "
-                    "video is composed of all the snapshots.");
+    def->tooltip = L("If enabled, a timelapse video will be generated for each print. "
+                     "After each layer is printed, the toolhead will move to the excess chute, "
+                     "and then a snapshot is taken with the chamber camera. "
+                     "All of these snapshots are composed into a timelapse video when printing completes. "
+                     "Since the melt filament may leak from the nozzle during the process of taking a snapshot, "
+                     "prime tower is required for nozzle priming.");
     def->mode = comSimple;
     def->set_default_value(new ConfigOptionBool(false));
 
@@ -2406,8 +2432,7 @@ void PrintConfigDef::init_fff_params()
     def = this->add("independent_support_layer_height", coBool);
     def->label = L("Independent support layer height");
     def->category = L("Support");
-    def->tooltip = L("Support layer uses layer height independent with object layer. This is to support custom support gap,"
-                   "but may cause extra filament switches if support is specified as different extruder with object");
+    def->tooltip = L("Support layer uses layer height independent with object layer. This is to support customizing z-gap and save print time.");
     def->mode = comAdvanced;
     def->set_default_value(new ConfigOptionBool(true));
 

src/libslic3r/PrintConfig.hpp

@@ -160,6 +160,7 @@ enum BedType {
     btPC = 0,
     btEP,
     btPEI,
+    btPTE,
     btCount
 };
 
@@ -185,6 +186,8 @@ static std::string bed_type_to_gcode_string(const BedType type)
     case btPEI:
         type_str = "high_temp_plate";
         break;
+    case btPTE:
+        type_str = "frosted_plate";
     default:
         type_str = "unknown";
         break;
@@ -204,6 +207,9 @@ static std::string get_bed_temp_key(const BedType type)
     if (type == btPEI)
         return "hot_plate_temp";
 
+    if (type == btPTE)
+        return "textured_plate_temp";
+
     return "";
 }
 
@@ -218,6 +224,9 @@ static std::string get_bed_temp_1st_layer_key(const BedType type)
     if (type == btPEI)
         return "hot_plate_temp_initial_layer";
 
+    if (type == btPTE)
+        return "textured_plate_temp_initial_layer";
+
     return "";
 }
 
@@ -786,9 +795,11 @@ PRINT_CONFIG_CLASS_DERIVED_DEFINE(
     ((ConfigOptionInts,               cool_plate_temp))
     ((ConfigOptionInts,               eng_plate_temp))
     ((ConfigOptionInts,               hot_plate_temp)) // hot is short for high temperature
+    ((ConfigOptionInts,               textured_plate_temp))
     ((ConfigOptionInts,               cool_plate_temp_initial_layer))
     ((ConfigOptionInts,               eng_plate_temp_initial_layer))
     ((ConfigOptionInts,               hot_plate_temp_initial_layer)) // hot is short for high temperature
+    ((ConfigOptionInts,               textured_plate_temp_initial_layer))
     ((ConfigOptionBools,              enable_overhang_bridge_fan))
     ((ConfigOptionInts,               overhang_fan_speed))
     ((ConfigOptionEnumsGeneric,       overhang_fan_threshold))

src/libslic3r/PrintObject.cpp

@@ -86,6 +86,7 @@ PrintObject::PrintObject(Print* print, ModelObject* model_object, const Transfor
 
 PrintObject::~PrintObject()
 {
+    BOOST_LOG_TRIVIAL(info) << __FUNCTION__ << boost::format(": this=%1%, m_shared_object %2%")%this%m_shared_object;
     if (m_shared_regions && -- m_shared_regions->m_ref_cnt == 0) delete m_shared_regions;
     clear_layers();
     clear_support_layers();
@@ -419,11 +420,17 @@ void PrintObject::generate_support_material()
             m_print->throw_if_canceled();
         } else {
             // BBS: pop a warning if objects have significant amount of overhangs but support material is not enabled
-            if (this->is_support_necessary()) {
+            SupportNecessaryType sntype = this->is_support_necessary();
+            if (sntype != NoNeedSupp) {
                 m_print->set_status(50, L("Checking support necessity"));
-
-                std::string warning_message = format(L("It seems object %s needs support to print. Please enable support generation."), this->model_object()->name);
-                this->active_step_add_warning(PrintStateBase::WarningLevel::CRITICAL, warning_message, PrintStateBase::SlicingNeedSupportOn);
+                if (sntype == SharpTail) {
+                    std::string warning_message = format(L("It seems object %s has completely floating regions. Please re-orient the object or enable support generation."),
+                                                         this->model_object()->name);
+                    this->active_step_add_warning(PrintStateBase::WarningLevel::CRITICAL, warning_message, PrintStateBase::SlicingNeedSupportOn);
+                } else {
+                    std::string warning_message = format(L("It seems object %s has large overhangs. Please enable support generation."), this->model_object()->name);
+                    this->active_step_add_warning(PrintStateBase::WarningLevel::CRITICAL, warning_message, PrintStateBase::SlicingNeedSupportOn);
+                }
             }
 
 #if 0
@@ -529,9 +536,11 @@ std::pair<FillAdaptive::OctreePtr, FillAdaptive::OctreePtr> PrintObject::prepare
 
 void PrintObject::clear_layers()
 {
-    for (Layer *l : m_layers)
-        delete l;
-    m_layers.clear();
+    if (!m_shared_object) {
+        for (Layer *l : m_layers)
+            delete l;
+        m_layers.clear();
+    }
 }
 
 Layer* PrintObject::add_layer(int id, coordf_t height, coordf_t print_z, coordf_t slice_z)
@@ -567,9 +576,11 @@ SupportLayer* PrintObject::get_support_layer_at_printz(coordf_t print_z, coordf_
 
 void PrintObject::clear_tree_support_layers()
 {
-    for (TreeSupportLayer* l : m_tree_support_layers)
-        delete l;
-    m_tree_support_layers.clear();
+    if (!m_shared_object) {
+        for (TreeSupportLayer* l : m_tree_support_layers)
+            delete l;
+        m_tree_support_layers.clear();
+    }
 }
 
 std::shared_ptr<TreeSupportData> PrintObject::alloc_tree_support_preview_cache()
@@ -596,9 +607,11 @@ TreeSupportLayer* PrintObject::add_tree_support_layer(int id, coordf_t height, c
 
 void PrintObject::clear_support_layers()
 {
-    for (Layer *l : m_support_layers)
-        delete l;
-    m_support_layers.clear();
+    if (!m_shared_object) {
+        for (Layer *l : m_support_layers)
+            delete l;
+        m_support_layers.clear();
+    }
 }
 
 SupportLayer* PrintObject::add_support_layer(int id, int interface_id, coordf_t height, coordf_t print_z)
@@ -2448,7 +2461,7 @@ template void PrintObject::remove_bridges_from_contacts<Polygons>(
     float max_bridge_length, bool break_bridge);
 
 
-bool PrintObject::is_support_necessary()
+SupportNecessaryType PrintObject::is_support_necessary()
 {
     static const double super_overhang_area_threshold = SQ(scale_(5.0));
 
@@ -2471,7 +2484,7 @@ bool PrintObject::is_support_necessary()
             for (const ExPolygon& expoly : layerm->raw_slices) {
                 // detect sharp tail
                 if (intersection_ex({ expoly }, lower_layer_offseted).empty())
-                    return true;
+                    return SharpTail;
             }
         }
 
@@ -2503,18 +2516,19 @@ bool PrintObject::is_support_necessary()
         double super_overhang_area = 0.0;
         for (Polygon& poly : super_overhang_polys) {
             bool is_ccw = poly.is_counter_clockwise();
+            double area_  = poly.area();
             if (is_ccw) {
-                if (super_overhang_area > super_overhang_area_threshold)
-                    return true;
-                super_overhang_area = poly.area();
+                if (area_ > super_overhang_area_threshold)
+                    return LargeOverhang;
+                super_overhang_area += area_;
             }
             else {
-                super_overhang_area -= poly.area();
+                super_overhang_area -= area_;
             }
         }
 
-        if (super_overhang_area > super_overhang_area_threshold)
-            return true;
+        //if (super_overhang_area > super_overhang_area_threshold)
+        //    return LargeOverhang;
 
         // 3. check overhang distance
         const double distance_threshold_scaled = extrusion_width_scaled * 2;
@@ -2529,10 +2543,10 @@ bool PrintObject::is_support_necessary()
             }),
             exceed_overhang.end());
         if (!exceed_overhang.empty())
-            return true;
+            return LargeOverhang;
     }
 
-    return false;
+    return NoNeedSupp;
 }
 
 static void project_triangles_to_slabs(ConstLayerPtrsAdaptor layers, const indexed_triangle_set &custom_facets, const Transform3f &tr, bool seam, std::vector<Polygons> &out)

src/libslic3r/PrintObjectSlice.cpp

@@ -474,9 +474,7 @@ std::string fix_slicing_errors(PrintObject* object, LayerPtrs &layers, const std
 
         if (layers[idx_layer]->slicing_errors) {
             buggy_layers.push_back(idx_layer);
-            //BBS
-            error_msg = L("Empty layers around bottom are replaced by nearest normal layers.");
-            }
+        }
         else
             break; // only detect empty layers near bed
     }
@@ -484,14 +482,15 @@ std::string fix_slicing_errors(PrintObject* object, LayerPtrs &layers, const std
     BOOST_LOG_TRIVIAL(debug) << "Slicing objects - fixing slicing errors in parallel - begin";
     tbb::parallel_for(
         tbb::blocked_range<size_t>(0, buggy_layers.size()),
-        [&layers, &throw_if_canceled, &buggy_layers](const tbb::blocked_range<size_t>& range) {
+        [&layers, &throw_if_canceled, &buggy_layers, &error_msg](const tbb::blocked_range<size_t>& range) {
             for (size_t buggy_layer_idx = range.begin(); buggy_layer_idx < range.end(); ++ buggy_layer_idx) {
                 throw_if_canceled();
                 size_t idx_layer = buggy_layers[buggy_layer_idx];
-                // BBS: only replace empty first layer
-                if (idx_layer > 0)
+                // BBS: only replace empty layers lower than 1mm
+                const coordf_t thresh_empty_layer_height = 1;
+                Layer* layer = layers[idx_layer];
+                if (layer->print_z>= thresh_empty_layer_height)
                     continue;
-                Layer *layer     = layers[idx_layer];
                 assert(layer->slicing_errors);
                 // Try to repair the layer surfaces by merging all contours and all holes from neighbor layers.
                 // BOOST_LOG_TRIVIAL(trace) << "Attempting to repair layer" << idx_layer;
@@ -500,42 +499,39 @@ std::string fix_slicing_errors(PrintObject* object, LayerPtrs &layers, const std
                     // Find the first valid layer below / above the current layer.
                     const Surfaces *upper_surfaces = nullptr;
                     const Surfaces *lower_surfaces = nullptr;
-                    //BBS: only repair first layer if the 2nd layer is Good
-                    for (size_t j = idx_layer + 1; j < /*layers.size()*/2; ++ j)
-                        if (! layers[j]->slicing_errors) {
+                    //BBS: only repair empty layers lowers than 1mm
+                    for (size_t j = idx_layer + 1; j < layers.size(); ++j) {
+                        if (!layers[j]->slicing_errors) {
                             upper_surfaces = &layers[j]->regions()[region_id]->slices.surfaces;
                             break;
                         }
-                    for (int j = /*int(idx_layer) -*/ 1; j >= 0; -- j)
-                        if (! layers[j]->slicing_errors) {
+                        if (layers[j]->print_z >= thresh_empty_layer_height) break;
+                    }
+                    for (int j = int(idx_layer) - 1; j >= 0; --j) {
+                        if (layers[j]->print_z >= thresh_empty_layer_height) continue;
+                        if (!layers[j]->slicing_errors) {
                             lower_surfaces = &layers[j]->regions()[region_id]->slices.surfaces;
                             break;
                         }
+                    }
                     // Collect outer contours and holes from the valid layers above & below.
-                    Polygons outer;
-                    outer.reserve(
+                    ExPolygons expolys;
+                    expolys.reserve(
                         ((upper_surfaces == nullptr) ? 0 : upper_surfaces->size()) +
                         ((lower_surfaces == nullptr) ? 0 : lower_surfaces->size()));
-                    size_t num_holes = 0;
                     if (upper_surfaces)
                         for (const auto &surface : *upper_surfaces) {
-                            outer.push_back(surface.expolygon.contour);
-                            num_holes += surface.expolygon.holes.size();
+                            expolys.emplace_back(surface.expolygon);
                         }
                     if (lower_surfaces)
                         for (const auto &surface : *lower_surfaces) {
-                            outer.push_back(surface.expolygon.contour);
-                            num_holes += surface.expolygon.holes.size();
+                            expolys.emplace_back(surface.expolygon);
                         }
-                    Polygons holes;
-                    holes.reserve(num_holes);
-                    if (upper_surfaces)
-                        for (const auto &surface : *upper_surfaces)
-                            polygons_append(holes, surface.expolygon.holes);
-                    if (lower_surfaces)
-                        for (const auto &surface : *lower_surfaces)
-                            polygons_append(holes, surface.expolygon.holes);
-                    layerm->slices.set(diff_ex(union_(outer), holes), stInternal);
+                    if (!expolys.empty()) {
+                        //BBS
+                        error_msg = L("Empty layers around bottom are replaced by nearest normal layers.");
+                        layerm->slices.set(union_ex(expolys), stInternal);
+                    }
                 }
                 // Update layer slices after repairing the single regions.
                 layer->make_slices();
@@ -555,8 +551,7 @@ std::string fix_slicing_errors(PrintObject* object, LayerPtrs &layers, const std
 
     //BBS
     if(error_msg.empty() && !buggy_layers.empty())
-        error_msg = L("The model has overlapping or self-intersecting facets. I tried to repair it, "
-            "however you might want to check the results or repair the input file and retry.");
+        error_msg = L("The model has too many empty layers.");
     return error_msg;
 }
 

src/libslic3r/ProjectTask.hpp

@@ -51,6 +51,7 @@ public:
     {
         profile_ = profile;
         prediction = 0;
+        weight = 0.0f;
     }
 
     BBLSliceInfo(const BBLSliceInfo& obj) {

src/libslic3r/Shape/TextShape.cpp

@@ -0,0 +1,203 @@
+#include "../libslic3r.h"
+#include "../Model.hpp"
+#include "../TriangleMesh.hpp"
+
+#include "TextShape.hpp"
+
+#include <string>
+#include <vector>
+
+#include "Standard_TypeDef.hxx"
+#include "STEPCAFControl_Reader.hxx"
+#include "BRepMesh_IncrementalMesh.hxx"
+#include "Interface_Static.hxx"
+#include "XCAFDoc_DocumentTool.hxx"
+#include "XCAFDoc_ShapeTool.hxx"
+#include "XCAFApp_Application.hxx"
+#include "TopoDS_Solid.hxx"
+#include "TopoDS_Compound.hxx"
+#include "TopoDS_Builder.hxx"
+#include "TopoDS.hxx"
+#include "TDataStd_Name.hxx"
+#include "BRepBuilderAPI_Transform.hxx"
+#include "TopExp_Explorer.hxx"
+#include "TopExp_Explorer.hxx"
+#include "BRep_Tool.hxx"
+#include "Font_BRepFont.hxx"
+#include "Font_BRepTextBuilder.hxx"
+#include "BRepPrimAPI_MakePrism.hxx"
+#include "Font_FontMgr.hxx"
+
+namespace Slic3r {
+
+std::vector<std::string> init_occt_fonts()
+{
+    std::vector<std::string> stdFontNames;
+
+    Handle(Font_FontMgr) aFontMgr = Font_FontMgr::GetInstance();
+    aFontMgr->InitFontDataBase();
+
+    TColStd_SequenceOfHAsciiString availFontNames;
+    aFontMgr->GetAvailableFontsNames(availFontNames);
+    stdFontNames.reserve(availFontNames.Size());
+
+    for (auto afn : availFontNames)
+        stdFontNames.push_back(afn->ToCString());
+
+    return stdFontNames;
+}
+
+static bool TextToBRep(const char* text, const char* font, const float theTextHeight, Font_FontAspect& theFontAspect, TopoDS_Shape& theShape)
+{
+    Standard_Integer anArgIt = 1;
+    Standard_CString aName = "text_shape";
+    Standard_CString aText = text;
+
+    Font_BRepFont           aFont;
+    //TCollection_AsciiString aFontName("Courier");
+    TCollection_AsciiString aFontName(font);
+    Standard_Real           aTextHeight = theTextHeight;
+    Font_FontAspect         aFontAspect = theFontAspect;
+    Standard_Boolean        anIsCompositeCurve = Standard_False;
+    gp_Ax3                  aPenAx3(gp::XOY());
+    gp_Dir                  aNormal(0.0, 0.0, 1.0);
+    gp_Dir                  aDirection(1.0, 0.0, 0.0);
+    gp_Pnt                  aPenLoc;
+
+    Graphic3d_HorizontalTextAlignment aHJustification = Graphic3d_HTA_LEFT;
+    Graphic3d_VerticalTextAlignment   aVJustification = Graphic3d_VTA_BOTTOM;
+    Font_StrictLevel aStrictLevel = Font_StrictLevel_Any;
+
+    aFont.SetCompositeCurveMode(anIsCompositeCurve);
+    if (!aFont.FindAndInit(aFontName.ToCString(), aFontAspect, aTextHeight, aStrictLevel))
+        return false;
+
+    aPenAx3 = gp_Ax3(aPenLoc, aNormal, aDirection);
+
+    Font_BRepTextBuilder aBuilder;
+    theShape = aBuilder.Perform(aFont, aText, aPenAx3, aHJustification, aVJustification);
+    return true;
+}
+
+static bool Prism(const TopoDS_Shape& theBase, const float thickness, TopoDS_Shape& theSolid)
+{
+    if (theBase.IsNull()) return false;
+
+    gp_Vec V(0.f, 0.f, thickness);
+    BRepPrimAPI_MakePrism* Prism = new BRepPrimAPI_MakePrism(theBase, V, Standard_False);
+
+    theSolid = Prism->Shape();
+    return true;
+}
+
+static void MakeMesh(TopoDS_Shape& theSolid, TriangleMesh& theMesh)
+{
+    const double STEP_TRANS_CHORD_ERROR = 0.005;
+    const double STEP_TRANS_ANGLE_RES = 1;
+
+    BRepMesh_IncrementalMesh mesh(theSolid, STEP_TRANS_CHORD_ERROR, false, STEP_TRANS_ANGLE_RES, true);
+    int aNbNodes = 0;
+    int aNbTriangles = 0;
+    for (TopExp_Explorer anExpSF(theSolid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
+        TopLoc_Location aLoc;
+        Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(anExpSF.Current()), aLoc);
+        if (!aTriangulation.IsNull()) {
+            aNbNodes += aTriangulation->NbNodes();
+            aNbTriangles += aTriangulation->NbTriangles();
+        }
+    }
+
+    stl_file stl;
+    stl.stats.type = inmemory;
+    stl.stats.number_of_facets = (uint32_t)aNbTriangles;
+    stl.stats.original_num_facets = stl.stats.number_of_facets;
+    stl_allocate(&stl);
+
+    std::vector<Vec3f> points;
+    points.reserve(aNbNodes);
+    //BBS: count faces missing triangulation
+    Standard_Integer aNbFacesNoTri = 0;
+    //BBS: fill temporary triangulation
+    Standard_Integer aNodeOffset = 0;
+    Standard_Integer aTriangleOffet = 0;
+    for (TopExp_Explorer anExpSF(theSolid, TopAbs_FACE); anExpSF.More(); anExpSF.Next()) {
+        const TopoDS_Shape& aFace = anExpSF.Current();
+        TopLoc_Location aLoc;
+        Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(TopoDS::Face(aFace), aLoc);
+        if (aTriangulation.IsNull()) {
+            ++aNbFacesNoTri;
+            continue;
+        }
+        //BBS: copy nodes
+        gp_Trsf aTrsf = aLoc.Transformation();
+        for (Standard_Integer aNodeIter = 1; aNodeIter <= aTriangulation->NbNodes(); ++aNodeIter) {
+            gp_Pnt aPnt = aTriangulation->Node(aNodeIter);
+            aPnt.Transform(aTrsf);
+            points.emplace_back(std::move(Vec3f(aPnt.X(), aPnt.Y(), aPnt.Z())));
+        }
+        //BBS: copy triangles
+        const TopAbs_Orientation anOrientation = anExpSF.Current().Orientation();
+        for (Standard_Integer aTriIter = 1; aTriIter <= aTriangulation->NbTriangles(); ++aTriIter) {
+            Poly_Triangle aTri = aTriangulation->Triangle(aTriIter);
+
+            Standard_Integer anId[3];
+            aTri.Get(anId[0], anId[1], anId[2]);
+            if (anOrientation == TopAbs_REVERSED) {
+                //BBS: swap 1, 2.
+                Standard_Integer aTmpIdx = anId[1];
+                anId[1] = anId[2];
+                anId[2] = aTmpIdx;
+            }
+            //BBS: Update nodes according to the offset.
+            anId[0] += aNodeOffset;
+            anId[1] += aNodeOffset;
+            anId[2] += aNodeOffset;
+            //BBS: save triangles facets
+            stl_facet facet;
+            facet.vertex[0] = points[anId[0] - 1].cast<float>();
+            facet.vertex[1] = points[anId[1] - 1].cast<float>();
+            facet.vertex[2] = points[anId[2] - 1].cast<float>();
+            facet.extra[0] = 0;
+            facet.extra[1] = 0;
+            stl_normal normal;
+            stl_calculate_normal(normal, &facet);
+            stl_normalize_vector(normal);
+            facet.normal = normal;
+            stl.facet_start[aTriangleOffet + aTriIter - 1] = facet;
+        }
+
+        aNodeOffset += aTriangulation->NbNodes();
+        aTriangleOffet += aTriangulation->NbTriangles();
+    }
+
+    theMesh.from_stl(stl);
+}
+
+void load_text_shape(const char*text, const char* font, const float text_height, const float thickness, bool is_bold, bool is_italic, TriangleMesh& text_mesh)
+{
+    Handle(Font_FontMgr) aFontMgr = Font_FontMgr::GetInstance();
+    if (aFontMgr->GetAvailableFonts().IsEmpty())
+        aFontMgr->InitFontDataBase();
+
+    TopoDS_Shape aTextBase;
+    Font_FontAspect aFontAspect = Font_FontAspect_UNDEFINED;
+    if (is_bold && is_italic)
+        aFontAspect = Font_FontAspect_BoldItalic;
+    else if (is_bold)
+        aFontAspect = Font_FontAspect_Bold;
+    else if (is_italic)
+        aFontAspect = Font_FontAspect_Italic;
+    else
+        aFontAspect = Font_FontAspect_Regular;
+
+    if (!TextToBRep(text, font, text_height, aFontAspect, aTextBase))
+        return;
+
+    TopoDS_Shape aTextShape;
+    if (!Prism(aTextBase, thickness, aTextShape))
+        return;
+
+    MakeMesh(aTextShape, text_mesh);
+}
+
+}; // namespace Slic3r

src/libslic3r/Shape/TextShape.hpp

@@ -0,0 +1,12 @@
+#ifndef slic3r_Text_Shape_hpp_
+#define slic3r_Text_Shape_hpp_
+
+namespace Slic3r {
+class TriangleMesh;
+
+extern std::vector<std::string> init_occt_fonts();
+extern void load_text_shape(const char* text, const char* font, const float text_height, const float thickness, bool is_bold, bool is_italic, TriangleMesh& text_mesh);
+
+}; // namespace Slic3r
+
+#endif // slic3r_Text_Shape_hpp_

src/libslic3r/ShortEdgeCollapse.cpp

@@ -0,0 +1,183 @@
+#include "ShortEdgeCollapse.hpp"
+#include "libslic3r/NormalUtils.hpp"
+
+#include <unordered_map>
+#include <unordered_set>
+#include <random>
+#include <algorithm>
+
+namespace Slic3r {
+
+void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) {
+    // whenever vertex is removed, its mapping is update to the index of vertex with wich it merged
+    std::vector<size_t> vertices_index_mapping(mesh.vertices.size());
+    for (size_t idx = 0; idx < vertices_index_mapping.size(); ++idx) {
+        vertices_index_mapping[idx] = idx;
+    }
+    // Algorithm uses get_final_index query to get the actual vertex index. The query also updates all mappings on the way, essentially flattening the mapping
+    std::vector<size_t> flatten_queue;
+    auto get_final_index = [&vertices_index_mapping, &flatten_queue](const size_t &orig_index) {
+        flatten_queue.clear();
+        size_t idx = orig_index;
+        while (vertices_index_mapping[idx] != idx) {
+            flatten_queue.push_back(idx);
+            idx = vertices_index_mapping[idx];
+        }
+        for (size_t i : flatten_queue) {
+            vertices_index_mapping[i] = idx;
+        }
+        return idx;
+
+    };
+
+    // if face is removed, mark it here
+    std::vector<bool> face_removal_flags(mesh.indices.size(), false);
+
+    std::vector<Vec3i> triangles_neighbors = its_face_neighbors_par(mesh);
+
+    // now compute vertices dot product - this is used during edge collapse,
+    // to determine which vertex to remove and which to keep;  We try to keep the one with larger angle, because it defines the shape "more".
+    // The min vertex dot product is lowest dot product of its normal with the normals of faces around it.
+    // the lower the dot product, the more we want to keep the vertex
+    // NOTE: This score is not updated, even though the decimation does change the mesh. It saves computation time, and there are no strong reasons to update.
+    std::vector<float> min_vertex_dot_product(mesh.vertices.size(), 1);
+    {
+        std::vector<Vec3f> face_normals = its_face_normals(mesh);
+        std::vector<Vec3f> vertex_normals = NormalUtils::create_normals(mesh);
+
+        for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) {
+            Vec3i t = mesh.indices[face_idx];
+            Vec3f n = face_normals[face_idx];
+            min_vertex_dot_product[t[0]] = std::min(min_vertex_dot_product[t[0]], n.dot(vertex_normals[t[0]]));
+            min_vertex_dot_product[t[1]] = std::min(min_vertex_dot_product[t[1]], n.dot(vertex_normals[t[1]]));
+            min_vertex_dot_product[t[2]] = std::min(min_vertex_dot_product[t[2]], n.dot(vertex_normals[t[2]]));
+        }
+    }
+
+    // lambda to remove face. It flags the face as removed, and updates neighbourhood info
+    auto remove_face = [&triangles_neighbors, &face_removal_flags](int face_idx, int other_face_idx) {
+        if (face_idx < 0) {
+            return;
+        }
+        face_removal_flags[face_idx] = true;
+        Vec3i neighbors = triangles_neighbors[face_idx];
+        int n_a = neighbors[0] != other_face_idx ? neighbors[0] : neighbors[1];
+        int n_b = neighbors[2] != other_face_idx ? neighbors[2] : neighbors[1];
+        if (n_a > 0)
+            for (int &n : triangles_neighbors[n_a]) {
+                if (n == face_idx) {
+                    n = n_b;
+                    break;
+                }
+            }
+        if (n_b > 0)
+            for (int &n : triangles_neighbors[n_b]) {
+                if (n == face_idx) {
+                    n = n_a;
+                    break;
+                }
+            }
+    };
+
+    std::mt19937_64 generator { 27644437 };// default constant seed! so that results are deterministic
+    std::vector<size_t> face_indices(mesh.indices.size());
+    for (size_t idx = 0; idx < face_indices.size(); ++idx) {
+        face_indices[idx] = idx;
+    }
+    //tmp face indices used only for swapping
+    std::vector<size_t> tmp_face_indices(mesh.indices.size());
+
+    float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all
+    float edge_len = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased
+
+    while (face_indices.size() > target_triangle_count) {
+        // simpple func to increase the edge len - if decimation ratio is low, it increases the len up to twice, if decimation ratio is high, increments are low
+        edge_len = edge_len * (1.0f + 1.0 - decimation_ratio);
+        float max_edge_len_squared = edge_len * edge_len;
+
+        //shuffle the faces and traverse in random order, this MASSIVELY improves the quality of the result
+        std::shuffle(face_indices.begin(), face_indices.end(), generator);
+
+        for (const size_t &face_idx : face_indices) {
+            if (face_removal_flags[face_idx]) {
+                // if face already removed from previous collapses, skip (each collapse removes two triangles [at least] )
+                continue;
+            }
+
+            // look at each edge if it is good candidate for collapse
+            for (size_t edge_idx = 0; edge_idx < 3; ++edge_idx) {
+                size_t vertex_index_keep = get_final_index(mesh.indices[face_idx][edge_idx]);
+                size_t vertex_index_remove = get_final_index(mesh.indices[face_idx][(edge_idx + 1) % 3]);
+                //check distance, skip long edges
+                if ((mesh.vertices[vertex_index_keep] - mesh.vertices[vertex_index_remove]).squaredNorm()
+                        > max_edge_len_squared) {
+                    continue;
+                }
+                // swap indexes if vertex_index_keep has higher dot product (we want to keep low dot product vertices)
+                if (min_vertex_dot_product[vertex_index_remove] < min_vertex_dot_product[vertex_index_keep]) {
+                    size_t tmp = vertex_index_keep;
+                    vertex_index_keep = vertex_index_remove;
+                    vertex_index_remove = tmp;
+                }
+
+                //remove vertex
+                {
+                    // map its index to the index of the kept vertex
+                    vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep];
+                }
+
+                int neighbor_to_remove_face_idx = triangles_neighbors[face_idx][edge_idx];
+                // remove faces
+                remove_face(face_idx, neighbor_to_remove_face_idx);
+                remove_face(neighbor_to_remove_face_idx, face_idx);
+
+                // break. this triangle is done
+                break;
+            }
+        }
+
+        // filter face_indices, remove those that have been collapsed
+        size_t prev_size = face_indices.size();
+        tmp_face_indices.clear();
+        for (size_t face_idx : face_indices) {
+            if (!face_removal_flags[face_idx]){
+                tmp_face_indices.push_back(face_idx);
+            }
+        }
+        face_indices.swap(tmp_face_indices);
+
+        decimation_ratio = float(prev_size - face_indices.size()) / float(prev_size);
+        //std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl;
+    }
+
+    //Extract the result mesh
+    std::unordered_map<size_t, size_t> final_vertices_mapping;
+    std::vector<Vec3f> final_vertices;
+    std::vector<Vec3i> final_indices;
+    final_indices.reserve(face_indices.size());
+    for (size_t idx : face_indices) {
+        Vec3i final_face;
+        for (size_t i = 0; i < 3; ++i) {
+            final_face[i] = get_final_index(mesh.indices[idx][i]);
+        }
+        if (final_face[0] == final_face[1] || final_face[1] == final_face[2] || final_face[2] == final_face[0]) {
+            continue; // discard degenerate triangles
+        }
+
+        for (size_t i = 0; i < 3; ++i) {
+            if (final_vertices_mapping.find(final_face[i]) == final_vertices_mapping.end()) {
+                final_vertices_mapping[final_face[i]] = final_vertices.size();
+                final_vertices.push_back(mesh.vertices[final_face[i]]);
+            }
+            final_face[i] = final_vertices_mapping[final_face[i]];
+        }
+
+        final_indices.push_back(final_face);
+    }
+
+    mesh.vertices = final_vertices;
+    mesh.indices = final_indices;
+}
+
+} //namespace Slic3r
+

src/libslic3r/ShortEdgeCollapse.hpp

@@ -0,0 +1,16 @@
+#ifndef SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_
+#define SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_
+
+#include "libslic3r/TriangleMesh.hpp"
+
+namespace Slic3r{
+
+// Decimates the model by collapsing short edges. It starts with very small edges and gradually increases the collapsible length,
+// until the target triangle count is reached (the algorithm will certainly undershoot the target count, result will have less triangles than target count)
+//  The algorithm does not check for triangle flipping, disconnections, self intersections or any other degeneration that can appear during mesh processing.
+void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count);
+
+}
+
+
+#endif /* SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_ */

src/libslic3r/TreeSupport.cpp

@@ -1421,6 +1421,15 @@ void TreeSupport::generate_toolpaths()
     bool obj_is_vertical = obj_size.x() < obj_size.y();
     int num_layers_to_change_infill_direction = int(HEIGHT_TO_SWITCH_INFILL_DIRECTION / object_config.layer_height.value);  // change direction every 30mm
 
+    std::shared_ptr<Fill> filler_interface = std::shared_ptr<Fill>(Fill::new_from_type(m_support_params.contact_fill_pattern));
+    std::shared_ptr<Fill> filler_Roof1stLayer = std::shared_ptr<Fill>(Fill::new_from_type(ipRectilinear));
+    std::shared_ptr<Fill> filler_support = std::shared_ptr<Fill>(Fill::new_from_type(m_support_params.base_fill_pattern));
+    filler_interface->set_bounding_box(bbox_object);
+    filler_Roof1stLayer->set_bounding_box(bbox_object);
+    filler_support->set_bounding_box(bbox_object);
+    filler_interface->angle = Geometry::deg2rad(object_config.support_angle.value + 90.);//(1 - obj_is_vertical) * M_PI_2;//((1-obj_is_vertical) + int(layer_id / num_layers_to_change_infill_direction)) * M_PI_2;;//layer_id % 2 ? 0 : M_PI_2;
+    filler_Roof1stLayer->angle = Geometry::deg2rad(object_config.support_angle.value + 90.);
+
     // generate tree support tool paths
     tbb::parallel_for(
         tbb::blocked_range<size_t>(m_raft_layers, m_object->tree_support_layer_count()),
@@ -1434,12 +1443,7 @@ void TreeSupport::generate_toolpaths()
 
                 TreeSupportLayer* ts_layer = m_object->get_tree_support_layer(layer_id);
                 Flow support_flow(support_extrusion_width, ts_layer->height, nozzle_diameter);
-                std::unique_ptr<Fill> filler_interface = std::unique_ptr<Fill>(Fill::new_from_type(m_support_params.contact_fill_pattern));
-                std::unique_ptr<Fill> filler_support   = std::unique_ptr<Fill>(Fill::new_from_type(m_support_params.base_fill_pattern));
-                filler_interface->set_bounding_box(bbox_object);
-                filler_support->set_bounding_box(bbox_object);
 
-                filler_interface->angle = Geometry::deg2rad(object_config.support_angle.value + 90.);//(1 - obj_is_vertical) * M_PI_2;//((1-obj_is_vertical) + int(layer_id / num_layers_to_change_infill_direction)) * M_PI_2;;//layer_id % 2 ? 0 : M_PI_2;
 
                 for (auto& area_group : ts_layer->area_groups) {
                     ExPolygon& poly = *area_group.first;
@@ -1465,9 +1469,9 @@ void TreeSupport::generate_toolpaths()
                         // roof_1st_layer
                         fill_params.density = interface_density;
                         // Note: spacing means the separation between two lines as if they are tightly extruded
-                        filler_interface->spacing = m_support_material_interface_flow.spacing();
-                        fill_expolygons_generate_paths(ts_layer->support_fills.entities, std::move(polys), filler_interface.get(), fill_params, erSupportMaterial,
-                                                       m_support_material_interface_flow);                        
+                        filler_Roof1stLayer->spacing = m_support_material_interface_flow.spacing();
+                        fill_expolygons_generate_paths(ts_layer->support_fills.entities, std::move(polys), filler_Roof1stLayer.get(), fill_params, erSupportMaterial,
+                                                       m_support_material_interface_flow);
                     } else if (area_group.second == TreeSupportLayer::FloorType) {
                         // floor_areas
                         fill_params.density = bottom_interface_density;
@@ -2670,7 +2674,7 @@ void TreeSupport::adjust_layer_heights(std::vector<std::vector<Node*>>& contact_
     for (int layer_nr = 1; layer_nr < contact_nodes.size(); layer_nr++) {
         std::vector<Node*>& curr_layer_nodes = contact_nodes[layer_nr];
         for (Node* node : curr_layer_nodes) {
-            if (node->support_roof_layers_below == top_intf_layers || node->support_floor_layers_above == bot_intf_layers) {
+            if (node->support_roof_layers_below >0 || node->support_floor_layers_above == bot_intf_layers) {
                 extremes.push_back(layer_nr);
                 break;
             }

src/libslic3r/TriangleMesh.cpp

@@ -9,7 +9,7 @@
 #include "Execution/ExecutionTBB.hpp"
 #include "Execution/ExecutionSeq.hpp"
 #include "Utils.hpp"
-
+#include "Format/STL.hpp"
 #include <libqhullcpp/Qhull.h>
 #include <libqhullcpp/QhullFacetList.h>
 #include <libqhullcpp/QhullVertexSet.h>
@@ -208,10 +208,10 @@ bool TriangleMesh::from_stl(stl_file& stl, bool repair)
     return true;
 }
 
-bool TriangleMesh::ReadSTLFile(const char* input_file, bool repair)
+bool TriangleMesh::ReadSTLFile(const char* input_file, bool repair, ImportstlProgressFn stlFn)
 { 
     stl_file stl;
-    if (! stl_open(&stl, input_file))
+    if (! stl_open(&stl, input_file, stlFn))
         return false;
     return from_stl(stl, repair);
 }

src/libslic3r/TriangleMesh.hpp

@@ -10,7 +10,7 @@
 #include "Point.hpp"
 #include "Polygon.hpp"
 #include "ExPolygon.hpp"
-
+#include "Format/STL.hpp"
 namespace Slic3r {
 
 class TriangleMesh;
@@ -94,7 +94,7 @@ public:
     explicit TriangleMesh(indexed_triangle_set &&M, const RepairedMeshErrors& repaired_errors = RepairedMeshErrors());
     void clear() { this->its.clear(); this->m_stats.clear(); }
     bool from_stl(stl_file& stl, bool repair = true);
-    bool ReadSTLFile(const char* input_file, bool repair = true);
+    bool ReadSTLFile(const char* input_file, bool repair = true, ImportstlProgressFn stlFn = nullptr);
     bool write_ascii(const char* output_file);
     bool write_binary(const char* output_file);
     float volume();
@@ -152,10 +152,14 @@ public:
 
     const TriangleMeshStats& stats() const { return m_stats; }
 
+    void set_init_shift(const Vec3d &offset) { m_init_shift = offset; }
+    Vec3d get_init_shift() const { return m_init_shift; }
+
     indexed_triangle_set its;
 
 private:
     TriangleMeshStats m_stats;
+    Vec3d m_init_shift {0.0, 0.0, 0.0};
 };
 
 // Index of face indices incident with a vertex index.

src/libslic3r/TriangleSelector.cpp

@@ -1427,7 +1427,7 @@ void TriangleSelector::get_facets(std::vector<indexed_triangle_set>& facets_per_
 {
     facets_per_type.clear();
 
-    for (int type = (int)EnforcerBlockerType::NONE; type < (int)EnforcerBlockerType::Extruder15; type++) {
+    for (int type = (int)EnforcerBlockerType::NONE; type <= (int)EnforcerBlockerType::ExtruderMax; type++) {
         facets_per_type.emplace_back();
         indexed_triangle_set& its = facets_per_type.back();
         std::vector<int> vertex_map(m_vertices.size(), -1);

src/libslic3r/TriangleSelector.hpp

@@ -7,12 +7,10 @@
 #include <cfloat>
 #include "Point.hpp"
 #include "TriangleMesh.hpp"
+#include "libslic3r/Model.hpp"
 
 namespace Slic3r {
 
-enum class EnforcerBlockerType : int8_t;
-
-
 // Following class holds information about selected triangles. It also has power
 // to recursively subdivide the triangles and make the selection finer.
 class TriangleSelector
@@ -275,7 +273,7 @@ public:
     std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> serialize() const;
 
     // Load serialized data. Assumes that correct mesh is loaded.
-    void deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>>& data, bool needs_reset = true, EnforcerBlockerType max_ebt = (EnforcerBlockerType)std::numeric_limits<int8_t>::max());
+    void deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>>& data, bool needs_reset = true, EnforcerBlockerType max_ebt = EnforcerBlockerType::ExtruderMax);
 
     // For all triangles, remove the flag indicating that the triangle was selected by seed fill.
     void seed_fill_unselect_all_triangles();
@@ -312,7 +310,7 @@ protected:
         void set_division(int sides_to_split, int special_side_idx);
 
         // Get/set current state.
-        void set_state(EnforcerBlockerType type) { assert(! is_split()); state = type; }
+        void set_state(EnforcerBlockerType type) { assert(!is_split()); state = type; }
         EnforcerBlockerType get_state() const { assert(! is_split()); return state; }
 
         // Set if the triangle has been selected or unselected by seed fill.

src/libslic3r/TriangleSetSampling.cpp

@@ -0,0 +1,71 @@
+#include "TriangleSetSampling.hpp"
+#include <map>
+#include <random>
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+
+namespace Slic3r {
+
+TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set) {
+    std::vector<double> triangles_area(triangle_set.indices.size());
+
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, triangle_set.indices.size()),
+            [&triangle_set, &triangles_area](
+                    tbb::blocked_range<size_t> r) {
+                for (size_t t_idx = r.begin(); t_idx < r.end(); ++t_idx) {
+                    const Vec3f &a = triangle_set.vertices[triangle_set.indices[t_idx].x()];
+                    const Vec3f &b = triangle_set.vertices[triangle_set.indices[t_idx].y()];
+                    const Vec3f &c = triangle_set.vertices[triangle_set.indices[t_idx].z()];
+                    double area = double(0.5 * (b - a).cross(c - a).norm());
+                    triangles_area[t_idx] = area;
+                }
+            });
+
+    std::map<double, size_t> area_sum_to_triangle_idx;
+    float area_sum = 0;
+    for (size_t t_idx = 0; t_idx < triangles_area.size(); ++t_idx) {
+        area_sum += triangles_area[t_idx];
+        area_sum_to_triangle_idx[area_sum] = t_idx;
+    }
+
+    std::mt19937_64 mersenne_engine { 27644437 };
+    // random numbers on interval [0, 1)
+    std::uniform_real_distribution<double> fdistribution;
+
+    auto get_random = [&fdistribution, &mersenne_engine]() {
+        return Vec3d { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) };
+    };
+
+    std::vector<Vec3d> random_samples(samples_count);
+    std::generate(random_samples.begin(), random_samples.end(), get_random);
+
+    TriangleSetSamples result;
+    result.total_area = area_sum;
+    result.positions.resize(samples_count);
+    result.normals.resize(samples_count);
+    result.triangle_indices.resize(samples_count);
+
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, samples_count),
+            [&triangle_set, &area_sum_to_triangle_idx, &area_sum, &random_samples, &result](
+                    tbb::blocked_range<size_t> r) {
+                for (size_t s_idx = r.begin(); s_idx < r.end(); ++s_idx) {
+                    double t_sample = random_samples[s_idx].x() * area_sum;
+                    size_t t_idx = area_sum_to_triangle_idx.upper_bound(t_sample)->second;
+
+                    double sq_u = std::sqrt(random_samples[s_idx].y());
+                    double v = random_samples[s_idx].z();
+
+                    Vec3f A = triangle_set.vertices[triangle_set.indices[t_idx].x()];
+                    Vec3f B = triangle_set.vertices[triangle_set.indices[t_idx].y()];
+                    Vec3f C = triangle_set.vertices[triangle_set.indices[t_idx].z()];
+
+                    result.positions[s_idx] = A * (1 - sq_u) + B * (sq_u * (1 - v)) + C * (v * sq_u);
+                    result.normals[s_idx] = ((B - A).cross(C - B)).normalized();
+                    result.triangle_indices[s_idx] = t_idx;
+                }
+            });
+
+    return result;
+}
+
+}

src/libslic3r/TriangleSetSampling.hpp

@@ -0,0 +1,20 @@
+#ifndef SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_
+#define SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_
+
+#include <admesh/stl.h>
+#include "libslic3r/Point.hpp"
+
+namespace Slic3r {
+
+struct TriangleSetSamples {
+    float total_area;
+    std::vector<Vec3f> positions;
+    std::vector<Vec3f> normals;
+    std::vector<size_t> triangle_indices;
+};
+
+TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set);
+
+}
+
+#endif /* SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_ */

src/libslic3r/Utils.hpp

@@ -283,6 +283,7 @@ inline typename CONTAINER_TYPE::value_type& next_value_modulo(typename CONTAINER
 }
 
 extern std::string xml_escape(std::string text, bool is_marked = false);
+extern std::string xml_unescape(std::string text);
 
 
 #if defined __GNUC__ && __GNUC__ < 5 && !defined __clang__

src/libslic3r/libslic3r.h

@@ -24,6 +24,13 @@
 #include <cmath>
 #include <type_traits>
 
+#ifdef _WIN32
+// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
+// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740
+// Thus it is recommended to use boost::container::deque instead.
+#include <boost/container/deque.hpp>
+#endif // _WIN32
+
 #include "Technologies.hpp"
 #include "Semver.hpp"
 
@@ -96,6 +103,16 @@ namespace Slic3r {
 
 extern Semver SEMVER;
 
+// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
+template<class T, class Allocator = std::allocator<T>>
+using deque =
+#ifdef _WIN32
+    // Use boost implementation, which allocates blocks of 512 bytes instead of blocks of 8 bytes.
+    boost::container::deque<T, Allocator>;
+#else  // _WIN32
+    std::deque<T, Allocator>;
+#endif // _WIN32
+
 template<typename T, typename Q>
 inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); }
 
@@ -341,6 +358,12 @@ public:
     inline bool   empty() const { return size() == 0; }
 };
 
+template<class T, class = FloatingOnly<T>>
+constexpr T NaN = std::numeric_limits<T>::quiet_NaN();
+
+constexpr float  NaNf = NaN<float>;
+constexpr double NaNd = NaN<double>;
+
 } // namespace Slic3r
 
 #endif

src/libslic3r/utils.cpp

@@ -1182,6 +1182,42 @@ std::string xml_escape(std::string text, bool is_marked/* = false*/)
     return text;
 }
 
+std::string xml_unescape(std::string s)
+{
+	std::string ret;
+	std::string::size_type i = 0;
+	std::string::size_type pos = 0;
+	while (i < s.size()) {
+		std::string rep;
+		if (s[i] == '&') {
+			if (s.substr(i, 4) == "&lt;") {
+				ret += s.substr(pos, i - pos) + "<";
+				i += 4;
+				pos = i;
+			}
+			else if (s.substr(i, 4) == "&gt;") {
+				ret += s.substr(pos, i - pos) + ">";
+				i += 4;
+				pos = i;
+			}
+			else if (s.substr(i, 5) == "&amp;") {
+				ret += s.substr(pos, i - pos) + "&";
+				i += 5;
+				pos = i;
+			}
+			else {
+				++i;
+			}
+		}
+		else {
+			++i;
+		}
+	}
+
+	ret += s.substr(pos);
+	return ret;
+}
+
 std::string format_memsize_MB(size_t n)
 {
     std::string out;