diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 5919143126..1ec0e1c8a0 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/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
diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp
index 6a7f825815..954f96fc3e 100644
--- a/src/libslic3r/GCode.cpp
+++ b/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.
diff --git a/src/libslic3r/GCode.hpp b/src/libslic3r/GCode.hpp
index 3aa59fb286..e8342d2547 100644
--- a/src/libslic3r/GCode.hpp
+++ b/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()
diff --git a/src/libslic3r/GCode/ModelVisibility.cpp b/src/libslic3r/GCode/ModelVisibility.cpp
new file mode 100644
index 0000000000..63f3beffc3
--- /dev/null
+++ b/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;
+}
+
+}
diff --git a/src/libslic3r/GCode/ModelVisibility.hpp b/src/libslic3r/GCode/ModelVisibility.hpp
new file mode 100644
index 0000000000..a6ce29fba2
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamAligned.cpp b/src/libslic3r/GCode/SeamAligned.cpp
new file mode 100644
index 0000000000..a2e22c6b33
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamAligned.hpp b/src/libslic3r/GCode/SeamAligned.hpp
new file mode 100644
index 0000000000..d3fbebab6c
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamChoice.cpp b/src/libslic3r/GCode/SeamChoice.cpp
new file mode 100644
index 0000000000..6db562e2a1
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamChoice.hpp b/src/libslic3r/GCode/SeamChoice.hpp
new file mode 100644
index 0000000000..d33933b71f
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamGeometry.cpp b/src/libslic3r/GCode/SeamGeometry.cpp
new file mode 100644
index 0000000000..2a90668215
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamGeometry.hpp b/src/libslic3r/GCode/SeamGeometry.hpp
new file mode 100644
index 0000000000..219fd1a4fe
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamPainting.cpp b/src/libslic3r/GCode/SeamPainting.cpp
new file mode 100644
index 0000000000..1f78eb5098
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamPainting.hpp b/src/libslic3r/GCode/SeamPainting.hpp
new file mode 100644
index 0000000000..a6fa6ff8d0
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamPerimeters.cpp b/src/libslic3r/GCode/SeamPerimeters.cpp
new file mode 100644
index 0000000000..e83b694afd
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamPerimeters.hpp b/src/libslic3r/GCode/SeamPerimeters.hpp
new file mode 100644
index 0000000000..20705754fd
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp
index ed0d480027..21b42584e0 100644
--- a/src/libslic3r/GCode/SeamPlacer.cpp
+++ b/src/libslic3r/GCode/SeamPlacer.cpp
@@ -1,1616 +1,323 @@
-///|/ Copyright (c) Prusa Research 2020 - 2023 Vojtěch Bubník @bubnikv, Lukáš Matěna @lukasmatena, Pavel Mikuš @Godrak
+///|/ Copyright (c) Prusa Research 2020 - 2023 Vojtěch Bubník @bubnikv, Lukáš Matěna @lukasmatena,
+/// Pavel Mikuš @Godrak
 ///|/ Copyright (c) SuperSlicer 2023 Remi Durand @supermerill
 ///|/
 ///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
 ///|/
+
+#include <boost/filesystem/operations.hpp>
+#include <boost/property_tree/json_parser.hpp>
 #include "SeamPlacer.hpp"
 
-#include "Color.hpp"
-#include "Polygon.hpp"
-#include "PrintConfig.hpp"
-#include "tbb/parallel_for.h"
-#include "tbb/blocked_range.h"
-#include "tbb/parallel_reduce.h"
-#include <boost/log/trivial.hpp>
-#include <random>
-#include <algorithm>
-#include <queue>
+#include "libslic3r/GCode/SeamAligned.hpp"
+#include "libslic3r/GCode/SeamRear.hpp"
+#include "libslic3r/GCode/SeamRandom.hpp"
+#include "libslic3r/GCode/ModelVisibility.hpp"
+#include "libslic3r/GCode/SeamGeometry.hpp"
 
