Bezier solver

core/bezier-solver.hpp

@@ -0,0 +1,94 @@
+
+#pragma once
+
+#include <cmath>
+#include "Vector2.hpp"
+
+// Parameters for iterative search of closest point on a cubic Bezier curve. Increase for higher precision.
+#define MSDFGEN_CUBIC_SEARCH_STARTS 4
+#define MSDFGEN_CUBIC_SEARCH_STEPS 4
+
+#define MSDFGEN_QUADRATIC_RATIO_LIMIT 1e8
+
+#ifndef MSDFGEN_CUBE_ROOT
+#define MSDFGEN_CUBE_ROOT(x) pow((x), 1/3.)
+#endif
+
+namespace msdfgen {
+
+/**
+ * Returns the parameter for the quadratic Bezier curve (P0, P1, P2) for the point closest to point P. May be outside the (0, 1) range.
+ * p = P-P0
+ * q = 2*P1-2*P0
+ * r = P2-2*P1+P0
+ */
+inline double quadraticNearPoint(const Vector2 p, const Vector2 q, const Vector2 r) {
+    double qq = q.squaredLength();
+    double rr = r.squaredLength();
+    if (qq >= MSDFGEN_QUADRATIC_RATIO_LIMIT*rr)
+        return dotProduct(p, q)/qq;
+    double norm = .5/rr;
+    double a = 3*norm*dotProduct(q, r);
+    double b = norm*(qq-2*dotProduct(p, r));
+    double c = norm*dotProduct(p, q);
+    double aa = a*a;
+    double g = 1/9.*(aa-3*b);
+    double h = 1/54.*(a*(aa+aa-9*b)-27*c);
+    double hh = h*h;
+    double ggg = g*g*g;
+    a *= 1/3.;
+    if (hh < ggg) {
+        double u = 1/3.*acos(h/sqrt(ggg));
+        g = -2*sqrt(g);
+        if (h >= 0) {
+            double t = g*cos(u)-a;
+            if (t >= 0)
+                return t;
+            return g*cos(u+2.0943951023931954923)-a; // 2.094 = PI*2/3
+        } else {
+            double t = g*cos(u+2.0943951023931954923)-a;
+            if (t <= 1)
+                return t;
+            return g*cos(u)-a;
+        }
+    }
+    double s = (h < 0 ? 1. : -1.)*MSDFGEN_CUBE_ROOT(fabs(h)+sqrt(hh-ggg));
+    return s+g/s-a;
+}
+
+/**
+ * Returns the parameter for the cubic Bezier curve (P0, P1, P2, P3) for the point closest to point P. Squared distance is provided as optional output parameter.
+ * p = P-P0
+ * q = 3*P1-3*P0
+ * r = 3*P2-6*P1+3*P0
+ * s = P3-3*P2+3*P1-P0
+ */
+inline double cubicNearPoint(const Vector2 p, const Vector2 q, const Vector2 r, const Vector2 s, double &squaredDistance) {
+    squaredDistance = p.squaredLength();
+    double bestT = 0;
+    for (int i = 0; i <= MSDFGEN_CUBIC_SEARCH_STARTS; ++i) {
+        double t = 1./MSDFGEN_CUBIC_SEARCH_STARTS*i;
+        Vector2 curP = p-(q+(r+s*t)*t)*t;
+        for (int step = 0; step < MSDFGEN_CUBIC_SEARCH_STEPS; ++step) {
+            Vector2 d0 = q+(r+r+3*s*t)*t;
+            Vector2 d1 = r+r+6*s*t;
+            t += dotProduct(curP, d0)/(d0.squaredLength()-dotProduct(curP, d1));
+            if (t <= 0 || t >= 1)
+                break;
+            curP = p-(q+(r+s*t)*t)*t;
+            double curSquaredDistance = curP.squaredLength();
+            if (curSquaredDistance < squaredDistance) {
+                squaredDistance = curSquaredDistance;
+                bestT = t;
+            }
+        }
+    }
+    return bestT;
+}
+
+inline double cubicNearPoint(const Vector2 p, const Vector2 q, const Vector2 r, const Vector2 s) {
+    double squaredDistance;
+    return cubicNearPoint(p, q, r, s, squaredDistance);
+}
+
+}

core/edge-segments.cpp

@@ -3,6 +3,9 @@
 
 #include "arithmetics.hpp"
 #include "equation-solver.h"
+#include "bezier-solver.hpp"
+
+#define MSDFGEN_USE_BEZIER_SOLVER
 
 namespace msdfgen {
 
@@ -171,6 +174,68 @@ SignedDistance LinearSegment::signedDistance(Point2 origin, double &param) const
     return SignedDistance(nonZeroSign(crossProduct(aq, ab))*endpointDistance, fabs(dotProduct(ab.normalize(), eq.normalize())));
 }
 
