Draw text from an atlas the CPU fills, composed by the GPU

Each size's symbols are drawn on the CPU into one 8-bit atlas, every pixel
the exact area its curves cover, and text is a quad per glyph reading it.
A 1080p page of 11,382 glyphs composes in 1.07 ms at 16 px (it took 525 ms
computing coverage per fragment), and the composed pixels are the CPU
raster's exactly (worst difference 0 for DejaVu's ASCII, 1 for Noto's
11,172 Hangul syllables).

Computing coverage on the GPU lost on both counts: a compute pass filling
the atlas took 5-13x longer than one Cortex-A76 core (35 against 2.7 ms at
16 px), and its faster sorting-network form computed whole rows wrongly
under this Mesa V3D. Coverage now has one implementation, in C, and the
shaders only compose. vega_symbol_raster takes a row stride to draw into
the atlas in place.

5 files changed, 239 insertions, 254 deletions

diff --git a/symbols/symbols-gl.c b/symbols/symbols-gl.c
index 0e78efe..c1ce64e 100644
--- a/symbols/symbols-gl.c
+++ b/symbols/symbols-gl.c
@@ -1,126 +1,34 @@
 #include "symbols-gl.h"
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 
-enum { BANDS = 16 };
+enum { ATLAS_MAX = 4096 };
 
-/* Coordinates travel doubled, so implied points stay whole: 2x of a point
- * is its raw i16 without the flag, and an implied point's is the sum of its
- * two neighbours' coordinates. */
 #define SHARED \
     "#version 310 es\n" \
     "precision highp float;\n" \
     "precision highp int;\n"
 
-static const char vertex_source[] = SHARED
-    "layout(location = 0) in vec4 place;   /* origin x, y in pixels; pixels per unit; unused */\n"
-    "layout(location = 1) in vec4 bounds;  /* left, bottom, right, top in units */\n"
-    "layout(location = 2) in uint bands;   /* its first band */\n"
+/* Text: a quad per glyph at whole pixels, each pixel one atlas texel. */
+static const char draw_vertex[] = SHARED
+    "layout(location = 0) in ivec4 rect; /* in the atlas */\n"
+    "layout(location = 1) in ivec2 at;   /* its top left on the target */\n"
     "uniform vec2 target;\n"
-    "out vec2 units;\n"
-    "flat out float scale;\n"
-    "flat out ivec2 rows;                 /* bottom and top in units */\n"
-    "flat out uint base;\n"
+    "out vec2 texel;\n"
     "void main() {\n"
-    "    vec2 corner = vec2(float(gl_VertexID & 1), float(gl_VertexID >> 1));\n"
-    "    vec2 pad = vec2(1.0 / place.z);\n"
-    "    units = mix(bounds.xy - pad, bounds.zw + pad, corner);\n"
-    "    vec2 pixel = place.xy + vec2(units.x, -units.y) * place.z;\n"
+    "    vec2 corner = vec2(float(gl_VertexID & 1), float(gl_VertexID >> 1)), size = vec2(rect.zw);\n"
+    "    vec2 pixel = vec2(at) + corner * size;\n"
+    "    texel = vec2(rect.xy) + corner * size;\n"
     "    gl_Position = vec4(pixel.x / target.x * 2.0 - 1.0, 1.0 - pixel.y / target.y * 2.0, 0.0, 1.0);\n"
-    "    scale = place.z;\n"
-    "    rows = ivec2(bounds.yw);\n"
-    "    base = bands;\n"
     "}\n";
 
-static const char fragment_source[] = SHARED
-    "layout(std430, binding = 0) readonly buffer Points { uint points[]; };\n"
-    "layout(std430, binding = 1) readonly buffer Bands { uint bands[]; };\n"
-    "uniform uint listed; /* where the lists of curves start in bands */\n"
-    "in vec2 units;\n"
-    "flat in float scale;\n"
-    "flat in ivec2 rows;\n"
-    "flat in uint base;\n"
+static const char draw_fragment[] = SHARED
+    "uniform highp sampler2D atlas;\n"
+    "in vec2 texel;\n"
     "out vec4 color;\n"
