Sepia

Sepia is a zero-dependency C++20 plotting framework: Xiaolin Wu anti-aliased lines, a built-in bitmap font, cache-aligned buffers and LTTB decimation render publication-quality 2D plots of a billion points in under a second.

Sepia works as a single header, no dependencies. Constraint behind every decision: render time follows output resolution, not input size.

01Pipeline

1. Ingest: copy into AlignedBuffer, or a zero-copy strided view over caller-owned data.
2. LTTB decimation (optional): N points collapse to target buckets, largest-triangle selection per bucket.
3. Coordinate transform: data space to pixel space, two FLOPs per point, precomputed scale.
4. Rasterize: Wu AA lines for series, Bresenham for axis/legend chrome.
5. PPM encode: P6 header plus raw RGB, no compression.

Memory

AlignedBuffer<T>: contiguous, 64-byte aligned (cache line, AVX-512 load), move-only over posix_memalign / _aligned_malloc. Copy constructor deleted: a deep copy of a megapoint buffer is the costliest user mistake.

Three ownership tiers:

A std::variant<Series, ExternalSeries> per entry dispatches at container level, not per point. ExternalSeries wraps DataView { ptr, count, stride }: strided reads into an existing SoA layout are free.

Per-frame scratch uses a bump Arena, one aligned slab, reset() per frame. Zero heap traffic in the render loop.

LTTB decimation

The algorithmic core. At N >> W thousands of segments collapse into one pixel column; every k-th point destroys spikes, outliers, envelope shape.

Largest-Triangle-Three-Buckets (Steinarsson): partition into target buckets; keep from each the point forming the largest triangle with the last kept point and the average of the next bucket. Area proxies visual information, so extrema survive.

// Inner loop: pick the point in this bucket with max triangle area
// against the last kept point (ax, ay) and next bucket's centroid.
f64 max_area = -1.0; usize max_idx = bucket_start;
f64 ax = xv[a], ay = yv[a];
for (usize j = bucket_start; j < bucket_end; ++j) {
    f64 area = std::abs(
        (ax - avg_x) * (yv[j] - ay) - (ax - xv[j]) * (avg_y - ay)
    ) * 0.5;   // 2D cross product: no sqrt, no trig, branchless-friendly
    if (area > max_area) { max_area = area; max_idx = j; }
}

Four multiplies, three subtracts, one abs. Autovectorizes under -O3 -mavx2 -mfma.

Lazy: runs inside render() on views, never mutating caller data, so re-rendering at another LOD target is free. Cost: one O(N) streaming pass plus O(target) rasterization.

Rasterization

Xiaolin Wu AA lines for series, Bresenham for axis and legend chrome. Wu splits column coverage across the two straddling pixels by fractional intercept: two set_pixel calls with computed alpha, ~2x Bresenham against 4-16x MSAA. Compositing stays integer, (src * a + dst * (255 - a)) / 255.

void clear(Color c) {
    u32* px = reinterpret_cast<u32*>(pixels_.data()); // buffer is 64B-aligned
    std::fill(px, px + usize(width_) * height_, pack(c)); // -> vectorized stores
}

Text without FreeType

FreeType would be the only dependency, and a heavy one, for 1% of the pixels. Glyphs are 5x7 bitmaps, one 5-bit row mask per u8.

struct Glyph { u8 rows[7]; };
// 'A' = 0x0E,0x11,0x11,0x1F,0x11,0x11,0x11  (each hex byte is one pixel row)
if (g.rows[row] & (0x10 >> col)) canvas.set_pixel(...);

35 bit tests per glyph. Integer scaling gives font sizes.

Rest

• Ticks: nice-numbers loop, normalize to [1,10), snap to 1/2/5/10. Decade path for log.
• Transform: scale precomputed per render; to_px_x/y is two inlined FLOPs. safe_log10 blocks -inf.
• Output: PPM (P6), 15-byte header plus raw RGB. convert out.ppm out.png.
• API: plot() returns a PlotCommand committing in its destructor. Scope exit is finalize, so no missing-finalize bug exists.

02Trade-offs

Through-line: optimize the O(N) stage, keep the O(pixels) stage simple.

03Benchmarks

AMD Ryzen 5 7600X, single-threaded, full render(), stress/stresstest.cpp. g++ -std=c++20 -O3 -mavx2 -mfma -ffast-math -fno-trapping-math -fno-math-errno -march=native

No-LTTB is linear from ~10⁵ (log-log slope 1). LTTB stays fixed-cost-dominated until ~10⁷, then same slope, far smaller constant: ~0.8 ns/point at 10⁹. Crossover is near-immediate, 4.8x at 50k points.

04Sharp edges

LTTB bucket boundaries. Index math (i-1) * bucket_size + 1 with float bucket_size needs three clamps: pathological N/target ratios yield an empty bucket, or push bucket_end past n-1. Failure mode is one duplicated or skipped point, visible only against a reference render.

-ffast-math vs. sentinels. BBox seeds bounds with +max/-max over strict </> chains, and the flag licenses no-NaN assumptions plus FP reassociation: "Debug fine, Release axis limits wrong." So BBox::empty() tests x_min > x_max, never the sentinels.

05Next

The constraint is how many points reach the rasterizer, not rasterizer speed. Highest leverage was an algorithm (LTTB), not a SIMD kernel: the 64-byte aligned, contiguous AlignedBuffer<T> layout is what lets that algorithm stream at nanoseconds per point, one cache line at a time, no pointer chasing.

Decimation is one independent pass per series: thread it, or SIMD the bucket scan. At ~1 ms per million points, re-running LTTB per zoom window makes a live viewport over a memory-mapped dataset viable.

06Further reading

• Downsampling Time Series for Visual Representation: Steinarsson's LTTB thesis.
• An Efficient Antialiasing Technique: Wu, SIGGRAPH '91.
• Sepia on GitHub: header, examples, stress harness.