plot
Header-only C++ SVG charts — heatmap, line, scatter — GR-style API + Solarized themes
Loading...
Searching...
No Matches
contour.hpp
1/* Copyright (c) 2026 Steven Varga, Toronto, ON, Canada
2 * MIT License — see LICENSE
3 *
4 * plot::contour(grid, opts...) — a deferred view (plot::view) drawing iso-contour
5 * lines of a row-major plot::mat<double> field via the marching-squares
6 * algorithm. plot::levels{n} sets the number of iso-levels (default 8); each
7 * level is coloured by its fraction along theme.gradient. Every cell that the
8 * iso-line crosses contributes one short <polyline> segment (linear interpolation
9 * on the crossed edges).
10 *
11 * Built on the #5 view pattern + impl::render_panel_into. Dependency-free.
12 */
13#ifndef PLOT_CONTOUR_HPP
14#define PLOT_CONTOUR_HPP
15
16#include <string>
17#include <vector>
18#include <tuple>
19#include <utility>
20#include <algorithm>
21#include <cmath>
22#include <cstddef>
23#include <cstdint>
24#include <limits>
25#include <type_traits>
26
27#include "tags.hpp"
28#include "meta.hpp"
29#include "attributes.hpp"
30#include "canvas.hpp"
31#include "theme.hpp"
32#include "gr.hpp"
33#include "heatmap.hpp" // plot::mat<T>
34#include "view.hpp"
35
36namespace plot {
37 // number of iso-levels for plot::contour (default 8).
38 struct levels { using value_type = tag::levels_t; std::size_t value; };
39}
40
41namespace plot::impl {
42 template <class... opt_t>
43 std::size_t levels_of(const opt_t&... opts){
44 using lv_t = typename arg::tpos<tag::levels_t, opt_t...>;
45 if constexpr( lv_t::present ){
46 auto tuple = std::forward_as_tuple(opts...);
47 std::size_t l = std::get<lv_t::position>(tuple).value;
48 return l ? l : 1;
49 }
50 return 8;
51 }
52
53 template <class... opt_t>
54 struct contour_view {
55 using value_type = tag::view_t;
56 std::vector<double> data; // owned copy, row-major
57 std::size_t rows = 0, cols = 0;
58 std::tuple<opt_t...> opts;
59
60 std::pair<std::size_t,std::size_t> natural() const {
61 return std::apply([](const auto&... o){
62 return impl::natural_size(std::size_t{640}, std::size_t{400}, o...); }, opts);
63 }
64 double at(std::size_t i, std::size_t j) const { return data[i*cols + j]; }
65
66 void draw_into(canvas_t& cv, float x, float y, float w, float h) const {
67 std::apply([&](const auto&... o){
68 const theme_t& th = impl::resolve_theme(o...);
69 auto [title, xl, yl] = impl::texts(o...);
70 const std::size_t nlev = impl::levels_of(o...);
71
72 double lo = std::numeric_limits<double>::infinity();
73 double hi = -std::numeric_limits<double>::infinity();
74 for(double d : data){ if(d<lo) lo=d; if(d>hi) hi=d; }
75 if( !(lo<=hi) ){ lo=0; hi=1; }
76 if( hi <= lo ) hi = lo + 1.0;
77
78 impl::scale_t sx = impl::make_scale(0.0, double(cols>0?cols-1:0), false, 10.0);
79 impl::scale_t sy = impl::make_scale(0.0, double(rows>0?rows-1:0), false, 10.0);
80 impl::render_panel_into(cv, x, y,
81 static_cast<std::size_t>(w), static_cast<std::size_t>(h),
82 th, sx, sy, title, xl, yl, {},
83 [&](impl::canvas_t& c, auto px, auto py){
84 using attribute_t = plot::attribute::element_t;
85 if( rows < 2 || cols < 2 ) return;
86 // linear interpolation of the crossing point between two
87 // grid corners (a,b) at field values (va,vb) for level L.
88 auto lerp = [](double a, double b, double va, double vb, double L){
89 double d = (vb - va);
90 double t = (std::abs(d) < 1e-12) ? 0.5 : (L - va)/d;
91 if( t < 0 ) t = 0; else if( t > 1 ) t = 1;
92 return a + t*(b - a);
93 };
94 for(std::size_t k=0;k<nlev;++k){
95 double frac = (nlev==1) ? 0.5 : double(k)/double(nlev-1);
96 double L = lo + frac*(hi-lo);
97 std::uint32_t col = impl::gradient3(th.gradient, frac);
98 attribute_t a; a.color = plot::attribute::color_t{ col };
99 a.stroke = plot::attribute::stroke_t{1.0f, 1.4f, {}, {}, {}};
100 for(std::size_t i=0;i+1<rows;++i)
101 for(std::size_t j=0;j+1<cols;++j){
102 // corner values: TL,TR,BR,BL of cell (i,j)
103 double v00 = at(i, j), v01 = at(i, j+1);
104 double v11 = at(i+1, j+1), v10 = at(i+1, j);
105 // marching-squares case index.
106 int cs = (v00>L?1:0) | (v01>L?2:0) | (v11>L?4:0) | (v10>L?8:0);
107 if( cs==0 || cs==15 ) continue;
108 // edge crossing points in grid coords (col=x, row=y):
109 // top edge between (i,j)-(i,j+1), right (i,j+1)-(i+1,j+1),
110 // bottom (i+1,j)-(i+1,j+1), left (i,j)-(i+1,j).
111 double tx = lerp(double(j), double(j+1), v00, v01, L), ty = double(i);
112 double rx = double(j+1), ry = lerp(double(i), double(i+1), v01, v11, L);
113 double bx = lerp(double(j), double(j+1), v10, v11, L), by = double(i+1);
114 double lx = double(j), ly = lerp(double(i), double(i+1), v00, v10, L);
115 auto seg = [&](double ax,double ay,double bx2,double by2){
116 std::vector<float> X{ px(ax), px(bx2) };
117 std::vector<float> Y{ py(ay), py(by2) };
118 c.poly_line(X, Y, a);
119 };
120 switch(cs){
121 case 1: case 14: seg(lx,ly, tx,ty); break;
122 case 2: case 13: seg(tx,ty, rx,ry); break;
123 case 3: case 12: seg(lx,ly, rx,ry); break;
124 case 4: case 11: seg(rx,ry, bx,by); break;
125 case 6: case 9: seg(tx,ty, bx,by); break;
126 case 7: case 8: seg(lx,ly, bx,by); break;
127 case 5: seg(lx,ly, tx,ty); seg(rx,ry, bx,by); break;
128 case 10: seg(tx,ty, rx,ry); seg(lx,ly, bx,by); break;
129 default: break;
130 }
131 }
132 }
133 });
134 }, opts);
135 }
136 };
137}
138
139namespace plot {
140 template <class T, class... opt_t,
141 class = std::enable_if_t<std::is_arithmetic_v<T>>>
142 impl::contour_view<opt_t...> contour(const mat<T>& grid, opt_t... opts){
143 impl::contour_view<opt_t...> v{ {}, grid.rows, grid.cols, std::make_tuple(opts...) };
144 v.data.assign(grid.begin(), grid.end());
145 return v;
146 }
147}
148#endif