Interpolation error correction WIP

Msdfgen.vcxproj

@@ -296,6 +296,7 @@
   </ItemDefinitionGroup>
   <ItemGroup>
     <ClInclude Include="core\arithmetics.hpp" />
+    <ClInclude Include="core\bitmap-interpolation.hpp" />
     <ClInclude Include="core\Bitmap.h" />
     <ClInclude Include="core\Bitmap.hpp" />
     <ClInclude Include="core\BitmapRef.hpp" />
@@ -308,6 +309,7 @@
     <ClInclude Include="core\EdgeHolder.h" />
     <ClInclude Include="core\equation-solver.h" />
     <ClInclude Include="core\estimate-sdf-error.h" />
+    <ClInclude Include="core\interpolation-error-correction.h" />
     <ClInclude Include="core\pixel-conversion.hpp" />
     <ClInclude Include="core\rasterization.h" />
     <ClInclude Include="core\render-sdf.h" />
@@ -334,6 +336,7 @@
     <ClCompile Include="core\EdgeHolder.cpp" />
     <ClCompile Include="core\equation-solver.cpp" />
     <ClCompile Include="core\estimate-sdf-error.cpp" />
+    <ClCompile Include="core\interpolation-error-correction.cpp" />
     <ClCompile Include="core\rasterization.cpp" />
     <ClCompile Include="core\render-sdf.cpp" />
     <ClCompile Include="core\save-bmp.cpp" />

Msdfgen.vcxproj.filters

@@ -105,6 +105,12 @@
     <ClInclude Include="core\estimate-sdf-error.h">
       <Filter>Core</Filter>
     </ClInclude>
+    <ClInclude Include="core\interpolation-error-correction.h">
+      <Filter>Core</Filter>
+    </ClInclude>
+    <ClInclude Include="core\bitmap-interpolation.hpp">
+      <Filter>Core</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="main.cpp">
@@ -179,6 +185,9 @@
     <ClCompile Include="core\estimate-sdf-error.cpp">
       <Filter>Core</Filter>
     </ClCompile>
+    <ClCompile Include="core\interpolation-error-correction.cpp">
+      <Filter>Core</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="Msdfgen.rc">

core/bitmap-interpolation.hpp

@@ -0,0 +1,25 @@
+
+#pragma once
+
+#include "arithmetics.hpp"
+#include "Vector2.h"
+#include "BitmapRef.hpp"
+
+namespace msdfgen {
+
+template <typename T, int N>
+static void interpolate(T *output, const BitmapConstRef<T, N> &bitmap, Point2 pos) {
+    pos -= .5;
+    int l = (int) floor(pos.x);
+    int b = (int) floor(pos.y);
+    int r = l+1;
+    int t = b+1;
+    double lr = pos.x-l;
+    double bt = pos.y-b;
+    l = clamp(l, bitmap.width-1), r = clamp(r, bitmap.width-1);
+    b = clamp(b, bitmap.height-1), t = clamp(t, bitmap.height-1);
+    for (int i = 0; i < N; ++i)
+        output[i] = mix(mix(bitmap(l, b)[i], bitmap(r, b)[i], lr), mix(bitmap(l, t)[i], bitmap(r, t)[i], lr), bt);
+}
+
+}

core/interpolation-error-correction.cpp

