diff --git a/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp b/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
index cb0580ef4..b57b8dc53 100644
--- a/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
+++ b/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
@@ -251,460 +251,6 @@ inline NfpResult<RawShape> nfpConvexOnly(const RawShape& sh,
     return {rsh, top_nfp};
 }
 
-template<class RawShape>
-NfpResult<RawShape> nfpSimpleSimple(const RawShape& cstationary,
-                                           const RawShape& cother)
-{
-
-    // Algorithms are from the original algorithm proposed in paper:
-    // https://eprints.soton.ac.uk/36850/1/CORMSIS-05-05.pdf
-
-    // /////////////////////////////////////////////////////////////////////////
-    // Algorithm 1: Obtaining the minkowski sum
-    // /////////////////////////////////////////////////////////////////////////
-
-    // I guess this is not a full minkowski sum of the two input polygons by
-    // definition. This yields a subset that is compatible with the next 2
-    // algorithms.
-
-    using Result = NfpResult<RawShape>;
-    using Vertex = TPoint<RawShape>;
-    using Coord = TCoord<Vertex>;
-    using Edge = _Segment<Vertex>;
-    namespace sl = shapelike;
-    using std::signbit;
-    using std::sort;
-    using std::vector;
-    using std::ref;
-    using std::reference_wrapper;
-
-    // TODO The original algorithms expects the stationary polygon in
-    // counter clockwise and the orbiter in clockwise order.
-    // So for preventing any further complication, I will make the input
-    // the way it should be, than make my way around the orientations.
-
-    // Reverse the stationary contour to counter clockwise
-    auto stcont = sl::contour(cstationary);
-    {
-        std::reverse(sl::begin(stcont), sl::end(stcont));
-        stcont.pop_back();
-        auto it = std::min_element(sl::begin(stcont), sl::end(stcont),
-                               [](const Vertex& v1, const Vertex& v2) {
-            return getY(v1) < getY(v2);
-        });
-        std::rotate(sl::begin(stcont), it, sl::end(stcont));
-        sl::addVertex(stcont, sl::front(stcont));
-    }
-    RawShape stationary;
-    sl::contour(stationary) = stcont;
-
-    // Reverse the orbiter contour to counter clockwise
-    auto orbcont = sl::contour(cother);
-    {
-        std::reverse(orbcont.begin(), orbcont.end());
-
-        // Step 1: Make the orbiter reverse oriented
-
-        orbcont.pop_back();
-        auto it = std::min_element(orbcont.begin(), orbcont.end(),
-                              [](const Vertex& v1, const Vertex& v2) {
-            return getY(v1) < getY(v2);
-        });
-
-        std::rotate(orbcont.begin(), it, orbcont.end());
-        orbcont.emplace_back(orbcont.front());
-
-        for(auto &v : orbcont) v = -v;
-
-    }
-
-    // Copy the orbiter (contour only), we will have to work on it
-    RawShape orbiter;
-    sl::contour(orbiter) = orbcont;
-
-    // An edge with additional data for marking it
-    struct MarkedEdge {
-        Edge e; Radians turn_angle = 0; bool is_turning_point = false;
-        MarkedEdge() = default;
-        MarkedEdge(const Edge& ed, Radians ta, bool tp):
-            e(ed), turn_angle(ta), is_turning_point(tp) {}
-
-        // debug
-        std::string label;
-    };
-
-    // Container for marked edges
-    using EdgeList = vector<MarkedEdge>;
-
-    EdgeList A, B;
-
-    // This is how an edge list is created from the polygons
-    auto fillEdgeList = [](EdgeList& L, const RawShape& ppoly, int dir) {
-        auto& poly = sl::contour(ppoly);
-
-        L.reserve(sl::contourVertexCount(poly));
-
-        if(dir > 0) {
-            auto it = poly.begin();
-            auto nextit = std::next(it);
-
-            double turn_angle = 0;
-            bool is_turn_point = false;
-
-            while(nextit != poly.end()) {
-                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
-                it++; nextit++;
-            }
-        } else {
-            auto it = sl::rbegin(poly);
-            auto nextit = std::next(it);
-
-            double turn_angle = 0;
-            bool is_turn_point = false;
-
-            while(nextit != sl::rend(poly)) {
-                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
-                it++; nextit++;
-            }
-        }
-
-        auto getTurnAngle = [](const Edge& e1, const Edge& e2) {
-            auto phi = e1.angleToXaxis();
-            auto phi_prev = e2.angleToXaxis();
-            auto turn_angle = phi-phi_prev;
-            if(turn_angle > Pi) turn_angle -= TwoPi;
-            if(turn_angle < -Pi) turn_angle += TwoPi;
-            return turn_angle;
-        };
-
-        auto eit = L.begin();
-        auto enext = std::next(eit);
-
-        eit->turn_angle = getTurnAngle(L.front().e, L.back().e);
-
-        while(enext != L.end()) {
-            enext->turn_angle = getTurnAngle( enext->e, eit->e);
-            eit->is_turning_point =
-                    signbit(enext->turn_angle) != signbit(eit->turn_angle);
-            ++eit; ++enext;
-        }
-
-        L.back().is_turning_point = signbit(L.back().turn_angle) !=
-                                    signbit(L.front().turn_angle);
-
-    };
-
-    // Step 2: Fill the edgelists
-    fillEdgeList(A, stationary, 1);
-    fillEdgeList(B, orbiter, 1);
-
-    int i = 1;
-    for(MarkedEdge& me : A) {
-        std::cout << "a" << i << ":\n\t"
-                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
-                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
-                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
-                  << std::endl;
-
-        me.label = "a"; me.label += std::to_string(i);
-        i++;
-    }
-
-    i = 1;
-    for(MarkedEdge& me : B) {
-        std::cout << "b" << i << ":\n\t"
-                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
-                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
-                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
-                  << std::endl;
-        me.label = "b"; me.label += std::to_string(i);
-        i++;
-    }
-
-    // A reference to a marked edge that also knows its container
-    struct MarkedEdgeRef {
-        reference_wrapper<MarkedEdge> eref;
-        reference_wrapper<vector<MarkedEdgeRef>> container;
-        Coord dir = 1;  // Direction modifier
-
-        inline Radians angleX() const { return eref.get().e.angleToXaxis(); }
-        inline const Edge& edge() const { return eref.get().e; }
-        inline Edge& edge() { return eref.get().e; }
-        inline bool isTurningPoint() const {
-            return eref.get().is_turning_point;
-        }
-        inline bool isFrom(const vector<MarkedEdgeRef>& cont ) {
-            return &(container.get()) == &cont;
-        }
-        inline bool eq(const MarkedEdgeRef& mr) {
-            return &(eref.get()) == &(mr.eref.get());
-        }
-
-        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
-                      reference_wrapper<vector<MarkedEdgeRef>> ec):
-            eref(er), container(ec), dir(1) {}
-
-        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
-                      reference_wrapper<vector<MarkedEdgeRef>> ec,
-                      Coord d):
-            eref(er), container(ec), dir(d) {}
-    };
-
-    using EdgeRefList = vector<MarkedEdgeRef>;
-
-    // Comparing two marked edges
-    auto sortfn = [](const MarkedEdgeRef& e1, const MarkedEdgeRef& e2) {
-        return e1.angleX() < e2.angleX();
-    };
-
-    EdgeRefList Aref, Bref;     // We create containers for the references
-    Aref.reserve(A.size()); Bref.reserve(B.size());
-
-    // Fill reference container for the stationary polygon
-    std::for_each(A.begin(), A.end(), [&Aref](MarkedEdge& me) {
-        Aref.emplace_back( ref(me), ref(Aref) );
-    });
-
-    // Fill reference container for the orbiting polygon
-    std::for_each(B.begin(), B.end(), [&Bref](MarkedEdge& me) {
-        Bref.emplace_back( ref(me), ref(Bref) );
-    });
-
-    auto mink = [sortfn] // the Mink(Q, R, direction) sub-procedure
-            (const EdgeRefList& Q, const EdgeRefList& R, bool positive)
-    {
-
-        // Step 1 "merge sort_list(Q) and sort_list(R) to form merge_list(Q,R)"
-        // Sort the containers of edge references and merge them.
-        // Q could be sorted only once and be reused here but we would still
-        // need to merge it with sorted(R).
-
-        EdgeRefList merged;
-        EdgeRefList S, seq;
-        merged.reserve(Q.size() + R.size());
-
-        merged.insert(merged.end(), R.begin(), R.end());
-        std::stable_sort(merged.begin(), merged.end(), sortfn);
-        merged.insert(merged.end(), Q.begin(), Q.end());
-        std::stable_sort(merged.begin(), merged.end(), sortfn);
-
-        // Step 2 "set i = 1, k = 1, direction = 1, s1 = q1"
-        // we don't use i, instead, q is an iterator into Q. k would be an index
-        // into the merged sequence but we use "it" as an iterator for that
-
-        // here we obtain references for the containers for later comparisons
-        const auto& Rcont = R.begin()->container.get();
-        const auto& Qcont = Q.begin()->container.get();
-
-        // Set the initial direction
-        Coord dir = 1;
-
-        // roughly i = 1 (so q = Q.begin()) and s1 = q1 so S[0] = q;
-        if(positive) {
-            auto q = Q.begin();
-            S.emplace_back(*q);
-
-            // Roughly step 3
-
-            std::cout << "merged size: " << merged.size() << std::endl;
-            auto mit = merged.begin();
-            for(bool finish = false; !finish && q != Q.end();) {
-                ++q; // "Set i = i + 1"
-
-                while(!finish && mit != merged.end()) {
-                    if(mit->isFrom(Rcont)) {
-                        auto s = *mit;
-                        s.dir = dir;
-                        S.emplace_back(s);
-                    }
-
-                    if(mit->eq(*q)) {
-                        S.emplace_back(*q);
-                        if(mit->isTurningPoint()) dir = -dir;
-                        if(q == Q.begin()) finish = true;
-                        break;
-                    }
-
-                    mit += dir;
-    //                __nfp::advance(mit, merged, dir > 0);
-                }
-            }
-        } else {
-            auto q = Q.rbegin();
-            S.emplace_back(*q);
-
-            // Roughly step 3
-
-            std::cout << "merged size: " << merged.size() << std::endl;
-            auto mit = merged.begin();
-            for(bool finish = false; !finish && q != Q.rend();) {
-                ++q; // "Set i = i + 1"
-
-                while(!finish && mit != merged.end()) {
-                    if(mit->isFrom(Rcont)) {
-                        auto s = *mit;
-                        s.dir = dir;
-                        S.emplace_back(s);
-                    }
-
-                    if(mit->eq(*q)) {
-                        S.emplace_back(*q);
-                        S.back().dir = -1;
-                        if(mit->isTurningPoint()) dir = -dir;
-                        if(q == Q.rbegin()) finish = true;
-                        break;
-                    }
-
-                    mit += dir;
-            //                __nfp::advance(mit, merged, dir > 0);
-                }
-            }
-        }
-
-
-        // Step 4:
-
-        // "Let starting edge r1 be in position si in sequence"
-        // whaaat? I guess this means the following:
-        auto it = S.begin();
-        while(!it->eq(*R.begin())) ++it;
-
-        // "Set j = 1, next = 2, direction = 1, seq1 = si"
-        // we don't use j, seq is expanded dynamically.
-        dir = 1;
-        auto next = std::next(R.begin()); seq.emplace_back(*it);
-
-        // Step 5:
-        // "If all si edges have been allocated to seqj" should mean that
-        // we loop until seq has equal size with S
-        auto send = it; //it == S.begin() ? it : std::prev(it);
-        while(it != S.end()) {
-            ++it; if(it == S.end()) it = S.begin();
-            if(it == send) break;
-
-            if(it->isFrom(Qcont)) {
-                seq.emplace_back(*it); // "If si is from Q, j = j + 1, seqj = si"
-
-                // "If si is a turning point in Q,
-                // direction = - direction, next = next + direction"
-                if(it->isTurningPoint()) {
-                    dir = -dir;
-                    next += dir;
-//                    __nfp::advance(next, R, dir > 0);
-                }
-            }
-
-            if(it->eq(*next) /*&& dir == next->dir*/) { // "If si = direction.rnext"
-                // "j = j + 1, seqj = si, next = next + direction"
-                seq.emplace_back(*it);
-                next += dir;
-//                __nfp::advance(next, R, dir > 0);
-            }
-        }
-
-        return seq;
-    };
-
-    std::vector<EdgeRefList> seqlist;
-    seqlist.reserve(Bref.size());
-
-    EdgeRefList Bslope = Bref;  // copy Bref, we will make a slope diagram
-
-    // make the slope diagram of B
-    std::sort(Bslope.begin(), Bslope.end(), sortfn);
-
-    auto slopeit = Bslope.begin(); // search for the first turning point
-    while(!slopeit->isTurningPoint() && slopeit != Bslope.end()) slopeit++;
-
-    if(slopeit == Bslope.end()) {
-        // no turning point means convex polygon.
-        seqlist.emplace_back(mink(Aref, Bref, true));
-    } else {
-        int dir = 1;
-
-        auto firstturn = Bref.begin();
-        while(!firstturn->eq(*slopeit)) ++firstturn;
-
-        assert(firstturn != Bref.end());
-
-        EdgeRefList bgroup; bgroup.reserve(Bref.size());
-        bgroup.emplace_back(*slopeit);
-
-        auto b_it = std::next(firstturn);
-        while(b_it != firstturn) {
-            if(b_it == Bref.end()) b_it = Bref.begin();
-
-            while(!slopeit->eq(*b_it)) {
-                __nfp::advance(slopeit, Bslope, dir > 0);
-            }
-
-            if(!slopeit->isTurningPoint()) {
-                bgroup.emplace_back(*slopeit);
-            } else {
-                if(!bgroup.empty()) {
-                    if(dir > 0) bgroup.emplace_back(*slopeit);
-                    for(auto& me : bgroup) {
-                        std::cout << me.eref.get().label << ", ";
-                    }
-                    std::cout << std::endl;
-                    seqlist.emplace_back(mink(Aref, bgroup, dir == 1 ? true : false));
-                    bgroup.clear();
-                    if(dir < 0) bgroup.emplace_back(*slopeit);
-                } else {
-                    bgroup.emplace_back(*slopeit);
-                }
-
-                dir *= -1;
-            }
-            ++b_it;
-        }
-    }
-
-//    while(it != Bref.end()) // This is step 3 and step 4 in one loop
-//        if(it->isTurningPoint()) {
-//            R = {R.last, it++};
-//            auto seq = mink(Q, R, orientation);
-
-//            // TODO step 6 (should be 5 shouldn't it?): linking edges from A
-//            // I don't get this step
-
-//            seqlist.insert(seqlist.end(), seq.begin(), seq.end());
-//            orientation = !orientation;
-//        } else ++it;
-
-//    if(seqlist.empty()) seqlist = mink(Q, {Bref.begin(), Bref.end()}, true);
-
-    // /////////////////////////////////////////////////////////////////////////
-    // Algorithm 2: breaking Minkowski sums into track line trips
-    // /////////////////////////////////////////////////////////////////////////
-
-
-    // /////////////////////////////////////////////////////////////////////////
-    // Algorithm 3: finding the boundary of the NFP from track line trips
-    // /////////////////////////////////////////////////////////////////////////
-
-
-    for(auto& seq : seqlist) {
-        std::cout << "seqlist size: " << seq.size() << std::endl;
-        for(auto& s : seq) {
-            std::cout << (s.dir > 0 ? "" : "-") << s.eref.get().label << ", ";
-        }
-        std::cout << std::endl;
-    }
-
-    auto& seq = seqlist.front();
-    RawShape rsh;
-    Vertex top_nfp;
-    std::vector<Edge> edgelist; edgelist.reserve(seq.size());
-    for(auto& s : seq) {
-        edgelist.emplace_back(s.eref.get().e);
-    }
-
-    __nfp::buildPolygon(edgelist, rsh, top_nfp);
-
-    return Result(rsh, top_nfp);
-}
-
 // Specializable NFP implementation class. Specialize it if you have a faster
 // or better NFP implementation
 template<class RawShape, NfpLevel nfptype>
