C cubic curve distance update

CMakeLists.txt

@@ -67,7 +67,7 @@ if(MSDFGEN_USE_VCPKG)
 endif()
 
 # Version is specified in vcpkg.json
-project(msdfgen VERSION ${MSDFGEN_VERSION} LANGUAGES CXX)
+project(msdfgen VERSION ${MSDFGEN_VERSION} LANGUAGES C CXX)
 
 if(MSDFGEN_DYNAMIC_RUNTIME)
     set(MSDFGEN_MSVC_RUNTIME "MultiThreaded$<$<CONFIG:Debug>:Debug>DLL")
@@ -87,11 +87,30 @@ if(MAX_WARNING_LEVEL)
     endif()
 endif()
 
+file(GLOB_RECURSE MSDFGEN_C_HEADERS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "c/*.h")
+file(GLOB_RECURSE MSDFGEN_C_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "c/*.c")
 file(GLOB_RECURSE MSDFGEN_CORE_HEADERS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "core/*.h" "core/*.hpp")
 file(GLOB_RECURSE MSDFGEN_CORE_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "core/*.cpp")
 file(GLOB_RECURSE MSDFGEN_EXT_HEADERS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "ext/*.h" "ext/*.hpp")
 file(GLOB_RECURSE MSDFGEN_EXT_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "ext/*.cpp" "lib/*.cpp")
 
+# C library
+add_library(msdfgen-c ${MSDFGEN_C_HEADERS} ${MSDFGEN_C_SOURCES})
+add_library(msdfgen::msdfgen-c ALIAS msdfgen-c)
+set_target_properties(msdfgen-c PROPERTIES PUBLIC_HEADER "${MSDFGEN_C_HEADERS}")
+set_property(TARGET msdfgen-c PROPERTY MSVC_RUNTIME_LIBRARY "${MSDFGEN_MSVC_RUNTIME}")
+target_compile_definitions(msdfgen-c PUBLIC
+    MSDFGEN_VERSION=${MSDFGEN_VERSION}
+    MSDFGEN_VERSION_MAJOR=${MSDFGEN_VERSION_MAJOR}
+    MSDFGEN_VERSION_MINOR=${MSDFGEN_VERSION_MINOR}
+    MSDFGEN_VERSION_REVISION=${MSDFGEN_VERSION_REVISION}
+    MSDFGEN_COPYRIGHT_YEAR=${MSDFGEN_COPYRIGHT_YEAR}
+)
+target_include_directories(msdfgen-c INTERFACE
+    $<INSTALL_INTERFACE:include/msdfgen>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/>
+)
+
 # Core library
 add_library(msdfgen-core "${CMAKE_CURRENT_SOURCE_DIR}/msdfgen.h" ${MSDFGEN_CORE_HEADERS} ${MSDFGEN_CORE_SOURCES})
 add_library(msdfgen::msdfgen-core ALIAS msdfgen-core)
@@ -166,8 +185,6 @@ if(NOT MSDFGEN_CORE_ONLY)
         PUBLIC
             $<INSTALL_INTERFACE:include/msdfgen>
             $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
-        PRIVATE
-            ${CMAKE_CURRENT_SOURCE_DIR}/include
     )
     set_property(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT msdfgen-ext)
 

c/CompiledShape.c