@@ -0,0 +1,194 @@
+
+#include "interpolation-error-correction.h"
+
+#include <vector>
+#include <utility>
+#include "arithmetics.hpp"
+#include "equation-solver.h"
+#include "bitmap-interpolation.hpp"
+
+#include <set> // TEMP
+
+#ifdef _DEBUG
+#define dPrintf printf
+#else
+#define dPrintf(...)
+#endif
+
+namespace msdfgen {
+
+// TODO EXPOSE
+static double simpleSignedDistance(const Shape &shape, const msdfgen::Point2 &p) {
+    double dummy;
+    msdfgen::SignedDistance minDistance;
+    for (const msdfgen::Contour &contour : shape.contours)
+        for (const msdfgen::EdgeHolder &edge : contour.edges) {
+            msdfgen::SignedDistance distance = edge->signedDistance(p, dummy);
+            if (distance < minDistance)
+                minDistance = distance;
+        }
+    return minDistance.distance;
+}
+
+static bool isHotspot(float am, float bm, float xm) {
+    return (am > .5f && bm > .5f && xm < .5f) || (am < .5f && bm < .5f && xm > .5f); // only at edge
+    return clamp(median(am, bm, xm)) != clamp(xm); // anywhere
+}
+
+static int findLinearChannelHotspots(double t[1], const float *a, const float *b, float dA, float dB) {
+    int found = 0;
+    double x = (double) dA/(dA-dB);
+    if (x > 0 && x < 1) {
+        float am = median(a[0], a[1], a[2]);
+        float bm = median(b[0], b[1], b[2]);
+        float xm = median(
+            mix(a[0], b[0], x),
+            mix(a[1], b[1], x),
+            mix(a[2], b[2], x)
+        );
+        if (isHotspot(am, bm, xm))
+            t[found++] = x;
+    }
+    return found;
+}
+
+static int findDiagonalChannelHotspots(double t[2], const float *a, const float *b, const float *c, const float *d, float dA, float dB, float dC, float dD) {
+    int found = 0;
+    double x[2];
+    int solutions = solveQuadratic(x, (dD-dC)-(dB-dA), dC+dB-2*dA, dA);
+    for (int i = 0; i < solutions; ++i)
+        if (x[i] > 0 && x[i] < 1) {
+            float am = median(a[0], a[1], a[2]);
+            float bm = median(b[0], b[1], b[2]);
+            float xm = median(
+                mix(mix(a[0], b[0], x[i]), mix(c[0], d[0], x[i]), x[i]),
+                mix(mix(a[1], b[1], x[i]), mix(c[1], d[1], x[i]), x[i]),
+                mix(mix(a[2], b[2], x[i]), mix(c[2], d[2], x[i]), x[i])
+            );
+            if (isHotspot(am, bm, xm))
+                t[found++] = x[i];
+        }
+    return found;
+}
+
+static int findLinearHotspots(double t[3], const float *a, const float *b) {
+    int found = 0;
+    found += findLinearChannelHotspots(t+found, a, b, a[1]-a[0], b[1]-b[0]);
+    found += findLinearChannelHotspots(t+found, a, b, a[2]-a[1], b[2]-b[1]);
+    found += findLinearChannelHotspots(t+found, a, b, a[0]-a[2], b[0]-b[2]);
+    return found;
+}
+
+static int findDiagonalHotspots(double t[6], const float *a, const float *b, const float *c, const float *d) {
+    int found = 0;
+    found += findDiagonalChannelHotspots(t+found, a, b, c, d, a[1]-a[0], b[1]-b[0], c[1]-c[0], d[1]-d[0]);
+    found += findDiagonalChannelHotspots(t+found, a, b, c, d, a[2]-a[1], b[2]-b[1], c[2]-c[1], d[2]-d[1]);
+    found += findDiagonalChannelHotspots(t+found, a, b, c, d, a[0]-a[2], b[0]-b[2], c[0]-c[2], d[0]-d[2]);
+    return found;
+}
+
+template <int N>
+void findHotspots(std::vector<Point2> &hotspots, const BitmapConstRef<float, N> &sdf) {
+    // All hotspots intersect either the horizontal, vertical, or diagonal line that connects neighboring texels
+    // Horizontal:
+    for (int y = 0; y < sdf.height; ++y) {
+        const float *left = sdf(0, y);
+        const float *right = sdf(1, y);
+        for (int x = 0; x < sdf.width-1; ++x) {
+            double t[3];
+            int found = findLinearHotspots(t, left, right);
+            for (int i = 0; i < found; ++i)
+                hotspots.push_back(Point2(x+.5+t[i], y+.5));
+            left += N, right += N;
+        }
+    }
+    // Vertical:
+    for (int y = 0; y < sdf.height-1; ++y) {
+        const float *bottom = sdf(0, y);
+        const float *top = sdf(0, y+1);
+        for (int x = 0; x < sdf.width; ++x) {
+            double t[3];
+            int found = findLinearHotspots(t, bottom, top);
+            for (int i = 0; i < found; ++i)
+                hotspots.push_back(Point2(x+.5, y+.5+t[i]));
+            bottom += N, top += N;
+        }
+    }
+    // Diagonal:
+    for (int y = 0; y < sdf.height-1; ++y) {
+        const float *lb = sdf(0, y);
+        const float *rb = sdf(1, y);
+        const float *lt = sdf(0, y+1);
+        const float *rt = sdf(1, y+1);
+        for (int x = 0; x < sdf.width-1; ++x) {
+            double t[6];
+            int found = 0;
+            found = findDiagonalHotspots(t, lb, rb, lt, rt);
+            for (int i = 0; i < found; ++i)
+                hotspots.push_back(Point2(x+.5+t[i], y+.5+t[i]));
+            found = findDiagonalHotspots(t, lt, rt, lb, rb);
+            for (int i = 0; i < found; ++i)
+                hotspots.push_back(Point2(x+.5+t[i], y+1.5-t[i]));
+            lb += N, rb += N, lt += N, rt += N;
+        }
+    }
+}
+
+template <int N>
+static void msdfInterpolationErrorCorrectionInner(const BitmapRef<float, N> &sdf, const Shape &shape, double range, const Vector2 &scale, const Vector2 &translate, bool overlapSupport) {
+    const float gray[] = { .5f, .5f, .5f, 1.f };
+    std::vector<Point2> hotspots;
+    findHotspots(hotspots, BitmapConstRef<float, N>(sdf));
+    dPrintf("Found %d hotspots (%.3fx texels)\n", (int) hotspots.size(), (double) hotspots.size()/(sdf.width*sdf.height));
+    //std::vector<std::pair<int, int> > artifacts;
+    //artifacts.reserve(hotspots.size());
+    std::set<std::pair<int, int> > artifacts;
+    for (std::vector<Point2>::const_iterator hotspot = hotspots.begin(); hotspot != hotspots.end(); ++hotspot) {
+        Point2 pos = *hotspot/scale-translate;
+        double distance = simpleSignedDistance(shape, pos);
+        float sd = float(distance/range+.5);
+
+        float *subject = sdf((int) hotspot->x, (int) hotspot->y);
+        float texel[N];
+        memcpy(texel, subject, N*sizeof(float));
+
+        float msd[N];
+        interpolate(msd, BitmapConstRef<float, N>(sdf), *hotspot);
+        float oldSsd = median(msd[0], msd[1], msd[2]);
+
+        float med = median(subject[0], subject[1], subject[2]);
+        subject[0] = med, subject[1] = med, subject[2] = med;
+
+        interpolate(msd, BitmapConstRef<float, N>(sdf), *hotspot);
+        float newSsd = median(msd[0], msd[1], msd[2]);
+
+        memcpy(subject, texel, N*sizeof(float));
+
+        dPrintf("Real sd = %f, old = %f, new = %f\n", sd, oldSsd, newSsd);
+        //memcpy(sdf((int) hotspot->x, (int) hotspot->y), gray, N*sizeof(float));
+
+        bool significant = fabsf(newSsd-sd) < fabsf(oldSsd-sd);
+        significant = (newSsd-.5f)*(oldSsd-.5f) < 0;
+        significant = fabsf(newSsd-sd) < fabsf(oldSsd-sd);
+
+        if (significant)
+            artifacts.insert(std::make_pair((int) hotspot->x, (int) hotspot->y));
+            //artifacts.push_back(std::make_pair((int) hotspot->x, (int) hotspot->y));
+    }
+    dPrintf("Found %d artifacts (%.2f%% hotspots, %.2f%% texels)\n", (int) artifacts.size(), 100.*artifacts.size()/hotspots.size(), 100.*artifacts.size()/(sdf.width*sdf.height));
+    for (std::set<std::pair<int, int> >::const_iterator artifact = artifacts.begin(); artifact != artifacts.end(); ++artifact) {
+        float *pixel = sdf(artifact->first, artifact->second);
+        float med = median(pixel[0], pixel[1], pixel[2]);
+        pixel[0] = med, pixel[1] = med, pixel[2] = med;
+    }
+}
+
+void msdfInterpolationErrorCorrection(const BitmapRef<float, 3> &sdf, const Shape &shape, double range, const Vector2 &scale, const Vector2 &translate, bool overlapSupport) {
+    msdfInterpolationErrorCorrectionInner(sdf, shape, range, scale, translate, overlapSupport);
+}
+
+void msdfInterpolationErrorCorrection(const BitmapRef<float, 4> &sdf, const Shape &shape, double range, const Vector2 &scale, const Vector2 &translate, bool overlapSupport) {
+    msdfInterpolationErrorCorrectionInner(sdf, shape, range, scale, translate, overlapSupport);
+}
+
+}

