diff options
author | tamasmeszaros <meszaros.q@gmail.com> | 2020-04-13 13:31:37 +0300 |
---|---|---|
committer | tamasmeszaros <meszaros.q@gmail.com> | 2020-04-23 20:05:32 +0300 |
commit | 247fca6d55318fbaa5c7424bf7605db7fb02e342 (patch) | |
tree | 828924c2f53b07b04260073ebee0c2928a0ba44b /src/libslic3r/MarchingSquares.hpp | |
parent | 6eb51a1cca4d99385f8b1795d5fc1063424c1635 (diff) |
Initial version of sl1 import with sla::Raster refactor.
Diffstat (limited to 'src/libslic3r/MarchingSquares.hpp')
-rw-r--r-- | src/libslic3r/MarchingSquares.hpp | 432 |
1 files changed, 432 insertions, 0 deletions
diff --git a/src/libslic3r/MarchingSquares.hpp b/src/libslic3r/MarchingSquares.hpp new file mode 100644 index 000000000..a5aeb0953 --- /dev/null +++ b/src/libslic3r/MarchingSquares.hpp @@ -0,0 +1,432 @@ +#ifndef MARCHINGSQUARES_HPP +#define MARCHINGSQUARES_HPP + +#include <type_traits> +#include <cstdint> +#include <vector> +#include <algorithm> +#include <cassert> + +namespace marchsq { + +// Marks a square in the grid +struct Coord { + size_t r = 0, c = 0; + + Coord() = default; + explicit Coord(size_t s) : r(s), c(s) {} + Coord(size_t _r, size_t _c): r(_r), c(_c) {} + + size_t seq(const Coord &res) const { return r * res.c + c; } + Coord& operator+=(const Coord& b) { r += b.r; c += b.c; return *this; } + Coord operator+(const Coord& b) const { Coord a = *this; a += b; return a; } +}; + +// Closed ring of cell coordinates +using Ring = std::vector<Coord>; + +// Specialize this struct to register a raster type for the Marching squares alg +template<class T, class Enable = void> struct _RasterTraits { + + // The type of pixel cell in the raster + using ValueType = typename T::ValueType; + + // Value at a given position + static ValueType get(const T &raster, size_t row, size_t col); + + // Number of rows and cols of the raster + static size_t rows(const T &raster); + static size_t cols(const T &raster); +}; + +// Specialize this to use parellel loops within the algorithm +template<class ExecutionPolicy, class Enable = void> struct _Loop { + template<class It, class Fn> static void for_each(It from, It to, Fn &&fn) + { + for (auto it = from; it < to; ++it) fn(*it, size_t(it - from)); + } +}; + +namespace __impl { + +template<class T> using RasterTraits = _RasterTraits<std::decay_t<T>>; +template<class T> using TRasterValue = typename RasterTraits<T>::ValueType; + +template<class T> TRasterValue<T> isoval(const T &raster, const Coord &crd) +{ + return RasterTraits<T>::get(raster, crd.r, crd.c); +} + +template<class T> size_t rows(const T &raster) +{ + return RasterTraits<T>::rows(raster); +} + +template<class T> size_t cols(const T &raster) +{ + return RasterTraits<T>::cols(raster); +} + +template<class ExecutionPolicy, class It, class Fn> +void for_each(ExecutionPolicy&& policy, It from, It to, Fn &&fn) +{ + _Loop<ExecutionPolicy>::for_each(from, to, fn); +} + +// Type of squares (tiles) depending on which vertices are inside an ROI +// The vertices would be marked a, b, c, d in counter clockwise order from the +// bottom left vertex of a square. +// d --- c +// | | +// | | +// a --- b +enum class SquareTag : uint8_t { +// 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 + none, a, b, ab, c, ac, bc, abc, d, ad, bd, abd, cd, acd, bcd, full +}; + +template<class E> constexpr std::underlying_type_t<E> _t(E e) noexcept +{ + return static_cast<std::underlying_type_t<E>>(e); +} + +enum class Dir: uint8_t { left, down, right, up, none}; + +static const constexpr Dir NEXT_CCW[] = { + /* 00 */ Dir::none, // SquareTag::none (empty square, nowhere to go) + /* 01 */ Dir::left, // SquareTag::a + /* 02 */ Dir::down, // SquareTag::b + /* 03 */ Dir::left, // SquareTag::ab + /* 04 */ Dir::right, // SquareTag::c + /* 05 */ Dir::none, // SquareTag::ac (ambiguous case) + /* 06 */ Dir::down, // SquareTag::bc + /* 07 */ Dir::left, // SquareTag::abc + /* 08 */ Dir::up, // SquareTag::d + /* 09 */ Dir::up, // SquareTag::ad + /* 10 */ Dir::none, // SquareTag::bd (ambiguous case) + /* 11 */ Dir::up, // SquareTag::abd + /* 12 */ Dir::right, // SquareTag::cd + /* 13 */ Dir::right, // SquareTag::acd + /* 14 */ Dir::down, // SquareTag::bcd + /* 15 */ Dir::none // SquareTag::full (full covered, nowhere to go) +}; + +static const constexpr uint8_t PREV_CCW[] = { + /* 00 */ 1 << _t(Dir::none), + /* 01 */ 1 << _t(Dir::up), + /* 02 */ 1 << _t(Dir::left), + /* 03 */ 1 << _t(Dir::left), + /* 04 */ 1 << _t(Dir::down), + /* 05 */ 1 << _t(Dir::up) | 1 << _t(Dir::down), + /* 06 */ 1 << _t(Dir::down), + /* 07 */ 1 << _t(Dir::down), + /* 08 */ 1 << _t(Dir::right), + /* 09 */ 1 << _t(Dir::up), + /* 10 */ 1 << _t(Dir::left) | 1 << _t(Dir::right), + /* 11 */ 1 << _t(Dir::left), + /* 12 */ 1 << _t(Dir::right), + /* 13 */ 1 << _t(Dir::up), + /* 14 */ 1 << _t(Dir::right), + /* 15 */ 1 << _t(Dir::none) +}; + +const constexpr uint8_t DIRMASKS[] = { + /*left: */ 0x01, /*down*/ 0x12, /*right */0x21, /*up*/ 0x10, /*none*/ 0x00 +}; + +inline Coord step(const Coord &crd, Dir d) +{ + uint8_t dd = DIRMASKS[uint8_t(d)]; + return {crd.r - 1 + (dd & 0x0f), crd.c - 1 + (dd >> 4)}; +} + +template<class Rst> class Grid { + const Rst * m_rst = nullptr; + Coord m_cellsize, m_res_1, m_window, m_gridsize; + std::vector<uint8_t> m_tags; // Assign tags to each square + + Coord rastercoord(const Coord &crd) const + { + return {crd.r * m_window.r, crd.c * m_window.c}; + } + + Coord bl(const Coord &crd) const { return tl(crd) + Coord{m_res_1.r, 0}; } + Coord br(const Coord &crd) const { return tl(crd) + Coord{m_res_1.r, m_res_1.c}; } + Coord tr(const Coord &crd) const { return tl(crd) + Coord{0, m_res_1.c}; } + Coord tl(const Coord &crd) const { return rastercoord(crd); } + + TRasterValue<Rst> bottomleft(const Coord &cell) const + { + return isoval(*m_rst, bl(cell)); + } + + TRasterValue<Rst> bottomright(const Coord &cell) const + { + return isoval(*m_rst, br(cell)); + } + + TRasterValue<Rst> topright(const Coord &cell) const + { + return isoval(*m_rst, tr(cell)); + } + + TRasterValue<Rst> topleft(const Coord &cell) const + { + return isoval(*m_rst, tl(cell)); + } + + // Calculate the tag for a cell (or square). The cell coordinates mark the + // top left vertex of a square in the raster. v is the isovalue + uint8_t get_tag_for_cell(const Coord &cell, TRasterValue<Rst> v) + { + uint8_t t = (bottomleft(cell) >= v) + + ((bottomright(cell) >= v) << 1) + + ((topright(cell) >= v) << 2) + + ((topleft(cell) >= v) << 3); + + assert(t < 16); + return t; + } + + // Get a cell coordinate from a sequential index + Coord coord(size_t i) const { return {i / m_gridsize.c, i % m_gridsize.c}; } + + size_t seq(const Coord &crd) const { return crd.seq(m_gridsize); } + + bool is_visited(size_t idx, Dir d = Dir::none) const + { + SquareTag t = get_tag(idx); + uint8_t ref = d == Dir::none ? PREV_CCW[_t(t)] : uint8_t(1 << _t(d)); + return t == SquareTag::full || t == SquareTag::none || + ((m_tags[idx] & 0xf0) >> 4) == ref; + } + + void set_visited(size_t idx, Dir d = Dir::none) + { + m_tags[idx] |= (1 << (_t(d)) << 4); + } + + bool is_ambiguous(size_t idx) const + { + SquareTag t = get_tag(idx); + return t == SquareTag::ac || t == SquareTag::bd; + } + + // Search for a new starting square + size_t search_start_cell(size_t i = 0) const + { + // Skip ambiguous tags as starting tags due to unknown previous + // direction. + while ((i < m_tags.size() && is_visited(i)) || is_ambiguous(i)) ++i; + + return i; + } + + SquareTag get_tag(size_t idx) const { return SquareTag(m_tags[idx] & 0x0f); } + + Dir next_dir(Dir prev, SquareTag tag) const + { + // Treat ambiguous cases as two separate regions in one square. + switch (tag) { + case SquareTag::ac: + switch (prev) { + case Dir::down: return Dir::right; + case Dir::up: return Dir::left; + default: assert(false); return Dir::none; + } + case SquareTag::bd: + switch (prev) { + case Dir::right: return Dir::up; + case Dir::left: return Dir::down; + default: assert(false); return Dir::none; + } + default: + return NEXT_CCW[uint8_t(tag)]; + } + + return Dir::none; + } + + struct CellIt { + Coord crd; Dir dir= Dir::none; const Rst *rst = nullptr; + TRasterValue<Rst> operator*() const { return isoval(*rst, crd); } + CellIt& operator++() { crd = step(crd, dir); return *this; } + CellIt operator++(int) { CellIt it = *this; ++(*this); return it; } + bool operator!=(const CellIt &it) { return crd.r != it.crd.r || crd.c != it.crd.c; } + + using value_type = TRasterValue<Rst>; + using pointer = TRasterValue<Rst> *; + using reference = TRasterValue<Rst> &; + using difference_type = long; + using iterator_category = std::forward_iterator_tag; + }; + + // Two cell iterators representing an edge of a square. This is then + // used for binary search for the first active pixel on the edge. + struct Edge { CellIt from, to; }; + + Edge edge(const Coord &ringvertex) + { + size_t idx = ringvertex.r; + Coord cell = coord(idx); + uint8_t tg = m_tags[ringvertex.r]; + SquareTag t = SquareTag(tg & 0x0f); + + switch (t) { + case SquareTag::a: + case SquareTag::ab: + case SquareTag::abc: + return {{tl(cell), Dir::down, m_rst}, {bl(cell)}}; + case SquareTag::b: + case SquareTag::bc: + case SquareTag::bcd: + return {{bl(cell), Dir::right, m_rst}, {br(cell)}}; + case SquareTag::c: + return {{br(cell), Dir::up, m_rst}, {tr(cell)}}; + case SquareTag::ac: + switch (Dir(ringvertex.c)) { + case Dir::left: return {{tl(cell), Dir::down, m_rst}, {bl(cell)}}; + case Dir::right: return {{br(cell), Dir::up, m_rst}, {tr(cell)}}; + default: assert(false); + } + case SquareTag::d: + case SquareTag::ad: + case SquareTag::abd: + return {{tr(cell), Dir::left, m_rst}, {tl(cell)}}; + case SquareTag::bd: + switch (Dir(ringvertex.c)) { + case Dir::down: return {{bl(cell), Dir::right, m_rst}, {br(cell)}}; + case Dir::up: return {{tr(cell), Dir::left, m_rst}, {tl(cell)}}; + default: assert(false); + } + case SquareTag::cd: + case SquareTag::acd: + return {{br(cell), Dir::up, m_rst}, {tr(cell)}}; + case SquareTag::full: + case SquareTag::none: { + Coord crd{tl(cell) + Coord{m_cellsize.r / 2, m_cellsize.c / 2}}; + return {{crd, Dir::none, m_rst}, crd}; + } + } + + return {}; + } + +public: + explicit Grid(const Rst &rst, const Coord &cellsz, const Coord &overlap) + : m_rst{&rst} + , m_cellsize{cellsz} + , m_res_1{m_cellsize.r - 1, m_cellsize.c - 1} + , m_window{overlap.r < cellsz.r ? cellsz.r - overlap.r : cellsz.r, + overlap.c < cellsz.c ? cellsz.c - overlap.c : cellsz.c} + , m_gridsize{(rows(rst) - overlap.r) / m_window.r, + (cols(rst) - overlap.c) / m_window.c} + , m_tags(m_gridsize.r * m_gridsize.c, 0) + {} + + // Go through the cells and mark them with the appropriate tag. + template<class ExecutionPolicy> + void tag_grid(ExecutionPolicy &&policy, TRasterValue<Rst> isoval) + { + // parallel for r + for_each (std::forward<ExecutionPolicy>(policy), + m_tags.begin(), m_tags.end(), + [this, isoval](uint8_t& tag, size_t idx) { + tag = get_tag_for_cell(coord(idx), isoval); + }); + } + + // Scan for the rings on the tagged grid. Each ring vertex stores the + // sequential index of the cell and the next direction (Dir). + // This info can be used later to calculate the exact raster coordinate. + std::vector<Ring> scan_rings() + { + std::vector<Ring> rings; + size_t startidx = 0; + while ((startidx = search_start_cell(startidx)) < m_tags.size()) { + Ring ring; + + size_t idx = startidx; + Dir prev = Dir::none, next = next_dir(prev, get_tag(idx)); + + while (next != Dir::none && !is_visited(idx, prev)) { + Coord ringvertex{idx, size_t(next)}; + ring.emplace_back(ringvertex); + set_visited(idx, prev); + + idx = seq(step(coord(idx), next)); + prev = next; + next = next_dir(next, get_tag(idx)); + } + + // To prevent infinite loops in case of degenerate input + if (next == Dir::none) m_tags[startidx] = _t(SquareTag::none); + + if (ring.size() > 1) { + ring.pop_back(); + rings.emplace_back(ring); + } + } + + return rings; + } + + // Calculate the exact raster position from the cells which store the + // sequantial index of the square and the next direction + template<class ExecutionPolicy> + void interpolate_rings(ExecutionPolicy && policy, + std::vector<Ring> &rings, + TRasterValue<Rst> isov) + { + for_each(std::forward<ExecutionPolicy>(policy), + rings.begin(), rings.end(), [this, isov] (Ring &ring, size_t) { + for (Coord &ringvertex : ring) { + Edge e = edge(ringvertex); + CellIt found = std::lower_bound(e.from, e.to, isov); + ringvertex = found.crd; + } + }); + } +}; + +template<class Raster, class ExecutionPolicy> +std::vector<marchsq::Ring> execute_with_policy(ExecutionPolicy && policy, + const Raster & raster, + TRasterValue<Raster> isoval, + Coord windowsize = {}) +{ + if (!rows(raster) || !cols(raster)) return {}; + + size_t ratio = cols(raster) / rows(raster); + + if (!windowsize.r) windowsize.r = 2; + if (!windowsize.c) + windowsize.c = std::max(size_t(2), windowsize.r * ratio); + + Coord overlap{1}; + + Grid<Raster> grid{raster, windowsize, overlap}; + + grid.tag_grid(std::forward<ExecutionPolicy>(policy), isoval); + std::vector<marchsq::Ring> rings = grid.scan_rings(); + grid.interpolate_rings(std::forward<ExecutionPolicy>(policy), rings, isoval); + + return rings; +} + +template<class Raster> +std::vector<marchsq::Ring> execute(const Raster &raster, + TRasterValue<Raster> isoval, + Coord windowsize = {}) +{ + return execute_with_policy(nullptr, raster, isoval, windowsize); +} + +} // namespace __impl + +using __impl::execute_with_policy; +using __impl::execute; + +} // namespace marchsq + +#endif // MARCHINGSQUARES_HPP |