Replace cairo curve flattener from Anitgrain Graphics

This also fixes an issue where non-closed curves were not dilated
properly.

8 files changed, 307 insertions, 42 deletions

diff --git a/gerbolyze/gerbolyze.py b/gerbolyze/gerbolyze.py
index 79bc69c..369212c 100755
--- a/gerbolyze/gerbolyze.py
+++ b/gerbolyze/gerbolyze.py
@@ -59,7 +59,8 @@ def vectorize(ctx, side, layer, exact, source, target, image, trace_space):
 @click.option('--bbox', help='Output file bounding box. Format: "w,h" to force [w] mm by [h] mm output canvas OR '
         '"x,y,w,h" to force [w] mm by [h] mm output canvas with its bottom left corner at the given input gerber '
         'coördinates. MUST MATCH --bbox GIVEN TO PREVIEW')
-@click.option('--dilate', default=0.1, help='Default dilation for subtraction operations in mm')
+@click.option('--dilate', default=0.1, type=float, help='Default dilation for subtraction operations in mm')
+@click.option('--curve-tolerance', type=float, help='Tolerance for curve flattening in mm')
 @click.option('--no-subtract', 'no_subtract', flag_value=True, help='Disable subtraction')
 @click.option('--subtract', help='Use user subtraction script from argument (see description above)')
 @click.option('--trace-space', type=float, default=0.1, help='passed through to svg-flatten')
@@ -70,7 +71,7 @@ def vectorize(ctx, side, layer, exact, source, target, image, trace_space):
 def paste(input_gerbers, output_gerbers,
         top, bottom, layer_top, layer_bottom,
         bbox,
-        dilate, no_subtract, subtract,
+        dilate, curve_tolerance, no_subtract, subtract,
         preserve_aspect_ratio,
         trace_space, vectorizer, vectorizer_map, exclude_groups):
     """ Render vector data and raster images from SVG file into gerbers. """
