Seam placement improvements.

Previously there was an algirithm that was fitting a curve through
the resulting seam points. This worked for smooth models, but was
failing for some very basic cases (such as a cylinder).

The new algorithm builds on the previously implemented visibility
algirithm but does not do the curve fitting anymore.

Now the code is more separated for the four seam placement
options (rear, random, aligned, nearest).

Nearest and random are handled as one would expect.
Aligned is handled in a more "brute force" manner
(trying multiple seams and picking the best one)
and rear smartly switches between two modes
(center straight line projection and max y).

src/libslic3r/CMakeLists.txt

@@ -177,6 +177,24 @@ set(SLIC3R_SOURCES
     GCode/SpiralVase.hpp
     GCode/SeamPlacer.cpp
     GCode/SeamPlacer.hpp
+    GCode/SeamChoice.cpp
+    GCode/SeamChoice.hpp
+    GCode/SeamPerimeters.cpp
+    GCode/SeamPerimeters.hpp
+    GCode/SeamShells.cpp
+    GCode/SeamShells.hpp
+    GCode/SeamGeometry.cpp
+    GCode/SeamGeometry.hpp
+    GCode/SeamAligned.cpp
+    GCode/SeamAligned.hpp
+    GCode/SeamRear.cpp
+    GCode/SeamRear.hpp
+    GCode/SeamRandom.cpp
+    GCode/SeamRandom.hpp
+    GCode/SeamPainting.cpp
+    GCode/SeamPainting.hpp
+    GCode/ModelVisibility.cpp
+    GCode/ModelVisibility.hpp
     GCode/SmoothPath.cpp
     GCode/SmoothPath.hpp
     GCode/ToolOrdering.cpp

src/libslic3r/GCode.cpp

@@ -1232,7 +1232,9 @@ void GCodeGenerator::_do_export(Print& print, GCodeOutputStream &file, Thumbnail
 
     // Collect custom seam data from all objects.
     std::function<void(void)> throw_if_canceled_func = [&print]() { print.throw_if_canceled();};
-    m_seam_placer.init(print, throw_if_canceled_func);
+
+    const Seams::Params params{Seams::Placer::get_params(print.full_print_config())};
+    m_seam_placer.init(print.objects(), params, throw_if_canceled_func);
 
     if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming)) {
         // Set initial extruder only after custom start G-code.
@@ -2953,7 +2955,7 @@ std::string GCodeGenerator::extrude_loop(const ExtrusionLoop &loop_src, const GC
     Point seam_point = this->last_position.has_value() ? *this->last_position : Point::Zero();
     if (!m_config.spiral_vase && comment_is_perimeter(description)) {
         assert(m_layer != nullptr);
-        seam_point = m_seam_placer.place_seam(m_layer, loop_src, m_config.external_perimeters_first, seam_point);
+        seam_point = m_seam_placer.place_seam(m_layer, loop_src, seam_point);
     }
     // Because the G-code export has 1um resolution, don't generate segments shorter than 1.5 microns,
     // thus empty path segments will not be produced by G-code export.

src/libslic3r/GCode.hpp

@@ -353,8 +353,7 @@ private:
     std::string     set_extruder(unsigned int extruder_id, double print_z);
     bool line_distancer_is_required(const std::vector<unsigned int>& extruder_ids);
 
-    // Cache for custom seam enforcers/blockers for each layer.
-    SeamPlacer                          m_seam_placer;
+    Seams::Placer                       m_seam_placer;
 
     /* Origin of print coordinates expressed in unscaled G-code coordinates.
        This affects the input arguments supplied to the extrude*() and travel_to()

src/libslic3r/GCode/ModelVisibility.cpp

@@ -0,0 +1,301 @@
+#include <boost/log/trivial.hpp>
+
+#include "libslic3r/ShortEdgeCollapse.hpp"
+#include "libslic3r/GCode/ModelVisibility.hpp"
+#include "libslic3r/AABBTreeIndirect.hpp"
+
+namespace Slic3r::ModelInfo {
+namespace Impl {
+
+CoordinateFunctor::CoordinateFunctor(const std::vector<Vec3f> *coords) : coordinates(coords) {}
+CoordinateFunctor::CoordinateFunctor() : coordinates(nullptr) {}
+
+const float &CoordinateFunctor::operator()(size_t idx, size_t dim) const {
+    return coordinates->operator[](idx)[dim];
+}
+
+
+template<typename T> int sgn(T val) {
+    return int(T(0) < val) - int(val < T(0));
+}
+
+/// Coordinate frame
+class Frame {
+public:
+    Frame() {
+        mX = Vec3f(1, 0, 0);
+        mY = Vec3f(0, 1, 0);
+        mZ = Vec3f(0, 0, 1);
+    }
+
+    Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) :
+            mX(x), mY(y), mZ(z) {
+    }
+
+    void set_from_z(const Vec3f &z) {
+        mZ = z.normalized();
+        Vec3f tmpZ = mZ;
+        Vec3f tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3f(0, 1, 0) : Vec3f(1, 0, 0);
+        mY = (tmpZ.cross(tmpX)).normalized();
+        mX = mY.cross(tmpZ);
+    }
+
+    Vec3f to_world(const Vec3f &a) const {
+        return a.x() * mX + a.y() * mY + a.z() * mZ;
+    }
+
+    Vec3f to_local(const Vec3f &a) const {
+        return Vec3f(mX.dot(a), mY.dot(a), mZ.dot(a));
+    }
+
+    const Vec3f& binormal() const {
+        return mX;
+    }
+
+    const Vec3f& tangent() const {
+        return mY;
+    }
+
+    const Vec3f& normal() const {
+        return mZ;
+    }
+
+private:
+    Vec3f mX, mY, mZ;
+};
+
+Vec3f sample_sphere_uniform(const Vec2f &samples) {
+    float term1 = 2.0f * float(PI) * samples.x();
+    float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
+    return {cos(term1) * term2, sin(term1) * term2,
+        1.0f - 2.0f * samples.y()};
+}
+
+Vec3f sample_hemisphere_uniform(const Vec2f &samples) {
+    float term1 = 2.0f * float(PI) * samples.x();
+    float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
+    return {cos(term1) * term2, sin(term1) * term2,
+        abs(1.0f - 2.0f * samples.y())};
+}
+
+Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) {
+    float term1 = 2.f * float(PI) * samples.x();
+    float term2 = pow(samples.y(), 1.f / (power + 1.f));
+    float term3 = sqrt(1.f - term2 * term2);
+
+    return Vec3f(cos(term1) * term3, sin(term1) * term3, term2);
+}
+
+std::vector<float> raycast_visibility(
+    const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
+    const indexed_triangle_set &triangles,
+    const TriangleSetSamples &samples,
+    size_t negative_volumes_start_index,
+    const Visibility::Params &params
+) {
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: raycast visibility of " << samples.positions.size() << " samples over " << triangles.indices.size()
+            << " triangles: end";
+
+    //prepare uniform samples of a hemisphere
+    float step_size = 1.0f / params.sqr_rays_per_sample_point;
+    std::vector<Vec3f> precomputed_sample_directions(
+            params.sqr_rays_per_sample_point * params.sqr_rays_per_sample_point);
+    for (size_t x_idx = 0; x_idx < params.sqr_rays_per_sample_point; ++x_idx) {
+        float sample_x = x_idx * step_size + step_size / 2.0;
+        for (size_t y_idx = 0; y_idx < params.sqr_rays_per_sample_point; ++y_idx) {
+            size_t dir_index = x_idx * params.sqr_rays_per_sample_point + y_idx;
+            float sample_y = y_idx * step_size + step_size / 2.0;
+            precomputed_sample_directions[dir_index] = sample_hemisphere_uniform( { sample_x, sample_y });
+        }
+    }
+
+    bool model_contains_negative_parts = negative_volumes_start_index < triangles.indices.size();
+
+    std::vector<float> result(samples.positions.size());
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, result.size()),
+            [&triangles, &precomputed_sample_directions, model_contains_negative_parts, negative_volumes_start_index,
+                    &raycasting_tree, &result, &samples, &params](tbb::blocked_range<size_t> r) {
+                // Maintaining hits memory outside of the loop, so it does not have to be reallocated for each query.
+                std::vector<igl::Hit> hits;
+                for (size_t s_idx = r.begin(); s_idx < r.end(); ++s_idx) {
+                    result[s_idx] = 1.0f;
+                    const float decrease_step = 1.0f
+                            / (params.sqr_rays_per_sample_point * params.sqr_rays_per_sample_point);
+
+                    const Vec3f &center = samples.positions[s_idx];
+                    const Vec3f &normal = samples.normals[s_idx];
+                    // apply the local direction via Frame struct - the local_dir is with respect to +Z being forward
+                    Frame f;
+                    f.set_from_z(normal);
+
+                    for (const auto &dir : precomputed_sample_directions) {
+                        Vec3f final_ray_dir = (f.to_world(dir));
+                        if (!model_contains_negative_parts) {
+                            igl::Hit hitpoint;
+                            // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and
+                            // direction.
+                            Vec3d final_ray_dir_d = final_ray_dir.cast<double>();
+                            Vec3d ray_origin_d = (center + normal * 0.01f).cast<double>(); // start above surface.
+                            bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices,
+                                    triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint);
+                            if (hit && its_face_normal(triangles, hitpoint.id).dot(final_ray_dir) <= 0) {
+                                result[s_idx] -= decrease_step;
+                            }
+                        } else { //TODO improve logic for order based boolean operations - consider order of volumes
+                            bool casting_from_negative_volume = samples.triangle_indices[s_idx]
+                                    >= negative_volumes_start_index;
+
+                            Vec3d ray_origin_d = (center + normal * 0.01f).cast<double>(); // start above surface.
+                            if (casting_from_negative_volume) { // if casting from negative volume face, invert direction, change start pos
+                                final_ray_dir = -1.0 * final_ray_dir;
+                                ray_origin_d = (center - normal * 0.01f).cast<double>();
+                            }
+                            Vec3d final_ray_dir_d = final_ray_dir.cast<double>();
+                            bool some_hit = AABBTreeIndirect::intersect_ray_all_hits(triangles.vertices,
+                                    triangles.indices, raycasting_tree,
+                                    ray_origin_d, final_ray_dir_d, hits);
+                            if (some_hit) {
+                                int counter = 0;
+                                // NOTE: iterating in reverse, from the last hit for one simple reason: We know the state of the ray at that point;
+                                //  It cannot be inside model, and it cannot be inside negative volume
+                                for (int hit_index = int(hits.size()) - 1; hit_index >= 0; --hit_index) {
+                                    Vec3f face_normal = its_face_normal(triangles, hits[hit_index].id);
+                                    if (hits[hit_index].id >= int(negative_volumes_start_index)) { //negative volume hit
+                                        counter -= sgn(face_normal.dot(final_ray_dir)); // if volume face aligns with ray dir, we are leaving negative space
+                                        // which in reverse hit analysis means, that we are entering negative space :) and vice versa
+                                    } else {
+                                        counter += sgn(face_normal.dot(final_ray_dir));
+                                    }
+                                }
+                                if (counter == 0) {
+                                    result[s_idx] -= decrease_step;
+                                }
+                            }
+                        }
+                    }
+                }
+            });
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: raycast visibility of " << samples.positions.size() << " samples over " << triangles.indices.size()
+            << " triangles: end";
+
+    return result;
+}
+}
+
+Visibility::Visibility(
+    const Transform3d &obj_transform,
+    const ModelVolumePtrs &volumes,
+    const Params &params,
+    const std::function<void(void)> &throw_if_canceled
+) {
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: gather occlusion meshes: start";
+    indexed_triangle_set triangle_set;
+    indexed_triangle_set negative_volumes_set;
+    //add all parts
+    for (const ModelVolume *model_volume : volumes) {
+        if (model_volume->type() == ModelVolumeType::MODEL_PART
+                || model_volume->type() == ModelVolumeType::NEGATIVE_VOLUME) {
+            auto model_transformation = model_volume->get_matrix();
+            indexed_triangle_set model_its = model_volume->mesh().its;
+            its_transform(model_its, model_transformation);
+            if (model_volume->type() == ModelVolumeType::MODEL_PART) {
+                its_merge(triangle_set, model_its);
+            } else {
+                its_merge(negative_volumes_set, model_its);
+            }
+        }
+    }
+    throw_if_canceled();
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: gather occlusion meshes: end";
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: decimate: start";
+    its_short_edge_collpase(triangle_set, params.fast_decimation_triangle_count_target);
+    its_short_edge_collpase(negative_volumes_set, params.fast_decimation_triangle_count_target);
+
+    size_t negative_volumes_start_index = triangle_set.indices.size();
+    its_merge(triangle_set, negative_volumes_set);
+    its_transform(triangle_set, obj_transform);
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: decimate: end";
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: Compute visibility sample points: start";
+
+    this->mesh_samples = sample_its_uniform_parallel(params.raycasting_visibility_samples_count,
+            triangle_set);
+    this->mesh_samples_coordinate_functor = Impl::CoordinateFunctor(&this->mesh_samples.positions);
+    this->mesh_samples_tree = KDTreeIndirect<3, float, Impl::CoordinateFunctor>(this->mesh_samples_coordinate_functor,
+            this->mesh_samples.positions.size());
+
+    // The following code determines search area for random visibility samples on the mesh when calculating visibility of each perimeter point
+    // number of random samples in the given radius (area) is approximately poisson distribution
+    // to compute ideal search radius (area), we use exponential distribution (complementary distr to poisson)
+    // parameters of exponential distribution to compute area that will have with probability="probability" more than given number of samples="samples"
+    float probability = 0.9f;
+    float samples = 4;
+    float density = params.raycasting_visibility_samples_count / this->mesh_samples.total_area;
+    // exponential probability distrubtion function is : f(x) = P(X > x) = e^(l*x) where l is the rate parameter (computed as 1/u where u is mean value)
+    // probability that sampled area A with S samples contains more than samples count:
+    //  P(S > samples in A) = e^-(samples/(density*A));   express A:
+    float search_area = samples / (-logf(probability) * density);
+    float search_radius = sqrt(search_area / PI);
+    this->mesh_samples_radius = search_radius;
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: Compute visiblity sample points: end";
+    throw_if_canceled();
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: Mesh sample raidus: " << this->mesh_samples_radius;
+
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: build AABB tree: start";
+    auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices,
+            triangle_set.indices);
+
+    throw_if_canceled();
+    BOOST_LOG_TRIVIAL(debug)
+    << "SeamPlacer: build AABB tree: end";
+    this->mesh_samples_visibility = Impl::raycast_visibility(raycasting_tree, triangle_set, this->mesh_samples,
+            negative_volumes_start_index, params);
+    throw_if_canceled();
+}
+
+float Visibility::calculate_point_visibility(const Vec3f &position) const {
+    std::vector<size_t> points = find_nearby_points(mesh_samples_tree, position, mesh_samples_radius);
+    if (points.empty()) {
+        return 1.0f;
+    }
+
+    auto compute_dist_to_plane = [](const Vec3f &position, const Vec3f &plane_origin,
+                                    const Vec3f &plane_normal) {
+        Vec3f orig_to_point = position - plane_origin;
+        return std::abs(orig_to_point.dot(plane_normal));
+    };
+
+    float total_weight = 0;
+    float total_visibility = 0;
+    for (size_t i = 0; i < points.size(); ++i) {
+        size_t sample_idx = points[i];
+
+        Vec3f sample_point = this->mesh_samples.positions[sample_idx];
+        Vec3f sample_normal = this->mesh_samples.normals[sample_idx];
+
+        float weight = mesh_samples_radius -
+            compute_dist_to_plane(position, sample_point, sample_normal);
+        weight += (mesh_samples_radius - (position - sample_point).norm());
+        total_visibility += weight * mesh_samples_visibility[sample_idx];
+        total_weight += weight;
+    }
+
+    return total_visibility / total_weight;
+}
+
+}

src/libslic3r/GCode/ModelVisibility.hpp

@@ -0,0 +1,51 @@
+#ifndef libslic3r_ModelVisibility_hpp_
+#define libslic3r_ModelVisibility_hpp_
+
+#include "libslic3r/KDTreeIndirect.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/Model.hpp"
+#include "libslic3r/TriangleSetSampling.hpp"
+
+namespace Slic3r::ModelInfo {
+namespace Impl {
+
+struct CoordinateFunctor
+{
+    const std::vector<Vec3f> *coordinates;
+    CoordinateFunctor(const std::vector<Vec3f> *coords);
+    CoordinateFunctor();
+
+    const float &operator()(size_t idx, size_t dim) const;
+};
+} // namespace Impl
+
+struct Visibility
+{
+    struct Params
+    {
+        // Number of samples generated on the mesh. There are
+        // sqr_rays_per_sample_point*sqr_rays_per_sample_point rays casted from each samples
+        size_t raycasting_visibility_samples_count{};
+        size_t fast_decimation_triangle_count_target{};
+        // square of number of rays per sample point
+        size_t sqr_rays_per_sample_point{};
+    };
+
+    Visibility(
+        const Transform3d &obj_transform,
+        const ModelVolumePtrs &volumes,
+        const Params &params,
+        const std::function<void(void)> &throw_if_canceled
+    );
+
+    TriangleSetSamples mesh_samples;
+    std::vector<float> mesh_samples_visibility;
+    Impl::CoordinateFunctor mesh_samples_coordinate_functor;
+    KDTreeIndirect<3, float, Impl::CoordinateFunctor> mesh_samples_tree{Impl::CoordinateFunctor{}};
+    float mesh_samples_radius;
+
+    float calculate_point_visibility(const Vec3f &position) const;
+};
+
+} // namespace Slic3r::ModelInfo
+#endif // libslic3r_ModelVisibility_hpp_

src/libslic3r/GCode/SeamAligned.cpp

@@ -0,0 +1,491 @@
+#include "libslic3r/GCode/SeamAligned.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
+#include "libslic3r/GCode/ModelVisibility.hpp"
+#include <fstream>
+
+namespace Slic3r::Seams::Aligned {
+using Perimeters::PointType;
+using Perimeters::PointClassification;
+
+namespace Impl {
+const Perimeters::Perimeter::PointTrees &pick_trees(
+    const Perimeters::Perimeter &perimeter, const PointType point_type
+) {
+    switch (point_type) {
+    case PointType::enforcer: return perimeter.enforced_points;
+    case PointType::blocker: return perimeter.blocked_points;
+    case PointType::common: return perimeter.common_points;
+    }
+    throw std::runtime_error("Point trees for point type do not exist.");
+}
+
+const Perimeters::Perimeter::OptionalPointTree &pick_tree(
+    const Perimeters::Perimeter::PointTrees &point_trees,
+    const PointClassification &point_classification
+) {
+    switch (point_classification) {
+    case PointClassification::overhang: return point_trees.overhanging_points;
+    case PointClassification::embedded: return point_trees.embedded_points;
+    case PointClassification::common: return point_trees.common_points;
+    }
+    throw std::runtime_error("Point tree for classification does not exist.");
+}
+
+unsigned point_value(PointType point_type, PointClassification point_classification) {
+    // Better be explicit than smart.
+    switch (point_type) {
+    case PointType::enforcer:
+        switch (point_classification) {
+        case PointClassification::embedded: return 9;
+        case PointClassification::common: return 8;
+        case PointClassification::overhang: return 7;
+        }
+    case PointType::common:
+        switch (point_classification) {
+        case PointClassification::embedded: return 6;
+        case PointClassification::common: return 5;
+        case PointClassification::overhang: return 4;
+        }
+    case PointType::blocker:
+        switch (point_classification) {
+        case PointClassification::embedded: return 3;
+        case PointClassification::common: return 2;
+        case PointClassification::overhang: return 1;
+        }
+    }
+    return 0;
+}
+
+SeamChoice pick_seam_option(const Perimeters::Perimeter &perimeter, const SeamOptions &options) {
+    const std::vector<PointType> &types{perimeter.point_types};
+    const std::vector<PointClassification> &classifications{perimeter.point_classifications};
+    const std::vector<Vec2d> &positions{perimeter.positions};
+
+    unsigned closeset_point_value =
+        point_value(types.at(options.closest), classifications[options.closest]);
+
+    if (options.snapped) {
+        unsigned snapped_point_value =
+            point_value(types.at(*options.snapped), classifications[*options.snapped]);
+        if (snapped_point_value >= closeset_point_value) {
+            const Vec2d position{positions.at(*options.snapped)};
+            return {*options.snapped, *options.snapped, position};
+        }
+    }
+
+    unsigned adjacent_point_value =
+        point_value(types.at(options.adjacent), classifications[options.adjacent]);
+    if (adjacent_point_value < closeset_point_value) {
+        const Vec2d position = positions[options.closest];
+        return {options.closest, options.closest, position};
+    }
+
+    const std::size_t next_index{options.adjacent_forward ? options.adjacent : options.closest};
+    const std::size_t previous_index{options.adjacent_forward ? options.closest : options.adjacent};
+    return {previous_index, next_index, options.on_edge};
+}
+
+std::optional<std::size_t> snap_to_angle(
+    const Vec2d &point,
+    const std::size_t search_start,
+    const Perimeters::Perimeter &perimeter,
+    const double max_detour
+) {
+    using Perimeters::AngleType;
+    const std::vector<Vec2d> &positions{perimeter.positions};
+    const std::vector<AngleType> &angle_types{perimeter.angle_types};
+
+    std::optional<std::size_t> match;
+    double min_distance{std::numeric_limits<double>::infinity()};
+    AngleType angle_type{AngleType::convex};
+
+    const auto visitor{[&](const std::size_t index) {
+        const double distance = (positions[index] - point).norm();
+        if (distance > max_detour) {
+            return true;
+        }
+        if (angle_types[index] == angle_type &&
+            distance < min_distance) {
+            match = index;
+            min_distance = distance;
+            return true;
+        }
+        return false;
+    }};
+    Geometry::visit_backward(search_start, positions.size(), visitor);
+    Geometry::visit_forward(search_start, positions.size(), visitor);
+    if (match) {
+        return match;
+    }
+
+    min_distance = std::numeric_limits<double>::infinity();
+    angle_type = AngleType::concave;
+
+    Geometry::visit_backward(search_start, positions.size(), visitor);
+    Geometry::visit_forward(search_start, positions.size(), visitor);
+
+    return match;
+}
+
+SeamOptions get_seam_options(
+    const Perimeters::Perimeter &perimeter,
+    const Vec2d &prefered_position,
+    const Perimeters::Perimeter::PointTree &points_tree,
+    const double max_detour
+) {
+    const std::vector<Vec2d> &positions{perimeter.positions};
+
+    const std::size_t closest{find_closest_point(points_tree, prefered_position.head<2>())};
+    std::size_t previous{closest == 0 ? positions.size() - 1 : closest - 1};
+    std::size_t next{closest == positions.size() - 1 ? 0 : closest + 1};
+
+    const Vec2d previous_adjacent_point{positions[previous]};
+    const Vec2d closest_point{positions[closest]};
+    const Vec2d next_adjacent_point{positions[next]};
+
+    const Linef previous_segment{previous_adjacent_point, closest_point};
+    const auto [previous_point, previous_distance] =
+        Geometry::distance_to_segment_squared(previous_segment, prefered_position);
+    const Linef next_segment{closest_point, next_adjacent_point};
+    const auto [next_point, next_distance] =
+        Geometry::distance_to_segment_squared(next_segment, prefered_position);
+
+    const bool adjacent_forward{next_distance < previous_distance};
+    const Vec2d nearest_point{adjacent_forward ? next_point : previous_point};
+    const std::size_t adjacent{adjacent_forward ? next : previous};
+
+    std::optional<std::size_t> snapped{
+        snap_to_angle(nearest_point.head<2>(), closest, perimeter, max_detour)};
+
+    return {
+        closest, adjacent, adjacent_forward, snapped, nearest_point,
+    };
+}
+
+std::optional<SeamChoice> LeastVisible::operator()(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) const {
+    std::optional<size_t> chosen_index;
+    double visibility{std::numeric_limits<double>::infinity()};
+
+    for (std::size_t i{0}; i < perimeter.positions.size(); ++i) {
+        if (perimeter.point_types[i] != point_type ||
+            perimeter.point_classifications[i] != point_classification) {
+            continue;
+        }
+        const Vec2d point{perimeter.positions[i]};
+        const double point_visibility{precalculated_visibility[i]};
+
+        if (point_visibility < visibility) {
+            visibility = point_visibility;
+            chosen_index = i;
+        }
+    }
+
+    if (chosen_index) {
+        return {{*chosen_index, *chosen_index, perimeter.positions[*chosen_index]}};
+    }
+    return std::nullopt;
+}
+
+std::optional<SeamChoice> Nearest::operator()(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) const {
+    const Perimeters::Perimeter::PointTrees &trees{pick_trees(perimeter, point_type)};
+    const Perimeters::Perimeter::OptionalPointTree &tree = pick_tree(trees, point_classification);
+    if (tree) {
+        const SeamOptions options{get_seam_options(perimeter, prefered_position, *tree, max_detour)};
+        return pick_seam_option(perimeter, options);
+    }
+    return std::nullopt;
+}
+} // namespace Impl
+
+double VisibilityCalculator::operator()(
+    const SeamChoice &choice, const Perimeters::Perimeter &perimeter
+) const {
+    double visibility = points_visibility.calculate_point_visibility(
+        to_3d(choice.position, perimeter.slice_z).cast<float>()
+    );
+
+    const double angle{
+        choice.previous_index == choice.next_index ? perimeter.angles[choice.previous_index] : 0.0};
+    visibility +=
+        get_angle_visibility_modifier(angle, convex_visibility_modifier, concave_visibility_modifier);
+    return visibility;
+}
+
+double VisibilityCalculator::get_angle_visibility_modifier(
+    double angle,
+    const double convex_visibility_modifier,
+    const double concave_visibility_modifier
+) {
+    const double weight_max{angle > 0 ? convex_visibility_modifier : concave_visibility_modifier};
+    angle = std::abs(angle);
+    const double right_angle{M_PI / 2.0};
+    if (angle > right_angle) {
+        return -weight_max;
+    }
+    const double angle_linear_weight{angle / right_angle};
+    // It is smooth and at angle 0 slope is equal to `angle linear weight`, at right angle the slope is 0 and value is equal to weight max.
+    const double angle_smooth_weight{angle / right_angle * weight_max + (right_angle - angle) / right_angle * angle_linear_weight};
+    return -angle_smooth_weight;
+}
+
+std::vector<Vec2d> extract_points(
+    const Perimeters::Perimeter &perimeter, const Perimeters::PointType point_type
+) {
+    std::vector<Vec2d> result;
+    for (std::size_t i{0}; i < perimeter.positions.size(); ++i) {
+        if (perimeter.point_types[i] == point_type) {
+            result.push_back(perimeter.positions[i]);
+        }
+    }
+    return result;
+}
+
+std::vector<Vec2d> get_starting_positions(const Shells::Shell<> &shell) {
+    const Perimeters::Perimeter &perimeter{shell.front().boundary};
+
+    std::vector<Vec2d> enforcers{extract_points(perimeter, Perimeters::PointType::enforcer)};
+    if (!enforcers.empty()) {
+        return enforcers;
+    }
+    std::vector<Vec2d> common{extract_points(perimeter, Perimeters::PointType::common)};
+    if (!common.empty()) {
+        return common;
+    }
+    return perimeter.positions;
+}
+
+struct LeastVisiblePoint
+{
+    SeamChoice choice;
+    double visibility;
+};
+
+struct SeamCandidate {
+    std::vector<SeamChoice> choices;
+    std::vector<double> visibilities;
+};
+
+SeamCandidate get_seam_candidate(
+    const Shells::Shell<> &shell,
+    const Vec2d &starting_position,
+    const SeamChoiceVisibility &visibility_calculator,
+    const Params &params,
+    const std::vector<std::vector<double>> &precalculated_visibility,
+    const std::vector<LeastVisiblePoint> &least_visible_points
+) {
+    using Perimeters::Perimeter, Perimeters::AngleType;
+
+    std::vector<double> choice_visibilities(shell.size(), 1.0);
+    std::vector<SeamChoice> choices{
+        Seams::get_shell_seam(shell, [&, previous_position{starting_position}](const Perimeter &perimeter, std::size_t slice_index) mutable {
+            SeamChoice candidate{Seams::choose_seam_point(
+                perimeter, Impl::Nearest{previous_position, params.max_detour}
+            )};
+            const bool is_too_far{
+                (candidate.position - previous_position).norm() > params.max_detour};
+            const LeastVisiblePoint &least_visible{least_visible_points[slice_index]};
+
+            const bool is_on_edge{
+                candidate.previous_index == candidate.next_index &&
+                perimeter.angle_types[candidate.next_index] != AngleType::smooth};
+
+            if (is_on_edge) {
+                choice_visibilities[slice_index] = precalculated_visibility[slice_index][candidate.previous_index];
+            } else {
+                choice_visibilities[slice_index] =
+                    visibility_calculator(candidate, perimeter);
+            }
+            const bool is_too_visible{
+                choice_visibilities[slice_index] >
+                least_visible.visibility + params.jump_visibility_threshold};
+            const bool can_be_on_edge{
+                !is_on_edge &&
+                perimeter.angle_types[least_visible.choice.next_index] != AngleType::smooth};
+            if (is_too_far || (can_be_on_edge  && is_too_visible)) {
+                candidate = least_visible.choice;
+            }
+            previous_position = candidate.position;
+            return candidate;
+        })};
+    return {std::move(choices), std::move(choice_visibilities)};
+}
+
+using ShellVertexVisibility = std::vector<std::vector<double>>;
+
+std::vector<ShellVertexVisibility> get_shells_vertex_visibility(
+    const Shells::Shells<> &shells, const SeamChoiceVisibility &visibility_calculator
+) {
+    std::vector<ShellVertexVisibility> result;
+
+    result.reserve(shells.size());
+    std::transform(
+        shells.begin(), shells.end(), std::back_inserter(result),
+        [](const Shells::Shell<> &shell) { return ShellVertexVisibility(shell.size()); }
+    );
+
+    Geometry::iterate_nested(shells, [&](const std::size_t shell_index, const std::size_t slice_index) {
+        const Shells::Shell<> &shell{shells[shell_index]};
+        const Shells::Slice<> &slice{shell[slice_index]};
+        const std::vector<Vec2d> &positions{slice.boundary.positions};
+
+        for (std::size_t point_index{0}; point_index < positions.size(); ++point_index) {
+            result[shell_index][slice_index].emplace_back(visibility_calculator(
+                SeamChoice{point_index, point_index, positions[point_index]}, slice.boundary
+            ));
+        }
+    });
+    return result;
+}
+
+using ShellLeastVisiblePoints = std::vector<LeastVisiblePoint>;
+
+std::vector<ShellLeastVisiblePoints> get_shells_least_visible_points(
+    const Shells::Shells<> &shells,
+    const std::vector<ShellVertexVisibility> &precalculated_visibility
+) {
+    std::vector<ShellLeastVisiblePoints> result;
+
+    result.reserve(shells.size());
+    std::transform(
+        shells.begin(), shells.end(), std::back_inserter(result),
+        [](const Shells::Shell<> &shell) { return ShellLeastVisiblePoints(shell.size()); }
+    );
+
+    Geometry::iterate_nested(shells, [&](const std::size_t shell_index, const std::size_t slice_index) {
+        const Shells::Shell<> &shell{shells[shell_index]};
+        const Shells::Slice<> &slice{shell[slice_index]};
+        const SeamChoice least_visibile{
+            Seams::choose_seam_point(slice.boundary, Impl::LeastVisible{precalculated_visibility[shell_index][slice_index]})};
+
+        const double visibility{precalculated_visibility[shell_index][slice_index][least_visibile.previous_index]};
+        result[shell_index][slice_index] = LeastVisiblePoint{least_visibile, visibility};
+    });
+    return result;
+}
+
+using ShellStartingPositions = std::vector<Vec2d>;
+
+std::vector<ShellStartingPositions> get_shells_starting_positions(
+    const Shells::Shells<> &shells
+) {
+    std::vector<ShellStartingPositions> result;
+    for (const Shells::Shell<> &shell : shells) {
+        std::vector<Vec2d> starting_positions{get_starting_positions(shell)};
+        result.push_back(std::move(starting_positions));
+    }
+    return result;
+}
+
+using ShellSeamCandidates = std::vector<SeamCandidate>;
+
+std::vector<ShellSeamCandidates> get_shells_seam_candidates(
+    const Shells::Shells<> &shells,
+    const std::vector<ShellStartingPositions> &starting_positions,
+    const SeamChoiceVisibility &visibility_calculator,
+    const std::vector<ShellVertexVisibility> &precalculated_visibility,
+    const std::vector<ShellLeastVisiblePoints> &least_visible_points,
+    const Params &params
+) {
+    std::vector<ShellSeamCandidates> result;
+
+    result.reserve(starting_positions.size());
+    std::transform(
+        starting_positions.begin(), starting_positions.end(), std::back_inserter(result),
+        [](const ShellStartingPositions &positions) { return ShellSeamCandidates(positions.size()); }
+    );
+
+    Geometry::iterate_nested(starting_positions, [&](const std::size_t shell_index, const std::size_t starting_position_index){
+        const Shells::Shell<> &shell{shells[shell_index]};
+        using Perimeters::Perimeter, Perimeters::AngleType;
+
+        result[shell_index][starting_position_index] = get_seam_candidate(
+            shell,
+            starting_positions[shell_index][starting_position_index],
+            visibility_calculator,
+            params,
+            precalculated_visibility[shell_index],
+            least_visible_points[shell_index]
+        );
+    });
+    return result;
+}
+
+std::vector<SeamChoice> get_shell_seam(
+    const Shells::Shell<> &shell,
+    std::vector<SeamCandidate> seam_candidates
+) {
+    std::vector<SeamChoice> seam;
+    double visibility{std::numeric_limits<double>::infinity()};
+
+    for (std::size_t i{0}; i < seam_candidates.size(); ++i) {
+        using Perimeters::Perimeter, Perimeters::AngleType;
+
+        SeamCandidate seam_candidate{seam_candidates[i]};
+
+        double seam_candidate_visibility{0.0};
+        for (std::size_t slice_index{}; slice_index < shell.size(); ++slice_index) {
+            seam_candidate_visibility += seam_candidate.visibilities[slice_index];
+        }
+
+        if (seam_candidate_visibility < visibility) {
+            seam = std::move(seam_candidate.choices);
+            visibility = seam_candidate_visibility;
+        }
+    }
+
+    return seam;
+}
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const SeamChoiceVisibility &visibility_calculator,
+    const Params &params
+) {
+    const std::vector<ShellVertexVisibility> precalculated_visibility{
+        get_shells_vertex_visibility(shells, visibility_calculator)};
+
+    const std::vector<ShellLeastVisiblePoints> least_visible_points{
+        get_shells_least_visible_points(shells, precalculated_visibility)
+    };
+
+    const std::vector<ShellStartingPositions> starting_positions{
+        get_shells_starting_positions(shells)
+    };
+
+    const std::vector<ShellSeamCandidates> seam_candidates{
+        get_shells_seam_candidates(
+            shells,
+            starting_positions,
+            visibility_calculator,
+            precalculated_visibility,
+            least_visible_points,
+            params
+        )
+    };
+
+    std::vector<std::vector<SeamPerimeterChoice>> layer_seams(get_layer_count(shells));
+
+    for (std::size_t shell_index{0}; shell_index < shells.size(); ++shell_index) {
+        Shells::Shell<> &shell{shells[shell_index]};
+
+        std::vector<SeamChoice> seam{
+            Aligned::get_shell_seam(shell, seam_candidates[shell_index])};
+
+        for (std::size_t perimeter_id{}; perimeter_id < shell.size(); ++perimeter_id) {
+            const SeamChoice &choice{seam[perimeter_id]};
+            Perimeters::Perimeter &perimeter{shell[perimeter_id].boundary};
+            layer_seams[shell[perimeter_id].layer_index].emplace_back(choice, std::move(perimeter));
+        }
+    }
+    return layer_seams;
+}
+
+} // namespace Slic3r::Seams::Aligned

src/libslic3r/GCode/SeamAligned.hpp

@@ -0,0 +1,96 @@
+#ifndef libslic3r_SeamAligned_hpp_
+#define libslic3r_SeamAligned_hpp_
+
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+#include "libslic3r/GCode/SeamChoice.hpp"
+#include "libslic3r/Point.hpp"
+
+namespace Slic3r::ModelInfo {
+    struct Visibility;
+}
+
+namespace Slic3r::Seams::Aligned {
+
+using SeamChoiceVisibility  = std::function<double(const SeamChoice &, const Perimeters::Perimeter &)>;
+
+namespace Impl {
+struct SeamOptions
+{
+    std::size_t closest;
+    std::size_t adjacent;
+    bool adjacent_forward;
+    std::optional<std::size_t> snapped;
+    Vec2d on_edge;
+};
+
+SeamChoice pick_seam_option(const Perimeters::Perimeter &perimeter, const SeamOptions &options);
+
+std::optional<std::size_t> snap_to_angle(
+    const Vec2d &point,
+    const std::size_t search_start,
+    const Perimeters::Perimeter &perimeter,
+    const double max_detour
+);
+
+SeamOptions get_seam_options(
+    const Perimeters::Perimeter &perimeter,
+    const Vec2d &prefered_position,
+    const Perimeters::Perimeter::PointTree &points_tree,
+    const double max_detour
+);
+
+struct Nearest
+{
+    Vec2d prefered_position;
+    double max_detour;
+
+    std::optional<SeamChoice> operator()(
+        const Perimeters::Perimeter &perimeter,
+        const Perimeters::PointType point_type,
+        const Perimeters::PointClassification point_classification
+    ) const;
+};
+
+struct LeastVisible
+{
+    std::optional<SeamChoice> operator()(
+        const Perimeters::Perimeter &perimeter,
+        const Perimeters::PointType point_type,
+        const Perimeters::PointClassification point_classification
+    ) const;
+
+    const std::vector<double> &precalculated_visibility;
+};
+}
+
+
+struct VisibilityCalculator
+{
+    const Slic3r::ModelInfo::Visibility &points_visibility;
+    double convex_visibility_modifier;
+    double concave_visibility_modifier;
+
+    double operator()(const SeamChoice &choice, const Perimeters::Perimeter &perimeter) const;
+
+private:
+    static double get_angle_visibility_modifier(
+        const double angle,
+        const double convex_visibility_modifier,
+        const double concave_visibility_modifier
+    );
+};
+
+struct Params {
+    double max_detour{};
+    double jump_visibility_threshold{};
+};
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const SeamChoiceVisibility& visibility_calculator,
+    const Params& params
+);
+
+} // namespace Slic3r::Seams::Aligned
+
+#endif // libslic3r_SeamAligned_hpp_

src/libslic3r/GCode/SeamChoice.cpp

@@ -0,0 +1,108 @@
+#include "libslic3r/GCode/SeamChoice.hpp"
+
+namespace Slic3r::Seams {
+std::optional<SeamChoice> maybe_choose_seam_point(
+    const Perimeters::Perimeter &perimeter, const SeamPicker &seam_picker
+) {
+    using Perimeters::PointType;
+    using Perimeters::PointClassification;
+
+    std::vector<PointType>
+        type_search_order{PointType::enforcer, PointType::common, PointType::blocker};
+    std::vector<PointClassification> classification_search_order{
+        PointClassification::embedded, PointClassification::common, PointClassification::overhang};
+    for (const PointType point_type : type_search_order) {
+        for (const PointClassification point_classification : classification_search_order) {
+            if (std::optional<SeamChoice> seam_choice{
+                    seam_picker(perimeter, point_type, point_classification)}) {
+                return seam_choice;
+            }
+        }
+    }
+
+    return std::nullopt;
+}
+
+SeamChoice choose_seam_point(const Perimeters::Perimeter &perimeter, const SeamPicker &seam_picker) {
+    using Perimeters::PointType;
+    using Perimeters::PointClassification;
+
+    std::optional<SeamChoice> seam_choice{maybe_choose_seam_point(perimeter, seam_picker)};
+
+    if (seam_choice) {
+        return *seam_choice;
+    }
+
+    // Failed to choose any reasonable point!
+    return SeamChoice{0, 0, perimeter.positions.front()};
+}
+
+std::optional<SeamChoice> choose_degenerate_seam_point(const Perimeters::Perimeter &perimeter) {
+    if (!perimeter.positions.empty()) {
+        return SeamChoice{0, 0, perimeter.positions.front()};
+    }
+    return std::nullopt;
+}
+
+std::optional<std::vector<SeamChoice>> maybe_get_shell_seam(
+    const Shells::Shell<> &shell,
+    const std::function<std::optional<SeamChoice>(const Perimeters::Perimeter &, std::size_t)> &chooser
+) {
+    std::vector<SeamChoice> result;
+    result.reserve(shell.size());
+    for (std::size_t i{0}; i < shell.size(); ++i) {
+        const Shells::Slice<> &slice{shell[i]};
+        if (slice.boundary.is_degenerate) {
+            if (std::optional<SeamChoice> seam_choice{
+                    choose_degenerate_seam_point(slice.boundary)}) {
+                result.push_back(*seam_choice);
+            } else {
+                result.emplace_back();
+            }
+        } else {
+            const std::optional<SeamChoice> choice{chooser(slice.boundary, i)};
+            if (!choice) {
+                return std::nullopt;
+            }
+            result.push_back(*choice);
+        }
+    }
+    return result;
+}
+
+std::vector<SeamChoice> get_shell_seam(
+    const Shells::Shell<> &shell,
+    const std::function<SeamChoice(const Perimeters::Perimeter &, std::size_t)> &chooser
+) {
+    std::optional<std::vector<SeamChoice>> seam{maybe_get_shell_seam(
+        shell,
+        [&](const Perimeters::Perimeter &perimeter, std::size_t slice_index) {
+            return chooser(perimeter, slice_index);
+        }
+    )};
+    if (!seam) {
+        // Should be unreachable as chooser always returns a SeamChoice!
+        return std::vector<SeamChoice>(shell.size());
+    }
+    return *seam;
+}
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const std::function<std::vector<SeamChoice>(const Shells::Shell<>&)> &get_shell_seam
+) {
+    std::vector<std::vector<SeamPerimeterChoice>> layer_seams(get_layer_count(shells));
+
+    for (Shells::Shell<> &shell : shells) {
+        std::vector<SeamChoice> seam{get_shell_seam(shell)};
+
+        for (std::size_t perimeter_id{}; perimeter_id < shell.size(); ++perimeter_id) {
+            const SeamChoice &choice{seam[perimeter_id]};
+            Perimeters::Perimeter &perimeter{shell[perimeter_id].boundary};
+            layer_seams[shell[perimeter_id].layer_index].emplace_back(choice, std::move(perimeter));
+        }
+    }
+
+    return layer_seams;
+}
+} // namespace Slic3r::Seams

src/libslic3r/GCode/SeamChoice.hpp

@@ -0,0 +1,72 @@
+#ifndef libslic3r_SeamChoice_hpp_
+#define libslic3r_SeamChoice_hpp_
+
+#include "libslic3r/Polygon.hpp"
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+
+namespace Slic3r::Seams {
+
+/**
+ * When previous_index == next_index, the point is at the point.
+ * Otherwise the point is at the edge.
+ */
+struct SeamChoice
+{
+    std::size_t previous_index{};
+    std::size_t next_index{};
+    Vec2d position{Vec2d::Zero()};
+};
+
+struct SeamPerimeterChoice
+{
+    SeamPerimeterChoice(const SeamChoice &choice, Perimeters::Perimeter &&perimeter)
+        : choice(choice)
+        , perimeter(std::move(perimeter))
+        , bounding_box(Polygon{Geometry::scaled(this->perimeter.positions)}.bounding_box()) {}
+
+    SeamChoice choice;
+    Perimeters::Perimeter perimeter;
+    BoundingBox bounding_box;
+};
+
+using SeamPicker = std::function<std::optional<
+    SeamChoice>(const Perimeters::Perimeter &, const Perimeters::PointType, const Perimeters::PointClassification)>;
+
+std::optional<SeamChoice> maybe_choose_seam_point(
+    const Perimeters::Perimeter &perimeter, const SeamPicker &seam_picker
+);
+
+/**
+ * Go throught points on perimeter and choose the best seam point closest to
+ * the prefered position.
+ *
+ * Points in the perimeter can be diveded into 3x3=9 categories. An example category is
+ * enforced overhanging point. These categories are searched in particualr order.
+ * For example enforced overhang will be always choosen over common embedded point, etc.
+ *
+ * A closest point is choosen from the first non-empty category.
+ */
+SeamChoice choose_seam_point(
+    const Perimeters::Perimeter &perimeter, const SeamPicker& seam_picker
+);
+
+std::optional<SeamChoice> choose_degenerate_seam_point(const Perimeters::Perimeter &perimeter);
+
+std::optional<std::vector<SeamChoice>> maybe_get_shell_seam(
+    const Shells::Shell<> &shell,
+    const std::function<std::optional<SeamChoice>(const Perimeters::Perimeter &, std::size_t)> &chooser
+);
+
+std::vector<SeamChoice> get_shell_seam(
+    const Shells::Shell<> &shell,
+    const std::function<SeamChoice(const Perimeters::Perimeter &, std::size_t)> &chooser
+);
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const std::function<std::vector<SeamChoice>(const Shells::Shell<> &)> &get_shell_seam
+);
+
+} // namespace Slic3r::Seams
+
+#endif // libslic3r_SeamChoice_hpp_

