Cubic curve distance optimization

core/Shape.cpp

@@ -3,6 +3,7 @@
 
 #include <cstdlib>
 #include "arithmetics.hpp"
+#include "convergent-curve-ordering.h"
 
 #define DECONVERGE_OVERSHOOT 1.11111111111111111 // moves control points slightly more than necessary to account for floating-point errors
 
@@ -79,8 +80,7 @@ void Shape::normalize() {
                 if (dotProduct(prevDir, curDir) < MSDFGEN_CORNER_DOT_EPSILON-1) {
                     double factor = DECONVERGE_OVERSHOOT*sqrt(1-(MSDFGEN_CORNER_DOT_EPSILON-1)*(MSDFGEN_CORNER_DOT_EPSILON-1))/(MSDFGEN_CORNER_DOT_EPSILON-1);
                     Vector2 axis = factor*(curDir-prevDir).normalize();
-                    // Determine curve ordering using third-order derivative (t = 0) of crossProduct((*prevEdge)->point(1-t)-p0, (*edge)->point(t)-p0) where p0 is the corner (*edge)->point(0)
-                    if (crossProduct((*prevEdge)->directionChange(1), (*edge)->direction(0))+crossProduct((*edge)->directionChange(0), (*prevEdge)->direction(1)) < 0)
+                    if (convergentCurveOrdering(*prevEdge, *edge) < 0)
                         axis = -axis;
                     deconvergeEdge(*prevEdge, 1, axis.getOrthogonal(true));
                     deconvergeEdge(*edge, 0, axis.getOrthogonal(false));

core/convergent-curve-ordering.cpp

@@ -0,0 +1,140 @@
+
+#include "convergent-curve-ordering.h"
+
+#include "arithmetics.hpp"
+#include "Vector2.hpp"
+
+/*
+ * For non-degenerate curves A(t), B(t) (ones where all control points are distinct) both originating at P = A(0) = B(0) = *corner,
+ * we are computing the limit of
+ *
+ *     sign(crossProduct( A(t / |A'(0)|) - P, B(t / |B'(0)|) - P ))
+ *
+ * for t -> 0 from 1. Of note is that the curves' parameter has to be normed by the first derivative at P,
+ * which ensures that the limit approaches P at the same rate along both curves - omitting this was the main error of earlier versions of deconverge.
+ *
+ * For degenerate cubic curves (ones where the first control point equals the origin point), the denominator |A'(0)| is zero,
+ * so to address that, we approach with the square root of t and use the derivative of A(sqrt(t)), which at t = 0 equals A''(0)/2
+ * Therefore, in these cases, we replace one factor of the cross product with A(sqrt(2*t / |A''(0)|)) - P
+ *
+ * The cross product results in a polynomial (in respect to t or t^2 in the degenerate case),
+ * the limit of sign of which at zero can be determined by the lowest order non-zero derivative,
+ * which equals to the sign of the first non-zero polynomial coefficient in the order of increasing exponents.
+ *
+ * The polynomial's constant and linear terms are zero, so the first derivative is definitely zero as well.
+ * The second derivative is assumed to be zero (or near zero) due to the curves being convergent - this is an input requirement
+ * (otherwise the correct result is the sign of the cross product of their directions at t = 0).
+ * Therefore, we skip the first and second derivatives.
+ */
+
+namespace msdfgen {
+
+static void simplifyDegenerateCurve(Point2 *controlPoints, int &order) {
+    if (order == 3 && (controlPoints[1] == controlPoints[0] || controlPoints[1] == controlPoints[3]) && (controlPoints[2] == controlPoints[0] || controlPoints[2] == controlPoints[3])) {
+        controlPoints[1] = controlPoints[3];
+        order = 1;
+    }
+    if (order == 2 && (controlPoints[1] == controlPoints[0] || controlPoints[1] == controlPoints[2])) {
+        controlPoints[1] = controlPoints[2];
+        order = 1;
+    }
+    if (order == 1 && controlPoints[0] == controlPoints[1])
+        order = 0;
+}
+
+int convergentCurveOrdering(const Point2 *corner, int controlPointsBefore, int controlPointsAfter) {
+    if (!(controlPointsBefore > 0 && controlPointsAfter > 0))
+        return 0;
+    Vector2 a1, a2, a3, b1, b2, b3;
+    a1 = *(corner-1)-*corner;
+    b1 = *(corner+1)-*corner;
+    if (controlPointsBefore >= 2)
+        a2 = *(corner-2)-*(corner-1)-a1;
+    if (controlPointsAfter >= 2)
+        b2 = *(corner+2)-*(corner+1)-b1;
+    if (controlPointsBefore >= 3) {
+        a3 = *(corner-3)-*(corner-2)-(*(corner-2)-*(corner-1))-a2;
+        a2 *= 3;
+    }
+    if (controlPointsAfter >= 3) {
+        b3 = *(corner+3)-*(corner+2)-(*(corner+2)-*(corner+1))-b2;
+        b2 *= 3;
+    }
+    a1 *= controlPointsBefore;
+    b1 *= controlPointsAfter;
+    // Non-degenerate case
+    if (a1 && b1) {
+        double as = a1.length();
+        double bs = b1.length();
+        // Third derivative
+        if (double d = as*crossProduct(a1, b2) + bs*crossProduct(a2, b1))
+            return sign(d);
+        // Fourth derivative
+        if (double d = as*as*crossProduct(a1, b3) + as*bs*crossProduct(a2, b2) + bs*bs*crossProduct(a3, b1))
+            return sign(d);
+        // Fifth derivative
+        if (double d = as*crossProduct(a2, b3) + bs*crossProduct(a3, b2))
+            return sign(d);
+        // Sixth derivative
+        return sign(crossProduct(a3, b3));
+    }
+    // Degenerate curve after corner (control point after corner equals corner)
+    int s = 1;
+    if (a1) { // !b1
+        // Swap aN <-> bN and handle in if (b1)
+        b1 = a1;
+        a1 = b2, b2 = a2, a2 = a1;
+        a1 = b3, b3 = a3, a3 = a1;
+        s = -1; // make sure to also flip output
+    }
+    // Degenerate curve before corner (control point before corner equals corner)
+    if (b1) { // !a1
+        // Two-and-a-half-th derivative
+        if (double d = crossProduct(a3, b1))
+            return s*sign(d);
+        // Third derivative
+        if (double d = crossProduct(a2, b2))
+            return s*sign(d);
+        // Three-and-a-half-th derivative
+        if (double d = crossProduct(a3, b2))
+            return s*sign(d);
+        // Fourth derivative
+        if (double d = crossProduct(a2, b3))
+            return s*sign(d);
+        // Four-and-a-half-th derivative
+        return s*sign(crossProduct(a3, b3));
+    }
+    // Degenerate curves on both sides of the corner (control point before and after corner equals corner)
+    { // !a1 && !b1
+        // Two-and-a-half-th derivative
+        if (double d = sqrt(a2.length())*crossProduct(a2, b3) + sqrt(b2.length())*crossProduct(a3, b2))
+            return sign(d);
+        // Third derivative
+        return sign(crossProduct(a3, b3));
+    }
+}
+
+int convergentCurveOrdering(const EdgeSegment *a, const EdgeSegment *b) {
+    Point2 controlPoints[12];
+    Point2 *corner = controlPoints+4;
+    Point2 *aCpTmp = controlPoints+8;
+    int aOrder = int(a->type());
+    int bOrder = int(b->type());
+    if (!(aOrder >= 1 && aOrder <= 3 && bOrder >= 1 && bOrder <= 3)) {
+        // Not implemented - only linear, quadratic, and cubic curves supported
+        return 0;
+    }
+    for (int i = 0; i <= aOrder; ++i)
+        aCpTmp[i] = a->controlPoints()[i];
+    for (int i = 0; i <= bOrder; ++i)
+        corner[i] = b->controlPoints()[i];
+    if (aCpTmp[aOrder] != *corner)
+        return 0;
+    simplifyDegenerateCurve(aCpTmp, aOrder);
+    simplifyDegenerateCurve(corner, bOrder);
+    for (int i = 0; i < aOrder; ++i)
+        corner[i-aOrder] = aCpTmp[i];
+    return convergentCurveOrdering(corner, aOrder, bOrder);
+}
+
+}

core/convergent-curve-ordering.h

@@ -0,0 +1,11 @@
+
+#pragma once
+
+#include "edge-segments.h"
+
+namespace msdfgen {
+
+/// For curves a, b converging at P = a->point(1) = b->point(0) with the same (opposite) direction, determines the relative ordering in which they exit P (i.e. whether a is to the left or right of b at the smallest positive radius around P)
+int convergentCurveOrdering(const EdgeSegment *a, const EdgeSegment *b);
+
+}

core/edge-segments.cpp

@@ -199,9 +199,9 @@ SignedDistance QuadraticSegment::signedDistance(Point2 origin, double &param) co
     double minDistance = nonZeroSign(crossProduct(epDir, qa))*qa.length(); // distance from A
     param = -dotProduct(qa, epDir)/dotProduct(epDir, epDir);
     {
-        epDir = direction(1);
         double distance = (p[2]-origin).length(); // distance from B
         if (distance < fabs(minDistance)) {
+            epDir = direction(1);
             minDistance = nonZeroSign(crossProduct(epDir, p[2]-origin))*distance;
             param = dotProduct(origin-p[1], epDir)/dotProduct(epDir, epDir);
         }
@@ -235,25 +235,31 @@ SignedDistance CubicSegment::signedDistance(Point2 origin, double &param) const
     double minDistance = nonZeroSign(crossProduct(epDir, qa))*qa.length(); // distance from A
     param = -dotProduct(qa, epDir)/dotProduct(epDir, epDir);
     {
-        epDir = direction(1);
         double distance = (p[3]-origin).length(); // distance from B
         if (distance < fabs(minDistance)) {
+            epDir = direction(1);
             minDistance = nonZeroSign(crossProduct(epDir, p[3]-origin))*distance;
             param = dotProduct(epDir-(p[3]-origin), epDir)/dotProduct(epDir, epDir);
         }
     }
     // Iterative minimum distance search
     for (int i = 0; i <= MSDFGEN_CUBIC_SEARCH_STARTS; ++i) {
-        double t = (double) i/MSDFGEN_CUBIC_SEARCH_STARTS;
+        double t = 1./MSDFGEN_CUBIC_SEARCH_STARTS*i;
         Vector2 qe = qa+3*t*ab+3*t*t*br+t*t*t*as;
-        for (int step = 0; step < MSDFGEN_CUBIC_SEARCH_STEPS; ++step) {
-            // Improve t
-            Vector2 d1 = 3*ab+6*t*br+3*t*t*as;
-            Vector2 d2 = 6*br+6*t*as;
-            t -= dotProduct(qe, d1)/(dotProduct(d1, d1)+dotProduct(qe, d2));
-            if (t <= 0 || t >= 1)
-                break;
-            qe = qa+3*t*ab+3*t*t*br+t*t*t*as;
+        Vector2 d1 = 3*ab+6*t*br+3*t*t*as;
+        Vector2 d2 = 6*br+6*t*as;
+        double improvedT = t-dotProduct(qe, d1)/(dotProduct(d1, d1)+dotProduct(qe, d2));
+        if (improvedT > 0 && improvedT < 1) {
+            int remainingSteps = MSDFGEN_CUBIC_SEARCH_STEPS;
+            do {
+                t = improvedT;
+                qe = qa+3*t*ab+3*t*t*br+t*t*t*as;
+                d1 = 3*ab+6*t*br+3*t*t*as;
+                if (!--remainingSteps)
+                    break;
+                d2 = 6*br+6*t*as;
+                improvedT = t-dotProduct(qe, d1)/(dotProduct(d1, d1)+dotProduct(qe, d2));
+            } while (improvedT > 0 && improvedT < 1);
             double distance = qe.length();
             if (distance < fabs(minDistance)) {
                 minDistance = nonZeroSign(crossProduct(d1, qe))*distance;