@@ -0,0 +1,79 @@
+
+#include "CompiledShape.h"
+
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+static MSDFGEN_Vector2 computeCornerVec(MSDFGEN_Vector2 normalizedDirection, MSDFGEN_Vector2 adjacentEdgeDirection1) {
+    return MSDFGEN_Vector2_normalizeOrZero(MSDFGEN_Vector2_sum(normalizedDirection, MSDFGEN_Vector2_normalizeOrZero(adjacentEdgeDirection1)));
+}
+
+void MSDFGEN_compileLinearEdge(MSDFGEN_CompiledLinearEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0) {
+    MSDFGEN_Vector2 p0p1 = MSDFGEN_Vector2_difference(point1, point0);
+    MSDFGEN_real p0p1SqLen = MSDFGEN_Vector2_squaredLength(p0p1);
+    dstEdgePtr->colorMask = colorMask;
+    dstEdgePtr->endpoint0 = point0;
+    dstEdgePtr->endpoint1 = point1;
+    dstEdgePtr->direction = MSDFGEN_Vector2_normalizeOrZero(p0p1);
+    dstEdgePtr->cornerVec0 = computeCornerVec(dstEdgePtr->direction, prevEdgeDirection1);
+    dstEdgePtr->cornerVec1 = computeCornerVec(dstEdgePtr->direction, nextEdgeDirection0);
+    dstEdgePtr->invDerivative.x = 0;
+    dstEdgePtr->invDerivative.y = 0;
+    if (p0p1SqLen != (MSDFGEN_real) 0) {
+        dstEdgePtr->invDerivative.x = -p0p1.x/p0p1SqLen;
+        dstEdgePtr->invDerivative.y = -p0p1.y/p0p1SqLen;
+    }
+}
+
+void MSDFGEN_compileQuadraticEdge(MSDFGEN_CompiledQuadraticEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 point2, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0) {
+    MSDFGEN_Vector2 p1p2 = MSDFGEN_Vector2_difference(point2, point1);
+    dstEdgePtr->colorMask = colorMask;
+    dstEdgePtr->endpoint0 = point0;
+    dstEdgePtr->endpoint1 = point2;
+    dstEdgePtr->derivative0 = MSDFGEN_Vector2_difference(point1, point0);
+    dstEdgePtr->derivative1 = MSDFGEN_Vector2_difference(p1p2, dstEdgePtr->derivative0);
+    dstEdgePtr->direction0 = dstEdgePtr->derivative0;
+    dstEdgePtr->direction1 = p1p2;
+    if (dstEdgePtr->direction0.x == (MSDFGEN_real) 0 && dstEdgePtr->direction0.y == (MSDFGEN_real) 0)
+        dstEdgePtr->direction0 = MSDFGEN_Vector2_difference(point2, point0);
+    if (dstEdgePtr->direction1.x == (MSDFGEN_real) 0 && dstEdgePtr->direction1.y == (MSDFGEN_real) 0)
+        dstEdgePtr->direction1 = MSDFGEN_Vector2_difference(point2, point0);
+    dstEdgePtr->direction0 = MSDFGEN_Vector2_normalizeOrZero(dstEdgePtr->direction0);
+    dstEdgePtr->direction1 = MSDFGEN_Vector2_normalizeOrZero(dstEdgePtr->direction1);
+    dstEdgePtr->cornerVec0 = computeCornerVec(dstEdgePtr->direction0, prevEdgeDirection1);
+    dstEdgePtr->cornerVec1 = computeCornerVec(dstEdgePtr->direction1, nextEdgeDirection0);
+}
+
+void MSDFGEN_compileCubicEdge(MSDFGEN_CompiledCubicEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 point2, MSDFGEN_Vector2 point3, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0) {
+    MSDFGEN_Vector2 p1p2 = MSDFGEN_Vector2_difference(point2, point1);
+    MSDFGEN_Vector2 p2p3 = MSDFGEN_Vector2_difference(point3, point2);
+    dstEdgePtr->colorMask = colorMask;
+    dstEdgePtr->endpoint0 = point0;
+    dstEdgePtr->endpoint1 = point3;
+    dstEdgePtr->derivative0 = MSDFGEN_Vector2_difference(point1, point0);
+    dstEdgePtr->derivative1 = MSDFGEN_Vector2_difference(p1p2, dstEdgePtr->derivative0);
+    dstEdgePtr->derivative2 = MSDFGEN_Vector2_difference(MSDFGEN_Vector2_difference(p2p3, p1p2), dstEdgePtr->derivative1);
+    dstEdgePtr->direction0 = dstEdgePtr->derivative0;
+    if (dstEdgePtr->direction0.x == (MSDFGEN_real) 0 && dstEdgePtr->direction0.y == (MSDFGEN_real) 0) {
+        dstEdgePtr->direction0 = MSDFGEN_Vector2_difference(point2, point0);
+        if (dstEdgePtr->direction0.x == (MSDFGEN_real) 0 && dstEdgePtr->direction0.y == (MSDFGEN_real) 0)
+            dstEdgePtr->direction0 = MSDFGEN_Vector2_difference(point3, point0);
+    }
+    dstEdgePtr->direction1 = p2p3;
+    if (dstEdgePtr->direction1.x == (MSDFGEN_real) 0 && dstEdgePtr->direction1.y == (MSDFGEN_real) 0) {
+        dstEdgePtr->direction1 = MSDFGEN_Vector2_difference(point3, point1);
+        if (dstEdgePtr->direction1.x == (MSDFGEN_real) 0 && dstEdgePtr->direction1.y == (MSDFGEN_real) 0)
+            dstEdgePtr->direction1 = MSDFGEN_Vector2_difference(point3, point0);
+    }
+    dstEdgePtr->direction0 = MSDFGEN_Vector2_normalizeOrZero(dstEdgePtr->direction0);
+    dstEdgePtr->direction1 = MSDFGEN_Vector2_normalizeOrZero(dstEdgePtr->direction1);
+    dstEdgePtr->cornerVec0 = computeCornerVec(dstEdgePtr->direction0, prevEdgeDirection1);
+    dstEdgePtr->cornerVec1 = computeCornerVec(dstEdgePtr->direction1, nextEdgeDirection0);
+}
+
+#ifdef __cplusplus
+}
+#endif

c/CompiledShape.h

