sdf_bitmap.cc

struct Vec2 {
  float x, y;
};

struct Segment {
  Vec2 a, b;
};

struct OutlineData {
  std::vector<Segment> segments;
  Vec2 pos{};
  Vec2 contour_start{};
};

// ── FT_Outline_Decompose callbacks ──────────────────────────────
// All coordinates arrive in 26.6 fixed-point; convert to float pixels.

static Vec2 FTVec(const FT_Vector* v) {
  return {v->x / 64.0f, v->y / 64.0f};
}

static void AddSeg(OutlineData* o, Vec2 to) {
  float dx = o->pos.x - to.x, dy = o->pos.y - to.y;
  if (dx * dx + dy * dy > 1e-12f) {
    o->segments.push_back({o->pos, to});
  }
  o->pos = to;
}

static int OLMoveTo(const FT_Vector* to, void* user) {
  auto* o = static_cast<OutlineData*>(user);
  AddSeg(o, o->contour_start);  // close previous contour
  o->pos = FTVec(to);
  o->contour_start = o->pos;
  return 0;
}

static int OLLineTo(const FT_Vector* to, void* user) {
  AddSeg(static_cast<OutlineData*>(user), FTVec(to));
  return 0;
}

// ── Adaptive Bézier flattening ──────────────────────────────────
// Subdivision tolerance: a control point within 0.125 px² of the
// chord midpoint is considered flat — roughly 1/3 pixel accuracy.

static void FlattenQuad(OutlineData* o, Vec2 p0, Vec2 p1, Vec2 p2, int d) {
  if (d > 8) {
    AddSeg(o, p2);
    return;
  }
  Vec2 mid = {(p0.x + p2.x) * 0.5f, (p0.y + p2.y) * 0.5f};
  float ex = p1.x - mid.x, ey = p1.y - mid.y;
  if (ex * ex + ey * ey < 0.125f) {
    AddSeg(o, p2);
    return;
  }
  Vec2 q0 = {(p0.x + p1.x) * .5f, (p0.y + p1.y) * .5f};
  Vec2 q1 = {(p1.x + p2.x) * .5f, (p1.y + p2.y) * .5f};
  Vec2 r = {(q0.x + q1.x) * .5f, (q0.y + q1.y) * .5f};
  FlattenQuad(o, p0, q0, r, d + 1);
  o->pos = r;
  FlattenQuad(o, r, q1, p2, d + 1);
}

static int OLConicTo(const FT_Vector* ctrl, const FT_Vector* to, void* user) {
  auto* o = static_cast<OutlineData*>(user);
  FlattenQuad(o, o->pos, FTVec(ctrl), FTVec(to), 0);
  return 0;
}

static void FlattenCubic(OutlineData* o, Vec2 p0, Vec2 p1, Vec2 p2, Vec2 p3,
                         int d) {
  if (d > 8) {
    AddSeg(o, p3);
    return;
  }
  Vec2 mid = {(p0.x + p3.x) * .5f, (p0.y + p3.y) * .5f};
  float d1 = (p1.x - mid.x) * (p1.x - mid.x) +
             (p1.y - mid.y) * (p1.y - mid.y);
  float d2 = (p2.x - mid.x) * (p2.x - mid.x) +
             (p2.y - mid.y) * (p2.y - mid.y);
  if (d1 + d2 < 0.25f) {
    AddSeg(o, p3);
    return;
  }
  Vec2 a = {(p0.x + p1.x) * .5f, (p0.y + p1.y) * .5f};
  Vec2 b = {(p1.x + p2.x) * .5f, (p1.y + p2.y) * .5f};
  Vec2 c = {(p2.x + p3.x) * .5f, (p2.y + p3.y) * .5f};
  Vec2 e = {(a.x + b.x) * .5f, (a.y + b.y) * .5f};
  Vec2 f = {(b.x + c.x) * .5f, (b.y + c.y) * .5f};
  Vec2 g = {(e.x + f.x) * .5f, (e.y + f.y) * .5f};
  FlattenCubic(o, p0, a, e, g, d + 1);
  o->pos = g;
  FlattenCubic(o, g, f, c, p3, d + 1);
}

