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. No FreeType, no libpng, no graphics stack.

How do you look at a hundred million samples, not summarize them? Every plotting stack I reached for from C++ had the same shape: serialize the data out of the process, hand it to Python or gnuplot or a browser, and wait. Serialization alone cost more than the analysis, and past a few million points the plotting layer became the bottleneck of the whole pipeline. The data was already sitting in memory, contiguous, hot in cache, and I was paying an IPC round-trip and a scripting-language rendering loop to draw lines through it.

The real question: what does a plotter look like if it's designed from the pixel buffer up to live inside a C++ process, with the same performance discipline as the code producing the data?

Existing native options either drag in a dependency tree (Cairo, FreeType, libpng, a GUI toolkit) or give up on rendering quality. So I built Sepia: a single-header, zero-dependency C++20 plotting framework for publication-quality 2D output. Everything (rasterization, anti-aliasing, font rendering, axis layout, image encoding) is in one sepia.hpp. Copy the header in and you're done. The name comes from Sepia officinalis, the cuttlefish: it releases ink, Sepia inks pixels.

The design constraint behind every decision below: render time should be a function of output resolution, not input size. A 700x450 plot has a fixed number of pixels; a billion input points shouldn't take 30,000x longer to draw than thirty thousand.

01How it works

The pipeline is straight-line: data ingestion -> (optional) LTTB decimation -> coordinate transform -> rasterization -> PPM encode.

Memory: cache-line aligned, move-only, ownership made explicit

All owned data lives in AlignedBuffer<T>, a move-only contiguous buffer aligned to 64 bytes: one cache line, and the natural alignment for AVX-512 loads should the hot loops get vectorized further.

template <typename T, usize Alignment = 64>
class AlignedBuffer {
    // posix_memalign / _aligned_malloc under the hood.
    // Move-only: copying megapoint buffers by accident should be
    // a compile error, not a profiler discovery.
    AlignedBuffer(const AlignedBuffer&) = delete;
    AlignedBuffer(AlignedBuffer&& o) noexcept
        : size_(o.size_), data_(o.data_) { o.data_ = nullptr; o.size_ = 0; }
    ...
};

Deleting the copy constructor matters: in a library whose point is huge datasets, an accidental deep copy is the most expensive mistake a user can make, so it doesn't compile.

Ownership is a three-tier API decision:

Internally these are a std::variant<Series, ExternalSeries> per plot entry: one dispatch at the container level rather than a virtual call per point. ExternalSeries wraps a DataView { ptr, count, stride }, so strided access into an existing struct-of-arrays layout costs nothing.

For per-frame scratch work (tick labels, layout rects) there's a small bump Arena allocator: one 64-byte-aligned slab, pointer-bump alloc, reset() between frames. No general-purpose heap traffic in the render loop.

LTTB decimation: the reason a billion points is fine

This is the algorithmic core. Rendering N points onto W pixels of plot width is wasted work when N >> W: thousands of segments collapse into the same pixel column. The naive fix (take every k-th point) destroys exactly the features you plot to see: spikes, outliers, envelope shape.

Sepia uses Largest-Triangle-Three-Buckets decimation, from Sveinn Steinarsson's thesis "Downsampling Time Series for Visual Representation": partition the series into target buckets; from each, keep the point forming the largest triangle with the previously selected point and the average of the next bucket. Maximizing triangle area is a proxy for retained visual information: it preferentially keeps extrema and inflections, what your eye reads.

// 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; }
}

Two properties worth noting. The area is a raw 2D cross product, four multiplies, three subtracts, one abs, no sqrt or division, so the compiler autovectorizes the bucket scan cleanly under -O3 -mavx2 -mfma. And the pass is lazy: it runs inside render() on views, never mutating the caller's data, so you can re-render at a different LOD target without touching source arrays.

The payoff: render cost O(N) with a tiny constant (one sequential streaming pass) followed by O(target) rasterization, instead of O(N) rasterization with all its per-segment fixed costs.

Rasterization: Xiaolin Wu, chosen deliberately

Sepia implements Xiaolin Wu's anti-aliased line algorithm for data series, with plain Bresenham as the fast path for axis/legend chrome where AA buys nothing:

void draw_line(f64 x0, f64 y0, f64 x1, f64 y1, Color c, f64 width, bool aa = true) {
    if (aa) draw_line_aa(x0, y0, x1, y1, c, width);   // Wu: data series
    else    draw_line_bresenham(x0, y0, x1, y1, c, width); // chrome: 1px, no blend
}