src/libslic3r/GCode/SeamGeometry.cpp

@@ -0,0 +1,357 @@
+#include "libslic3r/GCode/SeamGeometry.hpp"
+#include "KDTreeIndirect.hpp"
+#include "Layer.hpp"
+#include <fstream>
+#include <numeric>
+#include <oneapi/tbb/blocked_range.h>
+#include <oneapi/tbb/parallel_for.h>
+
+namespace Slic3r::Seams::Geometry {
+
+namespace MappingImpl {
+
+/**
+ * @brief Return 0, 1, ..., size - 1.
+ */
+std::vector<std::size_t> range(std::size_t size) {
+    std::vector<std::size_t> result(size);
+    std::iota(result.begin(), result.end(), 0);
+    return result;
+}
+
+/**
+ * @brief A link between lists.
+ */
+struct Link
+{
+    std::size_t bucket_id;
+    double weight;
+};
+
+/**
+ * @brief Get optional values. Replace any nullopt Links with new_bucket_id and increment new_bucket_id.
+ *
+ * @param links A list of optional links.
+ * @param new_bucket_id In-out parameter incremented on each nullopt replacement.
+ */
+std::vector<std::size_t> assign_buckets(
+    const std::vector<std::optional<Link>> &links, std::size_t &new_bucket_id
+) {
+    std::vector<std::size_t> result;
+    std::transform(
+        links.begin(), links.end(), std::back_inserter(result),
+        [&](const std::optional<Link> &link) {
+            if (link) {
+                return link->bucket_id;
+            }
+            return new_bucket_id++;
+        }
+    );
+    return result;
+}
+} // namespace MappingImpl
+
+Vec2d get_normal(const Vec2d &vector) { return Vec2d{vector.y(), -vector.x()}.normalized(); }
+
+Vec2d get_polygon_normal(
+    const std::vector<Vec2d> &points, const std::size_t index, const double min_arm_length
+) {
+    std::optional<std::size_t> previous_index;
+    std::optional<std::size_t> next_index;
+
+    visit_forward(index, points.size(), [&](const std::size_t index_candidate) {
+        if (index == index_candidate) {
+            return false;
+        }
+        const double distance{(points[index_candidate] - points[index]).norm()};
+        if (distance > min_arm_length) {
+            next_index = index_candidate;
+            return true;
+        }
+        return false;
+    });
+    visit_backward(index, points.size(), [&](const std::size_t index_candidate) {
+        const double distance{(points[index_candidate] - points[index]).norm()};
+        if (distance > min_arm_length) {
+            previous_index = index_candidate;
+            return true;
+        }
+        return false;
+    });
+
+    if (previous_index && next_index) {
+        const Vec2d previous_normal{
+            Geometry::get_normal(points.at(index) - points.at(*previous_index))};
+        const Vec2d next_normal{Geometry::get_normal(points.at(*next_index) - points.at(index))};
+        return (previous_normal + next_normal).normalized();
+    }
+    return Vec2d::Zero();
+}
+
+std::pair<Vec2d, double> distance_to_segment_squared(const Linef &segment, const Vec2d &point) {
+    Vec2d segment_point;
+    const double distance{line_alg::distance_to_squared(segment, point, &segment_point)};
+    return {segment_point, distance};
+}
+
+std::pair<Mapping, std::size_t> get_mapping(
+    const std::vector<std::size_t> &list_sizes, const MappingOperator &mapping_operator
+) {
+    using namespace MappingImpl;
+
+    std::vector<std::vector<std::size_t>> result;
+    result.reserve(list_sizes.size());
+    result.push_back(range(list_sizes.front()));
+
+    std::size_t new_bucket_id{result.back().size()};
+
+    for (std::size_t layer_index{0}; layer_index < list_sizes.size() - 1; ++layer_index) {
+        // Current layer is already assigned mapping.
+
+        // Links on the next layer to the current layer.
+        std::vector<std::optional<Link>> links(list_sizes[layer_index + 1]);
+
+        for (std::size_t item_index{0}; item_index < list_sizes[layer_index]; ++item_index) {
+            const MappingOperatorResult next_item{
+                mapping_operator(layer_index, item_index)};
+            if (next_item) {
+                const auto [index, weight] = *next_item;
+                const Link link{result.back()[item_index], weight};
+                if (!links[index] || links[index]->weight < link.weight) {
+                    links[index] = link;
+                }
+            }
+        }
+        result.push_back(assign_buckets(links, new_bucket_id));
+    }
+    return {result, new_bucket_id};
+}
+
+Extrusion::Extrusion(BoundingBox bounding_box, const double width, const ExPolygon &island_boundary)
+    : bounding_box(std::move(bounding_box)), width(width), island_boundary(island_boundary) {
+
+    this->island_boundary_bounding_boxes.push_back(island_boundary.contour.bounding_box());
+
+    std::transform(
+        this->island_boundary.holes.begin(), this->island_boundary.holes.end(),
+        std::back_inserter(this->island_boundary_bounding_boxes),
+        [](const Polygon &polygon) { return polygon.bounding_box(); }
+    );
+}
+
+Geometry::Extrusions get_external_perimeters(const Slic3r::Layer &layer, const LayerSlice &slice) {
+    std::vector<Geometry::Extrusion> result;
+    for (const LayerIsland &island : slice.islands) {
+        const LayerRegion &layer_region = *layer.get_region(island.perimeters.region());
+        for (const uint32_t perimeter_id : island.perimeters) {
+            const auto collection{static_cast<const ExtrusionEntityCollection *>(
+                layer_region.perimeters().entities[perimeter_id]
+            )};
+            for (const ExtrusionEntity *entity : *collection) {
+                if (entity->role().is_external_perimeter()) {
+                    const BoundingBox bounding_box{entity->as_polyline().points};
+                    const double width{layer_region.flow(FlowRole::frExternalPerimeter).width()};
+                    result.emplace_back(bounding_box, width, island.boundary);
+                }
+            }
+        }
+    }
+    return result;
+}
+
+std::vector<Extrusions> get_extrusions(tcb::span<const Slic3r::Layer *const> object_layers) {
+    std::vector<Extrusions> result;
+    result.reserve(object_layers.size());
+
+    for (const Slic3r::Layer *object_layer : object_layers) {
+        Extrusions extrusions;
+
+        for (const LayerSlice &slice : object_layer->lslices_ex) {
+            std::vector<Extrusion> external_perimeters{
+                get_external_perimeters(*object_layer, slice)};
+            for (Geometry::Extrusion &extrusion : external_perimeters) {
+                extrusions.push_back(std::move(extrusion));
+            }
+        }
+
+        result.push_back(std::move(extrusions));
+    }
+
+    return result;
+}
+
+std::vector<Vec2d> oversample_edge(const Vec2d &from, const Vec2d &to, const double max_distance) {
+    const double total_distance{(from - to).norm()};
+    const auto points_count{static_cast<std::size_t>(std::ceil(total_distance / max_distance)) + 1};
+    if (points_count < 3) {
+        return {};
+    }
+    const double step_size{total_distance / (points_count - 1)};
+    const Vec2d step_vector{step_size * (to - from).normalized()};
+    std::vector<Vec2d> result;
+    result.reserve(points_count - 2);
+    for (std::size_t i{1}; i < points_count - 1; ++i) {
+        result.push_back(from + i * step_vector);
+    }
+    return result;
+}
+
+void visit_forward(
+    const std::size_t start_index,
+    const std::size_t loop_size,
+    const std::function<bool(std::size_t)> &visitor
+) {
+    std::size_t last_index{loop_size - 1};
+    std::size_t index{start_index};
+    for (unsigned _{0}; _ < 30; ++_) { // Infinite loop prevention
+        if (visitor(index)) {
+            return;
+        }
+        index = index == last_index ? 0 : index + 1;
+    }
+    assert(false);
+}
+
+void visit_backward(
+    const std::size_t start_index,
+    const std::size_t loop_size,
+    const std::function<bool(std::size_t)> &visitor
+) {
+    std::size_t last_index{loop_size - 1};
+    std::size_t index{start_index == 0 ? loop_size - 1 : start_index - 1};
+    for (unsigned _{0}; _ < 30; ++_) { // Infinite loop prevention
+        if (visitor(index)) {
+            return;
+        }
+        index = index == 0 ? last_index : index - 1;
+    }
+    assert(false);
+}
+
+std::vector<Vec2d> unscaled(const Points &points) {
+    std::vector<Vec2d> result;
+    result.reserve(points.size());
+    using std::transform, std::begin, std::end, std::back_inserter;
+    transform(begin(points), end(points), back_inserter(result), [](const Point &point) {
+        return unscaled(point);
+    });
+    return result;
+}
+
+std::vector<Linef> unscaled(const Lines &lines) {
+    std::vector<Linef> result;
+    result.reserve(lines.size());
+    std::transform(lines.begin(), lines.end(), std::back_inserter(result), [](const Line &line) {
+        return Linef{unscaled(line.a), unscaled(line.b)};
+    });
+    return result;
+}
+
+Points scaled(const std::vector<Vec2d> &points) {
+    Points result;
+    for (const Vec2d &point : points) {
+        result.push_back(Slic3r::scaled(point));
+    }
+    return result;
+}
+
+std::vector<double> get_embedding_distances(
+    const std::vector<Vec2d> &points, const AABBTreeLines::LinesDistancer<Linef> &perimeter_distancer
+) {
+    std::vector<double> result;
+    result.reserve(points.size());
+    using std::transform, std::back_inserter;
+    transform(points.begin(), points.end(), back_inserter(result), [&](const Vec2d &point) {
+        const double distance{perimeter_distancer.distance_from_lines<true>(point)};
+        return distance < 0 ? -distance : 0.0;
+    });
+    return result;
+}
+
+std::vector<double> get_overhangs(
+    const std::vector<Vec2d> &points,
+    const AABBTreeLines::LinesDistancer<Linef> &previous_layer_perimeter_distancer,
+    const double layer_height
+) {
+    std::vector<double> result;
+    result.reserve(points.size());
+    using std::transform, std::back_inserter;
+    transform(points.begin(), points.end(), back_inserter(result), [&](const Vec2d &point) {
+        const double distance{previous_layer_perimeter_distancer.distance_from_lines<true>(point)};
+        return distance > 0 ? M_PI / 2 - std::atan(layer_height / distance) : 0.0;
+    });
+    return result;
+}
+
+// Measured from outside, convex is positive
+std::vector<double> get_vertex_angles(const std::vector<Vec2d> &points, const double min_arm_length) {
+    std::vector<double> result;
+    result.reserve(points.size());
+
+    for (std::size_t index{0}; index < points.size(); ++index) {
+        std::optional<std::size_t> previous_index;
+        std::optional<std::size_t> next_index;
+
+        visit_forward(index, points.size(), [&](const std::size_t index_candidate) {
+            if (index == index_candidate) {
+                return false;
+            }
+            const double distance{(points[index_candidate] - points[index]).norm()};
+            if (distance > min_arm_length) {
+                next_index = index_candidate;
+                return true;
+            }
+            return false;
+        });
+        visit_backward(index, points.size(), [&](const std::size_t index_candidate) {
+            const double distance{(points[index_candidate] - points[index]).norm()};
+            if (distance > min_arm_length) {
+                previous_index = index_candidate;
+                return true;
+            }
+            return false;
+        });
+
+        if (previous_index && next_index) {
+            const Vec2d &previous_point = points[*previous_index];
+            const Vec2d &point = points[index];
+            const Vec2d &next_point = points[*next_index];
+            result.push_back(-angle((point - previous_point), (next_point - point)));
+        } else {
+            result.push_back(0.0);
+        }
+    }
+
+    return result;
+}
+
+std::pair<std::size_t, double> pick_closest_bounding_box(
+    const BoundingBox &to, const BoundingBoxes &choose_from
+) {
+    double min_distance{std::numeric_limits<double>::infinity()};
+    std::size_t choosen_index{0};
+
+    for (std::size_t i{0}; i < choose_from.size(); ++i) {
+        const BoundingBox &candidate{choose_from[i]};
+        const double bb_max_distance{unscaled(Point{candidate.max - to.max}).norm()};
+        const double bb_min_distance{unscaled(Point{candidate.min - to.min}).norm()};
+        const double distance{std::max(bb_max_distance, bb_min_distance)};
+
+        if (distance < min_distance) {
+            choosen_index = i;
+            min_distance = distance;
+        }
+    }
+    return {choosen_index, min_distance};
+}
+
+Polygon to_polygon(const ExtrusionLoop &loop) {
+    Points loop_points{};
+    for (const ExtrusionPath &path : loop.paths) {
+        for (const Point &point : path.polyline.points) {
+            loop_points.push_back(point);
+        }
+    }
+    return Polygon{loop_points};
+}
+} // namespace Slic3r::Seams::Geometry