-    "\n"
-    "ivec3 point(uint i) { /* doubled x, doubled y, off */\n"
-    "    uint w = points[i];\n"
-    "    int x = int(w << 16u) >> 16, y = int(w) >> 16;\n"
-    "    return ivec3(x & ~1, y & ~1, x & 1);\n"
-    "}\n"
-    "int band(int doubled) { /* which of 16 bands a doubled y lies in */\n"
-    "    return clamp((doubled - 2 * rows.x) * 16 / (2 * (rows.y - rows.x) + 2), 0, 15);\n"
-    "}\n"
-    "int roots(float a, float b, float c, inout float ts[10], int n) {\n"
-    "    float r0 = -1.0, r1 = -1.0;\n"
-    "    if (a == 0.0) { if (b != 0.0) r0 = -c / b; }\n"
-    "    else {\n"
-    "        float d = b * b - 4.0 * a * c;\n"
-    "        if (d >= 0.0) {\n"
-    "            float q = -0.5 * (b + (b >= 0.0 ? sqrt(d) : -sqrt(d)));\n"
-    "            r0 = q / a;\n"
-    "            if (q != 0.0) r1 = c / q;\n"
-    "        }\n"
-    "    }\n"
-    "    if (r0 > 0.0 && r0 < 1.0) ts[n++] = r0;\n"
-    "    if (r1 > 0.0 && r1 < 1.0) ts[n++] = r1;\n"
-    "    return n;\n"
-    "}\n"
-    "/* As curve_area in symbols.c: the curve's share of the pixel [0, 1]^2. */\n"
-    "float area(vec2 p0, vec2 p1, vec2 p2) {\n"
-    "    if (max(max(p0.y, p1.y), p2.y) <= 0.0 || min(min(p0.y, p1.y), p2.y) >= 1.0 || min(min(p0.x, p1.x), p2.x) >= 1.0)\n"
-    "        return 0.0;\n"
-    "    if (max(max(p0.x, p1.x), p2.x) <= 0.0) return clamp(p2.y, 0.0, 1.0) - clamp(p0.y, 0.0, 1.0);\n"
-    "    vec2 a = p0 - 2.0 * p1 + p2, b = p0 - p1;\n"
-    "    float ts[10];\n"
-    "    ts[0] = 0.0;\n"
-    "    int n = 1;\n"
-    "    n = roots(a.y, -2.0 * b.y, p0.y, ts, n);\n"
-    "    n = roots(a.y, -2.0 * b.y, p0.y - 1.0, ts, n);\n"
-    "    n = roots(a.x, -2.0 * b.x, p0.x, ts, n);\n"
-    "    n = roots(a.x, -2.0 * b.x, p0.x - 1.0, ts, n);\n"
-    "    ts[n++] = 1.0;\n"
-    "    for (int i = 1; i < n; ++i)\n"
-    "        for (int j = i; j > 0 && ts[j] < ts[j - 1]; --j) { float t = ts[j]; ts[j] = ts[j - 1]; ts[j - 1] = t; }\n"
-    "    float u0 = 1.0 - p0.x, u1 = 2.0 * b.x, u2 = -a.x, v0 = -2.0 * b.y, v1 = 2.0 * a.y;\n"
-    "    float k1 = u0 * v0, k2 = (u0 * v1 + u1 * v0) / 2.0, k3 = (u1 * v1 + u2 * v0) / 3.0, k4 = u2 * v1 / 4.0, sum = 0.0;\n"
-    "    for (int i = 0; i + 1 < n; ++i) {\n"
-    "        float ta = ts[i], tb = ts[i + 1], t = (ta + tb) / 2.0;\n"
-    "        if (tb <= ta) continue;\n"
-    "        float y = (a.y * t - 2.0 * b.y) * t + p0.y, x = (a.x * t - 2.0 * b.x) * t + p0.x;\n"
-    "        if (y <= 0.0 || y >= 1.0 || x >= 1.0) continue;\n"
-    "        if (x <= 0.0) sum += (a.y * tb - 2.0 * b.y) * tb - (a.y * ta - 2.0 * b.y) * ta;\n"
-    "        else sum += (((k4 * tb + k3) * tb + k2) * tb + k1) * tb - (((k4 * ta + k3) * ta + k2) * ta + k1) * ta;\n"
-    "    }\n"
-    "    return sum;\n"
-    "}\n"
-    "void main() {\n"
-    "    vec2 corner = units - vec2(0.5 / scale); /* the pixel's lower left, in units */\n"
-    "    int b0 = band(2 * int(floor(corner.y))), b1 = band(2 * int(ceil(corner.y + 1.0 / scale)));\n"
-    "    float sum = 0.0;\n"
-    "    for (int b = b0; b <= b1; ++b) {\n"
-    "        uint at = bands[base + 2u * uint(b)], end = at + bands[base + 2u * uint(b) + 1u];\n"
-    "        for (; at < end; ++at) {\n"
-    "            uint k = bands[listed + at];\n"
-    "            ivec3 p = point(k), n = point(k + 1u);\n"
-    "            ivec2 s, c, e;\n"
-    "            if (p.z == 1) {\n"
-    "                ivec3 q = point(k - 1u);\n"
-    "                s = q.z == 1 ? (q.xy + p.xy) / 2 : q.xy;\n"
-    "                c = p.xy;\n"
-    "                e = n.z == 1 ? (p.xy + n.xy) / 2 : n.xy;\n"
-    "            } else {\n"
-    "                s = p.xy;\n"
-    "                c = (p.xy + n.xy) / 2;\n"
-    "                e = n.xy;\n"
-    "            }\n"
-    "            /* a curve over several of this pixel's bands counts in the first */\n"
-    "            if (max(band(min(min(s.y, c.y), e.y)), b0) != b) continue;\n"
-    "            sum += area((vec2(s) * 0.5 - corner) * scale, (vec2(c) * 0.5 - corner) * scale, (vec2(e) * 0.5 - corner) * scale);\n"
-    "        }\n"
-    "    }\n"
-    "    color = vec4(min(abs(sum), 1.0), 0.0, 0.0, 1.0);\n"
-    "}\n";
+    "void main() { color = vec4(texelFetch(atlas, ivec2(texel), 0).r, 0.0, 0.0, 1.0); }\n";
 
 static GLuint compile(GLenum type, const char *source, char *error, size_t error_size) {
     GLuint shader = glCreateShader(type);
@@ -136,125 +44,162 @@ static GLuint compile(GLenum type, const char *source, char *error, size_t error
     return shader;
 }
 
-/* Doubled coordinates make implied points whole; bands as the shader cuts them. */
-static int band(int doubled, int bottom, int top) {
-    int b = (doubled - 2 * bottom) * BANDS / (2 * (top - bottom) + 2);
-    return b < 0 ? 0 : b >= BANDS ? BANDS - 1 : b;
-}
-
-bool vega_gl_open(VegaGl *gl, const VegaSymbols *s, char *error, size_t error_size) {
+bool vega_gl_open(VegaGl *gl, char *error, size_t error_size) {
     *gl = (VegaGl){0};
-    GLuint vertex = compile(GL_VERTEX_SHADER, vertex_source, error, error_size);
-    GLuint fragment = vertex ? compile(GL_FRAGMENT_SHADER, fragment_source, error, error_size) : 0;
-    if (!fragment) return false;
-    gl->program = glCreateProgram();
-    glAttachShader(gl->program, vertex);
-    glAttachShader(gl->program, fragment);
-    glLinkProgram(gl->program);
+    GLuint vertex = compile(GL_VERTEX_SHADER, draw_vertex, error, error_size);
+    GLuint fragment = vertex ? compile(GL_FRAGMENT_SHADER, draw_fragment, error, error_size) : 0;
+    if (!fragment) {
+        if (vertex) glDeleteShader(vertex);
+        return false;
+    }
+    gl->draw = glCreateProgram();
+    glAttachShader(gl->draw, vertex);
+    glAttachShader(gl->draw, fragment);
+    glLinkProgram(gl->draw);
     glDeleteShader(vertex);
     glDeleteShader(fragment);
     GLint ok;
-    glGetProgramiv(gl->program, GL_LINK_STATUS, &ok);
+    glGetProgramiv(gl->draw, GL_LINK_STATUS, &ok);
     if (!ok) {
-        glGetProgramInfoLog(gl->program, (GLsizei)error_size, NULL, error);
+        glGetProgramInfoLog(gl->draw, (GLsizei)error_size, NULL, error);
+        vega_gl_close(gl);
         return false;
     }
-    /* The point table, as the file has it. */
-    glGenBuffers(1, &gl->points);
-    glBindBuffer(GL_SHADER_STORAGE_BUFFER, gl->points);
-    glBufferData(GL_SHADER_STORAGE_BUFFER, (GLsizeiptr)(s->points ? s->points : 1) * VEGA_POINT_BYTES,
-                 s->bytes + VEGA_SYMBOLS_HEADER + (size_t)s->symbols * VEGA_SYMBOL_BYTES, GL_STATIC_DRAW);
-    /* Bands: per symbol, 16 (start, count) pairs into one list of curves,
-     * each curve (by the point that begins it) in every band it reaches. */
-    size_t table = (size_t)s->symbols * BANDS * 2, listed = 0;
-    uint32_t *counts = calloc(table + 1, sizeof(*counts));
-    for (int pass = 0; pass < 2; ++pass) {
-        uint32_t *bands = pass ? malloc((table + listed + 1) * sizeof(*bands)) : NULL;
-        size_t next = 0;
-        for (uint32_t i = 0; i < s->symbols; ++i) {
-            VegaSymbol symbol = vega_symbol_at(s, i);
-            uint32_t *mine = pass ? bands + (size_t)i * BANDS * 2 : counts + (size_t)i * BANDS * 2;
-            if (pass)
-                for (int b = 0; b < BANDS; ++b) {
-                    mine[2 * b] = (uint32_t)next;
-                    next += counts[(size_t)i * BANDS * 2 + 2 * b + 1];
-                    mine[2 * b + 1] = 0;
-                }
-            VegaCurve c;
-            for (uint32_t k = symbol.first; k < symbol.first + symbol.count; ++k) {
-                if (!vega_symbol_curve(s, &symbol, k, &c)) continue;
-                float low = c.y0 < c.y1 ? c.y0 : c.y1, high = c.y0 > c.y1 ? c.y0 : c.y1;
-                low = c.y2 < low ? c.y2 : low;
-                high = c.y2 > high ? c.y2 : high;
-                int from = band((int)(2 * low), symbol.bottom, symbol.top), to = band((int)(2 * high), symbol.bottom, symbol.top);
-                for (int b = from; b <= to; ++b) {
-                    if (pass) bands[table + mine[2 * b] + mine[2 * b + 1]] = k;
-                    ++mine[2 * b + 1];
-                    listed += !pass;
-                }
-            }
-        }
-        if (pass) {
-            glGenBuffers(1, &gl->bands);
-            glBindBuffer(GL_SHADER_STORAGE_BUFFER, gl->bands);
-            glBufferData(GL_SHADER_STORAGE_BUFFER, (GLsizeiptr)((table + listed + 1) * sizeof(*bands)), bands, GL_STATIC_DRAW);
-            free(bands);
-        }
-    }
-    free(counts);
+    glUseProgram(gl->draw);
+    glUniform1i(glGetUniformLocation(gl->draw, "atlas"), 0);
+    gl->target = glGetUniformLocation(gl->draw, "target");
     glGenVertexArrays(1, &gl->array);
     glGenBuffers(1, &gl->instances);
     glBindVertexArray(gl->array);
     glBindBuffer(GL_ARRAY_BUFFER, gl->instances);
     glEnableVertexAttribArray(0);
     glEnableVertexAttribArray(1);
-    glEnableVertexAttribArray(2);
-    glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 36, (void *)0);
-    glVertexAttribPointer(1, 4, GL_FLOAT, GL_FALSE, 36, (void *)16);
-    glVertexAttribIPointer(2, 1, GL_UNSIGNED_INT, 36, (void *)32);
+    glVertexAttribIPointer(0, 4, GL_INT, 24, (void *)0);
+    glVertexAttribIPointer(1, 2, GL_INT, 24, (void *)16);
     glVertexAttribDivisor(0, 1);
     glVertexAttribDivisor(1, 1);
-    glVertexAttribDivisor(2, 1);
     glBindVertexArray(0);
-    glUseProgram(gl->program);
-    glUniform1ui(glGetUniformLocation(gl->program, "listed"), (GLuint)table);
-    gl->viewport = glGetUniformLocation(gl->program, "target");
     if (glGetError() != GL_NO_ERROR) {
-        snprintf(error, error_size, "GL error while uploading");
+        snprintf(error, error_size, "cannot set up the text program");
+        vega_gl_close(gl);
         return false;
     }
     return true;
 }
 
-void vega_gl_draw(VegaGl *gl, const VegaSymbols *s, const VegaGlyph *glyphs, uint32_t count, int width, int height) {
-    float *instances = malloc((size_t)(count ? count : 1) * 36);
+/* The pixels a symbol reaches at scale, one more on each side. */
+static VegaPlace place_of(const VegaSymbol *symbol, float scale) {
+    if (!symbol->count) return (VegaPlace){0};
+    int left = (int)floorf(symbol->left * scale) - 1, right = (int)ceilf(symbol->right * scale) + 1;
+    int top = (int)ceilf(symbol->top * scale) + 1, bottom = (int)floorf(symbol->bottom * scale) - 1;
+    return (VegaPlace){0, 0, right - left, top - bottom, left, top};
+}
+
+bool vega_gl_atlas(VegaGl *gl, const VegaSymbols *s, float scale, const uint32_t *wanted, uint32_t count, VegaAtlas *out,
+                   char *error, size_t error_size) {
+    (void)gl;
+    *out = (VegaAtlas){.scale = scale, .symbols = s->symbols};
+    if (!wanted) count = s->symbols;
+    out->places = calloc(s->symbols ? s->symbols : 1, sizeof(*out->places));
+    if (!out->places) {
+        snprintf(error, error_size, "out of memory");
+        return false;
+    }
+    /* Shelves, as wide as the square the rectangles would fill (a power of
+     * two), the next shelf starting under the tallest of the one before. */
+    uint64_t area = 0;
     for (uint32_t i = 0; i < count; ++i) {
-        VegaSymbol symbol = vega_symbol_at(s, glyphs[i].symbol);
-        float *f = instances + (size_t)i * 9;
-        f[0] = glyphs[i].x, f[1] = glyphs[i].y, f[2] = glyphs[i].scale, f[3] = 0;
-        f[4] = symbol.left, f[5] = symbol.bottom, f[6] = symbol.right, f[7] = symbol.top;
-        uint32_t base = glyphs[i].symbol * BANDS * 2;
-        memcpy(f + 8, &base, 4);
+        uint32_t index = wanted ? wanted[i] : i;
+        if (index >= s->symbols) continue;
+        VegaSymbol symbol = vega_symbol_at(s, index);
+        VegaPlace p = place_of(&symbol, scale);
+        area += (uint64_t)p.width * (uint64_t)p.height;
     }
-    glUseProgram(gl->program);
-    glUniform2f(gl->viewport, (float)width, (float)height);
-    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, gl->points);
-    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, gl->bands);
+    int width = 64;
+    while (width < ATLAS_MAX && (uint64_t)width * (uint64_t)width < area + area / 4) width *= 2;
+    int x = 0, y = 0, shelf = 0;
+    for (uint32_t i = 0; i < count; ++i) {
+        uint32_t index = wanted ? wanted[i] : i;
+        if (index >= s->symbols) continue;
+        VegaSymbol symbol = vega_symbol_at(s, index);
+        VegaPlace p = place_of(&symbol, scale);
+        if (!p.width) continue;
+        if (p.width > width) width = p.width;
+        if (x + p.width > width) { x = 0; y += shelf; shelf = 0; }
+        p.x = x, p.y = y;
+        x += p.width;
+        shelf = p.height > shelf ? p.height : shelf;
+        out->places[index] = p;
+    }
+    int height = y + shelf;
+    if (width > ATLAS_MAX || height > ATLAS_MAX) {
+        snprintf(error, error_size, "%u symbols at %.4f pixels a unit need %dx%d, more than one %d-pixel atlas", count,
+                 (double)scale, width, height, ATLAS_MAX);
+        vega_gl_atlas_close(out);
+        return false;
+    }
+    out->width = width, out->height = height ? height : 1;
+    unsigned char *pixels = calloc((size_t)out->width * (size_t)out->height, 1);
+    if (!pixels) {
+        snprintf(error, error_size, "out of memory");
+        vega_gl_atlas_close(out);
+        return false;
+    }
+    for (uint32_t i = 0; i < s->symbols; ++i) {
+        VegaPlace p = out->places[i];
+        if (!p.width) continue;
+        VegaSymbol symbol = vega_symbol_at(s, i);
+        vega_symbol_raster(s, &symbol, scale, p.left / scale, p.top / scale, p.width, p.height,
+                           pixels + (size_t)p.y * (size_t)out->width + (size_t)p.x, (size_t)out->width);
+    }
+    glGenTextures(1, &out->texture);
+    glBindTexture(GL_TEXTURE_2D, out->texture);
+    glTexStorage2D(GL_TEXTURE_2D, 1, GL_R8, out->width, out->height);
+    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
+    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, out->width, out->height, GL_RED, GL_UNSIGNED_BYTE, pixels);
+    free(pixels);
+    if (glGetError() != GL_NO_ERROR) {
+        snprintf(error, error_size, "cannot upload the atlas");
+        vega_gl_atlas_close(out);
+        return false;
+    }
+    return true;
+}
+
+void vega_gl_draw(VegaGl *gl, const VegaAtlas *atlas, const VegaGlyph *glyphs, uint32_t count, int width, int height) {
+    int32_t *instances = malloc((size_t)(count ? count : 1) * 6 * sizeof(*instances));
+    uint32_t n = 0;
+    for (uint32_t i = 0; instances && i < count; ++i) {
+        if (glyphs[i].symbol >= atlas->symbols) continue;
+        VegaPlace p = atlas->places[glyphs[i].symbol];
+        if (!p.width) continue;
+        int32_t *q = instances + (size_t)n++ * 6;
+        q[0] = p.x, q[1] = p.y, q[2] = p.width, q[3] = p.height;
+        q[4] = glyphs[i].x + p.left, q[5] = glyphs[i].y - p.top;
+    }
+    glUseProgram(gl->draw);
+    glUniform2f(gl->target, (float)width, (float)height);
+    glActiveTexture(GL_TEXTURE0);
+    glBindTexture(GL_TEXTURE_2D, atlas->texture);
     glBindVertexArray(gl->array);
     glBindBuffer(GL_ARRAY_BUFFER, gl->instances);
-    glBufferData(GL_ARRAY_BUFFER, (GLsizeiptr)count * 36, instances, GL_STREAM_DRAW);
+    glBufferData(GL_ARRAY_BUFFER, (GLsizeiptr)((size_t)n * 6 * sizeof(*instances)), instances, GL_STREAM_DRAW);
     glEnable(GL_BLEND);
     glBlendFunc(GL_ONE, GL_ONE);
-    glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, (GLsizei)count);
+    if (n) glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, (GLsizei)n);
     glBindVertexArray(0);
     free(instances);
 }
 
