Updated API of approximateSDF

README.md

@@ -113,25 +113,23 @@ in order to generate a distance field. Please note that all classes and function
 
 Example:
 ```c++
-#include "msdfgen.h"
-#include "msdfgen-ext.h"
+#include <msdfgen.h>
+#include <msdfgen-ext.h>
 
 using namespace msdfgen;
 
 int main() {
-    FreetypeHandle *ft = initializeFreetype();
-    if (ft) {
-        FontHandle *font = loadFont(ft, "C:\\Windows\\Fonts\\arialbd.ttf");
-        if (font) {
+    if (FreetypeHandle *ft = initializeFreetype()) {
+        if (FontHandle *font = loadFont(ft, "C:\\Windows\\Fonts\\arialbd.ttf")) {
             Shape shape;
-            if (loadGlyph(shape, font, 'A')) {
+            if (loadGlyph(shape, font, 'A', FONT_SCALING_EM_NORMALIZED)) {
                 shape.normalize();
                 //                      max. angle
                 edgeColoringSimple(shape, 3.0);
                 //          output width, height
                 Bitmap<float, 3> msdf(32, 32);
-                //                           scale, translation
-                SDFTransformation t(Projection(1.0, Vector2(4.0, 4.0)), Range(4.0));
+                //                            scale, translation (in em's)
+                SDFTransformation t(Projection(32.0, Vector2(0.125, 0.125)), Range(0.125));
                 generateMSDF(msdf, shape, t);
                 savePng(msdf, "output.png");
             }
@@ -141,7 +139,6 @@ int main() {
     }
     return 0;
 }
-
 ```
 
 ## Using a multi-channel distance field

core/MSDFErrorCorrection.cpp