src/libslic3r/GCode/SeamGeometry.hpp

@@ -0,0 +1,162 @@
+#ifndef libslic3r_SeamGeometry_hpp_
+#define libslic3r_SeamGeometry_hpp_
+
+#include "libslic3r/ExtrusionEntity.hpp"
+#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/AABBTreeLines.hpp"
+#include "libslic3r/Point.hpp"
+#include <oneapi/tbb/blocked_range.h>
+#include <oneapi/tbb/parallel_for.h>
+#include <vector>
+#include "tcbspan/span.hpp"
+
+namespace Slic3r {
+class Layer;
+}
+
+namespace Slic3r::Seams::Geometry {
+
+struct Extrusion
+{
+    Extrusion(BoundingBox bounding_box, const double width, const ExPolygon &island_boundary);
+
+    Extrusion(const Extrusion &) = delete;
+    Extrusion(Extrusion &&) = default;
+    Extrusion &operator=(const Extrusion &) = delete;
+    Extrusion &operator=(Extrusion &&) = delete;
+
+    BoundingBox bounding_box;
+    double width;
+    const ExPolygon &island_boundary;
+
+    // At index 0 there is the bounding box of contour. Rest are the bounding boxes of holes in order.
+    BoundingBoxes island_boundary_bounding_boxes;
+};
+
+using Extrusions = std::vector<Extrusion>;
+
+std::vector<Extrusions> get_extrusions(tcb::span<const Slic3r::Layer *const> object_layers);
+
+Vec2d get_polygon_normal(
+    const std::vector<Vec2d> &points, const std::size_t index, const double min_arm_length
+);
+
+Vec2d get_normal(const Vec2d &vector);
+
+std::pair<Vec2d, double> distance_to_segment_squared(const Linef &segment, const Vec2d &point);
+
+using Mapping = std::vector<std::vector<std::size_t>>;
+using MappingOperatorResult = std::optional<std::pair<std::size_t, double>>;
+using MappingOperator = std::function<MappingOperatorResult(std::size_t, std::size_t)>;
+
+/**
+ * @brief Indirectly map list of lists into buckets.
+ *
+ * Look for chains of items accross the lists.
+ * It may do this mapping: [[1, 2], [3, 4, 5], [6]] -> [[1, 4, 6], [2, 3], [5]].
+ * It depends on the weights provided by the mapping operator.
+ *
+ * Same bucket cannot be choosen for multiple items in any of the inner lists.
+ * Bucket is choosen **based on the weight** provided by the mapping operator. Multiple items from
+ * the same list may want to claim the same bucket. In that case, the item with the biggest weight
+ * wins the bucket. For example: [[1, 2], [3]] -> [[1, 3], [2]]
+ *
+ * @param list_sizes Vector of sizes of the original lists in a list.
+ * @param mapping_operator Operator that takes layer index and item index on that layer as input
+ * and returns the best fitting item index from the next layer, along with weight, representing how
+ * good the fit is. It may return nullopt if there is no good fit.
+ *
+ * @return Mapping [outter_list_index][inner_list_index] -> bucket id and the number of buckets.
+ */
+std::pair<Mapping, std::size_t> get_mapping(
+    const std::vector<std::size_t> &list_sizes, const MappingOperator &mapping_operator
+);
+
+std::vector<Vec2d> oversample_edge(const Vec2d &from, const Vec2d &to, const double max_distance);
+
+template <typename NestedVector>
+std::size_t get_flat_size(const NestedVector &nested_vector) {
+    return std::accumulate(
+        nested_vector.begin(), nested_vector.end(), std::size_t{0},
+        [](const std::size_t sum, const auto &vector) { return sum + vector.size(); }
+    );
+}
+
+template<typename NestedVector>
+std::vector<std::pair<std::size_t, std::size_t>> get_flat_index2indices_table(
+    const NestedVector &nested_vector
+) {
+    std::vector<std::pair<std::size_t, std::size_t>> result;
+    for (std::size_t parent_index{0}; parent_index < nested_vector.size(); ++parent_index) {
+        const auto &vector{nested_vector[parent_index]};
+        for (std::size_t nested_index{0}; nested_index < vector.size(); ++nested_index) {
+            result.push_back({parent_index, nested_index});
+        }
+    }
+    return result;
+}
+
+template <typename NestedVector>
+void iterate_nested(const NestedVector &nested_vector, const std::function<void(std::size_t, std::size_t)> &function) {
+    std::size_t flat_size{Geometry::get_flat_size(nested_vector)};
+    using Range = tbb::blocked_range<size_t>;
+    const Range range{0, flat_size};
+
+    std::vector<std::pair<std::size_t, std::size_t>> index_table{
+        get_flat_index2indices_table(nested_vector)};
+
+    // Iterate the shells as if it was flat.
+    tbb::parallel_for(range, [&](Range range) {
+        for (std::size_t index{range.begin()}; index < range.end(); ++index) {
+            const auto[parent_index, nested_index]{index_table[index]};
+            function(parent_index, nested_index);
+        }
+    });
+}
+
+void visit_forward(
+    const std::size_t start_index,
+    const std::size_t loop_size,
+    const std::function<bool(std::size_t)> &visitor
+);
+void visit_backward(
+    const std::size_t start_index,
+    const std::size_t loop_size,
+    const std::function<bool(std::size_t)> &visitor
+);
+
+std::vector<Vec2d> unscaled(const Points &points);
+
+std::vector<Linef> unscaled(const Lines &lines);
+
+Points scaled(const std::vector<Vec2d> &points);
+
+std::vector<double> get_embedding_distances(
+    const std::vector<Vec2d> &points, const AABBTreeLines::LinesDistancer<Linef> &perimeter_distancer
+);
+
+/**
+ * @brief Calculate overhang angle for each of the points over the previous layer perimeters.
+ *
+ * Larger angle <=> larger overhang. E.g. floating box has overhang = PI / 2.
+ *
+ * @returns Angles in radians <0, PI / 2>.
+ */
+std::vector<double> get_overhangs(
+    const std::vector<Vec2d> &points,
+    const AABBTreeLines::LinesDistancer<Linef> &previous_layer_perimeter_distancer,
+    const double layer_height
+);
+
+// Measured from outside, convex is positive
+std::vector<double> get_vertex_angles(const std::vector<Vec2d> &points, const double min_arm_length);
+
+std::pair<std::size_t, double> pick_closest_bounding_box(
+    const BoundingBox &to, const BoundingBoxes &choose_from
+);
+
+Polygon to_polygon(const ExtrusionLoop &loop);
+
+} // namespace Slic3r::Seams::Geometry
+
+#endif // libslic3r_SeamGeometry_hpp_