static int OLCubicTo(const FT_Vector* c1, const FT_Vector* c2,
                     const FT_Vector* to, void* user) {
  auto* o = static_cast<OutlineData*>(user);
  FlattenCubic(o, o->pos, FTVec(c1), FTVec(c2), FTVec(to), 0);
  return 0;
}

// ── Distance to line segment (squared) ──────────────────────────

static float DistToSegSq(Vec2 p, const Segment& s) {
  Vec2 ab = {s.b.x - s.a.x, s.b.y - s.a.y};
  Vec2 ap = {p.x - s.a.x, p.y - s.a.y};
  float len_sq = ab.x * ab.x + ab.y * ab.y;
  float t =
      len_sq > 1e-10f ? (ap.x * ab.x + ap.y * ab.y) / len_sq : 0.0f;
  t = std::clamp(t, 0.0f, 1.0f);
  float dx = p.x - (s.a.x + t * ab.x);
  float dy = p.y - (s.a.y + t * ab.y);
  return dx * dx + dy * dy;
}

static float MinDistSq(Vec2 p, const std::vector<Segment>& segs) {
  float best = 1e18f;
  for (const auto& s : segs) best = std::min(best, DistToSegSq(p, s));
  return best;
}

// ── Winding number ──────────────────────────────────────────────
// Non-zero winding rule: inside when winding != 0.

static int WindingNumber(Vec2 p, const std::vector<Segment>& segs) {
  int wn = 0;
  for (const auto& s : segs) {
    Vec2 a = s.a, b = s.b;
    if (a.y <= p.y) {
      if (b.y > p.y) {
        float cross = (b.x - a.x) * (p.y - a.y) - (p.x - a.x) * (b.y - a.y);
        if (cross > 0) wn++;
      }
    } else {
      if (b.y <= p.y) {
        float cross = (b.x - a.x) * (p.y - a.y) - (p.x - a.x) * (b.y - a.y);
        if (cross < 0) wn--;
      }
    }
  }
  return wn;
}

// ── Meijster EDT ────────────────────────────────────────────────

static void MeijsterEDT(const std::vector<bool>& on, int w, int h,
                         std::vector<float>& out_dsq) {
  const int INF = w + h;
  std::vector<int> g(w * h);

  for (int x = 0; x < w; x++) {
    g[x] = on[x] ? 0 : INF;
    for (int y = 1; y < h; y++) {
      int i = y * w + x;
      g[i] = on[i] ? 0 : g[(y - 1) * w + x] + 1;
    }
    for (int y = h - 2; y >= 0; y--) {
      int i = y * w + x;
      if (g[(y + 1) * w + x] + 1 < g[i]) g[i] = g[(y + 1) * w + x] + 1;
    }
  }

  out_dsq.resize(w * h);
  std::vector<int> s(w), t(w);

  for (int y = 0; y < h; y++) {
    auto f = [&](int x, int i) -> float {
      float gi = static_cast<float>(g[y * w + i]);
      return static_cast<float>((x - i) * (x - i)) + gi * gi;
    };
    auto sep_fn = [&](int i, int u) -> int {
      float gi = static_cast<float>(g[y * w + i]);
      float gu = static_cast<float>(g[y * w + u]);
      return static_cast<int>(std::floor(
          (static_cast<float>(u * u - i * i) + gu * gu - gi * gi) /
          (2.0f * static_cast<float>(u - i))));
    };

    int q = 0;
    s[0] = 0;
    t[0] = 0;
    for (int u = 1; u < w; u++) {
      while (q >= 0 && f(t[q], s[q]) >= f(t[q], u)) q--;
      if (q < 0) {
        q = 0;
        s[0] = u;
        t[0] = 0;
      } else {
        int wv = sep_fn(s[q], u) + 1;
        if (wv < w) {
          q++;
          s[q] = u;
          t[q] = wv;
        }
      }
    }
    for (int u = w - 1; u >= 0; u--) {
      out_dsq[y * w + u] = f(u, s[q]);
      if (u == t[q]) q--;
    }
  }
}

// ── Outline + exact distance SDF ────────────────────────────────