diff --git a/src/libnest2d/tests/test.cpp b/src/libnest2d/tests/test.cpp
index 1c4daf4af..3b0eae161 100644
--- a/src/libnest2d/tests/test.cpp
+++ b/src/libnest2d/tests/test.cpp
@@ -789,60 +789,6 @@ TEST(GeometryAlgorithms, nfpConvexConvex) {
 //    testNfp<NfpLevel::BOTH_CONCAVE, 1000>(nfp_concave_testdata);
 //}
 
-TEST(GeometryAlgorithms, nfpConcaveConcave) {
-    using namespace libnest2d;
-
-    Item stationary = {
-        {
-            {207, 76},
-            {194, 117},
-            {206, 117},
-            {206, 104},
-            {218, 104},
-            {231, 117},
-            {231, 130},
-            {244, 130},
-            {230, 92},
-            {220, 92},
-            {220, 84},
-            {239, 76},
-            {207, 76}
-        },
-        {}
-    };
-
-    Item orbiter = {
-        {
-            {78, 76},
-            {90, 89},
-            {76, 124},
-            {101, 124},
-            {101, 100},
-            {141, 113},
-            {141, 124},
-            {168, 124},
-            {158, 115},
-            {158, 104},
-            {121, 88},
-            {121, 76},
-            {78, 76}
-        },
-        {}
-    };
-
-    Rectangle r1(10, 10);
-    Rectangle r2(20, 20);
-    auto result = nfp::nfpSimpleSimple(stationary.transformedShape(),
-                                       orbiter.transformedShape());
-
-    svg::SVGWriter<PolygonImpl>::Config conf;
-    conf.mm_in_coord_units = 1;
-    svg::SVGWriter<PolygonImpl> wr(conf);
-    wr.writeItem(Item(result.first));
-    wr.save("simplesimple.svg");
-
-}
-
 TEST(GeometryAlgorithms, pointOnPolygonContour) {
     using namespace libnest2d;