src/libslic3r/GCode/SeamPainting.cpp

@@ -0,0 +1,48 @@
+#include "libslic3r/GCode/SeamPainting.hpp"
+
+namespace Slic3r::Seams::ModelInfo {
+Painting::Painting(const Transform3d &obj_transform, const ModelVolumePtrs &volumes) {
+    for (const ModelVolume *mv : volumes) {
+        if (mv->is_seam_painted()) {
+            auto model_transformation = obj_transform * mv->get_matrix();
+
+            indexed_triangle_set enforcers = mv->seam_facets
+                                                 .get_facets(*mv, EnforcerBlockerType::ENFORCER);
+            its_transform(enforcers, model_transformation);
+            its_merge(this->enforcers, enforcers);
+
+            indexed_triangle_set blockers = mv->seam_facets
+                                                .get_facets(*mv, EnforcerBlockerType::BLOCKER);
+            its_transform(blockers, model_transformation);
+            its_merge(this->blockers, blockers);
+        }
+    }
+
+    this->enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
+        this->enforcers.vertices, this->enforcers.indices
+    );
+    this->blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
+        this->blockers.vertices, this->blockers.indices
+    );
+}
+
+bool Painting::is_enforced(const Vec3f &position, float radius) const {
+    if (enforcers.empty()) {
+        return false;
+    }
+    float radius_sqr = radius * radius;
+    return AABBTreeIndirect::is_any_triangle_in_radius(
+        enforcers.vertices, enforcers.indices, enforcers_tree, position, radius_sqr
+    );
+}
+
+bool Painting::is_blocked(const Vec3f &position, float radius) const {
+    if (blockers.empty()) {
+        return false;
+    }
+    float radius_sqr = radius * radius;
+    return AABBTreeIndirect::is_any_triangle_in_radius(
+        blockers.vertices, blockers.indices, blockers_tree, position, radius_sqr
+    );
+}
+} // namespace Slic3r::Seams::ModelInfo

src/libslic3r/GCode/SeamPainting.hpp

@@ -0,0 +1,24 @@
+#ifndef libslic3r_GlobalModelInfo_hpp_
+#define libslic3r_GlobalModelInfo_hpp_
+
+#include "libslic3r/AABBTreeIndirect.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/Model.hpp"
+
+namespace Slic3r::Seams::ModelInfo {
+class Painting
+{
+public:
+    Painting(const Transform3d &obj_transform, const ModelVolumePtrs &volumes);
+
+    bool is_enforced(const Vec3f &position, float radius) const;
+    bool is_blocked(const Vec3f &position, float radius) const;
+
+private:
+    indexed_triangle_set enforcers;
+    indexed_triangle_set blockers;
+    AABBTreeIndirect::Tree<3, float> enforcers_tree;
+    AABBTreeIndirect::Tree<3, float> blockers_tree;
+};
+} // namespace Slic3r::Seams::ModelInfo
+#endif // libslic3r_GlobalModelInfo_hpp_

src/libslic3r/GCode/SeamPerimeters.cpp

@@ -0,0 +1,403 @@
+#include <fstream>
+
+#include "ClipperUtils.hpp"
+
+#include "libslic3r/Layer.hpp"
+
+#include "libslic3r/GCode/SeamGeometry.hpp"
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+
+namespace Slic3r::Seams::Perimeters::Impl {
+
+std::vector<Vec2d> oversample_painted(
+    const std::vector<Vec2d> &points,
+    const std::function<bool(Vec3f, double)> &is_painted,
+    const double slice_z,
+    const double max_distance
+) {
+    std::vector<Vec2d> result;
+
+    for (std::size_t index{0}; index < points.size(); ++index) {
+        const Vec2d &point{points[index]};
+
+        result.push_back(point);
+
+        const std::size_t next_index{index == points.size() - 1 ? 0 : index + 1};
+        const Vec2d &next_point{points[next_index]};
+        const float next_point_distance{static_cast<float>((point - next_point).norm())};
+        const Vec2d middle_point{(point + next_point) / 2.0};
+        Vec3f point3d{to_3d(middle_point, slice_z).cast<float>()};
+        if (is_painted(point3d, next_point_distance / 2.0)) {
+            for (const Vec2d &edge_point :
+                 Geometry::oversample_edge(point, next_point, max_distance)) {
+                result.push_back(edge_point);
+            }
+        }
+    }
+    return result;
+}
+
+std::pair<std::vector<Vec2d>, std::vector<PointType>> remove_redundant_points(
+    const std::vector<Vec2d> &points,
+    const std::vector<PointType> &point_types,
+    const double tolerance
+) {
+    std::vector<Vec2d> points_result;
+    std::vector<PointType> point_types_result;
+
+    auto range_start{points.begin()};
+
+    for (auto iterator{points.begin()}; iterator != points.end(); ++iterator) {
+        const std::int64_t index{std::distance(points.begin(), iterator)};
+        if (next(iterator) == points.end() || point_types[index] != point_types[index + 1]) {
+            std::vector<Vec2d> simplification_result;
+            douglas_peucker<double>(
+                range_start, next(iterator), std::back_inserter(simplification_result), tolerance,
+                [](const Vec2d &point) { return point; }
+            );
+
+            points_result.insert(
+                points_result.end(), simplification_result.begin(), simplification_result.end()
+            );
+            const std::vector<PointType>
+                point_types_to_add(simplification_result.size(), point_types[index]);
+            point_types_result.insert(
+                point_types_result.end(), point_types_to_add.begin(), point_types_to_add.end()
+            );
+
+            range_start = next(iterator);
+        }
+    }
+
+    return {points_result, point_types_result};
+}
+
+std::vector<PointType> get_point_types(
+    const std::vector<Vec2d> &positions,
+    const ModelInfo::Painting &painting,
+    const double slice_z,
+    const double painting_radius
+) {
+    std::vector<PointType> result;
+    result.reserve(positions.size());
+    using std::transform, std::back_inserter;
+    transform(
+        positions.begin(), positions.end(), back_inserter(result),
+        [&](const Vec2d &point) {
+            const Vec3f point3d{to_3d(point.cast<float>(), static_cast<float>(slice_z))};
+            if (painting.is_blocked(point3d, painting_radius)) {
+                return PointType::blocker;
+            }
+            if (painting.is_enforced(point3d, painting_radius)) {
+                return PointType::enforcer;
+            }
+            return PointType::common;
+        }
+    );
+    return result;
+}
+
+std::vector<PointClassification> classify_points(
+    const std::vector<double> &embeddings,
+    const std::optional<std::vector<double>> &overhangs,
+    const double overhang_threshold,
+    const double embedding_threshold
+) {
+    std::vector<PointClassification> result;
+    result.reserve(embeddings.size());
+    using std::transform, std::back_inserter;
+    transform(
+        embeddings.begin(), embeddings.end(), back_inserter(result),
+        [&, i = 0](const double embedding) mutable {
+            const unsigned index = i++;
+            if (overhangs && overhangs->operator[](index) > overhang_threshold) {
+                return PointClassification::overhang;
+            }
+            if (embedding > embedding_threshold) {
+                return PointClassification::embedded;
+            }
+            return PointClassification::common;
+        }
+    );
+    return result;
+}
+
+std::vector<AngleType> get_angle_types(
+    const std::vector<double> &angles, const double convex_threshold, const double concave_threshold
+) {
+    std::vector<AngleType> result;
+    using std::transform, std::back_inserter;
+    transform(angles.begin(), angles.end(), back_inserter(result), [&](const double angle) {
+        if (angle > convex_threshold) {
+            return AngleType::convex;
+        }
+        if (angle < -concave_threshold) {
+            return AngleType::concave;
+        }
+        return AngleType::smooth;
+    });
+    return result;
+}
+
+std::vector<AngleType> merge_angle_types(
+    const std::vector<AngleType> &angle_types,
+    const std::vector<AngleType> &smooth_angle_types,
+    const std::vector<Vec2d> &points,
+    const double min_arm_length
+) {
+    std::vector<AngleType> result;
+    result.reserve(angle_types.size());
+    for (std::size_t index{0}; index < angle_types.size(); ++index) {
+        const AngleType &angle_type{angle_types[index]};
+        const AngleType &smooth_angle_type{smooth_angle_types[index]};
+
+        AngleType resulting_type{angle_type};
+
+        if (smooth_angle_type != angle_type && smooth_angle_type != AngleType::smooth) {
+            resulting_type = smooth_angle_type;
+
+            // Check if there is a sharp angle in the vicinity. If so, do not use the smooth angle.
+            Geometry::visit_forward(index, angle_types.size(), [&](const std::size_t forward_index) {
+                const double distance{(points[forward_index] - points[index]).norm()};
+                if (distance > min_arm_length) {
+                    return true;
+                }
+                if (angle_types[forward_index] == smooth_angle_type) {
+                    resulting_type = angle_type;
+                }
+                return false;
+            });
+            Geometry::visit_backward(index, angle_types.size(), [&](const std::size_t backward_index) {
+                const double distance{(points[backward_index] - points[index]).norm()};
+                if (distance > min_arm_length) {
+                    return true;
+                }
+                if (angle_types[backward_index] == smooth_angle_type) {
+                    resulting_type = angle_type;
+                }
+                return false;
+            });
+        }
+        result.push_back(resulting_type);
+    }
+    return result;
+}
+
+} // namespace Slic3r::Seams::Perimeters::Impl
+
+namespace Slic3r::Seams::Perimeters {
+
+LayerInfos get_layer_infos(
+    tcb::span<const Slic3r::Layer* const> object_layers, const double elephant_foot_compensation
+) {
+    LayerInfos result(object_layers.size());
+
+    using Range = tbb::blocked_range<size_t>;
+    const Range range{0, object_layers.size()};
+    tbb::parallel_for(range, [&](Range range) {
+        for (std::size_t layer_index{range.begin()}; layer_index < range.end(); ++layer_index) {
+            result[layer_index] = LayerInfo::create(
+                *object_layers[layer_index], layer_index, elephant_foot_compensation
+            );
+        }
+    });
+    return result;
+}
+
+LayerInfo LayerInfo::create(
+    const Slic3r::Layer &object_layer,
+    const std::size_t index,
+    const double elephant_foot_compensation
+) {
+    AABBTreeLines::LinesDistancer<Linef> perimeter_distancer{
+        to_unscaled_linesf({object_layer.lslices})};
+
+    using PreviousLayerDistancer = std::optional<AABBTreeLines::LinesDistancer<Linef>>;
+    PreviousLayerDistancer previous_layer_perimeter_distancer;
+    if (object_layer.lower_layer != nullptr) {
+        previous_layer_perimeter_distancer = PreviousLayerDistancer{
+            to_unscaled_linesf(object_layer.lower_layer->lslices)};
+    }
+
+    return {
+        std::move(perimeter_distancer),
+        std::move(previous_layer_perimeter_distancer),
+        index,
+        object_layer.height,
+        object_layer.slice_z,
+        index == 0 ? elephant_foot_compensation : 0.0};
+}
+
+double Perimeter::IndexToCoord::operator()(const size_t index, size_t dim) const {
+    return positions[index][dim];
+}
+
+Perimeter::PointTrees get_kd_trees(
+    const PointType point_type,
+    const std::vector<PointType> &all_point_types,
+    const std::vector<PointClassification> &point_classifications,
+    const Perimeter::IndexToCoord &index_to_coord
+) {
+    std::vector<std::size_t> overhang_indexes;
+    std::vector<std::size_t> embedded_indexes;
+    std::vector<std::size_t> common_indexes;
+    for (std::size_t i{0}; i < all_point_types.size(); ++i) {
+        if (all_point_types[i] == point_type) {
+            switch (point_classifications[i]) {
+            case PointClassification::overhang: overhang_indexes.push_back(i); break;
+            case PointClassification::embedded: embedded_indexes.push_back(i); break;
+            case PointClassification::common: common_indexes.push_back(i); break;
+            }
+        }
+    }
+    Perimeter::PointTrees trees;
+    if (!overhang_indexes.empty()) {
+        trees.overhanging_points = Perimeter::PointTree{index_to_coord};
+        trees.overhanging_points->build(overhang_indexes);
+    }
+    if (!embedded_indexes.empty()) {
+        trees.embedded_points = Perimeter::PointTree{index_to_coord};
+        trees.embedded_points->build(embedded_indexes);
+    }
+    if (!common_indexes.empty()) {
+        trees.common_points = Perimeter::PointTree{index_to_coord};
+        trees.common_points->build(common_indexes);
+    }
+    return trees;
+}
+
+Perimeter::Perimeter(
+    const double slice_z,
+    const std::size_t layer_index,
+    std::vector<Vec2d> &&positions,
+    std::vector<double> &&angles,
+    std::vector<PointType> &&point_types,
+    std::vector<PointClassification> &&point_classifications,
+    std::vector<AngleType> &&angle_types
+)
+    : slice_z(slice_z)
+    , layer_index(layer_index)
+    , positions(std::move(positions))
+    , angles(std::move(angles))
+    , index_to_coord(IndexToCoord{tcb::span{this->positions}})
+    , point_types(std::move(point_types))
+    , point_classifications(std::move(point_classifications))
+    , angle_types(std::move(angle_types))
+    , enforced_points(get_kd_trees(
+          PointType::enforcer, this->point_types, this->point_classifications, this->index_to_coord
+      ))
+    , common_points(get_kd_trees(
+          PointType::common, this->point_types, this->point_classifications, this->index_to_coord
+      ))
+    , blocked_points(get_kd_trees(
+          PointType::blocker, this->point_types, this->point_classifications, this->index_to_coord
+      )) {}
+
+Perimeter Perimeter::create_degenerate(
+    std::vector<Vec2d> &&points, const double slice_z, const std::size_t layer_index
+) {
+    std::vector<PointType> point_types(points.size(), PointType::common);
+    std::vector<PointClassification>
+        point_classifications(points.size(), PointClassification::common);
+    std::vector<double> angles(points.size());
+    std::vector<AngleType> angle_types(points.size(), AngleType::smooth);
+    Perimeter perimeter{
+        slice_z,
+        layer_index,
+        std::move(points),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_types)};
+    perimeter.is_degenerate = true;
+    return perimeter;
+}
+
+Perimeter Perimeter::create(
+    const Polygon &polygon,
+    const ModelInfo::Painting &painting,
+    const LayerInfo &layer_info,
+    const PerimeterParams &params
+) {
+    if (polygon.size() < 3) {
+        return Perimeter::create_degenerate(
+            Geometry::unscaled(polygon.points), layer_info.slice_z, layer_info.index
+        );
+    }
+    std::vector<Vec2d> points;
+    if (layer_info.elephant_foot_compensation > 0) {
+        const Polygons expanded{expand(polygon, scaled(layer_info.elephant_foot_compensation))};
+        if (expanded.empty()) {
+            points = Geometry::unscaled(polygon.points);
+        } else {
+            points = Geometry::unscaled(expanded.front().points);
+        }
+    } else {
+        points = Geometry::unscaled(polygon.points);
+    }
+
+    auto is_painted{[&](const Vec3f &point, const double radius) {
+        return painting.is_enforced(point, radius) || painting.is_blocked(point, radius);
+    }};
+
+    std::vector<Vec2d> perimeter_points{
+        Impl::oversample_painted(points, is_painted, layer_info.slice_z, params.oversampling_max_distance)};
+
+    std::vector<PointType> point_types{
+        Impl::get_point_types(perimeter_points, painting, layer_info.slice_z, params.painting_radius)};
+
+    std::tie(perimeter_points, point_types) =
+        Impl::remove_redundant_points(perimeter_points, point_types, params.simplification_epsilon);
+
+    const std::vector<double> embeddings{
+        Geometry::get_embedding_distances(perimeter_points, layer_info.distancer)};
+    std::optional<std::vector<double>> overhangs;
+    if (layer_info.previous_distancer) {
+        overhangs = Geometry::get_overhangs(
+            perimeter_points, *layer_info.previous_distancer, layer_info.height
+        );
+    }
+    std::vector<PointClassification> point_classifications{
+        Impl::classify_points(embeddings, overhangs, params.overhang_threshold, params.embedding_threshold)};
+
+    std::vector<double> smooth_angles{Geometry::get_vertex_angles(perimeter_points, params.smooth_angle_arm_length)};
+    std::vector<double> angles{Geometry::get_vertex_angles(perimeter_points, params.sharp_angle_arm_length)};
+    std::vector<AngleType> angle_types{
+        Impl::get_angle_types(angles, params.convex_threshold, params.concave_threshold)};
+    std::vector<AngleType> smooth_angle_types{
+        Impl::get_angle_types(smooth_angles, params.convex_threshold, params.concave_threshold)};
+    angle_types = Impl::merge_angle_types(angle_types, smooth_angle_types, perimeter_points, params.smooth_angle_arm_length);
+
+    return Perimeter{
+        layer_info.slice_z,
+        layer_info.index,
+        std::move(perimeter_points),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_types)};
+}
+
+Shells::Shells<> create_perimeters(
+    const std::vector<Shells::Shell<Polygon>> &shells,
+    const std::vector<LayerInfo> &layer_infos,
+    const ModelInfo::Painting &painting,
+    const PerimeterParams &params
+) {
+    std::vector<Shells::Shell<>> result;
+    result.reserve(shells.size());
+    std::transform(
+        shells.begin(), shells.end(), std::back_inserter(result),
+        [](const Shells::Shell<Polygon> &shell) { return Shells::Shell<>(shell.size()); }
+    );
+
+    Geometry::iterate_nested(shells, [&](const std::size_t shell_index, const std::size_t polygon_index){
+        const Shells::Shell<Polygon> &shell{shells[shell_index]};
+        const Shells::Slice<Polygon>& slice{shell[polygon_index]};
+        const Polygon &polygon{slice.boundary};
+        const LayerInfo &layer_info{layer_infos[slice.layer_index]};
+        result[shell_index][polygon_index] = {Perimeter::create(polygon, painting, layer_info, params), slice.layer_index};
+    });
+    return result;
+}
+} // namespace Slic3r::Seams::Perimeter

