Improve precision of offset_cubic

(bzr r15291)

2 files changed, 195 insertions, 48 deletions

diff --git a/src/helper/geom-pathstroke.cpp b/src/helper/geom-pathstroke.cpp
index 10641692d..50627fd0b 100644
--- a/src/helper/geom-pathstroke.cpp
+++ b/src/helper/geom-pathstroke.cpp
@@ -1,6 +1,7 @@
 /* Authors:
  *   Liam P. White
  *   Tavmjong Bah
+ *   Alexander Brock
  *
  * Copyright (C) 2014-2015 Authors
  *
@@ -746,33 +747,36 @@ void get_cubic_data(Geom::CubicBezier const& bez, double time, double& len, doub
     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 @@ void offset_cubic(Geom::Path& p, Geom::CubicBezier const& bez, double width, dou
     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);
+    // 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 (err < tol) {
+    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 @@ void offset_cubic(Geom::Path& p, Geom::CubicBezier const& bez, double width, dou
         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 @@ void offset_quadratic(Geom::Path& p, Geom::QuadraticBezier const& bez, double wi
     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 @@ void offset_curve(Geom::Path& res, Geom::Curve const* current, double width)
             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 @@ void peak_cap(Geom::PathBuilder& res, Geom::Path const& with_dir, Geom::Path con
 
 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 @@ Geom::PathVector outline(Geom::Path const& input, double width, double miter, Li
     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 @@ Geom::Path half_outline(Geom::Path const& input, double width, double miter, Lin
     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 @@ Geom::Path half_outline(Geom::Path const& input, double width, double miter, Lin
         // 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 @@ Geom::Path half_outline(Geom::Path const& input, double width, double miter, Lin
         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 @@ void outline_join(Geom::Path &res, Geom::Path const& temp, Geom::Point in_tang,
 {
     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);
diff --git a/src/helper/geom-pathstroke.h b/src/helper/geom-pathstroke.h
index 6a753a6fd..4e4b8b91e 100644
--- a/src/helper/geom-pathstroke.h
+++ b/src/helper/geom-pathstroke.h
@@ -42,12 +42,19 @@ enum LineCapType {
  * @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 @@ Geom::PathVector outline(Geom::Path const& input, double width, double miter, Li
  * @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.