new sampling method that's not grid aligned
Haldean Brown
2 years ago
0 | 0 | #if 0 |
1 | 1 | set -eux |
2 | # g++ -ggdb -fsanitize=address \ | |
3 | g++ \ | |
4 | -D$1 -O3 -fopenmp -Wall -std=c++17 pntlsm.cpp -o pntlsm | |
5 | time ./pntlsm $2 | |
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 | 6 | exit |
7 | 7 | #else |
8 | 8 | |
9 | constexpr double pexp = 1.0; | |
10 | constexpr double ppmm = 72.0 / 25.4; // pixels per mm | |
11 | constexpr double left = 0; | |
12 | constexpr double top = 0; | |
13 | ||
14 | #include "lib/lodepng.h" | |
9 | #include "lib/array2d.hpp" | |
10 | #include "lib/pixels.hpp" | |
15 | 11 | #include "lib/lodepng.cpp" |
16 | ||
12 | #include "lib/hsvrgb.cpp" | |
13 | ||
14 | #include <algorithm> | |
17 | 15 | #include <fstream> |
18 | 16 | #include <iomanip> |
19 | 17 | #include <iostream> |
20 | 18 | #include <random> |
21 | 19 | #include <sstream> |
20 | #include <unordered_map> | |
22 | 21 | #include <vector> |
23 | 22 | |
24 | 23 | using namespace std; |
25 | 24 | |
26 | template<typename T> | |
27 | struct array2d { | |
28 | T *values; | |
29 | size_t w; | |
30 | size_t h; | |
31 | ||
32 | array2d() : array2d(0, 0) {}; | |
33 | ||
34 | array2d(const size_t w, const size_t h) noexcept : w(w), h(h) { | |
35 | values = (T*) calloc(w * h, sizeof(T)); | |
36 | } | |
37 | ||
38 | array2d(const array2d<T> &o) : w(o.w), h(o.h) { | |
39 | values = (T*) calloc(w * h, sizeof(T)); | |
40 | memcpy(values, o.values, w * h * sizeof(T)); | |
41 | } | |
42 | ||
43 | array2d(array2d<T> &&o) : w(o.w), h(o.h), values(o.values) { | |
44 | o.values = nullptr; | |
45 | } | |
46 | ||
47 | ~array2d() noexcept { | |
48 | if (values != nullptr) free(values); | |
49 | } | |
50 | ||
51 | array2d<T>& operator=(const array2d<T> &o) = delete; | |
52 | ||
53 | array2d<T>& operator=(array2d<T> &&o) { | |
54 | w = o.w; | |
55 | h = o.h; | |
56 | std::swap(values, o.values); | |
57 | return *this; | |
58 | } | |
59 | ||
60 | T& operator()(size_t i) noexcept { | |
61 | return values[i]; | |
62 | } | |
63 | ||
64 | const T& operator()(size_t i) const noexcept { | |
65 | return values[i]; | |
66 | } | |
67 | ||
68 | T& operator()(size_t i, size_t j) { | |
69 | if (__builtin_expect(i > w or j > h, 0)) { | |
70 | stringstream errmsg; | |
71 | errmsg << "bad access in array2d: i=" | |
72 | << i << " j=" << j << " w=" << w << " h=" << h; | |
73 | throw out_of_range(errmsg.str()); | |
74 | } | |
75 | return values[w * j + i]; | |
76 | } | |
77 | ||
78 | const T& operator()(size_t i, size_t j) const { | |
79 | if (__builtin_expect(i > w or j > h, 0)) { | |
80 | stringstream errmsg; | |
81 | errmsg << "bad access in array2d: i=" | |
82 | << i << " j=" << j << " w=" << w << " h=" << h; | |
83 | throw out_of_range(errmsg.str()); | |
84 | } | |
85 | return values[w * j + i]; | |
86 | } | |
87 | ||
88 | T* data() noexcept { | |
89 | return values; | |
90 | } | |
91 | ||
92 | const T* data() const noexcept { | |
93 | return values; | |
94 | } | |
95 | ||
96 | array2d<T> transposed() const { | |
97 | array2d<T> res(h, w); | |
98 | for (size_t i = 0; i < w; i++) { | |
99 | for (size_t j = 0; j < h; j++) { | |
100 | res(j, i) = operator()(i, j); | |
101 | } | |
102 | } | |
103 | return res; | |
104 | } | |
105 | ||
106 | }; | |
107 | ||
108 | bool to_png(string fname, const array2d<double> &grays) { | |
109 | vector<uint8_t> normed; | |
110 | const size_t n = grays.w * grays.h; | |
111 | normed.reserve(n); | |
112 | for (size_t i = 0; i < n; i++) { | |
113 | double clamped = grays(i); | |
114 | clamped = clamped > 1 ? 1 : clamped < 0 ? 0 : clamped; | |
115 | normed.push_back((uint8_t) (clamped * 255)); | |
116 | } | |
117 | auto err = lodepng::encode(fname, normed.data(), grays.w, grays.h, LCT_GREY); | |
118 | if (err) { | |
119 | cout << "failed to write png: " << err << endl; | |
120 | return false; | |
121 | } | |
122 | return true; | |
123 | } | |
124 | ||
125 | bool to_png(string fname, const array2d<uint8_t> &vals) { | |
126 | auto err = lodepng::encode(fname, vals.data(), vals.w, vals.h, LCT_GREY); | |
127 | if (err) { | |
128 | cout << "failed to write png: " << err << endl; | |
129 | return false; | |
130 | } | |
131 | return true; | |
132 | } | |
133 | ||
134 | template<typename T> | |
135 | struct pixels { | |
136 | array2d<T> r; | |
137 | array2d<T> g; | |
138 | array2d<T> b; | |
139 | size_t n; | |
140 | size_t w; | |
141 | size_t h; | |
142 | ||
143 | pixels(size_t w, size_t h) noexcept | |
144 | : r(w, h) | |
145 | , g(w, h) | |
146 | , b(w, h) | |
147 | , n(w * h) | |
148 | , w(w) | |
149 | , h(h) | |
150 | {} | |
151 | ||
152 | pixels() : pixels(0, 0) {}; | |
153 | pixels(const pixels&) = default; | |
154 | pixels(pixels&&) = default; | |
155 | pixels<T>& operator=(pixels<T> &&o) = default; | |
156 | }; | |
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; | |
29 | constexpr int seed = 0x54481; | |
30 | ||
31 | /* scale used for sample and sortimage */ | |
32 | constexpr double sscale = 2; | |
33 | constexpr double ppmm_sortimage = ppmm * sscale; | |
157 | 34 | |
158 | 35 | template<typename T> |
159 | 36 | constexpr double proportion(T v) { |
184 | 61 | } |
185 | 62 | |
186 | 63 | static pixels<uint8_t> image; |
64 | static pixels<uint8_t> sortimage; | |
187 | 65 | static array2d<double> probability; |
188 | 66 | static array2d<uint8_t> samples; |
67 | static vector<pair<double, double>> points; | |
189 | 68 | |
190 | 69 | bool load(char *fname) { |
191 | // 4 components per pixel, rgba | |
192 | vector<uint8_t> rawdata; | |
193 | uint32_t width, height; | |
194 | ||
195 | auto err = lodepng::decode(rawdata, width, height, fname); | |
196 | if (err) { | |
197 | cout << "couldn't load png: " << err << endl; | |
198 | return false; | |
199 | } | |
200 | if (rawdata.size() % 4 != 0) { | |
201 | cout << "bad size for rawdata " << rawdata.size() << endl; | |
202 | return false; | |
203 | } | |
204 | ||
205 | pixels<uint8_t> inpx(width, height); | |
206 | ||
207 | #pragma omp parallel for | |
208 | for (size_t i = 0; i < inpx.n; i++) { | |
209 | inpx.r(i) = rawdata[4 * i + 0]; | |
210 | inpx.g(i) = rawdata[4 * i + 1]; | |
211 | inpx.b(i) = rawdata[4 * i + 2]; | |
212 | } | |
213 | ||
214 | image = std::move(inpx); | |
215 | probability = array2d<double>(width, height); | |
216 | samples = array2d<uint8_t>(width, height); | |
70 | pixels<uint8_t> *loaded; | |
71 | bool ok; | |
72 | tie(loaded, ok) = from_png(fname); | |
73 | if (not ok) return false; | |
74 | image = std::move(*loaded); | |
75 | ||
76 | probability = array2d<double>(image.w, image.h); | |
77 | samples = array2d<uint8_t>(image.w, image.h); | |
217 | 78 | return true; |
218 | 79 | } |
219 | 80 | |
220 | void prob() { | |
81 | void genprob() { | |
221 | 82 | #pragma omp parallel for |
222 | 83 | for (size_t i = 0; i < probability.w; i++) { |
223 | 84 | for (size_t j = 0; j < probability.h; j++) { |
224 | 85 | probability(i, j) = P(image.r(i, j), image.g(i, j), image.b(i, j)); |
225 | } | |
226 | } | |
227 | } | |
228 | ||
229 | void sample() { | |
230 | std::random_device rd; | |
231 | std::mt19937 gen(rd()); | |
232 | std::uniform_real_distribution<double> pdist(0.0, 1.0); | |
233 | #pragma omp parallel for | |
234 | for (size_t i = 0; i < image.w; i++) { | |
235 | for (size_t j = 0; j < image.h; j++) { | |
236 | samples(i, j) = pdist(gen) < probability(i, j) ? 255 : 0; | |
237 | 86 | } |
238 | 87 | } |
239 | 88 | } |
262 | 111 | nc << "G1 Z1\n"; |
263 | 112 | } |
264 | 113 | |
265 | void savegcode() { | |
114 | void scatter_sample() { | |
115 | double xmax = image.w / ppmm; | |
116 | double ymax = image.h / ppmm; | |
117 | ||
118 | std::mt19937 gen(seed); | |
119 | std::uniform_real_distribution<double> xdist(0.0, xmax); | |
120 | std::uniform_real_distribution<double> ydist(0.0, ymax); | |
121 | std::uniform_real_distribution<double> pdist(0.0, 1.0); | |
122 | ||
123 | points.clear(); | |
124 | points.reserve(npoints); | |
125 | ||
126 | cout << "filling " << image.w << "x" << image.h | |
127 | << " (" << (image.w * image.h) << "px) with " | |
128 | << npoints << " points (" | |
129 | << ((npoints * 100.0) / ((double) image.w * image.h)) | |
130 | << "% fill)" | |
131 | << endl; | |
132 | ||
133 | size_t p = 0; | |
134 | while (p < npoints) { | |
135 | double x = xdist(gen); | |
136 | double y = ydist(gen); | |
137 | double xpx = x * ppmm; | |
138 | double ypx = y * ppmm; | |
139 | ||
140 | int xlo = floor(xpx); | |
141 | int xhi = ceil(xpx); | |
142 | int ylo = floor(ypx); | |
143 | int yhi = ceil(ypx); | |
144 | ||
145 | double s11 = probability(xlo, ylo); | |
146 | double s12 = probability(xlo, yhi); | |
147 | double s22 = probability(xhi, yhi); | |
148 | double s21 = probability(xhi, ylo); | |
149 | ||
150 | double y1 = (xhi - xpx) * s11 + (xpx - xlo) * s21; | |
151 | double y2 = (xhi - xpx) * s12 + (xpx - xlo) * s22; | |
152 | double prob = (yhi - ypx) * y1 + (ypx - ylo) * y2; | |
153 | ||
154 | if (pdist(gen) < prob) { | |
155 | points.push_back(make_pair(x, y)); | |
156 | p++; | |
157 | } | |
158 | } | |
159 | ||
160 | samples = array2d<uint8_t>(sscale * image.w, sscale * image.h); | |
161 | sortimage = pixels<uint8_t>(sscale * image.w, sscale * image.h); | |
162 | ||
163 | for (const auto &pair : points) { | |
164 | const double &xmm = pair.first; | |
165 | const double &ymm = pair.second; | |
166 | int xpx = xmm * ppmm_sortimage; | |
167 | int ypx = ymm * ppmm_sortimage; | |
168 | samples(xpx, ypx) = 255; | |
169 | } | |
170 | } | |
171 | ||
172 | void linefeeder_sample() { | |
173 | std::mt19937 gen(seed); | |
174 | std::uniform_real_distribution<double> pdist(0.0, 1.0); | |
175 | #pragma omp parallel for | |
176 | for (size_t i = 0; i < image.w; i++) { | |
177 | for (size_t j = 0; j < image.h; j++) { | |
178 | samples(i, j) = pdist(gen) < probability(i, j) ? 255 : 0; | |
179 | } | |
180 | } | |
181 | } | |
182 | ||
183 | void linefeeder() { | |
266 | 184 | ofstream nc("_pntlsm.nc"); |
267 | 185 | nc << setprecision(3) |
268 | 186 | << fixed |
304 | 222 | } |
305 | 223 | j++; |
306 | 224 | } |
307 | nc << "G0 Z1\n" | |
308 | << "G0 X0 Y0\n"; | |
225 | } | |
226 | ||
227 | typedef uint64_t bucketidx; | |
228 | ||
229 | bucketidx make_idx(double x, double y) { | |
230 | x /= bucketsize_mm; | |
231 | y /= bucketsize_mm; | |
232 | ||
233 | if (x < 0 || x > std::numeric_limits<uint32_t>::max()) { | |
234 | throw std::out_of_range("bad value for x"); | |
235 | } | |
236 | if (y < 0 || y > std::numeric_limits<uint32_t>::max()) { | |
237 | throw std::out_of_range("bad value for y"); | |
238 | } | |
239 | ||
240 | uint64_t xx = static_cast<uint64_t>(x); | |
241 | uint64_t yy = static_cast<uint64_t>(y); | |
242 | return (yy << 32) | xx; | |
243 | } | |
244 | ||
245 | void arcsort() { | |
246 | typedef vector<pair<double, double>> bucket; | |
247 | unordered_map<bucketidx, bucket> buckets; | |
248 | ||
249 | for (const auto &p : points) { | |
250 | buckets[make_idx(p.first, p.second)].push_back(p); | |
251 | } | |
252 | ||
253 | vector<bucketidx> sorted_bucket_indexes; | |
254 | sorted_bucket_indexes.reserve(buckets.size()); | |
255 | for (const auto &bentry : buckets) { | |
256 | sorted_bucket_indexes.push_back(bentry.first); | |
257 | } | |
258 | sort(sorted_bucket_indexes.begin(), sorted_bucket_indexes.end(), | |
259 | [](const auto &b1, const auto &b2) { | |
260 | uint64_t y1 = b1 & 0xFFFFFFFF00000000; | |
261 | uint64_t y2 = b2 & 0xFFFFFFFF00000000; | |
262 | if (y1 < y2) return true; | |
263 | if (y1 > y2) return false; | |
264 | if (y1 % 2) return b1 > b2; | |
265 | return b1 < b2; | |
266 | }); | |
267 | for (auto &bentry : buckets) { | |
268 | sort(bentry.second.begin(), bentry.second.end(), | |
269 | [](const auto &p1, const auto &p2) { | |
270 | if (p1.first < p2.first) return true; | |
271 | if (p1.first > p2.first) return false; | |
272 | return p1.second < p2.second; | |
273 | }); | |
274 | } | |
275 | ||
276 | size_t pi = 0; | |
277 | for (const auto &bidx : sorted_bucket_indexes) { | |
278 | const bucket &bucket = buckets.at(bidx); | |
279 | for (const auto &p : bucket) { | |
280 | double xpx = p.first * ppmm_sortimage; | |
281 | double ypx = p.second * ppmm_sortimage; | |
282 | double t = (double) pi++ / (double) npoints; | |
283 | uint8_t r, g, b; | |
284 | rainbow(t, r, g, b); | |
285 | sortimage.r(xpx, ypx) = r; | |
286 | sortimage.g(xpx, ypx) = g; | |
287 | sortimage.b(xpx, ypx) = b; | |
288 | } | |
289 | } | |
290 | ||
291 | ofstream nc("_pntlsm.nc"); | |
292 | nc << setprecision(3) << fixed | |
293 | << "G0 Z1\nG0 X0 Y0\n"; | |
294 | for (const auto &bidx : sorted_bucket_indexes) { | |
295 | const bucket &bucket = buckets.at(bidx); | |
296 | // put the conditional outside the loop for maximum loopage | |
297 | if (image.h > image.w) { | |
298 | for (const auto &p : bucket) { | |
299 | nc << "G0 X" << p.second << " Y" << p.first | |
300 | << "\nG1 Z0\nG1 Z1\n"; | |
301 | } | |
302 | } else { | |
303 | for (const auto &p : bucket) { | |
304 | nc << "G0 X" << p.first << " Y" << p.second | |
305 | << "\nG1 Z0\nG1 Z1\n"; | |
306 | } | |
307 | } | |
308 | } | |
309 | nc << "G0 Z1\nG0 X0 Y0\n"; | |
309 | 310 | } |
310 | 311 | |
311 | 312 | bool saveimages() { |
313 | 314 | cout << "failed to save prob image" << endl; |
314 | 315 | return false; |
315 | 316 | } |
317 | ||
316 | 318 | array2d<uint8_t> inverted = samples; |
317 | 319 | for (size_t i = 0; i < inverted.w * inverted.h; i++) { |
318 | 320 | inverted(i) = 255 - inverted(i); |
320 | 322 | if (not to_png("_pntlsm_samples.png", inverted)) { |
321 | 323 | cout << "failed to save sample image" << endl; |
322 | 324 | return false; |
325 | } | |
326 | ||
327 | if (sortimage) { | |
328 | if (not to_png("_pntlsm_output.png", sortimage)) { | |
329 | cout << "failed to save sortimage image" << endl; | |
330 | return false; | |
331 | } | |
323 | 332 | } |
324 | 333 | return true; |
325 | 334 | } |
332 | 341 | |
333 | 342 | if (not load(argv[1])) return EXIT_FAILURE; |
334 | 343 | |
335 | prob(); | |
336 | sample(); | |
344 | genprob(); | |
345 | ||
346 | #if defined(SCATTER) | |
347 | scatter_sample(); | |
348 | arcsort(); // must come before saveimages! | |
349 | #elif defined(FEEDER) | |
350 | linefeeder_sample(); | |
351 | linefeeder(); | |
352 | #else | |
353 | #error no sampling method specified | |
354 | #endif | |
337 | 355 | |
338 | 356 | if (not saveimages()) return EXIT_FAILURE; |
339 | savegcode(); | |
340 | ||
341 | 357 | return EXIT_SUCCESS; |
342 | 358 | } |
343 | 359 |