src/libslic3r/GCode/SeamPerimeters.hpp

@@ -0,0 +1,188 @@
+#ifndef libslic3r_SeamPerimeters_hpp_
+#define libslic3r_SeamPerimeters_hpp_
+
+#include <tcbspan/span.hpp>
+
+#include "libslic3r/GCode/SeamPainting.hpp"
+#include "libslic3r/KDTreeIndirect.hpp"
+
+#include "libslic3r/GCode/SeamShells.hpp"
+
+namespace Slic3r {
+    class Layer;
+}
+
+namespace Slic3r::Seams::ModelInfo {
+class Painting;
+}
+
+namespace Slic3r::Seams::Perimeters {
+enum class AngleType;
+enum class PointType;
+enum class PointClassification;
+struct Perimeter;
+struct PerimeterParams;
+
+struct LayerInfo
+{
+    static LayerInfo create(
+        const Slic3r::Layer &object_layer, std::size_t index, const double elephant_foot_compensation
+    );
+
+    AABBTreeLines::LinesDistancer<Linef> distancer;
+    std::optional<AABBTreeLines::LinesDistancer<Linef>> previous_distancer;
+    std::size_t index;
+    double height{};
+    double slice_z{};
+    double elephant_foot_compensation;
+};
+
+using LayerInfos = std::vector<LayerInfo>;
+
+/**
+ * @brief Construct LayerInfo for each of the provided layers.
+ */
+LayerInfos get_layer_infos(
+    tcb::span<const Slic3r::Layer* const> object_layers, const double elephant_foot_compensation
+);
+} // namespace Slic3r::Seams::Perimeters
+
+namespace Slic3r::Seams::Perimeters::Impl {
+
+
+/**
+ * @brief Split edges between points into multiple points if there is a painted point anywhere on
+ * the edge.
+ *
+ * The edge will be split by points no more than max_distance apart.
+ * Smaller max_distance -> more points.
+ *
+ * @return All the points (original and added) in order along the edges.
+ */
+std::vector<Vec2d> oversample_painted(
+    const std::vector<Vec2d> &points,
+    const std::function<bool(Vec3f, double)> &is_painted,
+    const double slice_z,
+    const double max_distance
+);
+
+/**
+ * @brief Call Duglas-Peucker for consecutive points of the same type.
+ *
+ * It never removes the first point and last point.
+ *
+ * @param tolerance Douglas-Peucker epsilon.
+ */
+std::pair<std::vector<Vec2d>, std::vector<PointType>> remove_redundant_points(
+    const std::vector<Vec2d> &points,
+    const std::vector<PointType> &point_types,
+    const double tolerance
+);
+
+} // namespace Slic3r::Seams::Perimeters::Impl
+
+namespace Slic3r::Seams::Perimeters {
+
+enum class AngleType { convex, concave, smooth };
+
+enum class PointType { enforcer, blocker, common };
+
+enum class PointClassification { overhang, embedded, common };
+
+struct PerimeterParams
+{
+    double elephant_foot_compensation{};
+    double oversampling_max_distance{};
+    double embedding_threshold{};
+    double overhang_threshold{};
+    double convex_threshold{};
+    double concave_threshold{};
+    double painting_radius{};
+    double simplification_epsilon{};
+    double smooth_angle_arm_length{};
+    double sharp_angle_arm_length{};
+};
+
+struct Perimeter
+{
+    struct IndexToCoord
+    {
+        double operator()(const size_t index, size_t dim) const;
+
+        tcb::span<const Vec2d> positions;
+    };
+
+    using PointTree = KDTreeIndirect<2, double, IndexToCoord>;
+    using OptionalPointTree = std::optional<PointTree>;
+
+    struct PointTrees
+    {
+        OptionalPointTree embedded_points;
+        OptionalPointTree common_points;
+        OptionalPointTree overhanging_points;
+    };
+
+    Perimeter() = default;
+
+    Perimeter(
+        const double slice_z,
+        const std::size_t layer_index,
+        std::vector<Vec2d> &&positions,
+        std::vector<double> &&angles,
+        std::vector<PointType> &&point_types,
+        std::vector<PointClassification> &&point_classifications,
+        std::vector<AngleType> &&angle_types
+    );
+
+    static Perimeter create(
+        const Polygon &polygon,
+        const ModelInfo::Painting &painting,
+        const LayerInfo &layer_info,
+        const PerimeterParams &params
+    );
+
+    static Perimeter create_degenerate(
+        std::vector<Vec2d> &&points, const double slice_z, const std::size_t layer_index
+    );
+
+    bool is_degenerate{false};
+    double slice_z{};
+    std::size_t layer_index{};
+    std::vector<Vec2d> positions{};
+    std::vector<double> angles{};
+    IndexToCoord index_to_coord{};
+    std::vector<PointType> point_types{};
+    std::vector<PointClassification> point_classifications{};
+    std::vector<AngleType> angle_types{};
+
+    PointTrees enforced_points{};
+    PointTrees common_points{};
+    PointTrees blocked_points{};
+};
+
+/**
+ * @brief Create a Perimeter for each polygon in each of the shells.
+ */
+Shells::Shells<Perimeter> create_perimeters(
+    const std::vector<Shells::Shell<Polygon>> &shells,
+    const std::vector<LayerInfo> &layer_infos,
+    const ModelInfo::Painting &painting,
+    const PerimeterParams &params
+);
+
+inline std::size_t get_layer_count(
+    const Shells::Shells<> &shells
+) {
+    std::size_t layer_count{0};
+    for (const Shells::Shell<> &shell : shells) {
+        for (const Shells::Slice<>& slice : shell) {
+            if (slice.layer_index >= layer_count) {
+                layer_count = slice.layer_index + 1;
+            }
+        }
+    }
+    return layer_count;
+}
+} // namespace Slic3r::Seams::Perimeters
+
+#endif // libslic3r_SeamPerimeters_hpp_

src/libslic3r/GCode/SeamPlacer.hpp

@@ -1,4 +1,5 @@
-///|/ Copyright (c) Prusa Research 2020 - 2022 Pavel Mikuš @Godrak, Lukáš Matěna @lukasmatena, Vojtěch Bubník @bubnikv
+///|/ Copyright (c) Prusa Research 2020 - 2022 Pavel Mikuš @Godrak, Lukáš Matěna @lukasmatena,
+/// Vojtěch Bubník @bubnikv
 ///|/
 ///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
 ///|/
@@ -10,159 +11,62 @@
 #include <memory>
 #include <atomic>
 
-#include "libslic3r/libslic3r.h"
-#include "libslic3r/ExtrusionEntity.hpp"
+#include "libslic3r/GCode/SeamAligned.hpp"
 #include "libslic3r/Polygon.hpp"
-#include "libslic3r/PrintConfig.hpp"
-#include "libslic3r/BoundingBox.hpp"
-#include "libslic3r/AABBTreeIndirect.hpp"
-#include "libslic3r/KDTreeIndirect.hpp"
+#include "libslic3r/Print.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+#include "libslic3r/GCode/SeamChoice.hpp"
+#include "libslic3r/GCode/ModelVisibility.hpp"
 
-namespace Slic3r {
+namespace Slic3r::Seams {
 
-class PrintObject;
-class ExtrusionLoop;
-class Print;
-class Layer;
-
-namespace EdgeGrid {
-class Grid;
-}
-
-namespace SeamPlacerImpl {
-
-
-struct GlobalModelInfo;
-struct SeamComparator;
-
-enum class EnforcedBlockedSeamPoint {
-    Blocked = 0,
-    Neutral = 1,
-    Enforced = 2,
+struct BoundedPerimeter {
+    Perimeters::Perimeter perimeter;
+    BoundingBox bounding_box;
 };
 
-// struct representing single perimeter loop
-struct Perimeter {
-    size_t start_index{};
-    size_t end_index{}; //inclusive!
-    size_t seam_index{};
-    float flow_width{};
+using ObjectSeams =
+    std::unordered_map<const PrintObject *, std::vector<std::vector<SeamPerimeterChoice>>>;
+using LayerPerimeters = std::vector<std::vector<BoundedPerimeter>>;
+using ObjectLayerPerimeters = std::unordered_map<const PrintObject *, LayerPerimeters>;
 
-    // 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
+struct Params
 {
-    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();
-    }
+    double max_nearest_detour;
+    double rear_project_threshold;
+    Aligned::Params aligned;
+    double max_distance{};
+    unsigned random_seed{};
+    double convex_visibility_modifier{};
+    double concave_visibility_modifier{};
+    Perimeters::PerimeterParams perimeter;
+    Slic3r::ModelInfo::Visibility::Params visibility;
+    SeamPosition seam_preference;
+    bool staggered_inner_seams;
 };
 
-class SeamPlacer {
+std::ostream& operator<<(std::ostream& os, const Params& params);
+
+class Placer
+{
 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;
-    static constexpr size_t fast_decimation_triangle_count_target = 16000;
-    //square of number of rays per sample point
-    static constexpr size_t sqr_rays_per_sample_point = 5;
+    static Params get_params(const DynamicPrintConfig &config);
 
-    // 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
-    static constexpr float overhang_angle_threshold = 50.0f * float(PI) / 180.0f;
+    void init(
+        SpanOfConstPtrs<PrintObject> objects,
+        const Params &params,
+        const std::function<void(void)> &throw_if_canceled
+    );
 
-    // 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);
-
-    Point place_seam(const Layer *layer, const ExtrusionLoop &loop, bool external_first, const Point &last_pos) const;
+    Point place_seam(const Layer *layer, const ExtrusionLoop &loop, const Point &last_pos) const;
 
 private:
-    void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info);
-    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;
+    Params params;
+    ObjectSeams seams_per_object;
+    ObjectLayerPerimeters perimeters_per_layer;
 };
 
-} // namespace Slic3r
+} // namespace Slic3r::Seams
 
 #endif // libslic3r_SeamPlacer_hpp_

src/libslic3r/GCode/SeamRandom.cpp

@@ -0,0 +1,136 @@
+#include <random>
+
+#include "libslic3r/GCode/SeamRandom.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
+
+namespace Slic3r::Seams::Random {
+using Perimeters::PointType;
+using Perimeters::PointClassification;
+
+namespace Impl {
+std::vector<PerimeterSegment> get_segments(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) {
+    const std::vector<Vec2d> &positions{perimeter.positions};
+    const std::vector<PointType> &point_types{perimeter.point_types};
+    const std::vector<PointClassification> &point_classifications{perimeter.point_classifications};
+
+    std::optional<double> current_begin;
+    std::optional<std::size_t> current_begin_index;
+    Vec2d previous_position{positions.front()};
+    double distance{0.0};
+    std::vector<PerimeterSegment> result;
+    for (std::size_t index{0}; index < positions.size(); ++index) {
+        distance += (positions[index] - previous_position).norm();
+        previous_position = positions[index];
+
+        if (point_types[index] == point_type &&
+            point_classifications[index] == point_classification) {
+            if (!current_begin) {
+                current_begin = distance;
+                current_begin_index = index;
+            }
+        } else {
+            if (current_begin) {
+                result.push_back(PerimeterSegment{*current_begin, distance, *current_begin_index});
+            }
+            current_begin = std::nullopt;
+            current_begin_index = std::nullopt;
+        }
+    }
+
+    if (current_begin) {
+        result.push_back(PerimeterSegment{*current_begin, distance, *current_begin_index});
+    }
+    return result;
+}
+
+PerimeterSegment pick_random_segment(
+    const std::vector<PerimeterSegment> &segments, std::mt19937 &random_engine
+) {
+    double length{0.0};
+    for (const PerimeterSegment &segment : segments) {
+        length += segment.length();
+    }
+
+    std::uniform_real_distribution<double> distribution{0.0, length};
+    double random_distance{distribution(random_engine)};
+
+    double distance{0.0};
+    return *std::find_if(segments.begin(), segments.end(), [&](const PerimeterSegment &segment) {
+        if (random_distance >= distance && random_distance <= distance + segment.length()) {
+            return true;
+        }
+        distance += segment.length();
+        return false;
+    });
+}
+
+SeamChoice pick_random_point(
+    const PerimeterSegment &segment, const Perimeters::Perimeter &perimeter, std::mt19937 &random_engine
+) {
+    const std::vector<Vec2d> &positions{perimeter.positions};
+
+    if (segment.length() < std::numeric_limits<double>::epsilon()) {
+        return {segment.begin_index, segment.begin_index, positions[segment.begin_index]};
+    }
+
+    std::uniform_real_distribution<double> distribution{0.0, segment.length()};
+    const double random_distance{distribution(random_engine)};
+
+    double distance{0.0};
+    std::size_t previous_index{segment.begin_index};
+    for (std::size_t index{segment.begin_index + 1}; index < perimeter.positions.size(); ++index) {
+        const Vec2d edge{positions[index] - positions[previous_index]};
+
+        if (distance + edge.norm() >= random_distance) {
+            if (random_distance - distance < std::numeric_limits<double>::epsilon()) {
+                index = previous_index;
+            } else if (distance + edge.norm() - random_distance < std::numeric_limits<double>::epsilon()) {
+                previous_index = index;
+            }
+
+            const double remaining_distance{random_distance - distance};
+            const Vec2d position{positions[previous_index] + remaining_distance * edge.normalized()};
+            return {previous_index, index, position};
+        }
+
+        distance += edge.norm();
+        previous_index = index;
+    }
+
+    // Should be unreachable.
+    return {segment.begin_index, segment.begin_index, positions[segment.begin_index]};
+}
+
+std::optional<SeamChoice> Random::operator()(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) const {
+    std::vector<PerimeterSegment> segments{
+        get_segments(perimeter, point_type, point_classification)};
+
+    if (!segments.empty()) {
+        const PerimeterSegment segment{pick_random_segment(segments, random_engine)};
+        return pick_random_point(segment, perimeter, random_engine);
+    }
+    return std::nullopt;
+}
+} // namespace Impl
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells, const unsigned fixed_seed
+) {
+    std::mt19937 random_engine{fixed_seed};
+    const Impl::Random random{random_engine};
+
+    return Seams::get_object_seams(std::move(shells), [&](const Shells::Shell<> &shell) {
+        return Seams::get_shell_seam(shell, [&](const Perimeters::Perimeter &perimeter, std::size_t) {
+            return Seams::choose_seam_point(perimeter, random);
+        });
+    });
+}
+} // namespace Slic3r::Seams::Random

src/libslic3r/GCode/SeamRandom.hpp

@@ -0,0 +1,29 @@
+#include "libslic3r/GCode/SeamChoice.hpp"
+#include <random>
+
+namespace Slic3r::Seams::Random {
+namespace Impl {
+struct PerimeterSegment
+{
+    double begin{};
+    double end{};
+    std::size_t begin_index{};
+
+    double length() const { return end - begin; }
+};
+
+struct Random
+{
+    std::mt19937 &random_engine;
+
+    std::optional<SeamChoice> operator()(
+        const Perimeters::Perimeter &perimeter,
+        const Perimeters::PointType point_type,
+        const Perimeters::PointClassification point_classification
+    ) const;
+};
+}
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells, const unsigned fixed_seed
+);
+}

src/libslic3r/GCode/SeamRear.cpp

@@ -0,0 +1,125 @@
+#include "libslic3r/GCode/SeamRear.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
+
+namespace Slic3r::Seams::Rear {
+using Perimeters::PointType;
+using Perimeters::PointClassification;
+
+namespace Impl {
+
+BoundingBoxf get_bounding_box(const Shells::Shell<> &shell) {
+    BoundingBoxf result;
+    for (const Shells::Slice<> &slice : shell) {
+        result.merge(BoundingBoxf{slice.boundary.positions});
+    }
+    return result;
+}
+
+std::optional<SeamChoice> get_rearest_point(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) {
+    double max_y{-std::numeric_limits<double>::infinity()};
+    std::optional<std::size_t> choosen_index;
+    for (std::size_t i{0}; i < perimeter.positions.size(); ++i) {
+        const Perimeters::PointType _point_type{perimeter.point_types[i]};
+        const Perimeters::PointClassification _point_classification{perimeter.point_classifications[i]};
+
+        if (point_type == _point_type && point_classification == _point_classification) {
+            const Vec2d &position{perimeter.positions[i]};
+            if (position.y() > max_y) {
+                max_y = position.y();
+                choosen_index = i;
+            }
+        }
+    }
+    if (choosen_index) {
+        return SeamChoice{*choosen_index, *choosen_index, perimeter.positions[*choosen_index]};
+    }
+    return std::nullopt;
+}
+
+std::optional<SeamChoice> StraightLine::operator()(
+    const Perimeters::Perimeter &perimeter,
+    const PointType point_type,
+    const PointClassification point_classification
+) const {
+    std::vector<PerimeterLine> possible_lines;
+    for (std::size_t i{0}; i < perimeter.positions.size() - 1; ++i) {
+        if (perimeter.point_types[i] != point_type) {
+            continue;
+        }
+        if (perimeter.point_classifications[i] != point_classification) {
+            continue;
+        }
+        if (perimeter.point_types[i + 1] != point_type) {
+            continue;
+        }
+        if (perimeter.point_classifications[i + 1] != point_classification) {
+            continue;
+        }
+        possible_lines.push_back(PerimeterLine{perimeter.positions[i], perimeter.positions[i+1], i, i + 1});
+    }
+    if (possible_lines.empty()) {
+        return std::nullopt;
+    }
+
+    const AABBTreeLines::LinesDistancer<PerimeterLine> possible_distancer{possible_lines};
+    const BoundingBoxf bounding_box{perimeter.positions};
+
+    const std::vector<std::pair<Vec2d, std::size_t>> intersections{
+        possible_distancer.intersections_with_line<true>(PerimeterLine{
+            this->prefered_position, Vec2d{this->prefered_position.x(), bounding_box.min.y()},
+            0, 0})};
+    if (!intersections.empty()) {
+        const auto[position, line_index]{intersections.front()};
+        if (position.y() < bounding_box.max.y() -
+                this->rear_project_threshold * (bounding_box.max.y() - bounding_box.min.y())) {
+            return std::nullopt;
+        }
+        const PerimeterLine &intersected_line{possible_lines[line_index]};
+        const SeamChoice intersected_choice{intersected_line.previous_index, intersected_line.next_index, position};
+        return intersected_choice;
+    }
+    return std::nullopt;
+}
+} // namespace Impl
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const double rear_project_threshold
+) {
+    double average_x_center{0.0};
+    std::size_t count{0};
+    for (const Shells::Shell<> &shell : shells) {
+        for (const Shells::Slice<> &slice : shell) {
+            if (slice.boundary.positions.empty()) {
+                continue;
+            }
+            const BoundingBoxf slice_bounding_box{slice.boundary.positions};
+            average_x_center += (slice_bounding_box.min.x() + slice_bounding_box.max.x()) / 2.0;
+            count++;
+        }
+    }
+    average_x_center /= count;
+    return Seams::get_object_seams(std::move(shells), [&](const Shells::Shell<> &shell) {
+        BoundingBoxf bounding_box{Impl::get_bounding_box(shell)};
+        const Vec2d back_center{average_x_center, bounding_box.max.y()};
+        std::optional<std::vector<SeamChoice>> straight_seam {
+            Seams::maybe_get_shell_seam(shell, [&](const Perimeters::Perimeter &perimeter, std::size_t) {
+                return Seams::maybe_choose_seam_point(
+                    perimeter,
+                    Impl::StraightLine{back_center, rear_project_threshold}
+                );
+            })
+        };
+        if (!straight_seam) {
+            return Seams::get_shell_seam(shell, [&](const Perimeters::Perimeter &perimeter, std::size_t) {
+                return Seams::choose_seam_point(perimeter, Impl::get_rearest_point);
+            });
+        }
+        return *straight_seam;
+    });
+}
+} // namespace Slic3r::Seams::Rear