core/interpolation-error-correction.h

@@ -0,0 +1,13 @@
+
+#pragma once
+
+#include "Vector2.h"
+#include "Shape.h"
+#include "BitmapRef.hpp"
+
+namespace msdfgen {
+
+void msdfInterpolationErrorCorrection(const BitmapRef<float, 3> &sdf, const Shape &shape, double range, const Vector2 &scale, const Vector2 &translate, bool overlapSupport = true);
+void msdfInterpolationErrorCorrection(const BitmapRef<float, 4> &sdf, const Shape &shape, double range, const Vector2 &scale, const Vector2 &translate, bool overlapSupport = true);
+
+}

core/render-sdf.cpp

@@ -3,25 +3,10 @@
 
 #include "arithmetics.hpp"
 #include "pixel-conversion.hpp"
+#include "bitmap-interpolation.hpp"
 
 namespace msdfgen {
 
-template <typename T, int N>
-static void sample(T *output, const BitmapConstRef<T, N> &bitmap, Point2 pos) {
-    double x = pos.x*bitmap.width-.5;
-    double y = pos.y*bitmap.height-.5;
-    int l = (int) floor(x);
-    int b = (int) floor(y);
-    int r = l+1;
-    int t = b+1;
-    double lr = x-l;
-    double bt = y-b;
-    l = clamp(l, bitmap.width-1), r = clamp(r, bitmap.width-1);
-    b = clamp(b, bitmap.height-1), t = clamp(t, bitmap.height-1);
-    for (int i = 0; i < N; ++i)
-        output[i] = mix(mix(bitmap(l, b)[i], bitmap(r, b)[i], lr), mix(bitmap(l, t)[i], bitmap(r, t)[i], lr), bt);
-}
-
 static float distVal(float dist, double pxRange) {
     if (!pxRange)
         return (float) (dist > .5f);
@@ -29,21 +14,23 @@ static float distVal(float dist, double pxRange) {
 }
 
 void renderSDF(const BitmapRef<float, 1> &output, const BitmapConstRef<float, 1> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd;
-            sample(&sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(&sd, sdf, scale*Point2(x+.5, y+.5));
             *output(x, y) = distVal(sd, pxRange);
         }
 }
 
 void renderSDF(const BitmapRef<float, 3> &output, const BitmapConstRef<float, 1> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd;
-            sample(&sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(&sd, sdf, scale*Point2(x+.5, y+.5));
             float v = distVal(sd, pxRange);
             output(x, y)[0] = v;
             output(x, y)[1] = v;
@@ -52,21 +39,23 @@ void renderSDF(const BitmapRef<float, 3> &output, const BitmapConstRef<float, 1>
 }
 
 void renderSDF(const BitmapRef<float, 1> &output, const BitmapConstRef<float, 3> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd[3];
-            sample(sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(sd, sdf, scale*Point2(x+.5, y+.5));
             *output(x, y) = distVal(median(sd[0], sd[1], sd[2]), pxRange);
         }
 }
 
 void renderSDF(const BitmapRef<float, 3> &output, const BitmapConstRef<float, 3> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd[3];
-            sample(sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(sd, sdf, scale*Point2(x+.5, y+.5));
             output(x, y)[0] = distVal(sd[0], pxRange);
             output(x, y)[1] = distVal(sd[1], pxRange);
             output(x, y)[2] = distVal(sd[2], pxRange);
@@ -74,21 +63,23 @@ void renderSDF(const BitmapRef<float, 3> &output, const BitmapConstRef<float, 3>
 }
 
 void renderSDF(const BitmapRef<float, 1> &output, const BitmapConstRef<float, 4> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd[4];
-            sample(sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(sd, sdf, scale*Point2(x+.5, y+.5));
             *output(x, y) = distVal(median(sd[0], sd[1], sd[2]), pxRange);
         }
 }
 
 void renderSDF(const BitmapRef<float, 4> &output, const BitmapConstRef<float, 4> &sdf, double pxRange) {
+    Vector2 scale((double) sdf.width/output.width, (double) sdf.height/output.height);
     pxRange *= (double) (output.width+output.height)/(sdf.width+sdf.height);
     for (int y = 0; y < output.height; ++y)
         for (int x = 0; x < output.width; ++x) {
             float sd[4];
-            sample(sd, sdf, Point2((x+.5)/output.width, (y+.5)/output.height));
+            interpolate(sd, sdf, scale*Point2(x+.5, y+.5));
             output(x, y)[0] = distVal(sd[0], pxRange);
             output(x, y)[1] = distVal(sd[1], pxRange);
             output(x, y)[2] = distVal(sd[2], pxRange);

main.cpp

@@ -668,6 +668,8 @@ int main(int argc, const char * const *argv) {
     if (suggestHelp)
         printf("Use -help for more information.\n");
 
+    edgeThreshold = 0;
+
     // Load input
     Vector2 svgDims;
     double glyphAdvance = 0;
@@ -878,11 +880,13 @@ int main(int argc, const char * const *argv) {
                 distanceSignCorrection(msdf, shape, scale, translate, fillRule);
                 if (edgeThreshold > 0)
                     msdfErrorCorrection(msdf, edgeThreshold/(scale*range));
+                msdfInterpolationErrorCorrection(msdf, shape, range, scale, translate, overlapSupport);
                 break;
             case MULTI_AND_TRUE:
                 distanceSignCorrection(mtsdf, shape, scale, translate, fillRule);
                 if (edgeThreshold > 0)
                     msdfErrorCorrection(mtsdf, edgeThreshold/(scale*range));
+                msdfInterpolationErrorCorrection(mtsdf, shape, range, scale, translate, overlapSupport);
                 break;
             default:;
         }

msdfgen.h

@@ -21,8 +21,10 @@
 #include "core/Shape.h"
 #include "core/BitmapRef.hpp"
 #include "core/Bitmap.h"
+#include "core/bitmap-interpolation.hpp"
 #include "core/pixel-conversion.hpp"
 #include "core/edge-coloring.h"
+#include "core/interpolation-error-correction.h"
 #include "core/render-sdf.h"
 #include "core/rasterization.h"
 #include "core/estimate-sdf-error.h"