Welcome to mirror list, hosted at ThFree Co, Russian Federation.

github.com/prusa3d/PrusaSlicer.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'src/libnest2d/tools/libnfporb/libnfporb.hpp')
-rw-r--r--src/libnest2d/tools/libnfporb/libnfporb.hpp1547
1 files changed, 1547 insertions, 0 deletions
diff --git a/src/libnest2d/tools/libnfporb/libnfporb.hpp b/src/libnest2d/tools/libnfporb/libnfporb.hpp
new file mode 100644
index 000000000..8cb34567e
--- /dev/null
+++ b/src/libnest2d/tools/libnfporb/libnfporb.hpp
@@ -0,0 +1,1547 @@
+#ifndef NFP_HPP_
+#define NFP_HPP_
+
+#include <iostream>
+#include <list>
+#include <string>
+#include <fstream>
+#include <streambuf>
+#include <vector>
+#include <set>
+#include <exception>
+#include <random>
+#include <limits>
+
+#if defined(_MSC_VER) && _MSC_VER <= 1800 || __cplusplus < 201103L
+ #define LIBNFP_NOEXCEPT
+ #define LIBNFP_CONSTEXPR
+#elif __cplusplus >= 201103L
+ #define LIBNFP_NOEXCEPT noexcept
+ #define LIBNFP_CONSTEXPR constexpr
+#endif
+
+#ifdef LIBNFP_USE_RATIONAL
+#include <boost/multiprecision/gmp.hpp>
+#include <boost/multiprecision/number.hpp>
+#endif
+#include <boost/geometry.hpp>
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/polygon.hpp>
+#include <boost/geometry/geometries/linestring.hpp>
+#include <boost/geometry/io/svg/svg_mapper.hpp>
+#include <boost/geometry/algorithms/intersects.hpp>
+#include <boost/geometry/geometries/register/point.hpp>
+
+#ifdef LIBNFP_USE_RATIONAL
+namespace bm = boost::multiprecision;
+#endif
+namespace bg = boost::geometry;
+namespace trans = boost::geometry::strategy::transform;
+
+
+namespace libnfporb {
+#ifdef NFP_DEBUG
+#define DEBUG_VAL(x) std::cerr << x << std::endl;
+#define DEBUG_MSG(title, value) std::cerr << title << ":" << value << std::endl;
+#else
+#define DEBUG_VAL(x)
+#define DEBUG_MSG(title, value)
+#endif
+
+using std::string;
+
+static LIBNFP_CONSTEXPR long double NFP_EPSILON=0.00000001;
+
+class LongDouble {
+private:
+ long double val_;
+public:
+ LongDouble() : val_(0) {
+ }
+
+ LongDouble(const long double& val) : val_(val) {
+ }
+
+ void setVal(const long double& v) {
+ val_ = v;
+ }
+
+ long double val() const {
+ return val_;
+ }
+
+ LongDouble operator/(const LongDouble& other) const {
+ return this->val_ / other.val_;
+ }
+
+ LongDouble operator*(const LongDouble& other) const {
+ return this->val_ * other.val_;
+ }
+
+ LongDouble operator-(const LongDouble& other) const {
+ return this->val_ - other.val_;
+ }
+
+ LongDouble operator-() const {
+ return this->val_ * -1;
+ }
+
+ LongDouble operator+(const LongDouble& other) const {
+ return this->val_ + other.val_;
+ }
+
+ void operator/=(const LongDouble& other) {
+ this->val_ = this->val_ / other.val_;
+ }
+
+ void operator*=(const LongDouble& other) {
+ this->val_ = this->val_ * other.val_;
+ }
+
+ void operator-=(const LongDouble& other) {
+ this->val_ = this->val_ - other.val_;
+ }
+
+ void operator+=(const LongDouble& other) {
+ this->val_ = this->val_ + other.val_;
+ }
+
+ bool operator==(const int& other) const {
+ return this->operator ==(static_cast<long double>(other));
+ }
+
+ bool operator==(const LongDouble& other) const {
+ return this->operator ==(other.val());
+ }
+
+ bool operator==(const long double& other) const {
+ return this->val() == other;
+ }
+
+ bool operator!=(const int& other) const {
+ return !this->operator ==(other);
+ }
+
+ bool operator!=(const LongDouble& other) const {
+ return !this->operator ==(other);
+ }
+
+ bool operator!=(const long double& other) const {
+ return !this->operator ==(other);
+ }
+
+ bool operator<(const int& other) const {
+ return this->operator <(static_cast<long double>(other));
+ }
+
+ bool operator<(const LongDouble& other) const {
+ return this->operator <(other.val());
+ }
+
+ bool operator<(const long double& other) const {
+ return this->val() < other;
+ }
+
+ bool operator>(const int& other) const {
+ return this->operator >(static_cast<long double>(other));
+ }
+
+ bool operator>(const LongDouble& other) const {
+ return this->operator >(other.val());
+ }
+
+ bool operator>(const long double& other) const {
+ return this->val() > other;
+ }
+
+ bool operator>=(const int& other) const {
+ return this->operator >=(static_cast<long double>(other));
+ }
+
+ bool operator>=(const LongDouble& other) const {
+ return this->operator >=(other.val());
+ }
+
+ bool operator>=(const long double& other) const {
+ return this->val() >= other;
+ }
+
+ bool operator<=(const int& other) const {
+ return this->operator <=(static_cast<long double>(other));
+ }
+
+ bool operator<=(const LongDouble& other) const {
+ return this->operator <=(other.val());
+ }
+
+ bool operator<=(const long double& other) const {
+ return this->val() <= other;
+ }
+};
+}
+
+
+namespace std {
+template<>
+ struct numeric_limits<libnfporb::LongDouble>
+ {
+ static const LIBNFP_CONSTEXPR bool is_specialized = true;
+
+ static const LIBNFP_CONSTEXPR long double
+ min() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::min(); }
+
+ static LIBNFP_CONSTEXPR long double
+ max() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::max(); }
+
+#if __cplusplus >= 201103L
+ static LIBNFP_CONSTEXPR long double
+ lowest() LIBNFP_NOEXCEPT { return -std::numeric_limits<long double>::lowest(); }
+#endif
+
+ static const LIBNFP_CONSTEXPR int digits = std::numeric_limits<long double>::digits;
+ static const LIBNFP_CONSTEXPR int digits10 = std::numeric_limits<long double>::digits10;
+#if __cplusplus >= 201103L
+ static const LIBNFP_CONSTEXPR int max_digits10
+ = std::numeric_limits<long double>::max_digits10;
+#endif
+ static const LIBNFP_CONSTEXPR bool is_signed = true;
+ static const LIBNFP_CONSTEXPR bool is_integer = false;
+ static const LIBNFP_CONSTEXPR bool is_exact = false;
+ static const LIBNFP_CONSTEXPR int radix = std::numeric_limits<long double>::radix;
+
+ static const LIBNFP_CONSTEXPR long double
+ epsilon() LIBNFP_NOEXCEPT { return libnfporb::NFP_EPSILON; }
+
+ static const LIBNFP_CONSTEXPR long double
+ round_error() LIBNFP_NOEXCEPT { return 0.5L; }
+
+ static const LIBNFP_CONSTEXPR int min_exponent = std::numeric_limits<long double>::min_exponent;
+ static const LIBNFP_CONSTEXPR int min_exponent10 = std::numeric_limits<long double>::min_exponent10;
+ static const LIBNFP_CONSTEXPR int max_exponent = std::numeric_limits<long double>::max_exponent;
+ static const LIBNFP_CONSTEXPR int max_exponent10 = std::numeric_limits<long double>::max_exponent10;
+
+
+ static const LIBNFP_CONSTEXPR bool has_infinity = std::numeric_limits<long double>::has_infinity;
+ static const LIBNFP_CONSTEXPR bool has_quiet_NaN = std::numeric_limits<long double>::has_quiet_NaN;
+ static const LIBNFP_CONSTEXPR bool has_signaling_NaN = has_quiet_NaN;
+ static const LIBNFP_CONSTEXPR float_denorm_style has_denorm
+ = std::numeric_limits<long double>::has_denorm;
+ static const LIBNFP_CONSTEXPR bool has_denorm_loss
+ = std::numeric_limits<long double>::has_denorm_loss;
+
+
+ static const LIBNFP_CONSTEXPR long double
+ infinity() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::infinity(); }
+
+ static const LIBNFP_CONSTEXPR long double
+ quiet_NaN() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::quiet_NaN(); }
+
+ static const LIBNFP_CONSTEXPR long double
+ signaling_NaN() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::signaling_NaN(); }
+
+
+ static const LIBNFP_CONSTEXPR long double
+ denorm_min() LIBNFP_NOEXCEPT { return std::numeric_limits<long double>::denorm_min(); }
+
+ static const LIBNFP_CONSTEXPR bool is_iec559
+ = has_infinity && has_quiet_NaN && has_denorm == denorm_present;
+
+ static const LIBNFP_CONSTEXPR bool is_bounded = true;
+ static const LIBNFP_CONSTEXPR bool is_modulo = false;
+
+ static const LIBNFP_CONSTEXPR bool traps = std::numeric_limits<long double>::traps;
+ static const LIBNFP_CONSTEXPR bool tinyness_before =
+ std::numeric_limits<long double>::tinyness_before;
+ static const LIBNFP_CONSTEXPR float_round_style round_style =
+ round_to_nearest;
+ };
+}
+
+namespace boost {
+namespace numeric {
+ template<>
+ struct raw_converter<boost::numeric::conversion_traits<double, libnfporb::LongDouble>>
+ {
+ typedef boost::numeric::conversion_traits<double, libnfporb::LongDouble>::result_type result_type ;
+ typedef boost::numeric::conversion_traits<double, libnfporb::LongDouble>::argument_type argument_type ;
+
+ static result_type low_level_convert ( argument_type s ) { return s.val() ; }
+ } ;
+}
+}
+
+namespace libnfporb {
+
+#ifndef LIBNFP_USE_RATIONAL
+typedef LongDouble coord_t;
+#else
+typedef bm::number<bm::gmp_rational, bm::et_off> rational_t;
+typedef rational_t coord_t;
+#endif
+
+bool equals(const LongDouble& lhs, const LongDouble& rhs);
+#ifdef LIBNFP_USE_RATIONAL
+bool equals(const rational_t& lhs, const rational_t& rhs);
+#endif
+bool equals(const long double& lhs, const long double& rhs);
+
+const coord_t MAX_COORD = 999999999999999999.0;
+const coord_t MIN_COORD = std::numeric_limits<coord_t>::min();
+
+class point_t {
+public:
+ point_t() : x_(0), y_(0) {
+ }
+ point_t(coord_t x, coord_t y) : x_(x), y_(y) {
+ }
+ bool marked_ = false;
+ coord_t x_;
+ coord_t y_;
+
+ point_t operator-(const point_t& other) const {
+ point_t result = *this;
+ bg::subtract_point(result, other);
+ return result;
+ }
+
+ point_t operator+(const point_t& other) const {
+ point_t result = *this;
+ bg::add_point(result, other);
+ return result;
+ }
+
+ bool operator==(const point_t& other) const {
+ return bg::equals(this, other);
+ }
+
+ bool operator!=(const point_t& other) const {
+ return !this->operator ==(other);
+ }
+
+ bool operator<(const point_t& other) const {
+ return boost::geometry::math::smaller(this->x_, other.x_) || (equals(this->x_, other.x_) && boost::geometry::math::smaller(this->y_, other.y_));
+ }
+};
+
+
+
+
+inline long double toLongDouble(const LongDouble& c) {
+ return c.val();
+}
+
+#ifdef LIBNFP_USE_RATIONAL
+inline long double toLongDouble(const rational_t& c) {
+ return bm::numerator(c).convert_to<long double>() / bm::denominator(c).convert_to<long double>();
+}
+#endif
+
+std::ostream& operator<<(std::ostream& os, const coord_t& p) {
+ os << toLongDouble(p);
+ return os;
+}
+
+std::istream& operator>>(std::istream& is, LongDouble& c) {
+ long double val;
+ is >> val;
+ c.setVal(val);
+ return is;
+}
+
+std::ostream& operator<<(std::ostream& os, const point_t& p) {
+ os << "{" << toLongDouble(p.x_) << "," << toLongDouble(p.y_) << "}";
+ return os;
+}
+const point_t INVALID_POINT = {MAX_COORD, MAX_COORD};
+
+typedef bg::model::segment<point_t> segment_t;
+}
+
+#ifdef LIBNFP_USE_RATIONAL
+inline long double acos(const libnfporb::rational_t& r) {
+ return acos(libnfporb::toLongDouble(r));
+}
+#endif
+
+inline long double acos(const libnfporb::LongDouble& ld) {
+ return acos(libnfporb::toLongDouble(ld));
+}
+
+#ifdef LIBNFP_USE_RATIONAL
+inline long double sqrt(const libnfporb::rational_t& r) {
+ return sqrt(libnfporb::toLongDouble(r));
+}
+#endif
+
+inline long double sqrt(const libnfporb::LongDouble& ld) {
+ return sqrt(libnfporb::toLongDouble(ld));
+}
+
+BOOST_GEOMETRY_REGISTER_POINT_2D(libnfporb::point_t, libnfporb::coord_t, cs::cartesian, x_, y_)
+
+
+namespace boost {
+namespace geometry {
+namespace math {
+namespace detail {
+
+template <>
+struct square_root<libnfporb::LongDouble>
+{
+ typedef libnfporb::LongDouble return_type;
+
+ static inline libnfporb::LongDouble apply(libnfporb::LongDouble const& a)
+ {
+ return std::sqrt(a.val());
+ }
+};
+
+#ifdef LIBNFP_USE_RATIONAL
+template <>
+struct square_root<libnfporb::rational_t>
+{
+ typedef libnfporb::rational_t return_type;
+
+ static inline libnfporb::rational_t apply(libnfporb::rational_t const& a)
+ {
+ return std::sqrt(libnfporb::toLongDouble(a));
+ }
+};
+#endif
+
+template<>
+struct abs<libnfporb::LongDouble>
+ {
+ static libnfporb::LongDouble apply(libnfporb::LongDouble const& value)
+ {
+ libnfporb::LongDouble const zero = libnfporb::LongDouble();
+ return value.val() < zero.val() ? -value.val() : value.val();
+ }
+ };
+
+template <>
+struct equals<libnfporb::LongDouble, false>
+{
+ template<typename Policy>
+ static inline bool apply(libnfporb::LongDouble const& lhs, libnfporb::LongDouble const& rhs, Policy const& policy)
+ {
+ if(lhs.val() == rhs.val())
+ return true;
+
+ return bg::math::detail::abs<libnfporb::LongDouble>::apply(lhs.val() - rhs.val()) <= policy.apply(lhs.val(), rhs.val()) * libnfporb::NFP_EPSILON;
+ }
+};
+
+template <>
+struct smaller<libnfporb::LongDouble>
+{
+ static inline bool apply(libnfporb::LongDouble const& lhs, libnfporb::LongDouble const& rhs)
+ {
+ if(lhs.val() == rhs.val() || bg::math::detail::abs<libnfporb::LongDouble>::apply(lhs.val() - rhs.val()) <= libnfporb::NFP_EPSILON * std::max(lhs.val(), rhs.val()))
+ return false;
+
+ return lhs < rhs;
+ }
+};
+}
+}
+}
+}
+
+namespace libnfporb {
+inline bool smaller(const LongDouble& lhs, const LongDouble& rhs) {
+ return boost::geometry::math::detail::smaller<LongDouble>::apply(lhs, rhs);
+}
+
+inline bool larger(const LongDouble& lhs, const LongDouble& rhs) {
+ return smaller(rhs, lhs);
+}
+
+bool equals(const LongDouble& lhs, const LongDouble& rhs) {
+ if(lhs.val() == rhs.val())
+ return true;
+
+ return bg::math::detail::abs<libnfporb::LongDouble>::apply(lhs.val() - rhs.val()) <= libnfporb::NFP_EPSILON * std::max(lhs.val(), rhs.val());
+}
+
+#ifdef LIBNFP_USE_RATIONAL
+inline bool smaller(const rational_t& lhs, const rational_t& rhs) {
+ return lhs < rhs;
+}
+
+inline bool larger(const rational_t& lhs, const rational_t& rhs) {
+ return smaller(rhs, lhs);
+}
+
+bool equals(const rational_t& lhs, const rational_t& rhs) {
+ return lhs == rhs;
+}
+#endif
+
+inline bool smaller(const long double& lhs, const long double& rhs) {
+ return lhs < rhs;
+}
+
+inline bool larger(const long double& lhs, const long double& rhs) {
+ return smaller(rhs, lhs);
+}
+
+
+bool equals(const long double& lhs, const long double& rhs) {
+ return lhs == rhs;
+}
+
+typedef bg::model::polygon<point_t, false, true> polygon_t;
+typedef std::vector<polygon_t::ring_type> nfp_t;
+typedef bg::model::linestring<point_t> linestring_t;
+
+typedef polygon_t::ring_type::size_type psize_t;
+
+typedef bg::model::d2::point_xy<long double> pointf_t;
+typedef bg::model::segment<pointf_t> segmentf_t;
+typedef bg::model::polygon<pointf_t, false, true> polygonf_t;
+
+polygonf_t::ring_type convert(const polygon_t::ring_type& r) {
+ polygonf_t::ring_type rf;
+ for(const auto& pt : r) {
+ rf.push_back(pointf_t(toLongDouble(pt.x_), toLongDouble(pt.y_)));
+ }
+ return rf;
+}
+
+polygonf_t convert(polygon_t p) {
+ polygonf_t pf;
+ pf.outer() = convert(p.outer());
+
+ for(const auto& r : p.inners()) {
+ pf.inners().push_back(convert(r));
+ }
+
+ return pf;
+}
+
+polygon_t nfpRingsToNfpPoly(const nfp_t& nfp) {
+ polygon_t nfppoly;
+ for (const auto& pt : nfp.front()) {
+ nfppoly.outer().push_back(pt);
+ }
+
+ for (size_t i = 1; i < nfp.size(); ++i) {
+ nfppoly.inners().push_back({});
+ for (const auto& pt : nfp[i]) {
+ nfppoly.inners().back().push_back(pt);
+ }
+ }
+
+ return nfppoly;
+}
+
+void write_svg(std::string const& filename,const std::vector<segment_t>& segments) {
+ std::ofstream svg(filename.c_str());
+
+ boost::geometry::svg_mapper<pointf_t> mapper(svg, 100, 100, "width=\"200mm\" height=\"200mm\" viewBox=\"-250 -250 500 500\"");
+ for(const auto& seg : segments) {
+ segmentf_t segf({toLongDouble(seg.first.x_), toLongDouble(seg.first.y_)}, {toLongDouble(seg.second.x_), toLongDouble(seg.second.y_)});
+ mapper.add(segf);
+ mapper.map(segf, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
+ }
+}
+
+void write_svg(std::string const& filename, const polygon_t& p, const polygon_t::ring_type& ring) {
+ std::ofstream svg(filename.c_str());
+
+ boost::geometry::svg_mapper<pointf_t> mapper(svg, 100, 100, "width=\"200mm\" height=\"200mm\" viewBox=\"-250 -250 500 500\"");
+ auto pf = convert(p);
+ auto rf = convert(ring);
+
+ mapper.add(pf);
+ mapper.map(pf, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
+ mapper.add(rf);
+ mapper.map(rf, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
+}
+
+void write_svg(std::string const& filename, std::vector<polygon_t> const& polygons) {
+ std::ofstream svg(filename.c_str());
+
+ boost::geometry::svg_mapper<pointf_t> mapper(svg, 100, 100, "width=\"200mm\" height=\"200mm\" viewBox=\"-250 -250 500 500\"");
+ for (auto p : polygons) {
+ auto pf = convert(p);
+ mapper.add(pf);
+ mapper.map(pf, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
+ }
+}
+
+void write_svg(std::string const& filename, std::vector<polygon_t> const& polygons, const nfp_t& nfp) {
+ polygon_t nfppoly;
+ for (const auto& pt : nfp.front()) {
+ nfppoly.outer().push_back(pt);
+ }
+
+ for (size_t i = 1; i < nfp.size(); ++i) {
+ nfppoly.inners().push_back({});
+ for (const auto& pt : nfp[i]) {
+ nfppoly.inners().back().push_back(pt);
+ }
+ }
+ std::ofstream svg(filename.c_str());
+
+ boost::geometry::svg_mapper<pointf_t> mapper(svg, 100, 100, "width=\"200mm\" height=\"200mm\" viewBox=\"-250 -250 500 500\"");
+ for (auto p : polygons) {
+ auto pf = convert(p);
+ mapper.add(pf);
+ mapper.map(pf, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
+ }
+ bg::correct(nfppoly);
+ auto nfpf = convert(nfppoly);
+ mapper.add(nfpf);
+ mapper.map(nfpf, "fill-opacity:0.5;fill:rgb(204,153,0);stroke:rgb(204,153,0);stroke-width:2");
+
+ for(auto& r: nfpf.inners()) {
+ if(r.size() == 1) {
+ mapper.add(r.front());
+ mapper.map(r.front(), "fill-opacity:0.5;fill:rgb(204,153,0);stroke:rgb(204,153,0);stroke-width:2");
+ } else if(r.size() == 2) {
+ segmentf_t seg(r.front(), *(r.begin()+1));
+ mapper.add(seg);
+ mapper.map(seg, "fill-opacity:0.5;fill:rgb(204,153,0);stroke:rgb(204,153,0);stroke-width:2");
+ }
+ }
+}
+
+std::ostream& operator<<(std::ostream& os, const segment_t& seg) {
+ os << "{" << seg.first << "," << seg.second << "}";
+ return os;
+}
+
+bool operator<(const segment_t& lhs, const segment_t& rhs) {
+ return lhs.first < rhs.first || ((lhs.first == rhs.first) && (lhs.second < rhs.second));
+}
+
+bool operator==(const segment_t& lhs, const segment_t& rhs) {
+ return (lhs.first == rhs.first && lhs.second == rhs.second) || (lhs.first == rhs.second && lhs.second == rhs.first);
+}
+
+bool operator!=(const segment_t& lhs, const segment_t& rhs) {
+ return !operator==(lhs,rhs);
+}
+
+enum Alignment {
+ LEFT,
+ RIGHT,
+ ON
+};
+
+point_t normalize(const point_t& pt) {
+ point_t norm = pt;
+ coord_t len = bg::length(segment_t{{0,0},pt});
+
+ if(len == 0.0L)
+ return {0,0};
+
+ norm.x_ /= len;
+ norm.y_ /= len;
+
+ return norm;
+}
+
+Alignment get_alignment(const segment_t& seg, const point_t& pt){
+ coord_t res = ((seg.second.x_ - seg.first.x_)*(pt.y_ - seg.first.y_)
+ - (seg.second.y_ - seg.first.y_)*(pt.x_ - seg.first.x_));
+
+ if(equals(res, 0)) {
+ return ON;
+ } else if(larger(res,0)) {
+ return LEFT;
+ } else {
+ return RIGHT;
+ }
+}
+
+long double get_inner_angle(const point_t& joint, const point_t& end1, const point_t& end2) {
+ coord_t dx21 = end1.x_-joint.x_;
+ coord_t dx31 = end2.x_-joint.x_;
+ coord_t dy21 = end1.y_-joint.y_;
+ coord_t dy31 = end2.y_-joint.y_;
+ coord_t m12 = sqrt((dx21*dx21 + dy21*dy21));
+ coord_t m13 = sqrt((dx31*dx31 + dy31*dy31));
+ if(m12 == 0.0L || m13 == 0.0L)
+ return 0;
+ return acos( (dx21*dx31 + dy21*dy31) / (m12 * m13) );
+}
+
+struct TouchingPoint {
+ enum Type {
+ VERTEX,
+ A_ON_B,
+ B_ON_A
+ };
+ Type type_;
+ psize_t A_;
+ psize_t B_;
+};
+
+struct TranslationVector {
+ point_t vector_;
+ segment_t edge_;
+ bool fromA_;
+ string name_;
+
+ bool operator<(const TranslationVector& other) const {
+ return this->vector_ < other.vector_ || ((this->vector_ == other.vector_) && (this->edge_ < other.edge_));
+ }
+};
+
+std::ostream& operator<<(std::ostream& os, const TranslationVector& tv) {
+ os << "{" << tv.edge_ << " -> " << tv.vector_ << "} = " << tv.name_;
+ return os;
+}
+
+
+void read_wkt_polygon(const string& filename, polygon_t& p) {
+ std::ifstream t(filename);
+
+ std::string str;
+ t.seekg(0, std::ios::end);
+ str.reserve(t.tellg());
+ t.seekg(0, std::ios::beg);
+
+ str.assign((std::istreambuf_iterator<char>(t)),
+ std::istreambuf_iterator<char>());
+
+ str.pop_back();
+ bg::read_wkt(str, p);
+ bg::correct(p);
+}
+
+std::vector<psize_t> find_minimum_y(const polygon_t& p) {
+ std::vector<psize_t> result;
+ coord_t min = MAX_COORD;
+ auto& po = p.outer();
+ for(psize_t i = 0; i < p.outer().size() - 1; ++i) {
+ if(smaller(po[i].y_, min)) {
+ result.clear();
+ min = po[i].y_;
+ result.push_back(i);
+ } else if (equals(po[i].y_, min)) {
+ result.push_back(i);
+ }
+ }
+ return result;
+}
+
+std::vector<psize_t> find_maximum_y(const polygon_t& p) {
+ std::vector<psize_t> result;
+ coord_t max = MIN_COORD;
+ auto& po = p.outer();
+ for(psize_t i = 0; i < p.outer().size() - 1; ++i) {
+ if(larger(po[i].y_, max)) {
+ result.clear();
+ max = po[i].y_;
+ result.push_back(i);
+ } else if (equals(po[i].y_, max)) {
+ result.push_back(i);
+ }
+ }
+ return result;
+}
+
+psize_t find_point(const polygon_t::ring_type& ring, const point_t& pt) {
+ for(psize_t i = 0; i < ring.size(); ++i) {
+ if(ring[i] == pt)
+ return i;
+ }
+ return std::numeric_limits<psize_t>::max();
+}
+
+std::vector<TouchingPoint> findTouchingPoints(const polygon_t::ring_type& ringA, const polygon_t::ring_type& ringB) {
+ std::vector<TouchingPoint> touchers;
+ for(psize_t i = 0; i < ringA.size() - 1; i++) {
+ psize_t nextI = i+1;
+ for(psize_t j = 0; j < ringB.size() - 1; j++) {
+ psize_t nextJ = j+1;
+ if(ringA[i] == ringB[j]) {
+ touchers.push_back({TouchingPoint::VERTEX, i, j});
+ } else if (ringA[nextI] != ringB[j] && bg::intersects(segment_t(ringA[i],ringA[nextI]), ringB[j])) {
+ touchers.push_back({TouchingPoint::B_ON_A, nextI, j});
+ } else if (ringB[nextJ] != ringA[i] && bg::intersects(segment_t(ringB[j],ringB[nextJ]), ringA[i])) {
+ touchers.push_back({TouchingPoint::A_ON_B, i, nextJ});
+ }
+ }
+ }
+ return touchers;
+}
+
+//TODO deduplicate code
+TranslationVector trimVector(const polygon_t::ring_type& rA, const polygon_t::ring_type& rB, const TranslationVector& tv) {
+ coord_t shortest = bg::length(tv.edge_);
+ TranslationVector trimmed = tv;
+ for(const auto& ptA : rA) {
+ point_t translated;
+ //for polygon A we invert the translation
+ trans::translate_transformer<coord_t, 2, 2> translate(-tv.vector_.x_, -tv.vector_.y_);
+ boost::geometry::transform(ptA, translated, translate);
+ linestring_t projection;
+ segment_t segproj(ptA, translated);
+ projection.push_back(ptA);
+ projection.push_back(translated);
+ std::vector<point_t> intersections;
+ bg::intersection(rB, projection, intersections);
+ if(bg::touches(projection, rB) && intersections.size() < 2) {
+ continue;
+ }
+
+ //find shortest intersection
+ coord_t len;
+ segment_t segi;
+ for(const auto& pti : intersections) {
+ segi = segment_t(ptA,pti);
+ len = bg::length(segi);
+ if(smaller(len, shortest)) {
+ trimmed.vector_ = ptA - pti;
+ trimmed.edge_ = segi;
+ shortest = len;
+ }
+ }
+ }
+
+ for(const auto& ptB : rB) {
+ point_t translated;
+
+ trans::translate_transformer<coord_t, 2, 2> translate(tv.vector_.x_, tv.vector_.y_);
+ boost::geometry::transform(ptB, translated, translate);
+ linestring_t projection;
+ segment_t segproj(ptB, translated);
+ projection.push_back(ptB);
+ projection.push_back(translated);
+ std::vector<point_t> intersections;
+ bg::intersection(rA, projection, intersections);
+ if(bg::touches(projection, rA) && intersections.size() < 2) {
+ continue;
+ }
+
+ //find shortest intersection
+ coord_t len;
+ segment_t segi;
+ for(const auto& pti : intersections) {
+
+ segi = segment_t(ptB,pti);
+ len = bg::length(segi);
+ if(smaller(len, shortest)) {
+ trimmed.vector_ = pti - ptB;
+ trimmed.edge_ = segi;
+ shortest = len;
+ }
+ }
+ }
+ return trimmed;
+}
+
+std::vector<TranslationVector> findFeasibleTranslationVectors(polygon_t::ring_type& ringA, polygon_t::ring_type& ringB, const std::vector<TouchingPoint>& touchers) {
+ //use a set to automatically filter duplicate vectors
+ std::vector<TranslationVector> potentialVectors;
+ std::vector<std::pair<segment_t,segment_t>> touchEdges;
+
+ for (psize_t i = 0; i < touchers.size(); i++) {
+ point_t& vertexA = ringA[touchers[i].A_];
+ vertexA.marked_ = true;
+
+ // adjacent A vertices
+ auto prevAindex = static_cast<signed long>(touchers[i].A_ - 1);
+ auto nextAindex = static_cast<signed long>(touchers[i].A_ + 1);
+
+ prevAindex = (prevAindex < 0) ? static_cast<signed long>(ringA.size() - 2) : prevAindex; // loop
+ nextAindex = (static_cast<psize_t>(nextAindex) >= ringA.size()) ? 1 : nextAindex; // loop
+
+ point_t& prevA = ringA[prevAindex];
+ point_t& nextA = ringA[nextAindex];
+
+ // adjacent B vertices
+ point_t& vertexB = ringB[touchers[i].B_];
+
+ auto prevBindex = static_cast<signed long>(touchers[i].B_ - 1);
+ auto nextBindex = static_cast<signed long>(touchers[i].B_ + 1);
+
+ prevBindex = (prevBindex < 0) ? static_cast<signed long>(ringB.size() - 2) : prevBindex; // loop
+ nextBindex = (static_cast<psize_t>(nextBindex) >= ringB.size()) ? 1 : nextBindex; // loop
+
+ point_t& prevB = ringB[prevBindex];
+ point_t& nextB = ringB[nextBindex];
+
+ if (touchers[i].type_ == TouchingPoint::VERTEX) {
+ segment_t a1 = { vertexA, nextA };
+ segment_t a2 = { vertexA, prevA };
+ segment_t b1 = { vertexB, nextB };
+ segment_t b2 = { vertexB, prevB };
+
+ //swap the segment elements so that always the first point is the touching point
+ //also make the second segment always a segment of ringB
+ touchEdges.push_back({a1, b1});
+ touchEdges.push_back({a1, b2});
+ touchEdges.push_back({a2, b1});
+ touchEdges.push_back({a2, b2});
+#ifdef NFP_DEBUG
+ write_svg("touchersV" + std::to_string(i) + ".svg", {a1,a2,b1,b2});
+#endif
+
+ //TODO test parallel edges for floating point stability
+ Alignment al;
+ //a1 and b1 meet at start vertex
+ al = get_alignment(a1, b1.second);
+ if(al == LEFT) {
+ potentialVectors.push_back({b1.first - b1.second, b1, false, "vertex1"});
+ } else if(al == RIGHT) {
+ potentialVectors.push_back({a1.second - a1.first, a1, true, "vertex2"});
+ } else {
+ potentialVectors.push_back({a1.second - a1.first, a1, true, "vertex3"});
+ }
+
+ //a1 and b2 meet at start and end
+ al = get_alignment(a1, b2.second);
+ if(al == LEFT) {
+ //no feasible translation
+ } else if(al == RIGHT) {
+ potentialVectors.push_back({a1.second - a1.first, a1, true, "vertex4"});
+ } else {
+ potentialVectors.push_back({a1.second - a1.first, a1, true, "vertex5"});
+ }
+
+ //a2 and b1 meet at end and start
+ al = get_alignment(a2, b1.second);
+ if(al == LEFT) {
+ //no feasible translation
+ } else if(al == RIGHT) {
+ potentialVectors.push_back({b1.first - b1.second, b1, false, "vertex6"});
+ } else {
+ potentialVectors.push_back({b1.first - b1.second, b1, false, "vertex7"});
+ }
+ } else if (touchers[i].type_ == TouchingPoint::B_ON_A) {
+ segment_t a1 = {vertexB, vertexA};
+ segment_t a2 = {vertexB, prevA};
+ segment_t b1 = {vertexB, prevB};
+ segment_t b2 = {vertexB, nextB};
+
+ touchEdges.push_back({a1, b1});
+ touchEdges.push_back({a1, b2});
+ touchEdges.push_back({a2, b1});
+ touchEdges.push_back({a2, b2});
+#ifdef NFP_DEBUG
+ write_svg("touchersB" + std::to_string(i) + ".svg", {a1,a2,b1,b2});
+#endif
+ potentialVectors.push_back({vertexA - vertexB, {vertexB, vertexA}, true, "bona"});
+ } else if (touchers[i].type_ == TouchingPoint::A_ON_B) {
+ //TODO testme
+ segment_t a1 = {vertexA, prevA};
+ segment_t a2 = {vertexA, nextA};
+ segment_t b1 = {vertexA, vertexB};
+ segment_t b2 = {vertexA, prevB};
+#ifdef NFP_DEBUG
+ write_svg("touchersA" + std::to_string(i) + ".svg", {a1,a2,b1,b2});
+#endif
+ touchEdges.push_back({a1, b1});
+ touchEdges.push_back({a2, b1});
+ touchEdges.push_back({a1, b2});
+ touchEdges.push_back({a2, b2});
+ potentialVectors.push_back({vertexA - vertexB, {vertexA, vertexB}, false, "aonb"});
+ }
+ }
+
+ //discard immediately intersecting translations
+ std::vector<TranslationVector> vectors;
+ for(const auto& v : potentialVectors) {
+ bool discarded = false;
+ for(const auto& sp : touchEdges) {
+ point_t normEdge = normalize(v.edge_.second - v.edge_.first);
+ point_t normFirst = normalize(sp.first.second - sp.first.first);
+ point_t normSecond = normalize(sp.second.second - sp.second.first);
+
+ Alignment a1 = get_alignment({{0,0},normEdge}, normFirst);
+ Alignment a2 = get_alignment({{0,0},normEdge}, normSecond);
+
+ if(a1 == a2 && a1 != ON) {
+ long double df = get_inner_angle({0,0},normEdge, normFirst);
+ long double ds = get_inner_angle({0,0},normEdge, normSecond);
+
+ point_t normIn = normalize(v.edge_.second - v.edge_.first);
+ if (equals(df, ds)) {
+ TranslationVector trimmed = trimVector(ringA,ringB, v);
+ polygon_t::ring_type translated;
+ trans::translate_transformer<coord_t, 2, 2> translate(trimmed.vector_.x_, trimmed.vector_.y_);
+ boost::geometry::transform(ringB, translated, translate);
+ if (!(bg::intersects(translated, ringA) && !bg::overlaps(translated, ringA) && !bg::covered_by(translated, ringA) && !bg::covered_by(ringA, translated))) {
+ discarded = true;
+ break;
+ }
+ } else {
+
+ if (normIn == normalize(v.vector_)) {
+ if (larger(ds, df)) {
+ discarded = true;
+ break;
+ }
+ } else {
+ if (smaller(ds, df)) {
+ discarded = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+ if(!discarded)
+ vectors.push_back(v);
+ }
+ return vectors;
+}
+
+bool find(const std::vector<TranslationVector>& h, const TranslationVector& tv) {
+ for(const auto& htv : h) {
+ if(htv.vector_ == tv.vector_)
+ return true;
+ }
+ return false;
+}
+
+TranslationVector getLongest(const std::vector<TranslationVector>& tvs) {
+ coord_t len;
+ coord_t maxLen = MIN_COORD;
+ TranslationVector longest;
+ longest.vector_ = INVALID_POINT;
+
+ for(auto& tv : tvs) {
+ len = bg::length(segment_t{{0,0},tv.vector_});
+ if(larger(len, maxLen)) {
+ maxLen = len;
+ longest = tv;
+ }
+ }
+ return longest;
+}
+
+TranslationVector selectNextTranslationVector(const polygon_t& pA, const polygon_t::ring_type& rA, const polygon_t::ring_type& rB, const std::vector<TranslationVector>& tvs, const std::vector<TranslationVector>& history) {
+ if(!history.empty()) {
+ TranslationVector last = history.back();
+ std::vector<TranslationVector> historyCopy = history;
+ if(historyCopy.size() >= 2) {
+ historyCopy.erase(historyCopy.end() - 1);
+ historyCopy.erase(historyCopy.end() - 1);
+ if(historyCopy.size() > 4) {
+ historyCopy.erase(historyCopy.begin(), historyCopy.end() - 4);
+ }
+
+ } else {
+ historyCopy.clear();
+ }
+ DEBUG_MSG("last", last);
+
+ psize_t laterI = std::numeric_limits<psize_t>::max();
+ point_t previous = rA[0];
+ point_t next;
+
+ if(last.fromA_) {
+ for (psize_t i = 1; i < rA.size() + 1; ++i) {
+ if (i >= rA.size())
+ next = rA[i % rA.size()];
+ else
+ next = rA[i];
+
+ segment_t candidate( previous, next );
+ if(candidate == last.edge_) {
+ laterI = i;
+ break;
+ }
+ previous = next;
+ }
+
+ if (laterI == std::numeric_limits<psize_t>::max()) {
+ point_t later;
+ if (last.vector_ == (last.edge_.second - last.edge_.first)) {
+ later = last.edge_.second;
+ } else {
+ later = last.edge_.first;
+ }
+
+ laterI = find_point(rA, later);
+ }
+ } else {
+ point_t later;
+ if (last.vector_ == (last.edge_.second - last.edge_.first)) {
+ later = last.edge_.second;
+ } else {
+ later = last.edge_.first;
+ }
+
+ laterI = find_point(rA, later);
+ }
+
+ if (laterI == std::numeric_limits<psize_t>::max()) {
+ throw std::runtime_error(
+ "Internal error: Can't find later point of last edge");
+ }
+
+ std::vector<segment_t> viableEdges;
+ previous = rA[laterI];
+ for(psize_t i = laterI + 1; i < rA.size() + laterI + 1; ++i) {
+ if(i >= rA.size())
+ next = rA[i % rA.size()];
+ else
+ next = rA[i];
+
+ viableEdges.push_back({previous, next});
+ previous = next;
+ }
+
+// auto rng = std::default_random_engine {};
+// std::shuffle(std::begin(viableEdges), std::end(viableEdges), rng);
+
+ //search with consulting the history to prevent oscillation
+ std::vector<TranslationVector> viableTrans;
+ for(const auto& ve: viableEdges) {
+ for(const auto& tv : tvs) {
+ if((tv.fromA_ && (normalize(tv.vector_) == normalize(ve.second - ve.first))) && (tv.edge_ != last.edge_ || tv.vector_.x_ != -last.vector_.x_ || tv.vector_.y_ != -last.vector_.y_) && !find(historyCopy, tv)) {
+ viableTrans.push_back(tv);
+ }
+ }
+ for (const auto& tv : tvs) {
+ if (!tv.fromA_) {
+ point_t later;
+ if (tv.vector_ == (tv.edge_.second - tv.edge_.first) && (tv.edge_ != last.edge_ || tv.vector_.x_ != -last.vector_.x_ || tv.vector_.y_ != -last.vector_.y_) && !find(historyCopy, tv)) {
+ later = tv.edge_.second;
+ } else if (tv.vector_ == (tv.edge_.first - tv.edge_.second)) {
+ later = tv.edge_.first;
+ } else
+ continue;
+
+ if (later == ve.first || later == ve.second) {
+ viableTrans.push_back(tv);
+ }
+ }
+ }
+ }
+
+ if(!viableTrans.empty())
+ return getLongest(viableTrans);
+
+ //search again without the history
+ for(const auto& ve: viableEdges) {
+ for(const auto& tv : tvs) {
+ if((tv.fromA_ && (normalize(tv.vector_) == normalize(ve.second - ve.first))) && (tv.edge_ != last.edge_ || tv.vector_.x_ != -last.vector_.x_ || tv.vector_.y_ != -last.vector_.y_)) {
+ viableTrans.push_back(tv);
+ }
+ }
+ for (const auto& tv : tvs) {
+ if (!tv.fromA_) {
+ point_t later;
+ if (tv.vector_ == (tv.edge_.second - tv.edge_.first) && (tv.edge_ != last.edge_ || tv.vector_.x_ != -last.vector_.x_ || tv.vector_.y_ != -last.vector_.y_)) {
+ later = tv.edge_.second;
+ } else if (tv.vector_ == (tv.edge_.first - tv.edge_.second)) {
+ later = tv.edge_.first;
+ } else
+ continue;
+
+ if (later == ve.first || later == ve.second) {
+ viableTrans.push_back(tv);
+ }
+ }
+ }
+ }
+ if(!viableTrans.empty())
+ return getLongest(viableTrans);
+
+ /*
+ //search again without the history and without checking last edge
+ for(const auto& ve: viableEdges) {
+ for(const auto& tv : tvs) {
+ if((tv.fromA_ && (normalize(tv.vector_) == normalize(ve.second - ve.first)))) {
+ return tv;
+ }
+ }
+ for (const auto& tv : tvs) {
+ if (!tv.fromA_) {
+ point_t later;
+ if (tv.vector_ == (tv.edge_.second - tv.edge_.first)) {
+ later = tv.edge_.second;
+ } else if (tv.vector_ == (tv.edge_.first - tv.edge_.second)) {
+ later = tv.edge_.first;
+ } else
+ continue;
+
+ if (later == ve.first || later == ve.second) {
+ return tv;
+ }
+ }
+ }
+ }*/
+
+ if(tvs.size() == 1)
+ return *tvs.begin();
+
+ TranslationVector tv;
+ tv.vector_ = INVALID_POINT;
+ return tv;
+ } else {
+ return getLongest(tvs);
+ }
+}
+
+bool inNfp(const point_t& pt, const nfp_t& nfp) {
+ for(const auto& r : nfp) {
+ if(bg::touches(pt, r))
+ return true;
+ }
+
+ return false;
+}
+
+enum SearchStartResult {
+ FIT,
+ FOUND,
+ NOT_FOUND
+};
+
+SearchStartResult searchStartTranslation(polygon_t::ring_type& rA, const polygon_t::ring_type& rB, const nfp_t& nfp,const bool& inside, point_t& result) {
+ for(psize_t i = 0; i < rA.size() - 1; i++) {
+ psize_t index;
+ if (i >= rA.size())
+ index = i % rA.size() + 1;
+ else
+ index = i;
+
+ auto& ptA = rA[index];
+
+ if(ptA.marked_)
+ continue;
+
+ ptA.marked_ = true;
+
+ for(const auto& ptB: rB) {
+ point_t testTranslation = ptA - ptB;
+ polygon_t::ring_type translated;
+ boost::geometry::transform(rB, translated, trans::translate_transformer<coord_t, 2, 2>(testTranslation.x_, testTranslation.y_));
+
+ //check if the translated rB is identical to rA
+ bool identical = false;
+ for(const auto& ptT: translated) {
+ identical = false;
+ for(const auto& ptA: rA) {
+ if(ptT == ptA) {
+ identical = true;
+ break;
+ }
+ }
+ if(!identical)
+ break;
+ }
+
+ if(identical) {
+ result = testTranslation;
+ return FIT;
+ }
+
+ bool bInside = false;
+ for(const auto& ptT: translated) {
+ if(bg::within(ptT, rA)) {
+ bInside = true;
+ break;
+ } else if(!bg::touches(ptT, rA)) {
+ bInside = false;
+ break;
+ }
+ }
+
+ if(((bInside && inside) || (!bInside && !inside)) && (!bg::overlaps(translated, rA) && !bg::covered_by(translated, rA) && !bg::covered_by(rA, translated)) && !inNfp(translated.front(), nfp)){
+ result = testTranslation;
+ return FOUND;
+ }
+
+ point_t nextPtA = rA[index + 1];
+ TranslationVector slideVector;
+ slideVector.vector_ = nextPtA - ptA;
+ slideVector.edge_ = {ptA, nextPtA};
+ slideVector.fromA_ = true;
+ TranslationVector trimmed = trimVector(rA, translated, slideVector);
+ polygon_t::ring_type translated2;
+ trans::translate_transformer<coord_t, 2, 2> trans(trimmed.vector_.x_, trimmed.vector_.y_);
+ boost::geometry::transform(translated, translated2, trans);
+
+ //check if the translated rB is identical to rA
+ identical = false;
+ for(const auto& ptT: translated) {
+ identical = false;
+ for(const auto& ptA: rA) {
+ if(ptT == ptA) {
+ identical = true;
+ break;
+ }
+ }
+ if(!identical)
+ break;
+ }
+
+ if(identical) {
+ result = trimmed.vector_ + testTranslation;
+ return FIT;
+ }
+
+ bInside = false;
+ for(const auto& ptT: translated2) {
+ if(bg::within(ptT, rA)) {
+ bInside = true;
+ break;
+ } else if(!bg::touches(ptT, rA)) {
+ bInside = false;
+ break;
+ }
+ }
+
+ if(((bInside && inside) || (!bInside && !inside)) && (!bg::overlaps(translated2, rA) && !bg::covered_by(translated2, rA) && !bg::covered_by(rA, translated2)) && !inNfp(translated2.front(), nfp)){
+ result = trimmed.vector_ + testTranslation;
+ return FOUND;
+ }
+ }
+ }
+ return NOT_FOUND;
+}
+
+enum SlideResult {
+ LOOP,
+ NO_LOOP,
+ NO_TRANSLATION
+};
+
+SlideResult slide(polygon_t& pA, polygon_t::ring_type& rA, polygon_t::ring_type& rB, nfp_t& nfp, const point_t& transB, bool inside) {
+ polygon_t::ring_type rifsB;
+ boost::geometry::transform(rB, rifsB, trans::translate_transformer<coord_t, 2, 2>(transB.x_, transB.y_));
+ rB = std::move(rifsB);
+
+#ifdef NFP_DEBUG
+ write_svg("ifs.svg", pA, rB);
+#endif
+
+ bool startAvailable = true;
+ psize_t cnt = 0;
+ point_t referenceStart = rB.front();
+ std::vector<TranslationVector> history;
+
+ //generate the nfp for the ring
+ while(startAvailable) {
+ DEBUG_VAL(cnt);
+ //use first point of rB as reference
+ nfp.back().push_back(rB.front());
+ if(cnt == 15)
+ std::cerr << "";
+
+ std::vector<TouchingPoint> touchers = findTouchingPoints(rA, rB);
+
+#ifdef NFP_DEBUG
+ DEBUG_MSG("touchers", touchers.size());
+ for(auto t : touchers) {
+ DEBUG_VAL(t.type_);
+ }
+#endif
+ if(touchers.empty()) {
+ throw std::runtime_error("Internal error: No touching points found");
+ }
+ std::vector<TranslationVector> transVectors = findFeasibleTranslationVectors(rA, rB, touchers);
+
+#ifdef NFP_DEBUG
+ DEBUG_MSG("collected vectors", transVectors.size());
+ for(auto pt : transVectors) {
+ DEBUG_VAL(pt);
+ }
+#endif
+
+ if(transVectors.empty()) {
+ return NO_LOOP;
+ }
+
+ TranslationVector next = selectNextTranslationVector(pA, rA, rB, transVectors, history);
+
+ if(next.vector_ == INVALID_POINT)
+ return NO_TRANSLATION;
+
+ DEBUG_MSG("next", next);
+
+ TranslationVector trimmed = trimVector(rA, rB, next);
+ DEBUG_MSG("trimmed", trimmed);
+
+ history.push_back(next);
+
+ polygon_t::ring_type nextRB;
+ boost::geometry::transform(rB, nextRB, trans::translate_transformer<coord_t, 2, 2>(trimmed.vector_.x_, trimmed.vector_.y_));
+ rB = std::move(nextRB);
+
+#ifdef NFP_DEBUG
+ write_svg("next" + std::to_string(cnt) + ".svg", pA,rB);
+#endif
+
+ ++cnt;
+ if(referenceStart == rB.front() || (inside && bg::touches(rB.front(), nfp.front()))) {
+ startAvailable = false;
+ }
+ }
+ return LOOP;
+}
+
+void removeCoLinear(polygon_t::ring_type& r) {
+ assert(r.size() > 2);
+ psize_t nextI;
+ psize_t prevI = 0;
+ segment_t segment(r[r.size() - 2], r[0]);
+ polygon_t::ring_type newR;
+
+ for (psize_t i = 1; i < r.size() + 1; ++i) {
+ if (i >= r.size())
+ nextI = i % r.size() + 1;
+ else
+ nextI = i;
+
+ if (get_alignment(segment, r[nextI]) != ON) {
+ newR.push_back(r[prevI]);
+ }
+ segment = {segment.second, r[nextI]};
+ prevI = nextI;
+ }
+
+ r = newR;
+}
+
+void removeCoLinear(polygon_t& p) {
+ removeCoLinear(p.outer());
+ for (auto& r : p.inners())
+ removeCoLinear(r);
+}
+
+nfp_t generateNFP(polygon_t& pA, polygon_t& pB, const bool checkValidity = true) {
+ removeCoLinear(pA);
+ removeCoLinear(pB);
+
+ if(checkValidity) {
+ std::string reason;
+ if(!bg::is_valid(pA, reason))
+ throw std::runtime_error("Polygon A is invalid: " + reason);
+
+ if(!bg::is_valid(pB, reason))
+ throw std::runtime_error("Polygon B is invalid: " + reason);
+ }
+
+ nfp_t nfp;
+
+#ifdef NFP_DEBUG
+ write_svg("start.svg", {pA, pB});
+#endif
+
+ DEBUG_VAL(bg::wkt(pA))
+ DEBUG_VAL(bg::wkt(pB));
+
+ //prevent double vertex connections at start because we might come back the same way we go which would end the nfp prematurely
+ std::vector<psize_t> ptyaminI = find_minimum_y(pA);
+ std::vector<psize_t> ptybmaxI = find_maximum_y(pB);
+
+ point_t pAstart;
+ point_t pBstart;
+
+ if(ptyaminI.size() > 1 || ptybmaxI.size() > 1) {
+ //find right-most of A and left-most of B to prevent double connection at start
+ coord_t maxX = MIN_COORD;
+ psize_t iRightMost = 0;
+ for(psize_t& ia : ptyaminI) {
+ const point_t& candidateA = pA.outer()[ia];
+ if(larger(candidateA.x_, maxX)) {
+ maxX = candidateA.x_;
+ iRightMost = ia;
+ }
+ }
+
+ coord_t minX = MAX_COORD;
+ psize_t iLeftMost = 0;
+ for(psize_t& ib : ptybmaxI) {
+ const point_t& candidateB = pB.outer()[ib];
+ if(smaller(candidateB.x_, minX)) {
+ minX = candidateB.x_;
+ iLeftMost = ib;
+ }
+ }
+ pAstart = pA.outer()[iRightMost];
+ pBstart = pB.outer()[iLeftMost];
+ } else {
+ pAstart = pA.outer()[ptyaminI.front()];
+ pBstart = pB.outer()[ptybmaxI.front()];
+ }
+
+ nfp.push_back({});
+ point_t transB = {pAstart - pBstart};
+
+ if(slide(pA, pA.outer(), pB.outer(), nfp, transB, false) != LOOP) {
+ throw std::runtime_error("Unable to complete outer nfp loop");
+ }
+
+ DEBUG_VAL("##### outer #####");
+ point_t startTrans;
+ while(true) {
+ SearchStartResult res = searchStartTranslation(pA.outer(), pB.outer(), nfp, false, startTrans);
+ if(res == FOUND) {
+ nfp.push_back({});
+ DEBUG_VAL("##### interlock start #####")
+ polygon_t::ring_type rifsB;
+ boost::geometry::transform(pB.outer(), rifsB, trans::translate_transformer<coord_t, 2, 2>(startTrans.x_, startTrans.y_));
+ if(inNfp(rifsB.front(), nfp)) {
+ continue;
+ }
+ SlideResult sres = slide(pA, pA.outer(), pB.outer(), nfp, startTrans, true);
+ if(sres != LOOP) {
+ if(sres == NO_TRANSLATION) {
+ //no initial slide found -> jiggsaw
+ if(!inNfp(pB.outer().front(),nfp)) {
+ nfp.push_back({});
+ nfp.back().push_back(pB.outer().front());
+ }
+ }
+ }
+ DEBUG_VAL("##### interlock end #####");
+ } else if(res == FIT) {
+ point_t reference = pB.outer().front();
+ point_t translated;
+ trans::translate_transformer<coord_t, 2, 2> translate(startTrans.x_, startTrans.y_);
+ boost::geometry::transform(reference, translated, translate);
+ if(!inNfp(translated,nfp)) {
+ nfp.push_back({});
+ nfp.back().push_back(translated);
+ }
+ break;
+ } else
+ break;
+ }
+
+
+ for(auto& rA : pA.inners()) {
+ while(true) {
+ SearchStartResult res = searchStartTranslation(rA, pB.outer(), nfp, true, startTrans);
+ if(res == FOUND) {
+ nfp.push_back({});
+ DEBUG_VAL("##### hole start #####");
+ slide(pA, rA, pB.outer(), nfp, startTrans, true);
+ DEBUG_VAL("##### hole end #####");
+ } else if(res == FIT) {
+ point_t reference = pB.outer().front();
+ point_t translated;
+ trans::translate_transformer<coord_t, 2, 2> translate(startTrans.x_, startTrans.y_);
+ boost::geometry::transform(reference, translated, translate);
+ if(!inNfp(translated,nfp)) {
+ nfp.push_back({});
+ nfp.back().push_back(translated);
+ }
+ break;
+ } else
+ break;
+ }
+ }
+
+#ifdef NFP_DEBUG
+ write_svg("nfp.svg", {pA,pB}, nfp);
+#endif
+
+ return nfp;
+}
+}
+#endif