src/libslic3r/GCode/SeamRear.hpp

@@ -0,0 +1,40 @@
+#ifndef libslic3r_SeamRear_hpp_
+#define libslic3r_SeamRear_hpp_
+
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+#include "libslic3r/GCode/SeamChoice.hpp"
+
+namespace Slic3r::Seams::Rear {
+namespace Impl {
+struct PerimeterLine
+{
+    Vec2d a;
+    Vec2d b;
+    std::size_t previous_index;
+    std::size_t next_index;
+
+    using Scalar = Vec2d::Scalar;
+    static const constexpr int Dim = 2;
+};
+
+struct StraightLine
+{
+    Vec2d prefered_position;
+    double rear_project_threshold;
+
+    std::optional<SeamChoice> operator()(
+        const Perimeters::Perimeter &perimeter,
+        const Perimeters::PointType point_type,
+        const Perimeters::PointClassification point_classification
+    ) const;
+};
+
+} // namespace Impl
+
+std::vector<std::vector<SeamPerimeterChoice>> get_object_seams(
+    Shells::Shells<> &&shells,
+    const double rear_project_threshold
+);
+} // namespace Slic3r::Seams::Rear
+
+#endif // libslic3r_SeamRear_hpp_

src/libslic3r/GCode/SeamShells.cpp

@@ -0,0 +1,99 @@
+#include "libslic3r/GCode/SeamShells.hpp"
+#include "libslic3r/ClipperUtils.hpp"
+#include <oneapi/tbb/blocked_range.h>
+#include <oneapi/tbb/parallel_for.h>
+
+namespace Slic3r::Seams::Shells::Impl {
+
+BoundedPolygons project_to_geometry(const Geometry::Extrusions &external_perimeters) {
+    BoundedPolygons result;
+    result.reserve(external_perimeters.size());
+
+    using std::transform, std::back_inserter;
+
+    transform(
+        external_perimeters.begin(), external_perimeters.end(), back_inserter(result),
+        [](const Geometry::Extrusion &external_perimeter) {
+            const auto [choosen_index, _]{Geometry::pick_closest_bounding_box(
+                external_perimeter.bounding_box,
+                external_perimeter.island_boundary_bounding_boxes
+            )};
+
+            const Polygon &adjacent_boundary{
+                choosen_index == 0 ? external_perimeter.island_boundary.contour :
+                                     external_perimeter.island_boundary.holes[choosen_index - 1]};
+            return BoundedPolygon{adjacent_boundary, external_perimeter.island_boundary_bounding_boxes[choosen_index]};
+        }
+    );
+    return result;
+}
+
+std::vector<BoundedPolygons> project_to_geometry(const std::vector<Geometry::Extrusions> &extrusions) {
+    std::vector<BoundedPolygons> result(extrusions.size());
+
+    for (std::size_t layer_index{0}; layer_index < extrusions.size(); ++layer_index) {
+        result[layer_index] = project_to_geometry(extrusions[layer_index]);
+    }
+
+    return result;
+}
+
+Shells<Polygon> map_to_shells(
+    std::vector<BoundedPolygons> &&layers, const Geometry::Mapping &mapping, const std::size_t shell_count
+) {
+    Shells<Polygon> result(shell_count);
+    for (std::size_t layer_index{0}; layer_index < layers.size(); ++layer_index) {
+        BoundedPolygons &perimeters{layers[layer_index]};
+        for (std::size_t perimeter_index{0}; perimeter_index < perimeters.size();
+             perimeter_index++) {
+            Polygon &perimeter{perimeters[perimeter_index].polygon};
+            result[mapping[layer_index][perimeter_index]].push_back(
+                Slice<Polygon>{std::move(perimeter), layer_index}
+            );
+        }
+    }
+    return result;
+}
+} // namespace Slic3r::Seams::Shells::Impl
+
+namespace Slic3r::Seams::Shells {
+Shells<Polygon> create_shells(
+    const std::vector<Geometry::Extrusions> &extrusions, const double max_distance
+) {
+    std::vector<Impl::BoundedPolygons> projected{Impl::project_to_geometry(extrusions)};
+
+    std::vector<std::size_t> layer_sizes;
+    layer_sizes.reserve(projected.size());
+    for (const Impl::BoundedPolygons &perimeters : projected) {
+        layer_sizes.push_back(perimeters.size());
+    }
+
+    const auto &[shell_mapping, shell_count]{Geometry::get_mapping(
+        layer_sizes,
+        [&](const std::size_t layer_index,
+            const std::size_t item_index) -> Geometry::MappingOperatorResult {
+            const Impl::BoundedPolygons &layer{projected[layer_index]};
+            const Impl::BoundedPolygons &next_layer{projected[layer_index + 1]};
+            if (next_layer.empty()) {
+                return std::nullopt;
+            }
+
+            BoundingBoxes next_layer_bounding_boxes;
+            for (const Impl::BoundedPolygon &bounded_polygon : next_layer) {
+                next_layer_bounding_boxes.emplace_back(bounded_polygon.bounding_box);
+            }
+
+            const auto [perimeter_index, distance] = Geometry::pick_closest_bounding_box(
+                layer[item_index].bounding_box, next_layer_bounding_boxes
+            );
+
+            if (distance > max_distance) {
+                return std::nullopt;
+            }
+            return std::pair{perimeter_index, 1.0 / distance};
+        }
+    )};
+
+    return Impl::map_to_shells(std::move(projected), shell_mapping, shell_count);
+}
+} // namespace Slic3r::Seams::Shells

src/libslic3r/GCode/SeamShells.hpp

@@ -0,0 +1,60 @@
+#ifndef libslic3r_SeamShells_hpp_
+#define libslic3r_SeamShells_hpp_
+
+#include <vector>
+#include <tcbspan/span.hpp>
+
+#include "libslic3r/Polygon.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
+
+namespace Slic3r {
+class Layer;
+}
+
+namespace Slic3r::Seams::Perimeters {
+struct Perimeter;
+}
+
+namespace Slic3r::Seams::Shells::Impl {
+
+struct BoundedPolygon {
+    Polygon polygon;
+    BoundingBox bounding_box;
+};
+
+using BoundedPolygons = std::vector<BoundedPolygon>;
+
+/**
+ * Project extrusion path to the original mesh.
+ *
+ * Takes the extrusion path and finds the closest polygon to it in
+ * the extruison island boundary.
+ *
+ * Then it expands the extrusion path so it roughly tracks the island boundary
+ * and check that all points in expanded extrusion path are within a reasonable
+ * distance (extrusion width) from the closes polygon.
+ *
+ * If the expanded extrusion path matches the boundary it returns the
+ * closeset polygon from the island boundary. Otherwise it returns
+ * the expanded extrusion.
+ */
+BoundedPolygons project_to_geometry(const Geometry::Extrusions &extrusions);
+}
+
+namespace Slic3r::Seams::Shells {
+template<typename T = Perimeters::Perimeter> struct Slice
+{
+    T boundary;
+    std::size_t layer_index;
+};
+
+template<typename T = Perimeters::Perimeter> using Shell = std::vector<Slice<T>>;
+
+template<typename T = Perimeters::Perimeter> using Shells = std::vector<Shell<T>>;
+
+Shells<Polygon> create_shells(
+    const std::vector<Geometry::Extrusions> &extrusions, const double max_distance
+);
+} // namespace Slic3r::Seams::Shells
+
+#endif // libslic3r_SeamShells_hpp_

src/libslic3r/Layer.cpp

