powerstroke: properly implement rounded joins. brain buster!!

(bzr r11111)

1 files changed, 81 insertions, 12 deletions

diff --git a/src/live_effects/lpe-powerstroke.cpp b/src/live_effects/lpe-powerstroke.cpp
index 257dba8ea..6cc2d2368 100644
--- a/src/live_effects/lpe-powerstroke.cpp
+++ b/src/live_effects/lpe-powerstroke.cpp
@@ -26,9 +26,10 @@
 #include <2geom/svg-path.h>
 #include <2geom/path-intersection.h>
 #include <2geom/crossing.h>
+#include <2geom/ellipse.h>
 
 namespace Geom {
-
+// should all be moved to 2geom at some point
 
 Point unitTangentAt( D2<SBasis> const & a, Coord t, unsigned n = 3) {
     std::vector<Point> derivs = a.valueAndDerivatives(t, n);
@@ -55,8 +56,48 @@ boost::optional<Point> intersection_point( Point const & origin_a, Point const &
     return boost::none;
 }
 
+Geom::CubicBezier sbasis_to_cubicbezier(Geom::D2<Geom::SBasis> const & sbasis_in)
+{
+    std::vector<Geom::Point> temp;
+    sbasis_to_bezier(temp, sbasis_in, 4);
+    return Geom::CubicBezier( temp );
 }
 
+/**
+ * document this!
+ * very quick: this finds the ellipse with minimum eccentricity
+   passing through point P and Q, with tangent PO at P and QO at Q
+   http://mathforum.org/kb/message.jspa?messageID=7471596&tstart=0
+ */
+static Ellipse find_ellipse(Point P, Point Q, Point O)
+{
+    Point p = P - O;
+    Point q = Q - O;
+    Coord K = 4 * dot(p,q) / (L2sq(p) + L2sq(q));
+
+    double cross = p[Y]*q[X] - p[X]*q[Y];
+    double a = -q[Y]/cross;
+    double b = q[X]/cross;
+    double c = (O[X]*q[Y] - O[Y]*q[X])/cross;
+
+    double d = p[Y]/cross;
+    double e = -p[X]/cross;
+    double f = (-O[X]*p[Y] + O[Y]*p[X])/cross;
+
+    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+    double A = (a*d*K+d*d+a*a);
+    double B = (a*e*K+b*d*K+2*d*e+2*a*b);
+    double C = (b*e*K+e*e+b*b);
+    double D = (a*f*K+c*d*K+2*d*f-2*d+2*a*c-2*a);
+    double E = (b*f*K+c*e*K+2*e*f-2*e+2*b*c-2*b);
+    double F = c*f*K+f*f-2*f+c*c-2*c+1;
+
+    return Ellipse(A, B, C, D, E, F);
+}
+
+
+} // namespace Geom
+
 namespace Inkscape {
 namespace LivePathEffect {
 
@@ -176,8 +217,16 @@ std::vector<discontinuity_data> find_discontinuities( Geom::Piecewise<Geom::D2<G
     for(unsigned i = 1; i < der.size(); i++) {
         if ( ! are_near(der[i-1].at1(), der[i].at0(), eps) ) {
             discontinuity_data data;
+
             data.der0 = der[i-1].at1();
             data.der1 = der[i].at0();
+            if ( Geom::are_near(data.der0.length(), 0) ) {
+                data.der0 = unitTangentAt(der[i-1], 1, 2);
+            }
+            if ( Geom::are_near(data.der1.length(), 0) ) {
+                data.der1 = unitTangentAt(der[i], 0, 2);
+            }
+
             double t = der.cuts[i];
             std::vector< double > rts = roots (x - t);  /// @todo this has multiple solutions for general strokewidth paths (generated by spiro interpolator...), ignore for now
             if (!rts.empty()) {
@@ -224,11 +273,36 @@ Geom::Path path_from_piecewise_fix_cusps( Geom::Piecewise<Geom::D2<Geom::SBasis>
             discontinuity_data cusp = cusps[cusp_i];
 
             switch (cusp_linecap) {
-            case LINECUSP_ROUND:  // properly bugged ^_^
-                pb.arcTo( abs(cusp.width), abs(cusp.width),
-                          angle_between(cusp.der0, cusp.der1), false, cusp.width < 0,
-                          B[i].at0() );
+            case LINECUSP_ROUND: {
+                if ( sign*cusp.width*angle_between(cusp.der0, cusp.der1) < 0.) {
+                    // we are on the outside: round corner
+                    /* for constant width paths, the rounding is a circular arc (rx == ry),
+                       for non-constant width paths, the rounding can be done with an ellipse but is hard and ambiguous.
+                       The elliptical arc should go through the discontinuity's start and end points (of course!)
+                       and also should match the discontinuity tangents at those start and end points.
+                       To resolve the ambiguity, the elliptical arc with minimal eccentricity should be chosen.
+                       A 2Geom method was created to do exactly this :)
+                       */
+
+                    Geom::Point tang1 = unitTangentAt(B[prev_i],1);
+                    Geom::Point tang2 = unitTangentAt(B[i],0);
+                    boost::optional<Geom::Point> O = intersection_point( B[prev_i].at1(), tang1,
+                                                                              B[i].at0(), tang2 );
+                    if (!O) {
+                        // no center found, i.e. 180 degrees round
+                       pb.lineTo(B[i].at0()); // default to bevel for too shallow cusp angles
+                       break;
+                    }
+
+                    Geom::Ellipse ellipse = find_ellipse(B[prev_i].at1(), B[i].at0(), *O);
+                    pb.arcTo( ellipse.ray(Geom::X), ellipse.ray(Geom::Y), ellipse.rot_angle(),
+                              false, cusp.width < 0, B[i].at0() );
+                } else {
+                    // we are on the inside, do a simple bevel to connect the paths
+                    pb.lineTo(B[i].at0()); // default to bevel for too shallow cusp angles
+                }
                 break;
+            }
 /*            case LINECUSP_NONE: {
                 if ( sign*cusp.width*angle_between(cusp.der0, cusp.der1) < 0.) {
                     // we are on the outside
@@ -250,15 +324,10 @@ Geom::Path path_from_piecewise_fix_cusps( Geom::Piecewise<Geom::D2<Geom::SBasis>
                     Geom::Point der2 = unitTangentAt(B[i],0);
 
                     Geom::D2<Geom::SBasis> newcurve1 = B[prev_i] * Geom::reflection(rot90(der1), B[prev_i].at1());
-                    newcurve1 = reverse(newcurve1);
-                    std::vector<Geom::Point> temp;
-                    sbasis_to_bezier(temp, newcurve1, 4);
-                    Geom::CubicBezier bzr1( temp );
+                    Geom::CubicBezier bzr1 = sbasis_to_cubicbezier( reverse(newcurve1) );
 
                     Geom::D2<Geom::SBasis> newcurve2 = B[i] * Geom::reflection(rot90(der2), B[i].at0());
-                    newcurve2 = reverse(newcurve2);
-                    sbasis_to_bezier(temp, newcurve2, 4);
-                    Geom::CubicBezier bzr2( temp );
+                    Geom::CubicBezier bzr2 = sbasis_to_cubicbezier( reverse(newcurve2) );
 
                     Geom::Crossings cross = crossings(bzr1, bzr2);
                     if (cross.empty()) {