Merge lp:~inkscape+alexander/inkscape/offset-cubic-improvement into lp:~inkscape.dev/inkscape/trunk

[Alexander Brock]

Improve precision of offset_cubic

I changed the behaviour of offset_cubic / half_outline:

1. Tolerance (convergence criterion) now depends on width to prevent problems with very thin lines where the old code would calculate a rough approximation.
See attachment for an example.
2. The old code does some constant adjustment to the length of the bezier handles and subdivides if the tolerance is exceeded at (t=0.5).
The new code does iterative adjustment of the length, evaluates the error at t = {0.1, 0.2, ..., 0.9} and subdivides at the point where the error is highest (if neccessary).

[Alexander Brock (brock-alexander)]

Comparison of old and new behaviour:

http://imgur.com/a/9ysxL

Path effect is "Join Type" with parameter "Line width" = 0.1

[Jabiertxof (jabiertxof)]

Hi Alexander. Could you try with this file?
https://inkscape.org/en/~jabiertxof/%E2%98%85drawing-to-fix-line-join

Seems a bit slow and also has a point to (0,0) on the climb.

[Jabiertxof (jabiertxof)]

As far I find the problem is the 0 width segment in the closing point of the mountain

m 101.00553,6.8925054 5.26676,5.3490096 c 0.49943,0.510326 0.49249,1.499268 0.21305,1.783664 l -2.61505,-2.091842 -0.51429,3.09763 -2.18499,-1.153436 -3.498951,2.210752 -1.15839,-4.660317 -1.87979,3.252214 h -2.87368 c -1.171269,0 -1.309008,-1.486387 -0.24475,-2.55064 1.85897,-2.006622 3.99243,-4.0518896 5.15181,-5.2370346 1.16533,-1.191091 3.19573,-1.1574 4.338271,0 z

[Jabiertxof (jabiertxof)]

Maybe you can add code to ignore it?

[Jabiertxof (jabiertxof)]

Ignore all degenerate curves

[Jabiertxof (jabiertxof)]

Geom::Curve::isDegenerate()

[lp:~inkscape+alexander/inkscape/offset-cubic-improvement]

15282. By Alexander Brock <email address hidden> on 2016-11-27

Increase tolerance to match old methods precision and improve speed

15283. By Alexander Brock <email address hidden> on 2016-11-27

Make tolerance of offset_curve a parameter

[Alexander Brock (brock-alexander)]

I added experimental code to ignore degenerate curves, but that changed nothing.

I found the place in the code where the "L 0 0 z" is created:

http://bazaar.launchpad.net/~inkscape+alexander/inkscape/offset-cubic-improvement/view/head:/src/helper/geom-pathstroke.cpp#L1116

content of temp:
M 15.835720153852229 8.14728407690592 C 15.835735083689977 8.14729920020721 15.835750013370388 8.147314323706892 15.835764942893457 8.147329447404967

content of temp2:
M 15.835674601501454 8.147237813173517 L 15.835720153852229 8.14728407690592

Tangents: (-0.7025234740496769, -0.7116605710654292), (0.7016071317351655, 0.7125639849854567)

content of temp after the call to outline_join:

M 15.835720153852229 8.14728407690592 C 15.835735083689977 8.14729920020721 15.835750013370388 8.147314323706892 15.835764942893457 8.147329447404967 L 0 0

I tried to figure out what exactly is happening in outline_join and now my head hurts. Since the problem has nothing to do with the offset_cubic improvements, can we merge this?

[Jabiertxof (jabiertxof)]

Hi Alexander, I just fix the problem as I try to tell in IRC. Also comment some cool improvements copied from IRC:

The Fix:
https://inkscape.org/en/paste/10473/
you need to apply to current trun geom-patstroke.cpp

The comments:
Just to know: is posibol remove overlaps when the offset in a path time value <= 0? Or is too complex. I could do with libarot in mayority cases but for example in the same SVG of climb when line join has a hight value, inside line cross outside line-
If this kind of work is done... simplify at ~80% the offset LPE, realy geom-pathstroke can do all work without boolops and without Livarot.

One cuestion about new tolerance.
¿Can it have a default value with similar process time than previous version? a deefault value in half_outline and in outline so you can use this parent and grandparent functions with tolerance
or with default value.

[lp:~inkscape+alexander/inkscape/offset-cubic-improvement]

15284. By Alexander Brock <email address hidden> on 2016-11-29

Fix instability caused by degenerate paths (patch by jabiertxof)

15285. By Alexander Brock <email address hidden> on 2016-11-29

Make tolerance a parameter on every level

15286. By Alexander Brock <email address hidden> on 2016-11-29

Handle negative width correctly in tolerance calculation

[Alexander Brock (brock-alexander)]

I made the tolerance a parameter of outline and half_outline. If not set (or set to some value <= 0) it is automatically calculated as 5% of the width parameter.

Default tolerance is 5% on my machine and runtime is 6 times slower than the old code in my benchmarks but it's not noticable on my machine in a debug build in the gui.

It creates a smaller number of segments than the old code and the accuracy is much better for small values of width, therefore I propose it for merging.

[Jabiertxof (jabiertxof)]

Thanks for working on this.

review: Approve (code + ux)