@@ -944,9 +944,11 @@ void Layer::sort_perimeters_into_islands(
 
     auto insert_into_island = [
         // Region where the perimeters, gap fills and fill expolygons are stored.
-        region_id, 
+        region_id,
         // Whether there are infills with different regions generated for this LayerSlice.
         has_multiple_regions,
+        // Layer split into surfaces
+        &slices,
         // Perimeters and gap fills to be sorted into islands.
         &perimeter_and_gapfill_ranges,
         // Infill regions to be sorted into islands.
@@ -959,6 +961,7 @@ void Layer::sort_perimeters_into_islands(
         lslices_ex[lslice_idx].islands.push_back({});
         LayerIsland &island = lslices_ex[lslice_idx].islands.back();
         island.perimeters = LayerExtrusionRange(region_id, perimeter_and_gapfill_ranges[source_slice_idx].first);
+        island.boundary = slices.surfaces[source_slice_idx].expolygon;
         island.thin_fills = perimeter_and_gapfill_ranges[source_slice_idx].second;
         if (ExPolygonRange fill_range = fill_expolygons_ranges[source_slice_idx]; ! fill_range.empty()) {
             if (has_multiple_regions) {

src/libslic3r/Layer.hpp

@@ -77,6 +77,8 @@ private:
     static constexpr const uint32_t fill_region_composite_id = std::numeric_limits<uint32_t>::max();
 
 public:
+    // Boundary of the LayerIsland before perimeter generation.
+    ExPolygon             boundary;
     // Perimeter extrusions in LayerRegion belonging to this island.
     LayerExtrusionRange   perimeters;
     // Thin fills of the same region as perimeters. Generated by classic perimeter generator, while Arachne puts them into perimeters.

src/libslic3r/ShortEdgeCollapse.cpp

@@ -9,11 +9,29 @@
 #include <unordered_set>
 #include <random>
 #include <algorithm>
+#include <boost/random/uniform_int_distribution.hpp>
 
 #include <ankerl/unordered_dense.h>
 
 namespace Slic3r {
 
+/**
+ * Simple implementation of Fisher-Yates algorithm using uniform int
+ * distribution from boost, ensurinng the result is the same
+ * accross platforms.
+ *
+ * DO NOT EXPECT IT TO BE PERFORMANT! Use it only when std::shuffle is
+ * not applicable.
+ */
+template<typename Range, typename UniformRandomNumberGenerator>
+void stable_shuffle(Range &range, UniformRandomNumberGenerator &generator) {
+    const int n{static_cast<int>(range.size())};
+    for (int i{0}; i < n - 2; ++i) {
+        int j{boost::random::uniform_int_distribution<int>{i, n-1}(generator)};
+        std::swap(range[i], range[j]);
+    }
+}
+
 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());
@@ -102,8 +120,8 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_
         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);
-        
+        stable_shuffle(face_indices, generator);
+
         int allowed_face_removals = int(face_indices.size()) - int(target_triangle_count);
         for (const size_t &face_idx : face_indices) {
             if (face_removal_flags[face_idx]) {

src/libslic3r/TriangleSetSampling.cpp

@@ -7,6 +7,7 @@
 #include <random>
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
+#include <boost/random/uniform_real_distribution.hpp>
 
 namespace Slic3r {
 
@@ -33,8 +34,9 @@ TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const index
     }
 
     std::mt19937_64 mersenne_engine { 27644437 };
+    // Use boost instead of std to ensure stability accross platforms!
     // random numbers on interval [0, 1)
-    std::uniform_real_distribution<double> fdistribution;
+    boost::random::uniform_real_distribution<double> fdistribution;
 
     auto get_random = [&fdistribution, &mersenne_engine]() {
         return Vec3d { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) };

tests/fff_print/CMakeLists.txt

@@ -14,6 +14,13 @@ add_executable(${_TEST_NAME}_tests
 	test_gaps.cpp
 	test_gcode.cpp
 	test_gcode_travels.cpp
+    test_seam_perimeters.cpp
+    test_seam_shells.cpp
+    test_seam_geometry.cpp
+    test_seam_aligned.cpp
+    test_seam_rear.cpp
+    test_seam_random.cpp
+    benchmark_seams.cpp
 	test_gcodefindreplace.cpp
 	test_gcodewriter.cpp
 	test_cancel_object.cpp
@@ -33,6 +40,7 @@ add_executable(${_TEST_NAME}_tests
 	)
 target_link_libraries(${_TEST_NAME}_tests test_common libslic3r)
 set_property(TARGET ${_TEST_NAME}_tests PROPERTY FOLDER "tests")
+target_compile_definitions(${_TEST_NAME}_tests PUBLIC CATCH_CONFIG_ENABLE_BENCHMARKING)
 
 if (WIN32)
     prusaslicer_copy_dlls(${_TEST_NAME}_tests)

tests/fff_print/benchmark_seams.cpp

@@ -0,0 +1,124 @@
+#include <catch2/catch.hpp>
+#include "test_data.hpp"
+
+#include "libslic3r/GCode/SeamGeometry.hpp"
+#include "libslic3r/GCode/SeamAligned.hpp"
+#include "libslic3r/GCode/SeamRear.hpp"
+#include "libslic3r/GCode/SeamRandom.hpp"
+
+TEST_CASE_METHOD(Slic3r::Test::SeamsFixture, "Seam benchmarks", "[Seams][.Benchmarks]") {
+    BENCHMARK_ADVANCED("Create extrusions benchy")(Catch::Benchmark::Chronometer meter) {
+        meter.measure([&] { return Slic3r::Seams::Geometry::get_extrusions(print_object->layers()); });
+    };
+
+    using namespace Slic3r::Seams;
+    BENCHMARK_ADVANCED("Create shells benchy")(Catch::Benchmark::Chronometer meter) {
+        meter.measure([&] { return Shells::create_shells(extrusions, params.max_distance); });
+    };
+
+
+    BENCHMARK_ADVANCED("Get layer infos benchy")(Catch::Benchmark::Chronometer meter) {
+        meter.measure([&] {
+            return Perimeters::get_layer_infos(
+                print_object->layers(), params.perimeter.elephant_foot_compensation
+            );
+        });
+    };
+
+    BENCHMARK_ADVANCED("Create perimeters benchy")(Catch::Benchmark::Chronometer meter) {
+        meter.measure([&] {
+            return Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter);
+        });
+    };
+
+    BENCHMARK_ADVANCED("Generate aligned seam benchy")(Catch::Benchmark::Chronometer meter) {
+        std::vector<Shells::Shells<>> inputs;
+        inputs.reserve(meter.runs());
+        std::generate_n(std::back_inserter(inputs), meter.runs(), [&]() {
+            return Perimeters::create_perimeters(
+                shell_polygons, layer_infos, painting, params.perimeter
+            );
+        });
+        meter.measure([&](const int i) {
+            return Aligned::get_object_seams(
+                std::move(inputs[i]), visibility_calculator, params.aligned
+            );
+        });
+    };
+
+    BENCHMARK_ADVANCED("Visibility constructor")(Catch::Benchmark::Chronometer meter) {
+        using Visibility = Slic3r::ModelInfo::Visibility;
+        std::vector<Catch::Benchmark::storage_for<Visibility>> storage(meter.runs());
+        meter.measure([&](const int i) {
+            storage[i].construct(transformation, volumes, params.visibility, []() {});
+        });
+    };
+
+    BENCHMARK_ADVANCED("Generate rear seam benchy")(Catch::Benchmark::Chronometer meter) {
+        std::vector<Shells::Shells<>> inputs;
+        inputs.reserve(meter.runs());
+        std::generate_n(std::back_inserter(inputs), meter.runs(), [&]() {
+            return create_perimeters(
+                shell_polygons, layer_infos, painting, params.perimeter
+            );
+        });
+        meter.measure([&](const int i) {
+            return Rear::get_object_seams(std::move(inputs[i]), params.rear_project_threshold);
+        });
+    };
+
+    BENCHMARK_ADVANCED("Generate random seam benchy")(Catch::Benchmark::Chronometer meter) {
+        std::vector<Shells::Shells<>> inputs;
+        inputs.reserve(meter.runs());
+        std::generate_n(std::back_inserter(inputs), meter.runs(), [&]() {
+            return Perimeters::create_perimeters(
+                shell_polygons, layer_infos, painting, params.perimeter
+            );
+        });
+        meter.measure([&](const int i) {
+            return Random::get_object_seams(std::move(inputs[i]), params.random_seed);
+        });
+    };
+
+    Placer placer;
+    BENCHMARK_ADVANCED("Init seam placer aligned")(Catch::Benchmark::Chronometer meter) {
+        meter.measure([&] {
+            return placer.init(print->objects(), params, [](){});
+        });
+    };
+
+    SECTION("Place seam"){
+        using namespace Slic3r;
+        Placer placer;
+        placer.init(print->objects(), params, [](){});
+        std::vector<std::pair<const Layer*, const ExtrusionLoop*>> loops;
+
+        const PrintObject* object{print->objects().front()};
+        for (const Layer* layer :object->layers()) {
+            for (const LayerSlice& lslice : layer->lslices_ex) {
+                for (const LayerIsland &island : lslice.islands) {
+                    const LayerRegion &layer_region = *layer->get_region(island.perimeters.region());
+                    for (uint32_t perimeter_id : island.perimeters) {
+                        const auto *entity_collection{static_cast<const ExtrusionEntityCollection*>(layer_region.perimeters().entities[perimeter_id])};
+                        if (entity_collection != nullptr) {
+                            for (const ExtrusionEntity *entity : *entity_collection) {
+                                const auto loop{static_cast<const ExtrusionLoop*>(entity)};
+                                if (loop == nullptr) {
+                                    continue;
+                                }
+                                loops.emplace_back(layer, loop);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        BENCHMARK_ADVANCED("Place seam benchy")(Catch::Benchmark::Chronometer meter) {
+            meter.measure([&] {
+                for (const auto &[layer, loop] : loops) {
+                    placer.place_seam(layer, *loop, {0, 0});
+                }
+            });
+        };
+    }
+}

tests/fff_print/test_data.hpp

@@ -2,12 +2,19 @@
 #define SLIC3R_TEST_DATA_HPP
 
 #include "libslic3r/Config.hpp"
+#include "libslic3r/Format/3mf.hpp"
+#include "libslic3r/GCode/ModelVisibility.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
+#include "libslic3r/GCode/SeamPerimeters.hpp"
 #include "libslic3r/Geometry.hpp"
 #include "libslic3r/Model.hpp"
 #include "libslic3r/Point.hpp"
 #include "libslic3r/Print.hpp"
 #include "libslic3r/TriangleMesh.hpp"
+#include "libslic3r/GCode/SeamPlacer.hpp"
+#include "libslic3r/GCode/SeamAligned.hpp"
 
+#include <filesystem>
 #include <unordered_map>
 
 namespace Slic3r { namespace Test {
@@ -39,14 +46,13 @@ enum class TestMesh {
 };
 
 // Neccessary for <c++17
-struct TestMeshHash {
-    std::size_t operator()(TestMesh tm) const {
-        return static_cast<std::size_t>(tm);
-    }
+struct TestMeshHash
+{
+    std::size_t operator()(TestMesh tm) const { return static_cast<std::size_t>(tm); }
 };
 
 /// Mesh enumeration to name mapping
-extern const std::unordered_map<TestMesh, const char*, TestMeshHash> mesh_names;
+extern const std::unordered_map<TestMesh, const char *, TestMeshHash> mesh_names;
 
 /// Port of Slic3r::Test::mesh
 /// Basic cubes/boxes should call TriangleMesh::make_cube() directly and rescale/translate it
@@ -56,35 +62,173 @@ TriangleMesh mesh(TestMesh m, Vec3d translate, Vec3d scale = Vec3d(1.0, 1.0, 1.0
 TriangleMesh mesh(TestMesh m, Vec3d translate, double scale = 1.0);
 
 /// Templated function to see if two values are equivalent (+/- epsilon)
-template <typename T>
-bool _equiv(const T& a, const T& b) { return std::abs(a - b) < EPSILON; }
+template<typename T> bool _equiv(const T &a, const T &b) { return std::abs(a - b) < EPSILON; }
 
-template <typename T>
-bool _equiv(const T& a, const T& b, double epsilon) { return abs(a - b) < epsilon; }
+template<typename T> bool _equiv(const T &a, const T &b, double epsilon) {
+    return abs(a - b) < epsilon;
+}
 
-Slic3r::Model model(const std::string& model_name, TriangleMesh&& _mesh);
-void init_print(std::vector<TriangleMesh> 		  &&meshes, Slic3r::Print &print, Slic3r::Model& model, const DynamicPrintConfig &config_in, bool comments = false, unsigned duplicate_count = 1);
-void init_print(std::initializer_list<TestMesh> 	meshes, Slic3r::Print &print, Slic3r::Model& model, const Slic3r::DynamicPrintConfig &config_in = Slic3r::DynamicPrintConfig::full_print_config(), bool comments = false, unsigned duplicate_count = 1);
-void init_print(std::initializer_list<TriangleMesh> meshes, Slic3r::Print &print, Slic3r::Model& model, const Slic3r::DynamicPrintConfig &config_in = Slic3r::DynamicPrintConfig::full_print_config(), bool comments = false, unsigned duplicate = 1);
-void init_print(std::initializer_list<TestMesh> 	meshes, Slic3r::Print &print, Slic3r::Model& model, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false, unsigned duplicate = 1);
-void init_print(std::initializer_list<TriangleMesh> meshes, Slic3r::Print &print, Slic3r::Model& model, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false, unsigned duplicate = 1);
+Slic3r::Model model(const std::string &model_name, TriangleMesh &&_mesh);
+void init_print(
+    std::vector<TriangleMesh> &&meshes,
+    Slic3r::Print &print,
+    Slic3r::Model &model,
+    const DynamicPrintConfig &config_in,
+    bool comments = false,
+    unsigned duplicate_count = 1
+);
+void init_print(
+    std::initializer_list<TestMesh> meshes,
+    Slic3r::Print &print,
+    Slic3r::Model &model,
+    const Slic3r::DynamicPrintConfig &config_in = Slic3r::DynamicPrintConfig::full_print_config(),
+    bool comments = false,
+    unsigned duplicate_count = 1
+);
+void init_print(
+    std::initializer_list<TriangleMesh> meshes,
+    Slic3r::Print &print,
+    Slic3r::Model &model,
+    const Slic3r::DynamicPrintConfig &config_in = Slic3r::DynamicPrintConfig::full_print_config(),
+    bool comments = false,
+    unsigned duplicate = 1
+);
+void init_print(
+    std::initializer_list<TestMesh> meshes,
+    Slic3r::Print &print,
+    Slic3r::Model &model,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false,
+    unsigned duplicate = 1
+);
+void init_print(
+    std::initializer_list<TriangleMesh> meshes,
+    Slic3r::Print &print,
+    Slic3r::Model &model,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false,
+    unsigned duplicate = 1
+);
 
-void init_and_process_print(std::initializer_list<TestMesh> 	meshes, Slic3r::Print &print, const DynamicPrintConfig& config, bool comments = false);
-void init_and_process_print(std::initializer_list<TriangleMesh> meshes, Slic3r::Print &print, const DynamicPrintConfig& config, bool comments = false);
-void init_and_process_print(std::initializer_list<TestMesh> 	meshes, Slic3r::Print &print, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false);
-void init_and_process_print(std::initializer_list<TriangleMesh> meshes, Slic3r::Print &print, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false);
+void init_and_process_print(
+    std::initializer_list<TestMesh> meshes,
+    Slic3r::Print &print,
+    const DynamicPrintConfig &config,
+    bool comments = false
+);
+void init_and_process_print(
+    std::initializer_list<TriangleMesh> meshes,
+    Slic3r::Print &print,
+    const DynamicPrintConfig &config,
+    bool comments = false
+);
+void init_and_process_print(
+    std::initializer_list<TestMesh> meshes,
+    Slic3r::Print &print,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false
+);
+void init_and_process_print(
+    std::initializer_list<TriangleMesh> meshes,
+    Slic3r::Print &print,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false
+);
 
-std::string gcode(Print& print);
+std::string gcode(Print &print);
 
-std::string slice(std::initializer_list<TestMesh> meshes, const DynamicPrintConfig &config, bool comments = false);
-std::string slice(std::initializer_list<TriangleMesh> meshes, const DynamicPrintConfig &config, bool comments = false);
-std::string slice(std::initializer_list<TestMesh> meshes, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false);
-std::string slice(std::initializer_list<TriangleMesh> meshes, std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items, bool comments = false);
+std::string slice(
+    std::initializer_list<TestMesh> meshes, const DynamicPrintConfig &config, bool comments = false
+);
+std::string slice(
+    std::initializer_list<TriangleMesh> meshes,
+    const DynamicPrintConfig &config,
+    bool comments = false
+);
+std::string slice(
+    std::initializer_list<TestMesh> meshes,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false
+);
+std::string slice(
+    std::initializer_list<TriangleMesh> meshes,
+    std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> config_items,
+    bool comments = false
+);
 
 bool contains(const std::string &data, const std::string &pattern);
 bool contains_regex(const std::string &data, const std::string &pattern);
 
-} } // namespace Slic3r::Test
+inline std::unique_ptr<Print> process_3mf(const std::filesystem::path &path) {
+    DynamicPrintConfig config;
+    auto print{std::make_unique<Print>()};
+    Model model;
 
+    ConfigSubstitutionContext context{ForwardCompatibilitySubstitutionRule::Disable};
+    load_3mf(path.string().c_str(), config, context, &model, false);
+
+    Slic3r::Test::init_print(std::vector<TriangleMesh>{}, *print, model, config);
+    print->process();
+
+    return print;
+}
+
+static std::map<std::string, std::unique_ptr<Print>> prints_3mfs;
+// Lazy getter, to avoid processing the 3mf multiple times, it already takes ages.
+inline Print *get_print(const std::filesystem::path &file_path) {
+    if (!prints_3mfs.count(file_path.string())) {
+        prints_3mfs[file_path.string()] = process_3mf(file_path.string());
+    }
+    return prints_3mfs[file_path.string()].get();
+}
+
+inline void serialize_seam(std::ostream &output, const std::vector<std::vector<Seams::SeamPerimeterChoice>> &seam) {
+    output << "x,y,z,layer_index" << std::endl;
+
+    for (const std::vector<Seams::SeamPerimeterChoice> &layer : seam) {
+        if (layer.empty()) {
+            continue;
+        }
+        const Seams::SeamPerimeterChoice &choice{layer.front()};
+
+        // clang-format off
+        output
+            << choice.choice.position.x() << ","
+            << choice.choice.position.y() << ","
+            << choice.perimeter.slice_z << ","
+            << choice.perimeter.layer_index << std::endl;
+        // clang-format on
+    }
+}
+
+struct SeamsFixture
+{
+    const std::filesystem::path file_3mf{
+        std::filesystem::path{TEST_DATA_DIR} / std::filesystem::path{"seam_test_object.3mf"}};
+    const Print *print{Test::get_print(file_3mf)};
+    const PrintObject *print_object{print->objects()[0]};
+
+    Seams::Params params{Seams::Placer::get_params(print->full_print_config())};
+
+    const Transform3d transformation{print_object->trafo_centered()};
+    const ModelVolumePtrs &volumes{print_object->model_object()->volumes};
+    Seams::ModelInfo::Painting painting{transformation, volumes};
+
+    const std::vector<Seams::Geometry::Extrusions> extrusions{
+        Seams::Geometry::get_extrusions(print_object->layers())};
+    const Seams::Perimeters::LayerInfos layer_infos{Seams::Perimeters::get_layer_infos(
+        print_object->layers(), params.perimeter.elephant_foot_compensation
+    )};
+    Seams::Shells::Shells<Polygon> shell_polygons{
+        Seams::Shells::create_shells(extrusions, params.max_distance)};
+
+    const std::size_t shell_index{15};
+
+    const ModelInfo::Visibility visibility{transformation, volumes, params.visibility, [](){}};
+    Seams::Aligned::VisibilityCalculator
+        visibility_calculator{visibility, params.convex_visibility_modifier, params.concave_visibility_modifier};
+};
+
+}} // namespace Slic3r::Test
 
 #endif // SLIC3R_TEST_DATA_HPP

tests/fff_print/test_seam_aligned.cpp

@@ -0,0 +1,159 @@
+#include <libslic3r/Point.hpp>
+#include <catch2/catch.hpp>
+#include <libslic3r/GCode/SeamAligned.hpp>
+#include "test_data.hpp"
+#include <fstream>
+
+using namespace Slic3r;
+using namespace Slic3r::Seams;
+
+constexpr bool debug_files{false};
+
+namespace AlignedTest {
+Perimeters::Perimeter get_perimeter() {
+    const double slice_z{1.0};
+    const std::size_t layer_index{};
+    std::vector<Vec2d> positions{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}, {0.0, 1.0}, {0.0, 0.5}};
+    std::vector<double> angles(positions.size(), -M_PI / 2.0);
+    angles[4] = 0.0;
+    std::vector<Perimeters::PointType> point_types(positions.size(), Perimeters::PointType::common);
+    std::vector<Perimeters::PointClassification>
+        point_classifications{positions.size(), Perimeters::PointClassification::common};
+    std::vector<Perimeters::AngleType> angle_type(positions.size(), Perimeters::AngleType::concave);
+    angle_type[4] = Perimeters::AngleType::smooth;
+
+    return {
+        slice_z,
+        layer_index,
+        std::move(positions),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_type)};
+}
+} // namespace AlignedTest
+
+TEST_CASE("Snap to angle", "[Seams][SeamAligned]") {
+    const Vec2d point{0.0, 0.4};
+    const std::size_t search_start{4};
+    const Perimeters::Perimeter perimeter{AlignedTest::get_perimeter()};
+
+    std::optional<std::size_t> snapped_to{
+        Aligned::Impl::snap_to_angle(point, search_start, perimeter, 0.5)};
+
+    REQUIRE(snapped_to);
+    CHECK(*snapped_to == 0);
+
+    snapped_to = Aligned::Impl::snap_to_angle(point, search_start, perimeter, 0.3);
+    REQUIRE(!snapped_to);
+}
+
+TEST_CASE("Get seam options", "[Seams][SeamAligned]") {
+    Perimeters::Perimeter perimeter{AlignedTest::get_perimeter()};
+    const Vec2d prefered_position{0.0, 0.3};
+
+    Aligned::Impl::SeamOptions options{Aligned::Impl::get_seam_options(
+        perimeter, prefered_position, *perimeter.common_points.common_points, 0.4
+    )};
+
+    CHECK(options.closest == 4);
+    CHECK(options.adjacent == 0);
+    CHECK((options.on_edge - Vec2d{0.0, 0.3}).norm() == Approx(0.0));
+    REQUIRE(options.snapped);
+    CHECK(options.snapped == 0);
+}
+
+struct PickSeamOptionFixture
+{
+    Perimeters::Perimeter perimeter{AlignedTest::get_perimeter()};
+
+    Aligned::Impl::SeamOptions options{
+        4,               // closest
+        0,               // adjacent
+        true,            // forward
+        false,           // snapped
+        Vec2d{0.0, 0.3}, // on_edge
+    };
+};
+
+TEST_CASE_METHOD(PickSeamOptionFixture, "Pick seam option", "[Seams][SeamAligned]") {
+    auto [previous_index, next_index, position]{pick_seam_option(perimeter, options)};
+    CHECK(previous_index == next_index);
+    CHECK((position - Vec2d{0.0, 0.0}).norm() == Approx(0.0));
+}
+
+TEST_CASE_METHOD(PickSeamOptionFixture, "Pick seam option picks enforcer", "[Seams][SeamAligned]") {
+    perimeter.point_types[4] = Perimeters::PointType::enforcer;
+
+    auto [previous_index, next_index, position]{pick_seam_option(perimeter, options)};
+    CHECK(previous_index == next_index);
+    CHECK((position - Vec2d{0.0, 0.5}).norm() == Approx(0.0));
+}
+
+TEST_CASE_METHOD(PickSeamOptionFixture, "Nearest point", "[Seams][SeamAligned]") {
+    const std::optional<SeamChoice> result{Aligned::Impl::Nearest{Vec2d{0.4, -0.1}, 0.2}(
+        perimeter, Perimeters::PointType::common, Perimeters::PointClassification::common
+    )};
+    CHECK(result->previous_index == 0);
+    CHECK(result->next_index == 1);
+    CHECK((result->position - Vec2d{0.4, 0.0}).norm() == Approx(0.0));
+}
+
+TEST_CASE_METHOD(PickSeamOptionFixture, "Least visible point", "[Seams][SeamAligned]") {
+    std::vector<double> precalculated_visibility{};
+    for (std::size_t i{0}; i < perimeter.positions.size(); ++i) {
+        precalculated_visibility.push_back(-static_cast<double>(i));
+    }
+    Aligned::Impl::LeastVisible least_visible{precalculated_visibility};
+    const std::optional<SeamChoice> result{least_visible(
+        perimeter, Perimeters::PointType::common, Perimeters::PointClassification::common
+    )};
+    CHECK(result->previous_index == 4);
+    CHECK(result->next_index == 4);
+    CHECK((result->position - Vec2d{0.0, 0.5}).norm() == Approx(0.0));
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Generate aligned seam", "[Seams][SeamAligned][Integration]") {
+    Shells::Shells<> perimeters{
+        Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+    Shells::Shells<> shell_perimeters;
+    shell_perimeters.push_back(std::move(perimeters[shell_index]));
+    const std::vector<std::vector<SeamPerimeterChoice>> seam{
+        Aligned::get_object_seams(std::move(shell_perimeters), visibility_calculator, params.aligned)};
+    REQUIRE(seam.size() == 125);
+
+    if constexpr (debug_files) {
+        std::ofstream csv{"aligned_seam.csv"};
+        Test::serialize_seam(csv, seam);
+    }
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Calculate visibility", "[Seams][SeamAligned][Integration]") {
+    if constexpr (debug_files) {
+        std::ofstream csv{"visibility.csv"};
+        csv << "x,y,z,visibility,total_visibility" << std::endl;
+        Shells::Shells<> perimeters{
+            Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+        for (const Shells::Shell<> &shell : perimeters) {
+            for (const Shells::Slice<> &slice : shell) {
+                for (std::size_t index{0}; index < slice.boundary.positions.size(); ++index) {
+                    const Vec2d &position{slice.boundary.positions[index]};
+                    const double point_visibility{visibility.calculate_point_visibility(
+                        to_3d(position.cast<float>(), slice.boundary.slice_z)
+                    )};
+                    const double total_visibility{
+                        visibility_calculator(SeamChoice{index, index, position}, slice.boundary)};
+
+                    // clang-format off
+                    csv <<
+                        position.x() << "," <<
+                        position.y() << "," <<
+                        slice.boundary.slice_z << "," <<
+                        point_visibility << "," <<
+                        total_visibility << std::endl;
+                    // clang-format on
+                }
+            }
+        }
+    }
+}

tests/fff_print/test_seam_geometry.cpp

@@ -0,0 +1,145 @@
+#include <libslic3r/Point.hpp>
+#include <catch2/catch.hpp>
+#include <libslic3r/GCode/SeamGeometry.hpp>
+#include <libslic3r/Geometry.hpp>
+
+using namespace Slic3r;
+
+TEST_CASE("Lists mapping", "[Seams][SeamGeometry]") {
+    // clang-format off
+    std::vector<std::vector<int>> list_of_lists{
+        {},
+        {7, 2, 3},
+        {9, 6, 3, 6, 7},
+        {1, 1, 3},
+        {1},
+        {3},
+        {1},
+        {},
+        {3}
+    };
+    // clang-format on
+
+    std::vector<std::size_t> sizes;
+    sizes.reserve(list_of_lists.size());
+    for (const std::vector<int> &list : list_of_lists) {
+        sizes.push_back(list.size());
+    }
+
+    const auto [mapping, bucket_cout]{Seams::Geometry::get_mapping(
+        sizes,
+        [&](const std::size_t layer_index,
+            const std::size_t item_index) -> Seams::Geometry::MappingOperatorResult {
+            unsigned max_diff{0};
+            std::optional<std::size_t> index;
+            const std::vector<int> &layer{list_of_lists[layer_index]};
+            const std::vector<int> &next_layer{list_of_lists[layer_index + 1]};
+            for (std::size_t i{0}; i < next_layer.size(); ++i) {
+                const long diff{std::abs(next_layer[i] - layer[item_index])};
+                if (diff > max_diff) {
+                    max_diff = diff;
+                    index = i;
+                }
+            }
+            if (!index) {
+                return std::nullopt;
+            }
+            return std::pair{*index, static_cast<double>(max_diff)};
+        }
+    )};
+
+    // clang-format off
+    CHECK(mapping == std::vector<std::vector<std::size_t>>{
+        {},
+        {0, 1, 2},
+        {1, 3, 0, 4, 5},
+        {1, 6, 7},
+        {7},
+        {7},
+        {7},
+        {},
+        {8}
+    });
+    // clang-format on
+}
+
+TEST_CASE("Vertex angle calculation counterclockwise", "[Seams][SeamGeometry]") {
+    std::vector<Vec2d> points{Vec2d{0, 0}, Vec2d{1, 0}, Vec2d{1, 1}, Vec2d{0, 1}};
+    std::vector<double> angles{Seams::Geometry::get_vertex_angles(points, 0.1)};
+
+    CHECK(angles.size() == 4);
+    for (const double angle : angles) {
+        CHECK(angle == Approx(-M_PI / 2));
+    }
+}
+
+TEST_CASE("Vertex angle calculation clockwise", "[Seams][SeamGeometry]") {
+    std::vector<Vec2d> points = {Vec2d{0, 0}, Vec2d{0, 1}, Vec2d{1, 1}, Vec2d{1, 0}};
+    std::vector<double> angles = Seams::Geometry::get_vertex_angles(points, 0.1);
+
+    CHECK(angles.size() == 4);
+    for (const double angle : angles) {
+        CHECK(angle == Approx(M_PI / 2));
+    }
+}
+
+TEST_CASE("Vertex angle calculation small convex", "[Seams][SeamGeometry]") {
+    std::vector<Vec2d> points = {Vec2d{0, 0}, Vec2d{-0.01, 1}, Vec2d{0, 2}, Vec2d{-2, 1}};
+    std::vector<double> angles = Seams::Geometry::get_vertex_angles(points, 0.1);
+
+    CHECK(angles.size() == 4);
+    CHECK(angles[1] > 0);
+    CHECK(angles[1] < 0.02);
+}
+
+TEST_CASE("Vertex angle calculation small concave", "[Seams][SeamGeometry]") {
+    std::vector<Vec2d> points = {Vec2d{0, 0}, Vec2d{0.01, 1}, Vec2d{0, 2}, Vec2d{-2, 1}};
+    std::vector<double> angles = Seams::Geometry::get_vertex_angles(points, 0.1);
+
+    CHECK(angles.size() == 4);
+    CHECK(angles[1] < 0);
+    CHECK(angles[1] > -0.02);
+}
+
+TEST_CASE("Vertex angle is rotation agnostic", "[Seams][SeamGeometry]") {
+    std::vector<Vec2d> points = {Vec2d{0, 0}, Vec2d{0.01, 1}, Vec2d{0, 2}, Vec2d{-2, 1}};
+    std::vector<double> angles = Seams::Geometry::get_vertex_angles(points, 0.1);
+
+    Points polygon_points;
+    using std::transform, std::back_inserter;
+    transform(points.begin(), points.end(), back_inserter(polygon_points), [](const Vec2d &point) {
+        return scaled(point);
+    });
+    Polygon polygon{polygon_points};
+    polygon.rotate(M_PI - Slic3r::Geometry::deg2rad(10.0));
+
+    std::vector<Vec2d> rotated_points;
+    using std::transform, std::back_inserter;
+    transform(
+        polygon.points.begin(), polygon.points.end(), back_inserter(rotated_points),
+        [](const Point &point) { return unscaled(point); }
+    );
+
+    std::vector<double> rotated_angles = Seams::Geometry::get_vertex_angles(points, 0.1);
+    CHECK(rotated_angles[1] == Approx(angles[1]));
+}
+
+TEST_CASE("Calculate overhangs", "[Seams][SeamGeometry]") {
+    const ExPolygon square{
+        scaled(Vec2d{0.0, 0.0}),
+        scaled(Vec2d{1.0, 0.0}),
+        scaled(Vec2d{1.0, 1.0}),
+        scaled(Vec2d{0.0, 1.0})
+    };
+    const std::vector<Vec2d> points{Seams::Geometry::unscaled(square.contour.points)};
+    ExPolygon previous_layer{square};
+    previous_layer.translate(scaled(Vec2d{-0.5, 0}));
+    AABBTreeLines::LinesDistancer<Linef> previous_layer_distancer{
+        to_unscaled_linesf({previous_layer})};
+    const std::vector<double> overhangs{
+        Seams::Geometry::get_overhangs(points, previous_layer_distancer, 0.5)};
+    REQUIRE(overhangs.size() == points.size());
+    CHECK_THAT(overhangs, Catch::Matchers::Approx(std::vector<double>{
+        0.0, M_PI / 4.0, M_PI / 4.0, 0.0
+    }));
+}

tests/fff_print/test_seam_perimeters.cpp

@@ -0,0 +1,180 @@
+#include "libslic3r/ClipperUtils.hpp"
+#include "libslic3r/GCode/SeamPerimeters.hpp"
+#include "libslic3r/Layer.hpp"
+#include "libslic3r/Point.hpp"
+#include <catch2/catch.hpp>
+#include <libslic3r/GCode/SeamGeometry.hpp>
+#include <libslic3r/Geometry.hpp>
+#include <fstream>
+
+#include "test_data.hpp"
+
+using namespace Slic3r;
+using namespace Slic3r::Seams;
+
+constexpr bool debug_files{false};
+
+const ExPolygon square{
+    scaled(Vec2d{0.0, 0.0}), scaled(Vec2d{1.0, 0.0}), scaled(Vec2d{1.0, 1.0}),
+    scaled(Vec2d{0.0, 1.0})};
+
+TEST_CASE("Oversample painted", "[Seams][SeamPerimeters]") {
+    auto is_painted{[](const Vec3f &position, float radius) {
+        return (position - Vec3f{0.5, 0.0, 1.0}).norm() < radius;
+    }};
+    std::vector<Vec2d> points{Perimeters::Impl::oversample_painted(
+        Seams::Geometry::unscaled(square.contour.points), is_painted, 1.0, 0.2
+    )};
+
+    REQUIRE(points.size() == 8);
+    CHECK((points[1] - Vec2d{0.2, 0.0}).norm() == Approx(0.0));
+
+    points = Perimeters::Impl::oversample_painted(
+        Seams::Geometry::unscaled(square.contour.points), is_painted, 1.0, 0.199
+    );
+    CHECK(points.size() == 9);
+}
+
+TEST_CASE("Remove redundant points", "[Seams][SeamPerimeters]") {
+    using Perimeters::PointType;
+    using Perimeters::PointClassification;
+
+    std::vector<Vec2d> points{{0.0, 0.0}, {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0},
+                              {3.0, 1.0}, {3.0, 2.0}, {0.0, 2.0}};
+    std::vector<PointType> point_types{PointType::common,
+                                       PointType::enforcer, // Should keep this.
+                                       PointType::enforcer, // Should keep this.
+                                       PointType::blocker,
+                                       PointType::blocker, // Should remove this.
+                                       PointType::blocker,  PointType::common};
+
+    const auto [resulting_points, resulting_point_types]{
+        Perimeters::Impl::remove_redundant_points(points, point_types, 0.1)};
+
+    REQUIRE(resulting_points.size() == 6);
+    REQUIRE(resulting_point_types.size() == 6);
+    CHECK((resulting_points[3] - Vec2d{3.0, 0.0}).norm() == Approx(0.0));
+    CHECK((resulting_points[4] - Vec2d{3.0, 2.0}).norm() == Approx(0.0));
+    CHECK(resulting_point_types[3] == PointType::blocker);
+    CHECK(resulting_point_types[4] == PointType::blocker);
+}
+
+TEST_CASE("Perimeter constructs KD trees", "[Seams][SeamPerimeters]") {
+    using Perimeters::PointType;
+    using Perimeters::PointClassification;
+    using Perimeters::AngleType;
+
+    std::vector<Vec2d> positions{Vec2d{0.0, 0.0}, Vec2d{1.0, 0.0}, Vec2d{1.0, 1.0}, Vec2d{0.0, 1.0}};
+    std::vector<double> angles(4, -M_PI / 2.0);
+    std::vector<PointType>
+        point_types{PointType::enforcer, PointType::blocker, PointType::common, PointType::common};
+    std::vector<PointClassification> point_classifications{
+        PointClassification::overhang, PointClassification::embedded, PointClassification::embedded,
+        PointClassification::common};
+    std::vector<AngleType>
+        angle_types{AngleType::convex, AngleType::concave, AngleType::smooth, AngleType::smooth};
+    Perimeters::Perimeter perimeter{
+        3.0,
+        2,
+        std::move(positions),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_types)};
+
+    CHECK(perimeter.enforced_points.overhanging_points);
+    CHECK(perimeter.blocked_points.embedded_points);
+    CHECK(perimeter.common_points.common_points);
+    CHECK(perimeter.common_points.embedded_points);
+}
+
+using std::filesystem::path;
+
+constexpr const char *to_string(Perimeters::PointType point_type) {
+    using Perimeters::PointType;
+
+    switch (point_type) {
+    case PointType::enforcer: return "enforcer";
+    case PointType::blocker: return "blocker";
+    case PointType::common: return "common";
+    }
+    throw std::runtime_error("Unreachable");
+}
+
+constexpr const char *to_string(Perimeters::PointClassification point_classification) {
+    using Perimeters::PointClassification;
+
+    switch (point_classification) {
+    case PointClassification::embedded: return "embedded";
+    case PointClassification::overhang: return "overhang";
+    case PointClassification::common: return "common";
+    }
+    throw std::runtime_error("Unreachable");
+}
+
+constexpr const char *to_string(Perimeters::AngleType angle_type) {
+    using Perimeters::AngleType;
+
+    switch (angle_type) {
+    case AngleType::convex: return "convex";
+    case AngleType::concave: return "concave";
+    case AngleType::smooth: return "smooth";
+    }
+    throw std::runtime_error("Unreachable");
+}
+
+void serialize_shell(std::ostream &output, const Shells::Shell<Perimeters::Perimeter> &shell) {
+    output << "x,y,z,point_type,point_classification,angle_type,layer_index,"
+              "point_index,distance,distance_to_previous,is_degenerate"
+           << std::endl;
+
+    for (std::size_t perimeter_index{0}; perimeter_index < shell.size(); ++perimeter_index) {
+        const Shells::Slice<> &slice{shell[perimeter_index]};
+        const Perimeters::Perimeter &perimeter{slice.boundary};
+        const std::vector<Vec2d> &points{perimeter.positions};
+
+        double total_distance{0.0};
+        for (std::size_t point_index{0}; point_index < perimeter.point_types.size(); ++point_index) {
+            const Vec3d point{to_3d(points[point_index], perimeter.slice_z)};
+            const Perimeters::PointType point_type{perimeter.point_types[point_index]};
+            const Perimeters::PointClassification point_classification{
+                perimeter.point_classifications[point_index]};
+            const Perimeters::AngleType angle_type{perimeter.angle_types[point_index]};
+            const std::size_t layer_index{slice.layer_index};
+            const std::size_t previous_index{point_index == 0 ? points.size() - 1 : point_index - 1};
+            const double distance_to_previous{(points[point_index] - points[previous_index]).norm()};
+            total_distance += point_index == 0 ? 0.0 : distance_to_previous;
+            const double distance{total_distance};
+            const bool is_degenerate{perimeter.is_degenerate};
+
+            // clang-format off
+                output
+                    << point.x() << ","
+                    << point.y() << ","
+                    << point.z() << ","
+                    << to_string(point_type) << ","
+                    << to_string(point_classification) << ","
+                    << to_string(angle_type) << ","
+                    << layer_index << ","
+                    << point_index << ","
+                    << distance << ","
+                    << distance_to_previous << ","
+                    << is_degenerate << std::endl;
+            // clang-format on
+        }
+    }
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Create perimeters", "[Seams][SeamPerimeters][Integration]") {
+    const Shells::Shells<> perimeters{
+        create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+
+    const Shells::Shell<> &shell{perimeters[shell_index]};
+
+    if constexpr (debug_files) {
+        std::ofstream csv{"perimeters.csv"};
+        serialize_shell(csv, shell);
+    }
+
+    REQUIRE(shell.size() == 54);
+}

tests/fff_print/test_seam_random.cpp

@@ -0,0 +1,98 @@
+#include <libslic3r/Point.hpp>
+#include <catch2/catch.hpp>
+#include <libslic3r/GCode/SeamRandom.hpp>
+#include "test_data.hpp"
+#include <fstream>
+
+using namespace Slic3r;
+using namespace Slic3r::Seams;
+
+constexpr bool debug_files{false};
+
+namespace RandomTest {
+Perimeters::Perimeter get_perimeter() {
+    const double slice_z{1.0};
+    const std::size_t layer_index{};
+    std::vector<Vec2d> positions{{0.0, 0.0}, {0.5, 0.0}, {1.0, 0.0}};
+    std::vector<double> angles(positions.size(), -M_PI / 2.0);
+    std::vector<Perimeters::PointType> point_types(positions.size(), Perimeters::PointType::common);
+    std::vector<Perimeters::PointClassification>
+        point_classifications{positions.size(), Perimeters::PointClassification::common};
+    std::vector<Perimeters::AngleType> angle_type(positions.size(), Perimeters::AngleType::concave);
+
+    return {
+        slice_z,
+        layer_index,
+        std::move(positions),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_type)};
+}
+} // namespace RandomTest
+
+double get_chi2_uniform(const std::vector<double> &data, double min, double max, const std::size_t bin_count) {
+    std::vector<std::size_t> bins(bin_count);
+    const double bin_size{(max - min) / bin_count};
+    const double expected_frequncy{static_cast<double>(data.size()) / bin_count};
+
+    for (const double value : data) {
+        auto bin{static_cast<int>(std::floor((value - min) / bin_size))};
+        bins[bin]++;
+    }
+
+    return std::accumulate(bins.begin(), bins.end(), 0.0, [&](const double total, const std::size_t count_in_bin){
+        return total + std::pow(static_cast<double>(count_in_bin - expected_frequncy), 2.0) / expected_frequncy;
+    });
+}
+
+TEST_CASE("Random is uniform", "[Seams][SeamRandom]") {
+    const int seed{42};
+    std::mt19937 random_engine{seed};
+    const Random::Impl::Random random{random_engine};
+    Perimeters::Perimeter perimeter{RandomTest::get_perimeter()};
+
+    std::vector<double> x_positions;
+    const std::size_t count{1001};
+    x_positions.reserve(count);
+    std::generate_n(std::back_inserter(x_positions), count, [&]() {
+        std::optional<SeamChoice> choice{
+            random(perimeter, Perimeters::PointType::common, Perimeters::PointClassification::common)};
+        return choice->position.x();
+    });
+    const std::size_t degrees_of_freedom{10};
+    const double critical{18.307}; // dof 10, significance 0.05
+
+    CHECK(get_chi2_uniform(x_positions, 0.0, 1.0, degrees_of_freedom + 1) < critical);
+}
+
+TEST_CASE("Random respects point type", "[Seams][SeamRandom]") {
+    const int seed{42};
+    std::mt19937 random_engine{seed};
+    const Random::Impl::Random random{random_engine};
+    Perimeters::Perimeter perimeter{RandomTest::get_perimeter()};
+    std::optional<SeamChoice> choice{
+        random(perimeter, Perimeters::PointType::common, Perimeters::PointClassification::common)};
+
+    REQUIRE(choice);
+    const std::size_t picked_index{choice->previous_index};
+    perimeter.point_types[picked_index] = Perimeters::PointType::blocker;
+    choice = random(perimeter, Perimeters::PointType::common, Perimeters::PointClassification::common);
+    REQUIRE(choice);
+    CHECK(choice->previous_index != picked_index);
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Generate random seam", "[Seams][SeamRandom][Integration]") {
+    Shells::Shells<> perimeters{
+        Seams::Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+    Shells::Shells<> shell_perimeters;
+    shell_perimeters.push_back(std::move(perimeters[shell_index]));
+    const std::vector<std::vector<SeamPerimeterChoice>> seam{
+        Random::get_object_seams(std::move(shell_perimeters), params.random_seed)};
+    REQUIRE(seam.size() == 125);
+
+    if constexpr (debug_files) {
+        std::ofstream csv{"random_seam.csv"};
+        Test::serialize_seam(csv, seam);
+    }
+}

tests/fff_print/test_seam_rear.cpp

@@ -0,0 +1,59 @@
+#include <libslic3r/Point.hpp>
+#include <catch2/catch.hpp>
+#include <libslic3r/GCode/SeamRear.hpp>
+#include "test_data.hpp"
+#include <fstream>
+
+using namespace Slic3r;
+using namespace Slic3r::Seams;
+
+constexpr bool debug_files{false};
+
+namespace RearTest {
+Perimeters::Perimeter get_perimeter() {
+    const double slice_z{1.0};
+    const std::size_t layer_index{};
+    std::vector<Vec2d> positions{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}, {0.5, 1.0}, {0.0, 1.0}};
+    std::vector<double> angles(positions.size(), -M_PI / 2.0);
+    angles[3] = 0.0;
+    std::vector<Perimeters::PointType> point_types(positions.size(), Perimeters::PointType::common);
+    std::vector<Perimeters::PointClassification>
+        point_classifications{positions.size(), Perimeters::PointClassification::common};
+    std::vector<Perimeters::AngleType> angle_type(positions.size(), Perimeters::AngleType::concave);
+    angle_type[3] = Perimeters::AngleType::smooth;
+
+    return {
+        slice_z,
+        layer_index,
+        std::move(positions),
+        std::move(angles),
+        std::move(point_types),
+        std::move(point_classifications),
+        std::move(angle_type)};
+}
+} // namespace RearTest
+
+TEST_CASE("StraightLine operator places seam point near the prefered position", "[Seams][SeamRear]") {
+    const Rear::Impl::StraightLine rearest{Vec2d{0.7, 2.0}};
+    std::optional<SeamChoice> choice{rearest(RearTest::get_perimeter(), Perimeters::PointType::common, Perimeters::PointClassification::common)};
+
+    REQUIRE(choice);
+    CHECK(scaled(choice->position) == scaled(Vec2d{0.7, 1.0}));
+    CHECK(choice->previous_index == 2);
+    CHECK(choice->next_index == 3);
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Generate rear seam", "[Seams][SeamRear][Integration]") {
+    Shells::Shells<> perimeters{
+        Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+    Shells::Shells<> shell_perimeters;
+    shell_perimeters.push_back(std::move(perimeters[shell_index]));
+    const std::vector<std::vector<SeamPerimeterChoice>> seam{
+        Rear::get_object_seams(std::move(shell_perimeters), params.rear_project_threshold)};
+    REQUIRE(seam.size() == 125);
+
+    if constexpr (debug_files) {
+        std::ofstream csv{"rear_seam.csv"};
+        Test::serialize_seam(csv, seam);
+    }
+}

tests/fff_print/test_seam_shells.cpp

@@ -0,0 +1,73 @@
+#include <catch2/catch.hpp>
+#include <filesystem>
+#include <fstream>
+#include "libslic3r/ClipperUtils.hpp"
+#include "libslic3r/GCode/SeamPainting.hpp"
+#include "test_data.hpp"
+
+#include "libslic3r/GCode/SeamShells.hpp"
+
+using namespace Slic3r;
+using namespace Slic3r::Seams;
+
+constexpr bool debug_files{false};
+
+struct ProjectionFixture
+{
+    Polygon extrusion_path{
+        Point{scaled(Vec2d{-1.0, -1.0})}, Point{scaled(Vec2d{1.0, -1.0})},
+        Point{scaled(Vec2d{1.0, 1.0})}, Point{scaled(Vec2d{-1.0, 1.0})}};
+
+    ExPolygon island_boundary;
+    Seams::Geometry::Extrusions extrusions;
+    double extrusion_width{0.2};
+
+    ProjectionFixture() {
+        extrusions.emplace_back(extrusion_path.bounding_box(), extrusion_width, island_boundary);
+    }
+};
+
+TEST_CASE_METHOD(ProjectionFixture, "Project to geometry matches", "[Seams][SeamShells]") {
+    Polygon boundary_polygon{extrusion_path};
+    // Add + 0.1 to check that boundary polygon has been picked.
+    boundary_polygon.scale(1.0 + extrusion_width / 2.0 + 0.1);
+    island_boundary.contour = boundary_polygon;
+
+    Shells::Impl::BoundedPolygons result{Shells::Impl::project_to_geometry(extrusions)};
+    REQUIRE(result.size() == 1);
+    REQUIRE(result[0].polygon.size() == 4);
+    // Boundary polygon is picked.
+    CHECK(result[0].polygon[0].x() == Approx(scaled(-(1.0 + extrusion_width / 2.0 + 0.1))));
+}
+
+void serialize_shells(
+    std::ostream &out, const Shells::Shells<Polygon> &shells, const double layer_height
+) {
+    out << "x,y,z,layer_index,slice_id,shell_id" << std::endl;
+    for (std::size_t shell_id{}; shell_id < shells.size(); ++shell_id) {
+        const Shells::Shell<Polygon> &shell{shells[shell_id]};
+        for (std::size_t slice_id{}; slice_id < shell.size(); ++slice_id) {
+            const Shells::Slice<Polygon> &slice{shell[slice_id]};
+            for (const Point &point : slice.boundary) {
+                // clang-format off
+                out
+                    << point.x() << ","
+                    << point.y() << ","
+                    << slice.layer_index * 1e6 * layer_height << ","
+                    << slice.layer_index << ","
+                    << slice_id << ","
+                    << shell_id << std::endl;
+                // clang-format on
+            }
+        }
+    }
+}
+
+TEST_CASE_METHOD(Test::SeamsFixture, "Create shells", "[Seams][SeamShells][Integration]") {
+    if constexpr (debug_files) {
+        std::ofstream csv{"shells.csv"};
+        serialize_shells(csv, shell_polygons, print->full_print_config().opt_float("layer_height"));
+    }
+
+    CHECK(shell_polygons.size() == 39);
+}