From 551ba962058aade6f952679f8bb14766ee5eb974 Mon Sep 17 00:00:00 2001
From: Lukas Matena <lukasmatena@seznam.cz>
Date: Wed, 27 Mar 2024 20:48:03 +0800
Subject: [PATCH] ENH: modify the multi-material segmentation and voronoi

This patch is cherry pick from Prusa, thanks to Prusa

Rework multi-material segmentation to work directly on the Voronoi diagram without creating a copy of it.

Previous algorithms assume that they can get an invalid Voronoi diagram. Because of that, during the multi-material segmentation, a copy of the Voronoi diagram was created, and there were several attempts to fix missing vertices and edges. But as it shows, this wasn't a good enough approach and sometimes led to several issues like bleeding layers.

After generalization, our approach for detection and repairs of invalid Voronoi diagrams from Arachne, we could assume that multi-material segmentation gets non-invalid Voronoi diagrams.
With this assumption, we reimplement multi-materials segmentation to work directly on the Voronoi diagram. That should make multi-material segmentation more stable.

So, this should fix several issues like bleeding layers. Also, memory consumption should decrease by a lot. Also, there should be some speedup of multi-materials segmentation.

Jira: none
Change-Id: I72aa6e1f9634d9ee8759aa469a0b39a36ace62f5
---
 src/ankerl/README.txt                         |    7 +
 src/ankerl/unordered_dense.h                  | 1584 +++++++++++++++++
 .../Arachne/SkeletalTrapezoidation.cpp        |  384 +---
 .../Arachne/SkeletalTrapezoidation.hpp        |   54 +-
 .../Arachne/utils/PolygonsSegmentIndex.hpp    |   19 +
 src/libslic3r/Arachne/utils/VoronoiUtils.cpp  |  251 ---
 src/libslic3r/Arachne/utils/VoronoiUtils.hpp  |   47 -
 src/libslic3r/CMakeLists.txt                  |    5 +-
 src/libslic3r/Geometry/MedialAxis.cpp         |    2 +-
 src/libslic3r/Geometry/Voronoi.cpp            |  354 ++++
 src/libslic3r/Geometry/Voronoi.hpp            |  186 +-
 src/libslic3r/Geometry/VoronoiOffset.cpp      |    6 -
 src/libslic3r/Geometry/VoronoiUtils.cpp       |  283 +++
 src/libslic3r/Geometry/VoronoiUtils.hpp       |  120 ++
 src/libslic3r/Geometry/VoronoiUtilsCgal.cpp   |  269 ++-
 src/libslic3r/Geometry/VoronoiUtilsCgal.hpp   |    9 +-
 src/libslic3r/MultiMaterialSegmentation.cpp   | 1376 +++++---------
 src/libslic3r/MultiMaterialSegmentation.hpp   |   30 +-
 18 files changed, 3304 insertions(+), 1682 deletions(-)
 create mode 100644 src/ankerl/README.txt
 create mode 100644 src/ankerl/unordered_dense.h
 delete mode 100644 src/libslic3r/Arachne/utils/VoronoiUtils.cpp
 delete mode 100644 src/libslic3r/Arachne/utils/VoronoiUtils.hpp
 create mode 100644 src/libslic3r/Geometry/Voronoi.cpp
 create mode 100644 src/libslic3r/Geometry/VoronoiUtils.cpp
 create mode 100644 src/libslic3r/Geometry/VoronoiUtils.hpp