@@ -129,10 +129,10 @@ void MSDFErrorCorrection::protectCorners(const Shape &shape) {
                 if (!(commonColor&(commonColor-1))) {
                     // Find the four texels that envelop the corner and mark them as protected.
                     Point2 p = transformation.project((*edge)->point(0));
-                    if (shape.inverseYAxis)
-                        p.y = stencil.height-p.y;
                     int l = (int) floor(p.x-.5);
                     int b = (int) floor(p.y-.5);
+                    if (shape.inverseYAxis)
+                        b = stencil.height-b-2;
                     int r = l+1;
                     int t = b+1;
                     // Check that the positions are within bounds.

core/approximate-sdf.cpp

@@ -0,0 +1,216 @@
+
+#include "approximate-sdf.h"
+
+#include <cmath>
+#include <queue>
+#include "arithmetics.hpp"
+
+#define ESTSDF_MAX_DIST 1e24f // Cannot be FLT_MAX because it might be divided by range, which could be < 1
+
+namespace msdfgen {
+
+void approximateSDF(const BitmapRef<float, 1> &output, const Shape &shape, const SDFTransformation &transformation) {
+    struct Entry {
+        float absDist;
+        int bitmapX, bitmapY;
+        Point2 nearPoint;
+
+        bool operator<(const Entry &other) const {
+            return absDist > other.absDist;
+        }
+    } entry;
+
+    float *firstRow = output.pixels;
+    ptrdiff_t stride = output.width;
+    if (shape.inverseYAxis) {
+        firstRow += (output.height-1)*stride;
+        stride = -stride;
+    }
+    #define ESTSDF_PIXEL_AT(x, y) ((firstRow+(y)*stride)[x])
+
+    for (float *p = output.pixels, *end = output.pixels+output.width*output.height; p < end; ++p)
+        *p = -ESTSDF_MAX_DIST;
+
+    Vector2 invScale = transformation.unprojectVector(Vector2(1));
+    DistanceMapping invDistanceMapping = transformation.distanceMapping.inverse();
+    float dLimit = float(max(fabs(invDistanceMapping(0)), fabs(invDistanceMapping(1))));
+    std::priority_queue<Entry> queue;
+    double x[3], y[3];
+    int dx[3], dy[3];
+
+    // Horizontal scanlines
+    for (int bitmapY = 0; bitmapY < output.height; ++bitmapY) {
+        float *row = firstRow+bitmapY*stride;
+        double y = transformation.unprojectY(bitmapY+.5);
+        entry.bitmapY = bitmapY;
+        for (std::vector<Contour>::const_iterator contour = shape.contours.begin(); contour != shape.contours.end(); ++contour) {
+            for (std::vector<EdgeHolder>::const_iterator edge = contour->edges.begin(); edge != contour->edges.end(); ++edge) {
+                int n = (*edge)->horizontalScanlineIntersections(x, dy, y);
+                for (int i = 0; i < n; ++i) {
+                    double bitmapX = transformation.projectX(x[i]);
+                    double bitmapX0 = floor(bitmapX-.5)+.5;
+                    double bitmapX1 = bitmapX0+1;
+                    if (bitmapX1 > 0 && bitmapX0 < output.width) {
+                        float sd0 = float(dy[i]*invScale.x*(bitmapX0-bitmapX));
+                        float sd1 = float(dy[i]*invScale.x*(bitmapX1-bitmapX));
+                        if (sd0 == 0.f) {
+                            if (sd1 == 0.f)
+                                continue;
+                            sd0 = -.000001f*float(sign(sd1));
+                        }
+                        if (sd1 == 0.f)
+                            sd1 = -.000001f*float(sign(sd0));
+                        if (bitmapX0 > 0) {
+                            entry.absDist = fabsf(sd0);
+                            entry.bitmapX = int(bitmapX0);
+                            float &sd = row[entry.bitmapX];
+                            if (entry.absDist < fabsf(sd)) {
+                                sd = sd0;
+                                entry.nearPoint = Point2(x[i], y);
+                                queue.push(entry);
+                            } else if (sd == -sd0)
+                                sd = -ESTSDF_MAX_DIST;
+                        }
+                        if (bitmapX1 < output.width) {
+                            entry.absDist = fabsf(sd1);
+                            entry.bitmapX = int(bitmapX1);
+                            float &sd = row[entry.bitmapX];
+                            if (entry.absDist < fabsf(sd)) {
+                                sd = sd1;
+                                entry.nearPoint = Point2(x[i], y);
+                                queue.push(entry);
+                            } else if (sd == -sd1)
+                                sd = -ESTSDF_MAX_DIST;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // Bake in distance signs
+    for (int y = 0; y < output.height; ++y) {
+        float *row = firstRow+y*stride;
+        int x = 0;
+        for (; x < output.width && row[x] == -ESTSDF_MAX_DIST; ++x);
+        if (x < output.width) {
+            bool flip = row[x] > 0;
+            if (flip) {
+                for (int i = 0; i < x; ++i)
+                    row[i] = ESTSDF_MAX_DIST;
+            }
+            for (; x < output.width; ++x) {
+                if (row[x] != -ESTSDF_MAX_DIST)
+                    flip = row[x] > 0;
+                else if (flip)
+                    row[x] = ESTSDF_MAX_DIST;
+            }
+        }
+    }
+
+    // Vertical scanlines
+    for (int bitmapX = 0; bitmapX < output.width; ++bitmapX) {
+        double x = transformation.unprojectX(bitmapX+.5);
+        entry.bitmapX = bitmapX;
+        for (std::vector<Contour>::const_iterator contour = shape.contours.begin(); contour != shape.contours.end(); ++contour) {
+            for (std::vector<EdgeHolder>::const_iterator edge = contour->edges.begin(); edge != contour->edges.end(); ++edge) {
+                int n = (*edge)->verticalScanlineIntersections(y, dx, x);
+                for (int i = 0; i < n; ++i) {
+                    double bitmapY = transformation.projectY(y[i]);
+                    double bitmapY0 = floor(bitmapY-.5)+.5;
+                    double bitmapY1 = bitmapY0+1;
+                    if (bitmapY0 > 0 && bitmapY1 < output.height) {
+                        float sd0 = float(dx[i]*invScale.y*(bitmapY-bitmapY0));
+                        float sd1 = float(dx[i]*invScale.y*(bitmapY-bitmapY1));
+                        if (sd0 == 0.f) {
+                            if (sd1 == 0.f)
+                                continue;
+                            sd0 = -.000001f*float(sign(sd1));
+                        }
+                        if (sd1 == 0.f)
+                            sd1 = -.000001f*float(sign(sd0));
+                        if (bitmapY0 > 0) {
+                            entry.absDist = fabsf(sd0);
+                            entry.bitmapY = int(bitmapY0);
+                            float &sd = ESTSDF_PIXEL_AT(bitmapX, entry.bitmapY);
+                            if (entry.absDist < fabsf(sd)) {
+                                sd = sd0;
+                                entry.nearPoint = Point2(x, y[i]);
+                                queue.push(entry);
+                            }
+                        }
+                        if (bitmapY1 < output.height) {
+                            entry.absDist = fabsf(sd1);
+                            entry.bitmapY = int(bitmapY1);
+                            float &sd = ESTSDF_PIXEL_AT(bitmapX, entry.bitmapY);
+                            if (entry.absDist < fabsf(sd)) {
+                                sd = sd1;
+                                entry.nearPoint = Point2(x, y[i]);
+                                queue.push(entry);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    if (queue.empty())
+        return;
+
+    while (!queue.empty()) {
+        Entry entry = queue.top();
+        queue.pop();
+        Entry newEntry = entry;
+        newEntry.bitmapX = entry.bitmapX-1;
+        if (newEntry.bitmapX >= 0) {
+            float &sd = ESTSDF_PIXEL_AT(newEntry.bitmapX, newEntry.bitmapY);
+            if (fabsf(sd) == ESTSDF_MAX_DIST) {
+                Point2 shapeCoord = transformation.unproject(Point2(newEntry.bitmapX+.5, newEntry.bitmapY+.5));
+                newEntry.absDist = float((shapeCoord-entry.nearPoint).length());
+                sd = float(sign(sd))*newEntry.absDist;
+                if (newEntry.absDist < dLimit)
+                    queue.push(newEntry);
+            }
+        }
+        newEntry.bitmapX = entry.bitmapX+1;
+        if (newEntry.bitmapX < output.width) {
+            float &sd = ESTSDF_PIXEL_AT(newEntry.bitmapX, newEntry.bitmapY);
+            if (fabsf(sd) == ESTSDF_MAX_DIST) {
+                Point2 shapeCoord = transformation.unproject(Point2(newEntry.bitmapX+.5, newEntry.bitmapY+.5));
+                newEntry.absDist = float((shapeCoord-entry.nearPoint).length());
+                sd = float(sign(sd))*newEntry.absDist;
+                if (newEntry.absDist < dLimit)
+                    queue.push(newEntry);
+            }
+        }
+        newEntry.bitmapX = entry.bitmapX;
+        newEntry.bitmapY = entry.bitmapY-1;
+        if (newEntry.bitmapY >= 0) {
+            float &sd = ESTSDF_PIXEL_AT(newEntry.bitmapX, newEntry.bitmapY);
+            if (fabsf(sd) == ESTSDF_MAX_DIST) {
+                Point2 shapeCoord = transformation.unproject(Point2(newEntry.bitmapX+.5, newEntry.bitmapY+.5));
+                newEntry.absDist = float((shapeCoord-entry.nearPoint).length());
+                sd = float(sign(sd))*newEntry.absDist;
+                if (newEntry.absDist < dLimit)
+                    queue.push(newEntry);
+            }
+        }
+        newEntry.bitmapY = entry.bitmapY+1;
+        if (newEntry.bitmapY < output.height) {
+            float &sd = ESTSDF_PIXEL_AT(newEntry.bitmapX, newEntry.bitmapY);
+            if (fabsf(sd) == ESTSDF_MAX_DIST) {
+                Point2 shapeCoord = transformation.unproject(Point2(newEntry.bitmapX+.5, newEntry.bitmapY+.5));
+                newEntry.absDist = float((shapeCoord-entry.nearPoint).length());
+                sd = float(sign(sd))*newEntry.absDist;
+                if (newEntry.absDist < dLimit)
+                    queue.push(newEntry);
+            }
+        }
+    }
+
+    for (float *p = output.pixels, *end = output.pixels+output.width*output.height; p < end; ++p)
+        *p = transformation.distanceMapping(*p);
+}
+
+}

core/approximate-sdf.h

@@ -0,0 +1,13 @@
+
+#pragma once
+
+#include "SDFTransformation.h"
+#include "Shape.h"
+#include "BitmapRef.hpp"
+
+namespace msdfgen {
+
+// Fast SDF approximation (out of range values not computed)
+void approximateSDF(const BitmapRef<float, 1> &output, const Shape &shape, const SDFTransformation &transformation);
+
+}

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 {
 
@@ -184,6 +187,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];
@@ -270,7 +335,21 @@ 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 {
+    return horizontalScanlineIntersections(x, dy, y);
+}
+
+int QuadraticSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
+    return horizontalScanlineIntersections(x, dy, y);
+}
+
+int CubicSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
+    return horizontalScanlineIntersections(x, dy, y);
+}
+
+int LinearSegment::horizontalScanlineIntersections(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);
         x[0] = mix(p[0].x, p[1].x, param);
@@ -280,7 +359,17 @@ int LinearSegment::scanlineIntersections(double x[3], int dy[3], double y) const
     return 0;
 }
 
-int QuadraticSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
+int LinearSegment::verticalScanlineIntersections(double y[3], int dx[3], double x) const {
+    if ((x >= p[0].x && x < p[1].x) || (x >= p[1].x && x < p[0].x)) {
+        double param = (x-p[0].x)/(p[1].x-p[0].x);
+        y[0] = mix(p[0].y, p[1].y, param);
+        dx[0] = sign(p[1].x-p[0].x);
+        return 1;
+    }
+    return 0;
+}
+
+int QuadraticSegment::horizontalScanlineIntersections(double x[3], int dy[3], double y) const {
     int total = 0;
     int nextDY = y > p[0].y ? 1 : -1;
     x[total] = p[0].x;
@@ -334,7 +423,61 @@ int QuadraticSegment::scanlineIntersections(double x[3], int dy[3], double y) co
     return total;
 }
 
-int CubicSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
+int QuadraticSegment::verticalScanlineIntersections(double y[3], int dx[3], double x) const {
+    int total = 0;
+    int nextDX = x > p[0].x ? 1 : -1;
+    y[total] = p[0].y;
+    if (p[0].x == x) {
+        if (p[0].x < p[1].x || (p[0].x == p[1].x && p[0].x < p[2].x))
+            dx[total++] = 1;
+        else
+            nextDX = 1;
+    }
+    {
+        Vector2 ab = p[1]-p[0];
+        Vector2 br = p[2]-p[1]-ab;
+        double t[2];
+        int solutions = solveQuadratic(t, br.x, 2*ab.x, p[0].x-x);
+        // Sort solutions
+        double tmp;
+        if (solutions >= 2 && t[0] > t[1])
+            tmp = t[0], t[0] = t[1], t[1] = tmp;
+        for (int i = 0; i < solutions && total < 2; ++i) {
+            if (t[i] >= 0 && t[i] <= 1) {
+                y[total] = p[0].y+2*t[i]*ab.y+t[i]*t[i]*br.y;
+                if (nextDX*(ab.x+t[i]*br.x) >= 0) {
+                    dx[total++] = nextDX;
+                    nextDX = -nextDX;
+                }
+            }
+        }
+    }
+    if (p[2].x == x) {
+        if (nextDX > 0 && total > 0) {
+            --total;
+            nextDX = -1;
+        }
+        if ((p[2].x < p[1].x || (p[2].x == p[1].x && p[2].x < p[0].x)) && total < 2) {
+            y[total] = p[2].y;
+            if (nextDX < 0) {
+                dx[total++] = -1;
+                nextDX = 1;
+            }
+        }
+    }
+    if (nextDX != (x >= p[2].x ? 1 : -1)) {
+        if (total > 0)
+            --total;
+        else {
+            if (fabs(p[2].x-x) < fabs(p[0].x-x))
+                y[total] = p[2].y;
+            dx[total++] = nextDX;
+        }
+    }
+    return total;
+}
+
+int CubicSegment::horizontalScanlineIntersections(double x[3], int dy[3], double y) const {
     int total = 0;
     int nextDY = y > p[0].y ? 1 : -1;
     x[total] = p[0].x;
@@ -396,6 +539,68 @@ int CubicSegment::scanlineIntersections(double x[3], int dy[3], double y) const
     return total;
 }
 
+int CubicSegment::verticalScanlineIntersections(double y[3], int dx[3], double x) const {
+    int total = 0;
+    int nextDX = x > p[0].x ? 1 : -1;
+    y[total] = p[0].y;
+    if (p[0].x == x) {
+        if (p[0].x < p[1].x || (p[0].x == p[1].x && (p[0].x < p[2].x || (p[0].x == p[2].x && p[0].x < p[3].x))))
+            dx[total++] = 1;
+        else
+            nextDX = 1;
+    }
+    {
+        Vector2 ab = p[1]-p[0];
+        Vector2 br = p[2]-p[1]-ab;
+        Vector2 as = (p[3]-p[2])-(p[2]-p[1])-br;
+        double t[3];
+        int solutions = solveCubic(t, as.x, 3*br.x, 3*ab.x, p[0].x-x);
+        // Sort solutions
+        double tmp;
+        if (solutions >= 2) {
+            if (t[0] > t[1])
+                tmp = t[0], t[0] = t[1], t[1] = tmp;
+            if (solutions >= 3 && t[1] > t[2]) {
+                tmp = t[1], t[1] = t[2], t[2] = tmp;
+                if (t[0] > t[1])
+                    tmp = t[0], t[0] = t[1], t[1] = tmp;
+            }
+        }
+        for (int i = 0; i < solutions && total < 3; ++i) {
+            if (t[i] >= 0 && t[i] <= 1) {
+                y[total] = p[0].y+3*t[i]*ab.y+3*t[i]*t[i]*br.y+t[i]*t[i]*t[i]*as.y;
+                if (nextDX*(ab.x+2*t[i]*br.x+t[i]*t[i]*as.x) >= 0) {
+                    dx[total++] = nextDX;
+                    nextDX = -nextDX;
+                }
+            }
+        }
+    }
+    if (p[3].x == x) {
+        if (nextDX > 0 && total > 0) {
+            --total;
+            nextDX = -1;
+        }
+        if ((p[3].x < p[2].x || (p[3].x == p[2].x && (p[3].x < p[1].x || (p[3].x == p[1].x && p[3].x < p[0].x)))) && total < 3) {
+            y[total] = p[3].y;
+            if (nextDX < 0) {
+                dx[total++] = -1;
+                nextDX = 1;
+            }
+        }
+    }
+    if (nextDX != (x >= p[3].x ? 1 : -1)) {
+        if (total > 0)
+            --total;
+        else {
+            if (fabs(p[3].x-x) < fabs(p[0].x-x))
+                y[total] = p[3].y;
+            dx[total++] = nextDX;
+        }
+    }
+    return total;
+}
+
 static void pointBounds(Point2 p, double &l, double &b, double &r, double &t) {
     if (p.x < l) l = p.x;
     if (p.y < b) b = p.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 {
 
@@ -41,6 +37,9 @@ public:
     virtual void distanceToPerpendicularDistance(SignedDistance &distance, Point2 origin, double param) const;
     /// Outputs a list of (at most three) intersections (their X coordinates) with an infinite horizontal scanline at y and returns how many there are.
     virtual int scanlineIntersections(double x[3], int dy[3], double y) const = 0;
+    virtual int horizontalScanlineIntersections(double x[3], int dy[3], double y) const = 0;
+    /// Outputs a list of (at most three) intersections (their Y coordinates) with an infinite vertical scanline at x and returns how many there are.
+    virtual int verticalScanlineIntersections(double y[3], int dx[3], double x) const = 0;
     /// Adjusts the bounding box to fit the edge segment.
     virtual void bound(double &l, double &b, double &r, double &t) const = 0;
 
@@ -75,6 +74,8 @@ public:
     double length() const;
     SignedDistance signedDistance(Point2 origin, double &param) const;
     int scanlineIntersections(double x[3], int dy[3], double y) const;
+    int horizontalScanlineIntersections(double x[3], int dy[3], double y) const;
+    int verticalScanlineIntersections(double y[3], int dx[3], double x) const;
     void bound(double &l, double &b, double &r, double &t) const;
 
     void reverse();
@@ -104,6 +105,8 @@ public:
     double length() const;
     SignedDistance signedDistance(Point2 origin, double &param) const;
     int scanlineIntersections(double x[3], int dy[3], double y) const;
+    int horizontalScanlineIntersections(double x[3], int dy[3], double y) const;
+    int verticalScanlineIntersections(double y[3], int dx[3], double x) const;
     void bound(double &l, double &b, double &r, double &t) const;
 
     void reverse();
@@ -134,6 +137,8 @@ public:
     Vector2 directionChange(double param) const;
     SignedDistance signedDistance(Point2 origin, double &param) const;
     int scanlineIntersections(double x[3], int dy[3], double y) const;
+    int horizontalScanlineIntersections(double x[3], int dy[3], double y) const;
+    int verticalScanlineIntersections(double y[3], int dx[3], double x) const;
     void bound(double &l, double &b, double &r, double &t) const;
 
     void reverse();

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)) {

msdfgen.h

@@ -33,6 +33,7 @@
 #include "core/msdf-error-correction.h"
 #include "core/render-sdf.h"
 #include "core/rasterization.h"
+#include "core/approximate-sdf.h"
 #include "core/sdf-error-estimation.h"
 #include "core/save-bmp.h"
 #include "core/save-tiff.h"