@@ -0,0 +1,61 @@
+
+#pragma once
+
+#include <stddef.h>
+#include "Vector2.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*typedef struct {
+    void *userPtr;
+    void *(*reallocPtr)(void *userPtr, void *memoryPtr, size_t memorySize);
+    void *(*freePtr)(void *userPtr, void *memoryPtr);
+} MSDFGEN_Allocator;*/
+
+typedef struct {
+    int colorMask;
+    MSDFGEN_Vector2 endpoint0, endpoint1;
+    MSDFGEN_Vector2 cornerVec0, cornerVec1;
+    MSDFGEN_Vector2 direction;
+    MSDFGEN_Vector2 invDerivative;
+} MSDFGEN_CompiledLinearEdge;
+
+typedef struct {
+    int colorMask;
+    MSDFGEN_Vector2 endpoint0, endpoint1;
+    MSDFGEN_Vector2 cornerVec0, cornerVec1;
+    MSDFGEN_Vector2 direction0, direction1;
+    MSDFGEN_Vector2 derivative0, derivative1;
+} MSDFGEN_CompiledQuadraticEdge;
+
+typedef struct {
+    int colorMask;
+    MSDFGEN_Vector2 endpoint0, endpoint1;
+    MSDFGEN_Vector2 cornerVec0, cornerVec1;
+    MSDFGEN_Vector2 direction0, direction1;
+    MSDFGEN_Vector2 derivative0, derivative1, derivative2;
+} MSDFGEN_CompiledCubicEdge;
+
+typedef struct {
+    MSDFGEN_CompiledLinearEdge *linearEdges;
+    MSDFGEN_CompiledQuadraticEdge *quadraticEdges;
+    MSDFGEN_CompiledCubicEdge *cubicEdges;
+    size_t nLinearEdges, nQuadraticEdges, nCubicEdges;
+} MSDFGEN_CompiledShape;
+
+/*void MSDFGEN_Shape_initialize(MSDFGEN_Shape *shapePtr);
+void MSDFGEN_Shape_clear(MSDFGEN_Shape *shapePtr, MSDFGEN_Allocator allocator);*/
+
+void MSDFGEN_compileLinearEdge(MSDFGEN_CompiledLinearEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0);
+void MSDFGEN_compileQuadraticEdge(MSDFGEN_CompiledQuadraticEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 point2, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0);
+void MSDFGEN_compileCubicEdge(MSDFGEN_CompiledCubicEdge *dstEdgePtr, int colorMask, MSDFGEN_Vector2 point0, MSDFGEN_Vector2 point1, MSDFGEN_Vector2 point2, MSDFGEN_Vector2 point3, MSDFGEN_Vector2 prevEdgeDirection1, MSDFGEN_Vector2 nextEdgeDirection0);
+
+/*void MSDFGEN_Shape_addLinearEdge(MSDFGEN_Shape *shapePtr, const MSDFGEN_Shape::LinearEdge *edgePtr, MSDFGEN_Allocator allocator);
+void MSDFGEN_Shape_addQuadraticEdge(MSDFGEN_Shape *shapePtr, const MSDFGEN_Shape::QuadraticEdge *edgePtr, MSDFGEN_Allocator allocator);
+void MSDFGEN_Shape_addCubicEdge(MSDFGEN_Shape *shapePtr, const MSDFGEN_Shape::CubicEdge *edgePtr, MSDFGEN_Allocator allocator);*/
+
+#ifdef __cplusplus
+}
+#endif

c/Vector2.h

@@ -0,0 +1,59 @@
+
+#pragma once
+
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef double MSDFGEN_real;
+
+typedef struct {
+    MSDFGEN_real x, y;
+} MSDFGEN_Vector2;
+
+inline MSDFGEN_real MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2 v) {
+    return v.x*v.x+v.y*v.y;
+}
+
+inline MSDFGEN_real MSDFGEN_Vector2_dot(MSDFGEN_Vector2 a, MSDFGEN_Vector2 b) {
+    return a.x*b.x+a.y*b.y;
+}
+
+inline MSDFGEN_real MSDFGEN_Vector2_cross(MSDFGEN_Vector2 a, MSDFGEN_Vector2 b) {
+    return a.x*b.y-a.y*b.x;
+}
+
+inline MSDFGEN_Vector2 MSDFGEN_Vector2_scale(MSDFGEN_real s, MSDFGEN_Vector2 v) {
+    v.x *= s;
+    v.y *= s;
+    return v;
+}
+
+inline MSDFGEN_Vector2 MSDFGEN_Vector2_sum(MSDFGEN_Vector2 a, MSDFGEN_Vector2 b) {
+    MSDFGEN_Vector2 result;
+    result.x = a.x+b.x;
+    result.y = a.y+b.y;
+    return result;
+}
+
+inline MSDFGEN_Vector2 MSDFGEN_Vector2_difference(MSDFGEN_Vector2 a, MSDFGEN_Vector2 b) {
+    MSDFGEN_Vector2 result;
+    result.x = a.x-b.x;
+    result.y = a.y-b.y;
+    return result;
+}
+
+inline MSDFGEN_Vector2 MSDFGEN_Vector2_normalizeOrZero(MSDFGEN_Vector2 v) {
+    MSDFGEN_real length = sqrt(MSDFGEN_Vector2_squaredLength(v));
+    if (length != (MSDFGEN_real) 0) {
+        v.x /= length;
+        v.y /= length;
+    }
+    return v;
+}
+
+#ifdef __cplusplus
+}
+#endif

c/equation-solver.c

