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/nfp_svgnest.hpp')
-rw-r--r--src/libnest2d/tools/nfp_svgnest.hpp1018
1 files changed, 1018 insertions, 0 deletions
diff --git a/src/libnest2d/tools/nfp_svgnest.hpp b/src/libnest2d/tools/nfp_svgnest.hpp
new file mode 100644
index 000000000..ac5700c10
--- /dev/null
+++ b/src/libnest2d/tools/nfp_svgnest.hpp
@@ -0,0 +1,1018 @@
+#ifndef NFP_SVGNEST_HPP
+#define NFP_SVGNEST_HPP
+
+#include <limits>
+#include <unordered_map>
+
+#include <libnest2d/geometry_traits_nfp.hpp>
+
+namespace libnest2d {
+
+namespace __svgnest {
+
+using std::sqrt;
+using std::min;
+using std::max;
+using std::abs;
+using std::isnan;
+
+//template<class Coord> struct _Scale {
+// static const BP2D_CONSTEXPR long long Value = 1000000;
+//};
+
+template<class S> struct _alg {
+ using Contour = TContour<S>;
+ using Point = TPoint<S>;
+ using iCoord = TCoord<Point>;
+ using Coord = double;
+ using Shapes = nfp::Shapes<S>;
+
+ static const Coord TOL;
+
+#define dNAN std::nan("")
+
+ struct Vector {
+ Coord x = 0.0, y = 0.0;
+ bool marked = false;
+ Vector() = default;
+ Vector(Coord X, Coord Y): x(X), y(Y) {}
+ Vector(const Point& p): x(Coord(getX(p))), y(Coord(getY(p))) {}
+ operator Point() const { return {iCoord(x), iCoord(y)}; }
+ Vector& operator=(const Point& p) {
+ x = getX(p), y = getY(p); return *this;
+ }
+ bool operator!=(const Vector& v) const {
+ return v.x != x || v.y != y;
+ }
+ Vector(std::initializer_list<Coord> il):
+ x(*il.begin()), y(*std::next(il.begin())) {}
+ };
+
+ static inline Coord x(const Point& p) { return Coord(getX(p)); }
+ static inline Coord y(const Point& p) { return Coord(getY(p)); }
+
+ static inline Coord x(const Vector& p) { return p.x; }
+ static inline Coord y(const Vector& p) { return p.y; }
+
+ class Cntr {
+ std::vector<Vector> v_;
+ public:
+ Cntr(const Contour& c) {
+ v_.reserve(c.size());
+ std::transform(c.begin(), c.end(), std::back_inserter(v_),
+ [](const Point& p) {
+ return Vector(double(x(p)) / 1e6, double(y(p)) / 1e6);
+ });
+ std::reverse(v_.begin(), v_.end());
+ v_.pop_back();
+ }
+ Cntr() = default;
+
+ Coord offsetx = 0;
+ Coord offsety = 0;
+ size_t size() const { return v_.size(); }
+ bool empty() const { return v_.empty(); }
+ typename std::vector<Vector>::const_iterator cbegin() const { return v_.cbegin(); }
+ typename std::vector<Vector>::const_iterator cend() const { return v_.cend(); }
+ typename std::vector<Vector>::iterator begin() { return v_.begin(); }
+ typename std::vector<Vector>::iterator end() { return v_.end(); }
+ Vector& operator[](size_t idx) { return v_[idx]; }
+ const Vector& operator[](size_t idx) const { return v_[idx]; }
+ template<class...Args>
+ void emplace_back(Args&&...args) {
+ v_.emplace_back(std::forward<Args>(args)...);
+ }
+ template<class...Args>
+ void push(Args&&...args) {
+ v_.emplace_back(std::forward<Args>(args)...);
+ }
+ void clear() { v_.clear(); }
+
+ operator Contour() const {
+ Contour cnt;
+ cnt.reserve(v_.size() + 1);
+ std::transform(v_.begin(), v_.end(), std::back_inserter(cnt),
+ [](const Vector& vertex) {
+ return Point(iCoord(vertex.x) * 1000000, iCoord(vertex.y) * 1000000);
+ });
+ if(!cnt.empty()) cnt.emplace_back(cnt.front());
+ S sh = shapelike::create<S>(cnt);
+
+// std::reverse(cnt.begin(), cnt.end());
+ return shapelike::getContour(sh);
+ }
+ };
+
+ inline static bool _almostEqual(Coord a, Coord b,
+ Coord tolerance = TOL)
+ {
+ return std::abs(a - b) < tolerance;
+ }
+
+ // returns true if p lies on the line segment defined by AB,
+ // but not at any endpoints may need work!
+ static bool _onSegment(const Vector& A, const Vector& B, const Vector& p) {
+
+ // vertical line
+ if(_almostEqual(A.x, B.x) && _almostEqual(p.x, A.x)) {
+ if(!_almostEqual(p.y, B.y) && !_almostEqual(p.y, A.y) &&
+ p.y < max(B.y, A.y) && p.y > min(B.y, A.y)){
+ return true;
+ }
+ else{
+ return false;
+ }
+ }
+
+ // horizontal line
+ if(_almostEqual(A.y, B.y) && _almostEqual(p.y, A.y)){
+ if(!_almostEqual(p.x, B.x) && !_almostEqual(p.x, A.x) &&
+ p.x < max(B.x, A.x) && p.x > min(B.x, A.x)){
+ return true;
+ }
+ else{
+ return false;
+ }
+ }
+
+ //range check
+ if((p.x < A.x && p.x < B.x) || (p.x > A.x && p.x > B.x) ||
+ (p.y < A.y && p.y < B.y) || (p.y > A.y && p.y > B.y))
+ return false;
+
+ // exclude end points
+ if((_almostEqual(p.x, A.x) && _almostEqual(p.y, A.y)) ||
+ (_almostEqual(p.x, B.x) && _almostEqual(p.y, B.y)))
+ return false;
+
+
+ double cross = (p.y - A.y) * (B.x - A.x) - (p.x - A.x) * (B.y - A.y);
+
+ if(abs(cross) > TOL) return false;
+
+ double dot = (p.x - A.x) * (B.x - A.x) + (p.y - A.y)*(B.y - A.y);
+
+ if(dot < 0 || _almostEqual(dot, 0)) return false;
+
+ double len2 = (B.x - A.x)*(B.x - A.x) + (B.y - A.y)*(B.y - A.y);
+
+ if(dot > len2 || _almostEqual(dot, len2)) return false;
+
+ return true;
+ }
+
+ // return true if point is in the polygon, false if outside, and null if exactly on a point or edge
+ static int pointInPolygon(const Vector& point, const Cntr& polygon) {
+ if(polygon.size() < 3){
+ return 0;
+ }
+
+ bool inside = false;
+ Coord offsetx = polygon.offsetx;
+ Coord offsety = polygon.offsety;
+
+ for (size_t i = 0, j = polygon.size() - 1; i < polygon.size(); j=i++) {
+ auto xi = polygon[i].x + offsetx;
+ auto yi = polygon[i].y + offsety;
+ auto xj = polygon[j].x + offsetx;
+ auto yj = polygon[j].y + offsety;
+
+ if(_almostEqual(xi, point.x) && _almostEqual(yi, point.y)){
+ return 0; // no result
+ }
+
+ if(_onSegment({xi, yi}, {xj, yj}, point)){
+ return 0; // exactly on the segment
+ }
+
+ if(_almostEqual(xi, xj) && _almostEqual(yi, yj)){ // ignore very small lines
+ continue;
+ }
+
+ bool intersect = ((yi > point.y) != (yj > point.y)) &&
+ (point.x < (xj - xi) * (point.y - yi) / (yj - yi) + xi);
+ if (intersect) inside = !inside;
+ }
+
+ return inside? 1 : -1;
+ }
+
+ static bool intersect(const Cntr& A, const Cntr& B){
+ Contour a = A, b = B;
+ return shapelike::intersects(shapelike::create<S>(a), shapelike::create<S>(b));
+ }
+
+ static Vector _normalizeVector(const Vector& v) {
+ if(_almostEqual(v.x*v.x + v.y*v.y, Coord(1))){
+ return Point(v); // given vector was already a unit vector
+ }
+ auto len = sqrt(v.x*v.x + v.y*v.y);
+ auto inverse = 1/len;
+
+ return { Coord(v.x*inverse), Coord(v.y*inverse) };
+ }
+
+ static double pointDistance( const Vector& p,
+ const Vector& s1,
+ const Vector& s2,
+ Vector normal,
+ bool infinite = false)
+ {
+ normal = _normalizeVector(normal);
+
+ Vector dir = {
+ normal.y,
+ -normal.x
+ };
+
+ auto pdot = p.x*dir.x + p.y*dir.y;
+ auto s1dot = s1.x*dir.x + s1.y*dir.y;
+ auto s2dot = s2.x*dir.x + s2.y*dir.y;
+
+ auto pdotnorm = p.x*normal.x + p.y*normal.y;
+ auto s1dotnorm = s1.x*normal.x + s1.y*normal.y;
+ auto s2dotnorm = s2.x*normal.x + s2.y*normal.y;
+
+ if(!infinite){
+ if (((pdot<s1dot || _almostEqual(pdot, s1dot)) &&
+ (pdot<s2dot || _almostEqual(pdot, s2dot))) ||
+ ((pdot>s1dot || _almostEqual(pdot, s1dot)) &&
+ (pdot>s2dot || _almostEqual(pdot, s2dot))))
+ {
+ // dot doesn't collide with segment,
+ // or lies directly on the vertex
+ return dNAN;
+ }
+ if ((_almostEqual(pdot, s1dot) && _almostEqual(pdot, s2dot)) &&
+ (pdotnorm>s1dotnorm && pdotnorm>s2dotnorm))
+ {
+ return min(pdotnorm - s1dotnorm, pdotnorm - s2dotnorm);
+ }
+ if ((_almostEqual(pdot, s1dot) && _almostEqual(pdot, s2dot)) &&
+ (pdotnorm<s1dotnorm && pdotnorm<s2dotnorm)){
+ return -min(s1dotnorm-pdotnorm, s2dotnorm-pdotnorm);
+ }
+ }
+
+ return -(pdotnorm - s1dotnorm + (s1dotnorm - s2dotnorm)*(s1dot - pdot)
+ / double(s1dot - s2dot));
+ }
+
+ static double segmentDistance( const Vector& A,
+ const Vector& B,
+ const Vector& E,
+ const Vector& F,
+ Vector direction)
+ {
+ Vector normal = {
+ direction.y,
+ -direction.x
+ };
+
+ Vector reverse = {
+ -direction.x,
+ -direction.y
+ };
+
+ auto dotA = A.x*normal.x + A.y*normal.y;
+ auto dotB = B.x*normal.x + B.y*normal.y;
+ auto dotE = E.x*normal.x + E.y*normal.y;
+ auto dotF = F.x*normal.x + F.y*normal.y;
+
+ auto crossA = A.x*direction.x + A.y*direction.y;
+ auto crossB = B.x*direction.x + B.y*direction.y;
+ auto crossE = E.x*direction.x + E.y*direction.y;
+ auto crossF = F.x*direction.x + F.y*direction.y;
+
+// auto crossABmin = min(crossA, crossB);
+// auto crossABmax = max(crossA, crossB);
+
+// auto crossEFmax = max(crossE, crossF);
+// auto crossEFmin = min(crossE, crossF);
+
+ auto ABmin = min(dotA, dotB);
+ auto ABmax = max(dotA, dotB);
+
+ auto EFmax = max(dotE, dotF);
+ auto EFmin = min(dotE, dotF);
+
+ // segments that will merely touch at one point
+ if(_almostEqual(ABmax, EFmin, TOL) || _almostEqual(ABmin, EFmax,TOL)) {
+ return dNAN;
+ }
+ // segments miss eachother completely
+ if(ABmax < EFmin || ABmin > EFmax){
+ return dNAN;
+ }
+
+ double overlap = 0;
+
+ if((ABmax > EFmax && ABmin < EFmin) || (EFmax > ABmax && EFmin < ABmin))
+ {
+ overlap = 1;
+ }
+ else{
+ auto minMax = min(ABmax, EFmax);
+ auto maxMin = max(ABmin, EFmin);
+
+ auto maxMax = max(ABmax, EFmax);
+ auto minMin = min(ABmin, EFmin);
+
+ overlap = (minMax-maxMin)/(maxMax-minMin);
+ }
+
+ auto crossABE = (E.y - A.y) * (B.x - A.x) - (E.x - A.x) * (B.y - A.y);
+ auto crossABF = (F.y - A.y) * (B.x - A.x) - (F.x - A.x) * (B.y - A.y);
+
+ // lines are colinear
+ if(_almostEqual(crossABE,0) && _almostEqual(crossABF,0)){
+
+ Vector ABnorm = {B.y-A.y, A.x-B.x};
+ Vector EFnorm = {F.y-E.y, E.x-F.x};
+
+ auto ABnormlength = sqrt(ABnorm.x*ABnorm.x + ABnorm.y*ABnorm.y);
+ ABnorm.x /= ABnormlength;
+ ABnorm.y /= ABnormlength;
+
+ auto EFnormlength = sqrt(EFnorm.x*EFnorm.x + EFnorm.y*EFnorm.y);
+ EFnorm.x /= EFnormlength;
+ EFnorm.y /= EFnormlength;
+
+ // segment normals must point in opposite directions
+ if(abs(ABnorm.y * EFnorm.x - ABnorm.x * EFnorm.y) < TOL &&
+ ABnorm.y * EFnorm.y + ABnorm.x * EFnorm.x < 0){
+ // normal of AB segment must point in same direction as
+ // given direction vector
+ auto normdot = ABnorm.y * direction.y + ABnorm.x * direction.x;
+ // the segments merely slide along eachother
+ if(_almostEqual(normdot,0, TOL)){
+ return dNAN;
+ }
+ if(normdot < 0){
+ return 0.0;
+ }
+ }
+ return dNAN;
+ }
+
+ std::vector<double> distances; distances.reserve(10);
+
+ // coincident points
+ if(_almostEqual(dotA, dotE)){
+ distances.emplace_back(crossA-crossE);
+ }
+ else if(_almostEqual(dotA, dotF)){
+ distances.emplace_back(crossA-crossF);
+ }
+ else if(dotA > EFmin && dotA < EFmax){
+ auto d = pointDistance(A,E,F,reverse);
+ if(!isnan(d) && _almostEqual(d, 0))
+ { // A currently touches EF, but AB is moving away from EF
+ auto dB = pointDistance(B,E,F,reverse,true);
+ if(dB < 0 || _almostEqual(dB*overlap,0)){
+ d = dNAN;
+ }
+ }
+ if(!isnan(d)){
+ distances.emplace_back(d);
+ }
+ }
+
+ if(_almostEqual(dotB, dotE)){
+ distances.emplace_back(crossB-crossE);
+ }
+ else if(_almostEqual(dotB, dotF)){
+ distances.emplace_back(crossB-crossF);
+ }
+ else if(dotB > EFmin && dotB < EFmax){
+ auto d = pointDistance(B,E,F,reverse);
+
+ if(!isnan(d) && _almostEqual(d, 0))
+ { // crossA>crossB A currently touches EF, but AB is moving away from EF
+ double dA = pointDistance(A,E,F,reverse,true);
+ if(dA < 0 || _almostEqual(dA*overlap,0)){
+ d = dNAN;
+ }
+ }
+ if(!isnan(d)){
+ distances.emplace_back(d);
+ }
+ }
+
+ if(dotE > ABmin && dotE < ABmax){
+ auto d = pointDistance(E,A,B,direction);
+ if(!isnan(d) && _almostEqual(d, 0))
+ { // crossF<crossE A currently touches EF, but AB is moving away from EF
+ double dF = pointDistance(F,A,B,direction, true);
+ if(dF < 0 || _almostEqual(dF*overlap,0)){
+ d = dNAN;
+ }
+ }
+ if(!isnan(d)){
+ distances.emplace_back(d);
+ }
+ }
+
+ if(dotF > ABmin && dotF < ABmax){
+ auto d = pointDistance(F,A,B,direction);
+ if(!isnan(d) && _almostEqual(d, 0))
+ { // && crossE<crossF A currently touches EF,
+ // but AB is moving away from EF
+ double dE = pointDistance(E,A,B,direction, true);
+ if(dE < 0 || _almostEqual(dE*overlap,0)){
+ d = dNAN;
+ }
+ }
+ if(!isnan(d)){
+ distances.emplace_back(d);
+ }
+ }
+
+ if(distances.empty()){
+ return dNAN;
+ }
+
+ return *std::min_element(distances.begin(), distances.end());
+ }
+
+ static double polygonSlideDistance( const Cntr& AA,
+ const Cntr& BB,
+ Vector direction,
+ bool ignoreNegative)
+ {
+// Vector A1, A2, B1, B2;
+ Cntr A = AA;
+ Cntr B = BB;
+
+ Coord Aoffsetx = A.offsetx;
+ Coord Boffsetx = B.offsetx;
+ Coord Aoffsety = A.offsety;
+ Coord Boffsety = B.offsety;
+
+ // close the loop for polygons
+ if(A[0] != A[A.size()-1]){
+ A.emplace_back(AA[0]);
+ }
+
+ if(B[0] != B[B.size()-1]){
+ B.emplace_back(BB[0]);
+ }
+
+ auto& edgeA = A;
+ auto& edgeB = B;
+
+ double distance = dNAN, d = dNAN;
+
+ Vector dir = _normalizeVector(direction);
+
+// Vector normal = {
+// dir.y,
+// -dir.x
+// };
+
+// Vector reverse = {
+// -dir.x,
+// -dir.y,
+// };
+
+ for(size_t i = 0; i < edgeB.size() - 1; i++){
+ for(size_t j = 0; j < edgeA.size() - 1; j++){
+ Vector A1 = {x(edgeA[j]) + Aoffsetx, y(edgeA[j]) + Aoffsety };
+ Vector A2 = {x(edgeA[j+1]) + Aoffsetx, y(edgeA[j+1]) + Aoffsety};
+ Vector B1 = {x(edgeB[i]) + Boffsetx, y(edgeB[i]) + Boffsety };
+ Vector B2 = {x(edgeB[i+1]) + Boffsetx, y(edgeB[i+1]) + Boffsety};
+
+ if((_almostEqual(A1.x, A2.x) && _almostEqual(A1.y, A2.y)) ||
+ (_almostEqual(B1.x, B2.x) && _almostEqual(B1.y, B2.y))){
+ continue; // ignore extremely small lines
+ }
+
+ d = segmentDistance(A1, A2, B1, B2, dir);
+
+ if(!isnan(d) && (isnan(distance) || d < distance)){
+ if(!ignoreNegative || d > 0 || _almostEqual(d, 0)){
+ distance = d;
+ }
+ }
+ }
+ }
+ return distance;
+ }
+
+ static double polygonProjectionDistance(const Cntr& AA,
+ const Cntr& BB,
+ Vector direction)
+ {
+ Cntr A = AA;
+ Cntr B = BB;
+
+ auto Boffsetx = B.offsetx;
+ auto Boffsety = B.offsety;
+ auto Aoffsetx = A.offsetx;
+ auto Aoffsety = A.offsety;
+
+ // close the loop for polygons
+ if(A[0] != A[A.size()-1]){
+ A.push(A[0]);
+ }
+
+ if(B[0] != B[B.size()-1]){
+ B.push(B[0]);
+ }
+
+ auto& edgeA = A;
+ auto& edgeB = B;
+
+ double distance = dNAN, d;
+// Vector p, s1, s2;
+
+ for(size_t i = 0; i < edgeB.size(); i++) {
+ // the shortest/most negative projection of B onto A
+ double minprojection = dNAN;
+ Vector minp;
+ for(size_t j = 0; j < edgeA.size() - 1; j++){
+ Vector p = {x(edgeB[i]) + Boffsetx, y(edgeB[i]) + Boffsety };
+ Vector s1 = {x(edgeA[j]) + Aoffsetx, y(edgeA[j]) + Aoffsety };
+ Vector s2 = {x(edgeA[j+1]) + Aoffsetx, y(edgeA[j+1]) + Aoffsety };
+
+ if(abs((s2.y-s1.y) * direction.x -
+ (s2.x-s1.x) * direction.y) < TOL) continue;
+
+ // project point, ignore edge boundaries
+ d = pointDistance(p, s1, s2, direction);
+
+ if(!isnan(d) && (isnan(minprojection) || d < minprojection)) {
+ minprojection = d;
+ minp = p;
+ }
+ }
+
+ if(!isnan(minprojection) && (isnan(distance) ||
+ minprojection > distance)){
+ distance = minprojection;
+ }
+ }
+
+ return distance;
+ }
+
+ static std::pair<bool, Vector> searchStartPoint(
+ const Cntr& AA, const Cntr& BB, bool inside, const std::vector<Cntr>& NFP = {})
+ {
+ // clone arrays
+ auto A = AA;
+ auto B = BB;
+
+// // close the loop for polygons
+// if(A[0] != A[A.size()-1]){
+// A.push(A[0]);
+// }
+
+// if(B[0] != B[B.size()-1]){
+// B.push(B[0]);
+// }
+
+ // returns true if point already exists in the given nfp
+ auto inNfp = [](const Vector& p, const std::vector<Cntr>& nfp){
+ if(nfp.empty()){
+ return false;
+ }
+
+ for(size_t i=0; i < nfp.size(); i++){
+ for(size_t j = 0; j< nfp[i].size(); j++){
+ if(_almostEqual(p.x, nfp[i][j].x) &&
+ _almostEqual(p.y, nfp[i][j].y)){
+ return true;
+ }
+ }
+ }
+
+ return false;
+ };
+
+ for(size_t i = 0; i < A.size() - 1; i++){
+ if(!A[i].marked) {
+ A[i].marked = true;
+ for(size_t j = 0; j < B.size(); j++){
+ B.offsetx = A[i].x - B[j].x;
+ B.offsety = A[i].y - B[j].y;
+
+ int Binside = 0;
+ for(size_t k = 0; k < B.size(); k++){
+ int inpoly = pointInPolygon({B[k].x + B.offsetx, B[k].y + B.offsety}, A);
+ if(inpoly != 0){
+ Binside = inpoly;
+ break;
+ }
+ }
+
+ if(Binside == 0){ // A and B are the same
+ return {false, {}};
+ }
+
+ auto startPoint = std::make_pair(true, Vector(B.offsetx, B.offsety));
+ if(((Binside && inside) || (!Binside && !inside)) &&
+ !intersect(A,B) && !inNfp(startPoint.second, NFP)){
+ return startPoint;
+ }
+
+ // slide B along vector
+ auto vx = A[i+1].x - A[i].x;
+ auto vy = A[i+1].y - A[i].y;
+
+ double d1 = polygonProjectionDistance(A,B,{vx, vy});
+ double d2 = polygonProjectionDistance(B,A,{-vx, -vy});
+
+ double d = dNAN;
+
+ // todo: clean this up
+ if(isnan(d1) && isnan(d2)){
+ // nothin
+ }
+ else if(isnan(d1)){
+ d = d2;
+ }
+ else if(isnan(d2)){
+ d = d1;
+ }
+ else{
+ d = min(d1,d2);
+ }
+
+ // only slide until no longer negative
+ // todo: clean this up
+ if(!isnan(d) && !_almostEqual(d,0) && d > 0){
+
+ }
+ else{
+ continue;
+ }
+
+ auto vd2 = vx*vx + vy*vy;
+
+ if(d*d < vd2 && !_almostEqual(d*d, vd2)){
+ auto vd = sqrt(vx*vx + vy*vy);
+ vx *= d/vd;
+ vy *= d/vd;
+ }
+
+ B.offsetx += vx;
+ B.offsety += vy;
+
+ for(size_t k = 0; k < B.size(); k++){
+ int inpoly = pointInPolygon({B[k].x + B.offsetx, B[k].y + B.offsety}, A);
+ if(inpoly != 0){
+ Binside = inpoly;
+ break;
+ }
+ }
+ startPoint = std::make_pair(true, Vector{B.offsetx, B.offsety});
+ if(((Binside && inside) || (!Binside && !inside)) &&
+ !intersect(A,B) && !inNfp(startPoint.second, NFP)){
+ return startPoint;
+ }
+ }
+ }
+ }
+
+ return {false, Vector(0, 0)};
+ }
+
+ static std::vector<Cntr> noFitPolygon(Cntr A,
+ Cntr B,
+ bool inside,
+ bool searchEdges)
+ {
+ if(A.size() < 3 || B.size() < 3) {
+ throw GeometryException(GeomErr::NFP);
+ return {};
+ }
+
+ A.offsetx = 0;
+ A.offsety = 0;
+
+ long i = 0, j = 0;
+
+ auto minA = y(A[0]);
+ long minAindex = 0;
+
+ auto maxB = y(B[0]);
+ long maxBindex = 0;
+
+ for(i = 1; i < A.size(); i++){
+ A[i].marked = false;
+ if(y(A[i]) < minA){
+ minA = y(A[i]);
+ minAindex = i;
+ }
+ }
+
+ for(i = 1; i < B.size(); i++){
+ B[i].marked = false;
+ if(y(B[i]) > maxB){
+ maxB = y(B[i]);
+ maxBindex = i;
+ }
+ }
+
+ std::pair<bool, Vector> startpoint;
+
+ if(!inside){
+ // shift B such that the bottom-most point of B is at the top-most
+ // point of A. This guarantees an initial placement with no
+ // intersections
+ startpoint = { true,
+ { x(A[minAindex]) - x(B[maxBindex]),
+ y(A[minAindex]) - y(B[maxBindex]) }
+ };
+ }
+ else {
+ // no reliable heuristic for inside
+ startpoint = searchStartPoint(A, B, true);
+ }
+
+ std::vector<Cntr> NFPlist;
+
+ struct Touch {
+ int type;
+ long A;
+ long B;
+ Touch(int t, long a, long b): type(t), A(a), B(b) {}
+ };
+
+ while(startpoint.first) {
+
+ B.offsetx = startpoint.second.x;
+ B.offsety = startpoint.second.y;
+
+ // maintain a list of touching points/edges
+ std::vector<Touch> touching;
+
+ struct V {
+ Coord x, y;
+ Vector *start, *end;
+ operator bool() {
+ return start != nullptr && end != nullptr;
+ }
+ operator Vector() const { return {x, y}; }
+ } prevvector = {0, 0, nullptr, nullptr};
+
+ Cntr NFP;
+ NFP.emplace_back(x(B[0]) + B.offsetx, y(B[0]) + B.offsety);
+
+ auto referencex = x(B[0]) + B.offsetx;
+ auto referencey = y(B[0]) + B.offsety;
+ auto startx = referencex;
+ auto starty = referencey;
+ unsigned counter = 0;
+
+ // sanity check, prevent infinite loop
+ while(counter < 10*(A.size() + B.size())){
+ touching.clear();
+
+ // find touching vertices/edges
+ for(i = 0; i < A.size(); i++){
+ long nexti = (i == A.size() - 1) ? 0 : i + 1;
+ for(j = 0; j < B.size(); j++){
+
+ long nextj = (j == B.size() - 1) ? 0 : j + 1;
+
+ if( _almostEqual(A[i].x, B[j].x+B.offsetx) &&
+ _almostEqual(A[i].y, B[j].y+B.offsety))
+ {
+ touching.emplace_back(0, i, j);
+ }
+ else if( _onSegment(
+ A[i], A[nexti],
+ { B[j].x+B.offsetx, B[j].y + B.offsety}) )
+ {
+ touching.emplace_back(1, nexti, j);
+ }
+ else if( _onSegment(
+ {B[j].x+B.offsetx, B[j].y + B.offsety},
+ {B[nextj].x+B.offsetx, B[nextj].y + B.offsety},
+ A[i]) )
+ {
+ touching.emplace_back(2, i, nextj);
+ }
+ }
+ }
+
+ // generate translation vectors from touching vertices/edges
+ std::vector<V> vectors;
+ for(i=0; i < touching.size(); i++){
+ auto& vertexA = A[touching[i].A];
+ vertexA.marked = true;
+
+ // adjacent A vertices
+ auto prevAindex = touching[i].A - 1;
+ auto nextAindex = touching[i].A + 1;
+
+ prevAindex = (prevAindex < 0) ? A.size() - 1 : prevAindex; // loop
+ nextAindex = (nextAindex >= A.size()) ? 0 : nextAindex; // loop
+
+ auto& prevA = A[prevAindex];
+ auto& nextA = A[nextAindex];
+
+ // adjacent B vertices
+ auto& vertexB = B[touching[i].B];
+
+ auto prevBindex = touching[i].B-1;
+ auto nextBindex = touching[i].B+1;
+
+ prevBindex = (prevBindex < 0) ? B.size() - 1 : prevBindex; // loop
+ nextBindex = (nextBindex >= B.size()) ? 0 : nextBindex; // loop
+
+ auto& prevB = B[prevBindex];
+ auto& nextB = B[nextBindex];
+
+ if(touching[i].type == 0){
+
+ V vA1 = {
+ prevA.x - vertexA.x,
+ prevA.y - vertexA.y,
+ &vertexA,
+ &prevA
+ };
+
+ V vA2 = {
+ nextA.x - vertexA.x,
+ nextA.y - vertexA.y,
+ &vertexA,
+ &nextA
+ };
+
+ // B vectors need to be inverted
+ V vB1 = {
+ vertexB.x - prevB.x,
+ vertexB.y - prevB.y,
+ &prevB,
+ &vertexB
+ };
+
+ V vB2 = {
+ vertexB.x - nextB.x,
+ vertexB.y - nextB.y,
+ &nextB,
+ &vertexB
+ };
+
+ vectors.emplace_back(vA1);
+ vectors.emplace_back(vA2);
+ vectors.emplace_back(vB1);
+ vectors.emplace_back(vB2);
+ }
+ else if(touching[i].type == 1){
+ vectors.emplace_back(V{
+ vertexA.x-(vertexB.x+B.offsetx),
+ vertexA.y-(vertexB.y+B.offsety),
+ &prevA,
+ &vertexA
+ });
+
+ vectors.emplace_back(V{
+ prevA.x-(vertexB.x+B.offsetx),
+ prevA.y-(vertexB.y+B.offsety),
+ &vertexA,
+ &prevA
+ });
+ }
+ else if(touching[i].type == 2){
+ vectors.emplace_back(V{
+ vertexA.x-(vertexB.x+B.offsetx),
+ vertexA.y-(vertexB.y+B.offsety),
+ &prevB,
+ &vertexB
+ });
+
+ vectors.emplace_back(V{
+ vertexA.x-(prevB.x+B.offsetx),
+ vertexA.y-(prevB.y+B.offsety),
+ &vertexB,
+ &prevB
+ });
+ }
+ }
+
+ // TODO: there should be a faster way to reject vectors that
+ // will cause immediate intersection. For now just check them all
+
+ V translate = {0, 0, nullptr, nullptr};
+ double maxd = 0;
+
+ for(i = 0; i < vectors.size(); i++) {
+ if(vectors[i].x == 0 && vectors[i].y == 0){
+ continue;
+ }
+
+ // if this vector points us back to where we came from, ignore it.
+ // ie cross product = 0, dot product < 0
+ if(prevvector && vectors[i].y * prevvector.y + vectors[i].x * prevvector.x < 0){
+
+ // compare magnitude with unit vectors
+ double vectorlength = sqrt(vectors[i].x*vectors[i].x+vectors[i].y*vectors[i].y);
+ Vector unitv = {Coord(vectors[i].x/vectorlength),
+ Coord(vectors[i].y/vectorlength)};
+
+ double prevlength = sqrt(prevvector.x*prevvector.x+prevvector.y*prevvector.y);
+ Vector prevunit = { prevvector.x/prevlength, prevvector.y/prevlength};
+
+ // we need to scale down to unit vectors to normalize vector length. Could also just do a tan here
+ if(abs(unitv.y * prevunit.x - unitv.x * prevunit.y) < 0.0001){
+ continue;
+ }
+ }
+
+ V vi = vectors[i];
+ double d = polygonSlideDistance(A, B, vi, true);
+ double vecd2 = vectors[i].x*vectors[i].x + vectors[i].y*vectors[i].y;
+
+ if(isnan(d) || d*d > vecd2){
+ double vecd = sqrt(vectors[i].x*vectors[i].x + vectors[i].y*vectors[i].y);
+ d = vecd;
+ }
+
+ if(!isnan(d) && d > maxd){
+ maxd = d;
+ translate = vectors[i];
+ }
+ }
+
+ if(!translate || _almostEqual(maxd, 0))
+ {
+ // didn't close the loop, something went wrong here
+ NFP.clear();
+ break;
+ }
+
+ translate.start->marked = true;
+ translate.end->marked = true;
+
+ prevvector = translate;
+
+ // trim
+ double vlength2 = translate.x*translate.x + translate.y*translate.y;
+ if(maxd*maxd < vlength2 && !_almostEqual(maxd*maxd, vlength2)){
+ double scale = sqrt((maxd*maxd)/vlength2);
+ translate.x *= scale;
+ translate.y *= scale;
+ }
+
+ referencex += translate.x;
+ referencey += translate.y;
+
+ if(_almostEqual(referencex, startx) &&
+ _almostEqual(referencey, starty)) {
+ // we've made a full loop
+ break;
+ }
+
+ // if A and B start on a touching horizontal line,
+ // the end point may not be the start point
+ bool looped = false;
+ if(NFP.size() > 0) {
+ for(i = 0; i < NFP.size() - 1; i++) {
+ if(_almostEqual(referencex, NFP[i].x) &&
+ _almostEqual(referencey, NFP[i].y)){
+ looped = true;
+ }
+ }
+ }
+
+ if(looped){
+ // we've made a full loop
+ break;
+ }
+
+ NFP.emplace_back(referencex, referencey);
+
+ B.offsetx += translate.x;
+ B.offsety += translate.y;
+
+ counter++;
+ }
+
+ if(NFP.size() > 0){
+ NFPlist.emplace_back(NFP);
+ }
+
+ if(!searchEdges){
+ // only get outer NFP or first inner NFP
+ break;
+ }
+
+ startpoint =
+ searchStartPoint(A, B, inside, NFPlist);
+
+ }
+
+ return NFPlist;
+ }
+};
+
+template<class S> const double _alg<S>::TOL = std::pow(10, -9);
+
+}
+}
+
+#endif // NFP_SVGNEST_HPP