+void vega_gl_atlas_close(VegaAtlas *atlas) {
+    if (atlas->texture) glDeleteTextures(1, &atlas->texture);
+    free(atlas->places);
+    *atlas = (VegaAtlas){0};
+}
+
 void vega_gl_close(VegaGl *gl) {
-    glDeleteBuffers(1, &gl->points);
-    glDeleteBuffers(1, &gl->bands);
-    glDeleteBuffers(1, &gl->instances);
-    glDeleteVertexArrays(1, &gl->array);
-    glDeleteProgram(gl->program);
+    if (gl->instances) glDeleteBuffers(1, &gl->instances);
+    if (gl->array) glDeleteVertexArrays(1, &gl->array);
+    if (gl->draw) glDeleteProgram(gl->draw);
     *gl = (VegaGl){0};
 }
diff --git a/symbols/symbols-gl.h b/symbols/symbols-gl.h
index 4c6db67..919a114 100644
--- a/symbols/symbols-gl.h
+++ b/symbols/symbols-gl.h
@@ -3,22 +3,41 @@
 #include "symbols.h"
 #include <GLES3/gl31.h>
 
-/* Symbols drawn by the GPU (OpenGL ES 3.1): the file's point table goes up
- * as it is, one storage buffer, and each pixel of a symbol's box sums the
- * exact area of the curves crossing its row, as vega_coverage does. Rows
- * are cut into bands at upload, each listing its curves, so a pixel reads
- * only those. Needs a current context. */
+/* Symbols on the GPU (OpenGL ES 3.1): the symbols wanted at one scale are
+ * drawn on the CPU into an atlas, each pixel the exact area its curves cover
+ * (vega_symbol_raster), and text is quads reading it. On the Pi 5 the CPU
+ * draws the atlas faster than the GPU's compute could, and the GPU composes
+ * a 1080p page of text in about a millisecond (docs/design/002). A scale that
+ * keeps changing (a zoom) draws its atlas again. Needs a current context. */
 typedef struct {
-    GLuint program, points, bands, instances, array;
-    GLint viewport;
+    GLuint draw, array, instances;
+    GLint target;
 } VegaGl;
 
