Updated API of approximateSDF

.gitignore

@@ -26,6 +26,7 @@ output.png
 render.png
 CMakeUserPresets.json
 out/
+/all-in-one/
 /build_xcode/
 /cmake-gen.bat
 /line-ending-check.bat

CHANGELOG.md

@@ -1,4 +1,23 @@
 
+## Version 1.12 (2024-05-18)
+
+- Added the possibility to specify asymmetrical distance range (`-arange`, `-apxrange`)
+- Added the ability to export the shape into an SVG file (`-exportsvg`)
+- Edge coloring no longer colors smooth contours white, which has been found suboptimal
+- Fixed incorrect scaling of font glyph coordinates. To preserve backwards compatibility, the user has to enable the fix with an explicit additional argument:
+    - `-emnormalize` in standalone, `FONT_SCALING_EM_NORMALIZED` in API for coordinates in ems
+    - `-noemnormalize` in standalone, `FONT_SCALING_NONE` in API for raw integer coordinates
+    - The default (backwards-compatible) behavior will change in a future version; a warning will be displayed if neither option is set
+- Added two new developer-friendly export image formats: RGBA and FL32
+- `-size` parameter renamed to `-dimensions` for clarity (old one will be kept for compatibility)
+- `generate*SDF` functions now combine projection and range into a single argument (`SDFTransformation`)
+- Conversion of floating point color values to 8-bit integers adjusted to match graphics hardware
+- Improved edge deconvergence procedure and made sure that calling `Shape::normalize` a second time has no effect
+- Fixed certain edge cases where Skia geometry preprocessing wouldn't make the geometry fully compliant
+- The term "perpendicular distance" now used instead of "pseudo-distance" (PSDF instead of PseudoSDF in API)
+- Fixed a bug in `savePng` where `fclose` could be called on null pointer
+- Minor code improvements
+
 ## Version 1.11 (2023-11-11)
 
 - Reworked SVG parser, which now supports multiple paths and other shapes - requires Skia

LICENSE.txt

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2016 - 2024 Viktor Chlumsky
+Copyright (c) 2014 - 2024 Viktor Chlumsky
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

build-release.bat

@@ -0,0 +1,96 @@
+@echo off
+setlocal
+pushd "%~dp0"
+
+if "%VCVARSALL%"=="" set "VCVARSALL=C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat"
+if "%VCTOOLSET%"=="" set "VCTOOLSET=VC143"
+
+if not exist "%VCVARSALL%" (
+    echo vcvarsall.bat not found, set VCVARSALL to its path
+    exit /b 1
+)
+
+if "%1"=="" (
+    echo No version specified, using current
+    set "version=current"
+) else (
+    set "version=%1"
+)
+set "builddir=%~dp0\build\release-%version%"
+
+if exist "%builddir%" (
+    echo Deleting contents of %builddir%
+    rd /s /q "%builddir%"
+)
+
+cmake --preset win32 -B "%builddir%\build-win32"
+cmake --preset win32-omp -B "%builddir%\build-win32-omp"
+cmake --preset win64 -B "%builddir%\build-win64"
+cmake --preset win64-omp -B "%builddir%\build-win64-omp"
+
+cmake --build "%builddir%\build-win32" --config Release
+cmake --build "%builddir%\build-win32-omp" --config Release
+cmake --build "%builddir%\build-win64" --config Release
+cmake --build "%builddir%\build-win64-omp" --config Release
+
+mkdir "%builddir%\rel-win32\msdfgen"
+mkdir "%builddir%\rel-win32-omp\msdfgen"
+mkdir "%builddir%\rel-win64\msdfgen"
+mkdir "%builddir%\rel-win64-omp\msdfgen"
+
+copy "%builddir%\build-win32\Release\msdfgen.exe" "%builddir%\rel-win32\msdfgen"
+copy "%builddir%\build-win32-omp\Release\msdfgen.exe" "%builddir%\rel-win32-omp\msdfgen"
+copy "%builddir%\build-win64\Release\msdfgen.exe" "%builddir%\rel-win64\msdfgen"
+copy "%builddir%\build-win64-omp\Release\msdfgen.exe" "%builddir%\rel-win64-omp\msdfgen"
+
+echo msdfgen.exe -defineshape "{ 1471,0; 1149,0; 1021,333; 435,333; 314,0; 0,0; 571,1466; 884,1466; # }{ 926,580; 724,1124; 526,580; # }" -size 16 16 -autoframe -testrender render.png 1024 1024 > "%builddir%\example.bat"
+
+copy "%builddir%\example.bat" "%builddir%\rel-win32\msdfgen"
+copy "%builddir%\example.bat" "%builddir%\rel-win32-omp\msdfgen"
+copy "%builddir%\example.bat" "%builddir%\rel-win64\msdfgen"
+copy "%builddir%\example.bat" "%builddir%\rel-win64-omp\msdfgen"
+
+call "%VCVARSALL%" x64
+
+set "omp32dll=%VCToolsRedistDir%\x86\Microsoft.%VCTOOLSET%.OPENMP\vcomp140.dll"
+set "omp64dll=%VCToolsRedistDir%\x64\Microsoft.%VCTOOLSET%.OPENMP\vcomp140.dll"
+
+if not exist "%omp32dll%" (
+    echo vcomp140.dll [x86] not found, make sure to set VCTOOLSET or update this script
+    exit /b 1
+)
+if not exist "%omp64dll%" (
+    echo vcomp140.dll [x64] not found, make sure to set VCTOOLSET or update this script
+    exit /b 1
+)
+
+copy "%omp32dll%" "%builddir%\rel-win32-omp\msdfgen"
+copy "%omp64dll%" "%builddir%\rel-win64-omp\msdfgen"
+
+if not exist "C:\Program Files\7-Zip\7z.exe" (
+    echo 7-Zip not found, you have to package it manually
+    exit /b 0
+)
+
+pushd "%builddir%\rel-win32"
+"C:\Program Files\7-Zip\7z.exe" a "..\msdfgen-%version%-win32.zip" msdfgen
+cd msdfgen
+call example.bat
+popd
+pushd "%builddir%\rel-win32-omp"
+"C:\Program Files\7-Zip\7z.exe" a "..\msdfgen-%version%-win32-openmp.zip" msdfgen
+cd msdfgen
+call example.bat
+popd
+pushd "%builddir%\rel-win64"
+"C:\Program Files\7-Zip\7z.exe" a "..\msdfgen-%version%-win64.zip" msdfgen
+cd msdfgen
+call example.bat
+popd
+pushd "%builddir%\rel-win64-omp"
+"C:\Program Files\7-Zip\7z.exe" a "..\msdfgen-%version%-win64-openmp.zip" msdfgen
+cd msdfgen
+call example.bat
+popd
+
+popd

core/MSDFErrorCorrection.cpp

@@ -89,6 +89,8 @@ public:
     bool protectedFlag;
     inline ShapeDistanceChecker(const BitmapConstRef<float, N> &sdf, const Shape &shape, const Projection &projection, DistanceMapping distanceMapping, double minImproveRatio) : distanceFinder(shape), sdf(sdf), distanceMapping(distanceMapping), minImproveRatio(minImproveRatio) {
         texelSize = projection.unprojectVector(Vector2(1));
+        if (shape.inverseYAxis)
+            texelSize.y = -texelSize.y;
     }
     inline ArtifactClassifier classifier(const Vector2 &direction, double span) {
         return ArtifactClassifier(this, direction, span);

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"