git.haldean.org plotter / 48ddd1e
lots of pntlsm stuff haldean a month ago
5 changed file(s) with 296 addition(s) and 20 deletion(s). Raw diff Collapse all Expand all
+6
-0
Makefile less more
0
1 pntlsm.exe: pntlsm.cpp lib/array2d.hpp lib/pixels.hpp lib/hsvrgb.cpp
2 g++ -ggdb -DRED -DSCATTER -O3 -fopenmp -ffast-math -Wall -std=gnu++17 $< -o $@
3
4 %.nc: art/%.png pntlsm.exe
5 ./pntlsm $< > $@
+129
-0
lib/array2d.hpp less more
0 #pragma once
1
2 #include <cstdio>
3 #include <cstdlib>
4 #include <cstring>
5 #include <iostream>
6 #include <stdexcept>
7 #include <string>
8 #include <vector>
9
10 #include "lodepng.h"
11
12 using namespace std;
13
14 template<typename T>
15 struct array2d {
16 T *values;
17 size_t w;
18 size_t h;
19
20 array2d() : array2d(0, 0) {};
21
22 array2d(const size_t w, const size_t h) noexcept : w(w), h(h) {
23 values = (T*) calloc(w * h, sizeof(T));
24 }
25
26 array2d(const array2d<T> &o) : w(o.w), h(o.h) {
27 values = (T*) calloc(w * h, sizeof(T));
28 memcpy(values, o.values, w * h * sizeof(T));
29 }
30
31 array2d(array2d<T> &&o) : w(o.w), h(o.h), values(o.values) {
32 o.values = nullptr;
33 }
34
35 ~array2d() noexcept {
36 if (values != nullptr) free(values);
37 }
38
39 array2d<T>& operator=(const array2d<T> &o) = delete;
40
41 array2d<T>& operator=(array2d<T> &&o) {
42 w = o.w;
43 h = o.h;
44 T *t = values;
45 values = o.values;
46 o.values = t;
47 return *this;
48 }
49
50 T& operator()(size_t i) noexcept {
51 return values[i];
52 }
53
54 const T& operator()(size_t i) const noexcept {
55 return values[i];
56 }
57
58 T& operator()(size_t i, size_t j) {
59 if (__builtin_expect(i > w or j > h, 0)) {
60 char *msg;
61 asprintf(&msg, "bad access in array2d: i=%llu j=%llu, w=%llu h=%llu",
62 i, j, w, h);
63 throw out_of_range(msg);
64 }
65 return values[w * j + i];
66 }
67
68 const T& operator()(size_t i, size_t j) const {
69 if (__builtin_expect(i > w or j > h, 0)) {
70 char *msg;
71 asprintf(&msg, "bad access in array2d: i=%llu j=%llu, w=%llu h=%llu",
72 i, j, w, h);
73 throw out_of_range(msg);
74 }
75 return values[w * j + i];
76 }
77
78 T* data() noexcept {
79 return values;
80 }
81
82 const T* data() const noexcept {
83 return values;
84 }
85
86 void fill(T val) {
87 const size_t n = w * h;
88 for (size_t i = 0; i < n; i++) {
89 values[i] = val;
90 }
91 }
92
93 array2d<T> transposed() const {
94 array2d<T> res(h, w);
95 for (size_t i = 0; i < w; i++) {
96 for (size_t j = 0; j < h; j++) {
97 res(j, i) = operator()(i, j);
98 }
99 }
100 return res;
101 }
102 };
103
104 bool to_png(string fname, const array2d<double> &grays) {
105 vector<uint8_t> normed;
106 const size_t n = grays.w * grays.h;
107 normed.reserve(n);
108 for (size_t i = 0; i < n; i++) {
109 double clamped = grays(i);
110 clamped = clamped > 1 ? 1 : clamped < 0 ? 0 : clamped;
111 normed.push_back((uint8_t) (clamped * 255));
112 }
113 auto err = lodepng::encode(fname, normed.data(), grays.w, grays.h, LCT_GREY);
114 if (err) {
115 cout << "failed to write png: " << err << endl;
116 return false;
117 }
118 return true;
119 }
120
121 bool to_png(string fname, const array2d<uint8_t> &vals) {
122 auto err = lodepng::encode(fname, vals.data(), vals.w, vals.h, LCT_GREY);
123 if (err) {
124 cout << "failed to write png: " << err << endl;
125 return false;
126 }
127 return true;
128 }
+66
-0
lib/hsvrgb.cpp less more
0 #include <cmath>
1
2 void
3 hsv2rgb(double h, double s, double v, double &r, double &g, double &b)
4 {
5 double hh, p, q, t, ff;
6 long i;
7
8 if(s <= 0.0) { // < is bogus, just shuts up warnings
9 r = v;
10 g = v;
11 b = v;
12 return;
13 }
14 hh = h;
15 if(hh >= 360.0) hh = 0.0;
16 hh /= 60.0;
17 i = (long)hh;
18 ff = hh - i;
19 p = v * (1.0 - s);
20 q = v * (1.0 - (s * ff));
21 t = v * (1.0 - (s * (1.0 - ff)));
22
23 switch(i) {
24 case 0:
25 r = v;
26 g = t;
27 b = p;
28 break;
29 case 1:
30 r = q;
31 g = v;
32 b = p;
33 break;
34 case 2:
35 r = p;
36 g = v;
37 b = t;
38 break;
39
40 case 3:
41 r = p;
42 g = q;
43 b = v;
44 break;
45 case 4:
46 r = t;
47 g = p;
48 b = v;
49 break;
50 case 5:
51 default:
52 r = v;
53 g = p;
54 b = q;
55 break;
56 }
57 }
58
59 void rainbow(double t, uint8_t &r, uint8_t &g, uint8_t &b) {
60 double rd, gd, bd;
61 hsv2rgb(360.0 * t, 1, 1, rd, gd, bd);
62 r = floor(255 * rd);
63 g = floor(255 * gd);
64 b = floor(255 * bd);
65 }
+79
-0
lib/pixels.hpp less more
0 #pragma once
1
2 #include "array2d.hpp"
3 #include "lodepng.h"
4
5 #include <vector>
6
7 using namespace std;
8
9 template<typename T>
10 struct pixels {
11 array2d<T> r;
12 array2d<T> g;
13 array2d<T> b;
14 size_t n;
15 size_t w;
16 size_t h;
17
18 pixels(size_t w, size_t h) noexcept
19 : r(w, h)
20 , g(w, h)
21 , b(w, h)
22 , n(w * h)
23 , w(w)
24 , h(h)
25 {}
26
27 pixels() : pixels(0, 0) {};
28 pixels(const pixels&) = delete;
29 pixels(pixels&&) = default;
30 pixels<T>& operator=(pixels<T> &&o) = default;
31
32 operator bool() const {
33 return n > 0;
34 }
35 };
36
37 pair<pixels<uint8_t>*, bool> from_png(char *fname) {
38 // 4 components per pixel, rgba
39 vector<uint8_t> rawdata;
40 uint32_t width, height;
41
42 auto err = lodepng::decode(rawdata, width, height, fname);
43 if (err) {
44 cout << "couldn't load png: " << err << endl;
45 return make_pair(nullptr, false);
46 }
47 if (rawdata.size() % 4 != 0) {
48 cout << "bad size for rawdata " << rawdata.size() << endl;
49 return make_pair(nullptr, false);
50 }
51
52 pixels<uint8_t> * const inpx = new pixels<uint8_t>(width, height);
53
54 #pragma omp parallel for
55 for (size_t i = 0; i < inpx->n; i++) {
56 inpx->r(i) = rawdata[4 * i + 0];
57 inpx->g(i) = rawdata[4 * i + 1];
58 inpx->b(i) = rawdata[4 * i + 2];
59 }
60
61 return make_pair(inpx, true);
62 }
63
64 bool to_png(string fname, const pixels<uint8_t> &img) {
65 uint8_t *data = (uint8_t*) malloc(3 * img.n);
66 for (size_t i = 0; i < img.n; i++) {
67 data[3 * i + 0] = img.r(i);
68 data[3 * i + 1] = img.g(i);
69 data[3 * i + 2] = img.b(i);
70 }
71
72 auto err = lodepng::encode(fname, data, img.w, img.h, LCT_RGB);
73 if (err) {
74 cout << "couldn't save pixel image: " << err << endl;
75 return false;
76 }
77 return true;
78 }
+16
-20
pntlsm.cpp less more
0 #if 0
1 set -eux
2 #g++ -ggdb -fsanitize=address \
3 g++ -ggdb \
4 -D$1 -D$2 -O3 -fopenmp -ffast-math -Wall -std=c++17 pntlsm.cpp -o pntlsm
5 time ./pntlsm $3
6 exit
7 #else
0 #define _GNU_SOURCE
81
92 #include "lib/array2d.hpp"
103 #include "lib/pixels.hpp"
2215
2316 using namespace std;
2417
25 constexpr double pexp = 1.0;
26 constexpr double ppmm = 72.0 / 25.4; // pixels per mm
27 constexpr size_t npoints = 65000;
28 constexpr double bucketsize_mm = 10;
18 constexpr double pexp = 1.8;
19 constexpr double ppmm = 150.0 / 25.4; // pixels per mm
20 constexpr double pointfill = 0.35; // number of output points per input pixel
21 constexpr double bucketsize_mm = 3;
2922 constexpr int seed = 0x54481;
3023
3124 /* scale used for sample and sortimage */
32 constexpr double sscale = 2;
25 constexpr double sscale = 1;
3326 constexpr double ppmm_sortimage = ppmm * sscale;
3427
3528 template<typename T>
112105 }
113106
114107 void scatter_sample() {
115 double xmax = image.w / ppmm;
116 double ymax = image.h / ppmm;
108 const double xmax = image.w / ppmm;
109 const double ymax = image.h / ppmm;
110 const size_t npoints = pointfill * image.n;
117111
118112 std::mt19937 gen(seed);
119113 std::uniform_real_distribution<double> xdist(0.0, xmax);
142136 int ylo = floor(ypx);
143137 int yhi = ceil(ypx);
144138
139 if (yhi >= (int64_t) probability.h) yhi = ylo;
140 if (xhi >= (int64_t) probability.h) xhi = xlo;
141
145142 double s11 = probability(xlo, ylo);
146143 double s12 = probability(xlo, yhi);
147144 double s22 = probability(xhi, yhi);
150147 double y1 = (xhi - xpx) * s11 + (xpx - xlo) * s21;
151148 double y2 = (xhi - xpx) * s12 + (xpx - xlo) * s22;
152149 double prob = (yhi - ypx) * y1 + (ypx - ylo) * y2;
150 prob = pow(prob, pexp);
153151
154152 if (pdist(gen) < prob) {
155153 points.push_back(make_pair(x, y));
181179 }
182180
183181 void linefeeder() {
184 ofstream nc("_pntlsm.nc");
182 ostream &nc = cout;
185183 nc << setprecision(3)
186184 << fixed
187185 << "G0 Z1\n"
279277 for (const auto &p : bucket) {
280278 double xpx = p.first * ppmm_sortimage;
281279 double ypx = p.second * ppmm_sortimage;
282 double t = (double) pi++ / (double) npoints;
280 double t = (double) pi++ / (double) points.size();
283281 uint8_t r, g, b;
284282 rainbow(t, r, g, b);
285283 sortimage.r(xpx, ypx) = r;
288286 }
289287 }
290288
291 ofstream nc("_pntlsm.nc");
289 ostream &nc = cout;
292290 nc << setprecision(3) << fixed
293291 << "G0 Z1\nG0 X0 Y0\n";
294292 for (const auto &bidx : sorted_bucket_indexes) {
356354 if (not saveimages()) return EXIT_FAILURE;
357355 return EXIT_SUCCESS;
358356 }
359
360 #endif