SDFGlyph GenerateSDFOutline(FT_Face face, uint32_t charcode, int pixel_size,
                            int spread) {
  FT_Set_Pixel_Sizes(face, 0, pixel_size);

  FT_UInt glyph_index = FT_Get_Char_Index(face, charcode);
  if (glyph_index == 0) return {};

  FT_Load_Glyph(face, glyph_index, FT_LOAD_NO_BITMAP);
  FT_GlyphSlot slot = face->glyph;

  // ── Decompose outline into flattened line segments ──
  OutlineData outline{};
  FT_Outline_Funcs funcs{};
  funcs.move_to = OLMoveTo;
  funcs.line_to = OLLineTo;
  funcs.conic_to = OLConicTo;
  funcs.cubic_to = OLCubicTo;
  FT_Outline_Decompose(&slot->outline, &funcs, &outline);
  AddSeg(&outline, outline.contour_start);  // close final contour

  if (outline.segments.empty()) return {};

  // ── Grid dimensions ──
  // Metrics are in 26.6 fixed-point.
  int bx = static_cast<int>(slot->metrics.horiBearingX) >> 6;
  int by = static_cast<int>(slot->metrics.horiBearingY) >> 6;
  int gw = (static_cast<int>(slot->metrics.width) + 63) >> 6;
  int gh = (static_cast<int>(slot->metrics.height) + 63) >> 6;

  int pad = spread + 1;
  int w = gw + 2 * pad;
  int h = gh + 2 * pad;

  // SDF texel (ix, iy) maps to glyph-space:
  //   gx = (bx - pad) + ix + 0.5
  //   gy = (by + pad) - iy - 0.5   (FreeType y-up → bitmap y-down)
  float ox = static_cast<float>(bx - pad);
  float oy = static_cast<float>(by + pad);

  // ── Build inside/outside grid via winding number ──
  std::vector<bool> inside(w * h, false);
  for (int iy = 0; iy < h; iy++) {
    for (int ix = 0; ix < w; ix++) {
      float gx = ox + static_cast<float>(ix) + 0.5f;
      float gy = oy - static_cast<float>(iy) - 0.5f;
      inside[iy * w + ix] = WindingNumber({gx, gy}, outline.segments) != 0;
    }
  }

  // ── Build boundary grid: pixels adjacent to an inside/outside transition ──
  std::vector<bool> boundary(w * h, false);
  for (int iy = 0; iy < h; iy++) {
    for (int ix = 0; ix < w; ix++) {
      int i = iy * w + ix;
      bool v = inside[i];
      if ((ix > 0 && inside[i - 1] != v) ||
          (ix < w - 1 && inside[i + 1] != v) ||
          (iy > 0 && inside[i - w] != v) ||
          (iy < h - 1 && inside[i + w] != v)) {
        boundary[i] = true;
      }
    }
  }

  // ── Meijster: unsigned distance² to nearest boundary pixel ──
  std::vector<float> meijster_dsq;
  MeijsterEDT(boundary, w, h, meijster_dsq);

  // ── Combine: exact near-field, Meijster far-field ──
  // Threshold: refine anything Meijster reports within (spread/2)².
  float refine_sq = static_cast<float>(spread * spread) * 0.25f;

  SDFGlyph result;
  result.width = w;
  result.height = h;
  result.bearing_x = bx - pad;
  result.bearing_y = by + pad;
  result.advance_x = static_cast<float>(slot->advance.x) / 64.0f;
  result.pixels.resize(w * h);

  float inv_spread = 1.0f / static_cast<float>(spread);
  for (int iy = 0; iy < h; iy++) {
    for (int ix = 0; ix < w; ix++) {
      int idx = iy * w + ix;
      float dsq;
      if (meijster_dsq[idx] <= refine_sq) {
        // Exact distance to nearest outline segment.
        float gx = ox + static_cast<float>(ix) + 0.5f;
        float gy = oy - static_cast<float>(iy) - 0.5f;
        dsq = MinDistSq({gx, gy}, outline.segments);
      } else {
        dsq = meijster_dsq[idx];
      }
      float d = std::sqrt(dsq);
      float sign = inside[idx] ? 1.0f : -1.0f;
      float norm = (sign * d) * inv_spread * 0.5f + 0.5f;
      result.pixels[idx] =
          static_cast<uint8_t>(std::clamp(norm, 0.0f, 1.0f) * 255.0f);
    }
  }
  return result;
}