=== modified file 'src/helper/geom-pathstroke.cpp'
--- src/helper/geom-pathstroke.cpp	2016-08-03 13:29:38 +0000
+++ src/helper/geom-pathstroke.cpp	2016-11-29 16:21:40 +0000
@@ -1,6 +1,7 @@
 /* Authors:
  *   Liam P. White
  *   Tavmjong Bah
+ *   Alexander Brock
  *
  * Copyright (C) 2014-2015 Authors
  *
@@ -746,33 +747,36 @@
     len = l;
 }
 
-void offset_cubic(Geom::Path& p, Geom::CubicBezier const& bez, double width, double tol, size_t levels)
-{
+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;
 
-    Geom::Point start_pos = bez.initialPoint();
-    Geom::Point end_pos = bez.finalPoint();
-
-    Geom::Point start_normal = Geom::rot90(bez.unitTangentAt(0));
-    Geom::Point end_normal = -Geom::rot90(Geom::unitTangentAt(Geom::reverse(bez.toSBasis()), 0.));
-
-    // offset the start and end control points out by the width
-    Geom::Point start_new = start_pos + start_normal*width;
-    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);
-
     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 / start_rad;
+        start_off += (width + width_correction) / start_rad;
     if (!Geom::are_near(end_rad, 0))
-        end_off += width / end_rad;
+        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;
     // --------
@@ -783,23 +787,137 @@
     mid2_new = Geom::Point(end_new[X] - mid2_new[X]/3., end_new[Y] - mid2_new[Y]/3.);
 
     // create the estimate curve
-    Geom::CubicBezier c = Geom::CubicBezier(start_new, mid1_new, mid2_new, end_new);
+    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) {
-        p.append(c);
+        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;
     }
 
-    // check the tolerance for our estimate to be a parallel curve
-    Geom::Point chk = c.pointAt(.5);
-    Geom::Point req = bez.pointAt(.5) + Geom::rot90(bez.unitTangentAt(.5))*width; // required accuracy
-
-    Geom::Point const diff = req - chk;
-    double const err = Geom::dot(diff, diff);
-
-    if (err < tol) {
+    // 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
         }
@@ -809,7 +927,7 @@
         return;
     } else {
         // split the curve in two
-        std::pair<Geom::CubicBezier, Geom::CubicBezier> s = bez.subdivide(.5);
+        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);
     }
@@ -827,9 +945,8 @@
     offset_cubic(p, cub, width, tol, levels);
 }
 
-void offset_curve(Geom::Path& res, Geom::Curve const* current, double width)
+void offset_curve(Geom::Path& res, Geom::Curve const* current, double width, double tolerance)
 {
-    double const tolerance = 0.0025;
     size_t levels = 8;
 
     if (current->isDegenerate()) return; // don't do anything
@@ -855,14 +972,14 @@
             default: {
                 Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(current->toSBasis(), tolerance);
                 for (size_t i = 0; i < sbasis_path.size(); ++i)
-                    offset_curve(res, &sbasis_path[i], width);
+                    offset_curve(res, &sbasis_path[i], width, tolerance);
                 break;
             }
         }
     } else {
         Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(current->toSBasis(), 0.1);
         for (size_t i = 0; i < sbasis_path.size(); ++i)
-            offset_curve(res, &sbasis_path[i], width);
+            offset_curve(res, &sbasis_path[i], width, tolerance);
     }
 }
 
@@ -902,13 +1019,19 @@
 
 namespace Inkscape {
 
-Geom::PathVector outline(Geom::Path const& input, double width, double miter, LineJoinType join, LineCapType butt)
+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);
-    Geom::Path against_dir = half_outline(input.reversed(), width/2., miter, join);
+    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);
@@ -947,8 +1070,21 @@
     return res.peek();
 }
 
-Geom::Path half_outline(Geom::Path const& input, double width, double miter, LineJoinType join)
+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;
 
@@ -963,12 +1099,13 @@
     res.start(start);
 
     // Do two curves at a time for efficiency, since the join function needs to know the outgoing curve as well
-    const size_t k = (input.back_closed().isDegenerate() && input.closed())
-            ?input.size_default()-1:input.size_default();
+    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);
+        offset_curve(temp, &input[u], width, tolerance);
 
         // on the first run through, there isn't a join
         if (u == 0) {
@@ -981,12 +1118,11 @@
         // odd number of paths
         if (u < k - 1) {
             temp.clear();
-            offset_curve(temp, &input[u+1], width);
+            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();
@@ -998,7 +1134,6 @@
         outline_join(temp, temp2, tang[0], tang[1], width, miter, join);
         res.erase(res.begin());
         res.erase_last();
-        //
         res.append(temp);
         res.close();
     }
@@ -1010,9 +1145,8 @@
 {
     if (res.size() == 0 || temp.size() == 0)
         return;
-
     Geom::Curve const& outgoing = temp.front();
-    if (Geom::are_near(res.finalPoint(), outgoing.initialPoint())) {
+    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);

=== modified file 'src/helper/geom-pathstroke.h'
--- src/helper/geom-pathstroke.h	2016-02-15 09:52:41 +0000
+++ src/helper/geom-pathstroke.h	2016-11-29 16:21:40 +0000
@@ -42,12 +42,19 @@
  * @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);
+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.
@@ -57,10 +64,16 @@
  * @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);
+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.