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: look at 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. The serialization alone cost more than the analysis did, and the moment the dataset crossed a few million points, the plotting layer became the bottleneck of the entire pipeline. That felt backwards. 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, then: 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 implemented in one sepia.hpp. You copy the header into your project and you're done. The name comes from Sepia officinalis, the cuttlefish: it releases ink, Sepia inks pixels.

The design constraint that shaped every decision below: render time should be a function of output resolution, not input size. A 700x450 plot has a fixed number of pixels; there is no reason a billion input points should 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. Each stage was chosen or written with a specific constraint in mind.

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 ever 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 is the important line. In a library whose whole point is huge datasets, an accidental deep copy is the single most expensive mistake a user can make, so it doesn't compile.

Ownership itself is a three-tier API decision, not an implementation detail:

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 data to see: spikes, outliers, envelope shape.

Sepia uses Largest-Triangle-Three-Buckets decimation, from Sveinn Steinarsson's thesis "Downsampling Time Series for Visual Representation". The mechanism in plain language: partition the series into target buckets; from each bucket, keep the single point that forms the largest triangle with the previously selected point and the average of the next bucket. Maximizing triangle area is a proxy for maximizing retained visual information: it preferentially keeps extrema and inflections, which is 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 square roots, no divisions inside the scan, 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; you can re-render at a different LOD target without touching the source arrays.

The payoff is a render cost that's O(N) with a very small constant (one sequential streaming pass, the prefetcher's favorite access pattern) followed by O(target) rasterization, instead of O(N) rasterization with all its per-segment fixed costs.

Rasterization: Xiaolin Wu, chosen deliberately

Line quality is where "homemade plotter" usually shows. 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. That's the whole trick: analytic coverage instead of multi-sample rendering, at roughly 2x Bresenham's cost rather than 4 to 16x for MSAA-style approaches.

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

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 have made 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, the whole printable-ASCII table built once into a static array:

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". It's not typography; it's legible, deterministic, and 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, which is what makes the log-log stress-test plot below self-hosting.
• 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. Encoding is a straight buffer walk with no compression library and no color-management surprises. 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 to the figure. 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 below are from the repo's reproducible stress harness (stress/stresstest.cpp); the results plot is itself rendered by Sepia, so the framework benchmarks and plots itself.

Read the shape, not just the endpoint. 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 essentially immediate; in the harness, 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, and every one of them was earned: for pathological N/target ratios the float truncation can produce an empty bucket or run bucket_end past n-1, and the failure mode isn't a crash. It's one silently duplicated or skipped point that you only notice as a subtle kink when diffing against a reference render.

The one that almost bit: -ffast-math vs. sentinel comparisons. The benchmark flags include -ffast-math, which 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. Dodging that class of bug is why BBox::empty() tests x_min > x_max rather than comparing against the sentinel constants, and why sentinel and limit 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 in the codebase 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. Layout and algorithm did the heavy lifting; the pixel code just had to not waste it.

That reframing points directly at the next version: since decimation is one independent streaming pass per series, it parallelizes trivially, either a thread per series or SIMD-friendly bucket scans within one (the cross-product inner loop is already reduction-shaped). 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. The stress harness can check both claims before a line of that version exists.

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.