+#ifdef MSDFGEN_USE_BEZIER_SOLVER
+
+SignedDistance QuadraticSegment::signedDistance(Point2 origin, double &param) const {
+    Vector2 ap = origin-p[0];
+    Vector2 bp = origin-p[2];
+    Vector2 q = 2*(p[1]-p[0]);
+    Vector2 r = p[2]-2*p[1]+p[0];
+    double aSqD = ap.squaredLength();
+    double bSqD = bp.squaredLength();
+    double t = quadraticNearPoint(ap, q, r);
+    if (t > 0 && t < 1) {
+        Vector2 tp = ap-(q+r*t)*t;
+        double tSqD = tp.squaredLength();
+        if (tSqD < aSqD && tSqD < bSqD) {
+            param = t;
+            return SignedDistance(nonZeroSign(crossProduct(tp, q+2*r*t))*sqrt(tSqD), 0);
+        }
+    }
+    if (bSqD < aSqD) {
+        Vector2 d = q+r+r;
+        if (!d)
+            d = p[2]-p[0];
+        param = dotProduct(bp, d)/d.squaredLength()+1;
+        return SignedDistance(nonZeroSign(crossProduct(bp, d))*sqrt(bSqD), dotProduct(bp.normalize(), d.normalize()));
+    }
+    if (!q)
+        q = p[2]-p[0];
+    param = dotProduct(ap, q)/q.squaredLength();
+    return SignedDistance(nonZeroSign(crossProduct(ap, q))*sqrt(aSqD), -dotProduct(ap.normalize(), q.normalize()));
+}
+
+SignedDistance CubicSegment::signedDistance(Point2 origin, double &param) const {
+    Vector2 ap = origin-p[0];
+    Vector2 bp = origin-p[3];
+    Vector2 q = 3*(p[1]-p[0]);
+    Vector2 r = 3*(p[2]-p[1])-q;
+    Vector2 s = p[3]-3*(p[2]-p[1])-p[0];
+    double aSqD = ap.squaredLength();
+    double bSqD = bp.squaredLength();
+    double tSqD;
+    double t = cubicNearPoint(ap, q, r, s, tSqD);
+    if (t > 0 && t < 1) {
+        if (tSqD < aSqD && tSqD < bSqD) {
+            param = t;
+            return SignedDistance(nonZeroSign(crossProduct(ap-(q+(r+s*t)*t)*t, q+(r+r+3*s*t)*t))*sqrt(tSqD), 0);
+        }
+    }
+    if (bSqD < aSqD) {
+        Vector2 d = q+r+r+3*s;
+        if (!d)
+            d = p[3]-p[1];
+        param = dotProduct(bp, d)/d.squaredLength()+1;
+        return SignedDistance(nonZeroSign(crossProduct(bp, d))*sqrt(bSqD), dotProduct(bp.normalize(), d.normalize()));
+    }
+    if (!q)
+        q = p[2]-p[0];
+    param = dotProduct(ap, q)/q.squaredLength();
+    return SignedDistance(nonZeroSign(crossProduct(ap, q))*sqrt(aSqD), -dotProduct(ap.normalize(), q.normalize()));
+}
+
+#else
+
 SignedDistance QuadraticSegment::signedDistance(Point2 origin, double &param) const {
     Vector2 qa = p[0]-origin;
     Vector2 ab = p[1]-p[0];
@@ -257,6 +322,8 @@ SignedDistance CubicSegment::signedDistance(Point2 origin, double &param) const
         return SignedDistance(minDistance, fabs(dotProduct(direction(1).normalize(), (p[3]-origin).normalize())));
 }
 
+#endif
+
 int LinearSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
     if ((y >= p[0].y && y < p[1].y) || (y >= p[1].y && y < p[0].y)) {
         double param = (y-p[0].y)/(p[1].y-p[0].y);

core/edge-segments.h

@@ -7,10 +7,6 @@
 
 namespace msdfgen {
 
-// Parameters for iterative search of closest point on a cubic Bezier curve. Increase for higher precision.
-#define MSDFGEN_CUBIC_SEARCH_STARTS 4
-#define MSDFGEN_CUBIC_SEARCH_STEPS 4
-
 /// An abstract edge segment.
 class EdgeSegment {
 

core/equation-solver.cpp

@@ -4,6 +4,10 @@
 #define _USE_MATH_DEFINES
 #include <cmath>
 
+#ifndef MSDFGEN_CUBE_ROOT
+#define MSDFGEN_CUBE_ROOT(x) pow((x), 1/3.)
+#endif
+
 namespace msdfgen {
 
 int solveQuadratic(double x[2], double a, double b, double c) {
@@ -49,7 +53,7 @@ static int solveCubicNormed(double x[3], double a, double b, double c) {
         x[2] = q*cos(1/3.*(t-2*M_PI))-a;
         return 3;
     } else {
-        double u = (r < 0 ? 1 : -1)*pow(fabs(r)+sqrt(r2-q3), 1/3.);
+        double u = (r < 0 ? 1 : -1)*MSDFGEN_CUBE_ROOT(fabs(r)+sqrt(r2-q3));
         double v = u == 0 ? 0 : q/u;
         x[0] = (u+v)-a;
         if (u == v || fabs(u-v) < 1e-12*fabs(u+v)) {