1 #ifndef MARCHINGSQUARES_HPP
2 #define MARCHINGSQUARES_HPP
3
4 #include <type_traits>
5 #include <cstdint>
6 #include <vector>
7 #include <algorithm>
8 #include <cassert>
9
10 namespace marchsq {
11
12 // Marks a square in the grid
13 struct Coord {
14 long r = 0, c = 0;
15
16 Coord() = default;
Coordmarchsq::Coord17 explicit Coord(long s) : r(s), c(s) {}
Coordmarchsq::Coord18 Coord(long _r, long _c): r(_r), c(_c) {}
19
seqmarchsq::Coord20 size_t seq(const Coord &res) const { return r * res.c + c; }
operator +=marchsq::Coord21 Coord& operator+=(const Coord& b) { r += b.r; c += b.c; return *this; }
operator +marchsq::Coord22 Coord operator+(const Coord& b) const { Coord a = *this; a += b; return a; }
23 };
24
25 // Closed ring of cell coordinates
26 using Ring = std::vector<Coord>;
27
28 // Specialize this struct to register a raster type for the Marching squares alg
29 template<class T, class Enable = void> struct _RasterTraits {
30
31 // The type of pixel cell in the raster
32 using ValueType = typename T::ValueType;
33
34 // Value at a given position
35 static ValueType get(const T &raster, size_t row, size_t col);
36
37 // Number of rows and cols of the raster
38 static size_t rows(const T &raster);
39 static size_t cols(const T &raster);
40 };
41
42 // Specialize this to use parellel loops within the algorithm
43 template<class ExecutionPolicy, class Enable = void> struct _Loop {
for_eachmarchsq::_Loop44 template<class It, class Fn> static void for_each(It from, It to, Fn &&fn)
45 {
46 for (auto it = from; it < to; ++it) fn(*it, size_t(it - from));
47 }
48 };
49
50 namespace __impl {
51
52 template<class T> using RasterTraits = _RasterTraits<std::decay_t<T>>;
53 template<class T> using TRasterValue = typename RasterTraits<T>::ValueType;
54
rows(const T & raster)55 template<class T> size_t rows(const T &raster)
56 {
57 return RasterTraits<T>::rows(raster);
58 }
59
cols(const T & raster)60 template<class T> size_t cols(const T &raster)
61 {
62 return RasterTraits<T>::cols(raster);
63 }
64
isoval(const T & rst,const Coord & crd)65 template<class T> TRasterValue<T> isoval(const T &rst, const Coord &crd)
66 {
67 return RasterTraits<T>::get(rst, crd.r, crd.c);
68 }
69
70 template<class ExecutionPolicy, class It, class Fn>
for_each(ExecutionPolicy && policy,It from,It to,Fn && fn)71 void for_each(ExecutionPolicy&& policy, It from, It to, Fn &&fn)
72 {
73 _Loop<ExecutionPolicy>::for_each(from, to, fn);
74 }
75
76 // Type of squares (tiles) depending on which vertices are inside an ROI
77 // The vertices would be marked a, b, c, d in counter clockwise order from the
78 // bottom left vertex of a square.
79 // d --- c
80 // | |
81 // | |
82 // a --- b
83 enum class SquareTag : uint8_t {
84 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
85 none, a, b, ab, c, ac, bc, abc, d, ad, bd, abd, cd, acd, bcd, full
86 };
87
_t(E e)88 template<class E> constexpr std::underlying_type_t<E> _t(E e) noexcept
89 {
90 return static_cast<std::underlying_type_t<E>>(e);
91 }
92
93 enum class Dir: uint8_t { left, down, right, up, none};
94
95 static const constexpr Dir NEXT_CCW[] = {
96 /* 00 */ Dir::none, // SquareTag::none (empty square, nowhere to go)
97 /* 01 */ Dir::left, // SquareTag::a
98 /* 02 */ Dir::down, // SquareTag::b
99 /* 03 */ Dir::left, // SquareTag::ab
100 /* 04 */ Dir::right, // SquareTag::c
101 /* 05 */ Dir::none, // SquareTag::ac (ambiguous case)
102 /* 06 */ Dir::down, // SquareTag::bc
103 /* 07 */ Dir::left, // SquareTag::abc
104 /* 08 */ Dir::up, // SquareTag::d
105 /* 09 */ Dir::up, // SquareTag::ad
106 /* 10 */ Dir::none, // SquareTag::bd (ambiguous case)
107 /* 11 */ Dir::up, // SquareTag::abd
108 /* 12 */ Dir::right, // SquareTag::cd
109 /* 13 */ Dir::right, // SquareTag::acd
110 /* 14 */ Dir::down, // SquareTag::bcd
111 /* 15 */ Dir::none // SquareTag::full (full covered, nowhere to go)
112 };
113
114 static const constexpr uint8_t PREV_CCW[] = {
115 /* 00 */ 1 << _t(Dir::none),
116 /* 01 */ 1 << _t(Dir::up),
117 /* 02 */ 1 << _t(Dir::left),
118 /* 03 */ 1 << _t(Dir::left),
119 /* 04 */ 1 << _t(Dir::down),
120 /* 05 */ 1 << _t(Dir::up) | 1 << _t(Dir::down),
121 /* 06 */ 1 << _t(Dir::down),
122 /* 07 */ 1 << _t(Dir::down),
123 /* 08 */ 1 << _t(Dir::right),
124 /* 09 */ 1 << _t(Dir::up),
125 /* 10 */ 1 << _t(Dir::left) | 1 << _t(Dir::right),
126 /* 11 */ 1 << _t(Dir::left),
127 /* 12 */ 1 << _t(Dir::right),
128 /* 13 */ 1 << _t(Dir::up),
129 /* 14 */ 1 << _t(Dir::right),
130 /* 15 */ 1 << _t(Dir::none)
131 };
132
133 const constexpr uint8_t DIRMASKS[] = {
134 /*left: */ 0x01, /*down*/ 0x12, /*right */0x21, /*up*/ 0x10, /*none*/ 0x00
135 };
136
step(const Coord & crd,Dir d)137 inline Coord step(const Coord &crd, Dir d)
138 {
139 uint8_t dd = DIRMASKS[uint8_t(d)];
140 return {crd.r - 1 + (dd & 0x0f), crd.c - 1 + (dd >> 4)};
141 }
142
143 template<class Rst> class Grid {
144 const Rst * m_rst = nullptr;
145 Coord m_cellsize, m_res_1, m_window, m_gridsize, m_grid_1;
146 std::vector<uint8_t> m_tags; // Assign tags to each square
147
rastercoord(const Coord & crd) const148 Coord rastercoord(const Coord &crd) const
149 {
150 return {(crd.r - 1) * m_window.r, (crd.c - 1) * m_window.c};
151 }
152
bl(const Coord & crd) const153 Coord bl(const Coord &crd) const { return tl(crd) + Coord{m_res_1.r, 0}; }
br(const Coord & crd) const154 Coord br(const Coord &crd) const { return tl(crd) + Coord{m_res_1.r, m_res_1.c}; }
tr(const Coord & crd) const155 Coord tr(const Coord &crd) const { return tl(crd) + Coord{0, m_res_1.c}; }
tl(const Coord & crd) const156 Coord tl(const Coord &crd) const { return rastercoord(crd); }
157
is_within(const Coord & crd)158 bool is_within(const Coord &crd)
159 {
160 long R = rows(*m_rst), C = cols(*m_rst);
161 return crd.r >= 0 && crd.r < R && crd.c >= 0 && crd.c < C;
162 };
163
164 // Calculate the tag for a cell (or square). The cell coordinates mark the
165 // top left vertex of a square in the raster. v is the isovalue
get_tag_for_cell(const Coord & cell,TRasterValue<Rst> v)166 uint8_t get_tag_for_cell(const Coord &cell, TRasterValue<Rst> v)
167 {
168 Coord sqr[] = {bl(cell), br(cell), tr(cell), tl(cell)};
169
170 uint8_t t = ((is_within(sqr[0]) && isoval(*m_rst, sqr[0]) >= v)) +
171 ((is_within(sqr[1]) && isoval(*m_rst, sqr[1]) >= v) << 1) +
172 ((is_within(sqr[2]) && isoval(*m_rst, sqr[2]) >= v) << 2) +
173 ((is_within(sqr[3]) && isoval(*m_rst, sqr[3]) >= v) << 3);
174
175 assert(t < 16);
176 return t;
177 }
178
179 // Get a cell coordinate from a sequential index
coord(size_t i) const180 Coord coord(size_t i) const
181 {
182 return {long(i) / m_gridsize.c, long(i) % m_gridsize.c};
183 }
184
seq(const Coord & crd) const185 size_t seq(const Coord &crd) const { return crd.seq(m_gridsize); }
186
is_visited(size_t idx,Dir d=Dir::none) const187 bool is_visited(size_t idx, Dir d = Dir::none) const
188 {
189 SquareTag t = get_tag(idx);
190 uint8_t ref = d == Dir::none ? PREV_CCW[_t(t)] : uint8_t(1 << _t(d));
191 return t == SquareTag::full || t == SquareTag::none ||
192 ((m_tags[idx] & 0xf0) >> 4) == ref;
193 }
194
set_visited(size_t idx,Dir d=Dir::none)195 void set_visited(size_t idx, Dir d = Dir::none)
196 {
197 m_tags[idx] |= (1 << (_t(d)) << 4);
198 }
199
is_ambiguous(size_t idx) const200 bool is_ambiguous(size_t idx) const
201 {
202 SquareTag t = get_tag(idx);
203 return t == SquareTag::ac || t == SquareTag::bd;
204 }
205
206 // Search for a new starting square
search_start_cell(size_t i=0) const207 size_t search_start_cell(size_t i = 0) const
208 {
209 // Skip ambiguous tags as starting tags due to unknown previous
210 // direction.
211 while ((i < m_tags.size()) && (is_visited(i) || is_ambiguous(i))) ++i;
212
213 return i;
214 }
215
get_tag(size_t idx) const216 SquareTag get_tag(size_t idx) const { return SquareTag(m_tags[idx] & 0x0f); }
217
next_dir(Dir prev,SquareTag tag) const218 Dir next_dir(Dir prev, SquareTag tag) const
219 {
220 // Treat ambiguous cases as two separate regions in one square.
221 switch (tag) {
222 case SquareTag::ac:
223 switch (prev) {
224 case Dir::down: return Dir::right;
225 case Dir::up: return Dir::left;
226 default: assert(false); return Dir::none;
227 }
228 case SquareTag::bd:
229 switch (prev) {
230 case Dir::right: return Dir::up;
231 case Dir::left: return Dir::down;
232 default: assert(false); return Dir::none;
233 }
234 default:
235 return NEXT_CCW[uint8_t(tag)];
236 }
237
238 return Dir::none;
239 }
240
241 struct CellIt {
242 Coord crd; Dir dir= Dir::none; const Rst *grid = nullptr;
243
operator *marchsq::__impl::Grid::CellIt244 TRasterValue<Rst> operator*() const { return isoval(*grid, crd); }
operator ++marchsq::__impl::Grid::CellIt245 CellIt& operator++() { crd = step(crd, dir); return *this; }
operator ++marchsq::__impl::Grid::CellIt246 CellIt operator++(int) { CellIt it = *this; ++(*this); return it; }
operator !=marchsq::__impl::Grid::CellIt247 bool operator!=(const CellIt &it) { return crd.r != it.crd.r || crd.c != it.crd.c; }
248
249 using value_type = TRasterValue<Rst>;
250 using pointer = TRasterValue<Rst> *;
251 using reference = TRasterValue<Rst> &;
252 using difference_type = long;
253 using iterator_category = std::forward_iterator_tag;
254 };
255
256 // Two cell iterators representing an edge of a square. This is then
257 // used for binary search for the first active pixel on the edge.
258 struct Edge { CellIt from, to; };
259
_edge(const Coord & ringvertex) const260 Edge _edge(const Coord &ringvertex) const
261 {
262 size_t idx = ringvertex.r;
263 Coord cell = coord(idx);
264 uint8_t tg = m_tags[ringvertex.r];
265 SquareTag t = SquareTag(tg & 0x0f);
266
267 switch (t) {
268 case SquareTag::a:
269 case SquareTag::ab:
270 case SquareTag::abc:
271 return {{tl(cell), Dir::down, m_rst}, {bl(cell)}};
272 case SquareTag::b:
273 case SquareTag::bc:
274 case SquareTag::bcd:
275 return {{bl(cell), Dir::right, m_rst}, {br(cell)}};
276 case SquareTag::c:
277 return {{br(cell), Dir::up, m_rst}, {tr(cell)}};
278 case SquareTag::ac:
279 switch (Dir(ringvertex.c)) {
280 case Dir::left: return {{tl(cell), Dir::down, m_rst}, {bl(cell)}};
281 case Dir::right: return {{br(cell), Dir::up, m_rst}, {tr(cell)}};
282 default: assert(false);
283 }
284 case SquareTag::d:
285 case SquareTag::ad:
286 case SquareTag::abd:
287 return {{tr(cell), Dir::left, m_rst}, {tl(cell)}};
288 case SquareTag::bd:
289 switch (Dir(ringvertex.c)) {
290 case Dir::down: return {{bl(cell), Dir::right, m_rst}, {br(cell)}};
291 case Dir::up: return {{tr(cell), Dir::left, m_rst}, {tl(cell)}};
292 default: assert(false);
293 }
294 case SquareTag::cd:
295 case SquareTag::acd:
296 return {{br(cell), Dir::up, m_rst}, {tr(cell)}};
297 case SquareTag::full:
298 case SquareTag::none: {
299 Coord crd{tl(cell) + Coord{m_cellsize.r / 2, m_cellsize.c / 2}};
300 return {{crd, Dir::none, m_rst}, crd};
301 }
302 }
303
304 return {};
305 }
306
edge(const Coord & ringvertex) const307 Edge edge(const Coord &ringvertex) const
308 {
309 const long R = rows(*m_rst), C = cols(*m_rst);
310 const long R_1 = R - 1, C_1 = C - 1;
311
312 Edge e = _edge(ringvertex);
313 e.to.dir = e.from.dir;
314 ++e.to;
315
316 e.from.crd.r = std::min(e.from.crd.r, R_1);
317 e.from.crd.r = std::max(e.from.crd.r, 0l);
318 e.from.crd.c = std::min(e.from.crd.c, C_1);
319 e.from.crd.c = std::max(e.from.crd.c, 0l);
320
321 e.to.crd.r = std::min(e.to.crd.r, R);
322 e.to.crd.r = std::max(e.to.crd.r, 0l);
323 e.to.crd.c = std::min(e.to.crd.c, C);
324 e.to.crd.c = std::max(e.to.crd.c, 0l);
325
326 return e;
327 }
328
329 public:
Grid(const Rst & rst,const Coord & cellsz,const Coord & overlap)330 explicit Grid(const Rst &rst, const Coord &cellsz, const Coord &overlap)
331 : m_rst{&rst}
332 , m_cellsize{cellsz}
333 , m_res_1{m_cellsize.r - 1, m_cellsize.c - 1}
334 , m_window{overlap.r < cellsz.r ? cellsz.r - overlap.r : cellsz.r,
335 overlap.c < cellsz.c ? cellsz.c - overlap.c : cellsz.c}
336 , m_gridsize{2 + (long(rows(rst)) - overlap.r) / m_window.r,
337 2 + (long(cols(rst)) - overlap.c) / m_window.c}
338 , m_tags(m_gridsize.r * m_gridsize.c, 0)
339 {}
340
341 // Go through the cells and mark them with the appropriate tag.
342 template<class ExecutionPolicy>
tag_grid(ExecutionPolicy && policy,TRasterValue<Rst> isoval)343 void tag_grid(ExecutionPolicy &&policy, TRasterValue<Rst> isoval)
344 {
345 // parallel for r
346 for_each (std::forward<ExecutionPolicy>(policy),
347 m_tags.begin(), m_tags.end(),
348 [this, isoval](uint8_t& tag, size_t idx) {
349 tag = get_tag_for_cell(coord(idx), isoval);
350 });
351 }
352
353 // Scan for the rings on the tagged grid. Each ring vertex stores the
354 // sequential index of the cell and the next direction (Dir).
355 // This info can be used later to calculate the exact raster coordinate.
scan_rings()356 std::vector<Ring> scan_rings()
357 {
358 std::vector<Ring> rings;
359 size_t startidx = 0;
360 while ((startidx = search_start_cell(startidx)) < m_tags.size()) {
361 Ring ring;
362
363 size_t idx = startidx;
364 Dir prev = Dir::none, next = next_dir(prev, get_tag(idx));
365
366 while (next != Dir::none && !is_visited(idx, prev)) {
367 Coord ringvertex{long(idx), long(next)};
368 ring.emplace_back(ringvertex);
369 set_visited(idx, prev);
370
371 idx = seq(step(coord(idx), next));
372 prev = next;
373 next = next_dir(next, get_tag(idx));
374 }
375
376 // To prevent infinite loops in case of degenerate input
377 if (next == Dir::none) m_tags[startidx] = _t(SquareTag::none);
378
379 if (ring.size() > 1) {
380 ring.pop_back();
381 rings.emplace_back(ring);
382 }
383 }
384
385 return rings;
386 }
387
388 // Calculate the exact raster position from the cells which store the
389 // sequantial index of the square and the next direction
390 template<class ExecutionPolicy>
interpolate_rings(ExecutionPolicy && policy,std::vector<Ring> & rings,TRasterValue<Rst> isov)391 void interpolate_rings(ExecutionPolicy && policy,
392 std::vector<Ring> &rings,
393 TRasterValue<Rst> isov)
394 {
395 for_each(std::forward<ExecutionPolicy>(policy),
396 rings.begin(), rings.end(), [this, isov] (Ring &ring, size_t)
397 {
398 for (Coord &ringvertex : ring) {
399 Edge e = edge(ringvertex);
400
401 CellIt found = std::lower_bound(e.from, e.to, isov);
402 ringvertex = found.crd;
403 }
404 });
405 }
406 };
407
408 template<class Raster, class ExecutionPolicy>
execute_with_policy(ExecutionPolicy && policy,const Raster & raster,TRasterValue<Raster> isoval,Coord windowsize={})409 std::vector<marchsq::Ring> execute_with_policy(ExecutionPolicy && policy,
410 const Raster & raster,
411 TRasterValue<Raster> isoval,
412 Coord windowsize = {})
413 {
414 if (!rows(raster) || !cols(raster)) return {};
415
416 size_t ratio = cols(raster) / rows(raster);
417
418 if (!windowsize.r) windowsize.r = 2;
419 if (!windowsize.c)
420 windowsize.c = std::max(2l, long(windowsize.r * ratio));
421
422 Coord overlap{1};
423
424 Grid<Raster> grid{raster, windowsize, overlap};
425
426 grid.tag_grid(std::forward<ExecutionPolicy>(policy), isoval);
427 std::vector<marchsq::Ring> rings = grid.scan_rings();
428 grid.interpolate_rings(std::forward<ExecutionPolicy>(policy), rings, isoval);
429
430 return rings;
431 }
432
433 template<class Raster>
execute(const Raster & raster,TRasterValue<Raster> isoval,Coord windowsize={})434 std::vector<marchsq::Ring> execute(const Raster &raster,
435 TRasterValue<Raster> isoval,
436 Coord windowsize = {})
437 {
438 return execute_with_policy(nullptr, raster, isoval, windowsize);
439 }
440
441 } // namespace __impl
442
443 using __impl::execute_with_policy;
444 using __impl::execute;
445
446 } // namespace marchsq
447
448 #endif // MARCHINGSQUARES_HPP
449