Wu's algorithm rides the true line and, per column, splits coverage between the two straddling pixels using the fractional part of the intercept: two set_pixel calls with computed alpha instead of supersampling. Analytic coverage instead of multi-sample rendering, at roughly 2x Bresenham's cost rather than 4 to 16x for MSAA.

Alpha compositing stays in integers, (src * a + dst * (255 - a)) / 255, no float conversion per pixel. The full-canvas clear() packs RGBA into a u32 and hands it to std::fill over a reinterpreted pixel buffer, which the compiler vectorizes:

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: a 5x7 font in bitmasks

Pulling in FreeType for axis labels would make it the project's only dependency, and a heavy one, for 1% of the pixels. Instead, glyphs are 5x7 bitmaps, each row a 5-bit mask in a 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(...);

Rendering a glyph is 35 bit tests. Integer scaling gives "font sizes." Not typography; legible, deterministic, free.

The rest of the machinery

• Tick generation uses the classic nice-numbers loop (normalize the rough step to [1,10), snap to 1/2/5/10), with a separate decade-based path for log scales.
• Coordinate transform (CoordTransform) precomputes scale factors once per render; the per-point to_px_x/y is two FLOPs and inlined. Log scales clamp through safe_log10 to avoid -inf poisoning the transform.
• Output is PPM (P6): a 15-byte header and raw RGB, no compression library. convert out.ppm out.png covers the last mile.
• The fluent API commits via RAII: figure.plot(...) returns a PlotCommand whose destructor commits the fully-configured entry. The builder chain (.color(...).width(...).label(...)) mutates a local entry; there's no "did you forget to call finalize()" failure mode, because scope exit is finalize.

02Trade-offs

The through-line: optimize the O(N) stage (decimation, memory layout, streaming access), keep the O(pixels) stage simple and correct.

03Benchmarks

Methodology: AMD Ryzen 5 7600X, g++ -std=c++20 -O3 -mavx2 -mfma -ffast-math -fno-trapping-math -fno-math-errno -march=native, single-threaded, timing the full render() (decimate + transform + rasterize) per dataset size. Numbers are from the repo's reproducible stress harness (stress/stresstest.cpp); the results plot is itself rendered by Sepia.

The no-LTTB curve is linear in N from ~10⁵ on (log-log slope 1), carrying per-segment transform, Wu coverage math, and blended pixel writes for every input point. The LTTB curve stays dominated by its fixed costs until ~10⁷, then settles onto the same slope with a far smaller constant: the cost of one streaming pass, around 0.8 ns/point at 10⁹. The crossover is nearly immediate; LTTB already wins 4.8x at 50k points, and below that the decimation pass costs nothing measurable. A billion points to a finished image in 823 ms, in-process, no serialization.

04Where it went sideways

LTTB bucket boundaries. The bucket index math ((i-1) * bucket_size + 1 with float bucket_size) has three separate clamp sites, each earned: for pathological N/target ratios float truncation can produce an empty bucket or run bucket_end past n-1. The failure mode isn't a crash, it's one silently duplicated or skipped point, noticeable only as a subtle kink when diffing against a reference render.

The one that almost bit: -ffast-math vs. sentinel comparisons. -ffast-math licenses the compiler to assume no NaNs and reassociate floating-point ops. The BBox bounds computation initializes to +max/-max sentinels and relies on strict </> chains; under aggressive math flags, code like this is exactly where "works in Debug, wrong axis limits in Release" bugs breed. That's why BBox::empty() tests x_min > x_max rather than comparing against the sentinel constants, and why sentinel logic stays out of code paths compiled with value-changing assumptions.

05What I'd do better next time

I started from the wrong mental model: "rasterization is the expensive part, optimize the drawing." The actual constraint is how many points survive to the rasterizer. The single highest-leverage decision is an algorithm choice (LTTB), not a SIMD kernel; the memory-layout work (alignment, move-only buffers, zero-copy views) is what lets that algorithm stream at nanoseconds per point.

That points directly at the next version: decimation is one independent streaming pass per series, so it parallelizes trivially, a thread per series or SIMD-friendly bucket scans within one. And because render cost is now decoupled from input size, interactive use stops being absurd: re-running LTTB per zoom window at ~1 ms per million points is fast enough for a live viewport over a memory-mapped dataset.

06Further reading

• Downsampling Time Series for Visual Representation: Sveinn Steinarsson's thesis introducing LTTB; the bucket/triangle derivation and comparisons against alternative decimators.
• An Efficient Antialiasing Technique: Xiaolin Wu, SIGGRAPH '91; analytic-coverage AA lines.
• Sepia on GitHub: single header, examples, and the stress harness used for every number in this post.