sdf.cc

// ── Meijster 2D Euclidean Distance Transform ────────────────────
// Computes squared Euclidean distance from every pixel to the nearest
// "on" pixel in the binary grid.  Two-pass, O(n) algorithm.

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);

  // Phase 1: column scan — g[y][x] = vertical distance to nearest on-pixel.
  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;
    }
  }

  // Phase 2: row scan — parabola envelope.
  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 = [&](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(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--;
    }
  }
}

// ── Nearest-edge propagation ────────────────────────────────────
// 8-connected BFS from edge pixels to propagate a Voronoi label so
// that every non-edge pixel knows which edge pixel is closest.

static void PropagateNearestEdge(const std::vector<bool>& is_edge, int w,
                                 int h, std::vector<int>& nearest) {
  nearest.assign(w * h, -1);
  std::deque<int> queue;
  for (int i = 0; i < w * h; i++) {
    if (is_edge[i]) {
      nearest[i] = i;
      queue.push_back(i);
    }
  }
  static constexpr int dx[] = {-1, 1, 0, 0, -1, -1, 1, 1};
  static constexpr int dy[] = {0, 0, -1, 1, -1, 1, -1, 1};
  while (!queue.empty()) {
    int idx = queue.front();
    queue.pop_front();
    int x = idx % w, y = idx / w;
    for (int d = 0; d < 8; d++) {
      int nx = x + dx[d], ny = y + dy[d];
      if (nx < 0 || nx >= w || ny < 0 || ny >= h) continue;
      int ni = ny * w + nx;
      if (nearest[ni] >= 0) continue;
      nearest[ni] = nearest[idx];
      queue.push_back(ni);
    }
  }
}

// ── Grayscale bitmap + Meijster SDF ─────────────────────────────

SDFGlyph GenerateSDFBitmap(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_DEFAULT);
  FT_Render_Glyph(face->glyph, FT_RENDER_MODE_NORMAL);

  FT_GlyphSlot slot = face->glyph;
  FT_Bitmap& bmp = slot->bitmap;
  int bw = static_cast<int>(bmp.width);
  int bh = static_cast<int>(bmp.rows);

  int pad = spread + 1;
  int w = bw + 2 * pad;
  int h = bh + 2 * pad;

  // ── Build padded coverage grid [0.0, 1.0] ──
  std::vector<float> coverage(w * h, 0.0f);
  for (int y = 0; y < bh; y++) {
    for (int x = 0; x < bw; x++) {
      coverage[(y + pad) * w + (x + pad)] =
          bmp.buffer[y * bmp.pitch + x] / 255.0f;
    }
  }

  // ── Classify pixels ──
  enum class Kind : uint8_t { Inside, Outside, Edge };
  std::vector<Kind> kind(w * h);
  constexpr float KEdgeLo = 1.0f / 255.0f;
  constexpr float KEdgeHi = 254.0f / 255.0f;
  for (int i = 0; i < w * h; i++) {
    if (coverage[i] <= KEdgeLo)
      kind[i] = Kind::Outside;
    else if (coverage[i] >= KEdgeHi)
      kind[i] = Kind::Inside;
    else
      kind[i] = Kind::Edge;
  }

  // ── Sub-pixel edge distance using coverage gradient ──
  // For edge pixels the gradient magnitude corrects for edge angle:
  //   d = (coverage - 0.5) / |∇coverage|
  // This is the core insight from Gustavson's EDTAA3 algorithm.
  std::vector<float> edge_dist(w * h, 0.0f);
  for (int y = 1; y < h - 1; y++) {
    for (int x = 1; x < w - 1; x++) {
      int i = y * w + x;
      if (kind[i] != Kind::Edge) continue;
      float gx = (coverage[i + 1] - coverage[i - 1]) * 0.5f;
      float gy = (coverage[i + w] - coverage[i - w]) * 0.5f;
      float glen = std::sqrt(gx * gx + gy * gy);
      if (glen > 0.001f) {
        edge_dist[i] = (coverage[i] - 0.5f) / glen;
      } else {
        edge_dist[i] = (coverage[i] - 0.5f) * 2.0f;
      }
      edge_dist[i] = std::clamp(edge_dist[i], -1.0f, 1.0f);
    }
  }

  // ── Meijster EDT from edge pixels ──
  std::vector<bool> is_edge(w * h, false);
  for (int i = 0; i < w * h; i++) is_edge[i] = (kind[i] == Kind::Edge);

  std::vector<float> meijster_dsq;
  MeijsterEDT(is_edge, w, h, meijster_dsq);

  // ── Propagate nearest-edge labels ──
  std::vector<int> nearest;
  PropagateNearestEdge(is_edge, w, h, nearest);

  // ── Combine into signed distance ──
  SDFGlyph result;
  result.width = w;
  result.height = h;
  result.bearing_x = slot->bitmap_left - pad;
  result.bearing_y = slot->bitmap_top + 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 i = 0; i < w * h; i++) {
    float d;
    if (kind[i] == Kind::Edge) {
      d = edge_dist[i];
    } else {
      float grid_dist = std::sqrt(meijster_dsq[i]);
      float sign = (kind[i] == Kind::Inside) ? 1.0f : -1.0f;
      float sub_pixel = 0.0f;
      if (nearest[i] >= 0) sub_pixel = edge_dist[nearest[i]];
      d = sign * (grid_dist + sign * sub_pixel);
    }
    float norm = d * inv_spread * 0.5f + 0.5f;
    result.pixels[i] =
        static_cast<uint8_t>(std::clamp(norm, 0.0f, 1.0f) * 255.0f);
  }
  return result;
}