diff --git a/core/approximate-sdf.cpp b/core/approximate-sdf.cpp
new file mode 100644
index 0000000..4a6f67f
--- /dev/null
+++ b/core/approximate-sdf.cpp
@@ -0,0 +1,221 @@
+
+#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 Projection &projection, double outerRange, double innerRange) {
+    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 = projection.unprojectVector(Vector2(1));
+    float dLimit = float(max(outerRange, innerRange));
+    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 = projection.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 = projection.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 = projection.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 = projection.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 = projection.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 = projection.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 = projection.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 = projection.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);
+            }
+        }
+    }
+
+    float rangeFactor = 1.f/float(outerRange+innerRange);
+    float zeroBias = rangeFactor*float(outerRange);
+    for (float *p = output.pixels, *end = output.pixels+output.width*output.height; p < end; ++p)
+        *p = rangeFactor**p+zeroBias;
+}
+
+void approximateSDF(const BitmapRef<float, 1> &output, const Shape &shape, const Projection &projection, double range) {
+    approximateSDF(output, shape, projection, .5*range, .5*range);
+}
+
+}
diff --git a/core/approximate-sdf.h b/core/approximate-sdf.h
new file mode 100644
index 0000000..8d4cf1a
--- /dev/null
+++ b/core/approximate-sdf.h
@@ -0,0 +1,14 @@
+
+#pragma once
+
+#include "Projection.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 Projection &projection, double outerRange, double innerRange);
+void approximateSDF(const BitmapRef<float, 1> &output, const Shape &shape, const Projection &projection, double range);
+
+}
diff --git a/core/edge-segments.cpp b/core/edge-segments.cpp
index 400c426..b72f0f9 100644
--- a/core/edge-segments.cpp
+++ b/core/edge-segments.cpp
@@ -338,6 +338,18 @@ SignedDistance CubicSegment::signedDistance(Point2 origin, double &param) const
 #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);
@@ -347,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;
@@ -401,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;
@@ -463,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;
diff --git a/core/edge-segments.h b/core/edge-segments.h
index 3880bec..950c407 100644
--- a/core/edge-segments.h
+++ b/core/edge-segments.h
@@ -37,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;
 
@@ -71,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();
@@ -100,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();
@@ -130,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();
diff --git a/msdfgen.h b/msdfgen.h
index ae909ef..a5393e6 100644
--- a/msdfgen.h
+++ b/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"