@@ -0,0 +1,76 @@
+
+#include "equation-solver.h"
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int MSDFGEN_solveQuadratic(double x[2], double a, double b, double c) {
+    // a == 0 -> linear equation
+    if (a == 0 || fabs(b) > 1e12*fabs(a)) {
+        // a == 0, b == 0 -> no solution
+        if (b == 0) {
+            if (c == 0)
+                return -1; // 0 == 0
+            return 0;
+        }
+        x[0] = -c/b;
+        return 1;
+    }
+    double dscr = b*b-4*a*c;
+    if (dscr > 0) {
+        dscr = sqrt(dscr);
+        x[0] = (-b+dscr)/(2*a);
+        x[1] = (-b-dscr)/(2*a);
+        return 2;
+    } else if (dscr == 0) {
+        x[0] = -b/(2*a);
+        return 1;
+    } else
+        return 0;
+}
+
+static int solveCubicNormed(double x[3], double a, double b, double c) {
+    double a2 = a*a;
+    double q = 1/9.*(a2-3*b);
+    double r = 1/54.*(a*(2*a2-9*b)+27*c);
+    double r2 = r*r;
+    double q3 = q*q*q;
+    a *= 1/3.;
+    if (r2 < q3) {
+        double t = r/sqrt(q3);
+        if (t < -1) t = -1;
+        if (t > 1) t = 1;
+        t = acos(t);
+        q = -2*sqrt(q);
+        x[0] = q*cos(1/3.*t)-a;
+        x[1] = q*cos(1/3.*(t+2*M_PI))-a;
+        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 v = u == 0 ? 0 : q/u;
+        x[0] = (u+v)-a;
+        if (u == v || fabs(u-v) < 1e-12*fabs(u+v)) {
+            x[1] = -.5*(u+v)-a;
+            return 2;
+        }
+        return 1;
+    }
+}
+
+int MSDFGEN_solveCubic(double x[3], double a, double b, double c, double d) {
+    if (a != 0) {
+        double bn = b/a;
+        if (fabs(bn) < 1e6) // Above this ratio, the numerical error gets larger than if we treated a as zero
+            return solveCubicNormed(x, bn, c/a, d/a);
+    }
+    return MSDFGEN_solveQuadratic(x, b, c, d);
+}
+
+#ifdef __cplusplus
+}
+#endif

c/equation-solver.h

@@ -0,0 +1,16 @@
+
+#pragma once
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// ax^2 + bx + c = 0
+int MSDFGEN_solveQuadratic(double x[2], double a, double b, double c);
+
+// ax^3 + bx^2 + cx + d = 0
+int MSDFGEN_solveCubic(double x[3], double a, double b, double c, double d);
+
+#ifdef __cplusplus
+}
+#endif

c/msdfgen-c.h

