Improved defines

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

all-in-one/.gitignore

@@ -0,0 +1,3 @@
+msdfgen.h
+msdfgen.cpp
+build/

all-in-one/CMakeLists.txt

@@ -0,0 +1,64 @@
+
+cmake_minimum_required(VERSION 3.15)
+
+option(MSDFGEN_USE_VCPKG "Use vcpkg package manager to link project dependencies" ON)
+option(MSDFGEN_USE_FREETYPE "Build with the FreeType library" ON)
+
+if(MSDFGEN_USE_VCPKG AND MSDFGEN_USE_FREETYPE)
+    # Make sure that vcpkg toolchain file is set
+    if(NOT CMAKE_TOOLCHAIN_FILE)
+        if(DEFINED ENV{VCPKG_ROOT})
+            set(CMAKE_TOOLCHAIN_FILE "$ENV{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake")
+        else()
+            message(SEND_ERROR "Vcpkg toolchain not configured. Either set VCPKG_ROOT environment variable or pass -DCMAKE_TOOLCHAIN_FILE=VCPKG_ROOT/scripts/buildsystems/vcpkg.cmake to cmake")
+        endif()
+    endif()
+    # Default to statically linked vcpkg triplet on Windows
+    if(WIN32 AND NOT VCPKG_TARGET_TRIPLET AND NOT MSDFGEN_DYNAMIC_RUNTIME)
+        if(CMAKE_GENERATOR_PLATFORM MATCHES "64$" AND NOT CMAKE_GENERATOR_PLATFORM STREQUAL "ARM64")
+            set(VCPKG_TARGET_TRIPLET "x64-windows-static")
+        elseif(CMAKE_GENERATOR_PLATFORM MATCHES "32$" OR CMAKE_GENERATOR_PLATFORM STREQUAL "x86")
+            set(VCPKG_TARGET_TRIPLET "x86-windows-static")
+        else()
+            if(CMAKE_GENERATOR_PLATFORM)
+                message(WARNING "Vcpkg triplet not explicitly specified and could not be deduced. Recommend using -DVCPKG_TARGET_TRIPLET=x64-windows-static or similar")
+            else()
+                message(WARNING "Vcpkg triplet not explicitly specified and could not be deduced. Recommend using -A to explicitly select platform (Win32 or x64)")
+            endif()
+        endif()
+    endif()
+endif()
+
+project(msdfgen LANGUAGES CXX)
+
+if(MAX_WARNING_LEVEL)
+    if (MSVC)
+        add_compile_options(/W4)
+    else()
+        add_compile_options(-Wall -Wextra -Wpedantic)
+    endif()
+endif()
+
+find_package(Python3 REQUIRED COMPONENTS Interpreter)
+if(MSDFGEN_USE_FREETYPE AND NOT TARGET Freetype::Freetype)
+    find_package(Freetype REQUIRED)
+endif()
+
+execute_process(COMMAND ${Python3_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/generate.py")
+
+add_library(msdfgen "${CMAKE_CURRENT_SOURCE_DIR}/msdfgen.h" "${CMAKE_CURRENT_SOURCE_DIR}/msdfgen.cpp")
+set_property(TARGET msdfgen PROPERTY MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>")
+target_compile_features(msdfgen PUBLIC cxx_std_11)
+if(MSDFGEN_USE_FREETYPE)
+    target_link_libraries(msdfgen PRIVATE Freetype::Freetype)
+else()
+    target_compile_definitions(msdfgen PUBLIC MSDFGEN_NO_FREETYPE)
+endif()
+
+add_executable(msdfgen-test "${CMAKE_CURRENT_SOURCE_DIR}/test.cpp")
+set_property(TARGET msdfgen-test PROPERTY MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>")
+target_link_libraries(msdfgen-test PRIVATE msdfgen)
+set_property(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT msdfgen-test)
+if(MSDFGEN_USE_FREETYPE)
+    target_link_libraries(msdfgen-test PRIVATE Freetype::Freetype)
+endif()

all-in-one/generate.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+
+import os
+import re
+
+rootDir = os.path.join(os.path.dirname(__file__), '..')
+
+sourceList = [
+    'core/base.h',
+    'core/arithmetics.hpp',
+    'core/equation-solver.h',
+    'core/equation-solver.cpp',
+    'core/Vector2.hpp',
+    'core/pixel-conversion.hpp',
+    'core/BitmapRef.hpp',
+    'core/Bitmap.h',
+    'core/Bitmap.hpp',
+    'core/Range.hpp',
+    'core/Projection.h',
+    'core/Projection.cpp',
+    'core/DistanceMapping.h',
+    'core/DistanceMapping.cpp',
+    'core/SDFTransformation.h',
+    'core/SignedDistance.hpp',
+    'core/Scanline.h',
+    'core/Scanline.cpp',
+    'core/EdgeColor.h',
+    'core/bezier-solver.hpp',
+    'core/edge-segments.h',
+    'core/edge-segments.cpp',
+    'core/EdgeHolder.h',
+    'core/EdgeHolder.cpp',
+    'core/Contour.h',
+    'core/Contour.cpp',
+    'core/Shape.h',
+    'core/Shape.cpp',
+    'core/bitmap-interpolation.hpp',
+    'core/edge-coloring.h',
+    'core/edge-coloring.cpp',
+    'core/edge-selectors.h',
+    'core/edge-selectors.cpp',
+    'core/contour-combiners.h',
+    'core/contour-combiners.cpp',
+    'core/ShapeDistanceFinder.h',
+    'core/ShapeDistanceFinder.hpp',
+    'core/approximate-sdf.h',
+    'core/approximate-sdf.cpp',
+    'core/generator-config.h',
+    'core/msdf-error-correction.h',
+    'core/msdf-error-correction.cpp',
+    'core/MSDFErrorCorrection.h',
+    'core/MSDFErrorCorrection.cpp',
+    'core/msdfgen.cpp',
+    'ext/import-font.h',
+    'ext/import-font.cpp',
+    'ext/resolve-shape-geometry.h',
+    'ext/resolve-shape-geometry.cpp',
+    'ext/import-svg.h',
+    'ext/import-svg.cpp',
+    'msdfgen.h'
+]
+
+header = """
+#pragma once
+
+#define MSDFGEN_USE_CPP11
+#ifndef MSDFGEN_NO_FREETYPE
+#define MSDFGEN_USE_FREETYPE
+#define MSDFGEN_DISABLE_VARIABLE_FONTS
+#endif
+#ifndef MSDFGEN_ENABLE_SVG
+#define MSDFGEN_DISABLE_SVG
+#endif
+
+#include <cstddef>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <vector>
+"""
+
+source = """
+#ifndef _CRT_SECURE_NO_WARNINGS
+#define _CRT_SECURE_NO_WARNINGS
+#endif
+
+#include "msdfgen.h"
+
+#include <queue>
+
+#ifdef MSDFGEN_USE_FREETYPE
+#include <ft2build.h>
+#include FT_FREETYPE_H
+#include FT_OUTLINE_H
+#ifndef MSDFGEN_DISABLE_VARIABLE_FONTS
+#include FT_MULTIPLE_MASTERS_H
+#endif
+#endif
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wshadow"
+#elif defined(_MSC_VER)
+#pragma warning(push)
+#pragma warning(disable : 4456 4457 4458 6246)
+#endif
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626433832795
+#endif
+"""
+
+namespaceStart = """
+#ifdef MSDFGEN_PARENT_NAMESPACE
+namespace MSDFGEN_PARENT_NAMESPACE {
+#endif
+"""
+
+namespaceEnd = """
+#ifdef MSDFGEN_PARENT_NAMESPACE
+} // namespace MSDFGEN_PARENT_NAMESPACE
+#endif
+"""
+
+sourceAppendix = """
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC diagnostic pop
+#elif defined(_MSC_VER)
+#pragma warning(pop)
+#endif
+"""
+
+header += namespaceStart
+source += namespaceStart
+
+with open(os.path.join(rootDir, 'LICENSE.txt'), 'r') as file:
+    license = file.read()
+license = '\n'.join([' * '+line for line in license.strip().split('\n')])
+
+for filename in sourceList:
+    with open(os.path.join(rootDir, *filename.split('/')), 'r') as file:
+        src = file.read()
+    src = '\n'.join([line for line in src.split('\n') if not re.match(r'^\s*#(include\s.*|pragma\s+once|define\s+(_CRT_SECURE_NO_WARNINGS|_USE_MATH_DEFINES))\s*$', line)])
+    if filename.startswith('ext/import-font.'):
+        src = '#ifdef MSDFGEN_USE_FREETYPE\n\n'+src+'\n\n#endif\n\n'
+    if filename.endswith('.h') or filename.endswith('.hpp'):
+        header += '\n\n'+src
+    if filename.endswith('.cpp'):
+        source += '\n\n'+src
+
+header = '\n'+re.sub(r'\n{3,}', '\n\n', re.sub(r'}\s*namespace\s+msdfgen\s*{', '', re.sub(r'\/\*[^\*].*?\*\/', '', header, flags=re.DOTALL))).strip()+'\n'
+source = '\n'+re.sub(r'\n{3,}', '\n\n', re.sub(r'}\s*namespace\s+msdfgen\s*{', '', re.sub(r'\/\*[^\*].*?\*\/', '', source, flags=re.DOTALL))).strip()+'\n'
+header += namespaceEnd
+source += namespaceEnd
+
+header = """
+/*
+ * MULTI-CHANNEL SIGNED DISTANCE FIELD GENERATOR
+ * ---------------------------------------------
+ * A utility by Viktor Chlumsky, (c) 2014 - 2024
+ * https://github.com/Chlumsky/msdfgen
+ * Published under the MIT license
+ *
+ * The technique used to generate multi-channel distance fields in this code
+ * was developed by Viktor Chlumsky in 2014 for his master's thesis,
+ * "Shape Decomposition for Multi-Channel Distance Fields". It provides improved
+ * quality of sharp corners in glyphs and other 2D shapes compared to monochrome
+ * distance fields. To reconstruct an image of the shape, apply the median of three
+ * operation on the triplet of sampled signed distance values.
+ *
+ */
+"""+header
+
+source = """
+/*
+ * MULTI-CHANNEL SIGNED DISTANCE FIELD GENERATOR
+ * ---------------------------------------------
+ * https://github.com/Chlumsky/msdfgen
+ *
+"""+license+"""
+ *
+ */
+"""+source+sourceAppendix
+
+with open(os.path.join(os.path.dirname(__file__), 'msdfgen.h'), 'w') as file:
+    file.write(header)
+with open(os.path.join(os.path.dirname(__file__), 'msdfgen.cpp'), 'w') as file:
+    file.write(source)

all-in-one/test.cpp

@@ -0,0 +1,52 @@
+
+#include <cstdio>
+#include "msdfgen.h"
+
+using namespace msdfgen;
+
+static bool saveRGBA(const char *filename, const BitmapConstRef<float, 4> &bitmap) {
+    if (FILE *f = fopen(filename, "wb")) {
+        char header[12];
+        header[0] = 'R', header[1] = 'G', header[2] = 'B', header[3] = 'A';
+        header[4] = char(bitmap.width>>24);
+        header[5] = char(bitmap.width>>16);
+        header[6] = char(bitmap.width>>8);
+        header[7] = char(bitmap.width);
+        header[8] = char(bitmap.height>>24);
+        header[9] = char(bitmap.height>>16);
+        header[10] = char(bitmap.height>>8);
+        header[11] = char(bitmap.height);
+        fwrite(header, 1, 12, f);
+        Bitmap<byte, 4> byteBitmap(bitmap.width, bitmap.height);
+        byte *d = (byte *) byteBitmap;
+        for (const float *s = bitmap.pixels, *end = bitmap.pixels+4*bitmap.width*bitmap.height; s < end; ++d, ++s)
+            *d = pixelFloatToByte(*s);
+        fwrite((const byte *) byteBitmap, 1, 4*bitmap.width*bitmap.height, f);
+        fclose(f);
+        return true;
+    }
+    return false;
+}
+
+int main() {
+#ifdef MSDFGEN_USE_FREETYPE
+    bool success = false;
+    if (FreetypeHandle *ft = initializeFreetype()) {
+        if (FontHandle *font = loadFont(ft, "C:\\Windows\\Fonts\\arialbd.ttf")) {
+            Shape shape;
+            if (loadGlyph(shape, font, 'A')) {
+                shape.normalize();
+                edgeColoringByDistance(shape, 3.0);
+                Bitmap<float, 4> msdf(32, 32);
+                generateMTSDF(msdf, shape, 4.0, 1.0, Vector2(4.0, 4.0));
+                success = saveRGBA("output.rgba", msdf);
+            }
+            destroyFont(font);
+        }
+        deinitializeFreetype(ft);
+    }
+    return success ? 0 : 1;
+#else
+    return 0;
+#endif
+}

all-in-one/vcpkg.json

@@ -0,0 +1,8 @@
+{
+    "$schema": "https://raw.githubusercontent.com/microsoft/vcpkg-tool/master/docs/vcpkg.schema.json",
+    "name": "msdfgen-all-in-one",
+    "version": "1.10.0",
+    "dependencies": [
+        "freetype"
+    ]
+}

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"