-#include "libslic3r/AABBTreeLines.hpp"
-#include "libslic3r/KDTreeIndirect.hpp"
-#include "libslic3r/ExtrusionEntity.hpp"
-#include "libslic3r/Print.hpp"
-#include "libslic3r/BoundingBox.hpp"
-#include "libslic3r/ClipperUtils.hpp"
-#include "libslic3r/Layer.hpp"
+namespace Slic3r::Seams {
 
-#include "libslic3r/Geometry/Curves.hpp"
-#include "libslic3r/ShortEdgeCollapse.hpp"
-#include "libslic3r/TriangleSetSampling.hpp"
+using ObjectShells = std::vector<std::pair<const PrintObject *, Shells::Shells<>>>;
+using ObjectPainting = std::map<const PrintObject*, ModelInfo::Painting>;
 
-#include "libslic3r/Utils.hpp"
+ObjectShells partition_to_shells(
+    SpanOfConstPtrs<PrintObject> objects,
+    const Params &params,
+    const ObjectPainting& object_painting,
+    const std::function<void(void)> &throw_if_canceled
+) {
+    ObjectShells result;
 
-//#define DEBUG_FILES
+    for (const PrintObject *print_object : objects) {
+        const ModelInfo::Painting &painting{object_painting.at(print_object)};
+        throw_if_canceled();
 
-#ifdef DEBUG_FILES
-#include <boost/nowide/cstdio.hpp>
-#include <SVG.hpp>
-#endif
+        const std::vector<Geometry::Extrusions> extrusions{
+            Geometry::get_extrusions(print_object->layers())};
+        const Perimeters::LayerInfos layer_infos{Perimeters::get_layer_infos(
+            print_object->layers(), params.perimeter.elephant_foot_compensation
+        )};
+        Shells::Shells<Polygon> shell_polygons{
+            Shells::create_shells(extrusions, params.max_distance)};
 
-namespace Slic3r {
-
-namespace SeamPlacerImpl {
-
-template<typename T> int sgn(T val) {
-    return int(T(0) < val) - int(val < T(0));
-}
-
-// base function: ((e^(((1)/(x^(2)+1)))-1)/(e-1))
-// checkout e.g. here: https://www.geogebra.org/calculator
-float gauss(float value, float mean_x_coord, float mean_value, float falloff_speed) {
-    float shifted = value - mean_x_coord;
-    float denominator = falloff_speed * shifted * shifted + 1.0f;
-    float exponent = 1.0f / denominator;
-    return mean_value * (std::exp(exponent) - 1.0f) / (std::exp(1.0f) - 1.0f);
-}
-
-float compute_angle_penalty(float ccw_angle) {
-    // This function is used:
-    // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x)))
-    // looks scary, but it is gaussian combined with sigmoid,
-    // so that concave points have much smaller penalty over convex ones
-    // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png
-    return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) +
-            1.0f / (2 + std::exp(-ccw_angle));
-}
-
-/// Coordinate frame
-class Frame {
-public:
-    Frame() {
-        mX = Vec3f(1, 0, 0);
-        mY = Vec3f(0, 1, 0);
-        mZ = Vec3f(0, 0, 1);
+        Shells::Shells<> perimeters{
+            Perimeters::create_perimeters(shell_polygons, layer_infos, painting, params.perimeter)};
+        throw_if_canceled();
+        result.emplace_back(print_object, std::move(perimeters));
     }
-
-    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) {
-    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 / SeamPlacer::sqr_rays_per_sample_point;
-    std::vector<Vec3f> precomputed_sample_directions(
-            SeamPlacer::sqr_rays_per_sample_point * SeamPlacer::sqr_rays_per_sample_point);
-    for (size_t x_idx = 0; x_idx < SeamPlacer::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 < SeamPlacer::sqr_rays_per_sample_point; ++y_idx) {
-            size_t dir_index = x_idx * SeamPlacer::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](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;
-                    constexpr float decrease_step = 1.0f
-                            / (SeamPlacer::sqr_rays_per_sample_point * SeamPlacer::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;
 }
 
-std::vector<float> calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector<float> &lengths,
-        float min_arm_length) {
-    std::vector<float> result(polygon.size());
-
-    if (polygon.size() == 1) {
-        result[0] = 0.0f;
-    }
-
-    size_t idx_prev = 0;
-    size_t idx_curr = 0;
-    size_t idx_next = 0;
-
-    float distance_to_prev = 0;
-    float distance_to_next = 0;
-
-    //push idx_prev far enough back as initialization
-    while (distance_to_prev < min_arm_length) {
-        idx_prev = Slic3r::prev_idx_modulo(idx_prev, polygon.size());
-        distance_to_prev += lengths[idx_prev];
-    }
-
-    for (size_t _i = 0; _i < polygon.size(); ++_i) {
-        // pull idx_prev to current as much as possible, while respecting the min_arm_length
-        while (distance_to_prev - lengths[idx_prev] > min_arm_length) {
-            distance_to_prev -= lengths[idx_prev];
-            idx_prev = Slic3r::next_idx_modulo(idx_prev, polygon.size());
-        }
-
-        //push idx_next forward as far as needed
-        while (distance_to_next < min_arm_length) {
-            distance_to_next += lengths[idx_next];
-            idx_next = Slic3r::next_idx_modulo(idx_next, polygon.size());
-        }
-
-        // Calculate angle between idx_prev, idx_curr, idx_next.
-        const Point &p0 = polygon.points[idx_prev];
-        const Point &p1 = polygon.points[idx_curr];
-        const Point &p2 = polygon.points[idx_next];
-        result[idx_curr] = float(angle(p1 - p0, p2 - p1));
-
-        // increase idx_curr by one
-        float curr_distance = lengths[idx_curr];
-        idx_curr++;
-        distance_to_prev += curr_distance;
-        distance_to_next -= curr_distance;
-    }
-
-    return result;
-}
-
-struct CoordinateFunctor {
-    const std::vector<Vec3f> *coordinates;
-    CoordinateFunctor(const std::vector<Vec3f> *coords) :
-            coordinates(coords) {
-    }
-    CoordinateFunctor() :
-            coordinates(nullptr) {
-    }
-
-    const float& operator()(size_t idx, size_t dim) const {
-        return coordinates->operator [](idx)[dim];
-    }
-};
-
-// structure to store global information about the model - occlusion hits, enforcers, blockers
-struct GlobalModelInfo {
-    TriangleSetSamples mesh_samples;
-    std::vector<float> mesh_samples_visibility;
-    CoordinateFunctor mesh_samples_coordinate_functor;
-    KDTreeIndirect<3, float, CoordinateFunctor> mesh_samples_tree { CoordinateFunctor { } };
-    float mesh_samples_radius;
-
-    indexed_triangle_set enforcers;
-    indexed_triangle_set blockers;
-    AABBTreeIndirect::Tree<3, float> enforcers_tree;
-    AABBTreeIndirect::Tree<3, float> blockers_tree;
-
-    bool is_enforced(const Vec3f &position, float radius) const {
-        if (enforcers.empty()) {
-            return false;
-        }
-        float radius_sqr = radius * radius;
-        return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices,
-                enforcers_tree, position, radius_sqr);
-    }
-
-    bool is_blocked(const Vec3f &position, float radius) const {
-        if (blockers.empty()) {
-            return false;
-        }
-        float radius_sqr = radius * radius;
-        return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices,
-                blockers_tree, position, radius_sqr);
-    }
-
-    float 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;
-
-    }
-
-#ifdef DEBUG_FILES
-    void debug_export(const indexed_triangle_set &obj_mesh) const {
-
-        indexed_triangle_set divided_mesh = obj_mesh;
-        Slic3r::CNumericLocalesSetter locales_setter;
-
-        {
-            auto filename = debug_out_path("visiblity.obj");
-            FILE *fp = boost::nowide::fopen(filename.c_str(), "w");
-            if (fp == nullptr) {
-                BOOST_LOG_TRIVIAL(error)
-                << "stl_write_obj: Couldn't open " << filename << " for writing";
-                return;
-            }
-
-            for (size_t i = 0; i < divided_mesh.vertices.size(); ++i) {
-                float visibility = calculate_point_visibility(divided_mesh.vertices[i]);
-                Vec3f color = value_to_rgbf(0.0f, 1.0f, visibility);
-                fprintf(fp, "v %f %f %f  %f %f %f\n",
-                        divided_mesh.vertices[i](0), divided_mesh.vertices[i](1), divided_mesh.vertices[i](2),
-                        color(0), color(1), color(2));
-            }
-            for (size_t i = 0; i < divided_mesh.indices.size(); ++i)
-                fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1,
-                        divided_mesh.indices[i][2] + 1);
-            fclose(fp);
-        }
-
-        {
-            auto filename = debug_out_path("visiblity_samples.obj");
-            FILE *fp = boost::nowide::fopen(filename.c_str(), "w");
-            if (fp == nullptr) {
-                BOOST_LOG_TRIVIAL(error)
-                << "stl_write_obj: Couldn't open " << filename << " for writing";
-                return;
-            }
-
-            for (size_t i = 0; i < mesh_samples.positions.size(); ++i) {
-                float visibility = mesh_samples_visibility[i];
-                Vec3f color = value_to_rgbf(0.0f, 1.0f, visibility);
-                fprintf(fp, "v %f %f %f  %f %f %f\n",
-                        mesh_samples.positions[i](0), mesh_samples.positions[i](1), mesh_samples.positions[i](2),
-                        color(0), color(1), color(2));
-            }
-            fclose(fp);
-        }
-
-    }
-#endif
-}
-;
-
-//Extract perimeter polygons of the given layer
-Polygons extract_perimeter_polygons(const Layer *layer, std::vector<const LayerRegion*> &corresponding_regions_out) {
-    Polygons polygons;
-    for (const LayerRegion *layer_region : layer->regions()) {
-        for (const ExtrusionEntity *ex_entity : layer_region->perimeters()) {
-            if (ex_entity->is_collection()) { //collection of inner, outer, and overhang perimeters
-                for (const ExtrusionEntity *perimeter : static_cast<const ExtrusionEntityCollection*>(ex_entity)->entities) {
-                    ExtrusionRole role = perimeter->role();
-                    if (perimeter->is_loop()) {
-                        for (const ExtrusionPath &path : static_cast<const ExtrusionLoop*>(perimeter)->paths) {
-                            if (path.role() == ExtrusionRole::ExternalPerimeter) {
-                                role = ExtrusionRole::ExternalPerimeter;
-                            }
-                        }
-                    }
-
-                    if (role == ExtrusionRole::ExternalPerimeter) {
-                        Points p;
-                        perimeter->collect_points(p);
-                        polygons.emplace_back(std::move(p));
-                        corresponding_regions_out.push_back(layer_region);
-                    }
-                }
-                if (polygons.empty()) {
-                    Points p;
-                    ex_entity->collect_points(p);
-                    polygons.emplace_back(std::move(p));
-                    corresponding_regions_out.push_back(layer_region);
-                }
-            } else {
-                Points p;
-                ex_entity->collect_points(p);
-                polygons.emplace_back(std::move(p));
-                corresponding_regions_out.push_back(layer_region);
-            }
-        }
-    }
-
-    if (polygons.empty()) { // If there are no perimeter polygons for whatever reason (disabled perimeters .. ) insert dummy point
-        // it is easier than checking everywhere if the layer is not emtpy, no seam will be placed to this layer anyway
-        polygons.emplace_back(Points{ { 0, 0 } });
-        corresponding_regions_out.push_back(nullptr);
-    }
-
-    return polygons;
-}
-
-// Insert SeamCandidates created from perimeter polygons in to the result vector.
-// Compute its type (Enfrocer,Blocker), angle, and position
-//each SeamCandidate also contains pointer to shared Perimeter structure representing the polygon
-// if Custom Seam modifiers are present, oversamples the polygon if necessary to better fit user intentions
-void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const LayerRegion *region,
-        const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) {
-    if (orig_polygon.size() == 0) {
-        return;
-    }
-    Polygon polygon = orig_polygon;
-    bool was_clockwise = polygon.make_counter_clockwise();
-    float angle_arm_len = region != nullptr ? region->flow(FlowRole::frExternalPerimeter).nozzle_diameter() : 0.5f;
-
-    std::vector<float> lengths { };
-    for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) {
-        lengths.push_back((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm());
-    }
-    lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.1));
-    std::vector<float> polygon_angles = calculate_polygon_angles_at_vertices(polygon, lengths,
-            angle_arm_len);
-
-    result.perimeters.push_back( { });
-    Perimeter &perimeter = result.perimeters.back();
-
-    std::queue<Vec3f> orig_polygon_points { };
-    for (size_t index = 0; index < polygon.size(); ++index) {
-        Vec2f unscaled_p = unscale(polygon[index]).cast<float>();
-        orig_polygon_points.emplace(unscaled_p.x(), unscaled_p.y(), z_coord);
-    }
-    Vec3f first = orig_polygon_points.front();
-    std::queue<Vec3f> oversampled_points { };
-    size_t orig_angle_index = 0;
-    perimeter.start_index = result.points.size();
-    perimeter.flow_width = region != nullptr ? region->flow(FlowRole::frExternalPerimeter).width() : 0.0f;
-    bool some_point_enforced = false;
-    while (!orig_polygon_points.empty() || !oversampled_points.empty()) {
-        EnforcedBlockedSeamPoint type = EnforcedBlockedSeamPoint::Neutral;
-        Vec3f position;
-        float local_ccw_angle = 0;
-        bool orig_point = false;
-        if (!oversampled_points.empty()) {
-            position = oversampled_points.front();
-            oversampled_points.pop();
-        } else {
-            position = orig_polygon_points.front();
-            orig_polygon_points.pop();
-            local_ccw_angle = was_clockwise ? -polygon_angles[orig_angle_index] : polygon_angles[orig_angle_index];
-            orig_angle_index++;
-            orig_point = true;
-        }
-
-        if (global_model_info.is_enforced(position, perimeter.flow_width)) {
-            type = EnforcedBlockedSeamPoint::Enforced;
-        }
-
-        if (global_model_info.is_blocked(position, perimeter.flow_width)) {
-            type = EnforcedBlockedSeamPoint::Blocked;
-        }
-        some_point_enforced = some_point_enforced || type == EnforcedBlockedSeamPoint::Enforced;
-
-        if (orig_point) {
-            Vec3f pos_of_next = orig_polygon_points.empty() ? first : orig_polygon_points.front();
-            float distance_to_next = (position - pos_of_next).norm();
-            if (global_model_info.is_enforced(position, distance_to_next)) {
-                Vec3f vec_to_next = (pos_of_next - position).normalized();
-                float step_size = SeamPlacer::enforcer_oversampling_distance;
-                float step = step_size;
-                while (step < distance_to_next) {
-                    oversampled_points.push(position + vec_to_next * step);
-                    step += step_size;
-                }
-            }
-        }
-
-        result.points.emplace_back(position, perimeter, local_ccw_angle, type);
-    }
-
-    perimeter.end_index = result.points.size();
-
-    if (some_point_enforced) {
-        // We will patches of enforced points (patch: continuous section of enforced points), choose
-        // the longest patch, and select the middle point or sharp point (depending on the angle)
-        // this point will have high priority on this perimeter
-        size_t perimeter_size = perimeter.end_index - perimeter.start_index;
-        const auto next_index = [&](size_t idx) {
-            return perimeter.start_index + Slic3r::next_idx_modulo(idx - perimeter.start_index, perimeter_size);
-        };
-
-        std::vector<size_t> patches_starts_ends;
-        for (size_t i = perimeter.start_index; i < perimeter.end_index; ++i) {
-            if (result.points[i].type != EnforcedBlockedSeamPoint::Enforced &&
-                    result.points[next_index(i)].type == EnforcedBlockedSeamPoint::Enforced) {
-                patches_starts_ends.push_back(next_index(i));
-            }
-            if (result.points[i].type == EnforcedBlockedSeamPoint::Enforced &&
-                    result.points[next_index(i)].type != EnforcedBlockedSeamPoint::Enforced) {
-                patches_starts_ends.push_back(next_index(i));
-            }
-        }
-        //if patches_starts_ends are empty, it means that the whole perimeter is enforced.. don't do anything in that case
-        if (!patches_starts_ends.empty()) {
-            //if the first point in the patches is not enforced, it marks a patch end. in that case, put it to the end and start on next
-            // to simplify the processing
-            assert(patches_starts_ends.size() % 2 == 0);
-            bool start_on_second = false;
-            if (result.points[patches_starts_ends[0]].type != EnforcedBlockedSeamPoint::Enforced) {
-                start_on_second = true;
-                patches_starts_ends.push_back(patches_starts_ends[0]);
-            }
-            //now pick the longest patch
-            std::pair<size_t, size_t> longest_patch { 0, 0 };
-            auto patch_len = [perimeter_size](const std::pair<size_t, size_t> &start_end) {
-                if (start_end.second < start_end.first) {
-                    return start_end.first + (perimeter_size - start_end.second);
-                } else {
-                    return start_end.second - start_end.first;
-                }
-            };
-            for (size_t patch_idx = start_on_second ? 1 : 0; patch_idx < patches_starts_ends.size(); patch_idx += 2) {
-                std::pair<size_t, size_t> current_patch { patches_starts_ends[patch_idx], patches_starts_ends[patch_idx
-                        + 1] };
-                if (patch_len(longest_patch) < patch_len(current_patch)) {
-                    longest_patch = current_patch;
-                }
-            }
-            std::vector<size_t> viable_points_indices;
-            std::vector<size_t> large_angle_points_indices;
-            for (size_t point_idx = longest_patch.first; point_idx != longest_patch.second;
-                    point_idx = next_index(point_idx)) {
-                viable_points_indices.push_back(point_idx);
-                if (std::abs(result.points[point_idx].local_ccw_angle)
-                        > SeamPlacer::sharp_angle_snapping_threshold) {
-                    large_angle_points_indices.push_back(point_idx);
-                }
-            }
-            assert(viable_points_indices.size() > 0);
-            if (large_angle_points_indices.empty()) {
-                size_t central_idx = viable_points_indices[viable_points_indices.size() / 2];
-                result.points[central_idx].central_enforcer = true;
-            } else {
-                size_t central_idx = large_angle_points_indices.size() / 2;
-                result.points[large_angle_points_indices[central_idx]].central_enforcer = true;
-            }
-        }
-    }
-
-}
-
-// Get index of previous and next perimeter point of the layer. Because SeamCandidates of all polygons of the given layer
-// are sequentially stored in the vector, each perimeter contains info about start and end index. These vales are used to
-// deduce index of previous and next neigbour in the corresponding perimeter.
-std::pair<size_t, size_t> find_previous_and_next_perimeter_point(const std::vector<SeamCandidate> &perimeter_points,
-        size_t point_index) {
-    const SeamCandidate &current = perimeter_points[point_index];
-    int prev = point_index - 1; //for majority of points, it is true that neighbours lie behind and in front of them in the vector
-    int next = point_index + 1;
-
-    if (point_index == current.perimeter.start_index) {
-        // if point_index is equal to start, it means that the previous neighbour is at the end
-        prev = current.perimeter.end_index;
-    }
-
-    if (point_index == current.perimeter.end_index - 1) {
-        // if point_index is equal to end, than next neighbour is at the start
-        next = current.perimeter.start_index;
-    }
-
-    assert(prev >= 0);
-    assert(next >= 0);
-    return {size_t(prev),size_t(next)};
-}
-
-// Computes all global model info - transforms object, performs raycasting
-void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po,
-        std::function<void(void)> throw_if_canceled) {
-    BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: gather occlusion meshes: start";
-    auto obj_transform = po->trafo_centered();
-    indexed_triangle_set triangle_set;
-    indexed_triangle_set negative_volumes_set;
-    //add all parts
-    for (const ModelVolume *model_volume : po->model_object()->volumes) {
-        if (model_volume->type() == ModelVolumeType::MODEL_PART
-                || model_volume->type() == ModelVolumeType::NEGATIVE_VOLUME) {
-            auto model_transformation = model_volume->get_matrix();
-            indexed_triangle_set model_its = model_volume->mesh().its;
-            its_transform(model_its, model_transformation);
-            if (model_volume->type() == ModelVolumeType::MODEL_PART) {
-                its_merge(triangle_set, model_its);
-            } else {
-                its_merge(negative_volumes_set, model_its);
-            }
-        }
-    }
-    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, SeamPlacer::fast_decimation_triangle_count_target);
-    its_short_edge_collpase(negative_volumes_set, SeamPlacer::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";
-
-    result.mesh_samples = sample_its_uniform_parallel(SeamPlacer::raycasting_visibility_samples_count,
-            triangle_set);
-    result.mesh_samples_coordinate_functor = CoordinateFunctor(&result.mesh_samples.positions);
-    result.mesh_samples_tree = KDTreeIndirect<3, float, CoordinateFunctor>(result.mesh_samples_coordinate_functor,
-            result.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 = SeamPlacer::raycasting_visibility_samples_count / result.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);
-    result.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: " << result.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";
-    result.mesh_samples_visibility = raycast_visibility(raycasting_tree, triangle_set, result.mesh_samples,
-            negative_volumes_start_index);
-    throw_if_canceled();
-#ifdef DEBUG_FILES
-    result.debug_export(triangle_set);
-#endif
-}
-
-void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) {
-    BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: start";
-
-    auto obj_transform = po->trafo_centered();
-
-    for (const ModelVolume *mv : po->model_object()->volumes) {
-        if (mv->is_seam_painted()) {
-            auto model_transformation = obj_transform * mv->get_matrix();
-
-            indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER);
-            its_transform(enforcers, model_transformation);
-            its_merge(result.enforcers, enforcers);
-
-            indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER);
-            its_transform(blockers, model_transformation);
-            its_merge(result.blockers, blockers);
-        }
-    }
-
-    result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices,
-            result.enforcers.indices);
-    result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices,
-            result.blockers.indices);
-
-    BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end";
-}
-
-struct SeamComparator {
-    SeamPosition setup;
-    float angle_importance;
-    explicit SeamComparator(SeamPosition setup) :
-            setup(setup) {
-        angle_importance =
-                setup == spNearest ? SeamPlacer::angle_importance_nearest : SeamPlacer::angle_importance_aligned;
-    }
-
-    // Standard comparator, must respect the requirements of comparators (e.g. give same result on same inputs) for sorting usage
-    // should return if a is better seamCandidate than b
-    bool is_first_better(const SeamCandidate &a, const SeamCandidate &b, const Vec2f &preffered_location = Vec2f { 0.0f,
-            0.0f }) const {
-        if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) {
-            return a.central_enforcer;
-        }
-
-        // Blockers/Enforcers discrimination, top priority
-        if (a.type != b.type) {
-            return a.type > b.type;
-        }
-
-        //avoid overhangs
-        if (a.overhang > 0.0f || b.overhang > 0.0f) {
-            return a.overhang < b.overhang;
-        }
-
-        // prefer hidden points (more than 0.5 mm inside)
-        if (a.embedded_distance < -0.5f && b.embedded_distance > -0.5f) {
-            return true;
-        }
-        if (b.embedded_distance < -0.5f && a.embedded_distance > -0.5f) {
-            return false;
-        }
-
-        if (setup == SeamPosition::spRear && a.position.y() != b.position.y()) {
-            return a.position.y() > b.position.y();
-        }
-
-        float distance_penalty_a = 0.0f;
-        float distance_penalty_b = 0.0f;
-        if (setup == spNearest) {
-            distance_penalty_a = 1.0f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f);
-            distance_penalty_b = 1.0f - gauss((b.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f);
-        }
-
-        // the penalites are kept close to range [0-1.x] however, it should not be relied upon
-        float penalty_a = a.overhang + a.visibility +
-                angle_importance * compute_angle_penalty(a.local_ccw_angle)
-                + distance_penalty_a;
-        float penalty_b = b.overhang + b.visibility +
-                angle_importance * compute_angle_penalty(b.local_ccw_angle)
-                + distance_penalty_b;
-
-        return penalty_a < penalty_b;
-    }
-
-    // Comparator used during alignment. If there is close potential aligned point, it is compared to the current
-    // seam point of the perimeter, to find out if the aligned point is not much worse than the current seam
-    // Also used by the random seam generator.
-    bool is_first_not_much_worse(const SeamCandidate &a, const SeamCandidate &b) const {
-        // Blockers/Enforcers discrimination, top priority
-        if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) {
-            // Prefer centers of enforcers.
-            return a.central_enforcer;
-        }
-
-        if (a.type == EnforcedBlockedSeamPoint::Enforced) {
-            return true;
-        }
-
-        if (a.type == EnforcedBlockedSeamPoint::Blocked) {
-            return false;
-        }
-
-        if (a.type != b.type) {
-            return a.type > b.type;
-        }
-
-        //avoid overhangs
-        if ((a.overhang > 0.0f || b.overhang > 0.0f)
-                && abs(a.overhang - b.overhang) > (0.1f * a.perimeter.flow_width)) {
-            return a.overhang < b.overhang;
-        }
-
-        // prefer hidden points (more than 0.5 mm inside)
-        if (a.embedded_distance < -0.5f && b.embedded_distance > -0.5f) {
-            return true;
-        }
-        if (b.embedded_distance < -0.5f && a.embedded_distance > -0.5f) {
-            return false;
-        }
-
-        if (setup == SeamPosition::spRandom) {
-            return true;
-        }
-
-        if (setup == SeamPosition::spRear) {
-            return a.position.y() + SeamPlacer::seam_align_score_tolerance * 5.0f > b.position.y();
-        }
-
-        float penalty_a = a.overhang + a.visibility
-                + angle_importance * compute_angle_penalty(a.local_ccw_angle);
-        float penalty_b = b.overhang + b.visibility +
-                angle_importance * compute_angle_penalty(b.local_ccw_angle);
-
-        return penalty_a <= penalty_b || penalty_a - penalty_b < SeamPlacer::seam_align_score_tolerance;
-    }
-
-    bool are_similar(const SeamCandidate &a, const SeamCandidate &b) const {
-        return is_first_not_much_worse(a, b) && is_first_not_much_worse(b, a);
-    }
-};
-
-#ifdef DEBUG_FILES
-void debug_export_points(const std::vector<PrintObjectSeamData::LayerSeams> &layers,
-        const BoundingBox &bounding_box, const SeamComparator &comparator) {
-    for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) {
-        std::string angles_file_name = debug_out_path(
-                ("angles_" + std::to_string(layer_idx) + ".svg").c_str());
-        SVG angles_svg { angles_file_name, bounding_box };
-        float min_vis = 0;
-        float max_vis = min_vis;
-
-        float min_weight = std::numeric_limits<float>::min();
-        float max_weight = min_weight;
-
-        for (const SeamCandidate &point : layers[layer_idx].points) {
-            Vec3i color = value_to_rgbi(-PI, PI, point.local_ccw_angle);
-            std::string fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + ","
-                    + std::to_string(color.z()) + ")";
-            angles_svg.draw(scaled(Vec2f(point.position.head<2>())), fill);
-            min_vis = std::min(min_vis, point.visibility);
-            max_vis = std::max(max_vis, point.visibility);
-
-            min_weight = std::min(min_weight, -compute_angle_penalty(point.local_ccw_angle));
-            max_weight = std::max(max_weight, -compute_angle_penalty(point.local_ccw_angle));
-
-        }
-
-        std::string visiblity_file_name = debug_out_path(
-                ("visibility_" + std::to_string(layer_idx) + ".svg").c_str());
-        SVG visibility_svg { visiblity_file_name, bounding_box };
-        std::string weights_file_name = debug_out_path(
-                ("weight_" + std::to_string(layer_idx) + ".svg").c_str());
-        SVG weight_svg { weights_file_name, bounding_box };
-        std::string overhangs_file_name = debug_out_path(
-                ("overhang_" + std::to_string(layer_idx) + ".svg").c_str());
-        SVG overhangs_svg { overhangs_file_name, bounding_box };
-
-        for (const SeamCandidate &point : layers[layer_idx].points) {
-            Vec3i color = value_to_rgbi(min_vis, max_vis, point.visibility);
-            std::string visibility_fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + ","
-                    + std::to_string(color.z()) + ")";
-            visibility_svg.draw(scaled(Vec2f(point.position.head<2>())), visibility_fill);
-
-            Vec3i weight_color = value_to_rgbi(min_weight, max_weight,
-                    -compute_angle_penalty(point.local_ccw_angle));
-            std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y())
-                    + ","
-                    + std::to_string(weight_color.z()) + ")";
-            weight_svg.draw(scaled(Vec2f(point.position.head<2>())), weight_fill);
-
-            Vec3i overhang_color = value_to_rgbi(-0.5, 0.5, std::clamp(point.overhang, -0.5f, 0.5f));
-            std::string overhang_fill = "rgb(" + std::to_string(overhang_color.x()) + ","
-                    + std::to_string(overhang_color.y())
-                    + ","
-                    + std::to_string(overhang_color.z()) + ")";
-            overhangs_svg.draw(scaled(Vec2f(point.position.head<2>())), overhang_fill);
-        }
-    }
-}
-#endif
-
-// Pick best seam point based on the given comparator
-void pick_seam_point(std::vector<SeamCandidate> &perimeter_points, size_t start_index,
-        const SeamComparator &comparator) {
-    size_t end_index = perimeter_points[start_index].perimeter.end_index;
-
-    size_t seam_index = start_index;
-    for (size_t index = start_index; index < end_index; ++index) {
-        if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index])) {
-            seam_index = index;
-        }
-    }
-    perimeter_points[start_index].perimeter.seam_index = seam_index;
-}
-
-size_t pick_nearest_seam_point_index(const std::vector<SeamCandidate> &perimeter_points, size_t start_index,
-        const Vec2f &preffered_location) {
-    size_t end_index = perimeter_points[start_index].perimeter.end_index;
-    SeamComparator comparator { spNearest };
-
-    size_t seam_index = start_index;
-    for (size_t index = start_index; index < end_index; ++index) {
-        if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index], preffered_location)) {
-            seam_index = index;
-        }
-    }
-    return seam_index;
-}
-
-// picks random seam point uniformly, respecting enforcers blockers and overhang avoidance.
-void pick_random_seam_point(const std::vector<SeamCandidate> &perimeter_points, size_t start_index) {
-    SeamComparator comparator { spRandom };
-
-    // algorithm keeps a list of viable points and their lengths. If it finds a point
-    // that is much better than the viable_example_index (e.g. better type, no overhang; see is_first_not_much_worse)
-    // then it throws away stored lists and starts from start
-    // in the end, the list should contain points with same type (Enforced > Neutral > Blocked) and also only those which are not
-    // big overhang.
-    size_t viable_example_index = start_index;
-    size_t end_index = perimeter_points[start_index].perimeter.end_index;
-    struct Viable {
-        // Candidate seam point index.
-        size_t index;
-        float edge_length;
-        Vec3f edge;
-    };
-    std::vector<Viable> viables;
-
-    const Vec3f pseudornd_seed = perimeter_points[viable_example_index].position;
-    float rand = std::abs(sin(pseudornd_seed.dot(Vec3f(12.9898f,78.233f, 133.3333f))) * 43758.5453f);
-    rand = rand - (int) rand;
-
-    for (size_t index = start_index; index < end_index; ++index) {
-        if (comparator.are_similar(perimeter_points[index], perimeter_points[viable_example_index])) {
-            // index ok, push info into viables
-            Vec3f edge_to_next { perimeter_points[index == end_index - 1 ? start_index : index + 1].position
-                    - perimeter_points[index].position };
-            float dist_to_next = edge_to_next.norm();
-            viables.push_back( { index, dist_to_next, edge_to_next });
-        } else if (comparator.is_first_not_much_worse(perimeter_points[viable_example_index],
-                perimeter_points[index])) {
-            // index is worse then viable_example_index, skip this point
-        } else {
-            // index is better than viable example index, update example, clear gathered info, start again
-            // clear up all gathered info, start from scratch, update example index
-            viable_example_index = index;
-            viables.clear();
-
-            Vec3f edge_to_next = (perimeter_points[index == end_index - 1 ? start_index : index + 1].position
-                    - perimeter_points[index].position);
-            float dist_to_next = edge_to_next.norm();
-            viables.push_back( { index, dist_to_next, edge_to_next });
-        }
-    }
-
-    // now pick random point from the stored options
-    float len_sum = std::accumulate(viables.begin(), viables.end(), 0.0f, [](const float acc, const Viable &v) {
-        return acc + v.edge_length;
-    });
-    float picked_len = len_sum * rand;
-
-    size_t point_idx = 0;
-    while (picked_len - viables[point_idx].edge_length > 0) {
-        picked_len = picked_len - viables[point_idx].edge_length;
-        point_idx++;
-    }
-
-    Perimeter &perimeter = perimeter_points[start_index].perimeter;
-    perimeter.seam_index = viables[point_idx].index;
-    perimeter.final_seam_position = perimeter_points[perimeter.seam_index].position
-            + viables[point_idx].edge.normalized() * picked_len;
-    perimeter.finalized = true;
-}
-
-} // namespace SeamPlacerImpl
-
-// Parallel process and extract each perimeter polygon of the given print object.
-// Gather SeamCandidates of each layer into vector and build KDtree over them
-// Store results in the SeamPlacer variables m_seam_per_object
-void SeamPlacer::gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info) {
-    using namespace SeamPlacerImpl;
-    PrintObjectSeamData &seam_data = m_seam_per_object.emplace(po, PrintObjectSeamData { }).first->second;
-    seam_data.layers.resize(po->layer_count());
-
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, po->layers().size()),
-            [po, &global_model_info, &seam_data]
-            (tbb::blocked_range<size_t> r) {
-                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
-                    PrintObjectSeamData::LayerSeams &layer_seams = seam_data.layers[layer_idx];
-                    const Layer *layer = po->get_layer(layer_idx);
-                    auto unscaled_z = layer->slice_z;
-                    std::vector<const LayerRegion*> regions;
-                    //NOTE corresponding region ptr may be null, if the layer has zero perimeters
-                    Polygons polygons = extract_perimeter_polygons(layer, regions);
-                    for (size_t poly_index = 0; poly_index < polygons.size(); ++poly_index) {
-                        process_perimeter_polygon(polygons[poly_index], unscaled_z,
-                                regions[poly_index], global_model_info, layer_seams);
-                    }
-                    auto functor = SeamCandidateCoordinateFunctor { layer_seams.points };
-                    seam_data.layers[layer_idx].points_tree =
-                            std::make_unique<PrintObjectSeamData::SeamCandidatesTree>(functor,
-                                    layer_seams.points.size());
-                }
-            }
-    );
-}
-
-void SeamPlacer::calculate_candidates_visibility(const PrintObject *po,
-        const SeamPlacerImpl::GlobalModelInfo &global_model_info) {
-    using namespace SeamPlacerImpl;
-
-    std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
-            [&layers, &global_model_info](tbb::blocked_range<size_t> r) {
-                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
-                    for (auto &perimeter_point : layers[layer_idx].points) {
-                        perimeter_point.visibility = global_model_info.calculate_point_visibility(
-                                perimeter_point.position);
-                    }
-                }
-            });
-}
-
-void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po) {
-    using namespace SeamPlacerImpl;
-    using PerimeterDistancer = AABBTreeLines::LinesDistancer<Linef>;
-
-    std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
-            [po, &layers](tbb::blocked_range<size_t> r) {
-                std::unique_ptr<PerimeterDistancer> prev_layer_distancer;
-                if (r.begin() > 0) { // previous layer exists
-                    prev_layer_distancer = std::make_unique<PerimeterDistancer>(to_unscaled_linesf(po->layers()[r.begin() - 1]->lslices));
-                }
-
-                for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
-                    size_t regions_with_perimeter = 0;
-                    for (const LayerRegion *region : po->layers()[layer_idx]->regions()) {
-                        if (region->perimeters().size() > 0) {
-                            regions_with_perimeter++;
-                        }
-                    };
-                    bool should_compute_layer_embedding = regions_with_perimeter > 1;
-                    std::unique_ptr<PerimeterDistancer> current_layer_distancer        = std::make_unique<PerimeterDistancer>(
-                        to_unscaled_linesf(po->layers()[layer_idx]->lslices));
-
-                    for (SeamCandidate &perimeter_point : layers[layer_idx].points) {
-                        Vec2f point = Vec2f { perimeter_point.position.head<2>() };
-                        if (prev_layer_distancer.get() != nullptr) {
-                            perimeter_point.overhang = prev_layer_distancer->distance_from_lines<true>(point.cast<double>())
-                                    + 0.6f * perimeter_point.perimeter.flow_width
-                                    - tan(SeamPlacer::overhang_angle_threshold)
-                                            * po->layers()[layer_idx]->height;
-                            perimeter_point.overhang =
-                                    perimeter_point.overhang < 0.0f ? 0.0f : perimeter_point.overhang;
-                        }
-
-                        if (should_compute_layer_embedding) { // search for embedded perimeter points (points hidden inside the print ,e.g. multimaterial join, best position for seam)
-                            perimeter_point.embedded_distance = current_layer_distancer->distance_from_lines<true>(point.cast<double>())
-                                    + 0.6f * perimeter_point.perimeter.flow_width;
-                        }
-                    }
-
-                    prev_layer_distancer.swap(current_layer_distancer);
-                }
-            }
+ObjectSeams precalculate_seams(
+    const Params &params,
+    ObjectShells &&seam_data,
+    const std::function<void(void)> &throw_if_canceled
+) {
+    ObjectSeams result;
+
+    for (auto &[print_object, shells] : seam_data) {
+        switch (params.seam_preference) {
+        case spAligned: {
+            const Transform3d transformation{print_object->trafo_centered()};
+            const ModelVolumePtrs &volumes{print_object->model_object()->volumes};
+
+            Slic3r::ModelInfo::Visibility
+                points_visibility{transformation, volumes, params.visibility, throw_if_canceled};
+            throw_if_canceled();
+            const Aligned::VisibilityCalculator visibility_calculator{
+                points_visibility, params.convex_visibility_modifier,
+                params.concave_visibility_modifier};
+
+            result[print_object] = Aligned::get_object_seams(
+                std::move(shells), visibility_calculator, params.aligned
             );
+            break;
         }
-
-// Estimates, if there is good seam point in the layer_idx which is close to last_point_pos
-// uses comparator.is_first_not_much_worse method to compare current seam with the closest point
-// (if current seam is too far away )
-// If the current chosen stream is close enough, it is stored in seam_string. returns true and updates last_point_pos
-// If the closest point is good enough to replace current chosen seam, it is stored in potential_string_seams, returns true and updates last_point_pos
-// Otherwise does nothing, returns false
-// Used by align_seam_points().
-std::optional<std::pair<size_t, size_t>> SeamPlacer::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 {
-    using namespace SeamPlacerImpl;
-    std::vector<size_t> nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position,
-            max_distance);
-
-    if (nearby_points_indices.empty()) {
-        return {};
-    }
-
-    size_t best_nearby_point_index = nearby_points_indices[0];
-    size_t nearest_point_index = nearby_points_indices[0];
-
-    // Now find best nearby point, nearest point, and corresponding indices
-    for (const size_t &nearby_point_index : nearby_points_indices) {
-        const SeamCandidate &point = layers[layer_idx].points[nearby_point_index];
-        if (point.perimeter.finalized) {
-            continue; // skip over finalized perimeters, try to find some that is not finalized
+        case spRear: {
+            result[print_object] = Rear::get_object_seams(std::move(shells), params.rear_project_threshold);
+            break;
         }
-        if (comparator.is_first_better(point, layers[layer_idx].points[best_nearby_point_index],
-                projected_position.head<2>())
-                || layers[layer_idx].points[best_nearby_point_index].perimeter.finalized) {
-            best_nearby_point_index = nearby_point_index;
+        case spRandom: {
+            result[print_object] = Random::get_object_seams(std::move(shells), params.random_seed);
+            break;
         }
-        if ((point.position - projected_position).squaredNorm()
-                < (layers[layer_idx].points[nearest_point_index].position - projected_position).squaredNorm()
-                || layers[layer_idx].points[nearest_point_index].perimeter.finalized) {
-            nearest_point_index = nearby_point_index;
+        case spNearest: {
+            throw std::runtime_error("Cannot precalculate seams for nearest position!");
         }
+        }
+        throw_if_canceled();
     }
-
-    const SeamCandidate &best_nearby_point = layers[layer_idx].points[best_nearby_point_index];
-    const SeamCandidate &nearest_point = layers[layer_idx].points[nearest_point_index];
-
-    if (nearest_point.perimeter.finalized) {
-        //all points are from already finalized perimeter, skip
-        return {};
-    }
-
-    //from the nearest_point, deduce index of seam in the next layer
-    const SeamCandidate &next_layer_seam = layers[layer_idx].points[nearest_point.perimeter.seam_index];
-
-    // First try to pick central enforcer if any present
-    if (next_layer_seam.central_enforcer
-            && (next_layer_seam.position - projected_position).squaredNorm()
-                    < sqr(3 * max_distance)) {
-        return {std::pair<size_t, size_t> {layer_idx, nearest_point.perimeter.seam_index}};
-    }
-
-    // First try to align the nearest, then try the best nearby
-    if (comparator.is_first_not_much_worse(nearest_point, next_layer_seam)) {
-        return {std::pair<size_t, size_t> {layer_idx, nearest_point_index}};
-    }
-    // If nearest point is not good enough, try it with the best nearby point.
-    if (comparator.is_first_not_much_worse(best_nearby_point, next_layer_seam)) {
-        return {std::pair<size_t, size_t> {layer_idx, best_nearby_point_index}};
-    }
-
-    return {};
+    return result;
 }
 
-std::vector<std::pair<size_t, size_t>> SeamPlacer::find_seam_string(const PrintObject *po,
-        std::pair<size_t, size_t> start_seam, const SeamPlacerImpl::SeamComparator &comparator) const {
-    const std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object.find(po)->second.layers;
-    int layer_idx = start_seam.first;
+Params Placer::get_params(const DynamicPrintConfig &config) {
+    Params params{};
 
-    //initialize searching for seam string - cluster of nearby seams on previous and next layers
-    int next_layer = layer_idx + 1;
-    int step = 1;
-    std::pair<size_t, size_t> prev_point_index = start_seam;
-    std::vector<std::pair<size_t, size_t>> seam_string { start_seam };
-
-    auto reverse_lookup_direction = [&]() {
-        step = -1;
-        prev_point_index = start_seam;
-        next_layer = layer_idx - 1;
-    };
-
-    while (next_layer >= 0) {
-        if (next_layer >= int(layers.size())) {
-            reverse_lookup_direction();
-            if (next_layer < 0) {
-                break;
-            }
-        }
-        float max_distance = SeamPlacer::seam_align_tolerable_dist_factor *
-                layers[start_seam.first].points[start_seam.second].perimeter.flow_width;
-        Vec3f prev_position = layers[prev_point_index.first].points[prev_point_index.second].position;
-        Vec3f projected_position = prev_position;
-        projected_position.z() = float(po->get_layer(next_layer)->slice_z);
-
-        std::optional<std::pair<size_t, size_t>> maybe_next_seam = find_next_seam_in_layer(layers, projected_position,
-                next_layer,
-                max_distance, comparator);
-
-        if (maybe_next_seam.has_value()) {
-            // For old macOS (pre 10.14), std::optional does not have .value() method, so the code is using operator*() instead.
-            seam_string.push_back(maybe_next_seam.operator*());
-            prev_point_index = seam_string.back();
-            //String added, prev_point_index updated
-        } else {
-            if (step == 1) {
-                reverse_lookup_direction();
-                if (next_layer < 0) {
-                    break;
-                }
-            } else {
-                break;
-            }
-        }
-        next_layer += step;
+    params.perimeter.elephant_foot_compensation = config.opt_float("elefant_foot_compensation");
+    if (config.opt_int("raft_layers") > 0) {
+        params.perimeter.elephant_foot_compensation = 0.0;
     }
-    return seam_string;
+    params.random_seed = 1653710332u;
+
+    params.aligned.max_detour = 1.0;
+    params.convex_visibility_modifier = 1.1;
+    params.concave_visibility_modifier = 0.9;
+    params.perimeter.overhang_threshold = Slic3r::Geometry::deg2rad(55.0);
+    params.perimeter.convex_threshold = Slic3r::Geometry::deg2rad(10.0);
+    params.perimeter.concave_threshold = Slic3r::Geometry::deg2rad(15.0);
+
+    params.seam_preference = config.opt_enum<SeamPosition>("seam_position");
+    params.staggered_inner_seams = config.opt_bool("staggered_inner_seams");
+
+    params.max_nearest_detour = 1.0;
+    params.rear_project_threshold = 0.05; // %
+    params.aligned.jump_visibility_threshold = 0.6;
+    params.max_distance = 5.0;
+    params.perimeter.oversampling_max_distance = 0.2;
+    params.perimeter.embedding_threshold = 0.5;
+    params.perimeter.painting_radius = 0.1;
+    params.perimeter.simplification_epsilon = 0.001;
+    params.perimeter.smooth_angle_arm_length = 0.5;
+    params.perimeter.sharp_angle_arm_length = 0.05;
+
+    params.visibility.raycasting_visibility_samples_count = 30000;
+    params.visibility.fast_decimation_triangle_count_target = 16000;
+    params.visibility.sqr_rays_per_sample_point = 5;
+
+    return params;
 }
 
-// clusters already chosen seam points into strings across multiple layers, and then
-// aligns the strings via polynomial fit
-// Does not change the positions of the SeamCandidates themselves, instead stores
-// the new aligned position into the shared Perimeter structure of each perimeter
-// Note that this position does not necesarilly lay on the perimeter.
-void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) {
-    using namespace SeamPlacerImpl;
+ObjectLayerPerimeters sort_to_layers(ObjectShells &&object_shells) {
+    ObjectLayerPerimeters result;
+    for (auto &[print_object, shells] : object_shells) {
+        const std::size_t layer_count{print_object->layer_count()};
+        result[print_object] = LayerPerimeters(layer_count);
 
-    // Prepares Debug files for writing.
-#ifdef DEBUG_FILES
-    Slic3r::CNumericLocalesSetter locales_setter;
-    auto clusters_f = debug_out_path("seam_clusters.obj");
-    FILE *clusters = boost::nowide::fopen(clusters_f.c_str(), "w");
-    if (clusters == nullptr) {
-        BOOST_LOG_TRIVIAL(error)
-        << "stl_write_obj: Couldn't open " << clusters_f << " for writing";
-        return;
-    }
-    auto aligned_f = debug_out_path("aligned_clusters.obj");
-    FILE *aligns = boost::nowide::fopen(aligned_f.c_str(), "w");
-    if (aligns == nullptr) {
-        BOOST_LOG_TRIVIAL(error)
-        << "stl_write_obj: Couldn't open " << clusters_f << " for writing";
-        return;
-    }
-#endif
-
-    //gather vector of all seams on the print_object - pair of layer_index and seam__index within that layer
-    const std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
-    std::vector<std::pair<size_t, size_t>> seams;
-    for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) {
-        const std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
-        size_t current_point_index = 0;
-        while (current_point_index < layer_perimeter_points.size()) {
-            seams.emplace_back(layer_idx, layer_perimeter_points[current_point_index].perimeter.seam_index);
-            current_point_index = layer_perimeter_points[current_point_index].perimeter.end_index;
+        for (Shells::Shell<> &shell : shells) {
+            for (Shells::Slice<> &slice : shell) {
+                const BoundingBox bounding_box{Geometry::scaled(slice.boundary.positions)};
+                result[print_object][slice.layer_index].push_back(
+                    BoundedPerimeter{std::move(slice.boundary), bounding_box}
+                );
+            }
         }
     }
+    return result;
+}
 
-    //sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice
-    std::stable_sort(seams.begin(), seams.end(),
-            [&comparator, &layers](const std::pair<size_t, size_t> &left,
-                    const std::pair<size_t, size_t> &right) {
-                return comparator.is_first_better(layers[left.first].points[left.second],
-                        layers[right.first].points[right.second]);
-            }
+void Placer::init(
+    SpanOfConstPtrs<PrintObject> objects,
+    const Params &params,
+    const std::function<void(void)> &throw_if_canceled
+) {
+    BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: init: start";
+
+    ObjectPainting object_painting;
+    for (const PrintObject *print_object : objects) {
+        const Transform3d transformation{print_object->trafo_centered()};
+        const ModelVolumePtrs &volumes{print_object->model_object()->volumes};
+        object_painting.emplace(print_object, ModelInfo::Painting{transformation, volumes});
+    }
+
+    ObjectShells seam_data{partition_to_shells(objects, params, object_painting, throw_if_canceled)};
+    this->params = params;
+
+    if (this->params.seam_preference != spNearest) {
+        this->seams_per_object =
+            precalculate_seams(params, std::move(seam_data), throw_if_canceled);
+    } else {
+        this->perimeters_per_layer = sort_to_layers(std::move(seam_data));
+    }
+
+    BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: init: end";
+}
+
+const SeamPerimeterChoice &choose_closest_seam(
+    const std::vector<SeamPerimeterChoice> &seams, const Polygon &loop_polygon
+) {
+    BoundingBoxes choose_from;
+    choose_from.reserve(seams.size());
+    for (const SeamPerimeterChoice &choice : seams) {
+        choose_from.push_back(choice.bounding_box);
+    }
+
+    const std::size_t choice_index{
+        Geometry::pick_closest_bounding_box(loop_polygon.bounding_box(), choose_from).first};
+
+    return seams[choice_index];
+}
+
+std::pair<std::size_t, Vec2d> project_to_extrusion_loop(
+    const SeamChoice &seam_choice, const Perimeters::Perimeter &perimeter, const Linesf &loop_lines
+) {
+    const AABBTreeLines::LinesDistancer<Linef> distancer{loop_lines};
+
+    const bool is_at_vertex{seam_choice.previous_index == seam_choice.next_index};
+    const Vec2d edge{
+        perimeter.positions[seam_choice.next_index] -
+        perimeter.positions[seam_choice.previous_index]};
+    const Vec2d normal{
+        is_at_vertex ?
+            Geometry::get_polygon_normal(perimeter.positions, seam_choice.previous_index, 0.1) :
+            Geometry::get_normal(edge)};
+
+    double depth{distancer.distance_from_lines<false>(seam_choice.position)};
+    const Vec2d final_position{seam_choice.position - normal * depth};
+
+    auto [_, loop_line_index, loop_point] = distancer.distance_from_lines_extra<false>(final_position
     );
+    return {loop_line_index, loop_point};
+}
 
-    //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment
-    // Keeping the vectors outside, so with a bit of luck they will not get reallocated after couple of for loop iterations.
-    std::vector<std::pair<size_t, size_t>> seam_string;
-    std::vector<std::pair<size_t, size_t>> alternative_seam_string;
-    std::vector<Vec2f> observations;
-    std::vector<float> observation_points;
-    std::vector<float> weights;
+std::optional<Vec2d> offset_along_loop_lines(
+    const Vec2d &point,
+    const std::size_t loop_line_index,
+    const Linesf &loop_lines,
+    const double offset
+) {
+    double distance{0};
+    Vec2d previous_point{point};
+    std::optional<Vec2d> offset_point;
+    Geometry::visit_forward(loop_line_index, loop_lines.size(), [&](std::size_t index) {
+        const Vec2d next_point{loop_lines[index].b};
+        const Vec2d edge{next_point - previous_point};
 
-    int global_index = 0;
-    while (global_index < int(seams.size())) {
-        size_t layer_idx = seams[global_index].first;
-        size_t seam_index = seams[global_index].second;
-        global_index++;
-        const std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
-        if (layer_perimeter_points[seam_index].perimeter.finalized) {
-            // This perimeter is already aligned, skip seam
-            continue;
-        } else {
-            seam_string = this->find_seam_string(po, { layer_idx, seam_index }, comparator);
-            size_t step_size = 1 + seam_string.size() / 20;
-            for (size_t alternative_start = 0; alternative_start < seam_string.size(); alternative_start += step_size) {
-                size_t start_layer_idx = seam_string[alternative_start].first;
-                size_t seam_idx =
-                        layers[start_layer_idx].points[seam_string[alternative_start].second].perimeter.seam_index;
-                alternative_seam_string = this->find_seam_string(po,
-                        std::pair<size_t, size_t>(start_layer_idx, seam_idx), comparator);
-                if (alternative_seam_string.size() > seam_string.size()) {
-                    seam_string = std::move(alternative_seam_string);
-                }
-            }
-            if (seam_string.size() < seam_align_minimum_string_seams) {
-                //string NOT long enough to be worth aligning, skip
-                continue;
-            }
+        if (distance + edge.norm() > offset) {
+            const double remaining_distance{offset - distance};
+            offset_point = previous_point + remaining_distance * edge.normalized();
+            return true;
+        }
 
-            // String is long enough, all string seams and potential string seams gathered, now do the alignment
-            //sort by layer index
-            std::sort(seam_string.begin(), seam_string.end(),
-                    [](const std::pair<size_t, size_t> &left, const std::pair<size_t, size_t> &right) {
-                        return left.first < right.first;
-                    });
+        distance += edge.norm();
+        previous_point = next_point;
 
-            //repeat the alignment for the current seam, since it could be skipped due to alternative path being aligned.
-            global_index--;
+        return false;
+    });
 
-            // gather all positions of seams and their weights
-            observations.resize(seam_string.size());
-            observation_points.resize(seam_string.size());
-            weights.resize(seam_string.size());
+    return offset_point;
+}
 
-            auto angle_3d = [](const Vec3f& a, const Vec3f& b){
-                return std::abs(acosf(a.normalized().dot(b.normalized())));
-            };
+double get_angle(const SeamChoice &seam_choice, const Perimeters::Perimeter &perimeter) {
+    const bool is_at_vertex{seam_choice.previous_index == seam_choice.next_index};
+    return is_at_vertex ? perimeter.angles[seam_choice.previous_index] : 0.0;
+}
 
-            auto angle_weight = [](float angle){
-                return 1.0f / (0.1f + compute_angle_penalty(angle));
-            };
+Point finalize_seam_position(
+    const Polygon &loop_polygon,
+    const SeamChoice &seam_choice,
+    const Perimeters::Perimeter &perimeter,
+    const double loop_width,
+    const bool do_staggering
+) {
+    const Linesf loop_lines{to_unscaled_linesf({ExPolygon{loop_polygon}})};
+    const auto [loop_line_index, loop_point]{
+        project_to_extrusion_loop(seam_choice, perimeter, loop_lines)};
 
-            //gather points positions and weights
-            float total_length = 0.0f;
-            Vec3f last_point_pos = layers[seam_string[0].first].points[seam_string[0].second].position;
-            for (size_t index = 0; index < seam_string.size(); ++index) {
-                const SeamCandidate &current = layers[seam_string[index].first].points[seam_string[index].second];
-                float layer_angle = 0.0f;
-                if (index > 0 && index < seam_string.size() - 1) {
-                    layer_angle = angle_3d(
-                                    current.position
-                                            - layers[seam_string[index - 1].first].points[seam_string[index - 1].second].position,
-                                    layers[seam_string[index + 1].first].points[seam_string[index + 1].second].position
-                                            - current.position
-                                            );
-                }
-                observations[index] = current.position.head<2>();
-                observation_points[index] = current.position.z();
-                weights[index] = angle_weight(current.local_ccw_angle);
-                float curling_influence = layer_angle > 2.0 * std::abs(current.local_ccw_angle) ? -0.8f : 1.0f;
-                if (current.type == EnforcedBlockedSeamPoint::Enforced) {
-                    curling_influence = 1.0f;
-                    weights[index] += 3.0f;
-                }
-                total_length += curling_influence * (last_point_pos - current.position).norm();
-                last_point_pos = current.position;
-            }
+    // ExtrusionRole::Perimeter is inner perimeter.
+    if (do_staggering) {
+        const double depth = (loop_point - seam_choice.position).norm() -
+            loop_width / 2.0;
+        const double angle{get_angle(seam_choice, perimeter)};
+        const double initial_offset{angle > 0 ? angle / 2.0 * depth : 0.0};
+        const double additional_offset{angle < 0 ? std::cos(angle / 2.0) * depth : depth};
 
-            if (comparator.setup == spRear) {
-                total_length *= 0.3f;
-            }
+        const double staggering_offset{initial_offset + additional_offset};
 
-            // Curve Fitting
-            size_t number_of_segments = std::max(size_t(1),
-                    size_t(std::max(0.0f,total_length) / SeamPlacer::seam_align_mm_per_segment));
-            auto curve = Geometry::fit_cubic_bspline(observations, observation_points, weights, number_of_segments);
+        std::optional<Vec2d> staggered_point{
+            offset_along_loop_lines(loop_point, loop_line_index, loop_lines, staggering_offset)};
 
-            // Do alignment - compute fitted point for each point in the string from its Z coord, and store the position into
-            // Perimeter structure of the point; also set flag aligned to true
-            for (size_t index = 0; index < seam_string.size(); ++index) {
-                const auto &pair = seam_string[index];
-                float t = std::min(1.0f, std::pow(std::abs(layers[pair.first].points[pair.second].local_ccw_angle)
-                        / SeamPlacer::sharp_angle_snapping_threshold, 3.0f));
-                if (layers[pair.first].points[pair.second].type == EnforcedBlockedSeamPoint::Enforced){
-                    t = std::max(0.4f, t);
-                }
-
-                Vec3f current_pos = layers[pair.first].points[pair.second].position;
-                Vec2f fitted_pos = curve.get_fitted_value(current_pos.z());
-
-                //interpolate between current and fitted position, prefer current pos for large weights.
-                Vec3f final_position = t * current_pos + (1.0f - t) * to_3d(fitted_pos, current_pos.z());
-
-                Perimeter &perimeter = layers[pair.first].points[pair.second].perimeter;
-                perimeter.seam_index = pair.second;
-                perimeter.final_seam_position = final_position;
-                perimeter.finalized = true;
-            }
-
-#ifdef DEBUG_FILES
-            auto randf = []() {
-                return float(rand()) / float(RAND_MAX);
-            };
-            Vec3f color { randf(), randf(), randf() };
-            for (size_t i = 0; i < seam_string.size(); ++i) {
-                auto orig_seam = layers[seam_string[i].first].points[seam_string[i].second];
-                fprintf(clusters, "v %f %f %f %f %f %f \n", orig_seam.position[0],
-                        orig_seam.position[1],
-                        orig_seam.position[2], color[0], color[1],
-                        color[2]);
-            }
-
-            color = Vec3f { randf(), randf(), randf() };
-            for (size_t i = 0; i < seam_string.size(); ++i) {
-                const Perimeter &perimeter = layers[seam_string[i].first].points[seam_string[i].second].perimeter;
-                fprintf(aligns, "v %f %f %f %f %f %f \n", perimeter.final_seam_position[0],
-                        perimeter.final_seam_position[1],
-                        perimeter.final_seam_position[2], color[0], color[1],
-                        color[2]);
-            }
-#endif
+        if (staggered_point) {
+            return scaled(*staggered_point);
         }
     }
 
-#ifdef DEBUG_FILES
-    fclose(clusters);
-    fclose(aligns);
-#endif
-
+    return scaled(loop_point);
 }
 
-void SeamPlacer::init(const Print &print, std::function<void(void)> throw_if_canceled_func) {
-    using namespace SeamPlacerImpl;
-    m_seam_per_object.clear();
-
-    for (const PrintObject *po : print.objects()) {
-        throw_if_canceled_func();
-        SeamPosition configured_seam_preference = po->config().seam_position.value;
-        SeamComparator comparator { configured_seam_preference };
-
-        {
-            GlobalModelInfo global_model_info { };
-            gather_enforcers_blockers(global_model_info, po);
-            throw_if_canceled_func();
-            if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) {
-                compute_global_occlusion(global_model_info, po, throw_if_canceled_func);
-            }
-            throw_if_canceled_func();
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: gather_seam_candidates: start";
-            gather_seam_candidates(po, global_model_info);
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: gather_seam_candidates: end";
-            throw_if_canceled_func();
-            if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) {
-                BOOST_LOG_TRIVIAL(debug)
-                << "SeamPlacer: calculate_candidates_visibility : start";
-                calculate_candidates_visibility(po, global_model_info);
-                BOOST_LOG_TRIVIAL(debug)
-                << "SeamPlacer: calculate_candidates_visibility : end";
-            }
-        } // destruction of global_model_info (large structure, no longer needed)
-        throw_if_canceled_func();
-        BOOST_LOG_TRIVIAL(debug)
-        << "SeamPlacer: calculate_overhangs and layer embdedding : start";
-        calculate_overhangs_and_layer_embedding(po);
-        BOOST_LOG_TRIVIAL(debug)
-        << "SeamPlacer: calculate_overhangs and layer embdedding: end";
-        throw_if_canceled_func();
-        if (configured_seam_preference != spNearest) { // For spNearest, the seam is picked in the place_seam method with actual nozzle position information
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: pick_seam_point : start";
-            //pick seam point
-            std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object[po].layers;
-            tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size()),
-                    [&layers, configured_seam_preference, comparator](tbb::blocked_range<size_t> r) {
-                        for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
-                            std::vector<SeamCandidate> &layer_perimeter_points = layers[layer_idx].points;
-                            for (size_t current = 0; current < layer_perimeter_points.size();
-                                    current = layer_perimeter_points[current].perimeter.end_index)
-                                if (configured_seam_preference == spRandom)
-                                    pick_random_seam_point(layer_perimeter_points, current);
-                                else
-                                    pick_seam_point(layer_perimeter_points, current, comparator);
-                        }
-                    });
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: pick_seam_point : end";
-        }
-        throw_if_canceled_func();
-        if (configured_seam_preference == spAligned || configured_seam_preference == spRear) {
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: align_seam_points : start";
-            align_seam_points(po, comparator);
-            BOOST_LOG_TRIVIAL(debug)
-            << "SeamPlacer: align_seam_points : end";
-        }
-
-#ifdef DEBUG_FILES
-        debug_export_points(m_seam_per_object[po].layers, po->bounding_box(), comparator);
-#endif
+std::pair<SeamChoice, std::size_t> place_seam_near(
+    const std::vector<BoundedPerimeter> &layer_perimeters,
+    const ExtrusionLoop &loop,
+    const Point &position,
+    const double max_detour
+) {
+    BoundingBoxes choose_from;
+    choose_from.reserve(layer_perimeters.size());
+    for (const BoundedPerimeter &perimeter : layer_perimeters) {
+        choose_from.push_back(perimeter.bounding_box);
     }
+
+    const Polygon loop_polygon{Geometry::to_polygon(loop)};
+
+    const std::size_t choice_index{
+        Geometry::pick_closest_bounding_box(loop_polygon.bounding_box(), choose_from).first};
+
+    Seams::Aligned::Impl::Nearest nearest{unscaled(position), max_detour};
+
+    const SeamChoice choice{Seams::choose_seam_point(layer_perimeters[choice_index].perimeter, nearest)};
+
+    return {choice, choice_index};
 }
 
-Point SeamPlacer::place_seam(const Layer *layer, const ExtrusionLoop &loop, bool external_first,
-        const Point &last_pos) const {
-    using namespace SeamPlacerImpl;
+Point Placer::place_seam(const Layer *layer, const ExtrusionLoop &loop, const Point &last_pos) const {
     const PrintObject *po = layer->object();
     // Must not be called with supprot layer.
-    assert(dynamic_cast<const SupportLayer*>(layer) == nullptr);
+    assert(dynamic_cast<const SupportLayer *>(layer) == nullptr);
     // Object layer IDs are incremented by the number of raft layers.
     assert(layer->id() >= po->slicing_parameters().raft_layers());
     const size_t layer_index = layer->id() - po->slicing_parameters().raft_layers();
-    const double unscaled_z = layer->slice_z;
 
-    auto get_next_loop_point = [loop](ExtrusionLoop::ClosestPathPoint current) {
-        current.segment_idx += 1;
-        if (current.segment_idx >= loop.paths[current.path_idx].polyline.points.size()) {
-            current.path_idx = next_idx_modulo(current.path_idx, loop.paths.size());
-            current.segment_idx = 0;
-        }
-        current.foot_pt = loop.paths[current.path_idx].polyline.points[current.segment_idx];
-        return current;
-    };
+    const Polygon loop_polygon{Geometry::to_polygon(loop)};
 
-    const PrintObjectSeamData::LayerSeams &layer_perimeters =
-            m_seam_per_object.find(layer->object())->second.layers[layer_index];
+    const bool do_staggering{this->params.staggered_inner_seams && loop.role() == ExtrusionRole::Perimeter};
+    const double loop_width{loop.paths.empty() ? 0.0 : loop.paths.front().width()};
 
-    // Find the closest perimeter in the SeamPlacer to this loop.
-    // Repeat search until two consecutive points of the loop are found, that result in the same closest_perimeter
-    // This is beacuse with arachne, T-Junctions may exist and sometimes the wrong perimeter was chosen
-    size_t closest_perimeter_point_index = 0;
-    { // local space for the closest_perimeter_point_index
-        Perimeter *closest_perimeter = nullptr;
-        ExtrusionLoop::ClosestPathPoint closest_point{0,0,loop.paths[0].polyline.points[0]};
-        size_t points_count = std::accumulate(loop.paths.begin(), loop.paths.end(), 0, [](size_t acc,const ExtrusionPath& p) {
-           return acc + p.polyline.points.size();
-        });
-        for (size_t i = 0; i < points_count; ++i) {
-            Vec2f unscaled_p = unscaled<float>(closest_point.foot_pt);
-            closest_perimeter_point_index = find_closest_point(*layer_perimeters.points_tree.get(),
-                    to_3d(unscaled_p, float(unscaled_z)));
-            if (closest_perimeter != &layer_perimeters.points[closest_perimeter_point_index].perimeter) {
-                closest_perimeter = &layer_perimeters.points[closest_perimeter_point_index].perimeter;
-                closest_point = get_next_loop_point(closest_point);
-            } else {
-                break;
-            }
-        }
-    }
-
-    Vec3f seam_position;
-    size_t seam_index;
-    if (const Perimeter &perimeter = layer_perimeters.points[closest_perimeter_point_index].perimeter;
-    perimeter.finalized) {
-        seam_position = perimeter.final_seam_position;
-        seam_index = perimeter.seam_index;
+    if (this->params.seam_preference == spNearest) {
+        const std::vector<BoundedPerimeter> &perimeters{this->perimeters_per_layer.at(po)[layer_index]};
+        const auto [seam_choice, perimeter_index] = place_seam_near(perimeters, loop, last_pos, this->params.max_nearest_detour);
+        return finalize_seam_position(loop_polygon, seam_choice, perimeters[perimeter_index].perimeter, loop_width, do_staggering);
     } else {
-        seam_index =
-                po->config().seam_position == spNearest ?
-                        pick_nearest_seam_point_index(layer_perimeters.points, perimeter.start_index,
-                                unscaled<float>(last_pos)) :
-                        perimeter.seam_index;
-        seam_position = layer_perimeters.points[seam_index].position;
+        const SeamPerimeterChoice &seam_perimeter_choice{choose_closest_seam(this->seams_per_object.at(po)[layer_index], loop_polygon)};
+        return finalize_seam_position(loop_polygon, seam_perimeter_choice.choice, seam_perimeter_choice.perimeter, loop_width, do_staggering);
     }
-
-    Point seam_point = Point::new_scale(seam_position.x(), seam_position.y());
-
-    if (loop.role() == ExtrusionRole::Perimeter) { //Hopefully inner perimeter
-        const SeamCandidate &perimeter_point = layer_perimeters.points[seam_index];
-        ExtrusionLoop::ClosestPathPoint projected_point = loop.get_closest_path_and_point(seam_point, false);
-        // determine depth of the seam point.
-        float depth = (float) unscale(Point(seam_point - projected_point.foot_pt)).norm();
-        float beta_angle = cos(perimeter_point.local_ccw_angle / 2.0f);
-        size_t index_of_prev =
-                seam_index == perimeter_point.perimeter.start_index ?
-                                                                      perimeter_point.perimeter.end_index - 1 :
-                                                                      seam_index - 1;
-        size_t index_of_next =
-                seam_index == perimeter_point.perimeter.end_index - 1 ?
-                                                                        perimeter_point.perimeter.start_index :
-                                                                        seam_index + 1;
-
-        if ((seam_position - perimeter_point.position).squaredNorm() < depth && // seam is on perimeter point
-                perimeter_point.local_ccw_angle < -EPSILON // In concave angles
-                        ) { // In this case, we are at internal perimeter, where the external perimeter has seam in concave angle. We want to align
-                // the internal seam into the concave corner, and not on the perpendicular projection on the closest edge (which is what the split_at function does)
-            Vec2f dir_to_middle =
-                    ((perimeter_point.position - layer_perimeters.points[index_of_prev].position).head<2>().normalized()
-                            + (perimeter_point.position - layer_perimeters.points[index_of_next].position).head<2>().normalized())
-                            * 0.5;
-            depth = 1.4142 * depth / beta_angle;
-            // There are some nice geometric identities in determination of the correct depth of new seam point.
-            //overshoot the target depth, in concave angles it will correctly snap to the corner; TODO: find out why such big overshoot is needed.
-            Vec2f final_pos = perimeter_point.position.head<2>() + depth * dir_to_middle;
-            projected_point = loop.get_closest_path_and_point(Point::new_scale(final_pos.x(), final_pos.y()), false);
-        } else { // not concave angle, in that case the nearest point is the good candidate
-            // but for staggering, we also need to recompute depth of the inner perimter, because in convex corners, the distance is larger than layer width
-            // we want the perpendicular depth, not distance to nearest point
-            depth = depth * beta_angle / 1.4142;
-        }
-
-        seam_point = projected_point.foot_pt;
-
-        //lastly, for internal perimeters, do the staggering if requested
-        if (po->config().staggered_inner_seams && loop.length() > 0.0) {
-            //fix depth, it is sometimes strongly underestimated
-            depth = std::max(loop.paths[projected_point.path_idx].width(), depth);
-
-            while (depth > 0.0f) {
-                auto next_point = get_next_loop_point(projected_point);
-                Vec2f a = unscale(projected_point.foot_pt).cast<float>();
-                Vec2f b = unscale(next_point.foot_pt).cast<float>();
-                float dist = (a - b).norm();
-                if (dist > depth) {
-                    Vec2f final_pos = a + (b - a) * depth / dist;
-                    next_point.foot_pt = Point::new_scale(final_pos.x(), final_pos.y());
-                }
-                depth -= dist;
-                projected_point = next_point;
-            }
-            seam_point = projected_point.foot_pt;
-        }
-    }
-
-    return seam_point;
 }
-
-} // namespace Slic3r
+} // namespace Slic3r::Seams
diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp
index c4b1edc00b..625170944a 100644
--- a/src/libslic3r/GCode/SeamPlacer.hpp
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamRandom.cpp b/src/libslic3r/GCode/SeamRandom.cpp
new file mode 100644
index 0000000000..27f01ca6be
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamRandom.hpp b/src/libslic3r/GCode/SeamRandom.hpp
new file mode 100644
index 0000000000..81c8dc4342
--- /dev/null
+++ b/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
+);
+}
diff --git a/src/libslic3r/GCode/SeamRear.cpp b/src/libslic3r/GCode/SeamRear.cpp
new file mode 100644
index 0000000000..d189b5ad08
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamRear.hpp b/src/libslic3r/GCode/SeamRear.hpp
new file mode 100644
index 0000000000..a313432780
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/GCode/SeamShells.cpp b/src/libslic3r/GCode/SeamShells.cpp
new file mode 100644
index 0000000000..f90ce0881c
--- /dev/null
+++ b/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
diff --git a/src/libslic3r/GCode/SeamShells.hpp b/src/libslic3r/GCode/SeamShells.hpp
new file mode 100644
index 0000000000..555332708e
--- /dev/null
+++ b/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_
diff --git a/src/libslic3r/Layer.cpp b/src/libslic3r/Layer.cpp
index e4b4d6539a..804999c796 100644
--- a/src/libslic3r/Layer.cpp
+++ b/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) {
diff --git a/src/libslic3r/Layer.hpp b/src/libslic3r/Layer.hpp
index 6b45c38925..0a2e457df9 100644
--- a/src/libslic3r/Layer.hpp
+++ b/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.
diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp
index 6917b24c7c..9089bae6a7 100644
--- a/src/libslic3r/ShortEdgeCollapse.cpp
+++ b/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]) {
diff --git a/src/libslic3r/TriangleSetSampling.cpp b/src/libslic3r/TriangleSetSampling.cpp
index 363df8659c..dd968fce77 100644
--- a/src/libslic3r/TriangleSetSampling.cpp
+++ b/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) };
diff --git a/tests/data/seam_test_object.3mf b/tests/data/seam_test_object.3mf
new file mode 100644
index 0000000000..c216b30725
Binary files /dev/null and b/tests/data/seam_test_object.3mf differ
diff --git a/tests/fff_print/CMakeLists.txt b/tests/fff_print/CMakeLists.txt
index da1d8b54c6..2872aa464b 100644
--- a/tests/fff_print/CMakeLists.txt
+++ b/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)
diff --git a/tests/fff_print/benchmark_seams.cpp b/tests/fff_print/benchmark_seams.cpp
new file mode 100644
index 0000000000..337c730a6d
--- /dev/null
+++ b/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});
+                }
+            });
+        };
+    }
+}
diff --git a/tests/fff_print/test_data.hpp b/tests/fff_print/test_data.hpp
index 12ff0f551d..61e1b92d42 100644
--- a/tests/fff_print/test_data.hpp
+++ b/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
diff --git a/tests/fff_print/test_seam_aligned.cpp b/tests/fff_print/test_seam_aligned.cpp
new file mode 100644
index 0000000000..66288c3ff5
--- /dev/null
+++ b/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
+                }
+            }
+        }
+    }
+}
diff --git a/tests/fff_print/test_seam_geometry.cpp b/tests/fff_print/test_seam_geometry.cpp
new file mode 100644
index 0000000000..a12fbdb3d6
--- /dev/null
+++ b/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
+    }));
+}
diff --git a/tests/fff_print/test_seam_perimeters.cpp b/tests/fff_print/test_seam_perimeters.cpp
new file mode 100644
index 0000000000..b97b59832f
--- /dev/null
+++ b/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);
+}
diff --git a/tests/fff_print/test_seam_random.cpp b/tests/fff_print/test_seam_random.cpp
new file mode 100644
index 0000000000..c877b39a70
--- /dev/null
+++ b/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);
+    }
+}
diff --git a/tests/fff_print/test_seam_rear.cpp b/tests/fff_print/test_seam_rear.cpp
new file mode 100644
index 0000000000..dda99b219e
--- /dev/null
+++ b/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);
+    }
+}
diff --git a/tests/fff_print/test_seam_shells.cpp b/tests/fff_print/test_seam_shells.cpp
new file mode 100644
index 0000000000..abbd4f3c45
--- /dev/null
+++ b/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);
+}