@@ -0,0 +1,85 @@
+
+#pragma once
+
+#include <stddef.h>
+#include "Vector2.h"
+#include "CompiledShape.h"
+
+#define MSDFGEN_CUBIC_SEARCH_STARTS 4
+#define MSDFGEN_CUBIC_SEARCH_STEPS 4
+
+#define MSDFGEN_ABTEST_ALT_CACHE 0
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+#if MSDFGEN_ABTEST_ALT_CACHE
+    MSDFGEN_Vector2 origin;
+#endif
+    MSDFGEN_real edgeDistance;
+} MSDFGEN_TrueDistanceEdgeCache;
+
+typedef struct {
+    MSDFGEN_Vector2 origin;
+    MSDFGEN_real minDistance;
+    MSDFGEN_TrueDistanceEdgeCache edgeData[1];
+} MSDFGEN_TrueDistanceCache;
+
+typedef struct {
+    MSDFGEN_Vector2 origin;
+    struct EdgeData {
+        MSDFGEN_real absDistance;
+        MSDFGEN_real domainDistance0, domainDistance1;
+        MSDFGEN_real perpendicularDistance0, perpendicularDistance1;
+    } edgeData[1];
+} MSDFGEN_PerpendicularDistanceCache;
+
+typedef struct {
+    struct {
+        MSDFGEN_real scale, translate;
+    } distanceMapping;
+    struct {
+        MSDFGEN_Vector2 scale, translate;
+    } projection;
+} MSDFGEN_SDFTransformation;
+
+typedef struct {
+    enum Mode {
+        DISABLED,
+        INDISCRIMINATE,
+        EDGE_PRIORITY,
+        EDGE_ONLY
+    } mode;
+    enum DistanceCheckMode {
+        DO_NOT_CHECK_DISTANCE,
+        CHECK_DISTANCE_AT_EDGE,
+        ALWAYS_CHECK_DISTANCE
+    } distanceCheckMode;
+    MSDFGEN_real minDeviationRatio;
+    MSDFGEN_real minImproveRatio;
+    void *buffer;
+} MSDFGEN_ErrorCorrectionSettings;
+
+size_t MSDFGEN_TrueDistanceCache_size(const MSDFGEN_CompiledShape *shapePtr);
+size_t MSDFGEN_PerpendicularDistanceCache_size(const MSDFGEN_CompiledShape *shapePtr);
+
+void MSDFGEN_TrueDistanceCache_initialize(MSDFGEN_TrueDistanceCache *cachePtr, const MSDFGEN_CompiledShape *shapePtr);
+void MSDFGEN_PerpendicularDistanceCache_initialize(MSDFGEN_PerpendicularDistanceCache *cachePtr, const MSDFGEN_CompiledShape *shapePtr);
+
+MSDFGEN_real MSDFGEN_signedDistance(const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_Vector2 origin, MSDFGEN_TrueDistanceCache *cachePtr);
+MSDFGEN_real MSDFGEN_perpendicularSignedDistance(const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_Vector2 origin, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+MSDFGEN_real MSDFGEN_multiSignedDistance(MSDFGEN_real dstMultiDistance[3], const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_Vector2 origin, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+
+void MSDFGEN_generateSDF(float *dstPixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_TrueDistanceCache *cachePtr);
+void MSDFGEN_generatePSDF(float *dstPixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+void MSDFGEN_generateMSDF(float *dstPixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+void MSDFGEN_generateMTSDF(float *dstPixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+
+void MSDFGEN_errorCorrectionMSDF(MSDFGEN_ErrorCorrectionSettings settings, float *pixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+void MSDFGEN_errorCorrectionMTSDF(MSDFGEN_ErrorCorrectionSettings settings, float *pixels, int width, int height, size_t rowStride, const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_SDFTransformation transformation, MSDFGEN_PerpendicularDistanceCache *cachePtr);
+
+#ifdef __cplusplus
+}
+#endif

c/msdfgen.c

@@ -0,0 +1,272 @@
+
+#include "msdfgen-c.h"
+
+#include <math.h>
+#include <float.h>
+#include "equation-solver.h"
+
+#define DISTANCE_DELTA_FACTOR 1.001
+
+//#define MSDFGEN_IF_CACHE_AND(...) // disable cache
+#define MSDFGEN_IF_CACHE_AND if
+
+extern long long MSDFGEN_PERFSTATS_CACHE_TESTS;
+extern long long MSDFGEN_PERFSTATS_CACHE_MISSES;
+
+long long MSDFGEN_PERFSTATS_CACHE_TESTS = 0;
+long long MSDFGEN_PERFSTATS_CACHE_MISSES = 0;
+
+#define MSDFGEN_PERFSTATS_CACHE_TEST() (++MSDFGEN_PERFSTATS_CACHE_TESTS)
+#define MSDFGEN_PERFSTATS_CACHE_MISS() (++MSDFGEN_PERFSTATS_CACHE_MISSES)
+
+//#define NO_CACHE_AND_SQUARE_DISTANCE
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+static int nonZeroSign(MSDFGEN_real x) {
+    return x > (MSDFGEN_real) 0 ? 1 : -1;
+}
+
+static int crossNonZeroSign(MSDFGEN_Vector2 a, MSDFGEN_Vector2 b) {
+    return a.x*b.y > b.x*a.y ? 1 : -1;
+}
+
+size_t MSDFGEN_TrueDistanceCache_size(const MSDFGEN_CompiledShape *shapePtr) {
+    return sizeof(MSDFGEN_TrueDistanceCache) + (shapePtr->nLinearEdges+shapePtr->nQuadraticEdges+shapePtr->nCubicEdges)*sizeof(MSDFGEN_TrueDistanceEdgeCache);
+}
+
+void MSDFGEN_TrueDistanceCache_initialize(MSDFGEN_TrueDistanceCache *cachePtr, const MSDFGEN_CompiledShape *shapePtr) {
+    MSDFGEN_TrueDistanceEdgeCache *cur, *end;
+    cachePtr->origin.x = 0;
+    cachePtr->origin.y = 0;
+    cachePtr->minDistance = .0625f*FLT_MAX;
+    for (cur = cachePtr->edgeData, end = cachePtr->edgeData+(shapePtr->nLinearEdges+shapePtr->nQuadraticEdges+shapePtr->nCubicEdges); cur < end; ++cur) {
+        #if MSDFGEN_ABTEST_ALT_CACHE
+            cur->origin.x = 0;
+            cur->origin.y = 0;
+        #else
+            cur->edgeDistance = 0;
+        #endif
+    }
+}
+
+#ifdef NO_CACHE_AND_SQUARE_DISTANCE
+
+// INCOMPLETE !!!!
+
+MSDFGEN_real MSDFGEN_signedDistance(const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_Vector2 origin, MSDFGEN_TrueDistanceCache *cachePtr) {
+    int sign = -1;
+    MSDFGEN_real sqd = FLT_MAX;
+    size_t i;
+
+    for (i = 0; i < shapePtr->nLinearEdges; ++i) {
+        MSDFGEN_Vector2 originP0 = MSDFGEN_Vector2_difference(shapePtr->linearEdges[i].endpoint0, origin);
+        MSDFGEN_real originP0SqD = MSDFGEN_Vector2_squaredLength(originP0);
+        MSDFGEN_real t = MSDFGEN_Vector2_dot(originP0, shapePtr->linearEdges[i].invDerivative);
+        if (originP0SqD < sqd) {
+            sqd = originP0SqD;
+            sign = crossNonZeroSign(shapePtr->linearEdges[i].cornerVec0, MSDFGEN_Vector2_sum(shapePtr->linearEdges[i].direction, originP0));
+        }
+        if (t > (MSDFGEN_real) 0 && t < (MSDFGEN_real) 1) {
+            MSDFGEN_real pd = MSDFGEN_Vector2_cross(shapePtr->linearEdges[i].direction, originP0);
+            MSDFGEN_real sqpd = pd*pd;
+            if (sqpd < sqd) {
+                sqd = sqpd;
+                sign = nonZeroSign(pd);
+            }
+        }
+    }
+
+    for (i = 0; i < shapePtr->nQuadraticEdges; ++i) {
+        MSDFGEN_Vector2 originP0 = MSDFGEN_Vector2_difference(shapePtr->quadraticEdges[i].endpoint0, origin);
+        MSDFGEN_real originP0SqD = MSDFGEN_Vector2_squaredLength(originP0);
+        MSDFGEN_real a = MSDFGEN_Vector2_squaredLength(shapePtr->quadraticEdges[i].derivative1);
+        MSDFGEN_real b = (MSDFGEN_real) 3*MSDFGEN_Vector2_dot(shapePtr->quadraticEdges[i].derivative0, shapePtr->quadraticEdges[i].derivative1);
+        MSDFGEN_real c = (MSDFGEN_real) 2*MSDFGEN_Vector2_squaredLength(shapePtr->quadraticEdges[i].derivative0) + MSDFGEN_Vector2_dot(originP0, shapePtr->quadraticEdges[i].derivative1);
+        MSDFGEN_real d = MSDFGEN_Vector2_dot(originP0, shapePtr->quadraticEdges[i].derivative0);
+        double t[3];
+        int solutions = MSDFGEN_solveCubic(t, a, b, c, d);
+        if (originP0SqD < sqd) {
+            sqd = originP0SqD;
+            sign = crossNonZeroSign(shapePtr->quadraticEdges[i].cornerVec0, MSDFGEN_Vector2_sum(shapePtr->quadraticEdges[i].direction0, originP0));
+        }
+        #define MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t) if (t > 0 && t < 1) { \
+            MSDFGEN_Vector2 originP = MSDFGEN_Vector2_sum(MSDFGEN_Vector2_sum(originP0, MSDFGEN_Vector2_scale((MSDFGEN_real) (2*t), shapePtr->quadraticEdges[i].derivative0)), MSDFGEN_Vector2_scale((MSDFGEN_real) (t*t), shapePtr->quadraticEdges[i].derivative1)); \
+            MSDFGEN_real originPSqD = MSDFGEN_Vector2_squaredLength(originP); \
+            if (originPSqD < sqd) { \
+                MSDFGEN_Vector2 direction = MSDFGEN_Vector2_sum(shapePtr->quadraticEdges[i].derivative0, MSDFGEN_Vector2_scale(t, shapePtr->quadraticEdges[i].derivative1)); \
+                sqd = originPSqD; \
+                sign = crossNonZeroSign(direction, originP); \
+            } \
+        }
+        if (solutions > 0) {
+            MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[0]);
+            if (solutions > 1) {
+                MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[1]);
+                if (solutions > 2)
+                    MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[2]);
+            }
+        }
+    }
+
+    return (MSDFGEN_real) sign*sqrt(sqd);
+}
+
+#else
+
+#if MSDFGEN_ABTEST_ALT_CACHE
+#define MSDFGEN_IF_TRUE_DISTANCE_UNCACHED MSDFGEN_PERFSTATS_CACHE_TEST(); MSDFGEN_IF_CACHE_AND(++edgeCache, edgeCache[-1].edgeDistance-DISTANCE_DELTA_FACTOR*sqrt(MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2_difference(origin, edgeCache[-1].origin))) <= minDistance)
+#else
+#define MSDFGEN_IF_TRUE_DISTANCE_UNCACHED MSDFGEN_PERFSTATS_CACHE_TEST(); MSDFGEN_IF_CACHE_AND((edgeCache++->edgeDistance -= cacheDelta) <= minDistance)
+#endif
+
+MSDFGEN_real MSDFGEN_signedDistance(const MSDFGEN_CompiledShape *shapePtr, MSDFGEN_Vector2 origin, MSDFGEN_TrueDistanceCache *cachePtr) {
+
+    MSDFGEN_real cacheDelta = DISTANCE_DELTA_FACTOR*sqrt(MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2_difference(origin, cachePtr->origin)));
+    MSDFGEN_real minDistance = cachePtr->minDistance+cacheDelta;
+    int distanceSign = -1;
+    MSDFGEN_TrueDistanceEdgeCache *edgeCache = cachePtr->edgeData;
+    size_t i;
+
+    for (i = 0; i < shapePtr->nLinearEdges; ++i) {
+        MSDFGEN_IF_TRUE_DISTANCE_UNCACHED {
+            MSDFGEN_Vector2 originP0 = MSDFGEN_Vector2_difference(shapePtr->linearEdges[i].endpoint0, origin);
+            MSDFGEN_real originP0Dist = sqrt(MSDFGEN_Vector2_squaredLength(originP0));
+            MSDFGEN_real t = MSDFGEN_Vector2_dot(originP0, shapePtr->linearEdges[i].invDerivative);
+            if (originP0Dist < minDistance) {
+                minDistance = originP0Dist;
+                distanceSign = crossNonZeroSign(shapePtr->linearEdges[i].cornerVec0, MSDFGEN_Vector2_sum(shapePtr->linearEdges[i].direction, originP0));
+            }
+            edgeCache[-1].edgeDistance = originP0Dist;
+            if (t > (MSDFGEN_real) 0 && t < (MSDFGEN_real) 1) {
+                MSDFGEN_real psd = MSDFGEN_Vector2_cross(shapePtr->linearEdges[i].direction, originP0);
+                MSDFGEN_real pd = fabs(psd);
+                if (pd < minDistance) {
+                    minDistance = pd;
+                    distanceSign = nonZeroSign(psd);
+                }
+                edgeCache[-1].edgeDistance = pd;
+            } else {
+                MSDFGEN_real originP1Dist = sqrt(MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2_difference(shapePtr->linearEdges[i].endpoint1, origin)));
+                if (originP1Dist < edgeCache[-1].edgeDistance)
+                    edgeCache[-1].edgeDistance = originP1Dist;
+            }
+            MSDFGEN_PERFSTATS_CACHE_MISS();
+            #if MSDFGEN_ABTEST_ALT_CACHE
+                edgeCache[-1].origin = origin;
+            #endif
+        }
+    }
+
+    for (i = 0; i < shapePtr->nQuadraticEdges; ++i) {
+        MSDFGEN_IF_TRUE_DISTANCE_UNCACHED {
+            MSDFGEN_Vector2 originP0 = MSDFGEN_Vector2_difference(shapePtr->quadraticEdges[i].endpoint0, origin);
+            MSDFGEN_real originP0Dist = sqrt(MSDFGEN_Vector2_squaredLength(originP0));
+            MSDFGEN_real originP1Dist = sqrt(MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2_difference(shapePtr->quadraticEdges[i].endpoint1, origin)));
+            MSDFGEN_real a = MSDFGEN_Vector2_squaredLength(shapePtr->quadraticEdges[i].derivative1);
+            MSDFGEN_real b = (MSDFGEN_real) 3*MSDFGEN_Vector2_dot(shapePtr->quadraticEdges[i].derivative0, shapePtr->quadraticEdges[i].derivative1);
+            MSDFGEN_real c = (MSDFGEN_real) 2*MSDFGEN_Vector2_squaredLength(shapePtr->quadraticEdges[i].derivative0) + MSDFGEN_Vector2_dot(originP0, shapePtr->quadraticEdges[i].derivative1);
+            MSDFGEN_real d = MSDFGEN_Vector2_dot(originP0, shapePtr->quadraticEdges[i].derivative0);
+            double t[3];
+            int solutions = MSDFGEN_solveCubic(t, a, b, c, d);
+            if (originP0Dist < minDistance) {
+                minDistance = originP0Dist;
+                distanceSign = crossNonZeroSign(shapePtr->quadraticEdges[i].cornerVec0, MSDFGEN_Vector2_sum(shapePtr->quadraticEdges[i].direction0, originP0));
+            }
+            edgeCache[-1].edgeDistance = originP0Dist;
+            if (originP1Dist < edgeCache[-1].edgeDistance)
+                edgeCache[-1].edgeDistance = originP1Dist;
+            #define MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t) if (t > 0 && t < 1) { \
+                MSDFGEN_Vector2 originP = MSDFGEN_Vector2_sum(MSDFGEN_Vector2_sum(originP0, MSDFGEN_Vector2_scale((MSDFGEN_real) (2*t), shapePtr->quadraticEdges[i].derivative0)), MSDFGEN_Vector2_scale((MSDFGEN_real) (t*t), shapePtr->quadraticEdges[i].derivative1)); \
+                MSDFGEN_real originPDist = sqrt(MSDFGEN_Vector2_squaredLength(originP)); \
+                if (originPDist < minDistance) { \
+                    MSDFGEN_Vector2 direction = MSDFGEN_Vector2_sum(shapePtr->quadraticEdges[i].derivative0, MSDFGEN_Vector2_scale(t, shapePtr->quadraticEdges[i].derivative1)); \
+                    minDistance = originPDist; \
+                    distanceSign = crossNonZeroSign(direction, originP); \
+                } \
+                if (originPDist < edgeCache[-1].edgeDistance) \
+                    edgeCache[-1].edgeDistance = originPDist; \
+            }
+            if (solutions > 0) {
+                MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[0]);
+                if (solutions > 1) {
+                    MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[1]);
+                    if (solutions > 2)
+                        MSDFGEN_SD_RESOLVE_QUADRATIC_SOLUTION(t[2]);
+                }
+            }
+            MSDFGEN_PERFSTATS_CACHE_MISS();
+            #if MSDFGEN_ABTEST_ALT_CACHE
+                edgeCache[-1].origin = origin;
+            #endif
+        }
+    }
+
+    for (i = 0; i < shapePtr->nCubicEdges; ++i) {
+        MSDFGEN_IF_TRUE_DISTANCE_UNCACHED {
+            int start, step;
+            MSDFGEN_Vector2 originP0 = MSDFGEN_Vector2_difference(shapePtr->cubicEdges[i].endpoint0, origin);
+            MSDFGEN_real originP0Dist = sqrt(MSDFGEN_Vector2_squaredLength(originP0));
+            MSDFGEN_real originP1Dist = sqrt(MSDFGEN_Vector2_squaredLength(MSDFGEN_Vector2_difference(shapePtr->cubicEdges[i].endpoint1, origin)));
+            if (originP0Dist < minDistance) {
+                minDistance = originP0Dist;
+                distanceSign = crossNonZeroSign(shapePtr->cubicEdges[i].cornerVec0, MSDFGEN_Vector2_sum(shapePtr->cubicEdges[i].direction0, originP0));
+            }
+            edgeCache[-1].edgeDistance = originP0Dist;
+            if (originP1Dist < edgeCache[-1].edgeDistance)
+                edgeCache[-1].edgeDistance = originP1Dist;
+            for (start = 0; start <= MSDFGEN_CUBIC_SEARCH_STARTS; ++start) {
+                int remainingSteps = MSDFGEN_CUBIC_SEARCH_STEPS;
+                MSDFGEN_real t, improvedT = (MSDFGEN_real) 1/(MSDFGEN_real) MSDFGEN_CUBIC_SEARCH_STARTS*(MSDFGEN_real) start;
+                MSDFGEN_Vector2 originP, derivative0, derivative1;
+                do {
+                    t = improvedT;
+                    originP = MSDFGEN_Vector2_sum(
+                        MSDFGEN_Vector2_sum(
+                            MSDFGEN_Vector2_sum(originP0, MSDFGEN_Vector2_scale((MSDFGEN_real) 3*t, shapePtr->cubicEdges[i].derivative0)),
+                            MSDFGEN_Vector2_scale((MSDFGEN_real) 3*(t*t), shapePtr->cubicEdges[i].derivative1)
+                        ), MSDFGEN_Vector2_scale(t*t*t, shapePtr->cubicEdges[i].derivative2)
+                    );
+                    derivative0 = MSDFGEN_Vector2_sum(
+                        MSDFGEN_Vector2_sum(
+                            MSDFGEN_Vector2_scale((MSDFGEN_real) 3, shapePtr->cubicEdges[i].derivative0),
+                            MSDFGEN_Vector2_scale((MSDFGEN_real) 6*t, shapePtr->cubicEdges[i].derivative1)
+                        ), MSDFGEN_Vector2_scale((MSDFGEN_real) 3*(t*t), shapePtr->cubicEdges[i].derivative2)
+                    );
+                    if (!remainingSteps--)
+                        break;
+                    derivative1 = MSDFGEN_Vector2_sum(
+                        MSDFGEN_Vector2_scale((MSDFGEN_real) 6, shapePtr->cubicEdges[i].derivative1),
+                        MSDFGEN_Vector2_scale((MSDFGEN_real) 6*t, shapePtr->cubicEdges[i].derivative2)
+                    );
+                    improvedT = t-MSDFGEN_Vector2_dot(originP, derivative0)/(MSDFGEN_Vector2_squaredLength(derivative0)+MSDFGEN_Vector2_dot(originP, derivative1));
+                } while (improvedT > (MSDFGEN_real) 0 && improvedT < (MSDFGEN_real) 1);
+                if (t > (MSDFGEN_real) 0 && t < (MSDFGEN_real) 1) {
+                    MSDFGEN_real originPDist = sqrt(MSDFGEN_Vector2_squaredLength(originP));
+                    if (originPDist < minDistance) {
+                        minDistance = originPDist;
+                        distanceSign = crossNonZeroSign(derivative0, originP);
+                    }
+                    if (originPDist < edgeCache[-1].edgeDistance)
+                        edgeCache[-1].edgeDistance = originPDist;
+                }
+            }
+            MSDFGEN_PERFSTATS_CACHE_MISS();
+            #if MSDFGEN_ABTEST_ALT_CACHE
+                edgeCache[-1].origin = origin;
+            #endif
+        }
+    }
+
+    cachePtr->origin = origin;
+    cachePtr->minDistance = minDistance;
+    return (MSDFGEN_real) distanceSign*minDistance;
+}
+
+#endif
+
+#ifdef __cplusplus
+}
+#endif

core/Shape.cpp

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

core/convergent-curve-ordering.cpp

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

core/convergent-curve-ordering.h

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

core/edge-segments.cpp

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