@@ -133,7 +134,7 @@ def paste(input_gerbers, output_gerbers,
             def do_dilate(layer, amount):
                 print('dilating', layer, 'by', amount)
                 outfile = tmpdir / f'dilated-{layer}-{amount}.gbr'
-                dilate_gerber(layers, layer, amount, bbox, tmpdir, outfile, units)
+                dilate_gerber(layers, layer, amount, bbox, tmpdir, outfile, units, curve_tolerance)
                 gbr = gerberex.read(str(outfile))
                 gbr.offset(bounds[0][0], bounds[1][0])
                 return gbr
@@ -157,7 +158,7 @@ def paste(input_gerbers, output_gerbers,
                 overlay_file = tmpdir / f'overlay-{side}-{layer}.gbr'
                 layer_arg = layer if target_layer is None else None # slightly confusing but trust me :)
                 svg_to_gerber(in_svg_or_png, overlay_file, layer_arg,
-                        trace_space, vectorizer, vectorizer_map, exclude_groups,
+                        trace_space, vectorizer, vectorizer_map, exclude_groups, curve_tolerance,
                         bounds_for_png=bounds, preserve_aspect_ratio=preserve_aspect_ratio)
 
                 overlay_grb = gerberex.read(str(overlay_file))
@@ -566,7 +567,7 @@ def create_template_from_svg(bounds, svg_data, extra_layers=DEFAULT_EXTRA_LAYERS
 # SVG/gerber import
 #==================
 
-def dilate_gerber(layers, layer_name, dilation, bbox, tmpdir, outfile, units):
+def dilate_gerber(layers, layer_name, dilation, bbox, tmpdir, outfile, units, curve_tolerance):
     if layer_name not in layers:
         raise ValueError(f'Cannot dilate layer {layer_name}: layer not found in input dir')
 
@@ -592,13 +593,13 @@ def dilate_gerber(layers, layer_name, dilation, bbox, tmpdir, outfile, units):
 
     # dilate & render back to gerber
     # TODO: the scale parameter is a hack. ideally we would fix svg-flatten to handle input units correctly.
-    svg_to_gerber(tmpfile, outfile, dilate=-dilation*72.0/25.4, dpi=72, scale=25.4/72.0)
+    svg_to_gerber(tmpfile, outfile, dilate=-dilation*72.0/25.4, dpi=72, scale=25.4/72.0, curve_tolerance=curve_tolerance)
 
 def svg_to_gerber(infile, outfile,
         layer=None, trace_space:'mm'=0.1,
         vectorizer=None, vectorizer_map=None,
         exclude_groups=None,
-        dilate=None,
+        dilate=None, curve_tolerance=None,
         dpi=None, scale=None, bounds_for_png=None,
         preserve_aspect_ratio=None,
         force_png=False, force_svg=False):
@@ -633,6 +634,9 @@ def svg_to_gerber(infile, outfile,
         args += ['--exclude-groups', exclude_groups]
     if dilate:
         args += ['--dilate', str(dilate)]
+    if curve_tolerance is not None:
+        print('applying curve tolerance', curve_tolerance)
+        args += ['--curve-tolerance', str(curve_tolerance)]
     if dpi:
         args += ['--usvg-dpi', str(dpi)]
     if scale:
diff --git a/svg-flatten/Makefile b/svg-flatten/Makefile
index 5f22a59..ca776d3 100644
--- a/svg-flatten/Makefile
+++ b/svg-flatten/Makefile
@@ -17,6 +17,7 @@ SOURCES := src/svg_color.cpp \
 	src/vec_core.cpp \
 	src/vec_grid.cpp \
 	src/main.cpp \
+	src/flatten.cpp \
 	src/out_svg.cpp \
 	src/out_gerber.cpp \
 	src/out_sexp.cpp \
diff --git a/svg-flatten/include/flatten.hpp b/svg-flatten/include/flatten.hpp
new file mode 100644
index 0000000..b620890
--- /dev/null
+++ b/svg-flatten/include/flatten.hpp
@@ -0,0 +1,27 @@
+
+#include "gerbolyze.hpp"
+
+namespace gerbolyze {
+    class curve4_div {
+        public:
+            curve4_div(double distance_tolerance=0.1, double angle_tolerance=0.0, double cusp_limit=0.0)
+                : m_distance_tolerance_square(0.25*distance_tolerance*distance_tolerance),
+                m_angle_tolerance(angle_tolerance),
+                m_cusp_limit(cusp_limit)
+                {
+                }
+
+            void run(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4);
+            const std::vector<d2p> &points() { return m_points; }
+
+        private:
+            void recursive_bezier(double x1, double y1, double x2, double y2, 
+                                  double x3, double y3, double x4, double y4,
+                                  unsigned level);
+            double m_cusp_limit;
+            double m_distance_tolerance_square;
+            double m_angle_tolerance;
+            std::vector<d2p> m_points;
+    };
+}
+
diff --git a/svg-flatten/src/flatten.cpp b/svg-flatten/src/flatten.cpp
new file mode 100644
index 0000000..e93f044
--- /dev/null
+++ b/svg-flatten/src/flatten.cpp
@@ -0,0 +1,231 @@
+/* Copied from Antigrain Graphics (AGG) v2.4 */
+/* Mirror: https://github.com/pelson/antigrain/blob/master/agg-2.4/src/agg_curves.cpp */
+
+#include <flatten.hpp>
+#include <cmath>
+
+using namespace gerbolyze;
+
+namespace gerbolyze {
+    const double curve_distance_epsilon                  = 1e-15;
+    const double curve_collinearity_epsilon              = 1e-15;
+    const double curve_angle_tolerance_epsilon           = 0.1;
+    constexpr unsigned curve_recursion_limit             = 20;
+}
+
+static inline double calc_sq_distance(double x1, double y1, double x2, double y2)
+{
+    double dx = x2-x1;
+    double dy = y2-y1;
+    return dx * dx + dy * dy;
+}
+
+void curve4_div::run(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4) {
+    m_points.clear();
+    m_points.emplace_back(d2p{x1, y1});
+    recursive_bezier(x1, y1, x2, y2, x3, y3, x4, y4, 0);
+    m_points.emplace_back(d2p{x4, y4});
+}
+
+void curve4_div::recursive_bezier(double x1, double y1,
+                                  double x2, double y2,
+                                  double x3, double y3,
+                                  double x4, double y4,
+                                  unsigned level)
+{
+    if(level > curve_recursion_limit) {
+        return;
+    }
+
+    double pi = M_PI;
+
+    // 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(d2p{x2, y2});
+                return;
+            }
+        } else {
+            if(d3 < m_distance_tolerance_square) {
+                m_points.emplace_back(d2p{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(d2p{x23, y23});
+                return;
+            }
+
+            // Angle Condition
+            //----------------------
+            da1 = fabs(atan2(y4 - y3, x4 - x3) - atan2(y3 - y2, x3 - x2));
+            if(da1 >= pi) da1 = 2*pi - da1;
+
+            if(da1 < m_angle_tolerance) {
+                m_points.emplace_back(d2p{x2, y2});
+                m_points.emplace_back(d2p{x3, y3});
+                return;
+            }
+
+            if(m_cusp_limit != 0.0) {
+                if(da1 > m_cusp_limit)
+                {
+                    m_points.emplace_back(d2p{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(d2p{x23, y23});
+                return;
+            }
+
+            // Angle Condition
+            //----------------------
+            da1 = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1));
+            if(da1 >= pi) da1 = 2*pi - da1;
+
+            if(da1 < m_angle_tolerance) {
+                m_points.emplace_back(d2p{x2, y2});
+                m_points.emplace_back(d2p{x3, y3});
+                return;
+            }
+
+            if(m_cusp_limit != 0.0) {
+                if(da1 > m_cusp_limit) {
+                    m_points.emplace_back(d2p{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(d2p{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 >= pi) da1 = 2*pi - da1;
+            if(da2 >= pi) da2 = 2*pi - da2;
+
+            if(da1 + da2 < m_angle_tolerance) {
+                // Finally we can stop the recursion
+                //----------------------
+                m_points.emplace_back(d2p{x23, y23});
+                return;
+            }
+
+            if(m_cusp_limit != 0.0) {
+                if(da1 > m_cusp_limit) {
+                    m_points.emplace_back(d2p{x2, y2});
+                    return;
+                }
+
+                if(da2 > m_cusp_limit) {
+                    m_points.emplace_back(d2p{x3, y3});
+                    return;
+                }
+            }
+        }
+        break;
+    }
+
+    // Continue subdivision
+    //----------------------
+    recursive_bezier(x1, y1, x12, y12, x123, y123, x1234, y1234, level + 1); 
+    recursive_bezier(x1234, y1234, x234, y234, x34, y34, x4, y4, level + 1); 
+}
+
diff --git a/svg-flatten/src/main.cpp b/svg-flatten/src/main.cpp
index 722b356..951b447 100644
--- a/svg-flatten/src/main.cpp
+++ b/svg-flatten/src/main.cpp
@@ -422,7 +422,7 @@ int main(int argc, char **argv) {
     SVGDocument doc;
     cerr << "Loading temporary file " << frob << endl;
     ifstream load_f(frob);
-    if (!doc.load(load_f)) {
+    if (!doc.load(load_f, "/tmp/debug.svg")) {
         cerr <<  "Error loading input file \"" << in_f_name << "\", exiting." << endl;
         return EXIT_FAILURE;
     }
diff --git a/svg-flatten/src/svg_geom.cpp b/svg-flatten/src/svg_geom.cpp
index 385e848..0ca66fe 100644
--- a/svg-flatten/src/svg_geom.cpp
+++ b/svg-flatten/src/svg_geom.cpp
@@ -20,6 +20,7 @@
 
 #include <cmath>
 #include <string>
+#include <iostream>
 #include <sstream>
 #include <queue>
 #include <assert.h>
@@ -112,7 +113,6 @@ static void dehole_polytree_worker(PolyNode &ptree, Paths &out, queue<PolyTree>
             out.push_back(nod->Contour);
 
         } else {
-
             /* Do not add children's children, those were handled in the recursive call above */
             Clipper c;
             c.AddPath(nod->Contour, ptSubject, /* closed= */ true);
diff --git a/svg-flatten/src/svg_import_util.cpp b/svg-flatten/src/svg_import_util.cpp
index efa7701..2193830 100644
--- a/svg-flatten/src/svg_import_util.cpp
+++ b/svg-flatten/src/svg_import_util.cpp
@@ -94,7 +94,7 @@ void gerbolyze::load_cairo_matrix_from_svg(const string &transform, cairo_matrix
 void gerbolyze::apply_cairo_transform_from_svg(cairo_t *cr, const string &transform) {
     cairo_matrix_t mat;
     load_cairo_matrix_from_svg(transform, mat);
-    cairo_transform(cr, &mat); /* or cairo_transform? */
+    cairo_transform(cr, &mat);
 }
 
 /* Cf. https://tools.ietf.org/html/rfc2397 */
diff --git a/svg-flatten/src/svg_path.cpp b/svg-flatten/src/svg_path.cpp
index ce2775d..f27f650 100644
--- a/svg-flatten/src/svg_path.cpp
+++ b/svg-flatten/src/svg_path.cpp
@@ -19,28 +19,26 @@
 #include <cmath>
 #include <assert.h>
 #include <iostream>
+#include <iomanip>
 #include <sstream>
 #include "cairo_clipper.hpp"
 #include "svg_import_defs.h"
 #include "svg_path.h"
+#include "flatten.hpp"
 
 using namespace std;
 
-static void clipper_add_cairo_path(cairo_t *cr, ClipperLib::Clipper &c, bool closed) {
-    ClipperLib::Paths in_poly;
-    ClipperLib::cairo::cairo_to_clipper(cr, in_poly, CAIRO_PRECISION, ClipperLib::cairo::tNone);
-    c.AddPaths(in_poly, ClipperLib::ptSubject, closed);
-}
-
-static pair<bool, bool> path_to_clipper_via_cairo(cairo_t *cr, ClipperLib::Clipper &c_stroke, ClipperLib::Clipper &c_fill, const pugi::char_t *path_data) {
+static pair<bool, bool> path_to_clipper_via_cairo(cairo_t *cr, ClipperLib::Clipper &c_stroke, ClipperLib::Clipper &c_fill, const pugi::char_t *path_data, double distance_tolerance_mm) {
     istringstream d(path_data);
 
     string cmd;
     double x, y, c1x, c1y, c2x, c2y;
 
+    ClipperLib::Path in_poly;
+    double scale = pow(10.0, CAIRO_PRECISION);
+
     bool first = true;
     bool has_closed = false;
-    bool path_is_empty = true;
     int num_subpaths = 0;
     while (!d.eof()) {
         d >> cmd;
@@ -48,24 +46,21 @@ static pair<bool, bool> path_to_clipper_via_cairo(cairo_t *cr, ClipperLib::Clipp
         assert(!first || cmd == "M");
         
         if (cmd == "Z") { /* Close path */
-            cairo_close_path(cr);
-            clipper_add_cairo_path(cr, c_stroke, /* closed= */ true);
-            clipper_add_cairo_path(cr, c_fill, /* closed= */ true);
+            c_stroke.AddPath(in_poly, ClipperLib::ptSubject, true);
+            c_fill.AddPath(in_poly, ClipperLib::ptSubject, true);
+
             has_closed = true;
-            cairo_new_path(cr);
-            path_is_empty = true;
+            in_poly.clear();
             num_subpaths += 1;
 
         } else if (cmd == "M") { /* Move to */
-            if (!first && !path_is_empty) {
-                cairo_close_path(cr);
-                clipper_add_cairo_path(cr, c_stroke, /* closed= */ false);
-                clipper_add_cairo_path(cr, c_fill, /* closed= */ true);
+            if (!first && !in_poly.empty()) {
+                c_stroke.AddPath(in_poly, ClipperLib::ptSubject, false);
+                c_fill.AddPath(in_poly, ClipperLib::ptSubject, true);
                 num_subpaths += 1;
+                in_poly.clear();
             }
 
-            cairo_new_path (cr);
-
             d >> x >> y;
             /* We need to transform all points ourselves here, and cannot use the transform feature of cairo_to_clipper:
              * Our transform may contain offsets, and clipper only passes its data into cairo's transform functions
@@ -75,17 +70,18 @@ static pair<bool, bool> path_to_clipper_via_cairo(cairo_t *cr, ClipperLib::Clipp
              */
             cairo_user_to_device(cr, &x, &y);
             assert (!d.fail());
-            path_is_empty = true;
-            cairo_move_to(cr, x, y);
+
+            in_poly.emplace_back(ClipperLib::IntPoint{(ClipperLib::cInt)round(x*scale), (ClipperLib::cInt)round(y*scale)});
 
         } else if (cmd == "L") { /* Line to */
             d >> x >> y;
             cairo_user_to_device(cr, &x, &y);
             assert (!d.fail());
-            cairo_line_to(cr, x, y);
-            path_is_empty = false;
+
+            in_poly.emplace_back(ClipperLib::IntPoint{(ClipperLib::cInt)round(x*scale), (ClipperLib::cInt)round(y*scale)});
 
         } else { /* Curve to */
+            double sx = x, sy = y;
             assert(cmd == "C");
             d >> c1x >> c1y; /* first control point */
             cairo_user_to_device(cr, &c1x, &c1y);
@@ -94,16 +90,20 @@ static pair<bool, bool> path_to_clipper_via_cairo(cairo_t *cr, ClipperLib::Clipp
             d >> x >> y; /* end point */
             cairo_user_to_device(cr, &x, &y);
             assert (!d.fail());
-            cairo_curve_to(cr, c1x, c1y, c2x, c2y, x, y);
-            path_is_empty = false;
+
+            gerbolyze::curve4_div c4div(distance_tolerance_mm);
+            c4div.run(sx, sy, c1x, c1y, c2x, c2y, x, y);
+
+            for (auto &pt : c4div.points()) {
+                in_poly.emplace_back(ClipperLib::IntPoint{(ClipperLib::cInt)round(pt[0]*scale), (ClipperLib::cInt)round(pt[1]*scale)});
+            }
         }
 
         first = false;
     }
-    if (!path_is_empty) {
-        cairo_close_path(cr);
-        clipper_add_cairo_path(cr, c_stroke, /* closed= */ false);
-        clipper_add_cairo_path(cr, c_fill, /* closed= */ true);
+    if (!in_poly.empty()) {
+        c_stroke.AddPath(in_poly, ClipperLib::ptSubject, false);
+        c_fill.AddPath(in_poly, ClipperLib::ptSubject, true);
         num_subpaths += 1;
     }
 
@@ -117,14 +117,13 @@ void gerbolyze::load_svg_path(cairo_t *cr, const pugi::xml_node &node, ClipperLi
     /* For open paths, clipper does not correctly remove self-intersections. Thus, we pass everything into
      * clipper twice: Once with all paths set to "closed" to compute fill areas, and once with correct
      * open/closed properties for stroke offsetting. */
-    cairo_set_tolerance (cr, curve_tolerance); /* FIXME make configurable, scale properly for units */
     cairo_set_fill_rule(cr, CAIRO_FILL_RULE_WINDING);
 
     ClipperLib::Clipper c_stroke;
     ClipperLib::Clipper c_fill;
     c_stroke.StrictlySimple(true);
     c_fill.StrictlySimple(true);
-    auto res = path_to_clipper_via_cairo(cr, c_stroke, c_fill, path_data);
+    auto res = path_to_clipper_via_cairo(cr, c_stroke, c_fill, path_data, curve_tolerance);
     bool has_closed = res.first, has_multiple = res.second;
 
     if (!has_closed && !has_multiple) {
@@ -145,9 +144,12 @@ void gerbolyze::load_svg_path(cairo_t *cr, const pugi::xml_node &node, ClipperLi
         auto le_min = -ClipperLib::loRange;
         auto le_max = ClipperLib::hiRange;
         ClipperLib::Path p = {{le_min, le_min}, {le_max, le_min}, {le_max, le_max}, {le_min, le_max}};
+
         c_stroke.AddPath(p, ClipperLib::ptClip, /* closed= */ true);
         c_stroke.Execute(ClipperLib::ctIntersection, ptree_stroke, fill_rule, ClipperLib::pftNonZero);
-        ptree_fill.Clear();
+
+        c_fill.AddPath(p, ClipperLib::ptClip, /* closed= */ true);
+        c_fill.Execute(ClipperLib::ctIntersection, ptree_fill, fill_rule, ClipperLib::pftNonZero);
 
     } else {
         /* We cannot clip the polygon here since that would produce incorrect results for our stroke. */