Adding upstream version 1.2.2.upstream/1.2.2upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

24 files changed, 4426 insertions, 0 deletions

diff --git a/src/helper/CMakeLists.txt b/src/helper/CMakeLists.txt
new file mode 100644
index 0000000..18baf7b
--- /dev/null
+++ b/src/helper/CMakeLists.txt
@@ -0,0 +1,46 @@
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+include(${CMAKE_SOURCE_DIR}/CMakeScripts/UseGlibMarshal.cmake)
+
+GLIB_MARSHAL(sp_marshal sp-marshal "${CMAKE_CURRENT_BINARY_DIR}/helper")
+
+set(sp_marshal_SRC
+	${CMAKE_CURRENT_BINARY_DIR}/sp-marshal.cpp
+	${CMAKE_CURRENT_BINARY_DIR}/sp-marshal.h
+)
+
+set(helper_SRC
+	geom.cpp
+	geom-nodetype.cpp
+	geom-pathstroke.cpp
+	geom-pathvector_nodesatellites.cpp
+	geom-nodesatellite.cpp
+	gettext.cpp
+	pixbuf-ops.cpp
+	png-write.cpp
+	stock-items.cpp
+	#units-test.cpp
+
+	# we generate this file and it's .h counter-part
+	${sp_marshal_SRC}
+
+
+	# -------
+	# Headers
+	geom-curves.h
+	geom-nodetype.h
+	geom-pathstroke.h
+	geom-pathvector_nodesatellites.h
+	geom-nodesatellite.h
+	geom.h
+	gettext.h
+	mathfns.h
+	pixbuf-ops.h
+	png-write.h
+	stock-items.h
+)
+
+set_source_files_properties(sp_marshal_SRC PROPERTIES GENERATED true)
+
+# add_inkscape_lib(helper_LIB "${helper_SRC}")
+add_inkscape_source("${helper_SRC}")
diff --git a/src/helper/README b/src/helper/README
new file mode 100644
index 0000000..d9a6b90
--- /dev/null
+++ b/src/helper/README
@@ -0,0 +1,8 @@
+
+
+This directory contains a variety of helper code.
+
+To do:
+
+* Merge with 'util'.
+* Move individual files to more appropriate directories.
diff --git a/src/helper/auto-connection.h b/src/helper/auto-connection.h
new file mode 100644
index 0000000..368db34
--- /dev/null
+++ b/src/helper/auto-connection.h
@@ -0,0 +1,42 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+
+#ifndef SEEN_AUTO_CONNECTION_H
+#define SEEN_AUTO_CONNECTION_H
+
+#include <sigc++/connection.h>
+
+namespace Inkscape {
+
+// class to simplify re-subsribing to connections; automates disconnecting
+
+class auto_connection {
+public:
+    auto_connection(const sigc::connection& c): _connection(c) {}
+
+    auto_connection() = default;
+
+    ~auto_connection() {
+        _connection.disconnect();
+    }
+
+    auto_connection(const auto_connection&) = delete;
+    auto_connection& operator = (const auto_connection&) = delete;
+
+    // re-assign
+    auto_connection& operator = (const sigc::connection& c) {
+        _connection.disconnect();
+        _connection = c;
+        return *this;
+    }
+
+    void disconnect() {
+        *this = sigc::connection();
+    }
+
+private:
+    sigc::connection _connection;
+};
+
+} // namespace
+
+#endif
diff --git a/src/helper/geom-curves.h b/src/helper/geom-curves.h
new file mode 100644
index 0000000..08403e2
--- /dev/null
+++ b/src/helper/geom-curves.h
@@ -0,0 +1,55 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef INKSCAPE_HELPER_GEOM_CURVES_H
+#define INKSCAPE_HELPER_GEOM_CURVES_H
+
+/**
+ * @file
+ * Specific curve type functions for Inkscape, not provided by lib2geom.
+ */
+/*
+ * Author:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2008-2009 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/line.h>
+#include <2geom/bezier-curve.h>
+
+/// \todo un-inline this function
+inline bool is_straight_curve(Geom::Curve const & c)
+{
+    if( dynamic_cast<Geom::LineSegment const*>(&c) )
+    {
+        return true;
+    }
+    // the curve can be a quad/cubic bezier, but could still be a perfect straight line
+    // if the control points are exactly on the line connecting the initial and final points.
+    Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&c);
+    if (curve) {
+        Geom::Line line(curve->initialPoint(), curve->finalPoint());
+        std::vector<Geom::Point> pts = curve->controlPoints();
+        for (unsigned i = 1; i < pts.size() - 1; ++i) {
+            if (!are_near(pts[i], line))
+                return false;
+        }
+        return true;
+    }
+
+    return false;
+}
+
+#endif  // INKSCAPE_HELPER_GEOM_CURVES_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-nodesatellite.cpp b/src/helper/geom-nodesatellite.cpp
new file mode 100644
index 0000000..fec75ef
--- /dev/null
+++ b/src/helper/geom-nodesatellite.cpp
@@ -0,0 +1,244 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * \brief NodeSatellite a per node holder of data.
+ *//*
+ * Authors:
+ * see git history
+ * 2015 Jabier Arraiza Cenoz<jabier.arraiza@marker.es>
+ *
+ * Copyright (C) 2018 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/curve.h>
+#include <2geom/nearest-time.h>
+#include <2geom/path-intersection.h>
+#include <2geom/ray.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <helper/geom-nodesatellite.h>
+#include <optional>
+// log cache
+#ifdef _WIN32
+#include <Windows.h>
+#else
+#include <sys/time.h>
+#include <ctime>
+#endif
+
+///@brief NodeSatellite a per node holder of data.
+NodeSatellite::NodeSatellite() = default;
+
+NodeSatellite::NodeSatellite(NodeSatelliteType nodesatellite_type)
+    : nodesatellite_type(nodesatellite_type)
+    , is_time(false)
+    , selected(false)
+    , has_mirror(false)
+    , hidden(true)
+    , amount(0.0)
+    , angle(0.0)
+    , steps(0)
+{}
+
+NodeSatellite::~NodeSatellite() = default;
+
+///Calculate the time in curve_in with a size of A
+//TODO: find a better place to it
+double timeAtArcLength(double const A, Geom::Curve const &curve_in)
+{
+    if ( A == 0 || curve_in.isDegenerate()) {
+        return 0;
+    }
+
+    Geom::D2<Geom::SBasis> d2_in = curve_in.toSBasis();
+    double t = 0;
+    double length_part = curve_in.length();
+    if (A >= length_part || curve_in.isLineSegment()) {
+        if (length_part != 0) {
+            t = A / length_part;
+        }
+    } else if (!curve_in.isLineSegment()) {
+        std::vector<double> t_roots = roots(Geom::arcLengthSb(d2_in) - A);
+        if (!t_roots.empty()) {
+            t = t_roots[0];
+        }
+    }
+    return t;
+}
+
+///Calculate the size in curve_in with a point at A
+//TODO: find a better place to it
+double arcLengthAt(double const A, Geom::Curve const &curve_in)
+{
+    if ( A == 0 || curve_in.isDegenerate()) {
+        return 0;
+    }
+
+    double s = 0;
+    double length_part = curve_in.length();
+    if (A > length_part || curve_in.isLineSegment()) {
+        s = (A * length_part);
+    } else if (!curve_in.isLineSegment()) {
+        Geom::Curve *curve = curve_in.portion(0.0, A);
+        s = curve->length();
+        delete curve;
+    }
+    return s;
+}
+
+/// Convert a arc radius of a fillet/chamfer to his nodesatellite length -point position where fillet/chamfer knot be on
+/// original curve
+double NodeSatellite::radToLen(double const A, Geom::Curve const &curve_in, Geom::Curve const &curve_out) const
+{
+    double len = 0;
+    Geom::D2<Geom::SBasis> d2_in = curve_in.toSBasis();
+    Geom::D2<Geom::SBasis> d2_out = curve_out.toSBasis();
+    Geom::Piecewise<Geom::D2<Geom::SBasis> > offset_curve0 = 
+        Geom::Piecewise<Geom::D2<Geom::SBasis> >(d2_in) +
+        rot90(unitVector(derivative(d2_in))) * (A);
+    Geom::Piecewise<Geom::D2<Geom::SBasis> > offset_curve1 =
+        Geom::Piecewise<Geom::D2<Geom::SBasis> >(d2_out) +
+        rot90(unitVector(derivative(d2_out))) * (A);
+    offset_curve0[0][0].normalize();
+    offset_curve0[0][1].normalize();
+    Geom::Path p0 = path_from_piecewise(offset_curve0, 0.1)[0];
+    offset_curve1[0][0].normalize();
+    offset_curve1[0][1].normalize();
+    Geom::Path p1 = path_from_piecewise(offset_curve1, 0.1)[0];
+    Geom::Crossings cs = Geom::crossings(p0, p1);
+    if (cs.size() > 0) {
+        Geom::Point cp = p0(cs[0].ta);
+        double p0pt = nearest_time(cp, curve_out);
+        len = arcLengthAt(p0pt, curve_out);
+    } else {
+        if (A > 0) {
+            len = radToLen(A * -1, curve_in, curve_out);
+        }
+    }
+    return len;
+}
+
+/// Convert a nodesatellite length -point position where fillet/chamfer knot be on original curve- to a arc radius of
+/// fillet/chamfer
+double NodeSatellite::lenToRad(double const A, Geom::Curve const &curve_in, Geom::Curve const &curve_out,
+                               NodeSatellite const previousNodeSatellite) const
+{
+    double time_in = (previousNodeSatellite).time(A, true, curve_in);
+    double time_out = timeAtArcLength(A, curve_out);
+    Geom::Point start_arc_point = curve_in.pointAt(time_in);
+    Geom::Point end_arc_point = curve_out.pointAt(time_out);
+    Geom::Curve *knot_curve1 = curve_in.portion(0, time_in);
+    Geom::Curve *knot_curve2 = curve_out.portion(time_out, 1);
+    Geom::CubicBezier const *cubic1 = dynamic_cast<Geom::CubicBezier const *>(&*knot_curve1);
+    Geom::Ray ray1(start_arc_point, curve_in.pointAt(1));
+    if (cubic1) {
+        ray1.setPoints((*cubic1)[2], start_arc_point);
+    }
+    Geom::CubicBezier const *cubic2 = dynamic_cast<Geom::CubicBezier const *>(&*knot_curve2);
+    Geom::Ray ray2(curve_out.pointAt(0), end_arc_point);
+    if (cubic2) {
+        ray2.setPoints(end_arc_point, (*cubic2)[1]);
+    }
+    bool ccw_toggle = cross(curve_in.pointAt(1) - start_arc_point,
+                           end_arc_point - start_arc_point) < 0;
+    double distance_arc =
+        Geom::distance(start_arc_point, middle_point(start_arc_point, end_arc_point));
+    double angle = angle_between(ray1, ray2, ccw_toggle);
+    double divisor = std::sin(angle / 2.0);
+    if (divisor > 0) {
+        return distance_arc / divisor;
+    }
+    return 0;
+}
+
+/// Get the time position of the nodesatellite in curve_in
+double NodeSatellite::time(Geom::Curve const &curve_in, bool inverse) const
+{
+    double t = amount;
+    if (!is_time) {
+        t = time(t, inverse, curve_in);
+    } else if (inverse) {
+        t = 1-t;
+    }
+    if (t > 1) {
+        t = 1;
+    }
+    return t;
+}
+
+///Get the time from a length A in other curve, a boolean inverse given to reverse time
+double NodeSatellite::time(double A, bool inverse, Geom::Curve const &curve_in) const
+{
+    if (A == 0 && inverse) {
+        return 1;
+    }
+    if (A == 0 && !inverse) {
+        return 0;
+    }
+    if (!inverse) {
+        return timeAtArcLength(A, curve_in);
+    }
+    double length_part = curve_in.length();
+    A = length_part - A;
+    return timeAtArcLength(A, curve_in);
+}
+
+/// Get the length of the nodesatellite in curve_in
+double NodeSatellite::arcDistance(Geom::Curve const &curve_in) const
+{
+    double s = amount;
+    if (is_time) {
+        s = arcLengthAt(s, curve_in);
+    }
+    return s;
+}
+
+/// Get the point position of the nodesatellite
+Geom::Point NodeSatellite::getPosition(Geom::Curve const &curve_in, bool inverse) const
+{
+    double t = time(curve_in, inverse);
+    return curve_in.pointAt(t);
+}
+
+/// Set the position of the nodesatellite from a given point P
+void NodeSatellite::setPosition(Geom::Point const p, Geom::Curve const &curve_in, bool inverse)
+{
+    Geom::Curve * curve = const_cast<Geom::Curve *>(&curve_in);
+    if (inverse) {
+        curve = curve->reverse();
+    }
+    double A = Geom::nearest_time(p, *curve);
+    if (!is_time) {
+        A = arcLengthAt(A, *curve);
+    }
+    amount = A;
+}
+
+/// Map a nodesatellite type with gchar
+void NodeSatellite::setNodeSatellitesType(gchar const *A)
+{
+    std::map<std::string, NodeSatelliteType> gchar_map_to_nodesatellite_type = boost::assign::map_list_of("F", FILLET)(
+        "IF", INVERSE_FILLET)("C", CHAMFER)("IC", INVERSE_CHAMFER)("KO", INVALID_SATELLITE);
+    std::map<std::string, NodeSatelliteType>::iterator it = gchar_map_to_nodesatellite_type.find(std::string(A));
+    if (it != gchar_map_to_nodesatellite_type.end()) {
+        nodesatellite_type = it->second;
+    }
+}
+
+/// Map a gchar with nodesatelliteType
+gchar const *NodeSatellite::getNodeSatellitesTypeGchar() const
+{
+    std::map<NodeSatelliteType, gchar const *> nodesatellite_type_to_gchar_map = boost::assign::map_list_of(
+        FILLET, "F")(INVERSE_FILLET, "IF")(CHAMFER, "C")(INVERSE_CHAMFER, "IC")(INVALID_SATELLITE, "KO");
+    return nodesatellite_type_to_gchar_map.at(nodesatellite_type);
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-nodesatellite.h b/src/helper/geom-nodesatellite.h
new file mode 100644
index 0000000..1bb6c54
--- /dev/null
+++ b/src/helper/geom-nodesatellite.h
@@ -0,0 +1,109 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * NodeSatellite -- a per node holder of data.
+ *//*
+ * Authors:
+ * see git history
+ * Jabier Arraiza Cenoz<jabier.arraiza@marker.es>
+ *
+ * Copyright (C) 2017 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifndef SEEN_SATELLITE_H
+#define SEEN_SATELLITE_H
+
+#include <map>
+#include <boost/assign.hpp>
+#include <2geom/sbasis-geometric.h>
+#include "util/enums.h"
+
+enum NodeSatelliteType
+{
+    FILLET = 0,       // Fillet
+    INVERSE_FILLET,   // Inverse Fillet
+    CHAMFER,          // Chamfer
+    INVERSE_CHAMFER,  // Inverse Chamfer
+    INVALID_SATELLITE // Invalid NodeSatellite
+};
+/**
+ * @brief NodeSatellite a per node holder of data.
+ */
+
+class NodeSatellite
+{
+public:
+    NodeSatellite();
+    NodeSatellite(NodeSatelliteType nodesatellite_type);
+
+    virtual ~NodeSatellite();
+    void setIsTime(bool set_is_time)
+    {
+        is_time = set_is_time;
+    }
+    void setSelected(bool set_selected)
+    {
+        selected = set_selected;
+    }
+    void setHasMirror(bool set_has_mirror)
+    {
+        has_mirror = set_has_mirror;
+    }
+    void setHidden(bool set_hidden)
+    {
+        hidden = set_hidden;
+    }
+    void setAmount(double set_amount)
+    {
+        amount = set_amount;
+    }
+    void setAngle(double set_angle)
+    {
+        angle = set_angle;
+    }
+    void setSteps(size_t set_steps)
+    {
+        steps = set_steps;
+    }
+    double lenToRad(double const A, Geom::Curve const &curve_in, Geom::Curve const &curve_out,
+                    NodeSatellite const previousNodeSatellite) const;
+    double radToLen(double const A, Geom::Curve const &curve_in,
+                    Geom::Curve const &curve_out) const;
+
+    double time(Geom::Curve const &curve_in, bool inverse = false) const;
+    double time(double A, bool inverse, Geom::Curve const &curve_in) const;
+    double arcDistance(Geom::Curve const &curve_in) const;
+
+    void setPosition(Geom::Point const p, Geom::Curve const &curve_in, bool inverse = false);
+    Geom::Point getPosition(Geom::Curve const &curve_in, bool inverse = false) const;
+
+    void setNodeSatellitesType(gchar const *A);
+    gchar const *getNodeSatellitesTypeGchar() const;
+    NodeSatelliteType nodesatellite_type;
+    // The value stored could be a time value of the nodesatellite in the curve or a length of distance to the node from
+    // the nodesatellite "is_time" tells us if it's a time or length value
+    bool is_time;
+    bool selected;
+    bool has_mirror;
+    bool hidden;
+    // in "amount" we store the time or distance used in the nodesatellite
+    double amount;
+    double angle;
+    size_t steps;
+};
+
+double timeAtArcLength(double const A, Geom::Curve const &curve_in);
+double arcLengthAt(double const A, Geom::Curve const &curve_in);
+
+#endif // SEEN_SATELLITE_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-nodetype.cpp b/src/helper/geom-nodetype.cpp
new file mode 100644
index 0000000..e04b08d
--- /dev/null
+++ b/src/helper/geom-nodetype.cpp
@@ -0,0 +1,59 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * Specific nodetype geometry functions for Inkscape, not provided my lib2geom.
+ *
+ * Author:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2008 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include "helper/geom-nodetype.h"
+
+#include <2geom/curve.h>
+
+namespace Geom {
+
+/*
+ * NOTE: THIS METHOD NEVER RETURNS "NODE_SYMM".
+ * Returns the nodetype between c_incoming and c_outgoing. Location of the node is
+ * at c_incoming.pointAt(1) == c_outgoing.pointAt(0). If these two are unequal, 
+ * the returned type is NODE_NONE.
+ * Comparison is based on the unitTangent, does not work for NODE_SYMM!
+ */
+NodeType get_nodetype(Curve const &c_incoming, Curve const &c_outgoing)
+{
+    if ( !are_near(c_incoming.pointAt(1), c_outgoing.pointAt(0)) )
+        return NODE_NONE;
+
+    Geom::Curve *crv = c_incoming.reverse();
+    Geom::Point deriv_1 = -crv->unitTangentAt(0);
+    delete crv;
+    Geom::Point deriv_2 = c_outgoing.unitTangentAt(0);
+    double this_angle_L2 = Geom::L2(deriv_1);
+    double next_angle_L2 = Geom::L2(deriv_2);
+    double both_angles_L2 = Geom::L2(deriv_1 + deriv_2);
+    if ( (this_angle_L2 > 1e-6) &&
+         (next_angle_L2 > 1e-6) &&
+         ((this_angle_L2 + next_angle_L2 - both_angles_L2) < 1e-3) )
+    {
+        return NODE_SMOOTH;
+    }
+
+    return NODE_CUSP;
+}
+
+} // end namespace Geom
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-nodetype.h b/src/helper/geom-nodetype.h
new file mode 100644
index 0000000..d7d4d4c
--- /dev/null
+++ b/src/helper/geom-nodetype.h
@@ -0,0 +1,57 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef INKSCAPE_HELPER_GEOM_NODETYPE_H
+#define INKSCAPE_HELPER_GEOM_NODETYPE_H
+
+/**
+ * @file
+ * Specific nodetype geometry functions for Inkscape, not provided my lib2geom.
+ */
+/*
+ * Author:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2008 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/forward.h>
+
+namespace Geom {
+
+/**
+ *  What kind of node is this?  This is the value for the node->type
+ *  field.  NodeType indicates the degree of continuity required for
+ *  the node.  I think that the corresponding integer indicates which
+ *  derivate is connected. (Thus 2 means that the node is continuous
+ *  to the second derivative, i.e. has matching endpoints and tangents)
+ */
+enum NodeType {
+/**  Discontinuous node, usually either start or endpoint of a path */
+    NODE_NONE,
+/**  This node continuously joins two segments, but the unit tangent is discontinuous.*/
+    NODE_CUSP,
+/**  This node continuously joins two segments, with continuous *unit* tangent. */
+    NODE_SMOOTH,
+/**  This node is symmetric. I.e. continuously joins two segments with continuous derivative */
+    NODE_SYMM
+};
+
+
+NodeType get_nodetype(Curve const &c_incoming, Curve const &c_outgoing);
+
+
+} // namespace Geom
+
+#endif  // INKSCAPE_HELPER_GEOM_NODETYPE_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-pathstroke.cpp b/src/helper/geom-pathstroke.cpp
new file mode 100644
index 0000000..228ea5c
--- /dev/null
+++ b/src/helper/geom-pathstroke.cpp
@@ -0,0 +1,1158 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/* Authors:
+ *   Liam P. White
+ *   Tavmjong Bah
+ *   Alexander Brock
+ *
+ * Copyright (C) 2014-2015 Authors
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <iomanip>
+#include <2geom/path-sink.h>
+#include <2geom/sbasis-to-bezier.h> // cubicbezierpath_from_sbasis
+#include <2geom/path-intersection.h>
+#include <2geom/circle.h>
+
+#include "helper/geom-pathstroke.h"
+
+namespace Geom {
+
+static Point intersection_point(Point origin_a, Point vector_a, Point origin_b, Point vector_b)
+{
+    Coord denom = cross(vector_a, vector_b);
+    if (!are_near(denom,0.)) {
+        Coord t = (cross(vector_b, origin_a) + cross(origin_b, vector_b)) / denom;
+        return origin_a + vector_a*t;
+    }
+    return Point(infinity(), infinity());
+}
+
+/**
+* Find circle that touches inside of the curve, with radius matching the curvature, at time value \c t.
+* Because this method internally uses unitTangentAt, t should be smaller than 1.0 (see unitTangentAt).
+*/
+static Circle touching_circle( D2<SBasis> const &curve, double t, double tol=0.01 )
+{
+    D2<SBasis> dM=derivative(curve);
+    if ( are_near(L2sq(dM(t)), tol) ) {
+        dM=derivative(dM);
+    }
+    if ( are_near(L2sq(dM(t)), tol) ) {   // try second time
+        dM=derivative(dM);
+    }
+    Piecewise<D2<SBasis> > unitv = unitVector(dM,tol);
+    Piecewise<SBasis> dMlength = dot(Piecewise<D2<SBasis> >(dM),unitv);
+    Piecewise<SBasis> k = cross(derivative(unitv),unitv);
+    k = divide(k,dMlength,tol,3);
+    double curv = k(t); // note that this value is signed
+
+    Geom::Point normal = unitTangentAt(curve, t).cw();
+    double radius = 1/curv;
+    Geom::Point center = curve(t) + radius*normal;
+    return Geom::Circle(center, fabs(radius));
+}
+
+
+// Area of triangle given three corner points
+static double area( Geom::Point a, Geom::Point b, Geom::Point c ) {
+
+    using Geom::X;
+    using Geom::Y;
+    return( 0.5 * fabs( ( a[X]*(b[Y]-c[Y]) + b[X]*(c[Y]-a[Y]) + c[X]*(a[Y]-b[Y]) ) ) );
+}
+
+// Alternative touching circle routine directly using Beziers. Works only at end points.
+static Circle touching_circle( CubicBezier const &curve, bool start ) {
+
+    double k = 0;
+    Geom::Point p;
+    Geom::Point normal;
+    if ( start ) {
+        double distance = Geom::distance( curve[1], curve[0] );
+        k = 4.0/3.0 * area( curve[0], curve[1], curve[2] ) /
+            (distance * distance * distance);
+        if( Geom::cross(curve[0]-curve[1], curve[1]-curve[2]) < 0 ) {
+            k = -k;
+        }
+        p = curve[0];
+        normal = Geom::Point(curve[1] - curve[0]).cw();
+        normal.normalize();
+        // std::cout << "Start k: " << k << " d: " << distance << " normal: " << normal << std::endl;
+    } else {
+        double distance = Geom::distance( curve[3], curve[2] );
+        k = 4.0/3.0 * area( curve[1], curve[2], curve[3] ) /
+            (distance * distance * distance);
+        if( Geom::cross(curve[1]-curve[2], curve[2]-curve[3]) < 0 ) {
+            k = -k;
+        }
+        p = curve[3];
+        normal = Geom::Point(curve[3] - curve[2]).cw();
+        normal.normalize();
+        // std::cout << "End   k: " << k << " d: " << distance << " normal: " << normal << std::endl;
+    }
+    if( k == 0 ) {
+        return Geom::Circle( Geom::Point(0,std::numeric_limits<float>::infinity()),
+                             std::numeric_limits<float>::infinity());
+    } else {
+        double radius = 1/k;
+        Geom::Point center = p + normal * radius;
+        return Geom::Circle( center, fabs(radius) );
+    }
+}
+}
+
+namespace {
+
+// Internal data structure
+
+struct join_data {
+    join_data(Geom::Path &_res, Geom::Path const&_outgoing, Geom::Point _in_tang, Geom::Point _out_tang, double _miter, double _width)
+        : res(_res), outgoing(_outgoing), in_tang(_in_tang), out_tang(_out_tang), miter(_miter), width(_width) {};
+
+    // contains the current path that is being built on
+    Geom::Path &res;
+
+    // contains the next curve to append
+    Geom::Path const& outgoing;
+
+    // input tangents
+    Geom::Point in_tang;
+    Geom::Point out_tang;
+
+    // line parameters
+    double miter;
+    double width; // half stroke width
+};
+
+// Join functions must append the outgoing path
+
+typedef void join_func(join_data jd);
+
+void bevel_join(join_data jd)
+{
+    jd.res.appendNew<Geom::LineSegment>(jd.outgoing.initialPoint());
+    jd.res.append(jd.outgoing);
+}
+
+void round_join(join_data jd)
+{
+    jd.res.appendNew<Geom::EllipticalArc>(jd.width, jd.width, 0, false, jd.width <= 0, jd.outgoing.initialPoint());
+    jd.res.append(jd.outgoing);
+}
+
+void miter_join_internal(join_data const &jd, bool clip)
+{
+    using namespace Geom;
+
+    Curve const& incoming = jd.res.back();
+    Curve const& outgoing = jd.outgoing.front();
+    Path &res = jd.res;
+    double width = jd.width, miter = jd.miter;
+
+    Point tang1 = jd.in_tang;
+    Point tang2 = jd.out_tang;
+    Point p = intersection_point(incoming.finalPoint(), tang1, outgoing.initialPoint(), tang2);
+
+    bool satisfied = false;
+    bool inc_ls = res.back_open().degreesOfFreedom() <= 4;
+
+    if (p.isFinite()) {
+        // check size of miter
+        Point point_on_path = incoming.finalPoint() + rot90(tang1)*width;
+        // SVG defines miter length as distance between inner intersection and outer intersection,
+        // which is twice the distance from p to point_on_path but width is half stroke width.
+        satisfied = distance(p, point_on_path) <= miter * width; 
+        if (satisfied) {
+            // miter OK, check to see if we can do a relocation
+            if (inc_ls) {
+                res.setFinal(p);
+            } else {
+                res.appendNew<LineSegment>(p);
+            }
+        } else if (clip) {
+            // std::cout << "  Clipping ------------ " << std::endl;
+            // miter needs clipping, find two points
+            Point bisector_versor = Line(point_on_path, p).versor();
+            Point point_limit = point_on_path + miter * width * bisector_versor;
+            // std::cout << "     bisector_versor: " << bisector_versor << std::endl;
+            // std::cout << "     point_limit: " << point_limit << std::endl;
+            Point p1 = intersection_point(incoming.finalPoint(), tang1, point_limit, bisector_versor.cw());
+            Point p2 = intersection_point(outgoing.initialPoint(), tang2, point_limit, bisector_versor.cw());
+            // std::cout << "     p1: " << p1 << std::endl;
+            // std::cout << "     p2: " << p2 << std::endl;
+            if (inc_ls) {
+                res.setFinal(p1);
+            } else {
+                res.appendNew<LineSegment>(p1);
+            }
+            res.appendNew<LineSegment>(p2);
+        }
+    }
+
+    res.appendNew<LineSegment>(outgoing.initialPoint());
+
+    // check if we can do another relocation
+    bool out_ls = outgoing.degreesOfFreedom() <= 4;
+
+    if ((satisfied || clip) && out_ls) {
+        res.setFinal(outgoing.finalPoint());
+    } else {
+        res.append(outgoing);
+    }
+
+    // either way, add the rest of the path
+    res.insert(res.end(), ++jd.outgoing.begin(), jd.outgoing.end());
+}
+
+void miter_join(join_data jd) { miter_join_internal(jd, false); }
+void miter_clip_join(join_data jd) { miter_join_internal(jd, true); }
+
+Geom::Point pick_solution(std::vector<Geom::ShapeIntersection> points, Geom::Point tang2, Geom::Point endPt)
+{
+    assert(points.size() == 2);
+    Geom::Point sol;
+    if ( dot(tang2, points[0].point() - endPt) > 0 ) {
+        // points[0] is bad, choose points[1]
+        sol = points[1];
+    } else if ( dot(tang2, points[1].point() - endPt) > 0 ) { // points[0] could be good, now check points[1]
+        // points[1] is bad, choose points[0]
+        sol = points[0];
+    } else {
+        // both points are good, choose nearest
+        sol = ( distanceSq(endPt, points[0].point()) < distanceSq(endPt, points[1].point()) )
+            ? points[0].point() : points[1].point();
+    }
+    return sol;
+}
+
+// Arcs line join. If two circles don't intersect, expand inner circle.
+Geom::Point expand_circle( Geom::Circle &inner_circle, Geom::Circle const &outer_circle, Geom::Point const &start_pt, Geom::Point const &start_tangent ) {
+    // std::cout << "expand_circle:" << std::endl;
+    // std::cout << "  outer_circle: radius: " << outer_circle.radius() << "  center: " << outer_circle.center() << std::endl;
+    // std::cout << "  start: point: " << start_pt << "  tangent: " << start_tangent << std::endl;
+
+    if( !(outer_circle.contains(start_pt) ) ) {
+        // std::cout << "  WARNING: Outer circle does not contain starting point!" << std::endl;
+        return Geom::Point(0,0);
+    }
+
+    Geom::Line secant1(start_pt, start_pt + start_tangent);
+    std::vector<Geom::ShapeIntersection> chord1_pts = outer_circle.intersect(secant1);
+    // std::cout << "  chord1: " << chord1_pts[0].point() << ", " << chord1_pts[1].point() << std::endl;
+    Geom::LineSegment chord1(chord1_pts[0].point(), chord1_pts[1].point());
+
+    Geom::Line bisector = make_bisector_line( chord1 );
+    std::vector<Geom::ShapeIntersection> chord2_pts = outer_circle.intersect(bisector);
+    // std::cout << "  chord2: " << chord2_pts[0].point() << ", " << chord2_pts[1].point() << std::endl;
+    Geom::LineSegment chord2(chord2_pts[0].point(), chord2_pts[1].point());
+
+    // Find D, point on chord2 and on circle closest to start point
+    Geom::Coord d0 = Geom::distance(chord2_pts[0].point(),start_pt);
+    Geom::Coord d1 = Geom::distance(chord2_pts[1].point(),start_pt);
+    // std::cout << "  d0: " << d0 << " d1: " << d1 << std::endl;
+    Geom::Point d = (d0 < d1) ? chord2_pts[0].point() : chord2_pts[1].point();
+    Geom::Line da(d,start_pt);
+
+    // Chord through start point and point D
+    std::vector<Geom::ShapeIntersection> chord3_pts =  outer_circle.intersect(da);
+    // std::cout << "  chord3: " << chord3_pts[0].point() << ", " << chord3_pts[1].point() << std::endl;
+
+    // Find farthest point on chord3 and on circle (could be more robust)
+    Geom::Coord d2 = Geom::distance(chord3_pts[0].point(),d);
+    Geom::Coord d3 = Geom::distance(chord3_pts[1].point(),d);
+    // std::cout << "  d2: " << d2 << " d3: " << d3 << std::endl;
+
+    // Find point P, the intersection of outer circle and new inner circle
+    Geom::Point p = (d2 > d3) ? chord3_pts[0].point() : chord3_pts[1].point();
+
+    // Find center of new circle: it is at the intersection of the bisector
+    // of the chord defined by the start point and point P and a line through
+    // the start point and parallel to the first bisector.
+    Geom::LineSegment chord4(start_pt,p);
+    Geom::Line bisector2 = make_bisector_line( chord4 );
+    Geom::Line diameter = make_parallel_line( start_pt, bisector );
+    std::vector<Geom::ShapeIntersection> center_new = bisector2.intersect( diameter );
+    // std::cout << "  center_new: " << center_new[0].point() << std::endl;
+    Geom::Coord r_new = Geom::distance( center_new[0].point(), start_pt );
+    // std::cout << "  r_new: " << r_new << std::endl;
+
+    inner_circle.setCenter( center_new[0].point() );
+    inner_circle.setRadius( r_new );
+    return p;
+}
+
+// Arcs line join. If two circles don't intersect, adjust both circles so they just touch.
+// Increase (decrease) the radius of circle 1 and decrease (increase) of circle 2 by the same amount keeping the given points and tangents fixed.
+Geom::Point adjust_circles( Geom::Circle &circle1, Geom::Circle &circle2, Geom::Point const &point1, Geom::Point const &point2, Geom::Point const &tan1, Geom::Point const &tan2 ) {
+
+    Geom::Point n1 = (circle1.center() - point1).normalized(); // Always points towards center
+    Geom::Point n2 = (circle2.center() - point2).normalized();
+    Geom::Point sum_n = n1 + n2;
+
+    double r1 = circle1.radius();
+    double r2 = circle2.radius();
+    double delta_r = r2 - r1;
+    Geom::Point c1 = circle1.center();
+    Geom::Point c2 = circle2.center();
+    Geom::Point delta_c = c2 - c1;
+
+    // std::cout << "adjust_circles:" << std::endl;
+    // std::cout << "    norm: " << n1 << "; " << n2 << std::endl;
+    // std::cout << "    sum_n: " << sum_n << std::endl;
+    // std::cout << "    delta_r: " << delta_r << std::endl;
+    // std::cout << "    delta_c: " << delta_c << std::endl;
+
+    // Quadratic equation
+    double A = 4 - sum_n.length() * sum_n.length();
+    double B = 4.0 * delta_r - 2.0 * Geom::dot( delta_c, sum_n );
+    double C = delta_r * delta_r - delta_c.length() * delta_c.length();
+
+    double s1 = 0;
+    double s2 = 0;
+
+    if( fabs(A) < 0.01 ) {
+        // std::cout << "    A near zero! $$$$$$$$$$$$$$$$$$" << std::endl;
+        if( B != 0 ) {
+            s1 = -C/B;
+            s2 = -s1;
+        }
+    } else {
+        s1 = (-B + sqrt(B*B - 4*A*C))/(2*A);
+        s2 = (-B - sqrt(B*B - 4*A*C))/(2*A);
+    }
+
+    double dr = (fabs(s1)<=fabs(s2)?s1:s2);
+
+    // std::cout << "    A: " << A << std::endl;
+    // std::cout << "    B: " << B << std::endl;
+    // std::cout << "    C: " << C << std::endl;
+    // std::cout << "    s1: " << s1 << " s2: " << s2 << " dr: " << dr << std::endl;
+
+    circle1 = Geom::Circle( c1 - dr*n1, r1-dr );
+    circle2 = Geom::Circle( c2 + dr*n2, r2+dr );
+
+    // std::cout << "    C1: " << circle1 << std::endl;
+    // std::cout << "    C2: " << circle2 << std::endl;
+    // std::cout << "    d': " << Geom::Point( circle1.center() - circle2.center() ).length() << std::endl;
+
+    Geom::Line bisector( circle1.center(), circle2.center() );
+    std::vector<Geom::ShapeIntersection> points = circle1.intersect( bisector );
+    Geom::Point p0 = points[0].point();
+    Geom::Point p1 = points[1].point();
+    // std::cout << "    points: " << p0 << "; " << p1 << std::endl;
+    if( std::abs( Geom::distance( p0, circle2.center() ) - circle2.radius() ) <
+        std::abs( Geom::distance( p1, circle2.center() ) - circle2.radius() ) ) {
+        return p0;
+    } else {
+        return p1;
+    }
+}
+
+void extrapolate_join_internal(join_data const &jd, int alternative)
+{
+    // std::cout << "\nextrapolate_join: entrance: alternative: " << alternative << std::endl;
+    using namespace Geom;
+
+    Geom::Path &res = jd.res;
+    Geom::Curve const& incoming = res.back();
+    Geom::Curve const& outgoing = jd.outgoing.front();
+    Geom::Point startPt = incoming.finalPoint();
+    Geom::Point endPt = outgoing.initialPoint();
+    Geom::Point tang1 = jd.in_tang;
+    Geom::Point tang2 = jd.out_tang;
+    // width is half stroke-width
+    double width = jd.width, miter = jd.miter;
+
+    // std::cout << "  startPt: " << startPt << "  endPt: " << endPt << std::endl;
+    // std::cout << "  tang1: " << tang1 << "  tang2: " << tang2 <<  std::endl;
+    // std::cout << "    dot product: " << Geom::dot( tang1, tang2 ) <<  std::endl;
+    // std::cout << "  width: " << width << std::endl;
+
+    // CIRCLE CALCULATION TESTING
+    Geom::Circle circle1 = touching_circle(Geom::reverse(incoming.toSBasis()), 0.);
+    Geom::Circle circle2 = touching_circle(outgoing.toSBasis(), 0);
+    // std::cout << "  circle1: " << circle1 << std::endl;
+    // std::cout << "  circle2: " << circle2 << std::endl;
+    if( Geom::CubicBezier const * in_bezier = dynamic_cast<Geom::CubicBezier const*>(&incoming) ) {
+        Geom::Circle circle_test1 = touching_circle(*in_bezier, false);
+        if( !Geom::are_near( circle1, circle_test1, 0.01 ) ) {
+            // std::cout << "  Circle 1 error!!!!!!!!!!!!!!!!!" << std::endl;
+            // std::cout << "           " << circle_test1 << std::endl;
+        }
+    }
+    if( Geom::CubicBezier const * out_bezier = dynamic_cast<Geom::CubicBezier const*>(&outgoing) ) {
+        Geom::Circle circle_test2 = touching_circle(*out_bezier, true);
+        if( !Geom::are_near( circle2, circle_test2, 0.01 ) ) {
+            // std::cout << "  Circle 2 error!!!!!!!!!!!!!!!!!" << std::endl;
+            // std::cout << "           " << circle_test2 << std::endl;
+        }
+    }
+    // END TESTING
+
+    Geom::Point center1 = circle1.center();
+    double side1 = tang1[Geom::X]*(startPt[Geom::Y]-center1[Geom::Y]) - tang1[Geom::Y]*(startPt[Geom::X]-center1[Geom::X]);
+    // std::cout << "  side1: " << side1 << std::endl;
+
+    bool inc_ls = !circle1.center().isFinite();
+    bool out_ls = !circle2.center().isFinite();
+
+    std::vector<Geom::ShapeIntersection> points;
+
+    Geom::EllipticalArc *arc1 = nullptr;
+    Geom::EllipticalArc *arc2 = nullptr;
+    Geom::LineSegment *seg1 = nullptr;
+    Geom::LineSegment *seg2 = nullptr;
+    Geom::Point sol;
+    Geom::Point p1;
+    Geom::Point p2;
+
+    if (!inc_ls && !out_ls) {
+        // std::cout << " two circles" << std::endl;
+
+        // See if tangent is backwards (radius < width/2 and circle is inside stroke).
+        Geom::Point node_on_path = startPt + Geom::rot90(tang1)*width;
+        // std::cout << "  node_on_path: " << node_on_path << "  -------------- " << std::endl;
+        bool b1 = false;
+        bool b2 = false;
+        if (circle1.radius() < width &&  distance( circle1.center(), node_on_path ) < width) {
+            b1 = true;
+        }
+        if (circle2.radius() < width &&  distance( circle2.center(), node_on_path ) < width) {
+            b2 = true;
+        }
+        // std::cout << "  b1: " << (b1?"true":"false")
+        //           << "  b2: " << (b2?"true":"false") << std::endl;
+
+        // Two circles
+        points = circle1.intersect(circle2);
+
+        if (points.size() != 2) {
+            // std::cout << "   Circles do not intersect, do backup" << std::endl;
+            switch (alternative) {
+                case 1:
+                {
+                    // Fallback to round if one path has radius smaller than half line width.
+                    if(b1 || b2) {
+                        // std::cout << "At one least path has radius smaller than half line width." << std::endl;
+                        return( round_join(jd) );
+                    }
+
+                    Point p;
+                    if( circle2.contains( startPt ) && !circle1.contains( endPt ) ) {
+                        // std::cout << "Expand circle1" << std::endl;
+                        p = expand_circle( circle1, circle2, startPt, tang1 );
+                        points.emplace_back( 0, 0, p );
+                        points.emplace_back( 0, 0, p );
+                    } else if( circle1.contains( endPt ) && !circle2.contains( startPt ) ) {
+                        // std::cout << "Expand circle2" << std::endl;
+                        p = expand_circle( circle2, circle1, endPt, tang2 );
+                        points.emplace_back( 0, 0, p );
+                        points.emplace_back( 0, 0, p );
+                    } else {
+                        // std::cout << "Either both points inside or both outside" << std::endl;
+                        return( miter_clip_join(jd) );
+                    }
+                    break;
+                    
+                }
+                case 2:
+                {
+                    // Fallback to round if one path has radius smaller than half line width.
+                    if(b1 || b2) {
+                        // std::cout << "At one least path has radius smaller than half line width." << std::endl;
+                        return( round_join(jd) );
+                    }
+
+                    if( ( circle2.contains( startPt ) && !circle1.contains( endPt ) ) ||
+                        ( circle1.contains( endPt ) && !circle2.contains( startPt ) ) ) {
+                        
+                        Geom::Point apex = adjust_circles( circle1, circle2, startPt, endPt, tang1, tang2 );
+                        points.emplace_back( 0, 0, apex );
+                        points.emplace_back( 0, 0, apex );
+                    } else {
+                        // std::cout << "Either both points inside or both outside" << std::endl;
+                        return( miter_clip_join(jd) );
+                    }
+                        
+                    break;
+                }
+                case 3:
+                    if( side1 > 0 ) {
+                        Geom::Line secant(startPt, startPt + tang1);
+                        points = circle2.intersect(secant);
+                        circle1.setRadius(std::numeric_limits<float>::infinity());
+                        circle1.setCenter(Geom::Point(0,std::numeric_limits<float>::infinity()));
+                    } else {
+                        Geom::Line secant(endPt, endPt + tang2);
+                        points = circle1.intersect(secant);
+                        circle2.setRadius(std::numeric_limits<float>::infinity());
+                        circle2.setCenter(Geom::Point(0,std::numeric_limits<float>::infinity()));
+                    }
+                    break;
+
+
+                case 0:
+                default:
+                    // Do nothing
+                    break;
+            }
+        }
+
+        if (points.size() == 2) {
+            sol = pick_solution(points, tang2, endPt);
+            if( circle1.radius() != std::numeric_limits<float>::infinity() ) {
+                arc1 = circle1.arc(startPt, 0.5*(startPt+sol), sol);
+            } else {
+                seg1 = new Geom::LineSegment(startPt, sol);
+            }
+            if( circle2.radius() != std::numeric_limits<float>::infinity() ) {
+                arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt);
+            } else {
+                seg2 = new Geom::LineSegment(sol, endPt);
+            }
+        }
+    } else if (inc_ls && !out_ls) {
+        // Line and circle
+        // std::cout << " line circle" << std::endl;
+        points = circle2.intersect(Line(incoming.initialPoint(), incoming.finalPoint()));
+        if (points.size() == 2) {
+            sol = pick_solution(points, tang2, endPt);
+            arc2 = circle2.arc(sol, 0.5*(sol+endPt), endPt);
+        }
+    } else if (!inc_ls && out_ls) {
+        // Circle and line
+        // std::cout << " circle line" << std::endl;
+        points = circle1.intersect(Line(outgoing.initialPoint(), outgoing.finalPoint()));
+        if (points.size() == 2) {
+            sol = pick_solution(points, tang2, endPt);
+            arc1 = circle1.arc(startPt, 0.5*(sol+startPt), sol);
+        }
+    }
+    if (points.size() != 2) {
+        // std::cout << " no solutions" << std::endl;
+        // no solutions available, fall back to miter
+        return miter_join(jd);
+    }
+
+    // We have a solution, thus sol is defined.
+    p1 = sol;
+    
+    // See if we need to clip. Miter length is measured along a circular arc that is tangent to the
+    // bisector of the incoming and out going angles and passes through the end point (sol) of the
+    // line join.
+
+    // Center of circle is intersection of a line orthogonal to bisector and a line bisecting
+    // a chord connecting the path end point (point_on_path) and the join end point (sol).
+    Geom::Point point_on_path = startPt + Geom::rot90(tang1)*width;
+    Geom::Line bisector = make_angle_bisector_line(startPt, point_on_path, endPt);
+    Geom::Line ortho = make_orthogonal_line(point_on_path, bisector); 
+
+    Geom::LineSegment chord(point_on_path, sol);
+    Geom::Line bisector_chord =  make_bisector_line(chord);
+
+    Geom::Line limit_line;
+    double miter_limit = width * miter;
+    bool clipped = false;
+
+    if (are_parallel(bisector_chord, ortho)) {
+        // No intersection (can happen if curvatures are equal but opposite)
+        if (Geom::distance(point_on_path, sol) > miter_limit) {
+            clipped = true;
+            Geom::Point temp = bisector.versor();
+            Geom::Point limit_point = point_on_path + miter_limit * temp; 
+            limit_line = make_parallel_line( limit_point, ortho );
+        }
+    } else {
+        Geom::Point center =
+            Geom::intersection_point( bisector_chord.pointAt(0), bisector_chord.versor(),
+                                      ortho.pointAt(0),          ortho.versor() );
+        Geom::Coord radius = distance(center, point_on_path);
+        Geom::Circle circle_center(center, radius);
+
+        double limit_angle = miter_limit / radius;
+
+        Geom::Ray start_ray(center, point_on_path);
+        Geom::Ray end_ray(center, sol);
+        Geom::Line limit_line(center, 0); // Angle set below
+
+        if (Geom::cross(start_ray.versor(), end_ray.versor()) < 0) {
+            limit_line.setAngle(start_ray.angle() - limit_angle);
+        } else {
+            limit_line.setAngle(start_ray.angle() + limit_angle);
+        }
+
+        Geom::EllipticalArc *arc_center = circle_center.arc(point_on_path, 0.5*(point_on_path + sol), sol);
+        if (arc_center && arc_center->sweepAngle() > limit_angle) {
+            // We need to clip
+            clipped = true;
+
+            if (!inc_ls) {
+                // Incoming circular
+                points = circle1.intersect(limit_line);
+                if (points.size() == 2) {
+                    p1 = pick_solution(points, tang2, endPt);
+                    delete arc1;
+                    arc1 = circle1.arc(startPt, 0.5*(p1+startPt), p1);
+                }
+            } else {
+                p1 = Geom::intersection_point(startPt, tang1, limit_line.pointAt(0), limit_line.versor());
+            }
+
+            if (!out_ls) {
+                // Outgoing circular
+                points = circle2.intersect(limit_line);
+                if (points.size() == 2) {
+                    p2 = pick_solution(points, tang1, endPt);
+                    delete arc2;
+                    arc2 = circle2.arc(p2, 0.5*(p2+endPt), endPt);
+                }
+            } else {
+                p2 = Geom::intersection_point(endPt, tang2, limit_line.pointAt(0), limit_line.versor());
+            }
+        }
+    }    
+
+    // Add initial
+    if (arc1) {
+        res.append(*arc1);
+    } else if (seg1 ) {
+        res.append(*seg1);
+    } else {
+        // Straight line segment: move last point
+        res.setFinal(p1);
+    }
+
+    if (clipped) {
+        res.appendNew<Geom::LineSegment>(p2);
+    }
+
+    // Add outgoing
+    if (arc2) {
+        res.append(*arc2);
+        res.append(outgoing);
+    } else if (seg2 ) {
+        res.append(*seg2);
+        res.append(outgoing);
+    } else {
+        // Straight line segment:
+        res.appendNew<Geom::LineSegment>(outgoing.finalPoint());
+    }
+
+    // add the rest of the path
+    res.insert(res.end(), ++jd.outgoing.begin(), jd.outgoing.end());
+
+    delete arc1;
+    delete arc2;
+    delete seg1;
+    delete seg2;
+}
+
+void extrapolate_join(     join_data jd) { extrapolate_join_internal(jd, 0); }
+void extrapolate_join_alt1(join_data jd) { extrapolate_join_internal(jd, 1); }
+void extrapolate_join_alt2(join_data jd) { extrapolate_join_internal(jd, 2); }
+void extrapolate_join_alt3(join_data jd) { extrapolate_join_internal(jd, 3); }
+
+
+void tangents(Geom::Point tang[2], Geom::Curve const& incoming, Geom::Curve const& outgoing)
+{
+    Geom::Point tang1 = Geom::unitTangentAt(reverse(incoming.toSBasis()), 0.);
+    Geom::Point tang2 = outgoing.unitTangentAt(0.);
+    tang[0] = tang1, tang[1] = tang2;
+}
+
+// Offsetting a line segment is mathematically stable and quick to do
+Geom::LineSegment offset_line(Geom::LineSegment const& l, double width)
+{
+    Geom::Point tang1 = Geom::rot90(l.unitTangentAt(0));
+    Geom::Point tang2 = Geom::rot90(unitTangentAt(reverse(l.toSBasis()), 0.));
+
+    Geom::Point start = l.initialPoint() + tang1 * width;
+    Geom::Point end = l.finalPoint() - tang2 * width;
+    
+    return Geom::LineSegment(start, end);
+}
+
+void get_cubic_data(Geom::CubicBezier const& bez, double time, double& len, double& rad)
+{
+    // get derivatives
+    std::vector<Geom::Point> derivs = bez.pointAndDerivatives(time, 3);
+
+    Geom::Point der1 = derivs[1]; // first deriv (tangent vector)
+    Geom::Point der2 = derivs[2]; // second deriv (tangent's tangent)
+    double l = Geom::L2(der1); // length
+
+    len = rad = 0;
+
+    // TODO: we might want to consider using Geom::touching_circle to determine the
+    // curvature radius here. Less code duplication, but slower
+
+    if (Geom::are_near(l, 0, 1e-4)) {
+        l = Geom::L2(der2);
+        Geom::Point der3 = derivs.at(3); // try second time
+        if (Geom::are_near(l, 0, 1e-4)) {
+            l = Geom::L2(der3);
+            if (Geom::are_near(l, 0)) {
+                return; // this isn't a segment...
+            }
+        rad = 1e8;
+        } else {
+            rad = -l * (Geom::dot(der2, der2) / Geom::cross(der2, der3));
+        }
+    } else {
+        rad = -l * (Geom::dot(der1, der1) / Geom::cross(der1, der2));
+    }
+    len = l;
+}
+
+double _offset_cubic_stable_sub(
+        Geom::CubicBezier const& bez,
+        Geom::CubicBezier& c,
+        const Geom::Point& start_normal,
+        const Geom::Point& end_normal,
+        const Geom::Point& start_new,
+        const Geom::Point& end_new,
+        const double start_rad,
+        const double end_rad,
+        const double start_len,
+        const double end_len,
+        const double width,
+        const double width_correction) {
+    using Geom::X;
+    using Geom::Y;
+
+    double start_off = 1, end_off = 1;
+    // correction of the lengths of the tangent to the offset
+    if (!Geom::are_near(start_rad, 0))
+        start_off += (width + width_correction) / start_rad;
+    if (!Geom::are_near(end_rad, 0))
+        end_off += (width + width_correction) / end_rad;
+
+    // We don't change the direction of the control points
+    if (start_off < 0) {
+        start_off = 0;
+    }
+    if (end_off < 0) {
+        end_off = 0;
+    }
+    start_off *= start_len;
+    end_off *= end_len;
+    // --------
+
+    Geom::Point mid1_new = start_normal.ccw()*start_off;
+    mid1_new = Geom::Point(start_new[X] + mid1_new[X]/3., start_new[Y] + mid1_new[Y]/3.);
+    Geom::Point mid2_new = end_normal.ccw()*end_off;
+    mid2_new = Geom::Point(end_new[X] - mid2_new[X]/3., end_new[Y] - mid2_new[Y]/3.);
+
+    // create the estimate curve
+    c = Geom::CubicBezier(start_new, mid1_new, mid2_new, end_new);
+
+    // check the tolerance for our estimate to be a parallel curve
+
+    double worst_residual = 0;
+    for (size_t ii = 3; ii <= 7; ii+=2) {
+        const double t = static_cast<double>(ii) / 10;
+        const Geom::Point req = bez.pointAt(t);
+        const Geom::Point chk = c.pointAt(c.nearestTime(req));
+        const double current_residual = (chk-req).length() - std::abs(width);
+        if (std::abs(current_residual) > std::abs(worst_residual)) {
+            worst_residual = current_residual;
+        }
+    }
+    return worst_residual;
+}
+
+void offset_cubic(Geom::Path& p, Geom::CubicBezier const& bez, double width, double tol, size_t levels)
+{
+    using Geom::X;
+    using Geom::Y;
+
+    const Geom::Point start_pos = bez.initialPoint();
+    const Geom::Point end_pos = bez.finalPoint();
+
+    const Geom::Point start_normal = Geom::rot90(bez.unitTangentAt(0));
+    const Geom::Point end_normal = -Geom::rot90(Geom::unitTangentAt(Geom::reverse(bez.toSBasis()), 0.));
+
+    // offset the start and end control points out by the width
+    const Geom::Point start_new = start_pos + start_normal*width;
+    const Geom::Point end_new = end_pos + end_normal*width;
+
+    // --------
+    double start_rad, end_rad;
+    double start_len, end_len; // tangent lengths
+    get_cubic_data(bez, 0, start_len, start_rad);
+    get_cubic_data(bez, 1, end_len, end_rad);
+
+
+    Geom::CubicBezier c;
+
+    double best_width_correction = 0;
+    double best_residual = _offset_cubic_stable_sub(
+                bez, c,
+                start_normal, end_normal,
+                start_new, end_new,
+                start_rad, end_rad,
+                start_len, end_len,
+                width, best_width_correction);
+    double stepsize = std::abs(width)/2;
+    bool seen_success = false;
+    double stepsize_threshold = 0;
+    // std::cout << "Residual from " << best_residual << " ";
+    size_t ii = 0;
+    for (; ii < 100 && stepsize > stepsize_threshold; ++ii) {
+        bool success = false;
+        const double width_correction = best_width_correction - (best_residual > 0 ? 1 : -1) * stepsize;
+        Geom::CubicBezier current_curve;
+        const double residual = _offset_cubic_stable_sub(
+                    bez, current_curve,
+                    start_normal, end_normal,
+                    start_new, end_new,
+                    start_rad, end_rad,
+                    start_len, end_len,
+                    width, width_correction);
+        if (std::abs(residual) < std::abs(best_residual)) {
+            best_residual = residual;
+            best_width_correction = width_correction;
+            c = current_curve;
+            success = true;
+            if (std::abs(best_residual) < tol/4) {
+                break;
+            }
+        }
+
+        if (success) {
+            if (!seen_success) {
+                seen_success = true;
+                //std::cout << "Stepsize factor: " << std::abs(width) / stepsize << std::endl;
+                stepsize_threshold = stepsize / 1000;
+            }
+        }
+        else {
+            stepsize /= 2;
+        }
+        if (std::abs(best_width_correction) >= std::abs(width)/2) {
+            //break; // Seems to prevent some numerical instabilities, not clear if useful
+        }
+    }
+
+    // reached maximum recursive depth
+    // don't bother with any more correction
+    if (levels == 0) {
+        try {
+            p.append(c);
+        }
+        catch (...) {
+            if ((p.finalPoint() - c.initialPoint()).length() < 1e-6) {
+                c.setInitial(p.finalPoint());
+            }
+            else {
+                auto line = Geom::LineSegment(p.finalPoint(), c.initialPoint());
+                p.append(line);
+            }
+            p.append(c);
+        }
+
+        return;
+    }
+
+    // We find the point on our new curve (c) for which the distance between
+    // (c) and (bez) differs the most from the desired distance (width).
+    double worst_err = std::abs(best_residual);
+    double worst_time = .5;
+    for (size_t ii = 1; ii <= 9; ++ii) {
+        const double t = static_cast<double>(ii) / 10;
+        const Geom::Point req = bez.pointAt(t);
+        // We use the exact solution with nearestTime because it is numerically
+        // much more stable than simply assuming that the point on (c) closest
+        // to bez.pointAt(t) is given by c.pointAt(t)
+        const Geom::Point chk = c.pointAt(c.nearestTime(req));
+
+        Geom::Point const diff = req - chk;
+        const double err = std::abs(diff.length() - std::abs(width));
+        if (err > worst_err) {
+            worst_err = err;
+            worst_time = t;
+        }
+    }
+
+    if (worst_err < tol) {
+        if (Geom::are_near(start_new, p.finalPoint())) {
+            p.setFinal(start_new); // if it isn't near, we throw
+        }
+
+        // we're good, curve is accurate enough
+        p.append(c);
+        return;
+    } else {
+        // split the curve in two
+        std::pair<Geom::CubicBezier, Geom::CubicBezier> s = bez.subdivide(worst_time);
+        offset_cubic(p, s.first, width, tol, levels - 1);
+        offset_cubic(p, s.second, width, tol, levels - 1);
+    }
+}
+
+void offset_quadratic(Geom::Path& p, Geom::QuadraticBezier const& bez, double width, double tol, size_t levels)
+{
+    // cheat
+    // it's faster
+    // seriously
+    std::vector<Geom::Point> points = bez.controlPoints();
+    Geom::Point b1 = points[0] + (2./3) * (points[1] - points[0]);
+    Geom::Point b2 = b1 + (1./3) * (points[2] - points[0]);
+    Geom::CubicBezier cub = Geom::CubicBezier(points[0], b1, b2, points[2]);
+    offset_cubic(p, cub, width, tol, levels);
+}
+
+void offset_curve(Geom::Path& res, Geom::Curve const* current, double width, double tolerance)
+{
+    size_t levels = 8;
+
+    if (current->isDegenerate()) return; // don't do anything
+
+    // TODO: we can handle SVGEllipticalArc here as well, do that!
+
+    if (Geom::BezierCurve const *b = dynamic_cast<Geom::BezierCurve const*>(current)) {
+        size_t order = b->order();
+        switch (order) {
+            case 1:
+                res.append(offset_line(static_cast<Geom::LineSegment const&>(*current), width));
+                break;
+            case 2: {
+                Geom::QuadraticBezier const& q = static_cast<Geom::QuadraticBezier const&>(*current);
+                offset_quadratic(res, q, width, tolerance, levels);
+                break;
+            }
+            case 3: {
+                Geom::CubicBezier const& cb = static_cast<Geom::CubicBezier const&>(*current);
+                offset_cubic(res, cb, width, tolerance, levels);
+                break;
+            }
+            default: {
+                Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(current->toSBasis(), tolerance);
+                for (const auto & i : sbasis_path)
+                    offset_curve(res, &i, width, tolerance);
+                break;
+            }
+        }
+    } else {
+        Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(current->toSBasis(), 0.1);
+        for (const auto & i : sbasis_path)
+            offset_curve(res, &i, width, tolerance);
+    }
+}
+
+typedef void cap_func(Geom::PathBuilder& res, Geom::Path const& with_dir, Geom::Path const& against_dir, double width);
+
+void flat_cap(Geom::PathBuilder& res, Geom::Path const&, Geom::Path const& against_dir, double)
+{
+    res.lineTo(against_dir.initialPoint());
+}
+
+void round_cap(Geom::PathBuilder& res, Geom::Path const&, Geom::Path const& against_dir, double width)
+{
+    res.arcTo(width / 2., width / 2., 0., true, false, against_dir.initialPoint());
+}
+
+void square_cap(Geom::PathBuilder& res, Geom::Path const& with_dir, Geom::Path const& against_dir, double width)
+{
+    width /= 2.;
+    Geom::Point normal_1 = -Geom::unitTangentAt(Geom::reverse(with_dir.back().toSBasis()), 0.);
+    Geom::Point normal_2 = -against_dir[0].unitTangentAt(0.);
+    res.lineTo(with_dir.finalPoint() + normal_1*width);
+    res.lineTo(against_dir.initialPoint() + normal_2*width);
+    res.lineTo(against_dir.initialPoint());
+}
+
+void peak_cap(Geom::PathBuilder& res, Geom::Path const& with_dir, Geom::Path const& against_dir, double width)
+{
+    width /= 2.;
+    Geom::Point normal_1 = -Geom::unitTangentAt(Geom::reverse(with_dir.back().toSBasis()), 0.);
+    Geom::Point normal_2 = -against_dir[0].unitTangentAt(0.);
+    Geom::Point midpoint = ((with_dir.finalPoint() + normal_1*width) + (against_dir.initialPoint() + normal_2*width)) * 0.5;
+    res.lineTo(midpoint);
+    res.lineTo(against_dir.initialPoint());
+}
+
+} // namespace
+
+namespace Inkscape {
+
+Geom::PathVector outline(
+        Geom::Path const& input,
+        double width,
+        double miter,
+        LineJoinType join,
+        LineCapType butt,
+        double tolerance)
+{
+    if (input.size() == 0) return Geom::PathVector(); // nope, don't even try
+
+    Geom::PathBuilder res;
+    Geom::Path with_dir = half_outline(input, width/2., miter, join, tolerance);
+    Geom::Path against_dir = half_outline(input.reversed(), width/2., miter, join, tolerance);
+    res.moveTo(with_dir[0].initialPoint());
+    res.append(with_dir);
+
+    cap_func *cf;
+    switch (butt) {
+        case BUTT_ROUND:
+            cf = &round_cap;
+            break;
+        case BUTT_SQUARE:
+            cf = &square_cap;
+            break;
+        case BUTT_PEAK:
+            cf = &peak_cap;
+            break;
+        default:
+            cf = &flat_cap;
+    }
+
+    // glue caps
+    if (!input.closed()) {
+        cf(res, with_dir, against_dir, width);
+    } else {
+        res.closePath();
+        res.moveTo(against_dir.initialPoint());
+    }
+
+    res.append(against_dir);
+
+    if (!input.closed()) {
+        cf(res, against_dir, with_dir, width);
+    }
+
+    res.closePath();
+    res.flush();
+    return res.peek();
+}
+
+Geom::Path half_outline(
+        Geom::Path const& input,
+        double width,
+        double miter,
+        LineJoinType join,
+        double tolerance)
+{
+    if (tolerance <= 0) {
+        if (std::abs(width) > 0) {
+            tolerance = 5.0 * (std::abs(width)/100);
+        }
+        else {
+            tolerance = 1e-4;
+        }
+    }
+    Geom::Path res;
+    if (input.size() == 0) return res;
+
+    Geom::Point tang1 = input[0].unitTangentAt(0);
+    Geom::Point start = input.initialPoint() + tang1 * width;
+    Geom::Path temp;
+    Geom::Point tang[2];
+
+    res.setStitching(true);
+    temp.setStitching(true);
+
+    res.start(start);
+
+    // Do two curves at a time for efficiency, since the join function needs to know the outgoing curve as well
+    const Geom::Curve &closingline = input.back_closed();
+    const size_t k = (are_near(closingline.initialPoint(), closingline.finalPoint()) && input.closed() )
+            ?input.size_open():input.size_default();
+    
+    for (size_t u = 0; u < k; u += 2) {
+        temp.clear();
+
+        offset_curve(temp, &input[u], width, tolerance);
+
+        // on the first run through, there isn't a join
+        if (u == 0) {
+            res.append(temp);
+        } else {
+            tangents(tang, input[u-1], input[u]);
+            outline_join(res, temp, tang[0], tang[1], width, miter, join);
+        }
+
+        // odd number of paths
+        if (u < k - 1) {
+            temp.clear();
+            offset_curve(temp, &input[u+1], width, tolerance);
+            tangents(tang, input[u], input[u+1]);
+            outline_join(res, temp, tang[0], tang[1], width, miter, join);
+        }
+    }
+    if (input.closed()) {
+        Geom::Curve const &c1 = res.back();
+        Geom::Curve const &c2 = res.front();
+        temp.clear();
+        temp.append(c1);
+        Geom::Path temp2;
+        temp2.append(c2);
+        tangents(tang, input.back(), input.front());
+        outline_join(temp, temp2, tang[0], tang[1], width, miter, join);
+        res.erase(res.begin());
+        res.erase_last();
+        res.append(temp);
+        res.close();
+    }
+    return res;
+}
+
+void outline_join(Geom::Path &res, Geom::Path const& temp, Geom::Point in_tang, Geom::Point out_tang, double width, double miter, Inkscape::LineJoinType join)
+{
+    if (res.size() == 0 || temp.size() == 0)
+        return;
+    Geom::Curve const& outgoing = temp.front();
+    if (Geom::are_near(res.finalPoint(), outgoing.initialPoint(), 0.01)) {
+        // if the points are /that/ close, just ignore this one
+        res.setFinal(temp.initialPoint());
+        res.append(temp);
+        return;
+    }
+
+    join_data jd(res, temp, in_tang, out_tang, miter, width);
+    if (!(Geom::cross(in_tang, out_tang) > 0)) {
+        join = Inkscape::JOIN_BEVEL;
+    }
+    join_func *jf;
+    switch (join) {
+        case Inkscape::JOIN_BEVEL:
+            jf = &bevel_join;
+            break;
+        case Inkscape::JOIN_ROUND:
+            jf = &round_join;
+            break;
+        case Inkscape::JOIN_EXTRAPOLATE:
+            jf = &extrapolate_join;
+            break;
+        case Inkscape::JOIN_EXTRAPOLATE1:
+            jf = &extrapolate_join_alt1;
+            break;
+        case Inkscape::JOIN_EXTRAPOLATE2:
+            jf = &extrapolate_join_alt2;
+            break;
+        case Inkscape::JOIN_EXTRAPOLATE3:
+            jf = &extrapolate_join_alt3;
+            break;
+        case Inkscape::JOIN_MITER_CLIP:
+            jf = &miter_clip_join;
+            break;
+        default:
+            jf = &miter_join;
+    }
+    jf(jd);
+ }
+
+} // namespace Inkscape
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8 :
diff --git a/src/helper/geom-pathstroke.h b/src/helper/geom-pathstroke.h
new file mode 100644
index 0000000..73d35b4
--- /dev/null
+++ b/src/helper/geom-pathstroke.h
@@ -0,0 +1,109 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef INKSCAPE_HELPER_PATH_STROKE_H
+#define INKSCAPE_HELPER_PATH_STROKE_H
+
+/* Authors:
+ *   Liam P. White
+ *   Tavmjong Bah
+ *
+ * Copyright (C) 2014-2015 Authors
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/path.h>
+#include <2geom/pathvector.h>
+
+namespace Inkscape {
+
+enum LineJoinType {
+    JOIN_BEVEL,
+    JOIN_ROUND,
+    JOIN_MITER,
+    JOIN_MITER_CLIP,
+    JOIN_EXTRAPOLATE,
+    JOIN_EXTRAPOLATE1,
+    JOIN_EXTRAPOLATE2,
+    JOIN_EXTRAPOLATE3,
+};
+
+enum LineCapType {
+    BUTT_FLAT,
+    BUTT_ROUND,
+    BUTT_SQUARE,
+    BUTT_PEAK, // This is not a line ending supported by the SVG standard.
+};
+
+/**
+ * Strokes the path given by @a input.
+ * Joins may behave oddly if the width is negative.
+ *
+ * @param[in] input Input path.
+ * @param[in] width Stroke width.
+ * @param[in] miter Miter limit. Only used when @a join is one of JOIN_MITER, JOIN_MITER_CLIP, and JOIN_EXTRAPOLATE.
+ * @param[in] join  Line join type used during offset. Member of LineJoinType enum.
+ * @param[in] cap   Line cap type used during stroking. Member of LineCapType enum.
+ * @param[in] tolerance Tolerance, values smaller than 0 lead to automatic tolerance depending on width.
+ *
+ * @return Stroked path.
+ *         If the input path is closed, the resultant vector will contain two paths.
+ *         Otherwise, there should be only one in the output.
+ */
+Geom::PathVector outline(
+        Geom::Path const& input,
+        double width,
+        double miter,
+        LineJoinType join = JOIN_BEVEL,
+        LineCapType cap = BUTT_FLAT,
+        double tolerance = -1);
+
+/**
+ * Offset the input path by @a width.
+ * Joins may behave oddly if the width is negative.
+ *
+ * @param[in] input Input path.
+ * @param[in] width Amount to offset.
+ * @param[in] miter Miter limit. Only used when @a join is one of JOIN_MITER, JOIN_MITER_CLIP, and JOIN_EXTRAPOLATE.
+ * @param[in] join  Line join type used during offset. Member of LineJoinType enum.
+ * @param[in] tolerance Tolerance, values smaller than 0 lead to automatic tolerance depending on width.
+ *
+ * @return Offsetted output.
+ */
+Geom::Path half_outline(
+        Geom::Path const& input,
+        double width,
+        double miter,
+        LineJoinType join = JOIN_BEVEL,
+        double tolerance = -1);
+
+/**
+ * Builds a join on the provided path.
+ * Joins may behave oddly if the width is negative.
+ *
+ * @param[inout] res      The path to build the join on. 
+ *                        The outgoing path (or a portion thereof) will be appended after the join is created.
+ *                        Previous segments may be modified as an optimization, beware!
+ *
+ * @param[in]    outgoing The segment to append on the outgoing portion of the join.
+ * @param[in]    in_tang  The end tangent to consider on the input path.
+ * @param[in]    out_tang The begin tangent to consider on the output path.
+ * @param[in]    width
+ * @param[in]    miter
+ * @param[in]    join
+ */
+void outline_join(Geom::Path &res, Geom::Path const& outgoing, Geom::Point in_tang, Geom::Point out_tang, double width, double miter, LineJoinType join);
+
+} // namespace Inkscape
+
+#endif // INKSCAPE_HELPER_PATH_STROKE_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8 :
diff --git a/src/helper/geom-pathvector_nodesatellites.cpp b/src/helper/geom-pathvector_nodesatellites.cpp
new file mode 100644
index 0000000..2fd4125
--- /dev/null
+++ b/src/helper/geom-pathvector_nodesatellites.cpp
@@ -0,0 +1,248 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * PathVectorNodeSatellites a class to manage nodesatellites -per node extra data- in a pathvector
+ *//*
+ * Authors: see git history
+ * Jabiertxof
+ * Nathan Hurst
+ * Johan Engelen
+ * Josh Andler
+ * suv
+ * Mc-
+ * Liam P. White
+ * Krzysztof Kosiński
+ * This code is in public domain
+ *
+ * Copyright (C) 2018 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <helper/geom-pathvector_nodesatellites.h>
+#include <helper/geom.h>
+
+#include "util/units.h"
+
+Geom::PathVector PathVectorNodeSatellites::getPathVector() const
+{
+    return _pathvector;
+}
+
+void PathVectorNodeSatellites::setPathVector(Geom::PathVector pathv)
+{
+    _pathvector = pathv;
+}
+
+NodeSatellites PathVectorNodeSatellites::getNodeSatellites()
+{
+    return _nodesatellites;
+}
+
+void PathVectorNodeSatellites::setNodeSatellites(NodeSatellites nodesatellites)
+{
+    _nodesatellites = nodesatellites;
+}
+
+size_t PathVectorNodeSatellites::getTotalNodeSatellites()
+{
+    size_t counter = 0;
+    for (auto &_nodesatellite : _nodesatellites) {
+        counter += _nodesatellite.size();
+    }
+    return counter;
+}
+
+std::pair<size_t, size_t> PathVectorNodeSatellites::getIndexData(size_t index)
+{
+    size_t counter = 0;
+    for (size_t i = 0; i < _nodesatellites.size(); ++i) {
+        for (size_t j = 0; j < _nodesatellites[i].size(); ++j) {
+            if (index == counter) {
+                return std::make_pair(i,j);
+            }
+            counter++;
+        }
+    }
+    return std::make_pair(0,0);
+}
+
+void PathVectorNodeSatellites::setSelected(std::vector<size_t> selected)
+{
+    size_t counter = 0;
+    for (auto &_nodesatellite : _nodesatellites) {
+        for (auto &j : _nodesatellite) {
+            if (find (selected.begin(), selected.end(), counter) != selected.end()) {
+                j.setSelected(true);
+            } else {
+                j.setSelected(false);
+            }
+            counter++;
+        }
+    }
+}
+
+void PathVectorNodeSatellites::updateSteps(size_t steps, bool apply_no_radius, bool apply_with_radius,
+                                           bool only_selected)
+{
+    for (auto &_nodesatellite : _nodesatellites) {
+        for (auto &j : _nodesatellite) {
+            if ((!apply_no_radius && j.amount == 0) ||
+                (!apply_with_radius && j.amount != 0)) 
+            {
+                continue;
+            }
+            if (only_selected) {
+                if (j.selected) {
+                    j.steps = steps;
+                }
+            } else {
+                j.steps = steps;
+            }
+        }
+    }
+}
+
+void PathVectorNodeSatellites::updateAmount(double radius, bool apply_no_radius, bool apply_with_radius,
+                                            bool only_selected, bool use_knot_distance, bool flexible)
+{
+    double power = 0;
+    if (!flexible) {
+        power = radius;
+    } else {
+        power = radius / 100;
+    }
+    for (size_t i = 0; i < _nodesatellites.size(); ++i) {
+        for (size_t j = 0; j < _nodesatellites[i].size(); ++j) {
+            std::optional<size_t> previous_index = std::nullopt;
+            if (j == 0 && _pathvector[i].closed()) {
+                previous_index = count_path_nodes(_pathvector[i]) - 1;
+            } else if (!_pathvector[i].closed() || j != 0) {
+                previous_index = j - 1;
+            }
+            if (!_pathvector[i].closed() && j == 0) {
+                _nodesatellites[i][j].amount = 0;
+                continue;
+            }
+            if (count_path_nodes(_pathvector[i]) == j) {
+                continue;
+            }
+            if ((!apply_no_radius && _nodesatellites[i][j].amount == 0) ||
+                (!apply_with_radius && _nodesatellites[i][j].amount != 0)) {
+                continue;
+            }
+
+            if (_nodesatellites[i][j].selected || !only_selected) {
+                if (!use_knot_distance && !flexible) {
+                    if (previous_index) {
+                        _nodesatellites[i][j].amount =
+                            _nodesatellites[i][j].radToLen(power, _pathvector[i][*previous_index], _pathvector[i][j]);
+                        if (power && !_nodesatellites[i][j].amount) {
+                            g_warning("Seems a too high radius value");
+                        }
+                    } else {
+                        _nodesatellites[i][j].amount = 0.0;
+                    }
+                } else {
+                    _nodesatellites[i][j].amount = power;
+                }
+            }
+        }
+    }
+}
+
+void PathVectorNodeSatellites::convertUnit(Glib::ustring in, Glib::ustring to, bool apply_no_radius,
+                                           bool apply_with_radius)
+{
+    for (size_t i = 0; i < _nodesatellites.size(); ++i) {
+        for (size_t j = 0; j < _nodesatellites[i].size(); ++j) {
+            if (!_pathvector[i].closed() && j == 0) {
+                _nodesatellites[i][j].amount = 0;
+                continue;
+            }
+            if (count_path_nodes(_pathvector[i]) == j) {
+                continue;
+            }
+            if ((!apply_no_radius && _nodesatellites[i][j].amount == 0) ||
+                (!apply_with_radius && _nodesatellites[i][j].amount != 0)) {
+                continue;
+            }
+            _nodesatellites[i][j].amount =
+                Inkscape::Util::Quantity::convert(_nodesatellites[i][j].amount, in.c_str(), to.c_str());
+        }
+    }
+}
+
+void PathVectorNodeSatellites::updateNodeSatelliteType(NodeSatelliteType nodesatellitetype, bool apply_no_radius,
+                                                       bool apply_with_radius, bool only_selected)
+{
+    for (size_t i = 0; i < _nodesatellites.size(); ++i) {
+        for (size_t j = 0; j < _nodesatellites[i].size(); ++j) {
+            if ((!apply_no_radius && _nodesatellites[i][j].amount == 0) ||
+                (!apply_with_radius && _nodesatellites[i][j].amount != 0)) {
+                continue;
+            }
+            if (count_path_nodes(_pathvector[i]) == j) {
+                if (!only_selected) {
+                    _nodesatellites[i][j].nodesatellite_type = nodesatellitetype;
+                }
+                continue;
+            }
+            if (only_selected) {
+                if (_nodesatellites[i][j].selected) {
+                    _nodesatellites[i][j].nodesatellite_type = nodesatellitetype;
+                }
+            } else {
+                _nodesatellites[i][j].nodesatellite_type = nodesatellitetype;
+            }
+        }
+    }
+}
+
+void PathVectorNodeSatellites::recalculateForNewPathVector(Geom::PathVector const pathv, NodeSatellite const S)
+{
+    // pathv && _pathvector came here:
+    // * with different number of nodes
+    // * without empty subpats
+    // * _pathvector and nodesatellites (old data) are paired
+    NodeSatellites nodesatellites;
+    bool found = false;
+    // TODO evaluate fix on nodes at same position
+    // size_t number_nodes = count_pathvector_nodes(pathv);
+    // size_t previous_number_nodes = getTotalNodeSatellites();
+    for (const auto & i : pathv) {
+        std::vector<NodeSatellite> path_nodesatellites;
+        size_t count = count_path_nodes(i);
+        for (size_t j = 0; j < count; j++) {
+            found = false;
+            for (size_t k = 0; k < _pathvector.size(); k++) {
+                size_t count2 = count_path_nodes(_pathvector[k]);
+                for (size_t l = 0; l < count2; l++) {
+                    if (Geom::are_near(_pathvector[k][l].initialPoint(),  i[j].initialPoint())) {
+                        path_nodesatellites.push_back(_nodesatellites[k][l]);
+                        found = true;
+                        break;
+                    }
+                }
+                if (found) {
+                    break;
+                }
+            }
+            if (!found) {
+                path_nodesatellites.push_back(S);
+            }
+        }
+        nodesatellites.push_back(path_nodesatellites);
+    }
+    setPathVector(pathv);
+    setNodeSatellites(nodesatellites);
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom-pathvector_nodesatellites.h b/src/helper/geom-pathvector_nodesatellites.h
new file mode 100644
index 0000000..4631b32
--- /dev/null
+++ b/src/helper/geom-pathvector_nodesatellites.h
@@ -0,0 +1,61 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * \brief PathVectorNodeSatellites a class to manage nodesatellites -per node extra data- in a pathvector
+ *//*
+ * Authors: see git history
+ * Jabiertxof
+ * Nathan Hurst
+ * Johan Engelen
+ * Josh Andler
+ * suv
+ * Mc-
+ * Liam P. White
+ * Krzysztof Kosiński
+ *
+ * Copyright (C) 2017 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifndef SEEN_PATHVECTORSATELLITES_H
+#define SEEN_PATHVECTORSATELLITES_H
+
+#include <2geom/path.h>
+#include <2geom/pathvector.h>
+#include <helper/geom-nodesatellite.h>
+
+typedef std::vector<std::vector<NodeSatellite>> NodeSatellites;
+///@brief PathVectorNodeSatellites a class to manage nodesatellites in a pathvector
+class PathVectorNodeSatellites
+{
+public:
+    Geom::PathVector getPathVector() const;
+    void setPathVector(Geom::PathVector pathv);
+    NodeSatellites getNodeSatellites();
+    void setNodeSatellites(NodeSatellites nodesatellites);
+    size_t getTotalNodeSatellites();
+    void setSelected(std::vector<size_t> selected);
+    void updateSteps(size_t steps, bool apply_no_radius, bool apply_with_radius, bool only_selected);
+    void updateAmount(double radius, bool apply_no_radius, bool apply_with_radius, bool only_selected, 
+                      bool use_knot_distance, bool flexible);
+    void convertUnit(Glib::ustring in, Glib::ustring to, bool apply_no_radius, bool apply_with_radius);
+    void updateNodeSatelliteType(NodeSatelliteType nodesatellitetype, bool apply_no_radius, bool apply_with_radius,
+                                 bool only_selected);
+    std::pair<size_t, size_t> getIndexData(size_t index);
+    void recalculateForNewPathVector(Geom::PathVector const pathv, NodeSatellite const S);
+
+private:
+    Geom::PathVector _pathvector;
+    NodeSatellites _nodesatellites;
+};
+
+#endif //SEEN_PATHVECTORSATELLITES_H
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom.cpp b/src/helper/geom.cpp
new file mode 100644
index 0000000..da2b63c
--- /dev/null
+++ b/src/helper/geom.cpp
@@ -0,0 +1,954 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * Specific geometry functions for Inkscape, not provided my lib2geom.
+ *
+ * Author:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2008 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <algorithm>
+#include "helper/geom.h"
+#include "helper/geom-curves.h"
+#include <2geom/curves.h>
+#include <2geom/sbasis-to-bezier.h>
+
+using Geom::X;
+using Geom::Y;
+
+//#################################################################################
+// BOUNDING BOX CALCULATIONS
+
+/* Fast bbox calculation */
+/* Thanks to Nathan Hurst for suggesting it */
+static void
+cubic_bbox (Geom::Coord x000, Geom::Coord y000, Geom::Coord x001, Geom::Coord y001, Geom::Coord x011, Geom::Coord y011, Geom::Coord x111, Geom::Coord y111, Geom::Rect &bbox)
+{
+    Geom::Coord a, b, c, D;
+
+    bbox[0].expandTo(x111);
+    bbox[1].expandTo(y111);
+
+    // It already contains (x000,y000) and (x111,y111)
+    // All points of the Bezier lie in the convex hull of (x000,y000), (x001,y001), (x011,y011) and (x111,y111)
+    // So, if it also contains (x001,y001) and (x011,y011) we don't have to compute anything else!
+    // Note that we compute it for the X and Y range separately to make it easier to use them below
+    bool containsXrange = bbox[0].contains(x001) && bbox[0].contains(x011);
+    bool containsYrange = bbox[1].contains(y001) && bbox[1].contains(y011);
+
+    /*
+     * xttt = s * (s * (s * x000 + t * x001) + t * (s * x001 + t * x011)) + t * (s * (s * x001 + t * x011) + t * (s * x011 + t * x111))
+     * xttt = s * (s2 * x000 + s * t * x001 + t * s * x001 + t2 * x011) + t * (s2 * x001 + s * t * x011 + t * s * x011 + t2 * x111)
+     * xttt = s * (s2 * x000 + 2 * st * x001 + t2 * x011) + t * (s2 * x001 + 2 * st * x011 + t2 * x111)
+     * xttt = s3 * x000 + 2 * s2t * x001 + st2 * x011 + s2t * x001 + 2st2 * x011 + t3 * x111
+     * xttt = s3 * x000 + 3s2t * x001 + 3st2 * x011 + t3 * x111
+     * xttt = s3 * x000 + (1 - s) 3s2 * x001 + (1 - s) * (1 - s) * 3s * x011 + (1 - s) * (1 - s) * (1 - s) * x111
+     * xttt = s3 * x000 + (3s2 - 3s3) * x001 + (3s - 6s2 + 3s3) * x011 + (1 - 2s + s2 - s + 2s2 - s3) * x111
+     * xttt = (x000 - 3 * x001 + 3 * x011 -     x111) * s3 +
+     *        (       3 * x001 - 6 * x011 + 3 * x111) * s2 +
+     *        (                  3 * x011 - 3 * x111) * s  +
+     *        (                                 x111)
+     * xttt' = (3 * x000 - 9 * x001 +  9 * x011 - 3 * x111) * s2 +
+     *         (           6 * x001 - 12 * x011 + 6 * x111) * s  +
+     *         (                       3 * x011 - 3 * x111)
+     */
+
+    if (!containsXrange) {
+        a = 3 * x000 - 9 * x001 + 9 * x011 - 3 * x111;
+        b = 6 * x001 - 12 * x011 + 6 * x111;
+        c = 3 * x011 - 3 * x111;
+
+        /*
+        * s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a;
+        */
+        if (fabs (a) < Geom::EPSILON) {
+            /* s = -c / b */
+            if (fabs (b) > Geom::EPSILON) {
+                double s;
+                s = -c / b;
+                if ((s > 0.0) && (s < 1.0)) {
+                    double t = 1.0 - s;
+                    double xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
+                    bbox[0].expandTo(xttt);
+                }
+            }
+        } else {
+            /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
+            D = b * b - 4 * a * c;
+            if (D >= 0.0) {
+                Geom::Coord d, s, t, xttt;
+                /* Have solution */
+                d = sqrt (D);
+                s = (-b + d) / (2 * a);
+                if ((s > 0.0) && (s < 1.0)) {
+                    t = 1.0 - s;
+                    xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
+                    bbox[0].expandTo(xttt);
+                }
+                s = (-b - d) / (2 * a);
+                if ((s > 0.0) && (s < 1.0)) {
+                    t = 1.0 - s;
+                    xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
+                    bbox[0].expandTo(xttt);
+                }
+            }
+        }
+    }
+
+    if (!containsYrange) {
+        a = 3 * y000 - 9 * y001 + 9 * y011 - 3 * y111;
+        b = 6 * y001 - 12 * y011 + 6 * y111;
+        c = 3 * y011 - 3 * y111;
+
+        if (fabs (a) < Geom::EPSILON) {
+            /* s = -c / b */
+            if (fabs (b) > Geom::EPSILON) {
+                double s;
+                s = -c / b;
+                if ((s > 0.0) && (s < 1.0)) {
+                    double t = 1.0 - s;
+                    double yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
+                    bbox[1].expandTo(yttt);
+                }
+            }
+        } else {
+            /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
+            D = b * b - 4 * a * c;
+            if (D >= 0.0) {
+                Geom::Coord d, s, t, yttt;
+                /* Have solution */
+                d = sqrt (D);
+                s = (-b + d) / (2 * a);
+                if ((s > 0.0) && (s < 1.0)) {
+                    t = 1.0 - s;
+                    yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
+                    bbox[1].expandTo(yttt);
+                }
+                s = (-b - d) / (2 * a);
+                if ((s > 0.0) && (s < 1.0)) {
+                    t = 1.0 - s;
+                    yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
+                    bbox[1].expandTo(yttt);
+                }
+            }
+        }
+    }
+}
+
+Geom::OptRect
+bounds_fast_transformed(Geom::PathVector const & pv, Geom::Affine const & t)
+{
+    return bounds_exact_transformed(pv, t); //use this as it is faster for now! :)
+//    return Geom::bounds_fast(pv * t);
+}
+
+Geom::OptRect
+bounds_exact_transformed(Geom::PathVector const & pv, Geom::Affine const & t)
+{
+    if (pv.empty())
+        return Geom::OptRect();
+
+    Geom::Point initial = pv.front().initialPoint() * t;
+    Geom::Rect bbox(initial, initial);        // obtain well defined bbox as starting point to unionWith
+
+    for (const auto & it : pv) {
+        bbox.expandTo(it.initialPoint() * t);
+
+        // don't loop including closing segment, since that segment can never increase the bbox
+        for (Geom::Path::const_iterator cit = it.begin(); cit != it.end_open(); ++cit) {
+            Geom::Curve const &c = *cit;
+
+            unsigned order = 0;
+            if (Geom::BezierCurve const* b = dynamic_cast<Geom::BezierCurve const*>(&c)) {
+                order = b->order();
+            }
+
+            if (order == 1) { // line segment
+                bbox.expandTo(c.finalPoint() * t);
+
+            // TODO: we can make the case for quadratics faster by degree elevating them to
+            // cubic and then taking the bbox of that.
+
+            } else if (order == 3) { // cubic bezier
+                Geom::CubicBezier const &cubic_bezier = static_cast<Geom::CubicBezier const&>(c);
+                Geom::Point c0 = cubic_bezier[0] * t;
+                Geom::Point c1 = cubic_bezier[1] * t;
+                Geom::Point c2 = cubic_bezier[2] * t;
+                Geom::Point c3 = cubic_bezier[3] * t;
+                cubic_bbox(c0[0], c0[1], c1[0], c1[1], c2[0], c2[1], c3[0], c3[1], bbox);
+            } else {
+                // should handle all not-so-easy curves:
+                Geom::Curve *ctemp = cit->transformed(t);
+                bbox.unionWith( ctemp->boundsExact());
+                delete ctemp;
+            }
+        }
+    }
+    //return Geom::bounds_exact(pv * t);
+    return bbox;
+}
+
+bool 
+pathv_similar(Geom::PathVector const &apv, Geom::PathVector const &bpv, double precission) 
+{
+    if (apv == bpv) {
+        return true;
+    }
+    size_t totala = apv.curveCount();
+    if (totala != bpv.curveCount()) {
+        return false;
+    }
+    std::vector<Geom::Coord> pos;
+    for (size_t i = 0; i < totala; i++) {
+        Geom::Point pointa = apv.pointAt(float(i)+0.2);
+        Geom::Point pointb = bpv.pointAt(float(i)+0.2);
+        Geom::Point pointc = apv.pointAt(float(i)+0.4);
+        Geom::Point pointd = bpv.pointAt(float(i)+0.4);
+        Geom::Point pointe = apv.pointAt(float(i));
+        Geom::Point pointf = bpv.pointAt(float(i));
+        if (!Geom::are_near(pointa[Geom::X], pointb[Geom::X], precission) ||
+            !Geom::are_near(pointa[Geom::Y], pointb[Geom::Y], precission) ||
+            !Geom::are_near(pointc[Geom::X], pointd[Geom::X], precission) ||
+            !Geom::are_near(pointc[Geom::Y], pointd[Geom::Y], precission) ||
+            !Geom::are_near(pointe[Geom::X], pointf[Geom::X], precission) ||
+            !Geom::are_near(pointe[Geom::Y], pointf[Geom::Y], precission)) 
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+static void
+geom_line_wind_distance (Geom::Coord x0, Geom::Coord y0, Geom::Coord x1, Geom::Coord y1, Geom::Point const &pt, int *wind, Geom::Coord *best)
+{
+    Geom::Coord Ax, Ay, Bx, By, Dx, Dy, s;
+    Geom::Coord dist2;
+
+    /* Find distance */
+    Ax = x0;
+    Ay = y0;
+    Bx = x1;
+    By = y1;
+    Dx = x1 - x0;
+    Dy = y1 - y0;
+    const Geom::Coord Px = pt[X];
+    const Geom::Coord Py = pt[Y];
+
+    if (best) {
+        s = ((Px - Ax) * Dx + (Py - Ay) * Dy) / (Dx * Dx + Dy * Dy);
+        if (s <= 0.0) {
+            dist2 = (Px - Ax) * (Px - Ax) + (Py - Ay) * (Py - Ay);
+        } else if (s >= 1.0) {
+            dist2 = (Px - Bx) * (Px - Bx) + (Py - By) * (Py - By);
+        } else {
+            Geom::Coord Qx, Qy;
+            Qx = Ax + s * Dx;
+            Qy = Ay + s * Dy;
+            dist2 = (Px - Qx) * (Px - Qx) + (Py - Qy) * (Py - Qy);
+        }
+
+        if (dist2 < (*best * *best)) *best = sqrt (dist2);
+    }
+
+    if (wind) {
+        /* Find wind */
+        if ((Ax >= Px) && (Bx >= Px)) return;
+        if ((Ay >= Py) && (By >= Py)) return;
+        if ((Ay < Py) && (By < Py)) return;
+        if (Ay == By) return;
+        /* Ctach upper y bound */
+        if (Ay == Py) {
+            if (Ax < Px) *wind -= 1;
+            return;
+        } else if (By == Py) {
+            if (Bx < Px) *wind += 1;
+            return;
+        } else {
+            Geom::Coord Qx;
+            /* Have to calculate intersection */
+            Qx = Ax + Dx * (Py - Ay) / Dy;
+            if (Qx < Px) {
+                *wind += (Dy > 0.0) ? 1 : -1;
+            }
+        }
+    }
+}
+
+static void
+geom_cubic_bbox_wind_distance (Geom::Coord x000, Geom::Coord y000,
+                 Geom::Coord x001, Geom::Coord y001,
+                 Geom::Coord x011, Geom::Coord y011,
+                 Geom::Coord x111, Geom::Coord y111,
+                 Geom::Point const &pt,
+                 Geom::Rect *bbox, int *wind, Geom::Coord *best,
+                 Geom::Coord tolerance)
+{
+    Geom::Coord x0, y0, x1, y1, len2;
+    int needdist, needwind;
+
+    const Geom::Coord Px = pt[X];
+    const Geom::Coord Py = pt[Y];
+
+    needdist = 0;
+    needwind = 0;
+
+    if (bbox) cubic_bbox (x000, y000, x001, y001, x011, y011, x111, y111, *bbox);
+
+    x0 = std::min (x000, x001);
+    x0 = std::min (x0, x011);
+    x0 = std::min (x0, x111);
+    y0 = std::min (y000, y001);
+    y0 = std::min (y0, y011);
+    y0 = std::min (y0, y111);
+    x1 = std::max (x000, x001);
+    x1 = std::max (x1, x011);
+    x1 = std::max (x1, x111);
+    y1 = std::max (y000, y001);
+    y1 = std::max (y1, y011);
+    y1 = std::max (y1, y111);
+
+    if (best) {
+        /* Quickly adjust to endpoints */
+        len2 = (x000 - Px) * (x000 - Px) + (y000 - Py) * (y000 - Py);
+        if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
+        len2 = (x111 - Px) * (x111 - Px) + (y111 - Py) * (y111 - Py);
+        if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
+
+        if (((x0 - Px) < *best) && ((y0 - Py) < *best) && ((Px - x1) < *best) && ((Py - y1) < *best)) {
+            /* Point is inside sloppy bbox */
+            /* Now we have to decide, whether subdivide */
+            /* fixme: (Lauris) */
+            if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
+                needdist = 1;
+            }
+        }
+    }
+    if (!needdist && wind) {
+        if ((y1 >= Py) && (y0 < Py) && (x0 < Px)) {
+            /* Possible intersection at the left */
+            /* Now we have to decide, whether subdivide */
+            /* fixme: (Lauris) */
+            if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
+                needwind = 1;
+            }
+        }
+    }
+
+    if (needdist || needwind) {
+        Geom::Coord x00t, x0tt, xttt, x1tt, x11t, x01t;
+        Geom::Coord y00t, y0tt, yttt, y1tt, y11t, y01t;
+        Geom::Coord s, t;
+
+        t = 0.5;
+        s = 1 - t;
+
+        x00t = s * x000 + t * x001;
+        x01t = s * x001 + t * x011;
+        x11t = s * x011 + t * x111;
+        x0tt = s * x00t + t * x01t;
+        x1tt = s * x01t + t * x11t;
+        xttt = s * x0tt + t * x1tt;
+
+        y00t = s * y000 + t * y001;
+        y01t = s * y001 + t * y011;
+        y11t = s * y011 + t * y111;
+        y0tt = s * y00t + t * y01t;
+        y1tt = s * y01t + t * y11t;
+        yttt = s * y0tt + t * y1tt;
+
+        geom_cubic_bbox_wind_distance (x000, y000, x00t, y00t, x0tt, y0tt, xttt, yttt, pt, nullptr, wind, best, tolerance);
+        geom_cubic_bbox_wind_distance (xttt, yttt, x1tt, y1tt, x11t, y11t, x111, y111, pt, nullptr, wind, best, tolerance);
+    } else {
+        geom_line_wind_distance (x000, y000, x111, y111, pt, wind, best);
+    }
+}
+
+static void
+geom_curve_bbox_wind_distance(Geom::Curve const & c, Geom::Affine const &m,
+                 Geom::Point const &pt,
+                 Geom::Rect *bbox, int *wind, Geom::Coord *dist,
+                 Geom::Coord tolerance, Geom::Rect const *viewbox,
+                 Geom::Point &p0) // pass p0 through as it represents the last endpoint added (the finalPoint of last curve)
+{
+    unsigned order = 0;
+    if (Geom::BezierCurve const* b = dynamic_cast<Geom::BezierCurve const*>(&c)) {
+        order = b->order();
+    }
+    if (order == 1) {
+        Geom::Point pe = c.finalPoint() * m;
+        if (bbox) {
+            bbox->expandTo(pe);
+        }
+        if (dist || wind) {
+            if (wind) { // we need to pick fill, so do what we're told
+                geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
+            } else { // only stroke is being picked; skip this segment if it's totally outside the viewbox
+                Geom::Rect swept(p0, pe);
+                if (!viewbox || swept.intersects(*viewbox))
+                    geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
+            }
+        }
+        p0 = pe;
+    }
+    else if (order == 3) {
+        Geom::CubicBezier const& cubic_bezier = static_cast<Geom::CubicBezier const&>(c);
+        Geom::Point p1 = cubic_bezier[1] * m;
+        Geom::Point p2 = cubic_bezier[2] * m;
+        Geom::Point p3 = cubic_bezier[3] * m;
+
+        // get approximate bbox from handles (convex hull property of beziers):
+        Geom::Rect swept(p0, p3);
+        swept.expandTo(p1);
+        swept.expandTo(p2);
+
+        if (!viewbox || swept.intersects(*viewbox)) { // we see this segment, so do full processing
+            geom_cubic_bbox_wind_distance ( p0[X], p0[Y],
+                                            p1[X], p1[Y],
+                                            p2[X], p2[Y],
+                                            p3[X], p3[Y],
+                                            pt,
+                                            bbox, wind, dist, tolerance);
+        } else {
+            if (wind) { // if we need fill, we can just pretend it's a straight line
+                geom_line_wind_distance (p0[X], p0[Y], p3[X], p3[Y], pt, wind, dist);
+            } else { // otherwise, skip it completely
+            }
+        }
+        p0 = p3;
+    } else { 
+        //this case handles sbasis as well as all other curve types
+        Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c.toSBasis(), 0.1);
+
+        //recurse to convert the new path resulting from the sbasis to svgd
+        for (const auto & iter : sbasis_path) {
+            geom_curve_bbox_wind_distance(iter, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
+        }
+    }
+}
+
+bool 
+pointInTriangle(Geom::Point const &p, Geom::Point const &p1, Geom::Point const &p2, Geom::Point const &p3)
+{
+    //http://totologic.blogspot.com.es/2014/01/accurate-point-in-triangle-test.html
+    using Geom::X;
+    using Geom::Y;
+    double denominator = (p1[X]*(p2[Y] - p3[Y]) + p1[Y]*(p3[X] - p2[X]) + p2[X]*p3[Y] - p2[Y]*p3[X]);
+    double t1 = (p[X]*(p3[Y] - p1[Y]) + p[Y]*(p1[X] - p3[X]) - p1[X]*p3[Y] + p1[Y]*p3[X]) / denominator;
+    double t2 = (p[X]*(p2[Y] - p1[Y]) + p[Y]*(p1[X] - p2[X]) - p1[X]*p2[Y] + p1[Y]*p2[X]) / -denominator;
+    double s = t1 + t2;
+
+    return 0 <= t1 && t1 <= 1 && 0 <= t2 && t2 <= 1 && s <= 1;
+}
+
+
+/* Calculates...
+   and returns ... in *wind and the distance to ... in *dist.
+   Returns bounding box in *bbox if bbox!=NULL.
+ */
+void
+pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, Geom::Affine const &m, Geom::Point const &pt,
+                         Geom::Rect *bbox, int *wind, Geom::Coord *dist,
+                         Geom::Coord tolerance, Geom::Rect const *viewbox)
+{
+    if (pathv.empty()) {
+        if (wind) *wind = 0;
+        if (dist) *dist = Geom::infinity();
+        return;
+    }
+
+    // remember last point of last curve
+    Geom::Point p0(0,0);
+
+    // remembering the start of subpath
+    Geom::Point p_start(0,0);
+    bool start_set = false;
+
+    for (const auto & it : pathv) {
+
+        if (start_set) { // this is a new subpath
+            if (wind && (p0 != p_start)) // for correct fill picking, each subpath must be closed
+                geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
+        }
+        p0 = it.initialPoint() * m;
+        p_start = p0;
+        start_set = true;
+        if (bbox) {
+            bbox->expandTo(p0);
+        }
+
+        // loop including closing segment if path is closed
+        for (Geom::Path::const_iterator cit = it.begin(); cit != it.end_default(); ++cit) {
+            geom_curve_bbox_wind_distance(*cit, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
+        }
+    }
+
+    if (start_set) { 
+        if (wind && (p0 != p_start)) // for correct picking, each subpath must be closed
+            geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
+    }
+}
+
+//#################################################################################
+
+/**
+ * Basic check on intersecting path vectors
+ */
+bool is_intersecting(Geom::PathVector const&a, Geom::PathVector const&b) {
+    for (auto &node : b.nodes()) {
+        if (a.winding(node)) {
+            return true;
+        }
+    }
+   for (auto &node : a.nodes()) {
+        if (b.winding(node)) {
+            return true;
+        }
+    }
+    return false;
+}
+
+/*
+ * Converts all segments in all paths to Geom::LineSegment or Geom::HLineSegment or
+ * Geom::VLineSegment or Geom::CubicBezier.
+ */
+Geom::PathVector
+pathv_to_linear_and_cubic_beziers( Geom::PathVector const &pathv )
+{
+    Geom::PathVector output;
+
+    for (const auto & pit : pathv) {
+        output.push_back( Geom::Path() );
+        output.back().setStitching(true);
+        output.back().start( pit.initialPoint() );
+
+        for (Geom::Path::const_iterator cit = pit.begin(); cit != pit.end_open(); ++cit) {
+            if (is_straight_curve(*cit)) {
+                Geom::LineSegment l(cit->initialPoint(), cit->finalPoint());
+                output.back().append(l);
+            } else {
+                Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&*cit);
+                if (curve && curve->order() == 3) {
+                    Geom::CubicBezier b((*curve)[0], (*curve)[1], (*curve)[2], (*curve)[3]);
+                    output.back().append(b);
+                } else {
+                    // convert all other curve types to cubicbeziers
+                    Geom::Path cubicbezier_path = Geom::cubicbezierpath_from_sbasis(cit->toSBasis(), 0.1);
+                    cubicbezier_path.close(false);
+                    output.back().append(cubicbezier_path);
+                }
+            }
+        }
+        
+        output.back().close( pit.closed() );
+    }
+    
+    return output;
+}
+
+/*
+ * Converts all segments in all paths to Geom::LineSegment.  There is an intermediate
+ * stage where some may be converted to beziers.  maxdisp is the maximum displacement from
+ * the line segment to the bezier curve; ** maxdisp is not used at this moment **.
+ *
+ * This is NOT a terribly fast method, but it should give a solution close to the one with the
+ * fewest points.
+ */
+Geom::PathVector
+pathv_to_linear( Geom::PathVector const &pathv, double /*maxdisp*/)
+{
+    Geom::PathVector output;
+    Geom::PathVector tmppath = pathv_to_linear_and_cubic_beziers(pathv);
+    
+    // Now all path segments are either already lines, or they are beziers.
+
+    for (const auto & pit : tmppath) {
+        output.push_back( Geom::Path() );
+        output.back().start( pit.initialPoint() );
+        output.back().close( pit.closed() );
+
+        for (Geom::Path::const_iterator cit = pit.begin(); cit != pit.end_open(); ++cit) {
+            if (is_straight_curve(*cit)) {
+                Geom::LineSegment ls(cit->initialPoint(), cit->finalPoint());
+                output.back().append(ls);
+            } 
+            else { /* all others must be Bezier curves */
+                Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&*cit);
+                std::vector<Geom::Point> bzrpoints = curve->controlPoints();
+                Geom::Point A = bzrpoints[0];
+                Geom::Point B = bzrpoints[1];
+                Geom::Point C = bzrpoints[2];
+                Geom::Point D = bzrpoints[3];
+                std::vector<Geom::Point> pointlist;
+                pointlist.push_back(A);
+                recursive_bezier4(
+                   A[X], A[Y], 
+                   B[X], B[Y], 
+                   C[X], C[Y], 
+                   D[X], D[Y],
+                   pointlist, 
+                   0);
+                pointlist.push_back(D);
+                Geom::Point r1 = pointlist[0];
+                for (unsigned int i=1; i<pointlist.size();i++){
+                   Geom::Point prev_r1 = r1;
+                   r1 = pointlist[i];
+                   Geom::LineSegment ls(prev_r1, r1);
+                   output.back().append(ls);
+                }
+                pointlist.clear();
+           }
+        }
+    }
+    
+    return output;
+}
+
+/*
+ * Converts all segments in all paths to Geom Cubic bezier.
+ * This is used in lattice2 LPE, maybe is better move the function to the effect
+ * But maybe could be usable by others, so i put here.
+ * The straight curve part is needed as is for the effect to work appropriately
+ */
+Geom::PathVector
+pathv_to_cubicbezier( Geom::PathVector const &pathv)
+{
+    Geom::PathVector output;
+    double cubicGap = 0.01;
+    for (const auto & pit : pathv) {
+        if (pit.empty()) {
+            continue;
+        }
+        output.push_back( Geom::Path() );
+        output.back().start( pit.initialPoint() );
+        output.back().close( pit.closed() );
+        bool end_open = false;
+        if (pit.closed()) {
+            const Geom::Curve &closingline = pit.back_closed();
+            if (!are_near(closingline.initialPoint(), closingline.finalPoint())) {
+                end_open = true;
+            }
+        }
+        Geom::Path pitCubic = (Geom::Path)pit;
+        if(end_open && pit.closed()){
+            pitCubic.close(false);
+            pitCubic.appendNew<Geom::LineSegment>( pitCubic.initialPoint() );
+            pitCubic.close(true);
+        }
+        for (Geom::Path::iterator cit = pitCubic.begin(); cit != pitCubic.end_open(); ++cit) {
+            if (is_straight_curve(*cit)) {
+                Geom::CubicBezier b(cit->initialPoint(), cit->pointAt(0.3334) + Geom::Point(cubicGap,cubicGap), cit->finalPoint(), cit->finalPoint());
+                output.back().append(b);
+            } else {
+                Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&*cit);
+                if (curve && curve->order() == 3) {
+                    Geom::CubicBezier b((*curve)[0], (*curve)[1], (*curve)[2], (*curve)[3]);
+                    output.back().append(b);
+                } else {
+                    // convert all other curve types to cubicbeziers
+                    Geom::Path cubicbezier_path = Geom::cubicbezierpath_from_sbasis(cit->toSBasis(), 0.1);
+                    output.back().append(cubicbezier_path);
+                }
+            }
+        }
+    }
+
+    return output;
+}
+
+//Study move to 2Geom
+size_t 
+count_pathvector_nodes(Geom::PathVector const &pathv) {
+    size_t tot = 0;
+    for (auto subpath : pathv) {
+        tot += count_path_nodes(subpath);
+    }
+    return tot;
+}
+size_t count_path_nodes(Geom::Path const &path)
+{
+    size_t tot = path.size_closed();
+    if (path.closed()) {
+        const Geom::Curve &closingline = path.back_closed();
+        // the closing line segment is always of type
+        // Geom::LineSegment.
+        if (are_near(closingline.initialPoint(), closingline.finalPoint())) {
+            // closingline.isDegenerate() did not work, because it only checks for
+            // *exact* zero length, which goes wrong for relative coordinates and
+            // rounding errors...
+            // the closing line segment has zero-length. So stop before that one!
+            tot -= 1;
+        }
+    }
+    return tot;
+}
+
+// The next routine is modified from curv4_div::recursive_bezier from file agg_curves.cpp
+//----------------------------------------------------------------------------
+// Anti-Grain Geometry (AGG) - Version 2.5
+// A high quality rendering engine for C++
+// Copyright (C) 2002-2006 Maxim Shemanarev
+// Contact: mcseem@antigrain.com
+//          mcseemagg@yahoo.com
+//          http://antigrain.com
+// 
+// AGG is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 2
+// of the License, or (at your option) any later version.
+// 
+// AGG is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with AGG; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 
+// MA 02110-1301, USA.
+//----------------------------------------------------------------------------
+void
+recursive_bezier4(const double x1, const double y1, 
+                  const double x2, const double y2, 
+                  const double x3, const double y3, 
+                  const double x4, const double y4,
+                  std::vector<Geom::Point> &m_points,
+                  int level)
+    {
+        // some of these should be parameters, but do it this way for now.
+        const double curve_collinearity_epsilon              = 1e-30;
+        const double curve_angle_tolerance_epsilon           = 0.01;
+        double       m_cusp_limit                            = 0.0;
+        double       m_angle_tolerance                       = 0.0;
+        double       m_approximation_scale                   = 1.0;
+        double       m_distance_tolerance_square = 0.5 / m_approximation_scale;
+        m_distance_tolerance_square *= m_distance_tolerance_square;
+        enum curve_recursion_limit_e { curve_recursion_limit = 32 };
+#define calc_sq_distance(A,B,C,D) ((A-C)*(A-C) + (B-D)*(B-D))
+
+        if(level > curve_recursion_limit) 
+        {
+            return;
+        }
+
+
+        // Calculate all the mid-points of the line segments
+        //----------------------
+        double x12   = (x1 + x2) / 2;
+        double y12   = (y1 + y2) / 2;
+        double x23   = (x2 + x3) / 2;
+        double y23   = (y2 + y3) / 2;
+        double x34   = (x3 + x4) / 2;
+        double y34   = (y3 + y4) / 2;
+        double x123  = (x12 + x23) / 2;
+        double y123  = (y12 + y23) / 2;
+        double x234  = (x23 + x34) / 2;
+        double y234  = (y23 + y34) / 2;
+        double x1234 = (x123 + x234) / 2;
+        double y1234 = (y123 + y234) / 2;
+
+
+        // Try to approximate the full cubic curve by a single straight line
+        //------------------
+        double dx = x4-x1;
+        double dy = y4-y1;
+
+        double d2 = fabs(((x2 - x4) * dy - (y2 - y4) * dx));
+        double d3 = fabs(((x3 - x4) * dy - (y3 - y4) * dx));
+        double da1, da2, k;
+
+        switch((int(d2 > curve_collinearity_epsilon) << 1) +
+                int(d3 > curve_collinearity_epsilon))
+        {
+        case 0:
+            // All collinear OR p1==p4
+            //----------------------
+            k = dx*dx + dy*dy;
+            if(k == 0)
+            {
+                d2 = calc_sq_distance(x1, y1, x2, y2);
+                d3 = calc_sq_distance(x4, y4, x3, y3);
+            }
+            else
+            {
+                k   = 1 / k;
+                da1 = x2 - x1;
+                da2 = y2 - y1;
+                d2  = k * (da1*dx + da2*dy);
+                da1 = x3 - x1;
+                da2 = y3 - y1;
+                d3  = k * (da1*dx + da2*dy);
+                if(d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1)
+                {
+                    // Simple collinear case, 1---2---3---4
+                    // We can leave just two endpoints
+                    return;
+                }
+                     if(d2 <= 0) d2 = calc_sq_distance(x2, y2, x1, y1);
+                else if(d2 >= 1) d2 = calc_sq_distance(x2, y2, x4, y4);
+                else             d2 = calc_sq_distance(x2, y2, x1 + d2*dx, y1 + d2*dy);
+
+                     if(d3 <= 0) d3 = calc_sq_distance(x3, y3, x1, y1);
+                else if(d3 >= 1) d3 = calc_sq_distance(x3, y3, x4, y4);
+                else             d3 = calc_sq_distance(x3, y3, x1 + d3*dx, y1 + d3*dy);
+            }
+            if(d2 > d3)
+            {
+                if(d2 < m_distance_tolerance_square)
+                {
+                    m_points.emplace_back(x2, y2);
+                    return;
+                }
+            }
+            else
+            {
+                if(d3 < m_distance_tolerance_square)
+                {
+                    m_points.emplace_back(x3, y3);
+                    return;
+                }
+            }
+            break;
+
+        case 1:
+            // p1,p2,p4 are collinear, p3 is significant
+            //----------------------
+            if(d3 * d3 <= m_distance_tolerance_square * (dx*dx + dy*dy))
+            {
+                if(m_angle_tolerance < curve_angle_tolerance_epsilon)
+                {
+                    m_points.emplace_back(x23, y23);
+                    return;
+                }
+
+                // Angle Condition
+                //----------------------
+                da1 = fabs(atan2(y4 - y3, x4 - x3) - atan2(y3 - y2, x3 - x2));
+                if(da1 >= M_PI) da1 = 2*M_PI - da1;
+
+                if(da1 < m_angle_tolerance)
+                {
+                    m_points.emplace_back(x2, y2);
+                    m_points.emplace_back(x3, y3);
+                    return;
+                }
+
+                if(m_cusp_limit != 0.0)
+                {
+                    if(da1 > m_cusp_limit)
+                    {
+                        m_points.emplace_back(x3, y3);
+                        return;
+                    }
+                }
+            }
+            break;
+
+        case 2:
+            // p1,p3,p4 are collinear, p2 is significant
+            //----------------------
+            if(d2 * d2 <= m_distance_tolerance_square * (dx*dx + dy*dy))
+            {
+                if(m_angle_tolerance < curve_angle_tolerance_epsilon)
+                {
+                    m_points.emplace_back(x23, y23);
+                    return;
+                }
+
+                // Angle Condition
+                //----------------------
+                da1 = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1));
+                if(da1 >= M_PI) da1 = 2*M_PI - da1;
+
+                if(da1 < m_angle_tolerance)
+                {
+                    m_points.emplace_back(x2, y2);
+                    m_points.emplace_back(x3, y3);
+                    return;
+                }
+
+                if(m_cusp_limit != 0.0)
+                {
+                    if(da1 > m_cusp_limit)
+                    {
+                        m_points.emplace_back(x2, y2);
+                        return;
+                    }
+                }
+            }
+            break;
+
+        case 3: 
+            // Regular case
+            //-----------------
+            if((d2 + d3)*(d2 + d3) <= m_distance_tolerance_square * (dx*dx + dy*dy))
+            {
+                // If the curvature doesn't exceed the distance_tolerance value
+                // we tend to finish subdivisions.
+                //----------------------
+                if(m_angle_tolerance < curve_angle_tolerance_epsilon)
+                {
+                    m_points.emplace_back(x23, y23);
+                    return;
+                }
+
+                // Angle & Cusp Condition
+                //----------------------
+                k   = atan2(y3 - y2, x3 - x2);
+                da1 = fabs(k - atan2(y2 - y1, x2 - x1));
+                da2 = fabs(atan2(y4 - y3, x4 - x3) - k);
+                if(da1 >= M_PI) da1 = 2*M_PI - da1;
+                if(da2 >= M_PI) da2 = 2*M_PI - da2;
+
+                if(da1 + da2 < m_angle_tolerance)
+                {
+                    // Finally we can stop the recursion
+                    //----------------------
+                    m_points.emplace_back(x23, y23);
+                    return;
+                }
+
+                if(m_cusp_limit != 0.0)
+                {
+                    if(da1 > m_cusp_limit)
+                    {
+                        m_points.emplace_back(x2, y2);
+                        return;
+                    }
+
+                    if(da2 > m_cusp_limit)
+                    {
+                        m_points.emplace_back(x3, y3);
+                        return;
+                    }
+                }
+            }
+            break;
+        }
+
+        // Continue subdivision
+        //----------------------
+        recursive_bezier4(x1, y1, x12, y12, x123, y123, x1234, y1234, m_points, level + 1); 
+        recursive_bezier4(x1234, y1234, x234, y234, x34, y34, x4, y4, m_points, level + 1); 
+}
+
+void 
+swap(Geom::Point &A, Geom::Point &B){
+    Geom::Point tmp = A;
+    A = B;
+    B = tmp;
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/geom.h b/src/helper/geom.h
new file mode 100644
index 0000000..ba12404
--- /dev/null
+++ b/src/helper/geom.h
@@ -0,0 +1,54 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef INKSCAPE_HELPER_GEOM_H
+#define INKSCAPE_HELPER_GEOM_H
+
+/**
+ * @file
+ * Specific geometry functions for Inkscape, not provided my lib2geom.
+ */
+/*
+ * Author:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2008 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/forward.h>
+#include <2geom/rect.h>
+#include <2geom/affine.h>
+
+Geom::OptRect bounds_fast_transformed(Geom::PathVector const & pv, Geom::Affine const & t);
+Geom::OptRect bounds_exact_transformed(Geom::PathVector const & pv, Geom::Affine const & t);
+
+void pathv_matrix_point_bbox_wind_distance ( Geom::PathVector const & pathv, Geom::Affine const &m, Geom::Point const &pt,
+                                             Geom::Rect *bbox, int *wind, Geom::Coord *dist,
+                                             Geom::Coord tolerance, Geom::Rect const *viewbox);
+
+bool is_intersecting(Geom::PathVector const&a, Geom::PathVector const&b);
+
+size_t count_pathvector_nodes(Geom::PathVector const &pathv );
+size_t count_path_nodes(Geom::Path const &path);
+bool pointInTriangle(Geom::Point const &p, Geom::Point const &p1, Geom::Point const &p2, Geom::Point const &p3);
+Geom::PathVector pathv_to_linear_and_cubic_beziers( Geom::PathVector const &pathv );
+Geom::PathVector pathv_to_linear( Geom::PathVector const &pathv, double maxdisp );
+Geom::PathVector pathv_to_cubicbezier( Geom::PathVector const &pathv);
+bool pathv_similar(Geom::PathVector const &apv, Geom::PathVector const &bpv, double precission = 0.001);
+void recursive_bezier4(const double x1, const double y1, const double x2, const double y2, 
+                       const double x3, const double y3, const double x4, const double y4,
+                       std::vector<Geom::Point> &pointlist,
+                       int level);
+void swap(Geom::Point &A, Geom::Point &B);
+#endif  // INKSCAPE_HELPER_GEOM_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/gettext.cpp b/src/helper/gettext.cpp
new file mode 100644
index 0000000..d3f8cc2
--- /dev/null
+++ b/src/helper/gettext.cpp
@@ -0,0 +1,89 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * helper functions for gettext
+ *//*
+ * Authors:
+ * see git history
+ *   Patrick Storz <eduard.braun2@gmx.de>
+ *
+ * Copyright (C) 2018 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"  // only include where actually required!
+#endif
+
+#ifdef _WIN32
+#include <windows.h>
+#endif
+
+#include "path-prefix.h"
+
+#include <string>
+#include <glibmm.h>
+#include <glibmm/i18n.h>
+
+namespace Inkscape {
+
+/** does all required gettext initialization and takes care of the respective locale directory paths */
+void initialize_gettext() {
+    auto localepath = Glib::getenv("INKSCAPE_LOCALEDIR");
+
+    if (localepath.empty()) {
+        localepath = Glib::build_filename(Glib::path_get_dirname(get_inkscape_datadir()), PACKAGE_LOCALE_DIR);
+    }
+
+    if (!Glib::file_test(localepath, Glib::FILE_TEST_IS_DIR)) {
+        localepath = PACKAGE_LOCALE_DIR_ABSOLUTE;
+    }
+
+#ifdef _WIN32
+    // obtain short path, bindtextdomain doesn't understand UTF-8
+    auto shortlocalepath = g_win32_locale_filename_from_utf8(localepath.c_str());
+    localepath = shortlocalepath;
+    g_free(shortlocalepath);
+#endif
+
+    bindtextdomain(GETTEXT_PACKAGE, localepath.c_str());
+
+    // common setup
+    bind_textdomain_codeset(GETTEXT_PACKAGE, "UTF-8");
+    textdomain(GETTEXT_PACKAGE);
+}
+
+/** set gettext codeset to UTF8 */
+void bind_textdomain_codeset_utf8() {
+    bind_textdomain_codeset(GETTEXT_PACKAGE, "UTF-8");
+}
+
+/** set gettext codeset to codeset of the system console
+ *   - on *nix this is typically the current locale,
+ *   - on Windows we don't care and simply use UTF8
+ *       any conversion would need to happen in our console wrappers (see winconsole.cpp) anyway
+ *       as we have no easy way of determining console encoding from inkscape/inkview.exe process;
+ *       for now do something even easier - switch console encoding to UTF8 and be done with it!
+ *     this also works nicely on MSYS consoles where UTF8 encoding is used by default, too  */
+void bind_textdomain_codeset_console() {
+#ifdef _WIN32
+    bind_textdomain_codeset(GETTEXT_PACKAGE, "UTF-8");
+#else
+    std::string charset;
+    Glib::get_charset(charset);
+    bind_textdomain_codeset(GETTEXT_PACKAGE, charset.c_str());
+#endif
+}
+
+}
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace .0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim:filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99:
diff --git a/src/helper/gettext.h b/src/helper/gettext.h
new file mode 100644
index 0000000..e5745ee
--- /dev/null
+++ b/src/helper/gettext.h
@@ -0,0 +1,33 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * \brief helper functions for gettext
+ *//*
+ * Authors:
+ * see git history
+ *   Patrick Storz <eduard.braun2@gmx.de>
+ *
+ * Copyright (C) 2017 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifndef SEEN_GETTEXT_HELPER_H
+#define SEEN_GETTEXT_HELPER_H
+
+namespace Inkscape {
+     void initialize_gettext();
+     void bind_textdomain_codeset_utf8();
+     void bind_textdomain_codeset_console();
+}
+
+#endif // SEEN_GETTEXT_HELPER_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace .0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim:filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99:
diff --git a/src/helper/mathfns.h b/src/helper/mathfns.h
new file mode 100644
index 0000000..e0e321c
--- /dev/null
+++ b/src/helper/mathfns.h
@@ -0,0 +1,83 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * ... some mathematical functions
+ *
+ * Authors:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *
+ * Copyright (C) 2007 Johan Engelen
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#ifndef SEEN_INKSCAPE_UTIL_MATHFNS_H
+#define SEEN_INKSCAPE_UTIL_MATHFNS_H
+
+#include <2geom/point.h>
+
+namespace Inkscape {
+
+namespace Util {
+
+/**
+ * Returns area in triangle given by points; may be negative.
+ */
+inline double
+triangle_area (Geom::Point p1, Geom::Point p2, Geom::Point p3)
+{
+    using Geom::X;
+    using Geom::Y;
+    return (p1[X]*p2[Y] + p1[Y]*p3[X] + p2[X]*p3[Y] - p2[Y]*p3[X] - p1[Y]*p2[X] - p1[X]*p3[Y]);
+}
+
+/**
+ * \return x rounded to the nearest multiple of c1 plus c0.
+ *
+ * \note
+ * If c1==0 (and c0 is finite), then returns +/-inf.  This makes grid spacing of zero
+ * mean "ignore the grid in this dimension".
+ */
+inline double round_to_nearest_multiple_plus(double x, double const c1, double const c0)
+{
+    return floor((x - c0) / c1 + .5) * c1 + c0;
+}
+
+/**
+ * \return x rounded to the lower multiple of c1 plus c0.
+ *
+ * \note
+ * If c1==0 (and c0 is finite), then returns +/-inf.  This makes grid spacing of zero
+ * mean "ignore the grid in this dimension".
+ */
+inline double round_to_lower_multiple_plus(double x, double const c1, double const c0 = 0)
+{
+    return floor((x - c0) / c1) * c1 + c0;
+}
+
+/**
+ * \return x rounded to the upper multiple of c1 plus c0.
+ *
+ * \note
+ * If c1==0 (and c0 is finite), then returns +/-inf.  This makes grid spacing of zero
+ * mean "ignore the grid in this dimension".
+ */
+inline double round_to_upper_multiple_plus(double x, double const c1, double const c0 = 0)
+{
+    return ceil((x - c0) / c1) * c1 + c0;
+}
+
+}
+
+}
+
+#endif
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/pixbuf-ops.cpp b/src/helper/pixbuf-ops.cpp
new file mode 100644
index 0000000..dbe9ff6
--- /dev/null
+++ b/src/helper/pixbuf-ops.cpp
@@ -0,0 +1,157 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * Helpers for SPItem -> gdk_pixbuf related stuff
+ *
+ * Authors:
+ *   John Cliff <simarilius@yahoo.com>
+ *   Jon A. Cruz <jon@joncruz.org>
+ *   Abhishek Sharma
+ *
+ * Copyright (C) 2008 John Cliff
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <2geom/transforms.h>
+
+#include "helper/png-write.h"
+#include "display/cairo-utils.h"
+#include "display/drawing.h"
+#include "display/drawing-context.h"
+#include "document.h"
+#include "object/sp-root.h"
+#include "object/sp-defs.h"
+#include "object/sp-use.h"
+#include "util/units.h"
+#include "inkscape.h"
+
+#include "helper/pixbuf-ops.h"
+
+#include <gdk/gdk.h>
+
+// TODO look for copy-n-paste duplication of this function:
+/**
+ * Hide all items that are not in list, recursively, skipping groups and defs.
+ */
+static void hide_other_items_recursively(SPObject *object, const std::vector<SPItem*> &items, unsigned dkey)
+{
+    SPItem *item = dynamic_cast<SPItem *>(object);
+    if (!item) {
+        // <defs>, <metadata>, etc.
+        return;
+    }
+
+    if (std::find (items.begin(), items.end(), item) != items.end()) {
+        // It's in our list, don't hide!
+        return;
+    }
+
+    if ( !dynamic_cast<SPRoot  *>(item) &&
+         !dynamic_cast<SPGroup *>(item) &&
+         !dynamic_cast<SPUse   *>(item) ) {
+        // Hide if not container or def.
+        item->invoke_hide(dkey);
+    }
+
+    for (auto& child: object->children) {
+        hide_other_items_recursively(&child, items, dkey);
+    }
+}
+
+
+/**
+    generates a bitmap from given items
+    the bitmap is stored in RAM and not written to file
+    @param document Inkscape document.
+    @param area     Export area in document units.
+    @param dpi      Resolution.
+    @param items    Vector of pointers to SPItems to export. Export all items if empty.
+    @param opaque   Set items opacity to 1 (used by Cairo renderer for filtered objects rendered as bitmaps).
+    @return The created GdkPixbuf structure or nullptr if rendering failed.
+*/
+Inkscape::Pixbuf *sp_generate_internal_bitmap(SPDocument *document,
+                                              Geom::Rect const &area,
+                                              double dpi,
+                                              std::vector<SPItem *> items,
+                                              bool opaque)
+{
+    // Geometry
+    if (area.hasZeroArea()) {
+        return nullptr;
+    }
+
+    Geom::Point origin = area.min();
+    double scale_factor = Inkscape::Util::Quantity::convert(dpi, "px", "in");
+    Geom::Affine affine = Geom::Translate(-origin) * Geom::Scale (scale_factor, scale_factor);
+
+    int width  = std::ceil(scale_factor * area.width());
+    int height = std::ceil(scale_factor * area.height());
+
+    // Document
+    document->ensureUpToDate();
+    unsigned dkey = SPItem::display_key_new(1);
+
+    // Drawing
+    Inkscape::Drawing drawing; // New drawing for offscreen rendering.
+    drawing.setExact(true); // Maximum quality for blurs.
+
+    /* Create ArenaItems and set transform */
+    Inkscape::DrawingItem *root = document->getRoot()->invoke_show( drawing, dkey, SP_ITEM_SHOW_DISPLAY);
+    root->setTransform(affine);
+    drawing.setRoot(root);
+
+    // Hide all items we don't want, instead of showing only requested items,
+    // because that would not work if the shown item references something in defs.
+    if (!items.empty()) {
+        hide_other_items_recursively(document->getRoot(), items, dkey);
+    }
+
+    Geom::IntRect final_area = Geom::IntRect::from_xywh(0, 0, width, height);
+    drawing.update(final_area);
+
+    if (opaque) {
+        // Required by sp_asbitmap_render().
+        for (auto item : items) {
+            if (item->get_arenaitem(dkey)) {
+                item->get_arenaitem(dkey)->setOpacity(1.0);
+            }
+        }
+    }
+
+    // Rendering
+    cairo_surface_t *surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, width, height);
+    Inkscape::Pixbuf* pixbuf = nullptr;
+
+    if (cairo_surface_status(surface) == CAIRO_STATUS_SUCCESS) {
+        Inkscape::DrawingContext dc(surface, Geom::Point(0,0));
+
+        // render items
+        drawing.render(dc, final_area, Inkscape::DrawingItem::RENDER_BYPASS_CACHE);
+
+        pixbuf = new Inkscape::Pixbuf(surface);
+
+    } else {
+
+        long long size =
+            (long long) height *
+            (long long) cairo_format_stride_for_width(CAIRO_FORMAT_ARGB32, width);
+        g_warning("sp_generate_internal_bitmap: not enough memory to create pixel buffer. Need %lld.", size);
+        cairo_surface_destroy(surface);
+    }
+
+    // Return to previous state.
+    document->getRoot()->invoke_hide(dkey);
+
+    return pixbuf;
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/pixbuf-ops.h b/src/helper/pixbuf-ops.h
new file mode 100644
index 0000000..91cd73b
--- /dev/null
+++ b/src/helper/pixbuf-ops.h
@@ -0,0 +1,27 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef __SP_PIXBUF_OPS_H__
+#define __SP_PIXBUF_OPS_H__
+
+/*
+ * Helpers for SPItem -> gdk_pixbuf related stuff
+ *
+ * Authors:
+ *   John Cliff <simarilius@yahoo.com>
+ *
+ * Copyright (C) 2008 John Cliff
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <glib.h>
+
+class SPDocument;
+class SPItem;
+namespace Inkscape { class Pixbuf; }
+
+Inkscape::Pixbuf *sp_generate_internal_bitmap(SPDocument *document,
+                                              Geom::Rect const &area,
+                                              double dpi,
+                                              std::vector<SPItem *> items = std::vector<SPItem*>(),
+                                              bool set_opaque = false);
+#endif
diff --git a/src/helper/png-write.cpp b/src/helper/png-write.cpp
new file mode 100644
index 0000000..88ca1a5
--- /dev/null
+++ b/src/helper/png-write.cpp
@@ -0,0 +1,508 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * PNG file format utilities
+ *
+ * Authors:
+ *   Lauris Kaplinski <lauris@kaplinski.com>
+ *   Whoever wrote this example in libpng documentation
+ *   Peter Bostrom
+ *   Jon A. Cruz <jon@joncruz.org>
+ *   Abhishek Sharma
+ *
+ * Copyright (C) 1999-2002 authors
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+
+#include <2geom/rect.h>
+#include <2geom/transforms.h>
+
+#include <png.h>
+
+#include "document.h"
+#include "inkscape.h"
+#include "png-write.h"
+#include "preferences.h"
+#include "rdf.h"
+
+#include "display/cairo-utils.h"
+#include "display/drawing-context.h"
+#include "display/drawing.h"
+
+#include "io/sys.h"
+
+#include "object/sp-defs.h"
+#include "object/sp-item.h"
+#include "object/sp-root.h"
+
+#include "ui/interface.h"
+#include "util/units.h"
+
+/* This is an example of how to use libpng to read and write PNG files.
+ * The file libpng.txt is much more verbose then this.  If you have not
+ * read it, do so first.  This was designed to be a starting point of an
+ * implementation.  This is not officially part of libpng, and therefore
+ * does not require a copyright notice.
+ *
+ * This file does not currently compile, because it is missing certain
+ * parts, like allocating memory to hold an image.  You will have to
+ * supply these parts to get it to compile.  For an example of a minimal
+ * working PNG reader/writer, see pngtest.c, included in this distribution.
+ */
+
+struct SPEBP {
+    unsigned long int width, height, sheight;
+    guint32 background;
+    Inkscape::Drawing *drawing; // it is assumed that all unneeded items are hidden
+    guchar *px;
+    unsigned (*status)(float, void *);
+    void *data;
+};
+
+/* write a png file */
+
+struct SPPNGBD {
+    guchar const *px;
+    int rowstride;
+};
+
+/**
+ * A simple wrapper to list png_text.
+ */
+class PngTextList {
+public:
+    PngTextList() : count(0), textItems(nullptr) {}
+    ~PngTextList();
+
+    void add(gchar const* key, gchar const* text);
+    gint getCount() {return count;}
+    png_text* getPtext() {return textItems;}
+
+private:
+    gint count;
+    png_text* textItems;
+};
+
+PngTextList::~PngTextList() {
+    for (gint i = 0; i < count; i++) {
+        if (textItems[i].key) {
+            g_free(textItems[i].key);
+        }
+        if (textItems[i].text) {
+            g_free(textItems[i].text);
+        }
+    }
+}
+
+void PngTextList::add(gchar const* key, gchar const* text)
+{
+    if (count < 0) {
+        count = 0;
+        textItems = nullptr;
+    }
+    png_text* tmp = (count > 0) ? g_try_renew(png_text, textItems, count + 1): g_try_new(png_text, 1);
+    if (tmp) {
+        textItems = tmp;
+        count++;
+
+        png_text* item = &(textItems[count - 1]);
+        item->compression = PNG_TEXT_COMPRESSION_NONE;
+        item->key = g_strdup(key);
+        item->text = g_strdup(text);
+        item->text_length = 0;
+#ifdef PNG_iTXt_SUPPORTED
+        item->itxt_length = 0;
+        item->lang = nullptr;
+        item->lang_key = nullptr;
+#endif // PNG_iTXt_SUPPORTED
+    } else {
+        g_warning("Unable to allocate array for %d PNG text data.", count);
+        textItems = nullptr;
+        count = 0;
+    }
+}
+
+static bool
+sp_png_write_rgba_striped(SPDocument *doc,
+                          gchar const *filename, unsigned long int width, unsigned long int height, double xdpi, double ydpi,
+                          int (* get_rows)(guchar const **rows, void **to_free, int row, int num_rows, void *data, int color_type, int bit_depth, int antialias),
+                          void *data, bool interlace, int color_type, int bit_depth, int zlib, int antialiasing)
+{
+    g_return_val_if_fail(filename != nullptr, false);
+    g_return_val_if_fail(data != nullptr, false);
+
+    struct SPEBP *ebp = (struct SPEBP *) data;
+    FILE *fp;
+    png_structp png_ptr;
+    png_infop info_ptr;
+    png_color_8 sig_bit;
+    png_uint_32 r;
+
+    /* open the file */
+
+    Inkscape::IO::dump_fopen_call(filename, "M");
+    fp = Inkscape::IO::fopen_utf8name(filename, "wb");
+    if(fp == nullptr) return false;
+
+    /* Create and initialize the png_struct with the desired error handler
+     * functions.  If you want to use the default stderr and longjump method,
+     * you can supply NULL for the last three parameters.  We also check that
+     * the library version is compatible with the one used at compile time,
+     * in case we are using dynamically linked libraries.  REQUIRED.
+     */
+    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
+
+    if (png_ptr == nullptr) {
+        fclose(fp);
+        return false;
+    }
+
+    /* Allocate/initialize the image information data.  REQUIRED */
+    info_ptr = png_create_info_struct(png_ptr);
+    if (info_ptr == nullptr) {
+        fclose(fp);
+        png_destroy_write_struct(&png_ptr, nullptr);
+        return false;
+    }
+
+    /* Set error handling.  REQUIRED if you aren't supplying your own
+     * error handling functions in the png_create_write_struct() call.
+     */
+    if (setjmp(png_jmpbuf(png_ptr))) {
+        // If we get here, we had a problem reading the file
+        fclose(fp);
+        png_destroy_write_struct(&png_ptr, &info_ptr);
+        return false;
+    }
+
+    /* set up the output control if you are using standard C streams */
+    png_init_io(png_ptr, fp);
+
+    /* Set the image information here.  Width and height are up to 2^31,
+     * bit_depth is one of 1, 2, 4, 8, or 16, but valid values also depend on
+     * the color_type selected. color_type is one of PNG_COLOR_TYPE_GRAY,
+     * PNG_COLOR_TYPE_GRAY_ALPHA, PNG_COLOR_TYPE_PALETTE, PNG_COLOR_TYPE_RGB,
+     * or PNG_COLOR_TYPE_RGB_ALPHA.  interlace is either PNG_INTERLACE_NONE or
+     * PNG_INTERLACE_ADAM7, and the compression_type and filter_type MUST
+     * currently be PNG_COMPRESSION_TYPE_BASE and PNG_FILTER_TYPE_BASE. REQUIRED
+     */
+
+    png_set_compression_level(png_ptr, zlib);
+
+    png_set_IHDR(png_ptr, info_ptr,
+                 width,
+                 height,
+                 bit_depth,
+                 color_type,
+                 interlace ? PNG_INTERLACE_ADAM7 : PNG_INTERLACE_NONE,
+                 PNG_COMPRESSION_TYPE_BASE,
+                 PNG_FILTER_TYPE_BASE);
+
+    if ((color_type&2) && bit_depth == 16) {
+        // otherwise, if we are dealing with a color image then
+        sig_bit.red = 8;
+        sig_bit.green = 8;
+        sig_bit.blue = 8;
+        // if the image has an alpha channel then
+        if (color_type&4)
+            sig_bit.alpha = 8;
+        png_set_sBIT(png_ptr, info_ptr, &sig_bit);
+    }
+
+    PngTextList textList;
+
+    textList.add("Software", "www.inkscape.org"); // Made by Inkscape comment
+    {
+        const gchar* pngToDc[] = {"Title", "title",
+                               "Author", "creator",
+                               "Description", "description",
+                               //"Copyright", "",
+                               "Creation Time", "date",
+                               //"Disclaimer", "",
+                               //"Warning", "",
+                               "Source", "source"
+                               //"Comment", ""
+        };
+        for (size_t i = 0; i < G_N_ELEMENTS(pngToDc); i += 2) {
+            struct rdf_work_entity_t * entity = rdf_find_entity ( pngToDc[i + 1] );
+            if (entity) {
+                gchar const* data = rdf_get_work_entity(doc, entity);
+                if (data && *data) {
+                    textList.add(pngToDc[i], data);
+                }
+            } else {
+                g_warning("Unable to find entity [%s]", pngToDc[i + 1]);
+            }
+        }
+
+
+        struct rdf_license_t *license =  rdf_get_license(doc);
+        if (license) {
+            if (license->name && license->uri) {
+                gchar* tmp = g_strdup_printf("%s %s", license->name, license->uri);
+                textList.add("Copyright", tmp);
+                g_free(tmp);
+            } else if (license->name) {
+                textList.add("Copyright", license->name);
+            } else if (license->uri) {
+                textList.add("Copyright", license->uri);
+            }
+        }
+    }
+    if (textList.getCount() > 0) {
+        png_set_text(png_ptr, info_ptr, textList.getPtext(), textList.getCount());
+    }
+
+    /* other optional chunks like cHRM, bKGD, tRNS, tIME, oFFs, pHYs, */
+    /* note that if sRGB is present the cHRM chunk must be ignored
+     * on read and must be written in accordance with the sRGB profile */
+    if(xdpi < 0.0254 ) xdpi = 0.0255;
+    if(ydpi < 0.0254 ) ydpi = 0.0255;
+
+    png_set_pHYs(png_ptr, info_ptr, unsigned(xdpi / 0.0254 ), unsigned(ydpi / 0.0254 ), PNG_RESOLUTION_METER);
+
+    /* Write the file header information.  REQUIRED */
+    png_write_info(png_ptr, info_ptr);
+
+    /* Once we write out the header, the compression type on the text
+     * chunks gets changed to PNG_TEXT_COMPRESSION_NONE_WR or
+     * PNG_TEXT_COMPRESSION_zTXt_WR, so it doesn't get written out again
+     * at the end.
+     */
+
+    /* set up the transformations you want.  Note that these are
+     * all optional.  Only call them if you want them.
+     */
+
+    /* --- CUT --- */
+
+    /* The easiest way to write the image (you may have a different memory
+     * layout, however, so choose what fits your needs best).  You need to
+     * use the first method if you aren't handling interlacing yourself.
+     */
+
+    png_bytep* row_pointers = new png_bytep[ebp->sheight];
+    int number_of_passes = interlace ? png_set_interlace_handling(png_ptr) : 1;
+
+    for(int i=0;i<number_of_passes; ++i){
+        r = 0;
+        while (r < static_cast<png_uint_32>(height)) {
+            void *to_free;
+            int n = get_rows((unsigned char const **) row_pointers, &to_free, r, height-r, data, color_type, bit_depth, antialiasing);
+            if (!n) break;
+            png_write_rows(png_ptr, row_pointers, n);
+            g_free(to_free);
+            r += n;
+        }
+    }
+
+    delete[] row_pointers;
+
+    /* You can write optional chunks like tEXt, zTXt, and tIME at the end
+     * as well.
+     */
+
+    /* It is REQUIRED to call this to finish writing the rest of the file */
+    png_write_end(png_ptr, info_ptr);
+
+    /* if you allocated any text comments, free them here */
+
+    /* clean up after the write, and free any memory allocated */
+    png_destroy_write_struct(&png_ptr, &info_ptr);
+
+    /* close the file */
+    fclose(fp);
+
+    /* that's it */
+    return true;
+}
+
+
+/**
+ *
+ */
+static int
+sp_export_get_rows(guchar const **rows, void **to_free, int row, int num_rows, void *data, int color_type, int bit_depth, int antialiasing)
+{
+    struct SPEBP *ebp = (struct SPEBP *) data;
+
+    if (ebp->status) {
+        if (!ebp->status((float) row / ebp->height, ebp->data)) return 0;
+    }
+
+    num_rows = MIN(num_rows, static_cast<int>(ebp->sheight));
+    num_rows = MIN(num_rows, static_cast<int>(ebp->height - row));
+
+    /* Set area of interest */
+    // bbox is now set to the entire image to prevent discontinuities
+    // in the image when blur is used (the borders may still be a bit
+    // off, but that's less noticeable).
+    Geom::IntRect bbox = Geom::IntRect::from_xywh(0, row, ebp->width, num_rows);
+
+    /* Update to renderable state */
+    ebp->drawing->update(bbox);
+
+    int stride = cairo_format_stride_for_width(CAIRO_FORMAT_ARGB32, ebp->width);
+    unsigned char *px = g_new(guchar, num_rows * stride);
+
+    cairo_surface_t *s = cairo_image_surface_create_for_data(
+        px, CAIRO_FORMAT_ARGB32, ebp->width, num_rows, stride);
+    Inkscape::DrawingContext dc(s, bbox.min());
+    dc.setSource(ebp->background);
+    dc.setOperator(CAIRO_OPERATOR_SOURCE);
+    dc.paint();
+    dc.setOperator(CAIRO_OPERATOR_OVER);
+
+    /* Render */
+    ebp->drawing->render(dc, bbox, 0, antialiasing);
+    cairo_surface_destroy(s);
+
+    // PNG stores data as unpremultiplied big-endian RGBA, which means
+    // it's identical to the GdkPixbuf format.
+    convert_pixels_argb32_to_pixbuf(px, ebp->width, num_rows, stride,
+                                    /* RGBA to ARGB with A=0 */ ebp->background >> 8);
+    
+    // If a custom bit depth or color type is asked, then convert rgb to grayscale, etc.
+    const guchar* new_data = pixbuf_to_png(rows, px, num_rows, ebp->width, stride, color_type, bit_depth);
+    *to_free = (void*) new_data;
+    free(px);
+
+    return num_rows;
+}
+
+/**
+ * Hide all items that are not listed in list, recursively, skipping groups and defs.
+ */
+static void hide_other_items_recursively(SPObject *o, const std::vector<SPItem*> &list, unsigned dkey)
+{
+    if ( SP_IS_ITEM(o)
+         && !SP_IS_DEFS(o)
+         && !SP_IS_ROOT(o)
+         && !SP_IS_GROUP(o)
+         && list.end()==find(list.begin(),list.end(),o))
+    {
+        SP_ITEM(o)->invoke_hide(dkey);
+    }
+
+    // recurse
+    if (list.end()==find(list.begin(),list.end(),o)) {
+        for (auto& child: o->children) {
+            hide_other_items_recursively(&child, list, dkey);
+        }
+    }
+}
+
+
+ExportResult sp_export_png_file(SPDocument *doc, gchar const *filename,
+                                double x0, double y0, double x1, double y1,
+                                unsigned long int width, unsigned long int height, double xdpi, double ydpi,
+                                unsigned long bgcolor,
+                                unsigned int (*status) (float, void *),
+                                void *data, bool force_overwrite,
+                                const std::vector<SPItem*> &items_only, bool interlace, int color_type, int bit_depth, int zlib, int antialiasing)
+{
+    return sp_export_png_file(doc, filename, Geom::Rect(Geom::Point(x0,y0),Geom::Point(x1,y1)),
+                              width, height, xdpi, ydpi, bgcolor, status, data, force_overwrite, items_only, interlace, color_type, bit_depth, zlib, antialiasing);
+}
+
+/**
+ * Export an area to a PNG file
+ *
+ * @param area Area in document coordinates
+ */
+ExportResult sp_export_png_file(SPDocument *doc, gchar const *filename,
+                                Geom::Rect const &area,
+                                unsigned long width, unsigned long height, double xdpi, double ydpi,
+                                unsigned long bgcolor,
+                                unsigned (*status)(float, void *),
+                                void *data, bool force_overwrite,
+                                const std::vector<SPItem*> &items_only, bool interlace, int color_type, int bit_depth, int zlib, int antialiasing)
+{
+    g_return_val_if_fail(doc != nullptr, EXPORT_ERROR);
+    g_return_val_if_fail(filename != nullptr, EXPORT_ERROR);
+    g_return_val_if_fail(width >= 1, EXPORT_ERROR);
+    g_return_val_if_fail(height >= 1, EXPORT_ERROR);
+    g_return_val_if_fail(!area.hasZeroArea(), EXPORT_ERROR);
+
+    if (!force_overwrite && !sp_ui_overwrite_file(filename)) {
+        // aborted overwrite
+	return EXPORT_ABORTED;
+    }
+
+    doc->ensureUpToDate();
+
+    /* Calculate translation by transforming to document coordinates (flipping Y)*/
+    Geom::Point translation = -area.min();
+
+    /*  This calculation is only valid when assumed that (x0,y0)= area.corner(0) and (x1,y1) = area.corner(2)
+     * 1) a[0] * x0 + a[2] * y1 + a[4] = 0.0
+     * 2) a[1] * x0 + a[3] * y1 + a[5] = 0.0
+     * 3) a[0] * x1 + a[2] * y1 + a[4] = width
+     * 4) a[1] * x0 + a[3] * y0 + a[5] = height
+     * 5) a[1] = 0.0;
+     * 6) a[2] = 0.0;
+     *
+     * (1,3) a[0] * x1 - a[0] * x0 = width
+     * a[0] = width / (x1 - x0)
+     * (2,4) a[3] * y0 - a[3] * y1 = height
+     * a[3] = height / (y0 - y1)
+     * (1) a[4] = -a[0] * x0
+     * (2) a[5] = -a[3] * y1
+     */
+
+    Geom::Affine const affine(Geom::Translate(translation)
+                            * Geom::Scale(width / area.width(),
+                                        height / area.height()));
+
+    struct SPEBP ebp;
+    ebp.width  = width;
+    ebp.height = height;
+    ebp.background = bgcolor;
+
+    /* Create new drawing */
+    Inkscape::Drawing drawing;
+    drawing.setExact(true); // export with maximum blur rendering quality
+    unsigned const dkey = SPItem::display_key_new(1);
+
+    // Create ArenaItems and set transform
+    drawing.setRoot(doc->getRoot()->invoke_show(drawing, dkey, SP_ITEM_SHOW_DISPLAY));
+    drawing.root()->setTransform(affine);
+    ebp.drawing = &drawing;
+
+    // We show all and then hide all items we don't want, instead of showing only requested items,
+    // because that would not work if the shown item references something in defs
+    if (!items_only.empty()) {
+        hide_other_items_recursively(doc->getRoot(), items_only, dkey);
+    }
+
+    ebp.status = status;
+    ebp.data   = data;
+
+    bool write_status = false;;
+
+    ebp.sheight = 64;
+    ebp.px = g_try_new(guchar, 4 * ebp.sheight * width);
+
+    if (ebp.px) {
+        write_status = sp_png_write_rgba_striped(doc, filename, width, height, xdpi, ydpi, sp_export_get_rows, &ebp, interlace, color_type, bit_depth, zlib, antialiasing);
+        g_free(ebp.px);
+    }
+
+    // Hide items, this releases arenaitem
+    doc->getRoot()->invoke_hide(dkey);
+
+    return write_status ? EXPORT_OK : EXPORT_ERROR;
+}
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/helper/png-write.h b/src/helper/png-write.h
new file mode 100644
index 0000000..6c79d20
--- /dev/null
+++ b/src/helper/png-write.h
@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+#ifndef SEEN_SP_PNG_WRITE_H
+#define SEEN_SP_PNG_WRITE_H
+
+/*
+ * PNG file format utilities
+ *
+ * Authors:
+ *   Lauris Kaplinski <lauris@kaplinski.com>
+ *   Peter Bostrom
+ *   Jon A. Cruz <jon@joncruz.org>
+ *
+ * Copyright (C) 1999-2002 Lauris Kaplinski
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#include <glib.h> // Only for gchar.
+#include <vector>
+
+#include <2geom/forward.h>
+
+class SPDocument;
+class SPItem;
+
+enum ExportResult {
+    EXPORT_ERROR = 0,
+    EXPORT_OK,
+    EXPORT_ABORTED
+};
+
+/**
+ * Export the given document as a Portable Network Graphics (PNG) file.
+ *
+ * @return EXPORT_OK if succeeded, EXPORT_ABORTED if no action was taken, EXPORT_ERROR (false) if an error occurred.
+ */
+ExportResult sp_export_png_file(SPDocument *doc, gchar const *filename,
+				double x0, double y0, double x1, double y1,
+				unsigned long int width, unsigned long int height, double xdpi, double ydpi,
+				unsigned long bgcolor,
+				unsigned int (*status) (float, void *), void *data, bool force_overwrite = false, const std::vector<SPItem*> &items_only = std::vector<SPItem*>(), 
+                                bool interlace = false, int color_type = 6, int bit_depth = 8, int zlib = 6, int antialiasing = 2);
+
+ExportResult sp_export_png_file(SPDocument *doc, gchar const *filename,
+				Geom::Rect const &area,
+				unsigned long int width, unsigned long int height, double xdpi, double ydpi,
+				unsigned long bgcolor,
+				unsigned int (*status) (float, void *), void *data, bool force_overwrite = false, const std::vector<SPItem*> &items_only = std::vector<SPItem*>(), 
+                                bool interlace = false, int color_type = 6, int bit_depth = 8, int zlib = 6, int antialiasing = 2);
+
+#endif // SEEN_SP_PNG_WRITE_H
diff --git a/src/helper/sp-marshal.list b/src/helper/sp-marshal.list
new file mode 100644
index 0000000..340d2c5
--- /dev/null
+++ b/src/helper/sp-marshal.list
@@ -0,0 +1,10 @@
+# SPDX-License-Identifier: GPL-2.0-or-later
+# marshallers for inkscape
+VOID:POINTER,UINT
+BOOLEAN:POINTER
+BOOLEAN:POINTER,UINT
+BOOLEAN:POINTER,POINTER
+INT:POINTER,POINTER
+DOUBLE:POINTER,UINT
+VOID:INT,INT
+VOID:STRING,STRING
diff --git a/src/helper/stock-items.cpp b/src/helper/stock-items.cpp
new file mode 100644
index 0000000..fbdf11c
--- /dev/null
+++ b/src/helper/stock-items.cpp
@@ -0,0 +1,243 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * Stock-items
+ *
+ * Stock Item management code
+ *
+ * Authors:
+ *  John Cliff <simarilius@yahoo.com>
+ *  Jon A. Cruz <jon@joncruz.org>
+ *   Abhishek Sharma
+ *
+ * Copyright 2004 John Cliff
+ *
+ *
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+
+#define noSP_SS_VERBOSE
+
+#include <cstring>
+#include <glibmm/fileutils.h>
+
+#include "path-prefix.h"
+
+#include <xml/repr.h>
+#include "inkscape.h"
+
+#include "io/sys.h"
+#include "io/resource.h"
+#include "stock-items.h"
+
+#include "object/sp-gradient.h"
+#include "object/sp-pattern.h"
+#include "object/sp-marker.h"
+#include "object/sp-defs.h"
+
+static SPDocument *load_paint_doc(char const *basename,
+                                  Inkscape::IO::Resource::Type type = Inkscape::IO::Resource::PAINT)
+{
+    using namespace Inkscape::IO::Resource;
+
+    for (Domain const domain : {SYSTEM, CREATE}) {
+        auto const filename = get_path_string(domain, type, basename);
+        if (Glib::file_test(filename, Glib::FILE_TEST_IS_REGULAR)) {
+            auto doc = SPDocument::createNewDoc(filename.c_str(), false);
+            if (doc) {
+                doc->ensureUpToDate();
+                return doc;
+            }
+        }
+    }
+
+    return nullptr;
+}
+
+// FIXME: these should be merged with the icon loading code so they
+// can share a common file/doc cache.  This function should just
+// take the dir to look in, and the file to check for, and cache
+// against that, rather than the existing copy/paste code seen here.
+
+static SPObject * sp_marker_load_from_svg(gchar const *name, SPDocument *current_doc)
+{
+    if (!current_doc) {
+        return nullptr;
+    }
+    /* Try to load from document */
+    static SPDocument *doc = load_paint_doc("markers.svg", Inkscape::IO::Resource::MARKERS);
+
+    if (doc) {
+        /* Get the marker we want */
+        SPObject *object = doc->getObjectById(name);
+        if (object && SP_IS_MARKER(object)) {
+            SPDefs *defs = current_doc->getDefs();
+            Inkscape::XML::Document *xml_doc = current_doc->getReprDoc();
+            Inkscape::XML::Node *mark_repr = object->getRepr()->duplicate(xml_doc);
+            defs->getRepr()->addChild(mark_repr, nullptr);
+            SPObject *cloned_item = current_doc->getObjectByRepr(mark_repr);
+            Inkscape::GC::release(mark_repr);
+            return cloned_item;
+        }
+    }
+    return nullptr;
+}
+
+
+static SPObject *
+sp_pattern_load_from_svg(gchar const *name, SPDocument *current_doc)
+{
+    if (!current_doc) {
+        return nullptr;
+    }
+    /* Try to load from document */
+    static SPDocument *doc = load_paint_doc("patterns.svg");
+
+    if (doc) {
+        /* Get the pattern we want */
+        SPObject *object = doc->getObjectById(name);
+        if (object && SP_IS_PATTERN(object)) {
+            SPDefs *defs = current_doc->getDefs();
+            Inkscape::XML::Document *xml_doc = current_doc->getReprDoc();
+            Inkscape::XML::Node *pat_repr = object->getRepr()->duplicate(xml_doc);
+            defs->getRepr()->addChild(pat_repr, nullptr);
+            Inkscape::GC::release(pat_repr);
+            return object;
+        }
+    }
+    return nullptr;
+}
+
+
+static SPObject *
+sp_gradient_load_from_svg(gchar const *name, SPDocument *current_doc)
+{
+    if (!current_doc) {
+        return nullptr;
+    }
+    /* Try to load from document */
+    static SPDocument *doc = load_paint_doc("gradients.svg");
+
+    if (doc) {
+        /* Get the gradient we want */
+        SPObject *object = doc->getObjectById(name);
+        if (object && SP_IS_GRADIENT(object)) {
+            SPDefs *defs = current_doc->getDefs();
+            Inkscape::XML::Document *xml_doc = current_doc->getReprDoc();
+            Inkscape::XML::Node *pat_repr = object->getRepr()->duplicate(xml_doc);
+            defs->getRepr()->addChild(pat_repr, nullptr);
+            Inkscape::GC::release(pat_repr);
+            return object;
+        }
+    }
+    return nullptr;
+}
+
+// get_stock_item returns a pointer to an instance of the desired stock object in the current doc
+// if necessary it will import the object. Copes with name clashes through use of the inkscape:stockid property
+// This should be set to be the same as the id in the library file.
+
+SPObject *get_stock_item(gchar const *urn, gboolean stock)
+{
+    g_assert(urn != nullptr);
+
+    /* check its an inkscape URN */
+    if (!strncmp (urn, "urn:inkscape:", 13)) {
+
+        gchar const *e = urn + 13;
+        int a = 0;
+        gchar * name = g_strdup(e);
+        gchar *name_p = name;
+        while (*name_p != ':' && *name_p != '\0'){
+            name_p++;
+            a++;
+        }
+
+        if (*name_p ==':') {
+            name_p++;
+        }
+
+        gchar * base = g_strndup(e, a);
+
+        SPDocument *doc = SP_ACTIVE_DOCUMENT;
+        SPDefs *defs = doc->getDefs();
+        if (!defs) {
+            g_free(base);
+            return nullptr;
+        }
+        SPObject *object = nullptr;
+        if (!strcmp(base, "marker") && !stock) {
+            for (auto& child: defs->children)
+            {
+                if (child.getRepr()->attribute("inkscape:stockid") &&
+                    !strcmp(name_p, child.getRepr()->attribute("inkscape:stockid")) &&
+                    SP_IS_MARKER(&child))
+                {
+                    object = &child;
+                }
+            }
+        }
+        else if (!strcmp(base,"pattern") && !stock)  {
+            for (auto& child: defs->children)
+            {
+                if (child.getRepr()->attribute("inkscape:stockid") &&
+                    !strcmp(name_p, child.getRepr()->attribute("inkscape:stockid")) &&
+                    SP_IS_PATTERN(&child))
+                {
+                    object = &child;
+                }
+            }
+        }
+        else if (!strcmp(base,"gradient") && !stock)  {
+            for (auto& child: defs->children)
+            {
+                if (child.getRepr()->attribute("inkscape:stockid") &&
+                    !strcmp(name_p, child.getRepr()->attribute("inkscape:stockid")) &&
+                    SP_IS_GRADIENT(&child))
+                {
+                    object = &child;
+                }
+            }
+        }
+
+        if (object == nullptr) {
+
+            if (!strcmp(base, "marker"))  {
+                object = sp_marker_load_from_svg(name_p, doc);
+            }
+            else if (!strcmp(base, "pattern"))  {
+                object = sp_pattern_load_from_svg(name_p, doc);
+            }
+            else if (!strcmp(base, "gradient"))  {
+                object = sp_gradient_load_from_svg(name_p, doc);
+            }
+        }
+
+        g_free(base);
+        g_free(name);
+
+        if (object) {
+            object->setAttribute("inkscape:isstock", "true");
+        }
+
+        return object;
+    }
+
+    else {
+
+        SPDocument *doc = SP_ACTIVE_DOCUMENT;
+        SPObject *object = doc->getObjectById(urn);
+
+        return object;
+    }
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
diff --git a/src/helper/stock-items.h b/src/helper/stock-items.h
new file mode 100644
index 0000000..0d1bb20
--- /dev/null
+++ b/src/helper/stock-items.h
@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/** @file
+ * TODO: insert short description here
+ *//*
+ * Authors:
+ * see git history
+ * John Cliff <simarilius@yahoo.com>
+ *
+ * Copyright (C) 2012 Authors
+ * Released under GNU GPL v2+, read the file 'COPYING' for more information.
+ */
+#ifndef SEEN_INK_STOCK_ITEMS_H
+#define SEEN_INK_STOCK_ITEMS_H
+
+#include <glib.h>
+
+class SPObject;
+
+SPObject *get_stock_item(gchar const *urn, gboolean stock=FALSE);
+
+#endif // SEEN_INK_STOCK_ITEMS_H