-/* A symbol placed: its origin in pixels from the target's top left, and
- * pixels per unit. */
-typedef struct { float x, y, scale; uint32_t symbol; } VegaGlyph;
+/* Where each symbol is in an atlas: its rectangle there (x, y, width, height)
+ * and the pixel offset of that rectangle's top left from the symbol's origin
+ * (right, up). A symbol with nothing to draw has an empty rectangle. */
+typedef struct { int32_t x, y, width, height, left, top; } VegaPlace;
 
-bool vega_gl_open(VegaGl *, const VegaSymbols *, char *error, size_t error_size);
+typedef struct {
+    GLuint texture; /* 8-bit coverage */
+    int width, height;
+    float scale;    /* pixels per unit */
+    uint32_t symbols;
+    VegaPlace *places; /* per symbol index */
+} VegaAtlas;
+
+/* A symbol at a pen position (whole pixels from the target's top left, on the
+ * baseline). */
+typedef struct { int32_t x, y; uint32_t symbol; } VegaGlyph;
+
+bool vega_gl_open(VegaGl *, char *error, size_t error_size);
+/* Draw symbols (their indices; all when wanted is NULL) at scale into a new
+ * atlas; false when they do not fit one texture. */
+bool vega_gl_atlas(VegaGl *, const VegaSymbols *, float scale, const uint32_t *wanted, uint32_t count, VegaAtlas *out,
+                   char *error, size_t error_size);
 /* Coverage added into the red channel of the bound framebuffer. */
