Improve precision of offset_cubic

(bzr r15280.1.1)

1 files changed, 149 insertions, 31 deletions

diff --git a/src/helper/geom-pathstroke.cpp b/src/helper/geom-pathstroke.cpp
index 10641692d..905984a2a 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 = 1; ii <= 9; ++ii) {
+        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);
     }
@@ -829,7 +947,7 @@ void offset_quadratic(Geom::Path& p, Geom::QuadraticBezier const& bez, double wi
 
 void offset_curve(Geom::Path& res, Geom::Curve const* current, double width)
 {
-    double const tolerance = 0.0025;
+    double const tolerance = 1.0 * (width / 100); // Tolerance is 1% of the given width.
     size_t levels = 8;
 
     if (current->isDegenerate()) return; // don't do anything