diff --git a/src/ankerl/README.txt b/src/ankerl/README.txt
new file mode 100644
index 000000000..0996bc28c
--- /dev/null
+++ b/src/ankerl/README.txt
@@ -0,0 +1,7 @@
+THIS DIRECTORY CONTAINS PIECES OF THE 
+ankerl::unordered_dense::{map, set}
+https://github.com/martinus/unordered_dense
+unordered_dense 3.1.1 10782bfc651c2bb75b11bf90491f50da122e5432
+SOURCE DISTRIBUTION.
+
+THIS IS NOT THE COMPLETE unordered_dense DISTRIBUTION. ONLY FILES NEEDED FOR COMPILING PRUSASLICER WERE PUT INTO THE PRUSASLICER SOURCE DISTRIBUTION.
diff --git a/src/ankerl/unordered_dense.h b/src/ankerl/unordered_dense.h
new file mode 100644
index 000000000..e294bdb4e
--- /dev/null
+++ b/src/ankerl/unordered_dense.h
@@ -0,0 +1,1584 @@
+///////////////////////// ankerl::unordered_dense::{map, set} /////////////////////////
+
+// A fast & densely stored hashmap and hashset based on robin-hood backward shift deletion.
+// Version 3.1.1
+// https://github.com/martinus/unordered_dense
+//
+// Licensed under the MIT License <http://opensource.org/licenses/MIT>.
+// SPDX-License-Identifier: MIT
+// Copyright (c) 2022-2023 Martin Leitner-Ankerl <martin.ankerl@gmail.com>
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#ifndef ANKERL_UNORDERED_DENSE_H
+#define ANKERL_UNORDERED_DENSE_H
+
+// see https://semver.org/spec/v2.0.0.html
+#define ANKERL_UNORDERED_DENSE_VERSION_MAJOR 3 // NOLINT(cppcoreguidelines-macro-usage) incompatible API changes
+#define ANKERL_UNORDERED_DENSE_VERSION_MINOR 1 // NOLINT(cppcoreguidelines-macro-usage) backwards compatible functionality
+#define ANKERL_UNORDERED_DENSE_VERSION_PATCH 1 // NOLINT(cppcoreguidelines-macro-usage) backwards compatible bug fixes
+
+// API versioning with inline namespace, see https://www.foonathan.net/2018/11/inline-namespaces/
+#define ANKERL_UNORDERED_DENSE_VERSION_CONCAT1(major, minor, patch) v##major##_##minor##_##patch
+#define ANKERL_UNORDERED_DENSE_VERSION_CONCAT(major, minor, patch) ANKERL_UNORDERED_DENSE_VERSION_CONCAT1(major, minor, patch)
+#define ANKERL_UNORDERED_DENSE_NAMESPACE   \
+    ANKERL_UNORDERED_DENSE_VERSION_CONCAT( \
+        ANKERL_UNORDERED_DENSE_VERSION_MAJOR, ANKERL_UNORDERED_DENSE_VERSION_MINOR, ANKERL_UNORDERED_DENSE_VERSION_PATCH)
+
+#if defined(_MSVC_LANG)
+#    define ANKERL_UNORDERED_DENSE_CPP_VERSION _MSVC_LANG
+#else
+#    define ANKERL_UNORDERED_DENSE_CPP_VERSION __cplusplus
+#endif
+
+#if defined(__GNUC__)
+// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
+#    define ANKERL_UNORDERED_DENSE_PACK(decl) decl __attribute__((__packed__))
+#elif defined(_MSC_VER)
+// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
+#    define ANKERL_UNORDERED_DENSE_PACK(decl) __pragma(pack(push, 1)) decl __pragma(pack(pop))
+#endif
+
+// exceptions
+#if defined(__cpp_exceptions) || defined(__EXCEPTIONS) || defined(_CPPUNWIND)
+#    define ANKERL_UNORDERED_DENSE_HAS_EXCEPTIONS() 1
+#else
+#    define ANKERL_UNORDERED_DENSE_HAS_EXCEPTIONS() 0
+#endif
+#ifdef _MSC_VER
+#    define ANKERL_UNORDERED_DENSE_NOINLINE __declspec(noinline)
+#else
+#    define ANKERL_UNORDERED_DENSE_NOINLINE __attribute__((noinline))
+#endif
+
+#if ANKERL_UNORDERED_DENSE_CPP_VERSION < 201703L
+#    error ankerl::unordered_dense requires C++17 or higher
+#else
+#    include <array>            // for array
+#    include <cstdint>          // for uint64_t, uint32_t, uint8_t, UINT64_C
+#    include <cstring>          // for size_t, memcpy, memset
+#    include <functional>       // for equal_to, hash
+#    include <initializer_list> // for initializer_list
+#    include <iterator>         // for pair, distance
+#    include <limits>           // for numeric_limits
+#    include <memory>           // for allocator, allocator_traits, shared_ptr
+#    include <stdexcept>        // for out_of_range
+#    include <string>           // for basic_string
+#    include <string_view>      // for basic_string_view, hash
+#    include <tuple>            // for forward_as_tuple
+#    include <type_traits>      // for enable_if_t, declval, conditional_t, ena...
+#    include <utility>          // for forward, exchange, pair, as_const, piece...
+#    include <vector>           // for vector
+#    if ANKERL_UNORDERED_DENSE_HAS_EXCEPTIONS() == 0
+#        include <cstdlib> // for abort
+#    endif
+
+#    define ANKERL_UNORDERED_DENSE_PMR 0 // NOLINT(cppcoreguidelines-macro-usage)
+#    if defined(__has_include)
+#        if __has_include(<memory_resource>)
+#            undef ANKERL_UNORDERED_DENSE_PMR
+#            define ANKERL_UNORDERED_DENSE_PMR 1 // NOLINT(cppcoreguidelines-macro-usage)
+#            define ANKERL_UNORDERED_DENSE_PMR_ALLOCATOR \
+                std::pmr::polymorphic_allocator // NOLINT(cppcoreguidelines-macro-usage)
+#            include <memory_resource>          // for polymorphic_allocator
+#        elif __has_include(<experimental/memory_resource>)
+#            undef ANKERL_UNORDERED_DENSE_PMR
+#            define ANKERL_UNORDERED_DENSE_PMR 1 // NOLINT(cppcoreguidelines-macro-usage)
+#            define ANKERL_UNORDERED_DENSE_PMR_ALLOCATOR \
+                std::experimental::pmr::polymorphic_allocator // NOLINT(cppcoreguidelines-macro-usage)
+#            include <experimental/memory_resource>           // for polymorphic_allocator
+#        endif
+#    endif
+
+#    if defined(_MSC_VER) && defined(_M_X64)
+#        include <intrin.h>
+#        pragma intrinsic(_umul128)
+#    endif
+
+#    if defined(__GNUC__) || defined(__INTEL_COMPILER) || defined(__clang__)
+#        define ANKERL_UNORDERED_DENSE_LIKELY(x) __builtin_expect(x, 1)   // NOLINT(cppcoreguidelines-macro-usage)
+#        define ANKERL_UNORDERED_DENSE_UNLIKELY(x) __builtin_expect(x, 0) // NOLINT(cppcoreguidelines-macro-usage)
+#    else
+#        define ANKERL_UNORDERED_DENSE_LIKELY(x) (x)   // NOLINT(cppcoreguidelines-macro-usage)
+#        define ANKERL_UNORDERED_DENSE_UNLIKELY(x) (x) // NOLINT(cppcoreguidelines-macro-usage)
+#    endif
+
+namespace ankerl::unordered_dense {
+inline namespace ANKERL_UNORDERED_DENSE_NAMESPACE {
+
+namespace detail {
+
+#    if ANKERL_UNORDERED_DENSE_HAS_EXCEPTIONS()
+
+// make sure this is not inlined as it is slow and dramatically enlarges code, thus making other
+// inlinings more difficult. Throws are also generally the slow path.
+[[noreturn]] inline ANKERL_UNORDERED_DENSE_NOINLINE void on_error_key_not_found() {
+    throw std::out_of_range("ankerl::unordered_dense::map::at(): key not found");
+}
+[[noreturn]] inline ANKERL_UNORDERED_DENSE_NOINLINE void on_error_bucket_overflow() {
+    throw std::overflow_error("ankerl::unordered_dense: reached max bucket size, cannot increase size");
+}
+[[noreturn]] inline ANKERL_UNORDERED_DENSE_NOINLINE void on_error_too_many_elements() {
+    throw std::out_of_range("ankerl::unordered_dense::map::replace(): too many elements");
+}
+
+#    else
+
+[[noreturn]] inline void on_error_key_not_found() {
+    abort();
+}
+[[noreturn]] inline void on_error_bucket_overflow() {
+    abort();
+}
+[[noreturn]] inline void on_error_too_many_elements() {
+    abort();
+}
+
+#    endif
+
+} // namespace detail
+
+// hash ///////////////////////////////////////////////////////////////////////
+
+// This is a stripped-down implementation of wyhash: https://github.com/wangyi-fudan/wyhash
+// No big-endian support (because different values on different machines don't matter),
+// hardcodes seed and the secret, reformattes the code, and clang-tidy fixes.
+namespace detail::wyhash {
+
+static inline void mum(uint64_t* a, uint64_t* b) {
+#    if defined(__SIZEOF_INT128__)
+    __uint128_t r = *a;
+    r *= *b;
+    *a = static_cast<uint64_t>(r);
+    *b = static_cast<uint64_t>(r >> 64U);
+#    elif defined(_MSC_VER) && defined(_M_X64)
+    *a = _umul128(*a, *b, b);
+#    else
+    uint64_t ha = *a >> 32U;
+    uint64_t hb = *b >> 32U;
+    uint64_t la = static_cast<uint32_t>(*a);
+    uint64_t lb = static_cast<uint32_t>(*b);
+    uint64_t hi{};
+    uint64_t lo{};
+    uint64_t rh = ha * hb;
+    uint64_t rm0 = ha * lb;
+    uint64_t rm1 = hb * la;
+    uint64_t rl = la * lb;
+    uint64_t t = rl + (rm0 << 32U);
+    auto c = static_cast<uint64_t>(t < rl);
+    lo = t + (rm1 << 32U);
+    c += static_cast<uint64_t>(lo < t);
+    hi = rh + (rm0 >> 32U) + (rm1 >> 32U) + c;
+    *a = lo;
+    *b = hi;
+#    endif
+}
+
+// multiply and xor mix function, aka MUM
+[[nodiscard]] static inline auto mix(uint64_t a, uint64_t b) -> uint64_t {
+    mum(&a, &b);
+    return a ^ b;
+}
+
+// read functions. WARNING: we don't care about endianness, so results are different on big endian!
+[[nodiscard]] static inline auto r8(const uint8_t* p) -> uint64_t {
+    uint64_t v{};
+    std::memcpy(&v, p, 8U);
+    return v;
+}
+
+[[nodiscard]] static inline auto r4(const uint8_t* p) -> uint64_t {
+    uint32_t v{};
+    std::memcpy(&v, p, 4);
+    return v;
+}
+
+// reads 1, 2, or 3 bytes
+[[nodiscard]] static inline auto r3(const uint8_t* p, size_t k) -> uint64_t {
+    return (static_cast<uint64_t>(p[0]) << 16U) | (static_cast<uint64_t>(p[k >> 1U]) << 8U) | p[k - 1];
+}
+
+[[maybe_unused]] [[nodiscard]] static inline auto hash(void const* key, size_t len) -> uint64_t {
+    static constexpr auto secret = std::array{UINT64_C(0xa0761d6478bd642f),
+                                              UINT64_C(0xe7037ed1a0b428db),
+                                              UINT64_C(0x8ebc6af09c88c6e3),
+                                              UINT64_C(0x589965cc75374cc3)};
+
+    auto const* p = static_cast<uint8_t const*>(key);
+    uint64_t seed = secret[0];
+    uint64_t a{};
+    uint64_t b{};
+    if (ANKERL_UNORDERED_DENSE_LIKELY(len <= 16)) {
+        if (ANKERL_UNORDERED_DENSE_LIKELY(len >= 4)) {
+            a = (r4(p) << 32U) | r4(p + ((len >> 3U) << 2U));
+            b = (r4(p + len - 4) << 32U) | r4(p + len - 4 - ((len >> 3U) << 2U));
+        } else if (ANKERL_UNORDERED_DENSE_LIKELY(len > 0)) {
+            a = r3(p, len);
+            b = 0;
+        } else {
+            a = 0;
+            b = 0;
+        }
+    } else {
+        size_t i = len;
+        if (ANKERL_UNORDERED_DENSE_UNLIKELY(i > 48)) {
+            uint64_t see1 = seed;
+            uint64_t see2 = seed;
+            do {
+                seed = mix(r8(p) ^ secret[1], r8(p + 8) ^ seed);
+                see1 = mix(r8(p + 16) ^ secret[2], r8(p + 24) ^ see1);
+                see2 = mix(r8(p + 32) ^ secret[3], r8(p + 40) ^ see2);
+                p += 48;
+                i -= 48;
+            } while (ANKERL_UNORDERED_DENSE_LIKELY(i > 48));
+            seed ^= see1 ^ see2;
+        }
+        while (ANKERL_UNORDERED_DENSE_UNLIKELY(i > 16)) {
+            seed = mix(r8(p) ^ secret[1], r8(p + 8) ^ seed);
+            i -= 16;
+            p += 16;
+        }
+        a = r8(p + i - 16);
+        b = r8(p + i - 8);
+    }
+
+    return mix(secret[1] ^ len, mix(a ^ secret[1], b ^ seed));
+}
+
+[[nodiscard]] static inline auto hash(uint64_t x) -> uint64_t {
+    return detail::wyhash::mix(x, UINT64_C(0x9E3779B97F4A7C15));
+}
+
+} // namespace detail::wyhash
+
+template <typename T, typename Enable = void>
+struct hash {
+    auto operator()(T const& obj) const noexcept(noexcept(std::declval<std::hash<T>>().operator()(std::declval<T const&>())))
+        -> uint64_t {
+        return std::hash<T>{}(obj);
+    }
+};
+
+template <typename CharT>
+struct hash<std::basic_string<CharT>> {
+    using is_avalanching = void;
+    auto operator()(std::basic_string<CharT> const& str) const noexcept -> uint64_t {
+        return detail::wyhash::hash(str.data(), sizeof(CharT) * str.size());
+    }
+};
+
+template <typename CharT>
+struct hash<std::basic_string_view<CharT>> {
+    using is_avalanching = void;
+    auto operator()(std::basic_string_view<CharT> const& sv) const noexcept -> uint64_t {
+        return detail::wyhash::hash(sv.data(), sizeof(CharT) * sv.size());
+    }
+};
+
+template <class T>
+struct hash<T*> {
+    using is_avalanching = void;
+    auto operator()(T* ptr) const noexcept -> uint64_t {
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        return detail::wyhash::hash(reinterpret_cast<uintptr_t>(ptr));
+    }
+};
+
+template <class T>
+struct hash<std::unique_ptr<T>> {
+    using is_avalanching = void;
+    auto operator()(std::unique_ptr<T> const& ptr) const noexcept -> uint64_t {
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        return detail::wyhash::hash(reinterpret_cast<uintptr_t>(ptr.get()));
+    }
+};
+
+template <class T>
+struct hash<std::shared_ptr<T>> {
+    using is_avalanching = void;
+    auto operator()(std::shared_ptr<T> const& ptr) const noexcept -> uint64_t {
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        return detail::wyhash::hash(reinterpret_cast<uintptr_t>(ptr.get()));
+    }
+};
+
+template <typename Enum>
+struct hash<Enum, typename std::enable_if<std::is_enum<Enum>::value>::type> {
+    using is_avalanching = void;
+    auto operator()(Enum e) const noexcept -> uint64_t {
+        using underlying = typename std::underlying_type_t<Enum>;
+        return detail::wyhash::hash(static_cast<underlying>(e));
+    }
+};
+
+// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
+#    define ANKERL_UNORDERED_DENSE_HASH_STATICCAST(T)                    \
+        template <>                                                      \
+        struct hash<T> {                                                 \
+            using is_avalanching = void;                                 \
+            auto operator()(T const& obj) const noexcept -> uint64_t {   \
+                return detail::wyhash::hash(static_cast<uint64_t>(obj)); \
+            }                                                            \
+        }
+
+#    if defined(__GNUC__) && !defined(__clang__)
+#        pragma GCC diagnostic push
+#        pragma GCC diagnostic ignored "-Wuseless-cast"
+#    endif
+// see https://en.cppreference.com/w/cpp/utility/hash
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(bool);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(char);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(signed char);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(unsigned char);
+#    if ANKERL_UNORDERED_DENSE_CPP_VERSION >= 202002L
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(char8_t);
+#    endif
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(char16_t);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(char32_t);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(wchar_t);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(short);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(unsigned short);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(int);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(unsigned int);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(long);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(long long);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(unsigned long);
+ANKERL_UNORDERED_DENSE_HASH_STATICCAST(unsigned long long);
+
+#    if defined(__GNUC__) && !defined(__clang__)
+#        pragma GCC diagnostic pop
+#    endif
+
+// bucket_type //////////////////////////////////////////////////////////
+
+namespace bucket_type {
+
+struct standard {
+    static constexpr uint32_t dist_inc = 1U << 8U;             // skip 1 byte fingerprint
+    static constexpr uint32_t fingerprint_mask = dist_inc - 1; // mask for 1 byte of fingerprint
+
+    uint32_t m_dist_and_fingerprint; // upper 3 byte: distance to original bucket. lower byte: fingerprint from hash
+    uint32_t m_value_idx;            // index into the m_values vector.
+};
+
+ANKERL_UNORDERED_DENSE_PACK(struct big {
+    static constexpr uint32_t dist_inc = 1U << 8U;             // skip 1 byte fingerprint
+    static constexpr uint32_t fingerprint_mask = dist_inc - 1; // mask for 1 byte of fingerprint
+
+    uint32_t m_dist_and_fingerprint; // upper 3 byte: distance to original bucket. lower byte: fingerprint from hash
+    size_t m_value_idx;              // index into the m_values vector.
+});
+
+} // namespace bucket_type
+
+namespace detail {
+
+struct nonesuch {};
+
+template <class Default, class AlwaysVoid, template <class...> class Op, class... Args>
+struct detector {
+    using value_t = std::false_type;
+    using type = Default;
+};
+
+template <class Default, template <class...> class Op, class... Args>
+struct detector<Default, std::void_t<Op<Args...>>, Op, Args...> {
+    using value_t = std::true_type;
+    using type = Op<Args...>;
+};
+
+template <template <class...> class Op, class... Args>
+using is_detected = typename detail::detector<detail::nonesuch, void, Op, Args...>::value_t;
+
+template <template <class...> class Op, class... Args>
+constexpr bool is_detected_v = is_detected<Op, Args...>::value;
+
+template <typename T>
+using detect_avalanching = typename T::is_avalanching;
+
+template <typename T>
+using detect_is_transparent = typename T::is_transparent;
+
+template <typename T>
+using detect_iterator = typename T::iterator;
+
+template <typename T>
+using detect_reserve = decltype(std::declval<T&>().reserve(size_t{}));
+
+// enable_if helpers
+
+template <typename Mapped>
+constexpr bool is_map_v = !std::is_void_v<Mapped>;
+
+// clang-format off
+template <typename Hash, typename KeyEqual>
+constexpr bool is_transparent_v = is_detected_v<detect_is_transparent, Hash>&& is_detected_v<detect_is_transparent, KeyEqual>;
+// clang-format on
+
+template <typename From, typename To1, typename To2>
+constexpr bool is_neither_convertible_v = !std::is_convertible_v<From, To1> && !std::is_convertible_v<From, To2>;
+
+template <typename T>
+constexpr bool has_reserve = is_detected_v<detect_reserve, T>;
+
+// base type for map has mapped_type
+template <class T>
+struct base_table_type_map {
+    using mapped_type = T;
+};
+
+// base type for set doesn't have mapped_type
+struct base_table_type_set {};
+
+// This is it, the table. Doubles as map and set, and uses `void` for T when its used as a set.
+template <class Key,
+          class T, // when void, treat it as a set.
+          class Hash,
+          class KeyEqual,
+          class AllocatorOrContainer,
+          class Bucket>
+class table : public std::conditional_t<is_map_v<T>, base_table_type_map<T>, base_table_type_set> {
+public:
+    using value_container_type = std::conditional_t<
+        is_detected_v<detect_iterator, AllocatorOrContainer>,
+        AllocatorOrContainer,
+        typename std::vector<typename std::conditional_t<is_map_v<T>, std::pair<Key, T>, Key>, AllocatorOrContainer>>;
+
+private:
+    using bucket_alloc =
+        typename std::allocator_traits<typename value_container_type::allocator_type>::template rebind_alloc<Bucket>;
+    using bucket_alloc_traits = std::allocator_traits<bucket_alloc>;
+
+    static constexpr uint8_t initial_shifts = 64 - 3; // 2^(64-m_shift) number of buckets
+    static constexpr float default_max_load_factor = 0.8F;
+
+public:
+    using key_type = Key;
+    using value_type = typename value_container_type::value_type;
+    using size_type = typename value_container_type::size_type;
+    using difference_type = typename value_container_type::difference_type;
+    using hasher = Hash;
+    using key_equal = KeyEqual;
+    using allocator_type = typename value_container_type::allocator_type;
+    using reference = typename value_container_type::reference;
+    using const_reference = typename value_container_type::const_reference;
+    using pointer = typename value_container_type::pointer;
+    using const_pointer = typename value_container_type::const_pointer;
+    using const_iterator = typename value_container_type::const_iterator;
+    using iterator = std::conditional_t<is_map_v<T>, typename value_container_type::iterator, const_iterator>;
+    using bucket_type = Bucket;
+
+private:
+    using value_idx_type = decltype(Bucket::m_value_idx);
+    using dist_and_fingerprint_type = decltype(Bucket::m_dist_and_fingerprint);
+
+    static_assert(std::is_trivially_destructible_v<Bucket>, "assert there's no need to call destructor / std::destroy");
+    static_assert(std::is_trivially_copyable_v<Bucket>, "assert we can just memset / memcpy");
+
+    value_container_type m_values{}; // Contains all the key-value pairs in one densely stored container. No holes.
+    typename std::allocator_traits<bucket_alloc>::pointer m_buckets{};
+    size_t m_num_buckets = 0;
+    size_t m_max_bucket_capacity = 0;
+    float m_max_load_factor = default_max_load_factor;
+    Hash m_hash{};
+    KeyEqual m_equal{};
+    uint8_t m_shifts = initial_shifts;
+
+    [[nodiscard]] auto next(value_idx_type bucket_idx) const -> value_idx_type {
+        return ANKERL_UNORDERED_DENSE_UNLIKELY(bucket_idx + 1U == m_num_buckets)
+                   ? 0
+                   : static_cast<value_idx_type>(bucket_idx + 1U);
+    }
+
+    // Helper to access bucket through pointer types
+    [[nodiscard]] static constexpr auto at(typename std::allocator_traits<bucket_alloc>::pointer bucket_ptr, size_t offset)
+        -> Bucket& {
+        return *(bucket_ptr + static_cast<typename std::allocator_traits<bucket_alloc>::difference_type>(offset));
+    }
+
+    // use the dist_inc and dist_dec functions so that uint16_t types work without warning
+    [[nodiscard]] static constexpr auto dist_inc(dist_and_fingerprint_type x) -> dist_and_fingerprint_type {
+        return static_cast<dist_and_fingerprint_type>(x + Bucket::dist_inc);
+    }
+
+    [[nodiscard]] static constexpr auto dist_dec(dist_and_fingerprint_type x) -> dist_and_fingerprint_type {
+        return static_cast<dist_and_fingerprint_type>(x - Bucket::dist_inc);
+    }
+
+    // The goal of mixed_hash is to always produce a high quality 64bit hash.
+    template <typename K>
+    [[nodiscard]] constexpr auto mixed_hash(K const& key) const -> uint64_t {
+        if constexpr (is_detected_v<detect_avalanching, Hash>) {
+            // we know that the hash is good because is_avalanching.
+            if constexpr (sizeof(decltype(m_hash(key))) < sizeof(uint64_t)) {
+                // 32bit hash and is_avalanching => multiply with a constant to avalanche bits upwards
+                return m_hash(key) * UINT64_C(0x9ddfea08eb382d69);
+            } else {
+                // 64bit and is_avalanching => only use the hash itself.
+                return m_hash(key);
+            }
+        } else {
+            // not is_avalanching => apply wyhash
+            return wyhash::hash(m_hash(key));
+        }
+    }
+
+    [[nodiscard]] constexpr auto dist_and_fingerprint_from_hash(uint64_t hash) const -> dist_and_fingerprint_type {
+        return Bucket::dist_inc | (static_cast<dist_and_fingerprint_type>(hash) & Bucket::fingerprint_mask);
+    }
+
+    [[nodiscard]] constexpr auto bucket_idx_from_hash(uint64_t hash) const -> value_idx_type {
+        return static_cast<value_idx_type>(hash >> m_shifts);
+    }
+
+    [[nodiscard]] static constexpr auto get_key(value_type const& vt) -> key_type const& {
+        if constexpr (is_map_v<T>) {
+            return vt.first;
+        } else {
+            return vt;
+        }
+    }
+
+    template <typename K>
+    [[nodiscard]] auto next_while_less(K const& key) const -> Bucket {
+        auto hash = mixed_hash(key);
+        auto dist_and_fingerprint = dist_and_fingerprint_from_hash(hash);
+        auto bucket_idx = bucket_idx_from_hash(hash);
+
+        while (dist_and_fingerprint < at(m_buckets, bucket_idx).m_dist_and_fingerprint) {
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+        }
+        return {dist_and_fingerprint, bucket_idx};
+    }
+
+    void place_and_shift_up(Bucket bucket, value_idx_type place) {
+        while (0 != at(m_buckets, place).m_dist_and_fingerprint) {
+            bucket = std::exchange(at(m_buckets, place), bucket);
+            bucket.m_dist_and_fingerprint = dist_inc(bucket.m_dist_and_fingerprint);
+            place = next(place);
+        }
+        at(m_buckets, place) = bucket;
+    }
+
+    [[nodiscard]] static constexpr auto calc_num_buckets(uint8_t shifts) -> size_t {
+        return (std::min)(max_bucket_count(), size_t{1} << (64U - shifts));
+    }
+
+    [[nodiscard]] constexpr auto calc_shifts_for_size(size_t s) const -> uint8_t {
+        auto shifts = initial_shifts;
+        while (shifts > 0 && static_cast<size_t>(static_cast<float>(calc_num_buckets(shifts)) * max_load_factor()) < s) {
+            --shifts;
+        }
+        return shifts;
+    }
+
+    // assumes m_values has data, m_buckets=m_buckets_end=nullptr, m_shifts is INITIAL_SHIFTS
+    void copy_buckets(table const& other) {
+        if (!empty()) {
+            m_shifts = other.m_shifts;
+            allocate_buckets_from_shift();
+            std::memcpy(m_buckets, other.m_buckets, sizeof(Bucket) * bucket_count());
+        }
+    }
+
+    /**
+     * True when no element can be added any more without increasing the size
+     */
+    [[nodiscard]] auto is_full() const -> bool {
+        return size() >= m_max_bucket_capacity;
+    }
+
+    void deallocate_buckets() {
+        auto ba = bucket_alloc(m_values.get_allocator());
+        if (nullptr != m_buckets) {
+            bucket_alloc_traits::deallocate(ba, m_buckets, bucket_count());
+        }
+        m_buckets = nullptr;
+        m_num_buckets = 0;
+        m_max_bucket_capacity = 0;
+    }
+
+    void allocate_buckets_from_shift() {
+        auto ba = bucket_alloc(m_values.get_allocator());
+        m_num_buckets = calc_num_buckets(m_shifts);
+        m_buckets = bucket_alloc_traits::allocate(ba, m_num_buckets);
+        if (m_num_buckets == max_bucket_count()) {
+            // reached the maximum, make sure we can use each bucket
+            m_max_bucket_capacity = max_bucket_count();
+        } else {
+            m_max_bucket_capacity = static_cast<value_idx_type>(static_cast<float>(m_num_buckets) * max_load_factor());
+        }
+    }
+
+    void clear_buckets() {
+        if (m_buckets != nullptr) {
+            std::memset(&*m_buckets, 0, sizeof(Bucket) * bucket_count());
+        }
+    }
+
+    void clear_and_fill_buckets_from_values() {
+        clear_buckets();
+        for (value_idx_type value_idx = 0, end_idx = static_cast<value_idx_type>(m_values.size()); value_idx < end_idx;
+             ++value_idx) {
+            auto const& key = get_key(m_values[value_idx]);
+            auto [dist_and_fingerprint, bucket] = next_while_less(key);
+
+            // we know for certain that key has not yet been inserted, so no need to check it.
+            place_and_shift_up({dist_and_fingerprint, value_idx}, bucket);
+        }
+    }
+
+    void increase_size() {
+        if (ANKERL_UNORDERED_DENSE_UNLIKELY(m_max_bucket_capacity == max_bucket_count())) {
+            on_error_bucket_overflow();
+        }
+        --m_shifts;
+        deallocate_buckets();
+        allocate_buckets_from_shift();
+        clear_and_fill_buckets_from_values();
+    }
+
+    void do_erase(value_idx_type bucket_idx) {
+        auto const value_idx_to_remove = at(m_buckets, bucket_idx).m_value_idx;
+
+        // shift down until either empty or an element with correct spot is found
+        auto next_bucket_idx = next(bucket_idx);
+        while (at(m_buckets, next_bucket_idx).m_dist_and_fingerprint >= Bucket::dist_inc * 2) {
+            at(m_buckets, bucket_idx) = {dist_dec(at(m_buckets, next_bucket_idx).m_dist_and_fingerprint),
+                                         at(m_buckets, next_bucket_idx).m_value_idx};
+            bucket_idx = std::exchange(next_bucket_idx, next(next_bucket_idx));
+        }
+        at(m_buckets, bucket_idx) = {};
+
+        // update m_values
+        if (value_idx_to_remove != m_values.size() - 1) {
+            // no luck, we'll have to replace the value with the last one and update the index accordingly
+            auto& val = m_values[value_idx_to_remove];
+            val = std::move(m_values.back());
+
+            // update the values_idx of the moved entry. No need to play the info game, just look until we find the values_idx
+            auto mh = mixed_hash(get_key(val));
+            bucket_idx = bucket_idx_from_hash(mh);
+
+            auto const values_idx_back = static_cast<value_idx_type>(m_values.size() - 1);
+            while (values_idx_back != at(m_buckets, bucket_idx).m_value_idx) {
+                bucket_idx = next(bucket_idx);
+            }
+            at(m_buckets, bucket_idx).m_value_idx = value_idx_to_remove;
+        }
+        m_values.pop_back();
+    }
+
+    template <typename K>
+    auto do_erase_key(K&& key) -> size_t {
+        if (empty()) {
+            return 0;
+        }
+
+        auto [dist_and_fingerprint, bucket_idx] = next_while_less(key);
+
+        while (dist_and_fingerprint == at(m_buckets, bucket_idx).m_dist_and_fingerprint &&
+               !m_equal(key, get_key(m_values[at(m_buckets, bucket_idx).m_value_idx]))) {
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+        }
+
+        if (dist_and_fingerprint != at(m_buckets, bucket_idx).m_dist_and_fingerprint) {
+            return 0;
+        }
+        do_erase(bucket_idx);
+        return 1;
+    }
+
+    template <class K, class M>
+    auto do_insert_or_assign(K&& key, M&& mapped) -> std::pair<iterator, bool> {
+        auto it_isinserted = try_emplace(std::forward<K>(key), std::forward<M>(mapped));
+        if (!it_isinserted.second) {
+            it_isinserted.first->second = std::forward<M>(mapped);
+        }
+        return it_isinserted;
+    }
+
+    template <typename K, typename... Args>
+    auto do_place_element(dist_and_fingerprint_type dist_and_fingerprint, value_idx_type bucket_idx, K&& key, Args&&... args)
+        -> std::pair<iterator, bool> {
+
+        // emplace the new value. If that throws an exception, no harm done; index is still in a valid state
+        m_values.emplace_back(std::piecewise_construct,
+                              std::forward_as_tuple(std::forward<K>(key)),
+                              std::forward_as_tuple(std::forward<Args>(args)...));
+
+        // place element and shift up until we find an empty spot
+        auto value_idx = static_cast<value_idx_type>(m_values.size() - 1);
+        place_and_shift_up({dist_and_fingerprint, value_idx}, bucket_idx);
+        return {begin() + static_cast<difference_type>(value_idx), true};
+    }
+
+    template <typename K, typename... Args>
+    auto do_try_emplace(K&& key, Args&&... args) -> std::pair<iterator, bool> {
+        if (ANKERL_UNORDERED_DENSE_UNLIKELY(is_full())) {
+            increase_size();
+        }
+
+        auto hash = mixed_hash(key);
+        auto dist_and_fingerprint = dist_and_fingerprint_from_hash(hash);
+        auto bucket_idx = bucket_idx_from_hash(hash);
+
+        while (true) {
+            auto* bucket = &at(m_buckets, bucket_idx);
+            if (dist_and_fingerprint == bucket->m_dist_and_fingerprint) {
+                if (m_equal(key, m_values[bucket->m_value_idx].first)) {
+                    return {begin() + static_cast<difference_type>(bucket->m_value_idx), false};
+                }
+            } else if (dist_and_fingerprint > bucket->m_dist_and_fingerprint) {
+                return do_place_element(dist_and_fingerprint, bucket_idx, std::forward<K>(key), std::forward<Args>(args)...);
+            }
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+        }
+    }
+
+    template <typename K>
+    auto do_find(K const& key) -> iterator {
+        if (ANKERL_UNORDERED_DENSE_UNLIKELY(empty())) {
+            return end();
+        }
+
+        auto mh = mixed_hash(key);
+        auto dist_and_fingerprint = dist_and_fingerprint_from_hash(mh);
+        auto bucket_idx = bucket_idx_from_hash(mh);
+        auto* bucket = &at(m_buckets, bucket_idx);
+
+        // unrolled loop. *Always* check a few directly, then enter the loop. This is faster.
+        if (dist_and_fingerprint == bucket->m_dist_and_fingerprint && m_equal(key, get_key(m_values[bucket->m_value_idx]))) {
+            return begin() + static_cast<difference_type>(bucket->m_value_idx);
+        }
+        dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+        bucket_idx = next(bucket_idx);
+        bucket = &at(m_buckets, bucket_idx);
+
+        if (dist_and_fingerprint == bucket->m_dist_and_fingerprint && m_equal(key, get_key(m_values[bucket->m_value_idx]))) {
+            return begin() + static_cast<difference_type>(bucket->m_value_idx);
+        }
+        dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+        bucket_idx = next(bucket_idx);
+        bucket = &at(m_buckets, bucket_idx);
+
+        while (true) {
+            if (dist_and_fingerprint == bucket->m_dist_and_fingerprint) {
+                if (m_equal(key, get_key(m_values[bucket->m_value_idx]))) {
+                    return begin() + static_cast<difference_type>(bucket->m_value_idx);
+                }
+            } else if (dist_and_fingerprint > bucket->m_dist_and_fingerprint) {
+                return end();
+            }
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+            bucket = &at(m_buckets, bucket_idx);
+        }
+    }
+
+    template <typename K>
+    auto do_find(K const& key) const -> const_iterator {
+        return const_cast<table*>(this)->do_find(key); // NOLINT(cppcoreguidelines-pro-type-const-cast)
+    }
+
+    template <typename K, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto do_at(K const& key) -> Q& {
+        if (auto it = find(key); ANKERL_UNORDERED_DENSE_LIKELY(end() != it)) {
+            return it->second;
+        }
+        on_error_key_not_found();
+    }
+
+    template <typename K, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto do_at(K const& key) const -> Q const& {
+        return const_cast<table*>(this)->at(key); // NOLINT(cppcoreguidelines-pro-type-const-cast)
+    }
+
+public:
+    table()
+        : table(0) {}
+
+    explicit table(size_t bucket_count,
+                   Hash const& hash = Hash(),
+                   KeyEqual const& equal = KeyEqual(),
+                   allocator_type const& alloc_or_container = allocator_type())
+        : m_values(alloc_or_container)
+        , m_hash(hash)
+        , m_equal(equal) {
+        if (0 != bucket_count) {
+            reserve(bucket_count);
+        }
+    }
+
+    table(size_t bucket_count, allocator_type const& alloc)
+        : table(bucket_count, Hash(), KeyEqual(), alloc) {}
+
+    table(size_t bucket_count, Hash const& hash, allocator_type const& alloc)
+        : table(bucket_count, hash, KeyEqual(), alloc) {}
+
+    explicit table(allocator_type const& alloc)
+        : table(0, Hash(), KeyEqual(), alloc) {}
+
+    template <class InputIt>
+    table(InputIt first,
+          InputIt last,
+          size_type bucket_count = 0,
+          Hash const& hash = Hash(),
+          KeyEqual const& equal = KeyEqual(),
+          allocator_type const& alloc = allocator_type())
+        : table(bucket_count, hash, equal, alloc) {
+        insert(first, last);
+    }
+
+    template <class InputIt>
+    table(InputIt first, InputIt last, size_type bucket_count, allocator_type const& alloc)
+        : table(first, last, bucket_count, Hash(), KeyEqual(), alloc) {}
+
+    template <class InputIt>
+    table(InputIt first, InputIt last, size_type bucket_count, Hash const& hash, allocator_type const& alloc)
+        : table(first, last, bucket_count, hash, KeyEqual(), alloc) {}
+
+    table(table const& other)
+        : table(other, other.m_values.get_allocator()) {}
+
+    table(table const& other, allocator_type const& alloc)
+        : m_values(other.m_values, alloc)
+        , m_max_load_factor(other.m_max_load_factor)
+        , m_hash(other.m_hash)
+        , m_equal(other.m_equal) {
+        copy_buckets(other);
+    }
+
+    table(table&& other) noexcept
+        : table(std::move(other), other.m_values.get_allocator()) {}
+
+    table(table&& other, allocator_type const& alloc) noexcept
+        : m_values(std::move(other.m_values), alloc)
+        , m_buckets(std::exchange(other.m_buckets, nullptr))
+        , m_num_buckets(std::exchange(other.m_num_buckets, 0))
+        , m_max_bucket_capacity(std::exchange(other.m_max_bucket_capacity, 0))
+        , m_max_load_factor(std::exchange(other.m_max_load_factor, default_max_load_factor))
+        , m_hash(std::exchange(other.m_hash, {}))
+        , m_equal(std::exchange(other.m_equal, {}))
+        , m_shifts(std::exchange(other.m_shifts, initial_shifts)) {
+        other.m_values.clear();
+    }
+
+    table(std::initializer_list<value_type> ilist,
+          size_t bucket_count = 0,
+          Hash const& hash = Hash(),
+          KeyEqual const& equal = KeyEqual(),
+          allocator_type const& alloc = allocator_type())
+        : table(bucket_count, hash, equal, alloc) {
+        insert(ilist);
+    }
+
+    table(std::initializer_list<value_type> ilist, size_type bucket_count, allocator_type const& alloc)
+        : table(ilist, bucket_count, Hash(), KeyEqual(), alloc) {}
+
+    table(std::initializer_list<value_type> init, size_type bucket_count, Hash const& hash, allocator_type const& alloc)
+        : table(init, bucket_count, hash, KeyEqual(), alloc) {}
+
+    ~table() {
+        if (nullptr != m_buckets) {
+            auto ba = bucket_alloc(m_values.get_allocator());
+            bucket_alloc_traits::deallocate(ba, m_buckets, bucket_count());
+        }
+    }
+
+    auto operator=(table const& other) -> table& {
+        if (&other != this) {
+            deallocate_buckets(); // deallocate before m_values is set (might have another allocator)
+            m_values = other.m_values;
+            m_max_load_factor = other.m_max_load_factor;
+            m_hash = other.m_hash;
+            m_equal = other.m_equal;
+            m_shifts = initial_shifts;
+            copy_buckets(other);
+        }
+        return *this;
+    }
+
+    auto operator=(table&& other) noexcept(
+        noexcept(std::is_nothrow_move_assignable_v<value_container_type>&& std::is_nothrow_move_assignable_v<Hash>&&
+                     std::is_nothrow_move_assignable_v<KeyEqual>)) -> table& {
+        if (&other != this) {
+            deallocate_buckets(); // deallocate before m_values is set (might have another allocator)
+            m_values = std::move(other.m_values);
+            m_buckets = std::exchange(other.m_buckets, nullptr);
+            m_num_buckets = std::exchange(other.m_num_buckets, 0);
+            m_max_bucket_capacity = std::exchange(other.m_max_bucket_capacity, 0);
+            m_max_load_factor = std::exchange(other.m_max_load_factor, default_max_load_factor);
+            m_hash = std::exchange(other.m_hash, {});
+            m_equal = std::exchange(other.m_equal, {});
+            m_shifts = std::exchange(other.m_shifts, initial_shifts);
+            other.m_values.clear();
+        }
+        return *this;
+    }
+
+    auto operator=(std::initializer_list<value_type> ilist) -> table& {
+        clear();
+        insert(ilist);
+        return *this;
+    }
+
+    auto get_allocator() const noexcept -> allocator_type {
+        return m_values.get_allocator();
+    }
+
+    // iterators //////////////////////////////////////////////////////////////
+
+    auto begin() noexcept -> iterator {
+        return m_values.begin();
+    }
+
+    auto begin() const noexcept -> const_iterator {
+        return m_values.begin();
+    }
+
+    auto cbegin() const noexcept -> const_iterator {
+        return m_values.cbegin();
+    }
+
+    auto end() noexcept -> iterator {
+        return m_values.end();
+    }
+
+    auto cend() const noexcept -> const_iterator {
+        return m_values.cend();
+    }
+
+    auto end() const noexcept -> const_iterator {
+        return m_values.end();
+    }
+
+    // capacity ///////////////////////////////////////////////////////////////
+
+    [[nodiscard]] auto empty() const noexcept -> bool {
+        return m_values.empty();
+    }
+
+    [[nodiscard]] auto size() const noexcept -> size_t {
+        return m_values.size();
+    }
+
+    [[nodiscard]] static constexpr auto max_size() noexcept -> size_t {
+        if constexpr ((std::numeric_limits<value_idx_type>::max)() == (std::numeric_limits<size_t>::max)()) {
+            return size_t{1} << (sizeof(value_idx_type) * 8 - 1);
+        } else {
+            return size_t{1} << (sizeof(value_idx_type) * 8);
+        }
+    }
+
+    // modifiers //////////////////////////////////////////////////////////////
+
+    void clear() {
+        m_values.clear();
+        clear_buckets();
+    }
+
+    auto insert(value_type const& value) -> std::pair<iterator, bool> {
+        return emplace(value);
+    }
+
+    auto insert(value_type&& value) -> std::pair<iterator, bool> {
+        return emplace(std::move(value));
+    }
+
+    template <class P, std::enable_if_t<std::is_constructible_v<value_type, P&&>, bool> = true>
+    auto insert(P&& value) -> std::pair<iterator, bool> {
+        return emplace(std::forward<P>(value));
+    }
+
+    auto insert(const_iterator /*hint*/, value_type const& value) -> iterator {
+        return insert(value).first;
+    }
+
+    auto insert(const_iterator /*hint*/, value_type&& value) -> iterator {
+        return insert(std::move(value)).first;
+    }
+
+    template <class P, std::enable_if_t<std::is_constructible_v<value_type, P&&>, bool> = true>
+    auto insert(const_iterator /*hint*/, P&& value) -> iterator {
+        return insert(std::forward<P>(value)).first;
+    }
+
+    template <class InputIt>
+    void insert(InputIt first, InputIt last) {
+        while (first != last) {
+            insert(*first);
+            ++first;
+        }
+    }
+
+    void insert(std::initializer_list<value_type> ilist) {
+        insert(ilist.begin(), ilist.end());
+    }
+
+    // nonstandard API: *this is emptied.
+    // Also see "A Standard flat_map" https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p0429r9.pdf
+    auto extract() && -> value_container_type {
+        return std::move(m_values);
+    }
+
+    // nonstandard API:
+    // Discards the internally held container and replaces it with the one passed. Erases non-unique elements.
+    auto replace(value_container_type&& container) {
+        if (ANKERL_UNORDERED_DENSE_UNLIKELY(container.size() > max_size())) {
+            on_error_too_many_elements();
+        }
+        auto shifts = calc_shifts_for_size(container.size());
+        if (0 == m_num_buckets || shifts < m_shifts || container.get_allocator() != m_values.get_allocator()) {
+            m_shifts = shifts;
+            deallocate_buckets();
+            allocate_buckets_from_shift();
+        }
+        clear_buckets();
+
+        m_values = std::move(container);
+
+        // can't use clear_and_fill_buckets_from_values() because container elements might not be unique
+        auto value_idx = value_idx_type{};
+
+        // loop until we reach the end of the container. duplicated entries will be replaced with back().
+        while (value_idx != static_cast<value_idx_type>(m_values.size())) {
+            auto const& key = get_key(m_values[value_idx]);
+
+            auto hash = mixed_hash(key);
+            auto dist_and_fingerprint = dist_and_fingerprint_from_hash(hash);
+            auto bucket_idx = bucket_idx_from_hash(hash);
+
+            bool key_found = false;
+            while (true) {
+                auto const& bucket = at(m_buckets, bucket_idx);
+                if (dist_and_fingerprint > bucket.m_dist_and_fingerprint) {
+                    break;
+                }
+                if (dist_and_fingerprint == bucket.m_dist_and_fingerprint &&
+                    m_equal(key, m_values[bucket.m_value_idx].first)) {
+                    key_found = true;
+                    break;
+                }
+                dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+                bucket_idx = next(bucket_idx);
+            }
+
+            if (key_found) {
+                if (value_idx != static_cast<value_idx_type>(m_values.size() - 1)) {
+                    m_values[value_idx] = std::move(m_values.back());
+                }
+                m_values.pop_back();
+            } else {
+                place_and_shift_up({dist_and_fingerprint, value_idx}, bucket_idx);
+                ++value_idx;
+            }
+        }
+    }
+
+    template <class M, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto insert_or_assign(Key const& key, M&& mapped) -> std::pair<iterator, bool> {
+        return do_insert_or_assign(key, std::forward<M>(mapped));
+    }
+
+    template <class M, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto insert_or_assign(Key&& key, M&& mapped) -> std::pair<iterator, bool> {
+        return do_insert_or_assign(std::move(key), std::forward<M>(mapped));
+    }
+
+    template <typename K,
+              typename M,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto insert_or_assign(K&& key, M&& mapped) -> std::pair<iterator, bool> {
+        return do_insert_or_assign(std::forward<K>(key), std::forward<M>(mapped));
+    }
+
+    template <class M, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto insert_or_assign(const_iterator /*hint*/, Key const& key, M&& mapped) -> iterator {
+        return do_insert_or_assign(key, std::forward<M>(mapped)).first;
+    }
+
+    template <class M, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto insert_or_assign(const_iterator /*hint*/, Key&& key, M&& mapped) -> iterator {
+        return do_insert_or_assign(std::move(key), std::forward<M>(mapped)).first;
+    }
+
+    template <typename K,
+              typename M,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto insert_or_assign(const_iterator /*hint*/, K&& key, M&& mapped) -> iterator {
+        return do_insert_or_assign(std::forward<K>(key), std::forward<M>(mapped)).first;
+    }
+
+    // Single arguments for unordered_set can be used without having to construct the value_type
+    template <class K,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<!is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto emplace(K&& key) -> std::pair<iterator, bool> {
+        if (is_full()) {
+            increase_size();
+        }
+
+        auto hash = mixed_hash(key);
+        auto dist_and_fingerprint = dist_and_fingerprint_from_hash(hash);
+        auto bucket_idx = bucket_idx_from_hash(hash);
+
+        while (dist_and_fingerprint <= at(m_buckets, bucket_idx).m_dist_and_fingerprint) {
+            if (dist_and_fingerprint == at(m_buckets, bucket_idx).m_dist_and_fingerprint &&
+                m_equal(key, m_values[at(m_buckets, bucket_idx).m_value_idx])) {
+                // found it, return without ever actually creating anything
+                return {begin() + static_cast<difference_type>(at(m_buckets, bucket_idx).m_value_idx), false};
+            }
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+        }
+
+        // value is new, insert element first, so when exception happens we are in a valid state
+        m_values.emplace_back(std::forward<K>(key));
+        // now place the bucket and shift up until we find an empty spot
+        auto value_idx = static_cast<value_idx_type>(m_values.size() - 1);
+        place_and_shift_up({dist_and_fingerprint, value_idx}, bucket_idx);
+        return {begin() + static_cast<difference_type>(value_idx), true};
+    }
+
+    template <class... Args>
+    auto emplace(Args&&... args) -> std::pair<iterator, bool> {
+        if (is_full()) {
+            increase_size();
+        }
+
+        // we have to instantiate the value_type to be able to access the key.
+        // 1. emplace_back the object so it is constructed. 2. If the key is already there, pop it later in the loop.
+        auto& key = get_key(m_values.emplace_back(std::forward<Args>(args)...));
+        auto hash = mixed_hash(key);
+        auto dist_and_fingerprint = dist_and_fingerprint_from_hash(hash);
+        auto bucket_idx = bucket_idx_from_hash(hash);
+
+        while (dist_and_fingerprint <= at(m_buckets, bucket_idx).m_dist_and_fingerprint) {
+            if (dist_and_fingerprint == at(m_buckets, bucket_idx).m_dist_and_fingerprint &&
+                m_equal(key, get_key(m_values[at(m_buckets, bucket_idx).m_value_idx]))) {
+                m_values.pop_back(); // value was already there, so get rid of it
+                return {begin() + static_cast<difference_type>(at(m_buckets, bucket_idx).m_value_idx), false};
+            }
+            dist_and_fingerprint = dist_inc(dist_and_fingerprint);
+            bucket_idx = next(bucket_idx);
+        }
+
+        // value is new, place the bucket and shift up until we find an empty spot
+        auto value_idx = static_cast<value_idx_type>(m_values.size() - 1);
+        place_and_shift_up({dist_and_fingerprint, value_idx}, bucket_idx);
+
+        return {begin() + static_cast<difference_type>(value_idx), true};
+    }
+
+    template <class... Args>
+    auto emplace_hint(const_iterator /*hint*/, Args&&... args) -> iterator {
+        return emplace(std::forward<Args>(args)...).first;
+    }
+
+    template <class... Args, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto try_emplace(Key const& key, Args&&... args) -> std::pair<iterator, bool> {
+        return do_try_emplace(key, std::forward<Args>(args)...);
+    }
+
+    template <class... Args, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto try_emplace(Key&& key, Args&&... args) -> std::pair<iterator, bool> {
+        return do_try_emplace(std::move(key), std::forward<Args>(args)...);
+    }
+
+    template <class... Args, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto try_emplace(const_iterator /*hint*/, Key const& key, Args&&... args) -> iterator {
+        return do_try_emplace(key, std::forward<Args>(args)...).first;
+    }
+
+    template <class... Args, typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto try_emplace(const_iterator /*hint*/, Key&& key, Args&&... args) -> iterator {
+        return do_try_emplace(std::move(key), std::forward<Args>(args)...).first;
+    }
+
+    template <
+        typename K,
+        typename... Args,
+        typename Q = T,
+        typename H = Hash,
+        typename KE = KeyEqual,
+        std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE> && is_neither_convertible_v<K&&, iterator, const_iterator>,
+                         bool> = true>
+    auto try_emplace(K&& key, Args&&... args) -> std::pair<iterator, bool> {
+        return do_try_emplace(std::forward<K>(key), std::forward<Args>(args)...);
+    }
+
+    template <
+        typename K,
+        typename... Args,
+        typename Q = T,
+        typename H = Hash,
+        typename KE = KeyEqual,
+        std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE> && is_neither_convertible_v<K&&, iterator, const_iterator>,
+                         bool> = true>
+    auto try_emplace(const_iterator /*hint*/, K&& key, Args&&... args) -> iterator {
+        return do_try_emplace(std::forward<K>(key), std::forward<Args>(args)...).first;
+    }
+
+    auto erase(iterator it) -> iterator {
+        auto hash = mixed_hash(get_key(*it));
+        auto bucket_idx = bucket_idx_from_hash(hash);
+
+        auto const value_idx_to_remove = static_cast<value_idx_type>(it - cbegin());
+        while (at(m_buckets, bucket_idx).m_value_idx != value_idx_to_remove) {
+            bucket_idx = next(bucket_idx);
+        }
+
+        do_erase(bucket_idx);
+        return begin() + static_cast<difference_type>(value_idx_to_remove);
+    }
+
+    template <typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto erase(const_iterator it) -> iterator {
+        return erase(begin() + (it - cbegin()));
+    }
+
+    auto erase(const_iterator first, const_iterator last) -> iterator {
+        auto const idx_first = first - cbegin();
+        auto const idx_last = last - cbegin();
+        auto const first_to_last = std::distance(first, last);
+        auto const last_to_end = std::distance(last, cend());
+
+        // remove elements from left to right which moves elements from the end back
+        auto const mid = idx_first + (std::min)(first_to_last, last_to_end);
+        auto idx = idx_first;
+        while (idx != mid) {
+            erase(begin() + idx);
+            ++idx;
+        }
+
+        // all elements from the right are moved, now remove the last element until all done
+        idx = idx_last;
+        while (idx != mid) {
+            --idx;
+            erase(begin() + idx);
+        }
+
+        return begin() + idx_first;
+    }
+
+    auto erase(Key const& key) -> size_t {
+        return do_erase_key(key);
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto erase(K&& key) -> size_t {
+        return do_erase_key(std::forward<K>(key));
+    }
+
+    void swap(table& other) noexcept(noexcept(std::is_nothrow_swappable_v<value_container_type>&&
+                                                  std::is_nothrow_swappable_v<Hash>&& std::is_nothrow_swappable_v<KeyEqual>)) {
+        using std::swap;
+        swap(other, *this);
+    }
+
+    // lookup /////////////////////////////////////////////////////////////////
+
+    template <typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto at(key_type const& key) -> Q& {
+        return do_at(key);
+    }
+
+    template <typename K,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto at(K const& key) -> Q& {
+        return do_at(key);
+    }
+
+    template <typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto at(key_type const& key) const -> Q const& {
+        return do_at(key);
+    }
+
+    template <typename K,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto at(K const& key) const -> Q const& {
+        return do_at(key);
+    }
+
+    template <typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto operator[](Key const& key) -> Q& {
+        return try_emplace(key).first->second;
+    }
+
+    template <typename Q = T, std::enable_if_t<is_map_v<Q>, bool> = true>
+    auto operator[](Key&& key) -> Q& {
+        return try_emplace(std::move(key)).first->second;
+    }
+
+    template <typename K,
+              typename Q = T,
+              typename H = Hash,
+              typename KE = KeyEqual,
+              std::enable_if_t<is_map_v<Q> && is_transparent_v<H, KE>, bool> = true>
+    auto operator[](K&& key) -> Q& {
+        return try_emplace(std::forward<K>(key)).first->second;
+    }
+
+    auto count(Key const& key) const -> size_t {
+        return find(key) == end() ? 0 : 1;
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto count(K const& key) const -> size_t {
+        return find(key) == end() ? 0 : 1;
+    }
+
+    auto find(Key const& key) -> iterator {
+        return do_find(key);
+    }
+
+    auto find(Key const& key) const -> const_iterator {
+        return do_find(key);
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto find(K const& key) -> iterator {
+        return do_find(key);
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto find(K const& key) const -> const_iterator {
+        return do_find(key);
+    }
+
+    auto contains(Key const& key) const -> bool {
+        return find(key) != end();
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto contains(K const& key) const -> bool {
+        return find(key) != end();
+    }
+
+    auto equal_range(Key const& key) -> std::pair<iterator, iterator> {
+        auto it = do_find(key);
+        return {it, it == end() ? end() : it + 1};
+    }
+
+    auto equal_range(const Key& key) const -> std::pair<const_iterator, const_iterator> {
+        auto it = do_find(key);
+        return {it, it == end() ? end() : it + 1};
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto equal_range(K const& key) -> std::pair<iterator, iterator> {
+        auto it = do_find(key);
+        return {it, it == end() ? end() : it + 1};
+    }
+
+    template <class K, class H = Hash, class KE = KeyEqual, std::enable_if_t<is_transparent_v<H, KE>, bool> = true>
+    auto equal_range(K const& key) const -> std::pair<const_iterator, const_iterator> {
+        auto it = do_find(key);
+        return {it, it == end() ? end() : it + 1};
+    }
+
+    // bucket interface ///////////////////////////////////////////////////////
+
+    auto bucket_count() const noexcept -> size_t { // NOLINT(modernize-use-nodiscard)
+        return m_num_buckets;
+    }
+
+    static constexpr auto max_bucket_count() noexcept -> size_t { // NOLINT(modernize-use-nodiscard)
+        return max_size();
+    }
+
+    // hash policy ////////////////////////////////////////////////////////////
+
+    [[nodiscard]] auto load_factor() const -> float {
+        return bucket_count() ? static_cast<float>(size()) / static_cast<float>(bucket_count()) : 0.0F;
+    }
+
+    [[nodiscard]] auto max_load_factor() const -> float {
+        return m_max_load_factor;
+    }
+
+    void max_load_factor(float ml) {
+        m_max_load_factor = ml;
+        if (m_num_buckets != max_bucket_count()) {
+            m_max_bucket_capacity = static_cast<value_idx_type>(static_cast<float>(bucket_count()) * max_load_factor());
+        }
+    }
+
+    void rehash(size_t count) {
+        count = (std::min)(count, max_size());
+        auto shifts = calc_shifts_for_size((std::max)(count, size()));
+        if (shifts != m_shifts) {
+            m_shifts = shifts;
+            deallocate_buckets();
+            m_values.shrink_to_fit();
+            allocate_buckets_from_shift();
+            clear_and_fill_buckets_from_values();
+        }
+    }
+
+    void reserve(size_t capa) {
+        capa = (std::min)(capa, max_size());
+        if constexpr (has_reserve<value_container_type>) {
+            // std::deque doesn't have reserve(). Make sure we only call when available
+            m_values.reserve(capa);
+        }
+        auto shifts = calc_shifts_for_size((std::max)(capa, size()));
+        if (0 == m_num_buckets || shifts < m_shifts) {
+            m_shifts = shifts;
+            deallocate_buckets();
+            allocate_buckets_from_shift();
+            clear_and_fill_buckets_from_values();
+        }
+    }
+
+    // observers //////////////////////////////////////////////////////////////
+
+    auto hash_function() const -> hasher {
+        return m_hash;
+    }
+
+    auto key_eq() const -> key_equal {
+        return m_equal;
+    }
+
+    // nonstandard API: expose the underlying values container
+    [[nodiscard]] auto values() const noexcept -> value_container_type const& {
+        return m_values;
+    }
+
+    // non-member functions ///////////////////////////////////////////////////
+
+    friend auto operator==(table const& a, table const& b) -> bool {
+        if (&a == &b) {
+            return true;
+        }
+        if (a.size() != b.size()) {
+            return false;
+        }
+        for (auto const& b_entry : b) {
+            auto it = a.find(get_key(b_entry));
+            if constexpr (is_map_v<T>) {
+                // map: check that key is here, then also check that value is the same
+                if (a.end() == it || !(b_entry.second == it->second)) {
+                    return false;
+                }
+            } else {
+                // set: only check that the key is here
+                if (a.end() == it) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+    friend auto operator!=(table const& a, table const& b) -> bool {
+        return !(a == b);
+    }
+};
+
+} // namespace detail
+
+template <class Key,
+          class T,
+          class Hash = hash<Key>,
+          class KeyEqual = std::equal_to<Key>,
+          class AllocatorOrContainer = std::allocator<std::pair<Key, T>>,
+          class Bucket = bucket_type::standard>
+using map = detail::table<Key, T, Hash, KeyEqual, AllocatorOrContainer, Bucket>;
+
+template <class Key,
+          class Hash = hash<Key>,
+          class KeyEqual = std::equal_to<Key>,
+          class AllocatorOrContainer = std::allocator<Key>,
+          class Bucket = bucket_type::standard>
+using set = detail::table<Key, void, Hash, KeyEqual, AllocatorOrContainer, Bucket>;
+
+#    if ANKERL_UNORDERED_DENSE_PMR
+
+namespace pmr {
+
+template <class Key,
+          class T,
+          class Hash = hash<Key>,
+          class KeyEqual = std::equal_to<Key>,
+          class Bucket = bucket_type::standard>
+using map = detail::table<Key, T, Hash, KeyEqual, ANKERL_UNORDERED_DENSE_PMR_ALLOCATOR<std::pair<Key, T>>, Bucket>;
+
+template <class Key, class Hash = hash<Key>, class KeyEqual = std::equal_to<Key>, class Bucket = bucket_type::standard>
+using set = detail::table<Key, void, Hash, KeyEqual, ANKERL_UNORDERED_DENSE_PMR_ALLOCATOR<Key>, Bucket>;
+
+} // namespace pmr
+
+#    endif
+
+// deduction guides ///////////////////////////////////////////////////////////
+
+// deduction guides for alias templates are only possible since C++20
+// see https://en.cppreference.com/w/cpp/language/class_template_argument_deduction
+
+} // namespace ANKERL_UNORDERED_DENSE_NAMESPACE
+} // namespace ankerl::unordered_dense
+
+// std extensions /////////////////////////////////////////////////////////////
+
+namespace std { // NOLINT(cert-dcl58-cpp)
+
+template <class Key, class T, class Hash, class KeyEqual, class AllocatorOrContainer, class Bucket, class Pred>
+// NOLINTNEXTLINE(cert-dcl58-cpp)
+auto erase_if(ankerl::unordered_dense::detail::table<Key, T, Hash, KeyEqual, AllocatorOrContainer, Bucket>& map, Pred pred)
+    -> size_t {
+    using map_t = ankerl::unordered_dense::detail::table<Key, T, Hash, KeyEqual, AllocatorOrContainer, Bucket>;
+
+    // going back to front because erase() invalidates the end iterator
+    auto const old_size = map.size();
+    auto idx = old_size;
+    while (idx) {
+        --idx;
+        auto it = map.begin() + static_cast<typename map_t::difference_type>(idx);
+        if (pred(*it)) {
+            map.erase(it);
+        }
+    }
+
+    return map.size() - old_size;
+}
+
+} // namespace std
+
+#endif
+#endif
diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
index 04dede546..5d085f333 100644
--- a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
+++ b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
@@ -5,14 +5,11 @@
 
 #include <stack>
 #include <functional>
-#include <unordered_set>
 #include <sstream>
 #include <queue>
 #include <functional>
 #include <boost/log/trivial.hpp>
 
-#include "utils/VoronoiUtils.hpp"
-
 #include "utils/linearAlg2D.hpp"
 #include "Utils.hpp"
 #include "SVG.hpp"
@@ -20,27 +17,10 @@
 #include "Geometry/VoronoiUtilsCgal.hpp"
 #include "../EdgeGrid.hpp"
 
+#include "Geometry/VoronoiUtils.hpp"
+
 #define SKELETAL_TRAPEZOIDATION_BEAD_SEARCH_MAX 1000 //A limit to how long it'll keep searching for adjacent beads. Increasing will re-use beadings more often (saving performance), but search longer for beading (costing performance).
 
-namespace boost::polygon {
-
-template<> struct geometry_concept<Slic3r::Arachne::PolygonsSegmentIndex>
-{
-    typedef segment_concept type;
-};
-
-template<> struct segment_traits<Slic3r::Arachne::PolygonsSegmentIndex>
-{
-    typedef coord_t          coordinate_type;
-    typedef Slic3r::Point    point_type;
-    static inline point_type get(const Slic3r::Arachne::PolygonsSegmentIndex &CSegment, direction_1d dir)
-    {
-        return dir.to_int() ? CSegment.p() : CSegment.next().p();
-    }
-};
-
-} // namespace boost::polygon
-
 namespace Slic3r::Arachne
 {
 
@@ -109,8 +89,7 @@ static void export_graph_to_svg(const std::string
 }
 #endif
 
-SkeletalTrapezoidation::node_t& SkeletalTrapezoidation::makeNode(vd_t::vertex_type& vd_node, Point p)
-{
+SkeletalTrapezoidation::node_t &SkeletalTrapezoidation::makeNode(const VD::vertex_type &vd_node, Point p) {
     auto he_node_it = vd_node_to_he_node.find(&vd_node);
     if (he_node_it == vd_node_to_he_node.end())
     {
@@ -125,8 +104,7 @@ SkeletalTrapezoidation::node_t& SkeletalTrapezoidation::makeNode(vd_t::vertex_ty
     }
 }
 
-void SkeletalTrapezoidation::transferEdge(Point from, Point to, vd_t::edge_type& vd_edge, edge_t*& prev_edge, Point& start_source_point, Point& end_source_point, const std::vector<Segment>& segments)
-{
+void SkeletalTrapezoidation::transferEdge(Point from, Point to, const VD::edge_type &vd_edge, edge_t *&prev_edge, Point &start_source_point, Point &end_source_point, const std::vector<Segment> &segments) {
     auto he_edge_it = vd_edge_to_he_edge.find(vd_edge.twin());
     if (he_edge_it != vd_edge_to_he_edge.end())
     { // Twin segment(s) have already been made
@@ -181,7 +159,7 @@ void SkeletalTrapezoidation::transferEdge(Point from, Point to, vd_t::edge_type&
     }
     else
     {
-        std::vector<Point> discretized = discretize(vd_edge, segments);
+        Points discretized = discretize(vd_edge, segments);
         assert(discretized.size() >= 2);
         if(discretized.size() < 2)
         {
@@ -236,45 +214,42 @@ void SkeletalTrapezoidation::transferEdge(Point from, Point to, vd_t::edge_type&
     }
 }
 
-std::vector<Point> SkeletalTrapezoidation::discretize(const vd_t::edge_type& vd_edge, const std::vector<Segment>& segments)
+Points SkeletalTrapezoidation::discretize(const VD::edge_type& vd_edge, const std::vector<Segment>& segments)
 {
+    assert(Geometry::VoronoiUtils::is_in_range<coord_t>(vd_edge));
+
     /*Terminology in this function assumes that the edge moves horizontally from
     left to right. This is not necessarily the case; the edge can go in any
     direction, but it helps to picture it in a certain direction in your head.*/
 
-    const vd_t::cell_type* left_cell = vd_edge.cell();
-    const vd_t::cell_type* right_cell = vd_edge.twin()->cell();
+    const VD::cell_type *left_cell  = vd_edge.cell();
+    const VD::cell_type *right_cell = vd_edge.twin()->cell();
 
-    assert(VoronoiUtils::p(vd_edge.vertex0()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge.vertex0()).x() >= std::numeric_limits<coord_t>::lowest());
-    assert(VoronoiUtils::p(vd_edge.vertex0()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge.vertex0()).y() >= std::numeric_limits<coord_t>::lowest());
-    assert(VoronoiUtils::p(vd_edge.vertex1()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge.vertex1()).x() >= std::numeric_limits<coord_t>::lowest());
-    assert(VoronoiUtils::p(vd_edge.vertex1()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge.vertex1()).y() >= std::numeric_limits<coord_t>::lowest());
+    Point start = Geometry::VoronoiUtils::to_point(vd_edge.vertex0()).cast<coord_t>();
+    Point end   = Geometry::VoronoiUtils::to_point(vd_edge.vertex1()).cast<coord_t>();
 
-    Point start = VoronoiUtils::p(vd_edge.vertex0()).cast<coord_t>();
-    Point end = VoronoiUtils::p(vd_edge.vertex1()).cast<coord_t>();
-    
     bool point_left = left_cell->contains_point();
     bool point_right = right_cell->contains_point();
     if ((!point_left && !point_right) || vd_edge.is_secondary()) // Source vert is directly connected to source segment
     {
-        return std::vector<Point>({ start, end });
+        return Points({ start, end });
     }
     else if (point_left != point_right) //This is a parabolic edge between a point and a line.
     {
-        Point p = VoronoiUtils::getSourcePoint(*(point_left ? left_cell : right_cell), segments);
-        const Segment& s = VoronoiUtils::getSourceSegment(*(point_left ? right_cell : left_cell), segments);
-        return VoronoiUtils::discretizeParabola(p, s, start, end, discretization_step_size, transitioning_angle);
+        Point          p = Geometry::VoronoiUtils::get_source_point(*(point_left ? left_cell : right_cell), segments.begin(), segments.end());
+        const Segment& s = Geometry::VoronoiUtils::get_source_segment(*(point_left ? right_cell : left_cell), segments.begin(), segments.end());
+        return Geometry::VoronoiUtils::discretize_parabola(p, s, start, end, discretization_step_size, transitioning_angle);
     }
     else //This is a straight edge between two points.
     {
         /*While the edge is straight, it is still discretized since the part
         becomes narrower between the two points. As such it may need different
         beadings along the way.*/
-        Point left_point = VoronoiUtils::getSourcePoint(*left_cell, segments);
-        Point right_point = VoronoiUtils::getSourcePoint(*right_cell, segments);
-        coord_t d = (right_point - left_point).cast<int64_t>().norm();
-        Point middle = (left_point + right_point) / 2;
-        Point x_axis_dir = Point(right_point - left_point).rotate_90_degree_ccw();
+        Point   left_point    = Geometry::VoronoiUtils::get_source_point(*left_cell, segments.begin(), segments.end());
+        Point   right_point   = Geometry::VoronoiUtils::get_source_point(*right_cell, segments.begin(), segments.end());
+        coord_t d             = (right_point - left_point).cast<int64_t>().norm();
+        Point   middle        = (left_point + right_point) / 2;
+        Point   x_axis_dir    = perp(Point(right_point - left_point));
         coord_t x_axis_length = x_axis_dir.cast<int64_t>().norm();
 
         const auto projected_x = [x_axis_dir, x_axis_length, middle](Point from) //Project a point on the edge.
@@ -311,7 +286,7 @@ std::vector<Point> SkeletalTrapezoidation::discretize(const vd_t::edge_type& vd_
         //Start generating points along the edge.
         Point a = start;
         Point b = end;
-        std::vector<Point> ret;
+        Points ret;
         ret.emplace_back(a);
 
         //Introduce an extra edge at the borders of the markings?
@@ -351,8 +326,7 @@ std::vector<Point> SkeletalTrapezoidation::discretize(const vd_t::edge_type& vd_
     }
 }
 
-bool SkeletalTrapezoidation::computePointCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::edge_type*& ending_vd_edge, const std::vector<Segment>& segments)
-{
+bool SkeletalTrapezoidation::computePointCellRange(const VD::cell_type &cell, Point &start_source_point, Point &end_source_point, const VD::edge_type *&starting_vd_edge, const VD::edge_type *&ending_vd_edge, const std::vector<Segment> &segments) {
     if (cell.incident_edge()->is_infinite())
         return false; //Infinite edges only occur outside of the polygon. Don't copy any part of this cell.
 
@@ -360,16 +334,16 @@ bool SkeletalTrapezoidation::computePointCellRange(vd_t::cell_type& cell, Point&
     // Copy whole cell into graph or not at all
 
     // If the cell.incident_edge()->vertex0() is far away so much that it doesn't even fit into Vec2i64, then there is no way that it will be inside the input polygon.
-    if (const vd_t::vertex_type &vert = *cell.incident_edge()->vertex0();
+    if (const VD::vertex_type &vert = *cell.incident_edge()->vertex0();
         vert.x() >= double(std::numeric_limits<int64_t>::max()) || vert.x() <= double(std::numeric_limits<int64_t>::lowest()) ||
         vert.y() >= double(std::numeric_limits<int64_t>::max()) || vert.y() <= double(std::numeric_limits<int64_t>::lowest()))
         return false; // Don't copy any part of this cell
 
-    const Point source_point = VoronoiUtils::getSourcePoint(cell, segments);
-    const PolygonsPointIndex source_point_index = VoronoiUtils::getSourcePointIndex(cell, segments);
-    Vec2i64 some_point = VoronoiUtils::p(cell.incident_edge()->vertex0());
+    const Point              source_point       = Geometry::VoronoiUtils::get_source_point(cell, segments.begin(), segments.end());
+    const PolygonsPointIndex source_point_index = Geometry::VoronoiUtils::get_source_point_index(cell, segments.begin(), segments.end());
+    Vec2i64                  some_point         = Geometry::VoronoiUtils::to_point(cell.incident_edge()->vertex0());
     if (some_point == source_point.cast<int64_t>())
-        some_point = VoronoiUtils::p(cell.incident_edge()->vertex1());
+        some_point = Geometry::VoronoiUtils::to_point(cell.incident_edge()->vertex1());
 
     //Test if the some_point is even inside the polygon.
     //The edge leading out of a polygon must have an endpoint that's not in the corner following the contour of the polygon at that vertex.
@@ -378,16 +352,16 @@ bool SkeletalTrapezoidation::computePointCellRange(vd_t::cell_type& cell, Point&
     if (!LinearAlg2D::isInsideCorner(source_point_index.prev().p(), source_point_index.p(), source_point_index.next().p(), some_point))
         return false; // Don't copy any part of this cell
 
-    vd_t::edge_type* vd_edge = cell.incident_edge();
+    const VD::edge_type* vd_edge = cell.incident_edge();
     do {
         assert(vd_edge->is_finite());
-        if (Vec2i64 p1 = VoronoiUtils::p(vd_edge->vertex1()); p1 == source_point.cast<int64_t>()) {
+        if (Vec2i64 p1 = Geometry::VoronoiUtils::to_point(vd_edge->vertex1()); p1 == source_point.cast<int64_t>()) {
             start_source_point = source_point;
             end_source_point = source_point;
             starting_vd_edge = vd_edge->next();
             ending_vd_edge = vd_edge;
         } else {
-            assert((VoronoiUtils::p(vd_edge->vertex0()) == source_point.cast<int64_t>() || !vd_edge->is_secondary()) && "point cells must end in the point! They cannot cross the point with an edge, because collinear edges are not allowed in the input.");
+            assert((Geometry::VoronoiUtils::to_point(vd_edge->vertex0()) == source_point.cast<int64_t>() || !vd_edge->is_secondary()) && "point cells must end in the point! They cannot cross the point with an edge, because collinear edges are not allowed in the input.");
         }
     }
     while (vd_edge = vd_edge->next(), vd_edge != cell.incident_edge());
@@ -396,47 +370,6 @@ bool SkeletalTrapezoidation::computePointCellRange(vd_t::cell_type& cell, Point&
     return true;
 }
 
-void SkeletalTrapezoidation::computeSegmentCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::edge_type*& ending_vd_edge, const std::vector<Segment>& segments)
-{
-    const Segment &source_segment = VoronoiUtils::getSourceSegment(cell, segments);
-    const Point    from           = source_segment.from();
-    const Point    to             = source_segment.to();
-
-    // Find starting edge
-    // Find end edge
-    bool seen_possible_start = false;
-    bool after_start = false;
-    bool ending_edge_is_set_before_start = false;
-    vd_t::edge_type* edge = cell.incident_edge();
-    do {
-        if (edge->is_infinite())
-            continue;
-
-        Vec2i64 v0 = VoronoiUtils::p(edge->vertex0());
-        Vec2i64 v1 = VoronoiUtils::p(edge->vertex1());
-
-        assert(!(v0 == to.cast<int64_t>() && v1 == from.cast<int64_t>() ));
-        if (v0 == to.cast<int64_t>()  && !after_start) { // Use the last edge which starts in source_segment.to
-            starting_vd_edge = edge;
-            seen_possible_start = true;
-        }
-        else if (seen_possible_start) {
-            after_start = true;
-        }
-
-        if (v1 == from.cast<int64_t>()  && (!ending_vd_edge || ending_edge_is_set_before_start)) {
-            ending_edge_is_set_before_start = !after_start;
-            ending_vd_edge = edge;
-        }
-    } while (edge = edge->next(), edge != cell.incident_edge());
-    
-    assert(starting_vd_edge && ending_vd_edge);
-    assert(starting_vd_edge != ending_vd_edge);
-    
-    start_source_point = source_segment.to();
-    end_source_point = source_segment.from();
-}
-
 SkeletalTrapezoidation::SkeletalTrapezoidation(const Polygons& polys, const BeadingStrategy& beading_strategy,
                                                double transitioning_angle, coord_t discretization_step_size,
                                                coord_t transition_filter_dist, coord_t allowed_filter_deviation,
@@ -451,128 +384,6 @@ SkeletalTrapezoidation::SkeletalTrapezoidation(const Polygons& polys, const Bead
     constructFromPolygons(polys);
 }
 
-static bool has_finite_edge_with_non_finite_vertex(const Geometry::VoronoiDiagram &voronoi_diagram)
-{
-    for (const VoronoiUtils::vd_t::edge_type &edge : voronoi_diagram.edges()) {
-        if (edge.is_finite()) {
-            assert(edge.vertex0() != nullptr && edge.vertex1() != nullptr);
-            if (edge.vertex0() == nullptr || edge.vertex1() == nullptr || !VoronoiUtils::is_finite(*edge.vertex0()) ||
-                !VoronoiUtils::is_finite(*edge.vertex1()))
-                return true;
-        }
-    }
-    return false;
-}
-
-static bool detect_missing_voronoi_vertex(const Geometry::VoronoiDiagram &voronoi_diagram, const std::vector<SkeletalTrapezoidation::Segment> &segments) {
-    if (has_finite_edge_with_non_finite_vertex(voronoi_diagram))
-        return true;
-
-    for (VoronoiUtils::vd_t::cell_type cell : voronoi_diagram.cells()) {
-        if (!cell.incident_edge())
-            continue; // There is no spoon
-
-        if (cell.contains_segment()) {
-            const SkeletalTrapezoidation::Segment &source_segment = VoronoiUtils::getSourceSegment(cell, segments);
-            const Point                            from           = source_segment.from();
-            const Point                            to             = source_segment.to();
-
-            // Find starting edge
-            // Find end edge
-            bool                           seen_possible_start             = false;
-            bool                           after_start                     = false;
-            bool                           ending_edge_is_set_before_start = false;
-            VoronoiUtils::vd_t::edge_type *starting_vd_edge                = nullptr;
-            VoronoiUtils::vd_t::edge_type *ending_vd_edge                  = nullptr;
-            VoronoiUtils::vd_t::edge_type *edge                            = cell.incident_edge();
-            do {
-                if (edge->is_infinite() || edge->vertex0() == nullptr || edge->vertex1() == nullptr || !VoronoiUtils::is_finite(*edge->vertex0()) || !VoronoiUtils::is_finite(*edge->vertex1()))
-                    continue;
-
-                Vec2i64 v0 = VoronoiUtils::p(edge->vertex0());
-                Vec2i64 v1 = VoronoiUtils::p(edge->vertex1());
-
-                assert(!(v0 == to.cast<int64_t>() && v1 == from.cast<int64_t>()));
-                if (v0 == to.cast<int64_t>() && !after_start) { // Use the last edge which starts in source_segment.to
-                    starting_vd_edge    = edge;
-                    seen_possible_start = true;
-                } else if (seen_possible_start) {
-                    after_start = true;
-                }
-
-                if (v1 == from.cast<int64_t>() && (!ending_vd_edge || ending_edge_is_set_before_start)) {
-                    ending_edge_is_set_before_start = !after_start;
-                    ending_vd_edge                  = edge;
-                }
-            } while (edge = edge->next(), edge != cell.incident_edge());
-
-            if (!starting_vd_edge || !ending_vd_edge || starting_vd_edge == ending_vd_edge)
-                return true;
-        }
-    }
-
-    return false;
-}
-
-static bool has_missing_twin_edge(const SkeletalTrapezoidationGraph &graph)
-{
-    for (const auto &edge : graph.edges)
-        if (edge.twin == nullptr)
-            return true;
-    return false;
-}
-
-inline static std::unordered_map<Point, Point, PointHash> try_to_fix_degenerated_voronoi_diagram_by_rotation(
-    Geometry::VoronoiDiagram                     &voronoi_diagram,
-    const Polygons                               &polys,
-    Polygons                                     &polys_rotated,
-    std::vector<SkeletalTrapezoidation::Segment> &segments,
-    const double                                  fix_angle)
-{
-    std::unordered_map<Point, Point, PointHash> vertex_mapping;
-    for (Polygon &poly : polys_rotated)
-        poly.rotate(fix_angle);
-
-    assert(polys_rotated.size() == polys.size());
-    for (size_t poly_idx = 0; poly_idx < polys.size(); ++poly_idx) {
-        assert(polys_rotated[poly_idx].size() == polys[poly_idx].size());
-        for (size_t point_idx = 0; point_idx < polys[poly_idx].size(); ++point_idx)
-            vertex_mapping.insert({polys_rotated[poly_idx][point_idx], polys[poly_idx][point_idx]});
-    }
-
-    segments.clear();
-    for (size_t poly_idx = 0; poly_idx < polys_rotated.size(); poly_idx++)
-        for (size_t point_idx = 0; point_idx < polys_rotated[poly_idx].size(); point_idx++)
-            segments.emplace_back(&polys_rotated, poly_idx, point_idx);
-
-    voronoi_diagram.clear();
-    construct_voronoi(segments.begin(), segments.end(), &voronoi_diagram);
-
-#ifdef ARACHNE_DEBUG_VORONOI
-    {
-        static int iRun = 0;
-        dump_voronoi_to_svg(debug_out_path("arachne_voronoi-diagram-rotated-%d.svg", iRun++).c_str(), voronoi_diagram, to_points(polys), to_lines(polys));
-    }
-#endif
-
-    assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(voronoi_diagram));
-
-    return vertex_mapping;
-}
-
-inline static void rotate_back_skeletal_trapezoidation_graph_after_fix(SkeletalTrapezoidationGraph                       &graph,
-                                                                       const double                                       fix_angle,
-                                                                       const std::unordered_map<Point, Point, PointHash> &vertex_mapping)
-{
-    for (STHalfEdgeNode &node : graph.nodes) {
-        // If a mapping exists between a rotated point and an original point, use this mapping. Otherwise, rotate a point in the opposite direction.
-        if (auto node_it = vertex_mapping.find(node.p); node_it != vertex_mapping.end())
-            node.p = node_it->second;
-        else
-            node.p.rotate(-fix_angle);
-    }
-}
-
 void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
 {
 #ifdef ARACHNE_DEBUG
@@ -604,8 +415,8 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
     }
 #endif
 
-    Geometry::VoronoiDiagram voronoi_diagram;
-    construct_voronoi(segments.begin(), segments.end(), &voronoi_diagram);
+    VD voronoi_diagram;
+    voronoi_diagram.construct_voronoi(segments.cbegin(), segments.cend());
 
 #ifdef ARACHNE_DEBUG_VORONOI
     {
@@ -614,126 +425,59 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
     }
 #endif
 
-    // Try to detect cases when some Voronoi vertex is missing and when
-    // the Voronoi diagram is not planar.
-    // When any Voronoi vertex is missing, or the Voronoi diagram is not
-    // planar, rotate the input polygon and try again.
-    const bool   has_missing_voronoi_vertex = detect_missing_voronoi_vertex(voronoi_diagram, segments);
-    // Detection of non-planar Voronoi diagram detects at least GH issues #8474, #8514 and #8446.
-    const bool   is_voronoi_diagram_planar  = Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram);
-    const double fix_angle                  = PI / 6;
-
-    std::unordered_map<Point, Point, PointHash> vertex_mapping;
-    // polys_copy is referenced through items stored in the std::vector segments.
-    Polygons                                    polys_copy = polys;
-    if (has_missing_voronoi_vertex || !is_voronoi_diagram_planar) {
-        if (has_missing_voronoi_vertex)
-            BOOST_LOG_TRIVIAL(warning) << "Detected missing Voronoi vertex, input polygons will be rotated back and forth.";
-        else if (!is_voronoi_diagram_planar)
-            BOOST_LOG_TRIVIAL(warning) << "Detected non-planar Voronoi diagram, input polygons will be rotated back and forth.";
-
-        vertex_mapping = try_to_fix_degenerated_voronoi_diagram_by_rotation(voronoi_diagram, polys, polys_copy, segments, fix_angle);
-
-        assert(!detect_missing_voronoi_vertex(voronoi_diagram, segments));
-        assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram));
-        if (detect_missing_voronoi_vertex(voronoi_diagram, segments))
-            BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex even after the rotation of input.";
-        else if (!Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram))
-            BOOST_LOG_TRIVIAL(error) << "Detected non-planar Voronoi diagram even after the rotation of input.";
-    }
-
-    bool degenerated_voronoi_diagram = has_missing_voronoi_vertex || !is_voronoi_diagram_planar;
-
-process_voronoi_diagram:
     assert(this->graph.edges.empty() && this->graph.nodes.empty() && this->vd_edge_to_he_edge.empty() && this->vd_node_to_he_node.empty());
-    for (vd_t::cell_type cell : voronoi_diagram.cells()) {
+    for (const VD::cell_type &cell : voronoi_diagram.cells()) {
         if (!cell.incident_edge())
             continue; // There is no spoon
 
-        Point start_source_point;
-        Point end_source_point;
-        vd_t::edge_type* starting_vonoroi_edge = nullptr;
-        vd_t::edge_type* ending_vonoroi_edge = nullptr;
+        Point                start_source_point;
+        Point                end_source_point;
+        const VD::edge_type *starting_voronoi_edge = nullptr;
+        const VD::edge_type *ending_voronoi_edge   = nullptr;
         // Compute and store result in above variables
-        
+
         if (cell.contains_point()) {
-            const bool keep_going = computePointCellRange(cell, start_source_point, end_source_point, starting_vonoroi_edge, ending_vonoroi_edge, segments);
+            const bool keep_going = computePointCellRange(cell, start_source_point, end_source_point, starting_voronoi_edge, ending_voronoi_edge, segments);
             if (!keep_going)
                 continue;
         } else {
             assert(cell.contains_segment());
-            computeSegmentCellRange(cell, start_source_point, end_source_point, starting_vonoroi_edge, ending_vonoroi_edge, segments);
+            Geometry::SegmentCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_segment_cell_range(cell, segments.cbegin(), segments.cend());
+            assert(cell_range.is_valid());
+            start_source_point    = cell_range.segment_start_point;
+            end_source_point      = cell_range.segment_end_point;
+            starting_voronoi_edge = cell_range.edge_begin;
+            ending_voronoi_edge   = cell_range.edge_end;
         }
-        
-        if (!starting_vonoroi_edge || !ending_vonoroi_edge) {
+
+        if (!starting_voronoi_edge || !ending_voronoi_edge) {
             assert(false && "Each cell should start / end in a polygon vertex");
             continue;
         }
-        
+
         // Copy start to end edge to graph
-        edge_t* prev_edge = nullptr;
-        assert(VoronoiUtils::p(starting_vonoroi_edge->vertex1()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(starting_vonoroi_edge->vertex1()).x() >= std::numeric_limits<coord_t>::lowest());
-        assert(VoronoiUtils::p(starting_vonoroi_edge->vertex1()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(starting_vonoroi_edge->vertex1()).y() >= std::numeric_limits<coord_t>::lowest());
-        transferEdge(start_source_point, VoronoiUtils::p(starting_vonoroi_edge->vertex1()).cast<coord_t>(), *starting_vonoroi_edge, prev_edge, start_source_point, end_source_point, segments);
-        node_t* starting_node = vd_node_to_he_node[starting_vonoroi_edge->vertex0()];
+        assert(Geometry::VoronoiUtils::is_in_range<coord_t>(*starting_voronoi_edge));
+        edge_t *prev_edge = nullptr;
+        transferEdge(start_source_point, Geometry::VoronoiUtils::to_point(starting_voronoi_edge->vertex1()).cast<coord_t>(), *starting_voronoi_edge, prev_edge, start_source_point, end_source_point, segments);
+        node_t *starting_node                    = vd_node_to_he_node[starting_voronoi_edge->vertex0()];
         starting_node->data.distance_to_boundary = 0;
 
         constexpr bool is_next_to_start_or_end = true;
         graph.makeRib(prev_edge, start_source_point, end_source_point, is_next_to_start_or_end);
-        for (vd_t::edge_type* vd_edge = starting_vonoroi_edge->next(); vd_edge != ending_vonoroi_edge; vd_edge = vd_edge->next()) {
+        for (const VD::edge_type* vd_edge = starting_voronoi_edge->next(); vd_edge != ending_voronoi_edge; vd_edge = vd_edge->next()) {
             assert(vd_edge->is_finite());
+            assert(Geometry::VoronoiUtils::is_in_range<coord_t>(*vd_edge));
 
-            assert(VoronoiUtils::p(vd_edge->vertex0()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge->vertex0()).x() >= std::numeric_limits<coord_t>::lowest());
-            assert(VoronoiUtils::p(vd_edge->vertex0()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge->vertex0()).y() >= std::numeric_limits<coord_t>::lowest());
-            assert(VoronoiUtils::p(vd_edge->vertex1()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge->vertex1()).x() >= std::numeric_limits<coord_t>::lowest());
-            assert(VoronoiUtils::p(vd_edge->vertex1()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(vd_edge->vertex1()).y() >= std::numeric_limits<coord_t>::lowest());
-
-            Point v1 = VoronoiUtils::p(vd_edge->vertex0()).cast<coord_t>();
-            Point v2 = VoronoiUtils::p(vd_edge->vertex1()).cast<coord_t>();
+            Point v1 = Geometry::VoronoiUtils::to_point(vd_edge->vertex0()).cast<coord_t>();
+            Point v2 = Geometry::VoronoiUtils::to_point(vd_edge->vertex1()).cast<coord_t>();
             transferEdge(v1, v2, *vd_edge, prev_edge, start_source_point, end_source_point, segments);
-
-            graph.makeRib(prev_edge, start_source_point, end_source_point, vd_edge->next() == ending_vonoroi_edge);
+            graph.makeRib(prev_edge, start_source_point, end_source_point, vd_edge->next() == ending_voronoi_edge);
         }
 
-        assert(VoronoiUtils::p(starting_vonoroi_edge->vertex0()).x() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(starting_vonoroi_edge->vertex0()).x() >= std::numeric_limits<coord_t>::lowest());
-        assert(VoronoiUtils::p(starting_vonoroi_edge->vertex0()).y() <= std::numeric_limits<coord_t>::max() && VoronoiUtils::p(starting_vonoroi_edge->vertex0()).y() >= std::numeric_limits<coord_t>::lowest());
-        transferEdge(VoronoiUtils::p(ending_vonoroi_edge->vertex0()).cast<coord_t>(), end_source_point, *ending_vonoroi_edge, prev_edge, start_source_point, end_source_point, segments);
+        transferEdge(Geometry::VoronoiUtils::to_point(ending_voronoi_edge->vertex0()).cast<coord_t>(), end_source_point, *ending_voronoi_edge, prev_edge, start_source_point, end_source_point, segments);
         prev_edge->to->data.distance_to_boundary = 0;
     }
 
-    // For some input polygons, as in GH issues #8474 and #8514 resulting Voronoi diagram is degenerated because it is not planar.
-    // When this degenerated Voronoi diagram is processed, the resulting half-edge structure contains some edges that don't have
-    // a twin edge. Based on this, we created a fast mechanism that detects those causes and tries to recompute the Voronoi
-    // diagram on slightly rotated input polygons that usually make the Voronoi generator generate a non-degenerated Voronoi diagram.
-    if (!degenerated_voronoi_diagram && has_missing_twin_edge(this->graph)) {
-        BOOST_LOG_TRIVIAL(warning) << "Detected degenerated Voronoi diagram, input polygons will be rotated back and forth.";
-        degenerated_voronoi_diagram = true;
-        vertex_mapping = try_to_fix_degenerated_voronoi_diagram_by_rotation(voronoi_diagram, polys, polys_copy, segments, fix_angle);
-
-        assert(!detect_missing_voronoi_vertex(voronoi_diagram, segments));
-        if (detect_missing_voronoi_vertex(voronoi_diagram, segments))
-            BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex after the rotation of input.";
-
-        assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(voronoi_diagram));
-
-        this->graph.edges.clear();
-        this->graph.nodes.clear();
-        this->vd_edge_to_he_edge.clear();
-        this->vd_node_to_he_node.clear();
-
-        goto process_voronoi_diagram;
-    }
-
-    if (degenerated_voronoi_diagram) {
-        assert(!has_missing_twin_edge(this->graph));
-
-        if (has_missing_twin_edge(this->graph))
-            BOOST_LOG_TRIVIAL(error) << "Detected degenerated Voronoi diagram even after the rotation of input.";
-    }
-
-    if (degenerated_voronoi_diagram)
-        rotate_back_skeletal_trapezoidation_graph_after_fix(this->graph, fix_angle, vertex_mapping);
-
 #ifdef ARACHNE_DEBUG
     assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(voronoi_diagram));
 #endif
@@ -742,7 +486,7 @@ process_voronoi_diagram:
 
     graph.collapseSmallEdges();
 
-    // Set [incident_edge] the the first possible edge that way we can iterate over all reachable edges from node.incident_edge,
+    // Set [incident_edge] the first possible edge that way we can iterate over all reachable edges from node.incident_edge,
     // without needing to iterate backward
     for (edge_t& edge : graph.edges)
         if (!edge.prev)
@@ -751,7 +495,7 @@ process_voronoi_diagram:
 
 void SkeletalTrapezoidation::separatePointyQuadEndNodes()
 {
-    std::unordered_set<node_t*> visited_nodes;
+    NodeSet visited_nodes;
     for (edge_t& edge : graph.edges)
     {
         if (edge.prev) 
@@ -2221,16 +1965,16 @@ void SkeletalTrapezoidation::addToolpathSegment(const ExtrusionJunction& from, c
 
 void SkeletalTrapezoidation::connectJunctions(ptr_vector_t<LineJunctions>& edge_junctions)
 {
-    std::unordered_set<edge_t*> unprocessed_quad_starts(graph.edges.size() * 5 / 2);
+    EdgeSet unprocessed_quad_starts(graph.edges.size() * 5 / 2);
     for (edge_t& edge : graph.edges)
     {
         if (!edge.prev)
         {
-            unprocessed_quad_starts.insert(&edge);
+            unprocessed_quad_starts.emplace(&edge);
         }
     }
 
-    std::unordered_set<edge_t*> passed_odd_edges;
+    EdgeSet passed_odd_edges;
 
     while (!unprocessed_quad_starts.empty())
     {
diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
index 819b71367..75ff94015 100644
--- a/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
+++ b/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
@@ -7,9 +7,10 @@
 #include <boost/polygon/voronoi.hpp>
 
 #include <memory> // smart pointers
-#include <unordered_map>
 #include <utility> // pair
 
+#include <ankerl/unordered_dense.h>
+
 #include "utils/HalfEdgeGraph.hpp"
 #include "utils/PolygonsSegmentIndex.hpp"
 #include "utils/ExtrusionJunction.hpp"
@@ -23,8 +24,9 @@
 //#define ARACHNE_DEBUG
 //#define ARACHNE_DEBUG_VORONOI
 
-namespace Slic3r::Arachne
-{
+namespace Slic3r::Arachne {
+
+using VD = Slic3r::Geometry::VoronoiDiagram;
 
 /*!
  * Main class of the dynamic beading strategies.
@@ -47,8 +49,6 @@ deposition modeling" by Kuipers et al.
  */
 class SkeletalTrapezoidation
 {
-    using pos_t = double;
-    using vd_t = boost::polygon::voronoi_diagram<pos_t>;
     using graph_t = SkeletalTrapezoidationGraph;
     using edge_t = STHalfEdge;
     using node_t = STHalfEdgeNode;
@@ -79,7 +79,11 @@ class SkeletalTrapezoidation
     const BeadingStrategy& beading_strategy;
 
 public:
-    using Segment = PolygonsSegmentIndex;
+    using Segment  = PolygonsSegmentIndex;
+    using NodeSet = ankerl::unordered_dense::set<node_t *>;
+    using EdgeSet = ankerl::unordered_dense::set<edge_t *>;
+    using EdgeMap = ankerl::unordered_dense::map<const VD::edge_type *, edge_t *>;
+    using NodeMap = ankerl::unordered_dense::map<const VD::vertex_type *, node_t *>;
 
     /*!
      * Construct a new trapezoidation problem to solve.
@@ -163,9 +167,9 @@ protected:
      * mapping each voronoi VD edge to the corresponding halfedge HE edge
      * In case the result segment is discretized, we map the VD edge to the *last* HE edge
      */
-    std::unordered_map<vd_t::edge_type*, edge_t*> vd_edge_to_he_edge;
-    std::unordered_map<vd_t::vertex_type*, node_t*> vd_node_to_he_node;
-    node_t& makeNode(vd_t::vertex_type& vd_node, Point p); //!< Get the node which the VD node maps to, or create a new mapping if there wasn't any yet.
+    EdgeMap vd_edge_to_he_edge;
+    NodeMap vd_node_to_he_node;
+    node_t &makeNode(const VD::vertex_type &vd_node, Point p); //!< Get the node which the VD node maps to, or create a new mapping if there wasn't any yet.
 
     /*!
      * (Eventual) returned 'polylines per index' result (from generateToolpaths):
@@ -176,7 +180,7 @@ protected:
      * Transfer an edge from the VD to the HE and perform discretization of parabolic edges (and vertex-vertex edges)
      * \p prev_edge serves as input and output. May be null as input.
      */
-    void transferEdge(Point from, Point to, vd_t::edge_type& vd_edge, edge_t*& prev_edge, Point& start_source_point, Point& end_source_point, const std::vector<Segment>& segments);
+    void transferEdge(Point from, Point to, const VD::edge_type &vd_edge, edge_t *&prev_edge, Point &start_source_point, Point &end_source_point, const std::vector<Segment> &segments);
 
     /*!
      * Discretize a Voronoi edge that represents the medial axis of a vertex-
@@ -203,7 +207,7 @@ protected:
      * \return A number of coordinates along the edge where the edge is broken
      * up into discrete pieces.
      */
-    std::vector<Point> discretize(const vd_t::edge_type& segment, const std::vector<Segment>& segments);
+    Points discretize(const VD::edge_type& segment, const std::vector<Segment>& segments);
 
     /*!
      * Compute the range of line segments that surround a cell of the skeletal
@@ -229,33 +233,7 @@ protected:
      * /return Whether the cell is inside of the polygon. If it's outside of the
      * polygon we should skip processing it altogether.
      */
-    bool computePointCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::edge_type*& ending_vd_edge, const std::vector<Segment>& segments);
-
-    /*!
-     * Compute the range of line segments that surround a cell of the skeletal
-     * graph that belongs to a line segment of the medial axis.
-     *
-     * This should only be used on cells that belong to a central line segment
-     * of the skeletal graph, e.g. trapezoid cells, not triangular cells.
-     *
-     * The resulting line segments is just the first and the last segment. They
-     * are linked to the neighboring segments, so you can iterate over the
-     * segments until you reach the last segment.
-     * \param cell The cell to compute the range of line segments for.
-     * \param[out] start_source_point The start point of the source segment of
-     * this cell.
-     * \param[out] end_source_point The end point of the source segment of this
-     * cell.
-     * \param[out] starting_vd_edge The edge of the Voronoi diagram where the
-     * loop around the cell starts.
-     * \param[out] ending_vd_edge The edge of the Voronoi diagram where the loop
-     * around the cell ends.
-     * \param points All vertices of the input Polygons.
-     * \param segments All edges of the input Polygons.
-     * /return Whether the cell is inside of the polygon. If it's outside of the
-     * polygon we should skip processing it altogether.
-     */
-    void computeSegmentCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::edge_type*& ending_vd_edge, const std::vector<Segment>& segments);
+    static bool computePointCellRange(const VD::cell_type &cell, Point &start_source_point, Point &end_source_point, const VD::edge_type *&starting_vd_edge, const VD::edge_type *&ending_vd_edge, const std::vector<Segment> &segments);
 
     /*!
      * For VD cells associated with an input polygon vertex, we need to separate the node at the end and start of the cell into two
diff --git a/src/libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp b/src/libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp
index 6eff3d62e..3258b41c7 100644
--- a/src/libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp
+++ b/src/libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp
@@ -27,5 +27,24 @@ public:
 
 } // namespace Slic3r::Arachne
 
+namespace boost::polygon {
+
+template<> struct geometry_concept<Slic3r::Arachne::PolygonsSegmentIndex>
+{
+    typedef segment_concept type;
+};
+
+template<> struct segment_traits<Slic3r::Arachne::PolygonsSegmentIndex>
+{
+    typedef coord_t       coordinate_type;
+    typedef Slic3r::Point point_type;
+
+    static inline point_type get(const Slic3r::Arachne::PolygonsSegmentIndex &CSegment, direction_1d dir)
+    {
+        return dir.to_int() ? CSegment.to() : CSegment.from();
+    }
+};
+
+} // namespace boost::polygon
 
 #endif//UTILS_POLYGONS_SEGMENT_INDEX_H
diff --git a/src/libslic3r/Arachne/utils/VoronoiUtils.cpp b/src/libslic3r/Arachne/utils/VoronoiUtils.cpp
deleted file mode 100644
index 82bd79523..000000000
--- a/src/libslic3r/Arachne/utils/VoronoiUtils.cpp
+++ /dev/null
@@ -1,251 +0,0 @@
-//Copyright (c) 2021 Ultimaker B.V.
-//CuraEngine is released under the terms of the AGPLv3 or higher.
-
-#include <stack>
-#include <optional>
-#include <boost/log/trivial.hpp>
-
-#include "linearAlg2D.hpp"
-#include "VoronoiUtils.hpp"
-
-namespace Slic3r::Arachne
-{
-
-Vec2i64 VoronoiUtils::p(const vd_t::vertex_type *node)
-{
-    const double x = node->x();
-    const double y = node->y();
-    assert(std::isfinite(x) && std::isfinite(y));
-    assert(x <= double(std::numeric_limits<int64_t>::max()) && x >= std::numeric_limits<int64_t>::lowest());
-    assert(y <= double(std::numeric_limits<int64_t>::max()) && y >= std::numeric_limits<int64_t>::lowest());
-    return {int64_t(x + 0.5 - (x < 0)), int64_t(y + 0.5 - (y < 0))}; // Round to the nearest integer coordinates.
-}
-
-Point VoronoiUtils::getSourcePoint(const vd_t::cell_type& cell, const std::vector<Segment>& segments)
-{
-    assert(cell.contains_point());
-    if(!cell.contains_point())
-        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source point!";
-
-    switch (cell.source_category()) {
-        case boost::polygon::SOURCE_CATEGORY_SINGLE_POINT:
-            assert(false && "Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!\n");
-            BOOST_LOG_TRIVIAL(error) << "Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!";
-            break;
-        case boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT:
-            assert(cell.source_index() < segments.size());
-            return segments[cell.source_index()].to();
-            break;
-        case boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT:
-            assert(cell.source_index() < segments.size());
-            return segments[cell.source_index()].from();
-            break;
-        default:
-            assert(false && "getSourcePoint should only be called on point cells!\n");
-            break;
-    }
-
-    assert(false && "cell.source_category() is equal to an invalid value!\n");
-    BOOST_LOG_TRIVIAL(error) << "cell.source_category() is equal to an invalid value!";
-    return {};
-}
-
-PolygonsPointIndex VoronoiUtils::getSourcePointIndex(const vd_t::cell_type& cell, const std::vector<Segment>& segments)
-{
-    assert(cell.contains_point());
-    if(!cell.contains_point())
-        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source point!";
-
-    assert(cell.source_category() != boost::polygon::SOURCE_CATEGORY_SINGLE_POINT);
-    switch (cell.source_category()) {
-        case boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT: {
-            assert(cell.source_index() < segments.size());
-            PolygonsPointIndex ret = segments[cell.source_index()];
-            ++ret;
-            return ret;
-            break;
-        }
-        case boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT: {
-            assert(cell.source_index() < segments.size());
-            return segments[cell.source_index()];
-            break;
-        }
-        default:
-            assert(false && "getSourcePoint should only be called on point cells!\n");
-            break;
-    }
-    PolygonsPointIndex ret = segments[cell.source_index()];
-    return ++ret;
-}
-
-const VoronoiUtils::Segment &VoronoiUtils::getSourceSegment(const vd_t::cell_type &cell, const std::vector<Segment> &segments)
-{
-    assert(cell.contains_segment());
-    if (!cell.contains_segment())
-        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source segment!";
-
-    return segments[cell.source_index()];
-}
-
-class PointMatrix
-{
-public:
-    double matrix[4];
-
-    PointMatrix()
-    {
-        matrix[0] = 1;
-        matrix[1] = 0;
-        matrix[2] = 0;
-        matrix[3] = 1;
-    }
-
-    PointMatrix(double rotation)
-    {
-        rotation = rotation / 180 * M_PI;
-        matrix[0] = cos(rotation);
-        matrix[1] = -sin(rotation);
-        matrix[2] = -matrix[1];
-        matrix[3] = matrix[0];
-    }
-
-    PointMatrix(const Point p)
-    {
-        matrix[0] = p.x();
-        matrix[1] = p.y();
-        double f = sqrt((matrix[0] * matrix[0]) + (matrix[1] * matrix[1]));
-        matrix[0] /= f;
-        matrix[1] /= f;
-        matrix[2] = -matrix[1];
-        matrix[3] = matrix[0];
-    }
-
-    static PointMatrix scale(double s)
-    {
-        PointMatrix ret;
-        ret.matrix[0] = s;
-        ret.matrix[3] = s;
-        return ret;
-    }
-
-    Point apply(const Point p) const
-    {
-        return Point(coord_t(p.x() * matrix[0] + p.y() * matrix[1]), coord_t(p.x() * matrix[2] + p.y() * matrix[3]));
-    }
-
-    Point unapply(const Point p) const
-    {
-        return Point(coord_t(p.x() * matrix[0] + p.y() * matrix[2]), coord_t(p.x() * matrix[1] + p.y() * matrix[3]));
-    }
-};
-std::vector<Point> VoronoiUtils::discretizeParabola(const Point& p, const Segment& segment, Point s, Point e, coord_t approximate_step_size, float transitioning_angle)
-{
-    std::vector<Point> discretized;
-    // x is distance of point projected on the segment ab
-    // xx is point projected on the segment ab
-    const Point a = segment.from();
-    const Point b = segment.to();
-    const Point ab = b - a;
-    const Point as = s - a;
-    const Point ae = e - a;
-    const coord_t ab_size = ab.cast<int64_t>().norm();
-    const coord_t sx = as.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
-    const coord_t ex = ae.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
-    const coord_t sxex = ex - sx;
-
-    assert((as.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
-    assert((ae.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
-
-    const Point ap = p - a;
-    const coord_t px = ap.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
-
-    assert((ap.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
-
-    Point pxx;
-    Line(a, b).distance_to_infinite_squared(p, &pxx);
-    const Point ppxx = pxx - p;
-    const coord_t d = ppxx.cast<int64_t>().norm();
-    const PointMatrix rot = PointMatrix(ppxx.rotate_90_degree_ccw());
-
-    if (d == 0)
-    {
-        discretized.emplace_back(s);
-        discretized.emplace_back(e);
-        return discretized;
-    }
-    
-    const float marking_bound = atan(transitioning_angle * 0.5);
-    int64_t msx = - marking_bound * int64_t(d); // projected marking_start
-    int64_t mex = marking_bound * int64_t(d); // projected marking_end
-
-    assert(msx <= std::numeric_limits<coord_t>::max());
-    assert(double(msx) * double(msx) <= double(std::numeric_limits<int64_t>::max()));
-    assert(mex <= std::numeric_limits<coord_t>::max());
-    assert(double(msx) * double(msx) / double(2 * d) + double(d / 2) <= std::numeric_limits<coord_t>::max());
-
-    const coord_t marking_start_end_h = msx * msx / (2 * d) + d / 2;
-    Point marking_start = rot.unapply(Point(coord_t(msx), marking_start_end_h)) + pxx;
-    Point marking_end = rot.unapply(Point(coord_t(mex), marking_start_end_h)) + pxx;
-    const int dir = (sx > ex) ? -1 : 1;
-    if (dir < 0)
-    {
-        std::swap(marking_start, marking_end);
-        std::swap(msx, mex);
-    }
-    
-    bool add_marking_start = msx * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && msx * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
-    bool add_marking_end = mex * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && mex * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
-
-    const Point apex = rot.unapply(Point(0, d / 2)) + pxx;
-    bool add_apex = int64_t(sx - px) * int64_t(dir) < 0 && int64_t(ex - px) * int64_t(dir) > 0;
-
-    assert(!(add_marking_start && add_marking_end) || add_apex);
-    if(add_marking_start && add_marking_end && !add_apex)
-    {
-        BOOST_LOG_TRIVIAL(warning) << "Failing to discretize parabola! Must add an apex or one of the endpoints.";
-    }
-    
-    const coord_t step_count = static_cast<coord_t>(static_cast<float>(std::abs(ex - sx)) / approximate_step_size + 0.5);
-    
-    discretized.emplace_back(s);
-    for (coord_t step = 1; step < step_count; step++)
-    {
-        assert(double(sxex) * double(step) <= double(std::numeric_limits<int64_t>::max()));
-        const int64_t x = int64_t(sx) + int64_t(sxex) * int64_t(step) / int64_t(step_count) - int64_t(px);
-        assert(double(x) * double(x) <= double(std::numeric_limits<int64_t>::max()));
-        assert(double(x) * double(x) / double(2 * d) + double(d / 2) <= double(std::numeric_limits<int64_t>::max()));
-        const int64_t y = int64_t(x) * int64_t(x) / int64_t(2 * d) + int64_t(d / 2);
-        
-        if (add_marking_start && msx * int64_t(dir) < int64_t(x) * int64_t(dir))
-        {
-            discretized.emplace_back(marking_start);
-            add_marking_start = false;
-        }
-        if (add_apex && int64_t(x) * int64_t(dir) > 0)
-        {
-            discretized.emplace_back(apex);
-            add_apex = false; // only add the apex just before the 
-        }
-        if (add_marking_end && mex * int64_t(dir) < int64_t(x) * int64_t(dir))
-        {
-            discretized.emplace_back(marking_end);
-            add_marking_end = false;
-        }
-        assert(x <= std::numeric_limits<coord_t>::max() && x >= std::numeric_limits<coord_t>::lowest());
-        assert(y <= std::numeric_limits<coord_t>::max() && y >= std::numeric_limits<coord_t>::lowest());
-        const Point result = rot.unapply(Point(x, y)) + pxx;
-        discretized.emplace_back(result);
-    }
-    if (add_apex)
-    {
-        discretized.emplace_back(apex);
-    }
-    if (add_marking_end)
-    {
-        discretized.emplace_back(marking_end);
-    }
-    discretized.emplace_back(e);
-    return discretized;
-}
-
-}//namespace Slic3r::Arachne
diff --git a/src/libslic3r/Arachne/utils/VoronoiUtils.hpp b/src/libslic3r/Arachne/utils/VoronoiUtils.hpp
deleted file mode 100644
index aa4693643..000000000
--- a/src/libslic3r/Arachne/utils/VoronoiUtils.hpp
+++ /dev/null
@@ -1,47 +0,0 @@
-//Copyright (c) 2020 Ultimaker B.V.
-//CuraEngine is released under the terms of the AGPLv3 or higher.
-
-
-#ifndef UTILS_VORONOI_UTILS_H
-#define UTILS_VORONOI_UTILS_H
-
-#include <vector>
-
-
-#include <boost/polygon/voronoi.hpp>
-
-#include "PolygonsSegmentIndex.hpp"
-
-namespace Slic3r::Arachne
-{
-
-/*!
- */
-class VoronoiUtils
-{
-public:
-    using Segment        = PolygonsSegmentIndex;
-    using voronoi_data_t = double;
-    using vd_t           = boost::polygon::voronoi_diagram<voronoi_data_t>;
-
-    static Point              getSourcePoint(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
-    static const Segment     &getSourceSegment(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
-    static PolygonsPointIndex getSourcePointIndex(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
-
-    static Vec2i64 p(const vd_t::vertex_type *node);
-
-    /*!
-     * Discretize a parabola based on (approximate) step size.
-     * The \p approximate_step_size is measured parallel to the \p source_segment, not along the parabola.
-     */
-    static std::vector<Point> discretizeParabola(const Point &source_point, const Segment &source_segment, Point start, Point end, coord_t approximate_step_size, float transitioning_angle);
-
-    static inline bool is_finite(const VoronoiUtils::vd_t::vertex_type &vertex)
-    {
-        return std::isfinite(vertex.x()) && std::isfinite(vertex.y());
-    }
-};
-
-} // namespace Slic3r::Arachne
-
-#endif // UTILS_VORONOI_UTILS_H
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 108b24f30..f93f233d2 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -171,9 +171,12 @@ set(lisbslic3r_sources
     Geometry/Curves.hpp
     Geometry/MedialAxis.cpp
     Geometry/MedialAxis.hpp
+	Geometry/Voronoi.cpp
     Geometry/Voronoi.hpp
     Geometry/VoronoiOffset.cpp
     Geometry/VoronoiOffset.hpp
+	Geometry/VoronoiUtils.hpp
+    Geometry/VoronoiUtils.cpp
     Geometry/VoronoiUtilsCgal.cpp
     Geometry/VoronoiUtilsCgal.hpp
     Geometry/VoronoiVisualUtils.hpp
@@ -391,8 +394,6 @@ set(lisbslic3r_sources
     Arachne/utils/PolygonsSegmentIndex.hpp
     Arachne/utils/PolylineStitcher.hpp
     Arachne/utils/PolylineStitcher.cpp
-    Arachne/utils/VoronoiUtils.hpp
-    Arachne/utils/VoronoiUtils.cpp
     Arachne/SkeletalTrapezoidation.hpp
     Arachne/SkeletalTrapezoidation.cpp
     Arachne/SkeletalTrapezoidationEdge.hpp
diff --git a/src/libslic3r/Geometry/MedialAxis.cpp b/src/libslic3r/Geometry/MedialAxis.cpp
index 39677c928..44eccd985 100644
--- a/src/libslic3r/Geometry/MedialAxis.cpp
+++ b/src/libslic3r/Geometry/MedialAxis.cpp
@@ -448,7 +448,7 @@ MedialAxis::MedialAxis(double min_width, double max_width, const ExPolygon &expo
 
 void MedialAxis::build(ThickPolylines* polylines)
 {
-    construct_voronoi(m_lines.begin(), m_lines.end(), &m_vd);
+    m_vd.construct_voronoi(m_lines.begin(), m_lines.end());
     Slic3r::Voronoi::annotate_inside_outside(m_vd, m_lines);
 //    static constexpr double threshold_alpha = M_PI / 12.; // 30 degrees
 //    std::vector<Vec2d> skeleton_edges = Slic3r::Voronoi::skeleton_edges_rough(vd, lines, threshold_alpha);
diff --git a/src/libslic3r/Geometry/Voronoi.cpp b/src/libslic3r/Geometry/Voronoi.cpp
new file mode 100644
index 000000000..e1df7322a
--- /dev/null
+++ b/src/libslic3r/Geometry/Voronoi.cpp
@@ -0,0 +1,354 @@
+#include "Voronoi.hpp"
+
+#include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
+#include "libslic3r/Geometry/VoronoiUtils.hpp"
+#include "libslic3r/Geometry/VoronoiUtilsCgal.hpp"
+#include "libslic3r/MultiMaterialSegmentation.hpp"
+
+#include <boost/log/trivial.hpp>
+
+namespace Slic3r::Geometry {
+
+using PolygonsSegmentIndexConstIt = std::vector<Arachne::PolygonsSegmentIndex>::const_iterator;
+using LinesIt                     = Lines::iterator;
+using ColoredLinesConstIt         = ColoredLines::const_iterator;
+
+// Explicit template instantiation.
+template void VoronoiDiagram::construct_voronoi(LinesIt, LinesIt, bool);
+template void VoronoiDiagram::construct_voronoi(ColoredLinesConstIt, ColoredLinesConstIt, bool);
+template void VoronoiDiagram::construct_voronoi(PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt, bool);
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    void>::type
+VoronoiDiagram::construct_voronoi(const SegmentIterator segment_begin, const SegmentIterator segment_end, const bool try_to_repair_if_needed) {
+    boost::polygon::construct_voronoi(segment_begin, segment_end, &m_voronoi_diagram);
+    if (try_to_repair_if_needed) {
+        if (m_issue_type = detect_known_issues(*this, segment_begin, segment_end); m_issue_type != IssueType::NO_ISSUE_DETECTED) {
+            if (m_issue_type == IssueType::MISSING_VORONOI_VERTEX) {
+                BOOST_LOG_TRIVIAL(warning) << "Detected missing Voronoi vertex, input polygons will be rotated back and forth.";
+            } else if (m_issue_type == IssueType::NON_PLANAR_VORONOI_DIAGRAM) {
+                BOOST_LOG_TRIVIAL(warning) << "Detected non-planar Voronoi diagram, input polygons will be rotated back and forth.";
+            } else if (m_issue_type == IssueType::VORONOI_EDGE_INTERSECTING_INPUT_SEGMENT) {
+                BOOST_LOG_TRIVIAL(warning) << "Detected Voronoi edge intersecting input segment, input polygons will be rotated back and forth.";
+            } else if (m_issue_type == IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX) {
+                BOOST_LOG_TRIVIAL(warning) << "Detected finite Voronoi vertex with non finite vertex, input polygons will be rotated back and forth.";
+            } else {
+                BOOST_LOG_TRIVIAL(error) << "Detected unknown Voronoi diagram issue, input polygons will be rotated back and forth.";
+            }
+
+            if (m_issue_type = try_to_repair_degenerated_voronoi_diagram(segment_begin, segment_end); m_issue_type != IssueType::NO_ISSUE_DETECTED) {
+                if (m_issue_type == IssueType::MISSING_VORONOI_VERTEX) {
+                    BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex even after the rotation of input.";
+                } else if (m_issue_type == IssueType::NON_PLANAR_VORONOI_DIAGRAM) {
+                    BOOST_LOG_TRIVIAL(error) << "Detected non-planar Voronoi diagram even after the rotation of input.";
+                } else if (m_issue_type == IssueType::VORONOI_EDGE_INTERSECTING_INPUT_SEGMENT) {
+                    BOOST_LOG_TRIVIAL(error) << "Detected Voronoi edge intersecting input segment even after the rotation of input.";
+                } else if (m_issue_type == IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX) {
+                    BOOST_LOG_TRIVIAL(error) << "Detected finite Voronoi vertex with non finite vertex even after the rotation of input.";
+                } else {
+                    BOOST_LOG_TRIVIAL(error) << "Detected unknown Voronoi diagram issue even after the rotation of input.";
+                }
+
+                m_state = State::REPAIR_UNSUCCESSFUL;
+            } else {
+                m_state = State::REPAIR_SUCCESSFUL;
+            }
+        } else {
+            m_state      = State::REPAIR_NOT_NEEDED;
+            m_issue_type = IssueType::NO_ISSUE_DETECTED;
+        }
+    } else {
+        m_state      = State::UNKNOWN;
+        m_issue_type = IssueType::UNKNOWN;
+    }
+}
+
+void VoronoiDiagram::clear()
+{
+    if (m_is_modified) {
+        m_vertices.clear();
+        m_edges.clear();
+        m_cells.clear();
+        m_is_modified = false;
+    } else {
+        m_voronoi_diagram.clear();
+    }
+
+    m_state      = State::UNKNOWN;
+    m_issue_type = IssueType::UNKNOWN;
+}
+
+void VoronoiDiagram::copy_to_local(voronoi_diagram_type &voronoi_diagram) {
+    m_edges.clear();
+    m_cells.clear();
+    m_vertices.clear();
+
+    // Copy Voronoi edges.
+    m_edges.reserve(voronoi_diagram.num_edges());
+    for (const edge_type &edge : voronoi_diagram.edges()) {
+        m_edges.emplace_back(edge.is_linear(), edge.is_primary());
+        m_edges.back().color(edge.color());
+    }
+
+    // Copy Voronoi cells.
+    m_cells.reserve(voronoi_diagram.num_cells());
+    for (const cell_type &cell : voronoi_diagram.cells()) {
+        m_cells.emplace_back(cell.source_index(), cell.source_category());
+        m_cells.back().color(cell.color());
+
+        if (cell.incident_edge()) {
+            size_t incident_edge_idx = cell.incident_edge() - voronoi_diagram.edges().data();
+            m_cells.back().incident_edge(&m_edges[incident_edge_idx]);
+        }
+    }
+
+    // Copy Voronoi vertices.
+    m_vertices.reserve(voronoi_diagram.num_vertices());
+    for (const vertex_type &vertex : voronoi_diagram.vertices()) {
+        m_vertices.emplace_back(vertex.x(), vertex.y());
+        m_vertices.back().color(vertex.color());
+
+        if (vertex.incident_edge()) {
+            size_t incident_edge_idx = vertex.incident_edge() - voronoi_diagram.edges().data();
+            m_vertices.back().incident_edge(&m_edges[incident_edge_idx]);
+        }
+    }
+
+    // Assign all pointers for each Voronoi edge.
+    for (const edge_type &old_edge : voronoi_diagram.edges()) {
+        size_t     edge_idx = &old_edge - voronoi_diagram.edges().data();
+        edge_type &new_edge = m_edges[edge_idx];
+
+        if (old_edge.cell()) {
+            size_t cell_idx = old_edge.cell() - voronoi_diagram.cells().data();
+            new_edge.cell(&m_cells[cell_idx]);
+        }
+
+        if (old_edge.vertex0()) {
+            size_t vertex0_idx = old_edge.vertex0() - voronoi_diagram.vertices().data();
+            new_edge.vertex0(&m_vertices[vertex0_idx]);
+        }
+
+        if (old_edge.twin()) {
+            size_t twin_edge_idx = old_edge.twin() - voronoi_diagram.edges().data();
+            new_edge.twin(&m_edges[twin_edge_idx]);
+        }
+
+        if (old_edge.next()) {
+            size_t next_edge_idx = old_edge.next() - voronoi_diagram.edges().data();
+            new_edge.next(&m_edges[next_edge_idx]);
+        }
+
+        if (old_edge.prev()) {
+            size_t prev_edge_idx = old_edge.prev() - voronoi_diagram.edges().data();
+            new_edge.prev(&m_edges[prev_edge_idx]);
+        }
+    }
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    VoronoiDiagram::IssueType>::type
+VoronoiDiagram::detect_known_issues(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end)
+{
+    if (has_finite_edge_with_non_finite_vertex(voronoi_diagram)) {
+        return IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX;
+    } else if (const IssueType cell_issue_type = detect_known_voronoi_cell_issues(voronoi_diagram, segment_begin, segment_end); cell_issue_type != IssueType::NO_ISSUE_DETECTED) {
+        return cell_issue_type;
+    } else if (!VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram, segment_begin, segment_end)) {
+        // Detection of non-planar Voronoi diagram detects at least GH issues #8474, #8514 and #8446.
+        return IssueType::NON_PLANAR_VORONOI_DIAGRAM;
+    }
+
+    return IssueType::NO_ISSUE_DETECTED;
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    VoronoiDiagram::IssueType>::type
+VoronoiDiagram::detect_known_voronoi_cell_issues(const VoronoiDiagram &voronoi_diagram,
+                                                 const SegmentIterator segment_begin,
+                                                 const SegmentIterator segment_end)
+{
+    using Segment          = typename std::iterator_traits<SegmentIterator>::value_type;
+    using Point            = typename boost::polygon::segment_point_type<Segment>::type;
+    using SegmentCellRange = SegmentCellRange<Point>;
+
+    for (VD::cell_type cell : voronoi_diagram.cells()) {
+        if (cell.is_degenerate() || !cell.contains_segment())
+            continue; // Skip degenerated cell that has no spoon. Also, skip a cell that doesn't contain a segment.
+
+        if (const SegmentCellRange cell_range = VoronoiUtils::compute_segment_cell_range(cell, segment_begin, segment_end); cell_range.is_valid()) {
+            // Detection if Voronoi edge is intersecting input segment.
+            // It detects this type of issue at least in GH issues #8446, #8474 and #8514.
+
+            const Segment &source_segment      = Geometry::VoronoiUtils::get_source_segment(cell, segment_begin, segment_end);
+            const Vec2d    source_segment_from = boost::polygon::segment_traits<Segment>::get(source_segment, boost::polygon::LOW).template cast<double>();
+            const Vec2d    source_segment_to   = boost::polygon::segment_traits<Segment>::get(source_segment, boost::polygon::HIGH).template cast<double>();
+            const Vec2d    source_segment_vec  = source_segment_to - source_segment_from;
+
+            // All Voronoi vertices must be on the left side of the source segment, otherwise the Voronoi diagram is invalid.
+            for (const VD::edge_type *edge = cell_range.edge_begin; edge != cell_range.edge_end; edge = edge->next()) {
+                if (edge->is_infinite()) {
+                    // When there is a missing Voronoi vertex, we may encounter an infinite Voronoi edge.
+                    // This happens, for example, in GH issue #8846.
+                    return IssueType::MISSING_VORONOI_VERTEX;
+                } else if (const Vec2d edge_v1(edge->vertex1()->x(), edge->vertex1()->y()); Slic3r::cross2(source_segment_vec, edge_v1 - source_segment_from) < 0) {
+                    return IssueType::VORONOI_EDGE_INTERSECTING_INPUT_SEGMENT;
+                }
+            }
+        } else {
+            // When there is a missing Voronoi vertex (especially at one of the endpoints of the input segment),
+            // the returned cell_range is marked as invalid.
+            // It detects this type of issue at least in GH issue #8846.
+            return IssueType::MISSING_VORONOI_VERTEX;
+        }
+    }
+
+    return IssueType::NO_ISSUE_DETECTED;
+}
+
+bool VoronoiDiagram::has_finite_edge_with_non_finite_vertex(const VoronoiDiagram &voronoi_diagram)
+{
+    for (const voronoi_diagram_type::edge_type &edge : voronoi_diagram.edges()) {
+        if (edge.is_finite()) {
+            assert(edge.vertex0() != nullptr && edge.vertex1() != nullptr);
+            if (edge.vertex0() == nullptr || edge.vertex1() == nullptr || !VoronoiUtils::is_finite(*edge.vertex0()) || !VoronoiUtils::is_finite(*edge.vertex1()))
+                return true;
+        }
+    }
+    return false;
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    VoronoiDiagram::IssueType>::type
+VoronoiDiagram::try_to_repair_degenerated_voronoi_diagram(const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    IssueType issue_type = m_issue_type;
+
+    const std::vector<double> fix_angles = {PI / 6, PI / 5, PI / 7, PI / 11};
+    for (const double fix_angle : fix_angles) {
+        issue_type = try_to_repair_degenerated_voronoi_diagram_by_rotation(segment_begin, segment_end, fix_angle);
+        if (issue_type == IssueType::NO_ISSUE_DETECTED) {
+            return issue_type;
+        }
+    }
+
+    return issue_type;
+}
+
+inline VD::vertex_type::color_type encode_input_segment_endpoint(const VD::cell_type::source_index_type cell_source_index, const boost::polygon::direction_1d dir)
+{
+    return (cell_source_index + 1) << 1 | (dir.to_int() ? 1 : 0);
+}
+
+template<typename SegmentIterator>
+inline typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type
+decode_input_segment_endpoint(const VD::vertex_type::color_type color, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    using SegmentType = typename std::iterator_traits<SegmentIterator>::value_type;
+    using PointType   = typename boost::polygon::segment_traits<SegmentType>::point_type;
+
+    const size_t          segment_idx  = (color >> 1) - 1;
+    const SegmentIterator segment_it   = segment_begin + segment_idx;
+    const PointType       source_point = boost::polygon::segment_traits<SegmentType>::get(*segment_it, ((color & 1) ? boost::polygon::HIGH :
+                                                                                                                      boost::polygon::LOW));
+    return source_point;
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    VoronoiDiagram::IssueType>::type
+VoronoiDiagram::try_to_repair_degenerated_voronoi_diagram_by_rotation(const SegmentIterator segment_begin,
+                                                                      const SegmentIterator segment_end,
+                                                                      const double          fix_angle)
+{
+    using SegmentType = typename std::iterator_traits<SegmentIterator>::value_type;
+    using PointType   = typename boost::polygon::segment_traits<SegmentType>::point_type;
+
+    // Copy all segments and rotate their vertices.
+    std::vector<VoronoiDiagram::Segment> segments_rotated;
+    segments_rotated.reserve(std::distance(segment_begin, segment_end));
+    for (auto segment_it = segment_begin; segment_it != segment_end; ++segment_it) {
+        PointType from = boost::polygon::segment_traits<SegmentType>::get(*segment_it, boost::polygon::LOW);
+        PointType to   = boost::polygon::segment_traits<SegmentType>::get(*segment_it, boost::polygon::HIGH);
+        segments_rotated.emplace_back(from.rotated(fix_angle), to.rotated(fix_angle));
+    }
+
+    VoronoiDiagram::voronoi_diagram_type voronoi_diagram_rotated;
+    boost::polygon::construct_voronoi(segments_rotated.begin(), segments_rotated.end(), &voronoi_diagram_rotated);
+
+    this->copy_to_local(voronoi_diagram_rotated);
+    const IssueType issue_type = detect_known_issues(*this, segments_rotated.begin(), segments_rotated.end());
+
+    // We want to remap all Voronoi vertices at the endpoints of input segments
+    // to ensure that Voronoi vertices at endpoints will be preserved after rotation.
+    // So we assign every Voronoi vertices color to map this Vertex into input segments.
+    for (cell_type cell : m_cells) {
+        if (cell.is_degenerate())
+            continue;
+
+        if (cell.contains_segment()) {
+            if (const SegmentCellRange cell_range = VoronoiUtils::compute_segment_cell_range(cell, segments_rotated.begin(), segments_rotated.end()); cell_range.is_valid()) {
+                if (cell_range.edge_end->vertex1()->color() == 0) {
+                    // Vertex 1 of edge_end points to the starting endpoint of the input segment (from() or line.a).
+                    VD::vertex_type::color_type color = encode_input_segment_endpoint(cell.source_index(), boost::polygon::LOW);
+                    cell_range.edge_end->vertex1()->color(color);
+                }
+
+                if (cell_range.edge_begin->vertex0()->color() == 0) {
+                    // Vertex 0 of edge_end points to the ending endpoint of the input segment (to() or line.b).
+                    VD::vertex_type::color_type color = encode_input_segment_endpoint(cell.source_index(), boost::polygon::HIGH);
+                    cell_range.edge_begin->vertex0()->color(color);
+                }
+            } else {
+                // This could happen when there is a missing Voronoi vertex even after rotation.
+                assert(cell_range.is_valid());
+            }
+        }
+
+        // FIXME @hejllukas: Implement mapping also for source points and not just for source segments.
+    }
+
+    // Rotate all Voronoi vertices back.
+    // When a Voronoi vertex can be mapped to the input segment endpoint, then we don't need to do rotation back.
+    for (vertex_type &vertex : m_vertices) {
+        if (vertex.color() == 0) {
+            // This vertex isn't mapped to any vertex, so we rotate it back.
+            vertex = VoronoiUtils::make_rotated_vertex(vertex, -fix_angle);
+        } else {
+            // This vertex can be mapped to the input segment endpoint.
+            PointType   endpoint = decode_input_segment_endpoint(vertex.color(), segment_begin, segment_end);
+            vertex_type endpoint_vertex{double(endpoint.x()), double(endpoint.y())};
+            endpoint_vertex.incident_edge(vertex.incident_edge());
+            endpoint_vertex.color(vertex.color());
+            vertex = endpoint_vertex;
+        }
+    }
+
+    // We have to clear all marked vertices because some algorithms expect that all vertices have a color equal to 0.
+    for (vertex_type &vertex : m_vertices)
+        vertex.color(0);
+
+    m_voronoi_diagram.clear();
+    m_is_modified = true;
+
+    return issue_type;
+}
+
+} // namespace Slic3r::Geometry
diff --git a/src/libslic3r/Geometry/Voronoi.hpp b/src/libslic3r/Geometry/Voronoi.hpp
index 5529750f3..ef1c62276 100644
--- a/src/libslic3r/Geometry/Voronoi.hpp
+++ b/src/libslic3r/Geometry/Voronoi.hpp
@@ -4,10 +4,8 @@
 #include "../Line.hpp"
 #include "../Polyline.hpp"
 
-#define BOOST_VORONOI_USE_GMP 1
-
 #ifdef _MSC_VER
-// Suppress warning C4146 in OpenVDB: unary minus operator applied to unsigned type, result still unsigned 
+// Suppress warning C4146 in OpenVDB: unary minus operator applied to unsigned type, result still unsigned
 #pragma warning(push)
 #pragma warning(disable : 4146)
 #endif // _MSC_VER
@@ -16,18 +14,182 @@
 #pragma warning(pop)
 #endif // _MSC_VER
 
-namespace Slic3r { 
+namespace Slic3r::Geometry {
 
-namespace Geometry {
-
-class VoronoiDiagram : public boost::polygon::voronoi_diagram<double> {
+class VoronoiDiagram
+{
 public:
-    typedef double                                          coord_type;
-    typedef boost::polygon::point_data<coordinate_type>     point_type;
-    typedef boost::polygon::segment_data<coordinate_type>   segment_type;
-    typedef boost::polygon::rectangle_data<coordinate_type> rect_type;
+    using coord_type           = double;
+    using voronoi_diagram_type = boost::polygon::voronoi_diagram<coord_type>;
+    using point_type           = boost::polygon::point_data<voronoi_diagram_type::coordinate_type>;
+    using segment_type         = boost::polygon::segment_data<voronoi_diagram_type::coordinate_type>;
+    using rect_type            = boost::polygon::rectangle_data<voronoi_diagram_type::coordinate_type>;
+
+    using coordinate_type = voronoi_diagram_type::coordinate_type;
+    using vertex_type     = voronoi_diagram_type::vertex_type;
+    using edge_type       = voronoi_diagram_type::edge_type;
+    using cell_type       = voronoi_diagram_type::cell_type;
+
+    using const_vertex_iterator = voronoi_diagram_type::const_vertex_iterator;
+    using const_edge_iterator   = voronoi_diagram_type::const_edge_iterator;
+    using const_cell_iterator   = voronoi_diagram_type::const_cell_iterator;
+
+    using vertex_container_type = voronoi_diagram_type::vertex_container_type;
+    using edge_container_type   = voronoi_diagram_type::edge_container_type;
+    using cell_container_type   = voronoi_diagram_type::cell_container_type;
+
+    enum class IssueType {
+        NO_ISSUE_DETECTED,
+        FINITE_EDGE_WITH_NON_FINITE_VERTEX,
+        MISSING_VORONOI_VERTEX,
+        NON_PLANAR_VORONOI_DIAGRAM,
+        VORONOI_EDGE_INTERSECTING_INPUT_SEGMENT,
+        UNKNOWN                                  // Repairs are disabled in the constructor.
+    };
+
+    enum class State {
+        REPAIR_NOT_NEEDED,   // The original Voronoi diagram doesn't have any issue.
+        REPAIR_SUCCESSFUL,   // The original Voronoi diagram has some issues, but it was repaired.
+        REPAIR_UNSUCCESSFUL, // The original Voronoi diagram has some issues, but it wasn't repaired.
+        UNKNOWN              // Repairs are disabled in the constructor.
+    };
+
+    VoronoiDiagram() = default;
+
+    virtual ~VoronoiDiagram() = default;
+
+    IssueType get_issue_type() const { return m_issue_type; }
+
+    State get_state() const { return m_state; }
+
+    bool is_valid() const { return m_state != State::REPAIR_UNSUCCESSFUL; }
+
+    void clear();
+
+    const vertex_container_type &vertices() const { return m_is_modified ? m_vertices : m_voronoi_diagram.vertices(); }
+
+    const edge_container_type &edges() const { return m_is_modified ? m_edges : m_voronoi_diagram.edges(); }
+
+    const cell_container_type &cells() const { return m_is_modified ? m_cells : m_voronoi_diagram.cells(); }
+
+    std::size_t num_vertices() const { return m_is_modified ? m_vertices.size() : m_voronoi_diagram.num_vertices(); }
+
+    std::size_t num_edges() const { return m_is_modified ? m_edges.size() : m_voronoi_diagram.num_edges(); }
+
+    std::size_t num_cells() const { return m_is_modified ? m_cells.size() : m_voronoi_diagram.num_cells(); }
+
+    template<typename SegmentIterator>
+    typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        void>::type
+    construct_voronoi(SegmentIterator segment_begin, SegmentIterator segment_end, bool try_to_repair_if_needed = true);
+
+    template<typename PointIterator>
+    typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_point_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<PointIterator>::value_type>::type>::type>::type,
+        void>::type
+    construct_voronoi(const PointIterator first, const PointIterator last)
+    {
+        boost::polygon::construct_voronoi(first, last, &m_voronoi_diagram);
+        m_state      = State::UNKNOWN;
+        m_issue_type = IssueType::UNKNOWN;
+    }
+
+    template<typename PointIterator, typename SegmentIterator>
+    typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_and<
+            typename boost::polygon::gtl_if<typename boost::polygon::is_point_concept<
+                typename boost::polygon::geometry_concept<typename std::iterator_traits<PointIterator>::value_type>::type>::type>::type,
+            typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<typename boost::polygon::geometry_concept<
+                typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type>::type,
+        void>::type
+    construct_voronoi(const PointIterator p_first, const PointIterator p_last, const SegmentIterator s_first, const SegmentIterator s_last)
+    {
+        boost::polygon::construct_voronoi(p_first, p_last, s_first, s_last, &m_voronoi_diagram);
+        m_state      = State::UNKNOWN;
+        m_issue_type = IssueType::UNKNOWN;
+    }
+
+    // Try to detect cases when some Voronoi vertex is missing, when the Voronoi diagram
+    // is not planar or some Voronoi edge is intersecting input segment.
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        IssueType>::type
+    detect_known_issues(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    template<typename SegmentIterator>
+    typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        VoronoiDiagram::IssueType>::type
+    try_to_repair_degenerated_voronoi_diagram_by_rotation(SegmentIterator segment_begin, SegmentIterator segment_end, double fix_angle);
+
+    template<typename SegmentIterator>
+    typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        VoronoiDiagram::IssueType>::type
+    try_to_repair_degenerated_voronoi_diagram(SegmentIterator segment_begin, SegmentIterator segment_end);
+
+private:
+    struct Segment
+    {
+        Point from;
+        Point to;
+
+        Segment() = delete;
+        explicit Segment(const Point &from, const Point &to) : from(from), to(to) {}
+    };
+
+    void copy_to_local(voronoi_diagram_type &voronoi_diagram);
+
+    // Detect issues related to Voronoi cells, or that can be detected by iterating over Voronoi cells.
+    // The first type of issue that can be detected is a missing Voronoi vertex, especially when it is
+    // missing at one of the endpoints of the input segment.
+    // The second type of issue that can be detected is a Voronoi edge that intersects the input segment.
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        IssueType>::type
+    detect_known_voronoi_cell_issues(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    static bool has_finite_edge_with_non_finite_vertex(const VoronoiDiagram &voronoi_diagram);
+
+    voronoi_diagram_type  m_voronoi_diagram;
+    vertex_container_type m_vertices;
+    edge_container_type   m_edges;
+    cell_container_type   m_cells;
+    bool                  m_is_modified = false;
+    State                 m_state       = State::UNKNOWN;
+    IssueType             m_issue_type  = IssueType::UNKNOWN;
+
+public:
+    using SegmentIt = std::vector<Slic3r::Geometry::VoronoiDiagram::Segment>::iterator;
+
+    friend struct boost::polygon::segment_traits<Slic3r::Geometry::VoronoiDiagram::Segment>;
 };
 
-} } // namespace Slicer::Geometry
+} // namespace Slic3r::Geometry
+
+namespace boost::polygon {
+template<> struct geometry_concept<Slic3r::Geometry::VoronoiDiagram::Segment>
+{
+    typedef segment_concept type;
+};
+
+template<> struct segment_traits<Slic3r::Geometry::VoronoiDiagram::Segment>
+{
+    using coordinate_type = coord_t;
+    using point_type      = Slic3r::Point;
+    using segment_type    = Slic3r::Geometry::VoronoiDiagram::Segment;
+
+    static inline point_type get(const segment_type &segment, direction_1d dir) { return dir.to_int() ? segment.to : segment.from; }
+};
+} // namespace boost::polygon
 
 #endif // slic3r_Geometry_Voronoi_hpp_
diff --git a/src/libslic3r/Geometry/VoronoiOffset.cpp b/src/libslic3r/Geometry/VoronoiOffset.cpp
index ed8d9c2b1..46105220a 100644
--- a/src/libslic3r/Geometry/VoronoiOffset.cpp
+++ b/src/libslic3r/Geometry/VoronoiOffset.cpp
@@ -782,9 +782,6 @@ void annotate_inside_outside(VD &vd, const Lines &lines)
 
     for (const VD::edge_type &edge : vd.edges())
         if (edge.vertex1() == nullptr) {
-            if (edge.vertex0() == nullptr)
-                continue;
-
             // Infinite Voronoi edge separating two Point sites or a Point site and a Segment site.
             // Infinite edge is always outside and it references at least one valid vertex.
             assert(edge.is_infinite());
@@ -891,9 +888,6 @@ void annotate_inside_outside(VD &vd, const Lines &lines)
     for (const VD::edge_type &edge : vd.edges()) {
         assert((edge_category(edge) == EdgeCategory::Unknown) == (edge_category(edge.twin()) == EdgeCategory::Unknown));
         if (edge_category(edge) == EdgeCategory::Unknown) {
-            if (!edge.is_finite())
-                continue;
-
             assert(edge.is_finite());
             const VD::cell_type &cell   = *edge.cell();
             const VD::cell_type &cell2  = *edge.twin()->cell();
diff --git a/src/libslic3r/Geometry/VoronoiUtils.cpp b/src/libslic3r/Geometry/VoronoiUtils.cpp
new file mode 100644
index 000000000..1e9436325
--- /dev/null
+++ b/src/libslic3r/Geometry/VoronoiUtils.cpp
@@ -0,0 +1,283 @@
+#include <boost/log/trivial.hpp>
+
+#include <Arachne/utils/PolygonsSegmentIndex.hpp>
+#include <MultiMaterialSegmentation.hpp>
+
+#include "VoronoiUtils.hpp"
+
+namespace Slic3r::Geometry {
+
+using PolygonsSegmentIndexConstIt = std::vector<Arachne::PolygonsSegmentIndex>::const_iterator;
+using LinesIt                     = Lines::iterator;
+using ColoredLinesIt              = ColoredLines::iterator;
+using ColoredLinesConstIt         = ColoredLines::const_iterator;
+
+// Explicit template instantiation.
+template LinesIt::reference VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &, LinesIt, LinesIt);
+template VD::SegmentIt::reference VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &, VD::SegmentIt, VD::SegmentIt);
+template ColoredLinesIt::reference VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &, ColoredLinesIt, ColoredLinesIt);
+template ColoredLinesConstIt::reference VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &, ColoredLinesConstIt, ColoredLinesConstIt);
+template PolygonsSegmentIndexConstIt::reference VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+template Point VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &, LinesIt, LinesIt);
+template Point VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &, VD::SegmentIt, VD::SegmentIt);
+template Point VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &, ColoredLinesIt, ColoredLinesIt);
+template Point VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &, ColoredLinesConstIt, ColoredLinesConstIt);
+template Point VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, LinesIt, LinesIt);
+template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, VD::SegmentIt, VD::SegmentIt);
+template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, ColoredLinesConstIt, ColoredLinesConstIt);
+template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+template Points VoronoiUtils::discretize_parabola(const Point &, const Arachne::PolygonsSegmentIndex &, const Point &, const Point &, coord_t, float);
+template Arachne::PolygonsPointIndex VoronoiUtils::get_source_point_index(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    typename std::iterator_traits<SegmentIterator>::reference>::type
+VoronoiUtils::get_source_segment(const VoronoiDiagram::cell_type &cell, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    if (!cell.contains_segment())
+        throw Slic3r::InvalidArgument("Voronoi cell doesn't contain a source segment!");
+
+    if (cell.source_index() >= size_t(std::distance(segment_begin, segment_end)))
+        throw Slic3r::OutOfRange("Voronoi cell source index is out of range!");
+
+    return *(segment_begin + cell.source_index());
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type
+VoronoiUtils::get_source_point(const VoronoiDiagram::cell_type &cell, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    using Segment = typename std::iterator_traits<SegmentIterator>::value_type;
+
+    if (!cell.contains_point())
+        throw Slic3r::InvalidArgument("Voronoi cell doesn't contain a source point!");
+
+    if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) {
+        assert(int(cell.source_index()) < std::distance(segment_begin, segment_end));
+        const SegmentIterator segment_it = segment_begin + cell.source_index();
+        return boost::polygon::segment_traits<Segment>::get(*segment_it, boost::polygon::LOW);
+    } else if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT) {
+        assert(int(cell.source_index()) < std::distance(segment_begin, segment_end));
+        const SegmentIterator segment_it = segment_begin + cell.source_index();
+        return boost::polygon::segment_traits<Segment>::get(*segment_it, boost::polygon::HIGH);
+    } else if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT) {
+        throw Slic3r::RuntimeError("Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!");
+    } else {
+        throw Slic3r::InvalidArgument("Function get_source_point() should only be called on point cells!");
+    }
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    Arachne::PolygonsPointIndex>::type
+VoronoiUtils::get_source_point_index(const VD::cell_type &cell, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    if (!cell.contains_point())
+        throw Slic3r::InvalidArgument("Voronoi cell doesn't contain a source point!");
+
+    if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) {
+        assert(int(cell.source_index()) < std::distance(segment_begin, segment_end));
+        const SegmentIterator segment_it = segment_begin + cell.source_index();
+        return (*segment_it);
+    } else if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT) {
+        assert(int(cell.source_index()) < std::distance(segment_begin, segment_end));
+        const SegmentIterator segment_it = segment_begin + cell.source_index();
+        return (*segment_it).next();
+    } else if (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT) {
+        throw Slic3r::RuntimeError("Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!");
+    } else {
+        throw Slic3r::InvalidArgument("Function get_source_point_index() should only be called on point cells!");
+    }
+}
+
+template<typename Segment>
+typename boost::polygon::enable_if<typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+    typename boost::polygon::geometry_concept<Segment>::type>::type>::type,
+    Points>::type
+VoronoiUtils::discretize_parabola(const Point &source_point, const Segment &source_segment, const Point &start, const Point &end, const coord_t approximate_step_size, float transitioning_angle)
+{
+    Points discretized;
+    // x is distance of point projected on the segment ab
+    // xx is point projected on the segment ab
+    const Point   a       = source_segment.from();
+    const Point   b       = source_segment.to();
+    const Point   ab      = b - a;
+    const Point   as      = start - a;
+    const Point   ae      = end - a;
+    const coord_t ab_size = ab.cast<int64_t>().norm();
+    const coord_t sx      = as.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+    const coord_t ex      = ae.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+    const coord_t sxex    = ex - sx;
+
+    const Point   ap = source_point - a;
+    const coord_t px = ap.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+
+    Point pxx;
+    Line(a, b).distance_to_infinite_squared(source_point, &pxx);
+    const Point   ppxx = pxx - source_point;
+    const coord_t d    = ppxx.cast<int64_t>().norm();
+
+    const Vec2d  rot           = perp(ppxx).cast<double>().normalized();
+    const double rot_cos_theta = rot.x();
+    const double rot_sin_theta = rot.y();
+
+    if (d == 0) {
+        discretized.emplace_back(start);
+        discretized.emplace_back(end);
+        return discretized;
+    }
+
+    const double marking_bound = atan(transitioning_angle * 0.5);
+    int64_t      msx           = -marking_bound * int64_t(d); // projected marking_start
+    int64_t      mex           = marking_bound * int64_t(d);  // projected marking_end
+
+    const coord_t marking_start_end_h = msx * msx / (2 * d) + d / 2;
+    Point         marking_start       = Point(coord_t(msx), marking_start_end_h).rotated(rot_cos_theta, rot_sin_theta) + pxx;
+    Point         marking_end         = Point(coord_t(mex), marking_start_end_h).rotated(rot_cos_theta, rot_sin_theta) + pxx;
+    const int     dir                 = (sx > ex) ? -1 : 1;
+    if (dir < 0) {
+        std::swap(marking_start, marking_end);
+        std::swap(msx, mex);
+    }
+
+    bool add_marking_start = msx * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && msx * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
+    bool add_marking_end   = mex * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && mex * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
+
+    const Point apex     = Point(0, d / 2).rotated(rot_cos_theta, rot_sin_theta) + pxx;
+    bool        add_apex = int64_t(sx - px) * int64_t(dir) < 0 && int64_t(ex - px) * int64_t(dir) > 0;
+
+    assert(!add_marking_start || !add_marking_end || add_apex);
+    if (add_marking_start && add_marking_end && !add_apex)
+        BOOST_LOG_TRIVIAL(warning) << "Failing to discretize parabola! Must add an apex or one of the endpoints.";
+
+    const coord_t step_count = lround(static_cast<double>(std::abs(ex - sx)) / approximate_step_size);
+    discretized.emplace_back(start);
+    for (coord_t step = 1; step < step_count; ++step) {
+        const int64_t x = int64_t(sx) + int64_t(sxex) * int64_t(step) / int64_t(step_count) - int64_t(px);
+        const int64_t y = int64_t(x) * int64_t(x) / int64_t(2 * d) + int64_t(d / 2);
+
+        if (add_marking_start && msx * int64_t(dir) < int64_t(x) * int64_t(dir)) {
+            discretized.emplace_back(marking_start);
+            add_marking_start = false;
+        }
+
+        if (add_apex && int64_t(x) * int64_t(dir) > 0) {
+            discretized.emplace_back(apex);
+            add_apex = false; // only add the apex just before the
+        }
+
+        if (add_marking_end && mex * int64_t(dir) < int64_t(x) * int64_t(dir)) {
+            discretized.emplace_back(marking_end);
+            add_marking_end = false;
+        }
+
+        assert(is_in_range<coord_t>(x) && is_in_range<coord_t>(y));
+        const Point result = Point(x, y).rotated(rot_cos_theta, rot_sin_theta) + pxx;
+        discretized.emplace_back(result);
+    }
+
+    if (add_apex)
+        discretized.emplace_back(apex);
+
+    if (add_marking_end)
+        discretized.emplace_back(marking_end);
+
+    discretized.emplace_back(end);
+    return discretized;
+}
+
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    Geometry::SegmentCellRange<
+        typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>>::type
+VoronoiUtils::compute_segment_cell_range(const VD::cell_type &cell, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    using Segment          = typename std::iterator_traits<SegmentIterator>::value_type;
+    using Point            = typename boost::polygon::segment_point_type<Segment>::type;
+    using SegmentCellRange = SegmentCellRange<Point>;
+
+    const Segment &source_segment = Geometry::VoronoiUtils::get_source_segment(cell, segment_begin, segment_end);
+    const Point    from           = boost::polygon::segment_traits<Segment>::get(source_segment, boost::polygon::LOW);
+    const Point    to             = boost::polygon::segment_traits<Segment>::get(source_segment, boost::polygon::HIGH);
+    const Vec2i64  from_i64       = from.template cast<int64_t>();
+    const Vec2i64  to_i64         = to.template cast<int64_t>();
+
+    // FIXME @hejllukas: Ensure that there is no infinite edge during iteration between edge_begin and edge_end.
+    SegmentCellRange cell_range(to, from);
+
+    // Find starting edge and end edge
+    bool                 seen_possible_start             = false;
+    bool                 after_start                     = false;
+    bool                 ending_edge_is_set_before_start = false;
+    const VD::edge_type *edge                            = cell.incident_edge();
+    do {
+        if (edge->is_infinite())
+            continue;
+
+        Vec2i64 v0 = Geometry::VoronoiUtils::to_point(edge->vertex0());
+        Vec2i64 v1 = Geometry::VoronoiUtils::to_point(edge->vertex1());
+        assert(v0 != to_i64 || v1 != from_i64);
+
+        if (v0 == to_i64 && !after_start) { // Use the last edge which starts in source_segment.to
+            cell_range.edge_begin = edge;
+            seen_possible_start   = true;
+        } else if (seen_possible_start) {
+            after_start = true;
+        }
+
+        if (v1 == from_i64 && (!cell_range.edge_end || ending_edge_is_set_before_start)) {
+            ending_edge_is_set_before_start = !after_start;
+            cell_range.edge_end             = edge;
+        }
+    } while (edge = edge->next(), edge != cell.incident_edge());
+
+    return cell_range;
+}
+
+Vec2i64 VoronoiUtils::to_point(const VD::vertex_type *vertex)
+{
+    assert(vertex != nullptr);
+    return VoronoiUtils::to_point(*vertex);
+}
+
+Vec2i64 VoronoiUtils::to_point(const VD::vertex_type &vertex)
+{
+    const double x = vertex.x(), y = vertex.y();
+
+    assert(std::isfinite(x) && std::isfinite(y));
+    assert(is_in_range<int64_t>(x) && is_in_range<int64_t>(y));
+
+    return {std::llround(x), std::llround(y)};
+}
+
+bool VoronoiUtils::is_finite(const VD::vertex_type &vertex)
+{
+    return std::isfinite(vertex.x()) && std::isfinite(vertex.y());
+}
+
+VD::vertex_type VoronoiUtils::make_rotated_vertex(VD::vertex_type &vertex, const double angle)
+{
+    const double cos_a = std::cos(angle);
+    const double sin_a = std::sin(angle);
+
+    const double rotated_x = (cos_a * vertex.x() - sin_a * vertex.y());
+    const double rotated_y = (cos_a * vertex.y() + sin_a * vertex.x());
+
+    VD::vertex_type rotated_vertex{rotated_x, rotated_y};
+    rotated_vertex.incident_edge(vertex.incident_edge());
+    rotated_vertex.color(vertex.color());
+
+    return rotated_vertex;
+}
+
+} // namespace Slic3r::Geometry
diff --git a/src/libslic3r/Geometry/VoronoiUtils.hpp b/src/libslic3r/Geometry/VoronoiUtils.hpp
new file mode 100644
index 000000000..bf6391467
--- /dev/null
+++ b/src/libslic3r/Geometry/VoronoiUtils.hpp
@@ -0,0 +1,120 @@
+#ifndef slic3r_VoronoiUtils_hpp_
+#define slic3r_VoronoiUtils_hpp_
+
+#include "libslic3r/Geometry/Voronoi.hpp"
+#include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
+
+using VD = Slic3r::Geometry::VoronoiDiagram;
+
+namespace Slic3r::Geometry {
+
+// Represent trapezoid Voronoi cell around segment.
+template<typename PT> struct SegmentCellRange
+{
+    const PT             segment_start_point;  // The start point of the source segment of this cell.
+    const PT             segment_end_point;    // The end point of the source segment of this cell.
+    const VD::edge_type *edge_begin = nullptr; // The edge of the Voronoi diagram where the loop around the cell starts.
+    const VD::edge_type *edge_end   = nullptr; // The edge of the Voronoi diagram where the loop around the cell ends.
+
+    SegmentCellRange() = delete;
+    explicit SegmentCellRange(const PT &segment_start_point, const PT &segment_end_point)
+        : segment_start_point(segment_start_point), segment_end_point(segment_end_point)
+    {}
+
+    bool is_valid() const { return edge_begin && edge_end && edge_begin != edge_end; }
+};
+
+class VoronoiUtils
+{
+public:
+    static Vec2i64 to_point(const VD::vertex_type *vertex);
+
+    static Vec2i64 to_point(const VD::vertex_type &vertex);
+
+    static bool is_finite(const VD::vertex_type &vertex);
+
+    static VD::vertex_type make_rotated_vertex(VD::vertex_type &vertex, double angle);
+
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        typename std::iterator_traits<SegmentIterator>::reference>::type
+    get_source_segment(const VD::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type
+    get_source_point(const VoronoiDiagram::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        Arachne::PolygonsPointIndex>::type
+    get_source_point_index(const VD::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    /**
+     * Discretize a parabola based on (approximate) step size.
+     *
+     * Adapted from CuraEngine VoronoiUtils::discretizeParabola by Tim Kuipers @BagelOrb and @Ghostkeeper.
+     *
+     * @param approximate_step_size is measured parallel to the source_segment, not along the parabola.
+     */
+    template<typename Segment>
+    static typename boost::polygon::enable_if<typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<Segment>::type>::type>::type,
+        Points>::type
+    discretize_parabola(const Point &source_point, const Segment &source_segment, const Point &start, const Point &end, coord_t approximate_step_size, float transitioning_angle);
+
+    /**
+     * Compute the range of line segments that surround a cell of the skeletal
+     * graph that belongs to a line segment of the medial axis.
+     *
+     * This should only be used on cells that belong to a central line segment
+     * of the skeletal graph, e.g. trapezoid cells, not triangular cells.
+     *
+     * The resulting line segments is just the first and the last segment. They
+     * are linked to the neighboring segments, so you can iterate over the
+     * segments until you reach the last segment.
+     *
+     * Adapted from CuraEngine VoronoiUtils::computePointCellRange by Tim Kuipers @BagelOrb,
+     * Jaime van Kessel @nallath, Remco Burema @rburema and @Ghostkeeper.
+     *
+     * @param cell The cell to compute the range of line segments for.
+     * @param segment_begin Begin iterator for all edges of the input Polygons.
+     * @param segment_end End iterator for all edges of the input Polygons.
+     * @return Range of line segments that surround the cell.
+     */
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        Geometry::SegmentCellRange<
+            typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>>::type
+    compute_segment_cell_range(const VD::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
+
+    template<typename T> static bool is_in_range(double value)
+    {
+        return double(std::numeric_limits<T>::lowest()) <= value && value <= double(std::numeric_limits<T>::max());
+    }
+
+    template<typename T> static bool is_in_range(const VD::vertex_type &vertex)
+    {
+        return VoronoiUtils::is_finite(vertex) && is_in_range<T>(vertex.x()) && is_in_range<T>(vertex.y());
+    }
+
+    template<typename T> static bool is_in_range(const VD::edge_type &edge)
+    {
+        if (edge.vertex0() == nullptr || edge.vertex1() == nullptr)
+            return false;
+
+        return is_in_range<T>(*edge.vertex0()) && is_in_range<T>(*edge.vertex1());
+    }
+};
+
+} // namespace Slic3r::Geometry
+
+#endif // slic3r_VoronoiUtils_hpp_
diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp
index 062a3b397..5817b33ec 100644
--- a/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp
+++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp
@@ -3,7 +3,9 @@
 #include <CGAL/Surface_sweep_2_algorithms.h>
 
 #include "libslic3r/Geometry/Voronoi.hpp"
-#include "libslic3r/Arachne/utils/VoronoiUtils.hpp"
+#include "libslic3r/Geometry/VoronoiUtils.hpp"
+#include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
+#include "libslic3r/MultiMaterialSegmentation.hpp"
 
 #include "VoronoiUtilsCgal.hpp"
 
@@ -11,18 +13,145 @@ using VD = Slic3r::Geometry::VoronoiDiagram;
 
 namespace Slic3r::Geometry {
 
-using CGAL_Point   = CGAL::Exact_predicates_exact_constructions_kernel::Point_2;
-using CGAL_Segment = CGAL::Arr_segment_traits_2<CGAL::Exact_predicates_exact_constructions_kernel>::Curve_2;
+using PolygonsSegmentIndexConstIt = std::vector<Arachne::PolygonsSegmentIndex>::const_iterator;
+using LinesIt                     = Lines::iterator;
+using ColoredLinesConstIt         = ColoredLines::const_iterator;
 
-inline static CGAL_Point to_cgal_point(const VD::vertex_type &pt) { return {pt.x(), pt.y()}; }
+// Explicit template instantiation.
+template bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VD &, LinesIt, LinesIt);
+template bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VD &, VD::SegmentIt, VD::SegmentIt);
+template bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VD &, ColoredLinesConstIt, ColoredLinesConstIt);
+template bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VD &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+
+// The tangent vector of the parabola is computed based on the Proof of the reflective property.
+// https://en.wikipedia.org/wiki/Parabola#Proof_of_the_reflective_property
+// https://math.stackexchange.com/q/2439647/2439663#comment5039739_2439663
+namespace impl {
+    using K   = CGAL::Simple_cartesian<double>;
+    using FK  = CGAL::Simple_cartesian<CGAL::Interval_nt_advanced>;
+    using EK  = CGAL::Simple_cartesian<CGAL::MP_Float>;
+    using C2E = CGAL::Cartesian_converter<K, EK>;
+    using C2F = CGAL::Cartesian_converter<K, FK>;
+    class Epick : public CGAL::Filtered_kernel_adaptor<CGAL::Type_equality_wrapper<K::Base<Epick>::Type, Epick>, true> {};
+
+    template<typename K>
+    inline typename K::Vector_2 calculate_parabolic_tangent_vector(
+        // Test point on the parabola, where the tangent will be calculated.
+        const typename K::Point_2 &p,
+        // Focus point of the parabola.
+        const typename K::Point_2 &f,
+        // Points of a directrix of the parabola.
+        const typename K::Point_2 &u,
+        const typename K::Point_2 &v,
+        // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent.
+        const typename K::Orientation &tangent_orientation)
+    {
+        using RT       = typename K::RT;
+        using Vector_2 = typename K::Vector_2;
+
+        const Vector_2 directrix_vec            = v - u;
+        const RT       directrix_vec_sqr_length = CGAL::scalar_product(directrix_vec, directrix_vec);
+        Vector_2       focus_vec                = (f - u) * directrix_vec_sqr_length - directrix_vec * CGAL::scalar_product(directrix_vec, p - u);
+        Vector_2       tangent_vec              = focus_vec.perpendicular(tangent_orientation);
+        return tangent_vec;
+    }
+
+    template<typename K> struct ParabolicTangentToSegmentOrientationPredicate
+    {
+        using Point_2     = typename K::Point_2;
+        using Vector_2    = typename K::Vector_2;
+        using Orientation = typename K::Orientation;
+        using result_type = typename K::Orientation;
+
+        result_type operator()(
+            // Test point on the parabola, where the tangent will be calculated.
+            const Point_2 &p,
+            // End of the linear segment (p, q), for which orientation towards the tangent to parabola will be evaluated.
+            const Point_2 &q,
+            // Focus point of the parabola.
+            const Point_2 &f,
+            // Points of a directrix of the parabola.
+            const Point_2 &u,
+            const Point_2 &v,
+            // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent.
+            const Orientation &tangent_orientation) const
+        {
+            assert(tangent_orientation == CGAL::Orientation::LEFT_TURN || tangent_orientation == CGAL::Orientation::RIGHT_TURN);
+
+            Vector_2 tangent_vec = calculate_parabolic_tangent_vector<K>(p, f, u, v, tangent_orientation);
+            Vector_2 linear_vec  = q - p;
+
+            return CGAL::sign(tangent_vec.x() * linear_vec.y() - tangent_vec.y() * linear_vec.x());
+        }
+    };
+
+    template<typename K> struct ParabolicTangentToParabolicTangentOrientationPredicate
+    {
+        using Point_2     = typename K::Point_2;
+        using Vector_2    = typename K::Vector_2;
+        using Orientation = typename K::Orientation;
+        using result_type = typename K::Orientation;
+
+        result_type operator()(
+            // Common point on both parabolas, where the tangent will be calculated.
+            const Point_2 &p,
+            // Focus point of the first parabola.
+            const Point_2 &f_0,
+            // Points of a directrix of the first parabola.
+            const Point_2 &u_0,
+            const Point_2 &v_0,
+            // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent.
+            const Orientation &tangent_orientation_0,
+            // Focus point of the second parabola.
+            const Point_2 &f_1,
+            // Points of a directrix of the second parabola.
+            const Point_2 &u_1,
+            const Point_2 &v_1,
+            // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent.
+            const Orientation &tangent_orientation_1) const
+        {
+            assert(tangent_orientation_0 == CGAL::Orientation::LEFT_TURN || tangent_orientation_0 == CGAL::Orientation::RIGHT_TURN);
+            assert(tangent_orientation_1 == CGAL::Orientation::LEFT_TURN || tangent_orientation_1 == CGAL::Orientation::RIGHT_TURN);
+
+            Vector_2 tangent_vec_0 = calculate_parabolic_tangent_vector<K>(p, f_0, u_0, v_0, tangent_orientation_0);
+            Vector_2 tangent_vec_1 = calculate_parabolic_tangent_vector<K>(p, f_1, u_1, v_1, tangent_orientation_1);
+
+            return CGAL::sign(tangent_vec_0.x() * tangent_vec_1.y() - tangent_vec_0.y() * tangent_vec_1.x());
+        }
+    };
+
+    using ParabolicTangentToSegmentOrientationPredicateFiltered = CGAL::Filtered_predicate<ParabolicTangentToSegmentOrientationPredicate<EK>, ParabolicTangentToSegmentOrientationPredicate<FK>, C2E, C2F>;
+    using ParabolicTangentToParabolicTangentOrientationPredicateFiltered = CGAL::Filtered_predicate<ParabolicTangentToParabolicTangentOrientationPredicate<EK>, ParabolicTangentToParabolicTangentOrientationPredicate<FK>, C2E, C2F>;
+} // namespace impl
+
+using ParabolicTangentToSegmentOrientation = impl::ParabolicTangentToSegmentOrientationPredicateFiltered;
+using ParabolicTangentToParabolicTangentOrientation = impl::ParabolicTangentToParabolicTangentOrientationPredicateFiltered;
+using CGAL_Point   = impl::K::Point_2;
+
+inline CGAL_Point to_cgal_point(const VD::vertex_type *pt) { return {pt->x(), pt->y()}; }
+inline CGAL_Point to_cgal_point(const Point &pt) { return {pt.x(), pt.y()}; }
+inline CGAL_Point to_cgal_point(const Vec2d &pt) { return {pt.x(), pt.y()}; }
+
+inline Linef make_linef(const VD::edge_type &edge)
+{
+    const VD::vertex_type *v0 = edge.vertex0();
+    const VD::vertex_type *v1 = edge.vertex1();
+    return {Vec2d(v0->x(), v0->y()), Vec2d(v1->x(), v1->y())};
+}
+
+[[maybe_unused]] inline bool is_equal(const VD::vertex_type &vertex_first, const VD::vertex_type &vertex_second) { return vertex_first.x() == vertex_second.x() && vertex_first.y() == vertex_second.y(); }
 
 // FIXME Lukas H.: Also includes parabolic segments.
 bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_diagram)
 {
+    using CGAL_E_Point   = CGAL::Exact_predicates_exact_constructions_kernel::Point_2;
+    using CGAL_E_Segment = CGAL::Arr_segment_traits_2<CGAL::Exact_predicates_exact_constructions_kernel>::Curve_2;
+    auto to_cgal_point   = [](const VD::vertex_type &pt) -> CGAL_E_Point { return {pt.x(), pt.y()}; };
+
     assert(std::all_of(voronoi_diagram.edges().cbegin(), voronoi_diagram.edges().cend(),
                        [](const VD::edge_type &edge) { return edge.color() == 0; }));
 
-    std::vector<CGAL_Segment> segments;
+    std::vector<CGAL_E_Segment> segments;
     segments.reserve(voronoi_diagram.num_edges());
 
     for (const VD::edge_type &edge : voronoi_diagram.edges()) {
@@ -30,7 +159,7 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_
             continue;
 
         if (edge.is_finite() && edge.is_linear() && edge.vertex0() != nullptr && edge.vertex1() != nullptr &&
-            Arachne::VoronoiUtils::is_finite(*edge.vertex0()) && Arachne::VoronoiUtils::is_finite(*edge.vertex1())) {
+            VoronoiUtils::is_finite(*edge.vertex0()) && VoronoiUtils::is_finite(*edge.vertex1())) {
             segments.emplace_back(to_cgal_point(*edge.vertex0()), to_cgal_point(*edge.vertex1()));
             edge.color(1);
             assert(edge.twin() != nullptr);
@@ -41,42 +170,136 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_
     for (const VD::edge_type &edge : voronoi_diagram.edges())
         edge.color(0);
 
-    std::vector<CGAL_Point> intersections_pt;
+    std::vector<CGAL_E_Point> intersections_pt;
     CGAL::compute_intersection_points(segments.begin(), segments.end(), std::back_inserter(intersections_pt));
     return intersections_pt.empty();
 }
 
-static bool check_if_three_vectors_are_ccw(const CGAL_Point &common_pt, const CGAL_Point &pt_1, const CGAL_Point &pt_2, const CGAL_Point &test_pt) {
-    CGAL::Orientation orientation = CGAL::orientation(common_pt, pt_1, pt_2);
+struct ParabolicSegment
+{
+    const Point focus;
+    const Line  directrix;
+    // Two points on the parabola;
+    const Linef segment;
+    // Indicate if focus point is on the left side or right side relative to parabolic segment endpoints.
+    const CGAL::Orientation is_focus_on_left;
+};
+
+template<typename SegmentIterator>
+inline static typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    ParabolicSegment>::type
+get_parabolic_segment(const VD::edge_type &edge, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    using Segment = typename std::iterator_traits<SegmentIterator>::value_type;
+    assert(edge.is_curved());
+
+    const VD::cell_type *left_cell  = edge.cell();
+    const VD::cell_type *right_cell = edge.twin()->cell();
+
+    const Point       focus_pt   = VoronoiUtils::get_source_point(*(left_cell->contains_point() ? left_cell : right_cell), segment_begin, segment_end);
+    const Segment    &directrix  = VoronoiUtils::get_source_segment(*(left_cell->contains_point() ? right_cell : left_cell), segment_begin, segment_end);
+    CGAL::Orientation focus_side = CGAL::opposite(CGAL::orientation(to_cgal_point(edge.vertex0()), to_cgal_point(edge.vertex1()), to_cgal_point(focus_pt)));
+
+    assert(focus_side == CGAL::Orientation::LEFT_TURN || focus_side == CGAL::Orientation::RIGHT_TURN);
+
+    const Point directrix_from = boost::polygon::segment_traits<Segment>::get(directrix, boost::polygon::LOW);
+    const Point directrix_to   = boost::polygon::segment_traits<Segment>::get(directrix, boost::polygon::HIGH);
+    return {focus_pt, Line(directrix_from, directrix_to), make_linef(edge), focus_side};
+}
+
+template<typename SegmentIterator>
+inline static typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    CGAL::Orientation>::type
+orientation_of_two_edges(const VD::edge_type &edge_a, const VD::edge_type &edge_b, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    assert(is_equal(*edge_a.vertex0(), *edge_b.vertex0()));
+    CGAL::Orientation orientation;
+    if (edge_a.is_linear() && edge_b.is_linear()) {
+        orientation = CGAL::orientation(to_cgal_point(edge_a.vertex0()), to_cgal_point(edge_a.vertex1()), to_cgal_point(edge_b.vertex1()));
+    } else if (edge_a.is_curved() && edge_b.is_curved()) {
+        const ParabolicSegment parabolic_a = get_parabolic_segment(edge_a, segment_begin, segment_end);
+        const ParabolicSegment parabolic_b = get_parabolic_segment(edge_b, segment_begin, segment_end);
+        orientation = ParabolicTangentToParabolicTangentOrientation{}(to_cgal_point(parabolic_a.segment.a),
+                                                                      to_cgal_point(parabolic_a.focus),
+                                                                      to_cgal_point(parabolic_a.directrix.a),
+                                                                      to_cgal_point(parabolic_a.directrix.b),
+                                                                      parabolic_a.is_focus_on_left,
+                                                                      to_cgal_point(parabolic_b.focus),
+                                                                      to_cgal_point(parabolic_b.directrix.a),
+                                                                      to_cgal_point(parabolic_b.directrix.b),
+                                                                      parabolic_b.is_focus_on_left);
+        return orientation;
+    } else {
+        assert(edge_a.is_curved() != edge_b.is_curved());
+
+        const VD::edge_type   &linear_edge    = edge_a.is_curved() ? edge_b : edge_a;
+        const VD::edge_type   &parabolic_edge = edge_a.is_curved() ? edge_a : edge_b;
+        const ParabolicSegment parabolic      = get_parabolic_segment(parabolic_edge, segment_begin, segment_end);
+        orientation = ParabolicTangentToSegmentOrientation{}(to_cgal_point(parabolic.segment.a), to_cgal_point(linear_edge.vertex1()),
+                                                             to_cgal_point(parabolic.focus),
+                                                             to_cgal_point(parabolic.directrix.a),
+                                                             to_cgal_point(parabolic.directrix.b),
+                                                             parabolic.is_focus_on_left);
+
+        if (edge_b.is_curved())
+            orientation = CGAL::opposite(orientation);
+    }
+
+    return orientation;
+}
+
+template<typename SegmentIterator>
+static typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    bool>::type
+check_if_three_edges_are_ccw(const VD::edge_type  &edge_first,
+                             const VD::edge_type  &edge_second,
+                             const VD::edge_type  &edge_third,
+                             const SegmentIterator segment_begin,
+                             const SegmentIterator segment_end)
+{
+    assert(is_equal(*edge_first.vertex0(), *edge_second.vertex0()) && is_equal(*edge_second.vertex0(), *edge_third.vertex0()));
+
+    CGAL::Orientation orientation = orientation_of_two_edges(edge_first, edge_second, segment_begin, segment_end);
     if (orientation == CGAL::Orientation::COLLINEAR) {
         // The first two edges are collinear, so the third edge must be on the right side on the first of them.
-        return CGAL::orientation(common_pt, pt_1, test_pt) == CGAL::Orientation::RIGHT_TURN;
+        return orientation_of_two_edges(edge_first, edge_third, segment_begin, segment_end) == CGAL::Orientation::RIGHT_TURN;
     } else if (orientation == CGAL::Orientation::LEFT_TURN) {
         // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is bellow PI.
         // So we need to check if test_pt isn't between them.
-        CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt);
-        CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt);
+        CGAL::Orientation orientation1 = orientation_of_two_edges(edge_first, edge_third, segment_begin, segment_end);
+        CGAL::Orientation orientation2 = orientation_of_two_edges(edge_second, edge_third, segment_begin, segment_end);
         return (orientation1 != CGAL::Orientation::LEFT_TURN || orientation2 != CGAL::Orientation::RIGHT_TURN);
     } else {
         assert(orientation == CGAL::Orientation::RIGHT_TURN);
         // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is upper PI.
         // So we need to check if test_pt is between them.
-        CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt);
-        CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt);
+        CGAL::Orientation orientation1 = orientation_of_two_edges(edge_first, edge_third, segment_begin, segment_end);
+        CGAL::Orientation orientation2 = orientation_of_two_edges(edge_second, edge_third, segment_begin, segment_end);
         return (orientation1 == CGAL::Orientation::RIGHT_TURN || orientation2 == CGAL::Orientation::LEFT_TURN);
     }
 }
 
-bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram)
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    bool>::type
+VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VD             &voronoi_diagram,
+                                                  const SegmentIterator segment_begin,
+                                                  const SegmentIterator segment_end)
 {
     for (const VD::vertex_type &vertex : voronoi_diagram.vertices()) {
         std::vector<const VD::edge_type *> edges;
         const VD::edge_type               *edge = vertex.incident_edge();
 
         do {
-            // FIXME Lukas H.: Also process parabolic segments.
-            if (edge->is_finite() && edge->is_linear() && edge->vertex0() != nullptr && edge->vertex1() != nullptr &&
-                Arachne::VoronoiUtils::is_finite(*edge->vertex0()) && Arachne::VoronoiUtils::is_finite(*edge->vertex1()))
+            if (edge->is_finite() && edge->vertex0() != nullptr && edge->vertex1() != nullptr && VoronoiUtils::is_finite(*edge->vertex0()) && VoronoiUtils::is_finite(*edge->vertex1()))
                 edges.emplace_back(edge);
 
             edge = edge->rot_next();
@@ -85,12 +308,11 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &vor
         // Checking for CCW make sense for three and more edges.
         if (edges.size() > 2) {
             for (auto edge_it = edges.begin() ; edge_it != edges.end(); ++edge_it) {
-                const Geometry::VoronoiDiagram::edge_type *prev_edge = edge_it == edges.begin() ? edges.back() : *std::prev(edge_it);
-                const Geometry::VoronoiDiagram::edge_type *curr_edge = *edge_it;
-                const Geometry::VoronoiDiagram::edge_type *next_edge = std::next(edge_it) == edges.end() ? edges.front() : *std::next(edge_it);
+                const VD::edge_type *prev_edge = edge_it == edges.begin() ? edges.back() : *std::prev(edge_it);
+                const VD::edge_type *curr_edge = *edge_it;
+                const VD::edge_type *next_edge = std::next(edge_it) == edges.end() ? edges.front() : *std::next(edge_it);
 
-                if (!check_if_three_vectors_are_ccw(to_cgal_point(*prev_edge->vertex0()), to_cgal_point(*prev_edge->vertex1()),
-                                                    to_cgal_point(*curr_edge->vertex1()), to_cgal_point(*next_edge->vertex1())))
+                if (!check_if_three_edges_are_ccw(*prev_edge, *curr_edge, *next_edge, segment_begin, segment_end))
                     return false;
             }
         }
@@ -99,5 +321,4 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &vor
     return true;
 }
 
-
 } // namespace Slic3r::Geometry
\ No newline at end of file
diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp
index 897891bd9..33ae84780 100644
--- a/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp
+++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp
@@ -2,6 +2,7 @@
 #define slic3r_VoronoiUtilsCgal_hpp_
 
 #include "Voronoi.hpp"
+#include "../Arachne/utils/PolygonsSegmentIndex.hpp"
 
 namespace Slic3r::Geometry {
 class VoronoiDiagram;
@@ -13,8 +14,12 @@ public:
     static bool is_voronoi_diagram_planar_intersection(const VoronoiDiagram &voronoi_diagram);
 
     // Check if the Voronoi diagram is planar using verification that all neighboring edges are ordered CCW for each vertex.
-    static bool is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram);
-
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        bool>::type
+    is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end);
 };
 } // namespace Slic3r::Geometry
 
diff --git a/src/libslic3r/MultiMaterialSegmentation.cpp b/src/libslic3r/MultiMaterialSegmentation.cpp
index a8bf9e6ff..45f464cf9 100644
--- a/src/libslic3r/MultiMaterialSegmentation.cpp
+++ b/src/libslic3r/MultiMaterialSegmentation.cpp
@@ -4,11 +4,11 @@
 #include "Layer.hpp"
 #include "Print.hpp"
 #include "Geometry/VoronoiVisualUtils.hpp"
+#include "Geometry/VoronoiUtils.hpp"
 #include "MutablePolygon.hpp"
 #include "format.hpp"
 
 #include <utility>
-#include <cfloat>
 #include <unordered_set>
 
 #include <boost/log/trivial.hpp>
@@ -16,36 +16,19 @@
 #include <mutex>
 #include <boost/thread/lock_guard.hpp>
 
-namespace Slic3r {
-struct ColoredLine {
-    Line line;
-    int color;
-    int poly_idx = -1;
-    int local_line_idx = -1;
-};
-}
+//#define MM_SEGMENTATION_DEBUG_GRAPH
+//#define MM_SEGMENTATION_DEBUG_REGIONS
+//#define MM_SEGMENTATION_DEBUG_INPUT
+//#define MM_SEGMENTATION_DEBUG_PAINTED_LINES
+//#define MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
 
-#include <boost/polygon/polygon.hpp>
-namespace boost::polygon {
-template <>
-struct geometry_concept<Slic3r::ColoredLine> { typedef segment_concept type; };
+#if defined(MM_SEGMENTATION_DEBUG_GRAPH) || defined(MM_SEGMENTATION_DEBUG_REGIONS) || \
+    defined(MM_SEGMENTATION_DEBUG_INPUT) || defined(MM_SEGMENTATION_DEBUG_PAINTED_LINES) || \
+    defined(MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS)
+#define MM_SEGMENTATION_DEBUG
+#endif
 
-template <>
-struct segment_traits<Slic3r::ColoredLine> {
-    typedef coord_t coordinate_type;
-    typedef Slic3r::Point point_type;
-
-    static inline point_type get(const Slic3r::ColoredLine& line, const direction_1d& dir) {
-        return dir.to_int() ? line.line.b : line.line.a;
-    }
-};
-}
-
-//#define MMU_SEGMENTATION_DEBUG_GRAPH
-//#define MMU_SEGMENTATION_DEBUG_REGIONS
-//#define MMU_SEGMENTATION_DEBUG_INPUT
-//#define MMU_SEGMENTATION_DEBUG_PAINTED_LINES
-//#define MMU_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
+//#define MM_SEGMENTATION_DEBUG_TOP_BOTTOM
 
 namespace Slic3r {
 bool is_equal(float left, float right, float eps = 1e-3) {
@@ -156,55 +139,31 @@ struct PaintedLineVisitor
     static inline const double                                                            append_threshold2 = Slic3r::sqr(append_threshold);
 };
 
-static Polygon colored_points_to_polygon(const std::vector<ColoredLine> &lines)
-{
-    Polygon out;
-    out.points.reserve(lines.size());
-    for (const ColoredLine &l : lines)
-        out.points.emplace_back(l.line.a);
-    return out;
-}
-
-static Polygons colored_points_to_polygon(const std::vector<std::vector<ColoredLine>> &lines)
-{
-    Polygons out;
-    out.reserve(lines.size());
-    for (const std::vector<ColoredLine> &l : lines)
-        out.emplace_back(colored_points_to_polygon(l));
-    return out;
+BoundingBox get_extents(const std::vector<ColoredLines> &colored_polygons) {
+    BoundingBox bbox;
+    for (const ColoredLines &colored_lines : colored_polygons) {
+        for (const ColoredLine &colored_line : colored_lines) {
+            bbox.merge(colored_line.line.a);
+            bbox.merge(colored_line.line.b);
+        }
+    }
+    return bbox;
 }
 
 // Flatten the vector of vectors into a vector.
-static inline std::vector<ColoredLine> to_lines(const std::vector<std::vector<ColoredLine>> &c_lines)
+static inline ColoredLines to_lines(const std::vector<ColoredLines> &c_lines)
 {
     size_t n_lines = 0;
     for (const auto &c_line : c_lines)
         n_lines += c_line.size();
-    std::vector<ColoredLine> lines;
+    ColoredLines lines;
     lines.reserve(n_lines);
     for (const auto &c_line : c_lines)
         lines.insert(lines.end(), c_line.begin(), c_line.end());
     return lines;
 }
 
-static bool vertex_equal_to_point(const Voronoi::VD::vertex_type &vertex, const Vec2d &ipt)
-{
-    // Convert ipt to doubles, force the 80bit FPU temporary to 64bit and then compare.
-    // This should work with any settings of math compiler switches and the C++ compiler
-    // shall understand the memcpies as type punning and it shall optimize them out.
-    using ulp_cmp_type = boost::polygon::detail::ulp_comparison<double>;
-    ulp_cmp_type ulp_cmp;
-    static constexpr int ULPS = boost::polygon::voronoi_diagram_traits<double>::vertex_equality_predicate_type::ULPS;
-    return ulp_cmp(vertex.x(), ipt.x(), ULPS) == ulp_cmp_type::EQUAL &&
-           ulp_cmp(vertex.y(), ipt.y(), ULPS) == ulp_cmp_type::EQUAL;
-}
-
-static inline bool vertex_equal_to_point(const Voronoi::VD::vertex_type *vertex, const Vec2d &ipt)
-{
-    return vertex_equal_to_point(*vertex, ipt);
-}
-
-static std::vector<std::pair<size_t, size_t>> get_segments(const std::vector<ColoredLine> &polygon)
+static std::vector<std::pair<size_t, size_t>> get_segments(const ColoredLines &polygon)
 {
     std::vector<std::pair<size_t, size_t>> segments;
 
@@ -229,16 +188,6 @@ static std::vector<std::pair<size_t, size_t>> get_segments(const std::vector<Col
     return segments;
 }
 
-static std::vector<std::vector<std::pair<size_t, size_t>>> get_all_segments(const std::vector<std::vector<ColoredLine>> &color_poly)
-{
-    std::vector<std::vector<std::pair<size_t, size_t>>> all_segments(color_poly.size());
-    for (size_t poly_idx = 0; poly_idx < color_poly.size(); ++poly_idx) {
-        const std::vector<ColoredLine> &c_polygon = color_poly[poly_idx];
-        all_segments[poly_idx]                    = get_segments(c_polygon);
-    }
-    return all_segments;
-}
-
 static std::vector<PaintedLine> filter_painted_lines(const Line &line_to_process, const size_t start_idx, const size_t end_idx, const std::vector<PaintedLine> &painted_lines)
 {
     const int                filter_eps_value = scale_(0.1f);
@@ -319,7 +268,7 @@ static std::vector<std::vector<PaintedLine>> post_process_painted_lines(const st
 }
 
 #ifndef NDEBUG
-static bool are_lines_connected(const std::vector<ColoredLine> &colored_lines)
+static bool are_lines_connected(const ColoredLines &colored_lines)
 {
     for (size_t line_idx = 1; line_idx < colored_lines.size(); ++line_idx)
         if (colored_lines[line_idx - 1].line.b != colored_lines[line_idx].line.a)
@@ -328,7 +277,7 @@ static bool are_lines_connected(const std::vector<ColoredLine> &colored_lines)
 }
 #endif
 
-static std::vector<ColoredLine> colorize_line(const Line &line_to_process,
+static ColoredLines colorize_line(const Line &line_to_process,
                                               const size_t              start_idx,
                                               const size_t              end_idx,
                                               const std::vector<PaintedLine> &painted_contour)
@@ -336,9 +285,9 @@ static std::vector<ColoredLine> colorize_line(const Line &line_to_process,
     assert(start_idx < painted_contour.size() && end_idx < painted_contour.size() && start_idx <= end_idx);
     assert(std::all_of(painted_contour.begin() + start_idx, painted_contour.begin() + end_idx + 1, [&painted_contour, &start_idx](const auto &p_line) { return painted_contour[start_idx].line_idx == p_line.line_idx; }));
 
-    const int                filter_eps_value = scale_(0.1f);
-    std::vector<ColoredLine> final_lines;
-    const PaintedLine       &first_line = painted_contour[start_idx];
+    const int          filter_eps_value = scale_(0.1f);
+    ColoredLines       final_lines;
+    const PaintedLine &first_line = painted_contour[start_idx];
     if (double dist_to_start = (first_line.projected_line.a - line_to_process.a).cast<double>().norm(); dist_to_start > filter_eps_value)
         final_lines.push_back({Line(line_to_process.a, first_line.projected_line.a), 0});
     final_lines.push_back({first_line.projected_line, first_line.color});
@@ -377,7 +326,7 @@ static std::vector<ColoredLine> colorize_line(const Line &line_to_process,
             if (line_1.line.length() <= scale_(0.2)) line_1.color = line_0.color;
     }
 
-    std::vector<ColoredLine> colored_lines_simple;
+    ColoredLines colored_lines_simple;
     colored_lines_simple.emplace_back(final_lines.front());
     for (size_t line_idx = 1; line_idx < final_lines.size(); ++line_idx) {
         const ColoredLine &line_0 = final_lines[line_idx];
@@ -405,7 +354,7 @@ static std::vector<ColoredLine> colorize_line(const Line &line_to_process,
     return final_lines;
 }
 
-static std::vector<ColoredLine> filter_colorized_polygon(std::vector<ColoredLine> &&new_lines) {
+static ColoredLines filter_colorized_polygon(ColoredLines &&new_lines) {
     for (size_t line_idx = 2; line_idx < new_lines.size(); ++line_idx) {
         const ColoredLine &line_0 = new_lines[line_idx - 2];
         ColoredLine       &line_1 = new_lines[line_idx - 1];
@@ -494,10 +443,10 @@ static std::vector<ColoredLine> filter_colorized_polygon(std::vector<ColoredLine
     return std::move(new_lines);
 }
 
-static std::vector<ColoredLine> colorize_contour(const EdgeGrid::Contour &contour, const std::vector<PaintedLine> &painted_contour) {
+static ColoredLines colorize_contour(const EdgeGrid::Contour &contour, const std::vector<PaintedLine> &painted_contour) {
     assert(painted_contour.empty() || std::all_of(painted_contour.begin(), painted_contour.end(), [&painted_contour](const auto &p_line) { return painted_contour.front().contour_idx == p_line.contour_idx; }));
 
-    std::vector<ColoredLine> colorized_contour;
+    ColoredLines colorized_contour;
     if (painted_contour.empty()) {
         // Appends contour with default color for lines before the first PaintedLine.
         colorized_contour.reserve(contour.num_segments());
@@ -534,297 +483,33 @@ static std::vector<ColoredLine> colorize_contour(const EdgeGrid::Contour &contou
     return filter_colorized_polygon(std::move(colorized_contour));
 }
 
-static std::vector<std::vector<ColoredLine>> colorize_contours(const std::vector<EdgeGrid::Contour> &contours, const std::vector<std::vector<PaintedLine>> &painted_contours)
+static std::vector<ColoredLines> colorize_contours(const std::vector<EdgeGrid::Contour> &contours, const std::vector<std::vector<PaintedLine>> &painted_contours)
 {
     assert(contours.size() == painted_contours.size());
-    std::vector<std::vector<ColoredLine>> colorized_contours(contours.size());
+    std::vector<ColoredLines> colorized_contours(contours.size());
     for (const std::vector<PaintedLine> &painted_contour : painted_contours) {
         size_t contour_idx              = &painted_contour - &painted_contours.front();
         colorized_contours[contour_idx] = colorize_contour(contours[contour_idx], painted_contours[contour_idx]);
     }
+
+    size_t poly_idx = 0;
+    for (ColoredLines &color_lines : colorized_contours) {
+        size_t line_idx = 0;
+        for (size_t color_line_idx = 0; color_line_idx < color_lines.size(); ++color_line_idx) {
+            color_lines[color_line_idx].poly_idx       = int(poly_idx);
+            color_lines[color_line_idx].local_line_idx = int(line_idx);
+            ++line_idx;
+        }
+        ++poly_idx;
+    }
+
     return colorized_contours;
 }
 
-using boost::polygon::voronoi_diagram;
-
-static inline Point mk_point(const Voronoi::VD::vertex_type *point) { return {coord_t(point->x()), coord_t(point->y())}; }
-
-static inline Point mk_point(const Voronoi::Internal::point_type &point) { return {coord_t(point.x()), coord_t(point.y())}; }
-
-static inline Point mk_point(const voronoi_diagram<double>::vertex_type &point) { return {coord_t(point.x()), coord_t(point.y())}; }
-
-static inline Point mk_point(const Vec2d &point) { return {coord_t(std::round(point.x())), coord_t(std::round(point.y()))}; }
-
-static inline Vec2d mk_vec2(const voronoi_diagram<double>::vertex_type *point) { return {point->x(), point->y()}; }
-
-struct MMU_Graph
-{
-    enum class ARC_TYPE { BORDER, NON_BORDER };
-
-    struct Arc
-    {
-        size_t   from_idx;
-        size_t   to_idx;
-        int      color;
-        ARC_TYPE type;
-
-        bool operator==(const Arc &rhs) const { return (from_idx == rhs.from_idx) && (to_idx == rhs.to_idx) && (color == rhs.color) && (type == rhs.type); }
-        bool operator!=(const Arc &rhs) const { return !operator==(rhs); }
-    };
-
-    struct Node
-    {
-        Vec2d             point;
-        std::list<size_t> arc_idxs;
-
-        void remove_edge(const size_t to_idx, MMU_Graph &graph)
-        {
-            for (auto arc_it = this->arc_idxs.begin(); arc_it != this->arc_idxs.end(); ++arc_it) {
-                MMU_Graph::Arc &arc = graph.arcs[*arc_it];
-                if (arc.to_idx == to_idx) {
-                    assert(arc.type != ARC_TYPE::BORDER);
-                    this->arc_idxs.erase(arc_it);
-                    break;
-                }
-            }
-        }
-    };
-
-    std::vector<MMU_Graph::Node> nodes;
-    std::vector<MMU_Graph::Arc>  arcs;
-    size_t                       all_border_points{};
-
-    std::vector<size_t> polygon_idx_offset;
-    std::vector<size_t> polygon_sizes;
-
-    void remove_edge(const size_t from_idx, const size_t to_idx)
-    {
-        nodes[from_idx].remove_edge(to_idx, *this);
-        nodes[to_idx].remove_edge(from_idx, *this);
-    }
-
-    [[nodiscard]] size_t get_global_index(const size_t poly_idx, const size_t point_idx) const { return polygon_idx_offset[poly_idx] + point_idx; }
-
-    void append_edge(const size_t &from_idx, const size_t &to_idx, int color = -1, ARC_TYPE type = ARC_TYPE::NON_BORDER)
-    {
-        // Don't append duplicate edges between the same nodes.
-        for (const size_t &arc_idx : this->nodes[from_idx].arc_idxs)
-            if (arcs[arc_idx].to_idx == to_idx)
-                return;
-        for (const size_t &arc_idx : this->nodes[to_idx].arc_idxs)
-            if (arcs[arc_idx].to_idx == from_idx)
-                return;
-
-        this->nodes[from_idx].arc_idxs.push_back(this->arcs.size());
-        this->arcs.push_back({from_idx, to_idx, color, type});
-
-        // Always insert only one directed arc for the input polygons.
-        // Two directed arcs in both directions are inserted if arcs aren't between points of the input polygons.
-        if (type == ARC_TYPE::NON_BORDER) {
-            this->nodes[to_idx].arc_idxs.push_back(this->arcs.size());
-            this->arcs.push_back({to_idx, from_idx, color, type});
-        }
-    }
-
-    // It assumes that between points of the input polygons is always only one directed arc,
-    // with the same direction as lines of the input polygon.
-    [[nodiscard]] MMU_Graph::Arc get_border_arc(size_t idx) const {
-        assert(idx < this->all_border_points);
-        return this->arcs[idx];
-    }
-
-    [[nodiscard]] size_t nodes_count() const { return this->nodes.size(); }
-
-    void remove_nodes_with_one_arc()
-    {
-        std::queue<size_t> update_queue;
-        for (const MMU_Graph::Node &node : this->nodes) {
-            size_t node_idx = &node - &this->nodes.front();
-            // Skip nodes that represent points of input polygons.
-            if (node.arc_idxs.size() == 1 && node_idx >= this->all_border_points)
-                update_queue.emplace(&node - &this->nodes.front());
-        }
-
-        while (!update_queue.empty()) {
-            size_t           node_from_idx = update_queue.front();
-            MMU_Graph::Node &node_from     = this->nodes[update_queue.front()];
-            update_queue.pop();
-            if (node_from.arc_idxs.empty())
-                continue;
-
-            assert(node_from.arc_idxs.size() == 1);
-            size_t           node_to_idx = arcs[node_from.arc_idxs.front()].to_idx;
-            MMU_Graph::Node &node_to     = this->nodes[node_to_idx];
-            this->remove_edge(node_from_idx, node_to_idx);
-            if (node_to.arc_idxs.size() == 1 && node_to_idx >= this->all_border_points)
-                update_queue.emplace(node_to_idx);
-        }
-    }
-
-    void add_contours(const std::vector<std::vector<ColoredLine>> &color_poly)
-    {
-        this->all_border_points = nodes.size();
-        this->polygon_sizes     = std::vector<size_t>(color_poly.size());
-        for (size_t polygon_idx = 0; polygon_idx < color_poly.size(); ++polygon_idx)
-            this->polygon_sizes[polygon_idx] = color_poly[polygon_idx].size();
-        this->polygon_idx_offset    = std::vector<size_t>(color_poly.size());
-        this->polygon_idx_offset[0] = 0;
-        for (size_t polygon_idx = 1; polygon_idx < color_poly.size(); ++polygon_idx) {
-            this->polygon_idx_offset[polygon_idx] = this->polygon_idx_offset[polygon_idx - 1] + color_poly[polygon_idx - 1].size();
-        }
-
-        size_t poly_idx = 0;
-        for (const std::vector<ColoredLine> &color_lines : color_poly) {
-            size_t line_idx = 0;
-            for (const ColoredLine &color_line : color_lines) {
-                size_t from_idx = this->get_global_index(poly_idx, line_idx);
-                size_t to_idx   = this->get_global_index(poly_idx, (line_idx + 1) % color_lines.size());
-                this->append_edge(from_idx, to_idx, color_line.color, ARC_TYPE::BORDER);
-                ++line_idx;
-            }
-            ++poly_idx;
-        }
-    }
-
-    // Nodes 0..all_border_points are only one with are on countour. Other vertexis are consider as not on coouter. So we check if base on attach index
-    inline bool is_vertex_on_contour(const Voronoi::VD::vertex_type *vertex) const
-    {
-        assert(vertex != nullptr);
-        return vertex->color() < this->all_border_points;
-    }
-
-    [[nodiscard]] inline bool is_edge_attach_to_contour(const voronoi_diagram<double>::const_edge_iterator &edge_iterator) const
-    {
-        return this->is_vertex_on_contour(edge_iterator->vertex0()) || this->is_vertex_on_contour(edge_iterator->vertex1());
-    }
-
-    [[nodiscard]] inline bool is_edge_connecting_two_contour_vertices(const voronoi_diagram<double>::const_edge_iterator &edge_iterator) const
-    {
-        return this->is_vertex_on_contour(edge_iterator->vertex0()) && this->is_vertex_on_contour(edge_iterator->vertex1());
-    }
-
-    // All Voronoi vertices are post-processes to merge very close vertices to single. Witch eliminates issues with intersection edges.
-    // Also, Voronoi vertices outside of the bounding of input polygons are throw away by marking them.
-    void append_voronoi_vertices(const Geometry::VoronoiDiagram &vd, const Polygons &color_poly_tmp, BoundingBox bbox) {
-        bbox.offset(SCALED_EPSILON);
-
-        struct CPoint
-        {
-            CPoint() = delete;
-            CPoint(const Vec2d &point, size_t contour_idx, size_t point_idx) : m_point_double(point), m_point(mk_point(point)), m_point_idx(point_idx), m_contour_idx(contour_idx) {}
-            CPoint(const Vec2d &point, size_t point_idx) : m_point_double(point), m_point(mk_point(point)), m_point_idx(point_idx), m_contour_idx(0) {}
-            const Vec2d m_point_double;
-            const Point m_point;
-            size_t      m_point_idx;
-            size_t      m_contour_idx;
-
-            [[nodiscard]] const Vec2d &point_double() const { return m_point_double; }
-            [[nodiscard]] const Point &point() const { return m_point; }
-            bool operator==(const CPoint &rhs) const { return this->m_point_double == rhs.m_point_double && this->m_contour_idx == rhs.m_contour_idx && this->m_point_idx == rhs.m_point_idx; }
-        };
-        struct CPointAccessor { const Point* operator()(const CPoint &pt) const { return &pt.point(); }};
-        typedef ClosestPointInRadiusLookup<CPoint, CPointAccessor> CPointLookupType;
-
-        CPointLookupType closest_voronoi_point(coord_t(SCALED_EPSILON));
-        CPointLookupType closest_contour_point(3 * coord_t(SCALED_EPSILON));
-        for (const Polygon &polygon : color_poly_tmp)
-            for (const Point &pt : polygon.points)
-                closest_contour_point.insert(CPoint(Vec2d(pt.x(), pt.y()), &polygon - &color_poly_tmp.front(), &pt - &polygon.points.front()));
-
-        for (const voronoi_diagram<double>::vertex_type &vertex : vd.vertices()) {
-            vertex.color(-1);
-            Vec2d vertex_point_double = Vec2d(vertex.x(), vertex.y());
-            Point vertex_point        = mk_point(vertex);
-
-            const Vec2d &first_point_double  = this->nodes[this->get_border_arc(vertex.incident_edge()->cell()->source_index()).from_idx].point;
-            const Vec2d &second_point_double = this->nodes[this->get_border_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx].point;
-
-            if (vertex_equal_to_point(&vertex, first_point_double)) {
-                assert(vertex.color() != vertex.incident_edge()->cell()->source_index());
-                assert(vertex.color() != vertex.incident_edge()->twin()->cell()->source_index());
-                vertex.color(this->get_border_arc(vertex.incident_edge()->cell()->source_index()).from_idx);
-            } else if (vertex_equal_to_point(&vertex, second_point_double)) {
-                assert(vertex.color() != vertex.incident_edge()->cell()->source_index());
-                assert(vertex.color() != vertex.incident_edge()->twin()->cell()->source_index());
-                vertex.color(this->get_border_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx);
-            } else if (bbox.contains(vertex_point)) {
-                if (auto [contour_pt, c_dist_sqr] = closest_contour_point.find(vertex_point); contour_pt != nullptr && c_dist_sqr < Slic3r::sqr(3 * SCALED_EPSILON)) {
-                    vertex.color(this->get_global_index(contour_pt->m_contour_idx, contour_pt->m_point_idx));
-                } else if (auto [voronoi_pt, v_dist_sqr] = closest_voronoi_point.find(vertex_point); voronoi_pt == nullptr || v_dist_sqr >= Slic3r::sqr(SCALED_EPSILON / 10.0)) {
-                    closest_voronoi_point.insert(CPoint(vertex_point_double, this->nodes_count()));
-                    vertex.color(this->nodes_count());
-                    this->nodes.push_back({vertex_point_double});
-                } else {
-                    // Boost Voronoi diagram generator sometimes creates two very closed points instead of one point.
-                    // For the example points (146872.99999999997, -146872.99999999997) and (146873, -146873), this example also included in Voronoi generator test cases.
-                    std::vector<std::pair<const CPoint *, double>> all_closes_c_points = closest_voronoi_point.find_all(vertex_point);
-                    int                                            merge_to_point      = -1;
-                    for (const std::pair<const CPoint *, double> &c_point : all_closes_c_points)
-                        if ((vertex_point_double - c_point.first->point_double()).squaredNorm() <= Slic3r::sqr(EPSILON)) {
-                            merge_to_point = int(c_point.first->m_point_idx);
-                            break;
-                        }
-
-                    if (merge_to_point != -1) {
-                        vertex.color(merge_to_point);
-                    } else {
-                        closest_voronoi_point.insert(CPoint(vertex_point_double, this->nodes_count()));
-                        vertex.color(this->nodes_count());
-                        this->nodes.push_back({vertex_point_double});
-                    }
-                }
-            }
-        }
-    }
-
-    void garbage_collect()
-    {
-        std::vector<int> nodes_map(this->nodes.size(), -1);
-        int              nodes_count = 0;
-        size_t           arcs_count  = 0;
-        for (const MMU_Graph::Node &node : this->nodes)
-            if (size_t node_idx = &node - &this->nodes.front(); !node.arc_idxs.empty()) {
-                nodes_map[node_idx] = nodes_count++;
-                arcs_count += node.arc_idxs.size();
-            }
-
-        std::vector<MMU_Graph::Node> new_nodes;
-        std::vector<MMU_Graph::Arc>  new_arcs;
-        new_nodes.reserve(nodes_count);
-        new_arcs.reserve(arcs_count);
-        for (const MMU_Graph::Node &node : this->nodes)
-            if (size_t node_idx = &node - &this->nodes.front(); nodes_map[node_idx] >= 0) {
-                new_nodes.push_back({node.point});
-                for (const size_t &arc_idx : node.arc_idxs) {
-                    const Arc &arc = this->arcs[arc_idx];
-                    new_nodes.back().arc_idxs.emplace_back(new_arcs.size());
-                    new_arcs.push_back({size_t(nodes_map[arc.from_idx]), size_t(nodes_map[arc.to_idx]), arc.color, arc.type});
-                }
-            }
-
-        this->nodes = std::move(new_nodes);
-        this->arcs  = std::move(new_arcs);
-    }
-};
-
-static inline void mark_processed(const voronoi_diagram<double>::const_edge_iterator &edge_iterator)
-{
-    edge_iterator->color(true);
-    edge_iterator->twin()->color(true);
-}
-
-// Return true, if "p" is closer to line.a, then line.b
-static inline bool is_point_closer_to_beginning_of_line(const Line &line, const Point &p)
-{
-    return (p - line.a).cast<double>().squaredNorm() < (p - line.b).cast<double>().squaredNorm();
-}
-
-static inline bool has_same_color(const ColoredLine &cl1, const ColoredLine &cl2) { return cl1.color == cl2.color; }
-
 // Determines if the line points from the point between two contour lines is pointing inside polygon or outside.
 static inline bool points_inside(const Line &contour_first, const Line &contour_second, const Point &new_point)
 {
-    // Used in points_inside for decision if line leading thought the common point of two lines is pointing inside polygon or outside
+    // TODO: Used in points_inside for decision if line leading thought the common point of two lines is pointing inside polygon or outside
     auto three_points_inward_normal = [](const Point &left, const Point &middle, const Point &right) -> Vec2d {
         assert(left != middle);
         assert(middle != right);
@@ -839,534 +524,309 @@ static inline bool points_inside(const Line &contour_first, const Line &contour_
     return side > 0.;
 }
 
-static inline bool line_intersection_with_epsilon(const Line &line_to_extend, const Line &other, Point *intersection)
-{
-    Line extended_line = line_to_extend;
-    extended_line.extend(15 * SCALED_EPSILON);
-    return extended_line.intersection(other, intersection);
+enum VD_ANNOTATION : Voronoi::VD::cell_type::color_type {
+    VERTEX_ON_CONTOUR = 1,
+    DELETED           = 2
+};
+
+#ifdef MM_SEGMENTATION_DEBUG_GRAPH
+static void export_graph_to_svg(const std::string &path, const Voronoi::VD& vd, const std::vector<ColoredLines>& colored_polygons) {
+    const coordf_t                 stroke_width = scaled<coordf_t>(0.05f);
+    const BoundingBox              bbox         = get_extents(colored_polygons);
+
+    SVG svg(path.c_str(), bbox);
+    for (const ColoredLines &colored_lines : colored_polygons)
+        for (const ColoredLine &colored_line : colored_lines)
+            svg.draw(colored_line.line, "black", stroke_width);
+
+    for (const Voronoi::VD::vertex_type &vertex : vd.vertices()) {
+        if (Geometry::VoronoiUtils::is_in_range<coord_t>(vertex)) {
+            if (const Point pt = Geometry::VoronoiUtils::to_point(&vertex).cast<coord_t>(); vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR) {
+                svg.draw(pt, "blue", coord_t(stroke_width));
+            } else if (vertex.color() != VD_ANNOTATION::DELETED) {
+                svg.draw(pt, "green", coord_t(stroke_width));
+            }
+        }
+    }
+
+    for (const Voronoi::VD::edge_type &edge : vd.edges()) {
+        if (edge.is_infinite() || !Geometry::VoronoiUtils::is_in_range<coord_t>(edge))
+            continue;
+
+        const Point from = Geometry::VoronoiUtils::to_point(edge.vertex0()).cast<coord_t>();
+        const Point to   = Geometry::VoronoiUtils::to_point(edge.vertex1()).cast<coord_t>();
+
+        if (edge.color() != VD_ANNOTATION::DELETED)
+            svg.draw(Line(from, to), "red", stroke_width);
+    }
+}
+#endif // MM_SEGMENTATION_DEBUG_GRAPH
+
+static size_t non_deleted_edge_count(const VD::vertex_type &vertex) {
+    size_t               non_deleted_edge_cnt = 0;
+    const VD::edge_type *edge                 = vertex.incident_edge();
+    do {
+        if (edge->color() != VD_ANNOTATION::DELETED)
+            ++non_deleted_edge_cnt;
+    } while (edge = edge->prev()->twin(), edge != vertex.incident_edge());
+
+    return non_deleted_edge_cnt;
 }
 
-// For every ColoredLine in lines_colored_out, assign the index of the polygon to which belongs and also the index of this line inside of the polygon.
-static inline void init_polygon_indices(const MMU_Graph                             &graph,
-                                        const std::vector<std::vector<ColoredLine>> &color_poly,
-                                        std::vector<ColoredLine>                    &lines_colored_out)
-{
-    size_t poly_idx = 0;
-    for (const std::vector<ColoredLine> &color_lines : color_poly) {
-        size_t line_idx = 0;
-        for (size_t color_line_idx = 0; color_line_idx < color_lines.size(); ++color_line_idx) {
-            size_t from_idx                            = graph.get_global_index(poly_idx, line_idx);
-            lines_colored_out[from_idx].poly_idx       = int(poly_idx);
-            lines_colored_out[from_idx].local_line_idx = int(line_idx);
-            ++line_idx;
-        }
-        ++poly_idx;
+static bool can_vertex_be_deleted(const VD::vertex_type &vertex) {
+    if (vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR || vertex.color() == VD_ANNOTATION::DELETED)
+        return false;
+
+    return non_deleted_edge_count(vertex) <= 1;
+}
+
+static void delete_vertex_deep(const VD::vertex_type &vertex) {
+    std::queue<const VD::vertex_type *> vertices_to_delete;
+    vertices_to_delete.emplace(&vertex);
+
+    while (!vertices_to_delete.empty()) {
+        const VD::vertex_type &vertex_to_delete = *vertices_to_delete.front();
+        vertices_to_delete.pop();
+        vertex_to_delete.color(VD_ANNOTATION::DELETED);
+
+        const VD::edge_type *edge = vertex_to_delete.incident_edge();
+        do {
+            edge->color(VD_ANNOTATION::DELETED);
+            edge->twin()->color(VD_ANNOTATION::DELETED);
+
+            if (edge->is_finite() && can_vertex_be_deleted(*edge->vertex1()))
+                vertices_to_delete.emplace(edge->vertex1());
+        } while (edge = edge->prev()->twin(), edge != vertex_to_delete.incident_edge());
     }
 }
 
-// Voronoi edges produced by Voronoi generator cloud have coordinates that don't fit inside coord_t (int32_t).
-// Because of that, this function tries to clip edges that have one endpoint of the edge inside the BoundingBox.
-static inline Line clip_finite_voronoi_edge(const Voronoi::VD::edge_type &edge, const BoundingBoxf &bbox)
-{
+static inline Vec2d mk_point_vec2d(const VD::vertex_type *point) {
+    assert(point != nullptr);
+    return {point->x(), point->y()};
+}
+
+static inline Vec2d mk_vector_vec2d(const VD::edge_type *edge) {
+    assert(edge != nullptr);
+    return mk_point_vec2d(edge->vertex1()) - mk_point_vec2d(edge->vertex0());
+}
+
+static inline Vec2d mk_flipped_vector_vec2d(const VD::edge_type *edge) {
+    assert(edge != nullptr);
+    return mk_point_vec2d(edge->vertex0()) - mk_point_vec2d(edge->vertex1());
+}
+
+static double edge_length(const VD::edge_type &edge) {
     assert(edge.is_finite());
-    Vec2d v0          = mk_vec2(edge.vertex0());
-    Vec2d v1          = mk_vec2(edge.vertex1());
-    bool  contains_v0 = bbox.contains(v0);
-    bool  contains_v1 = bbox.contains(v1);
-    if ((contains_v0 && contains_v1) || (!contains_v0 && !contains_v1))
-        return {mk_point(edge.vertex0()), mk_point(edge.vertex1())};
-
-    Vec2d vector = (v1 - v0).normalized() * bbox.size().norm();
-    if (!contains_v0)
-        v0 = (v1 - vector);
-    else
-        v1 = (v0 + vector);
-
-    return {v0.cast<coord_t>(), v1.cast<coord_t>()};
-}
-
-static MMU_Graph build_graph(size_t layer_idx, const std::vector<std::vector<ColoredLine>> &color_poly)
-{
-    Geometry::VoronoiDiagram vd;
-    std::vector<ColoredLine> lines_colored  = to_lines(color_poly);
-    const Polygons           color_poly_tmp = colored_points_to_polygon(color_poly);
-    const Points             points         = to_points(color_poly_tmp);
-    const Lines              lines          = to_lines(color_poly_tmp);
-
-    // The algorithm adds edges to the graph that are between two different colors.
-    // If a polygon is colored entirely with one color, we need to add at least one edge from that polygon artificially.
-    // Adding this edge is necessary for cases where the expolygon has an outer contour colored whole with one color
-    // and a hole colored with a different color. If an edge wasn't added to the graph,
-    // the entire expolygon would be colored with single random color instead of two different.
-    std::vector<bool>        force_edge_adding(color_poly.size());
-
-    // For each polygon, check if it is all colored with the same color. If it is, we need to force adding one edge to it.
-    for (const std::vector<ColoredLine> &c_poly : color_poly) {
-        bool force_edge = true;
-        for (const ColoredLine &c_line : c_poly)
-            if (c_line.color != c_poly.front().color) {
-                force_edge = false;
-                break;
-            }
-        force_edge_adding[&c_poly - &color_poly.front()] = force_edge;
-    }
-
-    boost::polygon::construct_voronoi(lines_colored.begin(), lines_colored.end(), &vd);
-    MMU_Graph graph;
-    graph.nodes.reserve(points.size() + vd.vertices().size());
-    for (const Point &point : points)
-        graph.nodes.push_back({Vec2d(double(point.x()), double(point.y()))});
-
-    graph.add_contours(color_poly);
-    init_polygon_indices(graph, color_poly, lines_colored);
-
-    assert(graph.nodes.size() == lines_colored.size());
-    BoundingBox bbox = get_extents(color_poly_tmp);
-    graph.append_voronoi_vertices(vd, color_poly_tmp, bbox);
-
-    auto get_prev_contour_line = [&lines_colored, &color_poly, &graph](const voronoi_diagram<double>::const_edge_iterator &edge_it) -> ColoredLine {
-        size_t contour_line_local_idx = lines_colored[edge_it->cell()->source_index()].local_line_idx;
-        size_t contour_line_size      = color_poly[lines_colored[edge_it->cell()->source_index()].poly_idx].size();
-        size_t contour_prev_idx       = graph.get_global_index(lines_colored[edge_it->cell()->source_index()].poly_idx,
-                                                         (contour_line_local_idx > 0) ? contour_line_local_idx - 1 : contour_line_size - 1);
-        return lines_colored[contour_prev_idx];
-    };
-
-    auto get_next_contour_line = [&lines_colored, &color_poly, &graph](const voronoi_diagram<double>::const_edge_iterator &edge_it) -> ColoredLine {
-        size_t contour_line_local_idx = lines_colored[edge_it->cell()->source_index()].local_line_idx;
-        size_t contour_line_size      = color_poly[lines_colored[edge_it->cell()->source_index()].poly_idx].size();
-        size_t contour_next_idx       = graph.get_global_index(lines_colored[edge_it->cell()->source_index()].poly_idx,
-                                                         (contour_line_local_idx + 1) % contour_line_size);
-        return lines_colored[contour_next_idx];
-    };
-
-    bbox.offset(scale_(10.));
-    const BoundingBoxf bbox_clip(bbox.min.cast<double>(), bbox.max.cast<double>());
-    const double       bbox_dim_max = double(std::max(bbox.size().x(), bbox.size().y()));
-
-    // Make a copy of the input segments with the double type.
-    std::vector<Voronoi::Internal::segment_type> segments;
-    for (const Line &line : lines)
-        segments.emplace_back(Voronoi::Internal::point_type(double(line.a(0)), double(line.a(1))),
-                              Voronoi::Internal::point_type(double(line.b(0)), double(line.b(1))));
-
-    for (auto edge_it = vd.edges().begin(); edge_it != vd.edges().end(); ++edge_it) {
-        // Skip second half-edge
-        if (edge_it->cell()->source_index() > edge_it->twin()->cell()->source_index() || edge_it->color())
-            continue;
-
-        if (edge_it->is_infinite() && (edge_it->vertex0() != nullptr || edge_it->vertex1() != nullptr)) {
-            // Infinite edge is leading through a point on the counter, but there are no Voronoi vertices.
-            // So we could fix this case by computing the intersection between the contour line and infinity edge.
-            std::vector<Voronoi::Internal::point_type> samples;
-            Voronoi::Internal::clip_infinite_edge(points, segments, *edge_it, bbox_dim_max, &samples);
-            if (samples.empty())
-                continue;
-
-            const Line         edge_line(mk_point(samples[0]), mk_point(samples[1]));
-            const ColoredLine &contour_line = lines_colored[edge_it->cell()->source_index()];
-            Point              contour_intersection;
-
-            if (line_intersection_with_epsilon(contour_line.line, edge_line, &contour_intersection)) {
-                const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_it->cell()->source_index());
-                const size_t          from_idx  = (edge_it->vertex1() != nullptr) ? edge_it->vertex1()->color() : edge_it->vertex0()->color();
-                size_t                to_idx    = ((contour_line.line.a - contour_intersection).cast<double>().squaredNorm() <
-                                 (contour_line.line.b - contour_intersection).cast<double>().squaredNorm()) ?
-                                                      graph_arc.from_idx :
-                                                      graph_arc.to_idx;
-                if (from_idx != to_idx && from_idx < graph.nodes_count() && to_idx < graph.nodes_count()) {
-                    graph.append_edge(from_idx, to_idx);
-                    mark_processed(edge_it);
-                }
-            }
-        } else if (edge_it->is_finite()) {
-            // Both points are on contour, so skip them. In cases of duplicate Voronoi vertices, skip edges between the same two points.
-            if (graph.is_edge_connecting_two_contour_vertices(edge_it) || (edge_it->vertex0()->color() == edge_it->vertex1()->color()))
-                continue;
-
-            const Line        edge_line         = clip_finite_voronoi_edge(*edge_it, bbox_clip);
-            const Line        contour_line      = lines_colored[edge_it->cell()->source_index()].line;
-            const ColoredLine colored_line      = lines_colored[edge_it->cell()->source_index()];
-            const ColoredLine contour_line_prev = get_prev_contour_line(edge_it);
-            const ColoredLine contour_line_next = get_next_contour_line(edge_it);
-
-            if (edge_it->vertex0()->color() >= graph.nodes_count() || edge_it->vertex1()->color() >= graph.nodes_count()) {
-                enum class Vertex { VERTEX0, VERTEX1 };
-                auto append_edge_if_intersects_with_contour = [&graph, &lines_colored, &edge_line, &contour_line](const voronoi_diagram<double>::const_edge_iterator &edge_iterator, const Vertex vertex) {
-                    Point intersection;
-                    Line  contour_line_twin = lines_colored[edge_iterator->twin()->cell()->source_index()].line;
-                    if (line_intersection_with_epsilon(contour_line_twin, edge_line, &intersection)) {
-                        const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_iterator->twin()->cell()->source_index());
-                        const size_t          to_idx_l  = is_point_closer_to_beginning_of_line(contour_line_twin, intersection) ? graph_arc.from_idx :
-                                                                                                                                  graph_arc.to_idx;
-                        graph.append_edge(vertex == Vertex::VERTEX0 ? edge_iterator->vertex0()->color() : edge_iterator->vertex1()->color(), to_idx_l);
-                    } else if (line_intersection_with_epsilon(contour_line, edge_line, &intersection)) {
-                        const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_iterator->cell()->source_index());
-                        const size_t to_idx_l = is_point_closer_to_beginning_of_line(contour_line, intersection) ? graph_arc.from_idx : graph_arc.to_idx;
-                        graph.append_edge(vertex == Vertex::VERTEX0 ? edge_iterator->vertex0()->color() : edge_iterator->vertex1()->color(), to_idx_l);
-                    }
-                    mark_processed(edge_iterator);
-                };
-
-                if (edge_it->vertex0()->color() < graph.nodes_count() && !graph.is_vertex_on_contour(edge_it->vertex0()))
-                    append_edge_if_intersects_with_contour(edge_it, Vertex::VERTEX0);
-
-                if (edge_it->vertex1()->color() < graph.nodes_count() && !graph.is_vertex_on_contour(edge_it->vertex1()))
-                    append_edge_if_intersects_with_contour(edge_it, Vertex::VERTEX1);
-            } else if (graph.is_edge_attach_to_contour(edge_it)) {
-                mark_processed(edge_it);
-                // Skip edges witch connection two points on a contour
-                if (graph.is_edge_connecting_two_contour_vertices(edge_it))
-                    continue;
-
-                const size_t from_idx = edge_it->vertex0()->color();
-                const size_t to_idx   = edge_it->vertex1()->color();
-                if (graph.is_vertex_on_contour(edge_it->vertex0())) {
-                    if (is_point_closer_to_beginning_of_line(contour_line, edge_line.a)) {
-                        if ((!has_same_color(contour_line_prev, colored_line) || force_edge_adding[colored_line.poly_idx]) && points_inside(contour_line_prev.line, contour_line, edge_line.b)) {
-                            graph.append_edge(from_idx, to_idx);
-                            force_edge_adding[colored_line.poly_idx] = false;
-                        }
-                    } else {
-                        if ((!has_same_color(contour_line_next, colored_line) || force_edge_adding[colored_line.poly_idx]) && points_inside(contour_line, contour_line_next.line, edge_line.b)) {
-                            graph.append_edge(from_idx, to_idx);
-                            force_edge_adding[colored_line.poly_idx] = false;
-                        }
-                    }
-                } else {
-                    assert(graph.is_vertex_on_contour(edge_it->vertex1()));
-                    if (is_point_closer_to_beginning_of_line(contour_line, edge_line.b)) {
-                        if ((!has_same_color(contour_line_prev, colored_line) || force_edge_adding[colored_line.poly_idx]) && points_inside(contour_line_prev.line, contour_line, edge_line.a)) {
-                            graph.append_edge(from_idx, to_idx);
-                            force_edge_adding[colored_line.poly_idx] = false;
-                        }
-                    } else {
-                        if ((!has_same_color(contour_line_next, colored_line) || force_edge_adding[colored_line.poly_idx]) && points_inside(contour_line, contour_line_next.line, edge_line.a)) {
-                            graph.append_edge(from_idx, to_idx);
-                            force_edge_adding[colored_line.poly_idx] = false;
-                        }
-                    }
-                }
-            } else if (Point intersection; line_intersection_with_epsilon(contour_line, edge_line, &intersection)) {
-                mark_processed(edge_it);
-                Vec2d real_v0_double = graph.nodes[edge_it->vertex0()->color()].point;
-                Vec2d real_v1_double = graph.nodes[edge_it->vertex1()->color()].point;
-                Point real_v0        = Point(coord_t(real_v0_double.x()), coord_t(real_v0_double.y()));
-                Point real_v1        = Point(coord_t(real_v1_double.x()), coord_t(real_v1_double.y()));
-
-                if (is_point_closer_to_beginning_of_line(contour_line, intersection)) {
-                    Line first_part(intersection, real_v0);
-                    Line second_part(intersection, real_v1);
-
-                    if (!has_same_color(contour_line_prev, colored_line)) {
-                        if (points_inside(contour_line_prev.line, contour_line, first_part.b))
-                            graph.append_edge(edge_it->vertex0()->color(), graph.get_border_arc(edge_it->cell()->source_index()).from_idx);
-
-                        if (points_inside(contour_line_prev.line, contour_line, second_part.b))
-                            graph.append_edge(edge_it->vertex1()->color(), graph.get_border_arc(edge_it->cell()->source_index()).from_idx);
-                    }
-                } else {
-                    const size_t int_point_idx    = graph.get_border_arc(edge_it->cell()->source_index()).to_idx;
-                    const Vec2d  int_point_double = graph.nodes[int_point_idx].point;
-                    const Point  int_point        = Point(coord_t(int_point_double.x()), coord_t(int_point_double.y()));
-
-                    const Line first_part(int_point, real_v0);
-                    const Line second_part(int_point, real_v1);
-
-                    if (!has_same_color(contour_line_next, colored_line)) {
-                        if (points_inside(contour_line, contour_line_next.line, first_part.b))
-                            graph.append_edge(edge_it->vertex0()->color(), int_point_idx);
-
-                        if (points_inside(contour_line, contour_line_next.line, second_part.b))
-                            graph.append_edge(edge_it->vertex1()->color(), int_point_idx);
-                    }
-                }
-            }
-        }
-    }
-
-    for (auto edge_it = vd.edges().begin(); edge_it != vd.edges().end(); ++edge_it) {
-        // Skip second half-edge and processed edges
-        if (edge_it->cell()->source_index() > edge_it->twin()->cell()->source_index() || edge_it->color())
-            continue;
-
-        if (edge_it->is_finite() && !bool(edge_it->color()) && edge_it->vertex0()->color() < graph.nodes_count() &&
-            edge_it->vertex1()->color() < graph.nodes_count()) {
-            // Skip cases, when the edge is between two same vertices, which is in cases two near vertices were merged together.
-            if (edge_it->vertex0()->color() == edge_it->vertex1()->color())
-                continue;
-
-            size_t from_idx = edge_it->vertex0()->color();
-            size_t to_idx   = edge_it->vertex1()->color();
-            graph.append_edge(from_idx, to_idx);
-        }
-        mark_processed(edge_it);
-    }
-
-    graph.remove_nodes_with_one_arc();
-    return graph;
-}
-
-static inline Polygon to_polygon(const std::vector<std::pair<size_t, Linef>> &id_to_lines)
-{
-    std::vector<Linef> lines;
-    for (auto id_to_line : id_to_lines)
-        lines.emplace_back(id_to_line.second);
-
-    Polygon poly_out;
-    poly_out.points.reserve(lines.size());
-    for (const Linef &line : lines)
-        poly_out.points.emplace_back(mk_point(line.a));
-    return poly_out;
-}
-
-
-static std::vector<std::vector<const MMU_Graph::Arc *>> get_all_next_arcs(
-    const MMU_Graph &graph,
-    std::vector<bool> &used_arcs,
-    const Linef &process_line,
-    const MMU_Graph::Arc &original_arc,
-    const int color)
-{
-    std::vector<std::vector<const MMU_Graph::Arc *>> all_next_arcs;
-    for (const size_t &arc_idx : graph.nodes[original_arc.to_idx].arc_idxs) {
-        std::vector<const MMU_Graph::Arc *> next_continue_arc;
-
-        const MMU_Graph::Arc &              arc = graph.arcs[arc_idx];
-        if (graph.nodes[arc.to_idx].point == process_line.a || used_arcs[arc_idx])
-            continue;
-
-        if (original_arc.type == MMU_Graph::ARC_TYPE::BORDER && original_arc.color != color)
-            continue;
-
-        if (arc.type == MMU_Graph::ARC_TYPE::BORDER && arc.color != color)
-            continue;
-
-        Vec2d arc_line = graph.nodes[arc.to_idx].point - graph.nodes[arc.from_idx].point;
-        next_continue_arc.emplace_back(&arc);
-        all_next_arcs.emplace_back(next_continue_arc);
-    }
-    return all_next_arcs;
-}
-
-// two points that are very close are considered as one point
-// std::vector contain the close points
-static std::vector<const MMU_Graph::Arc *> get_next_arc(
-    const MMU_Graph &graph,
-    std::vector<bool> &used_arcs,
-    const Linef &process_line,
-    const MMU_Graph::Arc &original_arc,
-    const int color)
-{
-    std::vector<const MMU_Graph::Arc *> res;
-
-    std::vector<std::vector<const MMU_Graph::Arc *>> all_next_arcs = get_all_next_arcs(graph, used_arcs, process_line, original_arc, color);
-    if (all_next_arcs.empty()) {
-        res.emplace_back(&original_arc);
-        return res;
-    }
-
-    std::vector<std::pair<std::vector<const MMU_Graph::Arc *>, double>> sorted_arcs;
-    for (auto next_arc : all_next_arcs) {
-        if (next_arc.empty())
-            continue;
-
-        Vec2d process_line_vec_n   = (process_line.a - process_line.b).normalized();
-        Vec2d neighbour_line_vec_n = (graph.nodes[next_arc.back()->to_idx].point - graph.nodes[next_arc.back()->from_idx].point).normalized();
-
-        double angle = ::acos(std::clamp(neighbour_line_vec_n.dot(process_line_vec_n), -1.0, 1.0));
-        if (Slic3r::cross2(neighbour_line_vec_n, process_line_vec_n) < 0.0)
-            angle = 2.0 * (double) PI - angle;
-
-        sorted_arcs.emplace_back(next_arc, angle);
-    }
-
-    std::sort(sorted_arcs.begin(), sorted_arcs.end(),
-        [](std::pair<std::vector<const MMU_Graph::Arc *>, double> &l, std::pair<std::vector<const MMU_Graph::Arc *>, double> &r) -> bool {
-            return l.second < r.second;
-        });
-
-    // Try to return left most edge witch is unused
-    for (auto &sorted_arc : sorted_arcs) {
-        if (size_t arc_idx = sorted_arc.first.back() - &graph.arcs.front(); !used_arcs[arc_idx])
-            return sorted_arc.first;
-    }
-
-    if (sorted_arcs.empty()) {
-        res.emplace_back(&original_arc);
-        return res;
-    }
-
-    return sorted_arcs.front().first;
-}
-
-static bool is_profile_self_interaction(Polygon poly)
-{
-    auto lines = poly.lines();
-    Point intersection;
-    for (int i = 0; i < lines.size(); ++i) {
-        for (int j = i + 2; j < std::min(lines.size(), lines.size() + i - 1); ++j) {
-            if (lines[i].intersection(lines[j], &intersection))
-                return true;
-        }
-    }
-    return false;
-}
-
-// Returns list of polygons and assigned colors.
-// It iterates through all nodes on the border between two different colors, and from this point,
-// start selection always left most edges for every node to construct CCW polygons.
-// Assumes that graph is planar (without self-intersection edges)
-static std::vector<ExPolygons> extract_colored_segments(const MMU_Graph &graph, const size_t num_extruders)
-{
-    std::vector<bool> used_arcs(graph.arcs.size(), false);
-    
-    auto all_arc_used = [&used_arcs](const MMU_Graph::Node &node) -> bool {
-        return std::all_of(node.arc_idxs.cbegin(), node.arc_idxs.cend(), [&used_arcs](const size_t &arc_idx) -> bool { return used_arcs[arc_idx]; });
-    };
-
-    std::vector<ExPolygons> expolygons_segments(num_extruders + 1);
-    for (size_t node_idx = 0; node_idx < graph.all_border_points; ++node_idx) {
-        const MMU_Graph::Node &node = graph.nodes[node_idx];
-
-        for (const size_t &arc_idx : node.arc_idxs) {
-            const MMU_Graph::Arc &arc = graph.arcs[arc_idx];
-            if (arc.type == MMU_Graph::ARC_TYPE::NON_BORDER || used_arcs[arc_idx])
-                continue;
-
-            Linef process_line(graph.nodes[arc.from_idx].point, graph.nodes[arc.to_idx].point);
-            used_arcs[arc_idx] = true;
-
-            std::vector<std::pair<size_t, Linef>> arc_id_to_face_lines;
-            arc_id_to_face_lines.emplace_back(std::make_pair(arc_idx, process_line));
-            Vec2d start_p = process_line.a;
-
-            Linef                 p_vec = process_line;
-            const MMU_Graph::Arc *p_arc = &arc;
-            bool                  flag  = false;
-            do {
-                std::vector<const MMU_Graph::Arc *> nexts = get_next_arc(graph, used_arcs, p_vec, *p_arc, arc.color);
-                for (auto next : nexts) {
-                    size_t next_arc_idx = next - &graph.arcs.front();
-                    if (used_arcs[next_arc_idx]) {
-                        flag = true;
-                        break;
-                    }
-                }
-
-                if (flag)
-                    break;
-
-                for (auto next : nexts) {
-                    size_t next_arc_idx = next - &graph.arcs.front();
-                    arc_id_to_face_lines.emplace_back(std::make_pair(next_arc_idx, Linef(graph.nodes[next->from_idx].point, graph.nodes[next->to_idx].point)));
-                    used_arcs[next_arc_idx] = true;
-                }
-
-                p_vec = Linef(graph.nodes[nexts.back()->from_idx].point, graph.nodes[nexts.back()->to_idx].point);
-                p_arc = nexts.back();
-
-            } while (graph.nodes[p_arc->to_idx].point != start_p || !all_arc_used(graph.nodes[p_arc->to_idx]));
-
-            if (Polygon poly = to_polygon(arc_id_to_face_lines); poly.is_counter_clockwise() && poly.is_valid()) {
-                expolygons_segments[arc.color].emplace_back(std::move(poly));
-            } else{
-                while (arc_id_to_face_lines.size() > 1)
-                {
-                    auto id_to_line = arc_id_to_face_lines.back();
-                    used_arcs[id_to_line.first] = false;
-                    arc_id_to_face_lines.pop_back();
-                    Linef add_line(arc_id_to_face_lines.back().second.b, arc_id_to_face_lines.front().second.a);
-                    arc_id_to_face_lines.emplace_back(std::make_pair(-1, add_line));
-                    Polygon poly = to_polygon(arc_id_to_face_lines);
-                    if (!is_profile_self_interaction(poly) && poly.is_counter_clockwise() && poly.is_valid()) {
-                        expolygons_segments[arc.color].emplace_back(std::move(poly));
-                        break;
-                    }
-                    arc_id_to_face_lines.pop_back();
-                }
-            }
-        }
-    }
-    return expolygons_segments;
+    return mk_vector_vec2d(&edge).norm();
 }
 
 // Used in remove_multiple_edges_in_vertices()
 // Returns length of edge with is connected to contour. To this length is include other edges with follows it if they are almost straight (with the
 // tolerance of 15) And also if node between two subsequent edges is connected only to these two edges.
-static inline double compute_edge_length(const MMU_Graph &graph, const size_t start_idx, const size_t &start_arc_idx)
+static inline double calc_total_edge_length(const VD::edge_type &starting_edge)
 {
-    assert(start_arc_idx < graph.arcs.size());
-    std::vector<bool> used_arcs(graph.arcs.size(), false);
-
-    used_arcs[start_arc_idx]                = true;
-    const MMU_Graph::Arc *arc               = &graph.arcs[start_arc_idx];
-    size_t                idx               = start_idx;
-    double                line_total_length = (graph.nodes[arc->to_idx].point - graph.nodes[idx].point).norm();
-    while (graph.nodes[arc->to_idx].arc_idxs.size() == 2) {
-        bool found = false;
-        for (const size_t &arc_idx : graph.nodes[arc->to_idx].arc_idxs) {
-            if (const MMU_Graph::Arc &arc_n = graph.arcs[arc_idx]; arc_n.type == MMU_Graph::ARC_TYPE::NON_BORDER && !used_arcs[arc_idx] && arc_n.to_idx != idx) {
-                Linef first_line(graph.nodes[idx].point, graph.nodes[arc->to_idx].point);
-                Linef second_line(graph.nodes[arc->to_idx].point, graph.nodes[arc_n.to_idx].point);
-
-                Vec2d  first_line_vec    = (first_line.a - first_line.b);
-                Vec2d  second_line_vec   = (second_line.b - second_line.a);
-                Vec2d  first_line_vec_n  = first_line_vec.normalized();
-                Vec2d  second_line_vec_n = second_line_vec.normalized();
-                double angle             = ::acos(std::clamp(first_line_vec_n.dot(second_line_vec_n), -1.0, 1.0));
-                if (Slic3r::cross2(first_line_vec_n, second_line_vec_n) < 0.0)
-                    angle = 2.0 * (double) PI - angle;
-
-                if (std::abs(angle - PI) >= (PI / 12))
-                    continue;
-
-                idx = arc->to_idx;
-                arc = &arc_n;
-
-                line_total_length += (graph.nodes[arc->to_idx].point - graph.nodes[idx].point).norm();
-                used_arcs[arc_idx] = true;
-                found      = true;
-                break;
-            }
-        }
-        if (!found)
+    double               total_edge_length = edge_length(starting_edge);
+    const VD::edge_type *prev              = &starting_edge;
+    do {
+        if (prev->is_finite() && non_deleted_edge_count(*prev->vertex1()) > 2)
             break;
-    }
 
-    return line_total_length;
+        bool                 found_next_edge = false;
+        const VD::edge_type *current         = prev->next();
+        do {
+            if (current->color() == VD_ANNOTATION::DELETED)
+                continue;
+
+            Vec2d  first_line_vec_n  = mk_flipped_vector_vec2d(prev).normalized();
+            Vec2d  second_line_vec_n = mk_vector_vec2d(current).normalized();
+            double angle             = ::acos(std::clamp(first_line_vec_n.dot(second_line_vec_n), -1.0, 1.0));
+            if (Slic3r::cross2(first_line_vec_n, second_line_vec_n) < 0.0)
+                angle = 2.0 * (double) PI - angle;
+
+            if (std::abs(angle - PI) >= (PI / 12))
+                continue;
+
+            prev               = current;
+            found_next_edge    = true;
+            total_edge_length += edge_length(*current);
+
+            break;
+        } while (current = current->prev()->twin(), current != prev->next());
+
+        if (!found_next_edge)
+            break;
+
+    } while (prev != &starting_edge);
+
+    return total_edge_length;
 }
 
-// Used for fixing double Voronoi edges for concave parts of the polygon.
-static void remove_multiple_edges_in_vertices(MMU_Graph &graph, const std::vector<std::vector<ColoredLine>> &color_poly)
+// When a Voronoi vertex has more than one Voronoi edge (for example, in concave parts of a polygon),
+// we leave just one Voronoi edge in the Voronoi vertex.
+// This Voronoi edge is selected based on a heuristic.
+static void remove_multiple_edges_in_vertex(const VD::vertex_type &vertex) {
+    if (non_deleted_edge_count(vertex) <= 1)
+        return;
+
+    std::vector<std::pair<const VD::edge_type *, double>> edges_to_check;
+    const VD::edge_type *edge = vertex.incident_edge();
+    do {
+        if (edge->color() == VD_ANNOTATION::DELETED)
+            continue;
+
+        edges_to_check.emplace_back(edge, calc_total_edge_length(*edge));
+    } while (edge = edge->prev()->twin(), edge != vertex.incident_edge());
+
+    std::sort(edges_to_check.begin(), edges_to_check.end(), [](const auto &l, const auto &r) -> bool {
+        return l.second > r.second;
+    });
+
+    while (edges_to_check.size() > 1) {
+        const VD::edge_type &edge_to_check = *edges_to_check.back().first;
+        edge_to_check.color(VD_ANNOTATION::DELETED);
+        edge_to_check.twin()->color(VD_ANNOTATION::DELETED);
+
+        if (const VD::vertex_type &vertex_to_delete = *edge_to_check.vertex1(); can_vertex_be_deleted(vertex_to_delete))
+            delete_vertex_deep(vertex_to_delete);
+
+        edges_to_check.pop_back();
+    }
+}
+
+// Returns list of ExPolygons for each extruder + 1 for default unpainted regions.
+// It iterates through all nodes on the border between two different colors, and from this point,
+// start selection always left most edges for every node to construct CCW polygons.
+static std::vector<ExPolygons> extract_colored_segments(const std::vector<ColoredLines> &colored_polygons,
+                                                        const size_t                     num_extruders,
+                                                        const size_t                     layer_idx)
 {
-    std::vector<std::vector<std::pair<size_t, size_t>>> colored_segments = get_all_segments(color_poly);
-    for (const std::vector<std::pair<size_t, size_t>> &colored_segment_p : colored_segments) {
-        size_t poly_idx = &colored_segment_p - &colored_segments.front();
-        for (const std::pair<size_t, size_t> &colored_segment : colored_segment_p) {
-            size_t first_idx  = graph.get_global_index(poly_idx, colored_segment.first);
-            size_t second_idx = graph.get_global_index(poly_idx, (colored_segment.second + 1) % graph.polygon_sizes[poly_idx]);
-            Linef  seg_line(graph.nodes[first_idx].point, graph.nodes[second_idx].point);
+    const ColoredLines colored_lines = to_lines(colored_polygons);
+    const BoundingBox  bbox          = get_extents(colored_polygons);
 
-            if (graph.nodes[first_idx].arc_idxs.size() >= 3) {
-                std::vector<std::pair<MMU_Graph::Arc *, double>> arc_to_check;
-                for (const size_t &arc_idx : graph.nodes[first_idx].arc_idxs) {
-                    MMU_Graph::Arc &n_arc = graph.arcs[arc_idx];
-                    if (n_arc.type == MMU_Graph::ARC_TYPE::NON_BORDER) {
-                        double total_len = compute_edge_length(graph, first_idx, arc_idx);
-                        arc_to_check.emplace_back(&n_arc, total_len);
-                    }
-                }
-                std::sort(arc_to_check.begin(), arc_to_check.end(),
-                          [](std::pair<MMU_Graph::Arc *, double> &l, std::pair<MMU_Graph::Arc *, double> &r) -> bool { return l.second > r.second; });
+    auto get_next_contour_line = [&colored_polygons](const ColoredLine &line) -> const ColoredLine & {
+        size_t contour_line_size = colored_polygons[line.poly_idx].size();
+        size_t contour_next_idx  = (line.local_line_idx + 1) % contour_line_size;
+        return colored_polygons[line.poly_idx][contour_next_idx];
+    };
 
-                while (arc_to_check.size() > 1) {
-                    graph.remove_edge(first_idx, arc_to_check.back().first->to_idx);
-                    arc_to_check.pop_back();
-                }
-            }
+    Voronoi::VD vd;
+    vd.construct_voronoi(colored_lines.begin(), colored_lines.end());
+
+    // First, mark each Voronoi vertex on the input polygon to prevent it from being deleted later.
+    for (const Voronoi::VD::cell_type &cell : vd.cells()) {
+        if (cell.is_degenerate() || !cell.contains_segment())
+            continue;
+
+        if (const Geometry::SegmentCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_segment_cell_range(cell, colored_lines.begin(), colored_lines.end()); cell_range.is_valid())
+            cell_range.edge_begin->vertex0()->color(VD_ANNOTATION::VERTEX_ON_CONTOUR);
+    }
+
+    // Second, remove all Voronoi vertices that are outside the bounding box of input polygons.
+    // Such Voronoi vertices are definitely not inside of input polygons, so we don't care about them.
+    for (const Voronoi::VD::vertex_type &vertex : vd.vertices()) {
+        if (vertex.color() == VD_ANNOTATION::DELETED || vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR)
+            continue;
+
+        if (!Geometry::VoronoiUtils::is_in_range<coord_t>(vertex) || !bbox.contains(Geometry::VoronoiUtils::to_point(vertex).cast<coord_t>()))
+            delete_vertex_deep(vertex);
+    }
+
+    // Third, remove all Voronoi edges that are infinite.
+    for (const Voronoi::VD::edge_type &edge : vd.edges()) {
+        if (edge.color() != VD_ANNOTATION::DELETED && edge.is_infinite()) {
+            edge.color(VD_ANNOTATION::DELETED);
+            edge.twin()->color(VD_ANNOTATION::DELETED);
+
+            if (edge.vertex0() != nullptr && can_vertex_be_deleted(*edge.vertex0()))
+                delete_vertex_deep(*edge.vertex0());
+
+            if (edge.vertex1() != nullptr && can_vertex_be_deleted(*edge.vertex1()))
+                delete_vertex_deep(*edge.vertex1());
         }
     }
+
+    // Fourth, remove all edges that point outward from the input polygon.
+    for (Voronoi::VD::cell_type cell : vd.cells()) {
+        if (cell.is_degenerate() || !cell.contains_segment())
+            continue;
+
+        if (const Geometry::SegmentCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_segment_cell_range(cell, colored_lines.begin(), colored_lines.end()); cell_range.is_valid()) {
+            const ColoredLine &current_line = Geometry::VoronoiUtils::get_source_segment(cell, colored_lines.begin(), colored_lines.end());
+            const ColoredLine &next_line = get_next_contour_line(current_line);
+
+            const VD::edge_type *edge = cell_range.edge_begin;
+            do {
+                if (edge->color() == VD_ANNOTATION::DELETED)
+                    continue;
+
+                if (!points_inside(current_line.line, next_line.line, Geometry::VoronoiUtils::to_point(edge->vertex1()).cast<coord_t>())) {
+                    edge->color(VD_ANNOTATION::DELETED);
+                    edge->twin()->color(VD_ANNOTATION::DELETED);
+                    delete_vertex_deep(*edge->vertex1());
+                }
+            } while (edge = edge->prev()->twin(), edge != cell_range.edge_begin);
+        }
+    }
+
+    // Fifth, if a Voronoi vertex has more than one Voronoi edge, remove all but one of them based on heuristics.
+    for (const Voronoi::VD::vertex_type &vertex : vd.vertices()) {
+        if (vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR)
+            remove_multiple_edges_in_vertex(vertex);
+    }
+
+#ifdef MM_SEGMENTATION_DEBUG_GRAPH
+    {
+        static int iRun = 0;
+        export_graph_to_svg(debug_out_path("mm-graph-%d-%d.svg", layer_idx, iRun++), vd, colored_polygons);
+    }
+#endif // MM_SEGMENTATION_DEBUG_GRAPH
+
+    // Sixth, extract the colored segments from the annotated Voronoi diagram.
+    std::vector<ExPolygons> segmented_expolygons_per_extruder(num_extruders + 1);
+    for (const Voronoi::VD::cell_type &cell : vd.cells()) {
+        if (cell.is_degenerate() || !cell.contains_segment())
+            continue;
+
+        if (const Geometry::SegmentCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_segment_cell_range(cell, colored_lines.begin(), colored_lines.end()); cell_range.is_valid()) {
+            if (cell_range.edge_begin->vertex0()->color() != VD_ANNOTATION::VERTEX_ON_CONTOUR)
+                continue;
+
+            const ColoredLine source_segment = Geometry::VoronoiUtils::get_source_segment(cell, colored_lines.begin(), colored_lines.end());
+
+            Polygon segmented_polygon;
+            segmented_polygon.points.emplace_back(source_segment.line.b);
+
+            // We have ensured that each segmented_polygon have to start at edge_begin->vertex0() and end at edge_end->vertex1().
+            const VD::edge_type *edge = cell_range.edge_begin;
+            do {
+                if (edge->color() == VD_ANNOTATION::DELETED)
+                    continue;
+
+                const VD::vertex_type &next_vertex = *edge->vertex1();
+                segmented_polygon.points.emplace_back(Geometry::VoronoiUtils::to_point(next_vertex).cast<coord_t>());
+                edge->color(VD_ANNOTATION::DELETED);
+
+                if (next_vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR || next_vertex.color() == VD_ANNOTATION::DELETED) {
+                    assert(next_vertex.color() == VD_ANNOTATION::VERTEX_ON_CONTOUR);
+                    break;
+                }
+
+                edge = edge->twin();
+            } while (edge = edge->twin()->next(), edge != cell_range.edge_begin);
+
+            if (edge->vertex1() != cell_range.edge_end->vertex1())
+                continue;
+
+            cell_range.edge_begin->vertex0()->color(VD_ANNOTATION::DELETED);
+            segmented_expolygons_per_extruder[source_segment.color].emplace_back(std::move(segmented_polygon));
+        }
+    }
+
+    // Merge all polygons together for each extruder
+    for (auto &segmented_expolygons : segmented_expolygons_per_extruder)
+        segmented_expolygons = union_ex(segmented_expolygons);
+
+    return segmented_expolygons_per_extruder;
 }
 
 static void cut_segmented_layers(const std::vector<ExPolygons>        &input_expolygons,
@@ -1375,12 +835,12 @@ static void cut_segmented_layers(const std::vector<ExPolygons>        &input_exp
                                  const float                           interlocking_depth,
                                  const std::function<void()>          &throw_on_cancel_callback)
 {
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - cutting segmented layers in parallel - begin";
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, segmented_regions.size()),
-                      [&segmented_regions, &input_expolygons, &cut_width, &interlocking_depth, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - cutting segmented layers in parallel - begin";
+    const float interlocking_cut_width = interlocking_depth > 0.f ? std::max(cut_width - interlocking_depth, 0.f) : 0.f;
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, segmented_regions.size()),[&segmented_regions, &input_expolygons, &cut_width, &interlocking_cut_width, &throw_on_cancel_callback](const tbb::blocked_range<size_t>& range) {
         for (size_t layer_idx = range.begin(); layer_idx < range.end(); ++layer_idx) {
             throw_on_cancel_callback();
-            const float  region_cut_width = ((layer_idx % 2 == 0) && (interlocking_depth != 0.f)) ? interlocking_depth : cut_width;
+            const float  region_cut_width       = (layer_idx % 2 == 0 && interlocking_cut_width > 0.f) ? interlocking_cut_width : cut_width;
             const size_t num_extruders_plus_one = segmented_regions[layer_idx].size();
             if (region_cut_width > 0.f) {
                 std::vector<ExPolygons> segmented_regions_cuts(num_extruders_plus_one); // Indexed by extruder_id
@@ -1391,7 +851,7 @@ static void cut_segmented_layers(const std::vector<ExPolygons>        &input_exp
             }
         }
     }); // end of parallel_for
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - cutting segmented layers in parallel - end";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - cutting segmented layers in parallel - end";
 }
 
 static bool is_volume_sinking(const indexed_triangle_set &its, const Transform3d &trafo)
@@ -1431,9 +891,9 @@ static inline std::vector<std::vector<ExPolygons>> mmu_segmentation_top_and_bott
     std::vector<float> zs = zs_from_layers(layers);
     Transform3d        object_trafo = print_object.trafo_centered();
 
-#ifdef MMU_SEGMENTATION_DEBUG_TOP_BOTTOM
+#ifdef MM_SEGMENTATION_DEBUG_TOP_BOTTOM
     static int iRun = 0;
-#endif // NDEBUG
+#endif // MM_SEGMENTATION_DEBUG_TOP_BOTTOM
 
     if (max_top_layers > 0 || max_bottom_layers > 0) {
         for (const ModelVolume *mv : print_object.model_object()->volumes)
@@ -1441,12 +901,12 @@ static inline std::vector<std::vector<ExPolygons>> mmu_segmentation_top_and_bott
                 const Transform3d volume_trafo = object_trafo * mv->get_matrix();
                 for (size_t extruder_idx = 0; extruder_idx < num_extruders; ++ extruder_idx) {
                     const indexed_triangle_set painted = mv->mmu_segmentation_facets.get_facets_strict(*mv, EnforcerBlockerType(extruder_idx));
-#ifdef MMU_SEGMENTATION_DEBUG_TOP_BOTTOM
+#ifdef MM_SEGMENTATION_DEBUG_TOP_BOTTOM
                     {
                         static int iRun = 0;
                         its_write_obj(painted, debug_out_path("mm-painted-patch-%d-%d.obj", iRun ++, extruder_idx).c_str());
                     }
-#endif // MMU_SEGMENTATION_DEBUG_TOP_BOTTOM
+#endif // MM_SEGMENTATION_DEBUG_TOP_BOTTOM
                     if (! painted.indices.empty()) {
                         std::vector<Polygons> top, bottom;
                         if (!zs.empty() && is_volume_sinking(painted, volume_trafo)) {
@@ -1501,7 +961,7 @@ static inline std::vector<std::vector<ExPolygons>> mmu_segmentation_top_and_bott
     filter_out_small_polygons(top_raw, Slic3r::sqr(scale_(0.1f)));
     filter_out_small_polygons(bottom_raw, Slic3r::sqr(scale_(0.1f)));
 
-#ifdef MMU_SEGMENTATION_DEBUG_TOP_BOTTOM
+#ifdef MM_SEGMENTATION_DEBUG_TOP_BOTTOM
     {
         const char* colors[] = { "aqua", "black", "blue", "fuchsia", "gray", "green", "lime", "maroon", "navy", "olive", "purple", "red", "silver", "teal", "yellow" };
         static int iRun = 0;
@@ -1523,7 +983,7 @@ static inline std::vector<std::vector<ExPolygons>> mmu_segmentation_top_and_bott
         }
         ++ iRun;
     }
-#endif // MMU_SEGMENTATION_DEBUG_TOP_BOTTOM
+#endif // MM_SEGMENTATION_DEBUG_TOP_BOTTOM
 
     // When the upper surface of an object is occluded, it should no longer be considered the upper surface
     {
@@ -1708,7 +1168,7 @@ static std::vector<std::vector<ExPolygons>> merge_segmented_layers(
     segmented_regions_merged.assign(num_layers, std::vector<ExPolygons>(num_extruders));
     assert(num_extruders + 1 == top_and_bottom_layers.size());
 
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - merging segmented layers in parallel - begin";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - merging segmented layers in parallel - begin";
     tbb::parallel_for(tbb::blocked_range<size_t>(0, num_layers), [&segmented_regions, &top_and_bottom_layers, &segmented_regions_merged, &num_extruders, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
         for (size_t layer_idx = range.begin(); layer_idx < range.end(); ++layer_idx) {
             assert(segmented_regions[layer_idx].size() == num_extruders + 1);
@@ -1735,12 +1195,12 @@ static std::vector<std::vector<ExPolygons>> merge_segmented_layers(
             }
         }
     }); // end of parallel_for
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - merging segmented layers in parallel - end";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - merging segmented layers in parallel - end";
 
     return segmented_regions_merged;
 }
 
-#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
+#ifdef MM_SEGMENTATION_DEBUG_REGIONS
 static void export_regions_to_svg(const std::string &path, const std::vector<ExPolygons> &regions, const ExPolygons &lslices)
 {
     const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "yellow"};
@@ -1752,35 +1212,15 @@ static void export_regions_to_svg(const std::string &path, const std::vector<ExP
     svg.draw_outline(lslices, "green", "lime", stroke_width);
     for (const ExPolygons &by_extruder : regions) {
         size_t extrude_idx = &by_extruder - &regions.front();
-        if (extrude_idx >= 0 && extrude_idx < int(colors.size()))
-            svg.draw(by_extruder, colors[extrude_idx], stroke_width);
+        if (extrude_idx < int(colors.size()))
+            svg.draw(by_extruder, colors[extrude_idx]);
         else
-            svg.draw(by_extruder, "black", stroke_width);
+            svg.draw(by_extruder, "black");
     }
 }
-#endif // MMU_SEGMENTATION_DEBUG_REGIONS
+#endif // MM_SEGMENTATION_DEBUG_REGIONS
 
-#ifdef MMU_SEGMENTATION_DEBUG_GRAPH
-static void export_graph_to_svg(const std::string &path, const MMU_Graph &graph, const ExPolygons &lslices)
-{
-    const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "green", "yellow"};
-    coordf_t                       stroke_width = scale_(0.05);
-    BoundingBox                    bbox         = get_extents(lslices);
-    bbox.offset(scale_(1.));
-    ::Slic3r::SVG svg(path.c_str(), bbox);
-    for (const MMU_Graph::Node &node : graph.nodes)
-        for (const size_t &arc_idx : node.arc_idxs) {
-            const MMU_Graph::Arc &arc = graph.arcs[arc_idx];
-            Line arc_line(mk_point(node.point), mk_point(graph.nodes[arc.to_idx].point));
-            if (arc.type == MMU_Graph::ARC_TYPE::BORDER && arc.color >= 0 && arc.color < int(colors.size()))
-                svg.draw(arc_line, colors[arc.color], stroke_width);
-            else
-                svg.draw(arc_line, "black", stroke_width);
-        }
-}
-#endif // MMU_SEGMENTATION_DEBUG_GRAPH
-
-#ifdef MMU_SEGMENTATION_DEBUG_INPUT
+#ifdef MM_SEGMENTATION_DEBUG_INPUT
 void export_processed_input_expolygons_to_svg(const std::string &path, const LayerRegionPtrs &regions, const ExPolygons &processed_input_expolygons)
 {
     coordf_t    stroke_width = scale_(0.05);
@@ -1790,13 +1230,14 @@ void export_processed_input_expolygons_to_svg(const std::string &path, const Lay
     ::Slic3r::SVG svg(path.c_str(), bbox);
 
     for (LayerRegion *region : regions)
-        svg.draw_outline(region->slices.surfaces, "blue", "cyan", stroke_width);
+        for (const Surface &surface : region->slices.surfaces)
+            svg.draw_outline(surface, "blue", "cyan", stroke_width);
 
     svg.draw_outline(processed_input_expolygons, "red", "pink", stroke_width);
 }
-#endif // MMU_SEGMENTATION_DEBUG_INPUT
+#endif // MM_SEGMENTATION_DEBUG_INPUT
 
-#ifdef MMU_SEGMENTATION_DEBUG_PAINTED_LINES
+#ifdef MM_SEGMENTATION_DEBUG_PAINTED_LINES
 static void export_painted_lines_to_svg(const std::string &path, const std::vector<std::vector<PaintedLine>> &all_painted_lines, const ExPolygons &lslices)
 {
     const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "yellow"};
@@ -1812,10 +1253,10 @@ static void export_painted_lines_to_svg(const std::string &path, const std::vect
         for (const PaintedLine &painted_line : painted_lines)
             svg.draw(painted_line.projected_line, painted_line.color < int(colors.size()) ? colors[painted_line.color] : "black", stroke_width);
 }
-#endif // MMU_SEGMENTATION_DEBUG_PAINTED_LINES
+#endif // MM_SEGMENTATION_DEBUG_PAINTED_LINES
 
-#ifdef MMU_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
-static void export_colorized_polygons_to_svg(const std::string &path, const std::vector<std::vector<ColoredLine>> &colorized_polygons, const ExPolygons &lslices)
+#ifdef MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
+static void export_colorized_polygons_to_svg(const std::string &path, const std::vector<ColoredLines> &colorized_polygons, const ExPolygons &lslices)
 {
     const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "green", "yellow"};
     coordf_t                       stroke_width = scale_(0.05);
@@ -1823,19 +1264,19 @@ static void export_colorized_polygons_to_svg(const std::string &path, const std:
     bbox.offset(scale_(1.));
     ::Slic3r::SVG svg(path.c_str(), bbox);
 
-    for (const std::vector<ColoredLine> &colorized_polygon : colorized_polygons)
+    for (const ColoredLines &colorized_polygon : colorized_polygons)
         for (const ColoredLine &colorized_line : colorized_polygon)
             svg.draw(colorized_line.line, colorized_line.color < int(colors.size())? colors[colorized_line.color] : "black", stroke_width);
 }
-#endif // MMU_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
+#endif // MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
 
 // Check if all ColoredLine representing a single layer uses the same color.
-static bool has_layer_only_one_color(const std::vector<std::vector<ColoredLine>> &colored_polygons)
+static bool has_layer_only_one_color(const std::vector<ColoredLines> &colored_polygons)
 {
     assert(!colored_polygons.empty());
     assert(!colored_polygons.front().empty());
     int first_line_color = colored_polygons.front().front().color;
-    for (const std::vector<ColoredLine> &colored_polygon : colored_polygons)
+    for (const ColoredLines &colored_polygon : colored_polygons)
         for (const ColoredLine &colored_line : colored_polygon)
             if (first_line_color != colored_line.color)
                 return false;
@@ -1857,8 +1298,12 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
 
     throw_on_cancel_callback();
 
+#ifdef MM_SEGMENTATION_DEBUG
+    static int iRun = 0;
+#endif // MM_SEGMENTATION_DEBUG
+
     // Merge all regions and remove small holes
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - slices preparation in parallel - begin";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - slices preparation in parallel - begin";
     tbb::parallel_for(tbb::blocked_range<size_t>(0, num_layers), [&layers, &input_expolygons, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
         for (size_t layer_idx = range.begin(); layer_idx < range.end(); ++layer_idx) {
             throw_on_cancel_callback();
@@ -1878,17 +1323,14 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
             // Such close points sometimes caused that the Voronoi diagram has self-intersecting edges around these vertices.
             // This consequently leads to issues with the extraction of colored segments by function extract_colored_segments.
             // Calling expolygons_simplify fixed these issues.
-            input_expolygons[layer_idx] = remove_duplicates(expolygons_simplify(offset_ex(ex_polygons, -10.f * float(SCALED_EPSILON)), 5 * SCALED_EPSILON), scaled<coord_t>(0.01), PI / 6);
+            input_expolygons[layer_idx] = remove_duplicates(expolygons_simplify(offset_ex(ex_polygons, -10.f * float(SCALED_EPSILON)), 5 * SCALED_EPSILON), scaled<coord_t>(0.01), PI/6);
 
-#ifdef MMU_SEGMENTATION_DEBUG_INPUT
-            {
-                static int iRun = 0;
-                export_processed_input_expolygons_to_svg(debug_out_path("mm-input-%d-%d.svg", layer_idx, iRun++), layers[layer_idx]->regions(), input_expolygons[layer_idx]);
-            }
-#endif // MMU_SEGMENTATION_DEBUG_INPUT
+#ifdef MM_SEGMENTATION_DEBUG_INPUT
+            export_processed_input_expolygons_to_svg(debug_out_path("mm-input-%d-%d.svg", layer_idx, iRun), layers[layer_idx]->regions(), input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_INPUT
         }
     }); // end of parallel_for
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - slices preparation in parallel - end";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - slices preparation in parallel - end";
 
     std::vector<BoundingBox> layer_bboxes(num_layers);
     for (size_t layer_idx = 0; layer_idx < num_layers; ++layer_idx) {
@@ -1911,7 +1353,7 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
         edge_grids[layer_idx].create(input_expolygons[layer_idx], coord_t(scale_(10.)));
     }
 
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - projection of painted triangles - begin";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - projection of painted triangles - begin";
     for (const ModelVolume *mv : print_object.model_object()->volumes) {
         tbb::parallel_for(tbb::blocked_range<size_t>(1, num_extruders + 1), [&mv, &print_object, &layers, &edge_grids, &painted_lines, &painted_lines_mutex, &input_expolygons, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
             for (size_t extruder_idx = range.begin(); extruder_idx < range.end(); ++extruder_idx) {
@@ -2002,39 +1444,30 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
             }
         }); // end of parallel_for
     }
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - projection of painted triangles - end";
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - painted layers count: "
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - projection of painted triangles - end";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - painted layers count: "
                              << std::count_if(painted_lines.begin(), painted_lines.end(), [](const std::vector<PaintedLine> &pl) { return !pl.empty(); });
 
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - layers segmentation in parallel - begin";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - layers segmentation in parallel - begin";
     tbb::parallel_for(tbb::blocked_range<size_t>(0, num_layers), [&edge_grids, &input_expolygons, &painted_lines, &segmented_regions, &num_extruders, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
         for (size_t layer_idx = range.begin(); layer_idx < range.end(); ++layer_idx) {
             throw_on_cancel_callback();
             if (!painted_lines[layer_idx].empty()) {
-#ifdef MMU_SEGMENTATION_DEBUG_PAINTED_LINES
-                {
-                    static int iRun = 0;
-                    export_painted_lines_to_svg(debug_out_path("mm-painted-lines-%d-%d.svg", layer_idx, iRun++), {painted_lines[layer_idx]}, input_expolygons[layer_idx]);
-                }
-#endif // MMU_SEGMENTATION_DEBUG_PAINTED_LINES
+#ifdef MM_SEGMENTATION_DEBUG_PAINTED_LINES
+                export_painted_lines_to_svg(debug_out_path("0-mm-painted-lines-%d-%d.svg", layer_idx, iRun), {painted_lines[layer_idx]}, input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_PAINTED_LINES
 
                 std::vector<std::vector<PaintedLine>> post_processed_painted_lines = post_process_painted_lines(edge_grids[layer_idx].contours(), std::move(painted_lines[layer_idx]));
 
-#ifdef MMU_SEGMENTATION_DEBUG_PAINTED_LINES
-                {
-                    static int iRun = 0;
-                    export_painted_lines_to_svg(debug_out_path("mm-painted-lines-post-processed-%d-%d.svg", layer_idx, iRun++), post_processed_painted_lines, input_expolygons[layer_idx]);
-                }
-#endif // MMU_SEGMENTATION_DEBUG_PAINTED_LINES
+#ifdef MM_SEGMENTATION_DEBUG_PAINTED_LINES
+                export_painted_lines_to_svg(debug_out_path("1-mm-painted-lines-post-processed-%d-%d.svg", layer_idx, iRun), post_processed_painted_lines, input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_PAINTED_LINES
 
-                std::vector<std::vector<ColoredLine>> color_poly = colorize_contours(edge_grids[layer_idx].contours(), post_processed_painted_lines);
+                std::vector<ColoredLines> color_poly = colorize_contours(edge_grids[layer_idx].contours(), post_processed_painted_lines);
 
-#ifdef MMU_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
-                {
-                    static int iRun = 0;
-                    export_colorized_polygons_to_svg(debug_out_path("mm-colorized_polygons-%d-%d.svg", layer_idx, iRun++), color_poly, input_expolygons[layer_idx]);
-                }
-#endif // MMU_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
+#ifdef MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
+                export_colorized_polygons_to_svg(debug_out_path("2-mm-colorized_polygons-%d-%d.svg", layer_idx, iRun), color_poly, input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_COLORIZED_POLYGONS
 
                 assert(!color_poly.empty());
                 assert(!color_poly.front().empty());
@@ -2042,33 +1475,19 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
                     // If the whole layer is painted using the same color, it is not needed to construct a Voronoi diagram for the segmentation of this layer.
                     segmented_regions[layer_idx][size_t(color_poly.front().front().color)] = input_expolygons[layer_idx];
                 } else {
-                    MMU_Graph graph = build_graph(layer_idx, color_poly);
-                    remove_multiple_edges_in_vertices(graph, color_poly);
-                    graph.remove_nodes_with_one_arc();
-
-#ifdef MMU_SEGMENTATION_DEBUG_GRAPH
-                    {
-                        static int iRun = 0;
-                        export_graph_to_svg(debug_out_path("mm-graph-final-%d-%d.svg", layer_idx, iRun++), graph, input_expolygons[layer_idx]);
-                    }
-#endif // MMU_SEGMENTATION_DEBUG_GRAPH
-
-                    segmented_regions[layer_idx] = extract_colored_segments(graph, num_extruders);
+                    segmented_regions[layer_idx] = extract_colored_segments(color_poly, num_extruders, layer_idx);
                 }
 
-#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
-                {
-                    static int iRun = 0;
-                    export_regions_to_svg(debug_out_path("mm-regions-sides-%d-%d.svg", layer_idx, iRun++), segmented_regions[layer_idx], input_expolygons[layer_idx]);
-                }
-#endif // MMU_SEGMENTATION_DEBUG_REGIONS
+#ifdef MM_SEGMENTATION_DEBUG_REGIONS
+                export_regions_to_svg(debug_out_path("3-mm-regions-sides-%d-%d.svg", layer_idx, iRun), segmented_regions[layer_idx], input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_REGIONS
             }
         }
     }); // end of parallel_for
-    BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - layers segmentation in parallel - end";
+    BOOST_LOG_TRIVIAL(debug) << "MM segmentation - layers segmentation in parallel - end";
     throw_on_cancel_callback();
 
-    if (auto max_width = print_object.config().mmu_segmented_region_max_width, interlocking_depth = print_object.config().mmu_segmented_region_interlocking_depth; max_width > 0.f || interlocking_depth > 0.f) {
+    if (auto max_width = print_object.config().mmu_segmented_region_max_width, interlocking_depth = print_object.config().mmu_segmented_region_interlocking_depth; max_width > 0.f) {
         cut_segmented_layers(input_expolygons, segmented_regions, float(scale_(max_width)), float(scale_(interlocking_depth)), throw_on_cancel_callback);
         throw_on_cancel_callback();
     }
@@ -2080,13 +1499,14 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
     std::vector<std::vector<ExPolygons>> segmented_regions_merged = merge_segmented_layers(segmented_regions, std::move(top_and_bottom_layers), num_extruders, throw_on_cancel_callback);
     throw_on_cancel_callback();
 
-#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
-    {
-        static int iRun = 0;
-        for (size_t layer_idx = 0; layer_idx < print_object.layers().size(); ++layer_idx)
-            export_regions_to_svg(debug_out_path("mm-regions-merged-%d-%d.svg", layer_idx, iRun++), segmented_regions_merged[layer_idx], input_expolygons[layer_idx]);
-    }
-#endif // MMU_SEGMENTATION_DEBUG_REGIONS
+#ifdef MM_SEGMENTATION_DEBUG_REGIONS
+    for (size_t layer_idx = 0; layer_idx < print_object.layers().size(); ++layer_idx)
+        export_regions_to_svg(debug_out_path("4-mm-regions-merged-%d-%d.svg", layer_idx, iRun), segmented_regions_merged[layer_idx], input_expolygons[layer_idx]);
+#endif // MM_SEGMENTATION_DEBUG_REGIONS
+
+#ifdef MM_SEGMENTATION_DEBUG
+    ++iRun;
+#endif // MM_SEGMENTATION_DEBUG
 
     return segmented_regions_merged;
 }
diff --git a/src/libslic3r/MultiMaterialSegmentation.hpp b/src/libslic3r/MultiMaterialSegmentation.hpp
index 4efdc6951..91d0f298b 100644
--- a/src/libslic3r/MultiMaterialSegmentation.hpp
+++ b/src/libslic3r/MultiMaterialSegmentation.hpp
@@ -6,13 +6,41 @@
 
 namespace Slic3r {
 
-
 class PrintObject;
 class ExPolygon;
+using ExPolygons = std::vector<ExPolygon>;
+
+struct ColoredLine
+{
+    Line line;
+    int  color;
+    int  poly_idx       = -1;
+    int  local_line_idx = -1;
+};
+
+using ColoredLines = std::vector<ColoredLine>;
 
 // Returns MMU segmentation based on painting in MMU segmentation gizmo
 std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(const PrintObject &print_object, const std::function<void()> &throw_on_cancel_callback);
 
 } // namespace Slic3r
 
+namespace boost::polygon {
+template<> struct geometry_concept<Slic3r::ColoredLine>
+{
+    typedef segment_concept type;
+};
+
+template<> struct segment_traits<Slic3r::ColoredLine>
+{
+    typedef coord_t       coordinate_type;
+    typedef Slic3r::Point point_type;
+
+    static inline point_type get(const Slic3r::ColoredLine &line, const direction_1d &dir)
+    {
+        return dir.to_int() ? line.line.b : line.line.a;
+    }
+};
+} // namespace boost::polygon
+
 #endif // slic3r_MultiMaterialSegmentation_hpp_