-void vega_gl_draw(VegaGl *, const VegaSymbols *, const VegaGlyph *, uint32_t count, int width, int height);
+void vega_gl_draw(VegaGl *, const VegaAtlas *, const VegaGlyph *, uint32_t count, int width, int height);
+void vega_gl_atlas_close(VegaAtlas *);
 void vega_gl_close(VegaGl *);
 #endif
diff --git a/symbols/symbols-tool.c b/symbols/symbols-tool.c
index a94feb3..b003978 100644
--- a/symbols/symbols-tool.c
+++ b/symbols/symbols-tool.c
@@ -667,7 +667,7 @@ static int compare(int argc, char **argv) {
             if (need > raster_size && !(raster = realloc(raster, raster_size = need))) fail("out of memory");
             double started = seconds();
             vega_symbol_raster(&s, &symbol, scale, face->glyph->bitmap_left / scale, face->glyph->bitmap_top / scale,
-                               (int)b->width, (int)b->rows, raster);
+                               (int)b->width, (int)b->rows, raster, b->width);
             ours_seconds += seconds() - started;
             for (unsigned j = 0; j < b->rows; ++j)
                 for (unsigned i2 = 0; i2 < b->width; ++i2) {
@@ -730,7 +730,7 @@ static int render(int argc, char **argv) {
     int width = (int)ceil((symbol.advance > symbol.right ? symbol.advance : symbol.right) * scale) + 2;
     int height = (int)ceil(pixels * 1.25) + 2;
     unsigned char *raster = malloc((size_t)width * (size_t)height);
-    vega_symbol_raster(&s, &symbol, scale, -1 / scale, (pixels + 1) / scale, width, height, raster);
+    vega_symbol_raster(&s, &symbol, scale, -1 / scale, (pixels + 1) / scale, width, height, raster, (size_t)width);
     FILE *f = fopen(argv[5], "wb");
     if (!f) fail("cannot write");
     fprintf(f, "P5\n%d %d\n255\n", width, height);
@@ -775,9 +775,10 @@ static GLuint gl_target(int width, int height) {
     return texture;
 }
 
-/* Every symbol in a grid of cells, drawn by the GPU and by vega_coverage:
- * how far apart; then pages of text (glyphs cycling through the file)
- * drawn again and again: how long a page takes. */
+/* At each size: every symbol drawn into an atlas and composed by the GPU at a
+ * pen position, read back against vega_symbol_raster over the same pixels;
+ * then pages of text (glyphs cycling through the file) composed again and
+ * again. */
 static int gpu(int argc, char **argv) {
     if (argc < 5 || strcmp(argv[3], "--sizes")) fail("usage: vega-symbols gpu SYMBOLS.vsym --sizes 12,16,32,64 [--frames N]");
     int frames = argc > 6 && !strcmp(argv[5], "--frames") ? atoi(argv[6]) : 20;
@@ -788,7 +789,7 @@ static int gpu(int argc, char **argv) {
     printf("{\"renderer\":\"%s\"}\n", (const char *)glGetString(GL_RENDERER));
     VegaGl gl;
     double started = seconds();
-    if (!vega_gl_open(&gl, &s, error, sizeof(error))) fail(error);
+    if (!vega_gl_open(&gl, error, sizeof(error))) fail(error);
     glFinish();
     double open_ms = (seconds() - started) * 1e3;
     for (const char *p = argv[4]; *p;) {
@@ -796,72 +797,91 @@ static int gpu(int argc, char **argv) {
         long size = strtol(p, &end, 10);
         p = *end == ',' ? end + 1 : end;
         if (size <= 0 || size > 512) fail("--sizes");
-        double scale = (double)size / s.units;
-        /* Cells that hold any symbol, so none reaches into another's. */
+        float scale = (float)size / s.units;
+        VegaAtlas atlas;
+        started = seconds();
+        if (!vega_gl_atlas(&gl, &s, scale, NULL, 0, &atlas, error, sizeof(error))) {
+            printf("{\"pixels_per_em\":%ld,\"error\":\"%s\"}\n", size, error);
+            continue;
+        }
+        glFinish();
+        double fill_ms = (seconds() - started) * 1e3;
+        /* Every symbol composed at a pen position of its own cell, read back
+         * and set against vega_symbol_raster at the same pixels. */
         int left = 0, bottom = 0, right = 0, top = 0;
         for (uint32_t i = 0; i < s.symbols; ++i) {
-            VegaSymbol symbol = vega_symbol_at(&s, i);
-            left = symbol.left < left ? symbol.left : left;
-            bottom = symbol.bottom < bottom ? symbol.bottom : bottom;
-            right = symbol.right > right ? symbol.right : right;
-            top = symbol.top > top ? symbol.top : top;
+            VegaPlace place = atlas.places[i];
+            if (!place.width) continue;
+            left = place.left < left ? place.left : left;
+            top = place.top > top ? place.top : top;
+            right = place.left + place.width > right ? place.left + place.width : right;
+            bottom = place.top - place.height < bottom ? place.top - place.height : bottom;
         }
-        int cell = (int)ceil((right - left > top - bottom ? right - left : top - bottom) * scale) + 4;
-        int across = 4096 / cell, rows = (int)((s.symbols + (uint32_t)across - 1) / (uint32_t)across);
-        if (rows * cell > 4096) rows = 4096 / cell;
-        int width = across * cell, height = rows * cell;
+        int cell_width = right - left + 2, cell_height = top - bottom + 2;
+        int across = 2048 / cell_width, rows = (int)((s.symbols + (uint32_t)across - 1) / (uint32_t)across);
+        if (rows * cell_height > 4096) rows = 4096 / cell_height;
+        int grid_width = across * cell_width, grid_height = rows * cell_height;
         uint32_t shown = s.symbols < (uint32_t)(across * rows) ? s.symbols : (uint32_t)(across * rows);
-        VegaGlyph *glyphs = malloc(shown * sizeof(*glyphs));
+        VegaGlyph *cells = malloc((shown ? shown : 1) * sizeof(*cells));
         for (uint32_t i = 0; i < shown; ++i)
-            glyphs[i] = (VegaGlyph){(float)(i % (uint32_t)across * (uint32_t)cell) + 2 - (float)(left * scale),
-                                    (float)(i / (uint32_t)across * (uint32_t)cell) + 2 + (float)(top * scale), (float)scale, i};
-        GLuint texture = gl_target(width, height);
-        vega_gl_draw(&gl, &s, glyphs, shown, width, height);
-        unsigned char *drawn = malloc((size_t)width * (size_t)height), *reference = malloc((size_t)cell * (size_t)cell);
+            cells[i] = (VegaGlyph){(int)(i % (uint32_t)across) * cell_width + 1 - left,
+                                   (int)(i / (uint32_t)across) * cell_height + 1 + top, i};
+        GLuint grid = gl_target(grid_width, grid_height);
+        vega_gl_draw(&gl, &atlas, cells, shown, grid_width, grid_height);
+        unsigned char *drawn = malloc((size_t)grid_width * (size_t)grid_height), *reference = NULL;
+        size_t reference_size = 0;
         glPixelStorei(GL_PACK_ALIGNMENT, 1);
-        glReadPixels(0, 0, width, height, GL_RED, GL_UNSIGNED_BYTE, drawn);
+        glReadPixels(0, 0, grid_width, grid_height, GL_RED, GL_UNSIGNED_BYTE, drawn);
+        glDeleteTextures(1, &grid);
         double total = 0;
         uint64_t pixels = 0;
         int worst = 0;
+        size_t need = (size_t)cell_width * (size_t)cell_height;
+        if (!(reference = malloc(need))) fail("out of memory");
+        reference_size = need;
         for (uint32_t i = 0; i < shown; ++i) {
             VegaSymbol symbol = vega_symbol_at(&s, i);
-            int x0 = (int)(i % (uint32_t)across) * cell, y0 = (int)(i / (uint32_t)across) * cell;
-            vega_symbol_raster(&s, &symbol, scale, (x0 - glyphs[i].x) / scale, (glyphs[i].y - y0) / scale, cell, cell, reference);
-            for (int j = 0; j < cell; ++j)
-                for (int k = 0; k < cell; ++k) {
+            int x0 = (int)(i % (uint32_t)across) * cell_width, y0 = (int)(i / (uint32_t)across) * cell_height;
+            /* the cell's pixels, relative to the pen: from (x0 - pen x, pen y - y0) */
+            vega_symbol_raster(&s, &symbol, scale, (x0 - cells[i].x) / scale, (cells[i].y - y0) / scale, cell_width,
+                               cell_height, reference, (size_t)cell_width);
+            for (int j = 0; j < cell_height; ++j)
+                for (int k = 0; k < cell_width; ++k) {
                     /* GL's rows run bottom up */
-                    int gpu_value = drawn[(size_t)(height - 1 - (y0 + j)) * (size_t)width + (size_t)(x0 + k)];
-                    int difference = abs(gpu_value - reference[j * cell + k]);
+                    int gpu_value = drawn[(size_t)(grid_height - 1 - (y0 + j)) * (size_t)grid_width + (size_t)(x0 + k)];
+                    int difference = abs(gpu_value - reference[(size_t)j * (size_t)cell_width + (size_t)k]);
                     total += difference;
                     worst = difference > worst ? difference : worst;
                     ++pixels;
                 }
         }
-        /* Pages: glyphs at their advances in lines, a full target's worth. */
+        (void)reference_size;
+        free(cells);
+        /* Pages: glyphs at their advances in lines, a 1920x1080 target's worth. */
         int page_width = 1920, page_height = 1080;
         uint32_t page = 0, capacity = 1 << 16;
         VegaGlyph *text = malloc(capacity * sizeof(*text));
-        float x = 0, line = (float)size * 1.2f;
-        for (float y = line; y + (float)size * 0.3f < (float)page_height && page < capacity; ) {
-            VegaSymbol symbol = vega_symbol_at(&s, page % s.symbols);
-            float advance = (float)(symbol.advance * scale);
-            if (x + advance > (float)page_width) { x = 0; y += line; continue; }
-            text[page++] = (VegaGlyph){x, y, (float)scale, page % s.symbols};
+        int x = 0, line = (int)(size * 1.2);
+        for (int y = line; y + size / 3 < page_height && page < capacity;) {
+            uint32_t index = page % s.symbols;
+            int advance = (int)lroundf(vega_symbol_at(&s, index).advance * scale);
+            if (x + advance > page_width) { x = 0; y += line; continue; }
+            text[page++] = (VegaGlyph){x, y, index};
             x += advance;
         }
-        glDeleteTextures(1, &texture);
-        texture = gl_target(page_width, page_height);
-        vega_gl_draw(&gl, &s, text, page, page_width, page_height);
+        GLuint texture = gl_target(page_width, page_height);
+        vega_gl_draw(&gl, &atlas, text, page, page_width, page_height);
         glFinish();
         started = seconds();
-        for (int f = 0; f < frames; ++f) vega_gl_draw(&gl, &s, text, page, page_width, page_height);
+        for (int f = 0; f < frames; ++f) vega_gl_draw(&gl, &atlas, text, page, page_width, page_height);
         glFinish();
         double page_ms = (seconds() - started) / frames * 1e3;
-        printf("{\"pixels_per_em\":%ld,\"symbols\":%u,\"mean_difference\":%.3f,\"worst\":%d,\"page_glyphs\":%u,"
-               "\"page_ms\":%.2f,\"upload_ms\":%.1f}\n", size, shown, pixels ? total / (double)pixels : 0, worst, page, page_ms, open_ms);
+        printf("{\"pixels_per_em\":%ld,\"symbols\":%u,\"mean_difference\":%.3f,\"worst\":%d,\"atlas\":\"%dx%d\","
+               "\"fill_ms\":%.1f,\"page_glyphs\":%u,\"page_ms\":%.2f,\"upload_ms\":%.1f}\n", size, shown,
+               pixels ? total / (double)pixels : 0, worst, atlas.width, atlas.height, fill_ms, page, page_ms, open_ms);
         glDeleteTextures(1, &texture);
+        vega_gl_atlas_close(&atlas);
         free(text);
-        free(glyphs);
         free(drawn);
         free(reference);
     }
diff --git a/symbols/symbols.c b/symbols/symbols.c
index cedd68f..6f355f1 100644
--- a/symbols/symbols.c
+++ b/symbols/symbols.c
@@ -173,13 +173,13 @@ double vega_coverage(const VegaCurve *curves, uint32_t count, double scale, doub
 }
 
 void vega_symbol_raster(const VegaSymbols *s, const VegaSymbol *symbol, double scale, double left, double top, int width,
-                        int height, unsigned char *out) {
+                        int height, unsigned char *out, size_t stride) {
     VegaCurve *curves = malloc((symbol->count ? symbol->count : 1) * sizeof(*curves));
     uint32_t count = curves ? vega_symbol_curves(s, symbol, curves) : 0;
     for (int j = 0; j < height; ++j)
         for (int i = 0; i < width; ++i) {
             double coverage = vega_coverage(curves, count, scale, left + (i + 0.5) / scale, top - (j + 0.5) / scale);
-            out[(size_t)j * (size_t)width + (size_t)i] = (unsigned char)lround(coverage * 255);
+            out[(size_t)j * stride + (size_t)i] = (unsigned char)lround(coverage * 255);
         }
     free(curves);
 }
diff --git a/symbols/symbols.h b/symbols/symbols.h
index 998ff3f..20fbd22 100644
--- a/symbols/symbols.h
+++ b/symbols/symbols.h
@@ -47,8 +47,9 @@ uint32_t vega_symbol_curves(const VegaSymbols *, const VegaSymbol *, VegaCurve *
 /* The coverage (0 to 1) of the pixel whose center is (x, y) in units, at
  * scale pixels per unit: the exact area the curves' outline covers. */
 double vega_coverage(const VegaCurve *, uint32_t count, double scale, double x, double y);
-/* width x height 8-bit coverage, rows top first; pixel (i, j) has its center
- * at (left + (i + 0.5) / scale, top - (j + 0.5) / scale) in units. */
+/* width x height 8-bit coverage, rows top first, stride bytes apart; pixel
+ * (i, j) has its center at (left + (i + 0.5) / scale, top - (j + 0.5) / scale)
+ * in units. */
 void vega_symbol_raster(const VegaSymbols *, const VegaSymbol *, double scale, double left, double top, int width,
-                        int height, unsigned char *out);
+                        int height, unsigned char *out, size_t stride);
 #endif