diff --git a/src/libslic3r/Format/SL1.cpp b/src/libslic3r/Format/SL1.cpp index 3973ab27af..5796c108e5 100644 --- a/src/libslic3r/Format/SL1.cpp +++ b/src/libslic3r/Format/SL1.cpp @@ -305,8 +305,8 @@ ConfigSubstitutions import_sla_archive( std::function progr) { // Ensure minimum window size for marching squares - windowsize.x() = std::max(2, windowsize.x()); - windowsize.y() = std::max(2, windowsize.y()); + windowsize.x() = std::max(1, windowsize.x()); + windowsize.y() = std::max(1, windowsize.y()); std::string exclude_entries{"thumbnail"}; ArchiveData arch = extract_sla_archive(zipfname, exclude_entries); diff --git a/src/libslic3r/MarchingSquares.hpp b/src/libslic3r/MarchingSquares.hpp index 8370c8cef9..7eb0b25806 100644 --- a/src/libslic3r/MarchingSquares.hpp +++ b/src/libslic3r/MarchingSquares.hpp @@ -1,49 +1,160 @@ #ifndef MARCHINGSQUARES_HPP #define MARCHINGSQUARES_HPP +#include "Execution/ExecutionTBB.hpp" #include #include #include #include #include +// Marching squares +// +// This algorithm generates 2D contour rings for a 2D scalar field (height +// map). See https://en.wikipedia.org/wiki/Marching_squares for details. +// +// This algorithm uses square cells (tiles) with corner tags that are set +// depending on which vertices are inside an ROI (higher than the selected +// isovalue threshold). We label the vertices counter-clockwise, and for ease +// of extraction from the m_tags bitmap we put them into a four bit bitmap in +// little-endian row-column order ; +// +// ^ u (up) ----> column "c" direction +// | +// a ----- d +// | | | +// l <--| |--> r | increasing row "r" direction +// (left) | | (right) \/ +// b ----- c +// | +// \/ d (down) +// +// The m_tags bitmap has the top-left "a" corner of each cell packed +// into uint32_t values stored in little-endian order, with the "a" tag for +// the first cell (r=0,c=0) in the LSB bit0 of the first 32bit word of the +// first row. The grid has a one cell border round the raster area that is +// always clear so lines never extend out of the grid. +// +// From this grid of points we can get the tags bitmap for the four +// corners of each cell with "a" in the LSB bit0 and "c" in bit3. Note the +// bits b,c,d are also the "a" bits of the cell's next-colum/row neighbour +// cells. Cells in the right-most column get their cleared right-side border +// tags by wrapping around the grid and getting the cleared values +// from the left-side border. Note the grid bitmap always includes at least +// one extra cleared boarder bit for the wrap-around of the last tag in the +// last cell. Cells in the bottom row get their cleared bottom-side border +// tags by leaving them clear because their grid index exceeds the +// grid size. +// +// These tag bits are then translated into up/down/left/right cell-exit +// directions indicating the cell edge(s) that any lines will cross and exit +// the cell. As you walk counter-clock-wise around the cell corners in the +// order a->b->c->d, any edge from a unset-corner-tag to an set-corner-tag is +// where a line will cross, exiting the cell into the next cell. +// +// Note any set-to-unset-transition edge is where a line will enter a cell, +// but for marching the lines through the cells we only need to keep the next +// cell exit directions. The exit directions are stored as 4bit +// left/down/right/up flags packed into uint32_t values in m_dirs, in +// counter-clockwise little-endian order. +// +// To build the line-rings we scan through the m_dirs directions for the +// cells until we find one with a set-edge. We then start a line and "march" +// it through the cells following the directions in each cell, recording each +// edge the line crosses that will become a point in the ring. +// +// Note that the outside-edge of the grid extends one cell past the raster +// grid, and this outer edge is always cleared. This means no line will cross +// out of the grid, and they will all eventually connect back to their +// starting cell to become a ring around an island of set tag bits, possibly +// containing rings around lakes of unset tag bits within them. +// +// As we march through the cells and build the lines, we clear the direction +// bits to indicate the edge-transition has been processed. After completing +// a ring we return to scanning for the next non-zeroed starting cell for the +// next ring. Eventually all rings will be found and completed, and the +// m_dirs map will be completely cleared. +// +// Note there are two ambigous cases when the diagonally opposite corner tags +// are set. In these cases two lines will enter and exit the cell. To avoid +// messy lines, we consistently choose the exit direction based on what edge +// the line entered the cell from. This also means when we start scanning for +// the next ring we should start checking at the same cell the last ring +// started from, because it could have a second ring passing through it. +// namespace marchsq { -// Marks a square in the grid -struct Coord { +// Marks a square or point in grid or raster coordintes. +struct Coord +{ long r = 0, c = 0; - + Coord() = default; explicit Coord(long s) : r(s), c(s) {} - Coord(long _r, long _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; } + Coord(long _r, long _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; + } + bool operator==(const Coord& o) const { return (r == o.r) && (c == o.c); } + bool operator!=(const Coord& o) const { return !(*this == o); } }; +inline std::ostream& operator<<(std::ostream& os, const Coord& o) { return os << "(r=" << o.r << ",c=" << o.c << ")"; } + // Closed ring of cell coordinates using Ring = std::vector; -// Specialize this struct to register a raster type for the Marching squares alg -template struct _RasterTraits { - +inline std::ostream& operator<<(std::ostream& os, const Ring& r) +{ + os << "[" << r.size() << "]:"; + for (Coord c : r) + os << " " << c; + return os << "\n"; +} + +// Specialize this struct to register a raster type for MarchingSquares. +template 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); - + 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); + static size_t rows(const T& raster); + static size_t cols(const T& raster); }; // Specialize this to use parellel loops within the algorithm -template struct _Loop { - template static void for_each(It from, It to, Fn &&fn) +template struct _Loop +{ + template static void for_each_idx(It from, It to, Fn&& fn) { - for (auto it = from; it < to; ++it) fn(*it, size_t(it - from)); + for (auto it = from; it < to; ++it) + fn(*it, size_t(it - from)); + } +}; + +// Add Specialization for using ExecutionTBB for parallel loops. +using namespace Slic3r; +template<> struct _Loop +{ + template static void for_each_idx(It from, It to, Fn&& fn) + { + execution::for_each( + ex_tbb, size_t(0), size_t(to - from), [&from, &fn](size_t i) { fn(from[i], i); }, execution::max_concurrency(ex_tbb)); } }; @@ -52,388 +163,481 @@ namespace __impl { template using RasterTraits = _RasterTraits>; template using TRasterValue = typename RasterTraits::ValueType; -template size_t rows(const T &raster) +template size_t rows(const T& raster) { return RasterTraits::rows(raster); } + +template size_t cols(const T& raster) { return RasterTraits::cols(raster); } + +template TRasterValue isoval(const T& rst, const Coord& crd) { return RasterTraits::get(rst, crd.r, crd.c); } + +template void for_each_idx(ExecutionPolicy&& policy, It from, It to, Fn&& fn) { - return RasterTraits::rows(raster); + _Loop::for_each_idx(from, to, fn); } -template size_t cols(const T &raster) -{ - return RasterTraits::cols(raster); -} +template constexpr std::underlying_type_t _t(E e) noexcept { return static_cast>(e); } -template TRasterValue isoval(const T &rst, const Coord &crd) -{ - return RasterTraits::get(rst, crd.r, crd.c); -} - -template -void for_each(ExecutionPolicy&& policy, It from, It to, Fn &&fn) -{ - _Loop::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 +// A set of next direction flags for a cell to indicate what sides(s) +// any yet-to-be-marched line(s) will cross to leave the cell. +enum class Dir : uint8_t { + none = 0b0000, + left = 0b0001, /* exit through a-b edge */ + down = 0b0010, /* exit through b-c edge */ + right = 0b0100, /* exit through c-d edge */ + leftright = left | right, + up = 0b1000, /* exit through d-a edge */ + updown = up | down, + all = 0b1111 }; -template constexpr std::underlying_type_t _t(E e) noexcept +inline std::ostream& operator<<(std::ostream& os, const Dir& d) { return os << ".x##^#X#####"[_t(d)]; } + +// This maps square tag column/row order bitmaps to the next +// direction(s) a line will exit the cell. The directions ensure the +// line-rings circle counter-clock-wise around the clusters of set tags. +static const constexpr std::array NEXT_CCW = [] { + auto map = decltype(NEXT_CCW){}; + for (uint32_t cbda = 0; cbda < 16; cbda++) { + const uint32_t dcba = (cbda & 0b0001) | ((cbda >> 1) & 0b0110) | ((cbda << 2) & 0b1000); + const uint32_t adcb = (dcba >> 1) | ((dcba << 3) & 0b1000); + map[cbda] = ~dcba & adcb; + } + return map; +}(); + +// Step a point in a direction, optionally by n steps. +inline void step(Coord& crd, const Dir d, const long n = 1) { - return static_cast>(e); + switch (d) { + case Dir::left: crd.c -= n; break; + case Dir::down: crd.r += n; break; + case Dir::right: crd.c += n; break; + case Dir::up: crd.r -= n; break; + } } -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) +template class Grid { - uint8_t dd = DIRMASKS[uint8_t(d)]; - return {crd.r - 1 + (dd & 0x0f), crd.c - 1 + (dd >> 4)}; -} + const Rst* m_rst = nullptr; + Coord m_window, m_rastsize, m_gridsize; + size_t m_gridlen; // The number of cells in the grid. + std::vector m_tags; // bit-packed squaretags for each corner. + std::vector m_dirs; // bit-packed next directions for each cell. -template class Grid { - const Rst * m_rst = nullptr; - Coord m_cellsize, m_res_1, m_window, m_gridsize, m_grid_1; - std::vector m_tags; // Assign tags to each square - - Coord rastercoord(const Coord &crd) const + // Convert a grid coordinate point into raster coordinates. + inline Coord rastercoord(const Coord& crd) const { + // Note the -1 offset for the grid border around the raster area. return {(crd.r - 1) * m_window.r, (crd.c - 1) * 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); } - - bool is_within(const Coord &crd) - { - long R = rows(*m_rst), C = cols(*m_rst); - return crd.r >= 0 && crd.r < R && crd.c >= 0 && crd.c < C; - }; + // Get the 4 corners of a grid coordinate cell in raster coordinates. + Coord bl(const Coord& crd) const { return tl(crd) + Coord{m_window.r, 0}; } + Coord br(const Coord& crd) const { return tl(crd) + Coord{m_window.r, m_window.c}; } + Coord tr(const Coord& crd) const { return tl(crd) + Coord{0, m_window.c}; } + Coord tl(const Coord& crd) const { return rastercoord(crd); } - // 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 v) - { - Coord sqr[] = {bl(cell), br(cell), tr(cell), tl(cell)}; - - uint8_t t = ((is_within(sqr[0]) && isoval(*m_rst, sqr[0]) >= v)) + - ((is_within(sqr[1]) && isoval(*m_rst, sqr[1]) >= v) << 1) + - ((is_within(sqr[2]) && isoval(*m_rst, sqr[2]) >= v) << 2) + - ((is_within(sqr[3]) && isoval(*m_rst, sqr[3]) >= v) << 3); - - assert(t < 16); - return t; - } - - // Get a cell coordinate from a sequential index - Coord coord(size_t i) const - { - return {long(i) / m_gridsize.c, long(i) % m_gridsize.c}; - } + // Test if a raster coordinate point is within the raster area. + inline bool is_within(const Coord& crd) const { return crd.r >= 0 && crd.r < m_rastsize.r && crd.c >= 0 && crd.c < m_rastsize.c; } - size_t seq(const Coord &crd) const { return crd.seq(m_gridsize); } - - bool is_visited(size_t idx, Dir d = Dir::none) const + // Get a block of 32 tags for a block index from the raster isovals. + uint32_t get_tags_block32(const size_t bidx, const TRasterValue v) 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; + Coord gcrd = coord(bidx * 32); // position in grid coordinates of the block. + Coord rcrd = rastercoord(gcrd); // position in raster coordinates. + uint32_t tags = 0; + // Set a bit for each corner that has osoval > v. + for (uint32_t b = 1; b; b <<= 1) { + if (is_within(rcrd) && isoval(*m_rst, rcrd) > v) + tags |= b; + gcrd.c += 1; + rcrd.c += m_window.c; + // If we hit the end of the row, start on the next row. + if (gcrd.c >= m_gridsize.c) { + gcrd = Coord(gcrd.r + 1, 0); + rcrd = rastercoord(gcrd); } - 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; + return tags; } - - struct CellIt { - Coord crd; Dir dir= Dir::none; const Rst *grid = nullptr; - - TRasterValue operator*() const { return isoval(*grid, 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; } - + + // Get a block of 8 directions for a block index from the tags. + uint32_t get_dirs_block8(const size_t bidx) const + { + size_t gidx = bidx * 8; + uint32_t dirs = 0; + + // Get the next 9 top-row tags at this grid index into the bottom 16 + // bits, and the 9 bottom-row tags into the top 16 bits. + uint32_t tags9 = get_tags9(gidx) | (get_tags9(gidx + m_gridsize.c) << 16); + // Skip generating dirs if the tags are all 1's or all 0's. + if ((tags9 != 0) && (tags9 != 0x01ff01ff)) { + for (auto s = 0; s < 32; s += 4) { + uint8_t tags = (tags9 & 0b11) | ((tags9 >> 14) & 0b1100); + dirs |= NEXT_CCW[tags] << s; + tags9 >>= 1; + } + } + return dirs; + } + + // Get the next 9 corner tags on a row for building a dirs block at a grid index. + uint32_t get_tags9(size_t gidx) const + { + uint32_t tags = 0; // the tags value. + size_t i = gidx / 32; // the tags block index + int o = gidx % 32; // the tags block offset + + if (gidx < m_gridlen) { + // get the next 9 tags in the row. + tags = (m_tags[i] >> o) & 0x1ff; + // Some of the tags are in the next tags block. + if ((o > (32 - 9)) && ((i + 1) < m_tags.size())) { + tags |= (m_tags[i + 1] << (32 - o)) & 0x1ff; + } + } + return tags; + } + + // Clear directions in a cell's dirs for a grid index. + void clr_dirs(const size_t gidx, const Dir d) + { + assert(gidx < m_gridlen); + size_t i = gidx / 8; // the dirs block index + int o = (gidx % 8); // the dirs block offset + + m_dirs[i] &= ~(static_cast(d) << (o * 4)); + } + + // Get directions in a cell's dirs from the m_dirs store. + Dir get_dirs(const size_t gidx) const + { + assert(gidx < m_gridlen); + size_t i = gidx / 8; // the dirs block index + int o = (gidx % 8); // the dirs block offset + + return Dir((m_dirs[i] >> (o * 4)) & 0b1111); + } + + // Get a cell coordinate from a sequential grid index. + Coord coord(size_t i) const { return {long(i) / m_gridsize.c, long(i) % m_gridsize.c}; } + + // Get a sequential index from a cell coordinate. + size_t seq(const Coord& crd) const { return crd.seq(m_gridsize); } + + // Step a sequential grid index in a direction. + size_t stepidx(const size_t idx, const Dir d) const + { + assert(idx < m_gridlen); + switch (d) { + case Dir::left: return idx - 1; + case Dir::down: return idx + m_gridsize.c; + case Dir::right: return idx + 1; + case Dir::up: return idx - m_gridsize.c; + default: return idx; + } + } + + // Search for a new starting square. + size_t search_start_cell(size_t gidx = 0) const + { + size_t i = gidx / 8; // The dirs block index. + int o = gidx % 8; // The dirs block offset. + + // Skip cells without any unvisited edges. + while (i < m_dirs.size()) { + if (!m_dirs[i]) { + // Whole block is clear, advance to the next block; + i++; + o = 0; + } else { + // Block not clear, find the next non-zero tags in the block. Note + // all dirs before gidx are cleared, so any set bits must be in + // block offsets >= o, not before. + for (uint32_t m = 0b1111 << (o * 4); !(m_dirs[i] & m); m <<= 4) + o++; + break; + } + } + return i * 8 + o; + } + + // Get the next direction for a cell index after the prev direction. + Dir next_dir(size_t idx, Dir prev = Dir::all) const + { + Dir next = get_dirs(idx); + + // Treat ambiguous cases as two separate regions in one square. If + // there are two possible next directions, pick based on the prev + // direction. If prev=all we are starting a new line so pick the one + // that leads "forwards" (right or down) in the search. + switch (next) { + case Dir::leftright: + // We must be coming from up, down, or starting a new line. + assert(prev == Dir::up || prev == Dir::down || prev == Dir::all); + return (prev == Dir::up) ? Dir::left : Dir::right; + case Dir::updown: + // We must be coming from left, right, or starting a new line. + assert(prev == Dir::left || prev == Dir::right || prev == Dir::all); + return (prev == Dir::right) ? Dir::up : Dir::down; + default: + // Next must be a single direction or none to stop. + assert(next == Dir::none || next == Dir::left || next == Dir::down || next == Dir::right || next == Dir::up); + return next; + } + } + + struct CellIt + { + using iterator_category = std::random_access_iterator_tag; using value_type = TRasterValue; - using pointer = TRasterValue *; - using reference = TRasterValue &; + using pointer = TRasterValue*; + using reference = TRasterValue&; using difference_type = long; - using iterator_category = std::forward_iterator_tag; + + Coord crd; + Dir dir = Dir::none; + const Rst* grid = nullptr; + + inline TRasterValue operator*() const { return isoval(*grid, crd); } + inline TRasterValue& operator[](long n) const { return *((*this) + n); } + inline CellIt& operator++() + { + step(crd, dir); + return *this; + } + inline CellIt operator++(int) + { + CellIt o = *this; + ++(*this); + return o; + } + inline CellIt& operator--() + { + step(crd, dir, -1); + return *this; + } + inline CellIt operator--(int) + { + CellIt o = *this; + --(*this); + return o; + } + inline CellIt& operator+=(long n) + { + step(crd, dir, n); + return *this; + } + inline CellIt& operator-=(long n) + { + step(crd, dir, -n); + return *this; + } + inline CellIt operator+(long n) const + { + CellIt o = *this; + o += n; + return o; + } + inline CellIt operator-(long n) const + { + CellIt o = *this; + o -= n; + return o; + } + inline long operator-(const CellIt& o) const + { + switch (dir) { + case Dir::left: return o.crd.c - crd.c; + case Dir::down: return crd.r - o.crd.r; + case Dir::right: return crd.c - o.crd.c; + case Dir::up: return o.crd.r - crd.r; + default: return 0; + } + } + inline bool operator==(const CellIt& o) const { return crd == o.crd; } + inline bool operator!=(const CellIt& o) const { return crd != o.crd; } + inline bool operator<(const CellIt& o) const { return (*this - o) < 0; } + inline bool operator>(const CellIt& o) const { return (*this - o) > 0; } + inline bool operator<=(const CellIt& o) const { return (*this - o) <= 0; } + inline bool operator>=(const CellIt& o) const { return (*this - o) >= 0; } }; - + // 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) const + struct Edge { - 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}}; + CellIt from, to; + }; + + Edge _edge(const Coord& ringvertex) const + { + size_t idx = ringvertex.r; + Coord cell = coord(idx); + Dir d = Dir(ringvertex.c); + + switch (d) { + case Dir::left: return {{tl(cell), Dir::down, m_rst}, {bl(cell)}}; + case Dir::down: return {{bl(cell), Dir::right, m_rst}, {br(cell)}}; + case Dir::right: return {{br(cell), Dir::up, m_rst}, {tr(cell)}}; + case Dir::up: return {{tr(cell), Dir::left, m_rst}, {tl(cell)}}; + default: assert(false); } - } - - return {}; + return {}; } - - Edge edge(const Coord &ringvertex) const + + Edge edge(const Coord& ringvertex) const { - const long R = rows(*m_rst), C = cols(*m_rst); - const long R_1 = R - 1, C_1 = C - 1; - - Edge e = _edge(ringvertex); + Edge e = _edge(ringvertex); e.to.dir = e.from.dir; ++e.to; - - e.from.crd.r = std::min(e.from.crd.r, R_1); - e.from.crd.r = std::max(e.from.crd.r, 0l); - e.from.crd.c = std::min(e.from.crd.c, C_1); - e.from.crd.c = std::max(e.from.crd.c, 0l); - - e.to.crd.r = std::min(e.to.crd.r, R); - e.to.crd.r = std::max(e.to.crd.r, 0l); - e.to.crd.c = std::min(e.to.crd.c, C); - e.to.crd.c = std::max(e.to.crd.c, 0l); - + + e.from.crd.r = std::clamp(e.from.crd.r, 0l, m_rastsize.r - 1); + e.from.crd.c = std::clamp(e.from.crd.c, 0l, m_rastsize.c - 1); + e.to.crd.r = std::clamp(e.to.crd.r, 0l, m_rastsize.r); + e.to.crd.c = std::clamp(e.to.crd.c, 0l, m_rastsize.c); return e; } - -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{2 + (long(rows(rst)) - overlap.r) / m_window.r, - 2 + (long(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 - void tag_grid(ExecutionPolicy &&policy, TRasterValue isoval) - { - // parallel for r - for_each (std::forward(policy), - m_tags.begin(), m_tags.end(), - [this, isoval](uint8_t& tag, size_t idx) { - tag = get_tag_for_cell(coord(idx), isoval); - }); + + void interpolate_edge(Coord& ecrd, TRasterValue isoval) const + { + // The ecrd must have a grid index in r and a direction in c. + assert((0 <= ecrd.r) && (ecrd.r < m_gridlen)); + assert((ecrd.c == long(Dir::left)) || (ecrd.c == long(Dir::down)) || (ecrd.c == long(Dir::right)) || (ecrd.c == long(Dir::up))); + Edge e = edge(ecrd); + ecrd = std::lower_bound(e.from, e.to, isoval).crd; + // Shift bottom and right side points "out" by one to account for + // raster pixel width. Note "dir" is the direction of interpolation + // along the cell edge, not the next move direction. + if (e.from.dir == Dir::up) + ecrd.r += 1; + else if (e.from.dir == Dir::left) + ecrd.c += 1; } - - // 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. + +public: + explicit Grid(const Rst& rst, const Coord& window) + : m_rst{&rst} + , m_window{window} + , m_rastsize{static_cast(rows(rst)), static_cast(cols(rst))} + , m_gridsize{2 + m_rastsize.r / m_window.r, 2 + m_rastsize.c / m_window.c} + , m_gridlen{static_cast(m_gridsize.r * m_gridsize.c)} + , m_tags(m_gridlen / 32 + 1, 0) // 1 bit per tag means 32 per uint32_t. + , m_dirs(m_gridlen / 8 + 1, 0) // 4 bits per cell means 32/4=8 per uint32_t. + {} + + // Go through the cells getting their tags and dirs. + template void tag_grid(ExecutionPolicy&& policy, TRasterValue isoval) + { + // Get all the tags. parallel? + for_each_idx(std::forward(policy), m_tags.begin(), m_tags.end(), + [this, isoval](uint32_t& tag_block, size_t bidx) { tag_block = get_tags_block32(bidx, isoval); }); + // streamtags(std::cerr); + // Get all the dirs. parallel? + for_each_idx(std::forward(policy), m_dirs.begin(), m_dirs.end(), + [this](uint32_t& dirs_block, size_t bidx) { dirs_block = get_dirs_block8(bidx); }); + // streamdirs(std::cerr); + } + + // Scan for the rings on the tagged grid. Each ring vertex uses the Coord + // to store the sequential cell index (idx in r) and next direction (Dir + // in c) for the next point of the ring. This info can be used later to + // calculate the exact raster coordinate of the point. std::vector scan_rings() { std::vector rings; - size_t startidx = 0; - while ((startidx = search_start_cell(startidx)) < m_tags.size()) { + size_t startidx = 0; + while ((startidx = search_start_cell(startidx)) < m_gridlen) { 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)) { + + size_t idx = startidx; + Dir next = next_dir(idx); + do { + // We should never touch a cell with no remaining directions + // until we get back to the start cell. + assert(next != Dir::none); Coord ringvertex{long(idx), long(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(); + clr_dirs(idx, next); + idx = stepidx(idx, next); + next = next_dir(idx, next); + assert(idx < m_gridlen); + } while (idx != startidx); + // The start cell on returning should either have its directions + // cleared or have the second ambiguous direction. + assert(next == Dir::none || next == Dir::left || next == Dir::up); + if (ring.size() > 1) 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 - void interpolate_rings(ExecutionPolicy && policy, - std::vector &rings, - TRasterValue isov) + // sequential index of the square and the next direction + template void interpolate_rings(ExecutionPolicy&& policy, std::vector& rings, TRasterValue isov) { - for_each(std::forward(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; - } + for_each_idx(std::forward(policy), rings.begin(), rings.end(), [this, isov](Ring& ring, size_t) { + for (Coord& e : ring) + interpolate_edge(e, isov); }); } + + std::ostream& streamtags(std::ostream& os) + { + os << " :"; + for (auto c = 0; c < m_gridsize.c; c++) + os << (c % 10); + os << "\n"; + for (auto r = 0; r < m_gridsize.r; r++) { + os << std::setw(3) << r << ":"; + for (auto c = 0; c < m_gridsize.c; c++) + os << ((get_tags9(seq(Coord(r, c))) & 1) ? "H" : "."); + os << "\n"; + } + return os; + } + + std::ostream& streamdirs(std::ostream& os) + { + os << " :"; + for (auto c = 0; c < m_gridsize.c; c++) + os << (c % 10); + os << "\n"; + for (auto r = 0; r < m_gridsize.r; r++) { + os << std::setw(3) << r << ":"; + for (auto c = 0; c < m_gridsize.c; c++) + os << get_dirs(seq(Coord(r, c))); + os << std::endl; + } + return os; + } }; template -std::vector execute_with_policy(ExecutionPolicy && policy, - const Raster & raster, +std::vector execute_with_policy(ExecutionPolicy&& policy, + const Raster& raster, TRasterValue isoval, - Coord windowsize = {}) + Coord windowsize = {}) { - if (!rows(raster) || !cols(raster)) return {}; - + if (!rows(raster) || !cols(raster)) + return {}; + size_t ratio = cols(raster) / rows(raster); - - if (!windowsize.r) windowsize.r = 2; + + if (!windowsize.r) + windowsize.r = 2; if (!windowsize.c) windowsize.c = std::max(2l, long(windowsize.r * ratio)); - - Coord overlap{1}; - - Grid grid{raster, windowsize, overlap}; - + + Grid grid{raster, windowsize}; + grid.tag_grid(std::forward(policy), isoval); std::vector rings = grid.scan_rings(); grid.interpolate_rings(std::forward(policy), rings, isoval); - + return rings; } -template -std::vector execute(const Raster &raster, - TRasterValue isoval, - Coord windowsize = {}) +template std::vector execute(const Raster& raster, TRasterValue isoval, Coord windowsize = {}) { return execute_with_policy(nullptr, raster, isoval, windowsize); } diff --git a/src/libslic3r/SLA/RasterToPolygons.cpp b/src/libslic3r/SLA/RasterToPolygons.cpp index 79c21fe726..e05a721cd8 100644 --- a/src/libslic3r/SLA/RasterToPolygons.cpp +++ b/src/libslic3r/SLA/RasterToPolygons.cpp @@ -10,13 +10,13 @@ namespace marchsq { // Specialize this struct to register a raster type for the Marching squares alg template<> struct _RasterTraits { using Rst = Slic3r::sla::RasterGrayscaleAA; - + // The type of pixel cell in the raster using ValueType = uint8_t; - + // Value at a given position static uint8_t get(const Rst &rst, size_t row, size_t col) { return rst.read_pixel(col, row); } - + // Number of rows and cols of the raster static size_t rows(const Rst &rst) { return rst.resolution().height_px; } static size_t cols(const Rst &rst) { return rst.resolution().width_px; } @@ -34,57 +34,55 @@ template void foreach_vertex(ExPolygon &poly, Fn &&fn) } ExPolygons raster_to_polygons(const RasterGrayscaleAA &rst, Vec2i32 windowsize) -{ +{ size_t rows = rst.resolution().height_px, cols = rst.resolution().width_px; - + if (rows < 2 || cols < 2) return {}; - + Polygons polys; - long w_rows = std::max(2l, long(windowsize.y())); - long w_cols = std::max(2l, long(windowsize.x())); - + long w_rows = std::max(1l, long(windowsize.y())); + long w_cols = std::max(1l, long(windowsize.x())); + std::vector rings = marchsq::execute(rst, 128, {w_rows, w_cols}); - + polys.reserve(rings.size()); - + auto pxd = rst.pixel_dimensions(); - pxd.w_mm = (rst.resolution().width_px * pxd.w_mm) / (rst.resolution().width_px - 1); - pxd.h_mm = (rst.resolution().height_px * pxd.h_mm) / (rst.resolution().height_px - 1); - + for (const marchsq::Ring &ring : rings) { Polygon poly; Points &pts = poly.points; pts.reserve(ring.size()); - + for (const marchsq::Coord &crd : ring) pts.emplace_back(scaled(crd.c * pxd.w_mm), scaled(crd.r * pxd.h_mm)); - + polys.emplace_back(poly); } - + // reverse the raster transformations ExPolygons unioned = union_ex(polys); coord_t width = scaled(cols * pxd.h_mm), height = scaled(rows * pxd.w_mm); - + auto tr = rst.trafo(); for (ExPolygon &expoly : unioned) { if (tr.mirror_y) foreach_vertex(expoly, [height](Point &p) {p.y() = height - p.y(); }); - + if (tr.mirror_x) foreach_vertex(expoly, [width](Point &p) {p.x() = width - p.x(); }); - + expoly.translate(-tr.center_x, -tr.center_y); - + if (tr.flipXY) foreach_vertex(expoly, [](Point &p) { std::swap(p.x(), p.y()); }); - + if ((tr.mirror_x + tr.mirror_y + tr.flipXY) % 2) { expoly.contour.reverse(); for (auto &h : expoly.holes) h.reverse(); } } - + return unioned; } diff --git a/src/libslic3r/SLA/RasterToPolygons.hpp b/src/libslic3r/SLA/RasterToPolygons.hpp index e7e554fbb7..5320887c8f 100644 --- a/src/libslic3r/SLA/RasterToPolygons.hpp +++ b/src/libslic3r/SLA/RasterToPolygons.hpp @@ -8,7 +8,7 @@ namespace sla { class RasterGrayscaleAA; -ExPolygons raster_to_polygons(const RasterGrayscaleAA &rst, Vec2i32 windowsize = {2, 2}); +ExPolygons raster_to_polygons(const RasterGrayscaleAA &rst, Vec2i32 windowsize = {1, 1}); }} // namespace Slic3r::sla diff --git a/src/libslic3r/StreamUtils.hpp b/src/libslic3r/StreamUtils.hpp new file mode 100644 index 0000000000..78a9ac7335 --- /dev/null +++ b/src/libslic3r/StreamUtils.hpp @@ -0,0 +1,60 @@ +#ifndef slic3r_StreamUtils_hpp_ +#define slic3r_StreamUtils_hpp_ + +#include "Point.hpp" +#include "libslic3r.h" +#include "Polygon.hpp" +#include "Polyline.hpp" +#include "ExPolygon.hpp" +#include +#include + +namespace Slic3r { + +inline std::ostream& operator<<(std::ostream& os, const Points& pts) +{ + os << "[" << pts.size() << "]:"; + for (Point p : pts) + os << " (" << p << ")"; + os << "\n"; + return os; +} + +inline std::ostream& operator<<(std::ostream& os, const MultiPoint& mpts) +{ + os << "Multipoint" << mpts.points; + return os; +} + +inline std::ostream& operator<<(std::ostream& os, const Polygon& poly) +{ + os << "Polygon" << poly.points; + return os; +} + +inline std::ostream& operator<<(std::ostream& os, const Polygons& polys) +{ + os << "Polygons[" << polys.size() << "]:" << "\n"; + for (Polygon p : polys) + os << " " << p; + return os; +} + +inline std::ostream& operator<<(std::ostream& os, const ExPolygon& epoly) +{ + os << "ExPolygon:\n"; + os << " contour: " << epoly.contour; + os << " holes: " << epoly.holes; + return os; +} + +inline std::ostream& operator<<(std::ostream& os, const ExPolygons& epolys) +{ + os << "ExPolygons[" << epolys.size() << "]:" << "\n"; + for (ExPolygon p : epolys) + os << " " << p; + return os; +} + +} // namespace Slic3r +#endif // slic3r_StreamUtils_hpp_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 175b6e3a16..b97967fe34 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,6 +8,7 @@ add_library(Catch2 INTERFACE) list (APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake/modules/Catch2) target_include_directories(Catch2 INTERFACE ${CMAKE_CURRENT_LIST_DIR}) add_library(Catch2::Catch2 ALIAS Catch2) +target_compile_definitions(Catch2 INTERFACE -DCATCH_CONFIG_ENABLE_BENCHMARKING) if (APPLE) # OSX builds targeting OSX 10.9 do not support new std::uncought_exception() # see https://github.com/catchorg/Catch2/issues/1218 diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt index 2f71808cb7..ea01e27e99 100644 --- a/tests/libslic3r/CMakeLists.txt +++ b/tests/libslic3r/CMakeLists.txt @@ -1,44 +1,44 @@ get_filename_component(_TEST_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) -add_executable(${_TEST_NAME}_tests - ${_TEST_NAME}_tests.cpp - test_3mf.cpp - test_aabbindirect.cpp - test_clipper_offset.cpp - test_clipper_utils.cpp - test_config.cpp - test_elephant_foot_compensation.cpp - test_geometry.cpp - test_placeholder_parser.cpp - test_polygon.cpp - test_mutable_polygon.cpp - test_mutable_priority_queue.cpp - test_stl.cpp - test_meshboolean.cpp - # test_marchingsquares.cpp - test_timeutils.cpp - test_voronoi.cpp +add_executable(${_TEST_NAME}_tests + ${_TEST_NAME}_tests.cpp + test_3mf.cpp + test_aabbindirect.cpp + test_clipper_offset.cpp + test_clipper_utils.cpp + test_config.cpp + test_elephant_foot_compensation.cpp + test_geometry.cpp + test_placeholder_parser.cpp + test_polygon.cpp + test_mutable_polygon.cpp + test_mutable_priority_queue.cpp + test_stl.cpp + test_meshboolean.cpp + test_marchingsquares.cpp + test_timeutils.cpp + test_voronoi.cpp test_optimizers.cpp # test_png_io.cpp test_indexed_triangle_set.cpp ../libnest2d/printer_parts.cpp - ) + ) if (TARGET OpenVDB::openvdb) target_sources(${_TEST_NAME}_tests PRIVATE test_hollowing.cpp) endif() - + target_link_libraries(${_TEST_NAME}_tests test_common libslic3r) set_property(TARGET ${_TEST_NAME}_tests PROPERTY FOLDER "tests") if (WIN32) - if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") - orcaslicer_copy_dlls(COPY_DLLS "Debug" "d" output_dlls_Debug) - elseif("${CMAKE_BUILD_TYPE}" STREQUAL "RelWithDebInfo") - orcaslicer_copy_dlls(COPY_DLLS "RelWithDebInfo" "" output_dlls_Release) - else() - orcaslicer_copy_dlls(COPY_DLLS "Release" "" output_dlls_Release) - endif() + if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") + orcaslicer_copy_dlls(COPY_DLLS "Debug" "d" output_dlls_Debug) + elseif("${CMAKE_BUILD_TYPE}" STREQUAL "RelWithDebInfo") + orcaslicer_copy_dlls(COPY_DLLS "RelWithDebInfo" "" output_dlls_Release) + else() + orcaslicer_copy_dlls(COPY_DLLS "Release" "" output_dlls_Release) + endif() endif() catch_discover_tests(${_TEST_NAME}_tests TEST_PREFIX "${_TEST_NAME}: ") diff --git a/tests/libslic3r/test_clipper_utils.cpp b/tests/libslic3r/test_clipper_utils.cpp index 86afdcbcc4..6bbad8ef72 100644 --- a/tests/libslic3r/test_clipper_utils.cpp +++ b/tests/libslic3r/test_clipper_utils.cpp @@ -11,16 +11,16 @@ using namespace Slic3r; SCENARIO("Various Clipper operations - xs/t/11_clipper.t", "[ClipperUtils]") { - // CCW oriented contour - Slic3r::Polygon square{ { 200, 100 }, {200, 200}, {100, 200}, {100, 100} }; - // CW oriented contour - Slic3r::Polygon hole_in_square{ { 160, 140 }, { 140, 140 }, { 140, 160 }, { 160, 160 } }; - Slic3r::ExPolygon square_with_hole(square, hole_in_square); - GIVEN("square_with_hole") { + // CCW oriented contour + Slic3r::Polygon square{ { 200, 100 }, {200, 200}, {100, 200}, {100, 100} }; + // CW oriented contour + Slic3r::Polygon hole_in_square{ { 160, 140 }, { 140, 140 }, { 140, 160 }, { 160, 160 } }; + Slic3r::ExPolygon square_with_hole(square, hole_in_square); + GIVEN("square_with_hole") { WHEN("offset") { Polygons result = Slic3r::offset(square_with_hole, 5.f); THEN("offset matches") { - REQUIRE(result == Polygons { + REQUIRE(result == Polygons { { { 205, 205 }, { 95, 205 }, { 95, 95 }, { 205, 95 }, }, { { 155, 145 }, { 145, 145 }, { 145, 155 }, { 155, 155 } } }); } @@ -28,7 +28,7 @@ SCENARIO("Various Clipper operations - xs/t/11_clipper.t", "[ClipperUtils]") { WHEN("offset_ex") { ExPolygons result = Slic3r::offset_ex(square_with_hole, 5.f); THEN("offset matches") { - REQUIRE(result == ExPolygons { { + REQUIRE(result == ExPolygons { { { { 205, 205 }, { 95, 205 }, { 95, 95 }, { 205, 95 }, }, { { 145, 145 }, { 145, 155 }, { 155, 155 }, { 155, 145 } } } } ); } @@ -66,7 +66,7 @@ SCENARIO("Various Clipper operations - xs/t/11_clipper.t", "[ClipperUtils]") { GIVEN("polyline") { Polyline polyline { { 50, 150 }, { 300, 150 } }; WHEN("intersection_pl") { - Polylines result = Slic3r::intersection_pl({ polyline }, { square, hole_in_square }); + Polylines result = Slic3r::intersection_pl({ polyline }, Polygons{ square, hole_in_square }); THEN("correct number of result lines") { REQUIRE(result.size() == 2); } @@ -93,82 +93,82 @@ SCENARIO("Various Clipper operations - xs/t/11_clipper.t", "[ClipperUtils]") { } } } - GIVEN("Clipper bug #96 / Slic3r issue #2028") { - Slic3r::Polyline subject{ - { 44735000, 31936670 }, { 55270000, 31936670 }, { 55270000, 25270000 }, { 74730000, 25270000 }, { 74730000, 44730000 }, { 68063296, 44730000 }, { 68063296, 55270000 }, { 74730000, 55270000 }, - { 74730000, 74730000 }, { 55270000, 74730000 }, { 55270000, 68063296 }, { 44730000, 68063296 }, { 44730000, 74730000 }, { 25270000, 74730000 }, { 25270000, 55270000 }, { 31936670, 55270000 }, - { 31936670, 44730000 }, { 25270000, 44730000 }, { 25270000, 25270000 }, { 44730000, 25270000 }, { 44730000, 31936670 } }; - Slic3r::Polygon clip { {75200000, 45200000}, {54800000, 45200000}, {54800000, 24800000}, {75200000, 24800000} }; - Slic3r::Polylines result = Slic3r::intersection_pl(subject, { clip }); - THEN("intersection_pl - result is not empty") { - REQUIRE(result.size() == 1); - } - } - GIVEN("Clipper bug #122") { - Slic3r::Polyline subject { { 1975, 1975 }, { 25, 1975 }, { 25, 25 }, { 1975, 25 }, { 1975, 1975 } }; - Slic3r::Polygons clip { { { 2025, 2025 }, { -25, 2025 } , { -25, -25 }, { 2025, -25 } }, - { { 525, 525 }, { 525, 1475 }, { 1475, 1475 }, { 1475, 525 } } }; - Slic3r::Polylines result = Slic3r::intersection_pl({ subject }, clip); - THEN("intersection_pl - result is not empty") { - REQUIRE(result.size() == 1); - REQUIRE(result.front().points.size() == 5); - } - } - GIVEN("Clipper bug #126") { - Slic3r::Polyline subject { { 200000, 19799999 }, { 200000, 200000 }, { 24304692, 200000 }, { 15102879, 17506106 }, { 13883200, 19799999 }, { 200000, 19799999 } }; - Slic3r::Polygon clip { { 15257205, 18493894 }, { 14350057, 20200000 }, { -200000, 20200000 }, { -200000, -200000 }, { 25196917, -200000 } }; - Slic3r::Polylines result = Slic3r::intersection_pl(subject, { clip }); - THEN("intersection_pl - result is not empty") { - REQUIRE(result.size() == 1); - } - THEN("intersection_pl - result has same length as subject polyline") { - REQUIRE(result.front().length() == Approx(subject.length())); - } - } + GIVEN("Clipper bug #96 / Slic3r issue #2028") { + Slic3r::Polyline subject{ + { 44735000, 31936670 }, { 55270000, 31936670 }, { 55270000, 25270000 }, { 74730000, 25270000 }, { 74730000, 44730000 }, { 68063296, 44730000 }, { 68063296, 55270000 }, { 74730000, 55270000 }, + { 74730000, 74730000 }, { 55270000, 74730000 }, { 55270000, 68063296 }, { 44730000, 68063296 }, { 44730000, 74730000 }, { 25270000, 74730000 }, { 25270000, 55270000 }, { 31936670, 55270000 }, + { 31936670, 44730000 }, { 25270000, 44730000 }, { 25270000, 25270000 }, { 44730000, 25270000 }, { 44730000, 31936670 } }; + Slic3r::Polygon clip { {75200000, 45200000}, {54800000, 45200000}, {54800000, 24800000}, {75200000, 24800000} }; + Slic3r::Polylines result = Slic3r::intersection_pl(subject, Polygons{ clip }); + THEN("intersection_pl - result is not empty") { + REQUIRE(result.size() == 1); + } + } + GIVEN("Clipper bug #122") { + Slic3r::Polyline subject { { 1975, 1975 }, { 25, 1975 }, { 25, 25 }, { 1975, 25 }, { 1975, 1975 } }; + Slic3r::Polygons clip { { { 2025, 2025 }, { -25, 2025 } , { -25, -25 }, { 2025, -25 } }, + { { 525, 525 }, { 525, 1475 }, { 1475, 1475 }, { 1475, 525 } } }; + Slic3r::Polylines result = Slic3r::intersection_pl({ subject }, clip); + THEN("intersection_pl - result is not empty") { + REQUIRE(result.size() == 1); + REQUIRE(result.front().points.size() == 5); + } + } + GIVEN("Clipper bug #126") { + Slic3r::Polyline subject { { 200000, 19799999 }, { 200000, 200000 }, { 24304692, 200000 }, { 15102879, 17506106 }, { 13883200, 19799999 }, { 200000, 19799999 } }; + Slic3r::Polygon clip { { 15257205, 18493894 }, { 14350057, 20200000 }, { -200000, 20200000 }, { -200000, -200000 }, { 25196917, -200000 } }; + Slic3r::Polylines result = Slic3r::intersection_pl(subject, Polygons{ clip }); + THEN("intersection_pl - result is not empty") { + REQUIRE(result.size() == 1); + } + THEN("intersection_pl - result has same length as subject polyline") { + REQUIRE(result.front().length() == Approx(subject.length())); + } + } #if 0 - { - # Clipper does not preserve polyline orientation - my $polyline = Slic3r::Polyline->new([50, 150], [300, 150]); - my $result = Slic3r::Geometry::Clipper::intersection_pl([$polyline], [$square]); - is scalar(@$result), 1, 'intersection_pl - correct number of result lines'; - is_deeply $result->[0]->pp, [[100, 150], [200, 150]], 'clipped line orientation is preserved'; - } - { - # Clipper does not preserve polyline orientation - my $polyline = Slic3r::Polyline->new([300, 150], [50, 150]); - my $result = Slic3r::Geometry::Clipper::intersection_pl([$polyline], [$square]); - is scalar(@$result), 1, 'intersection_pl - correct number of result lines'; - is_deeply $result->[0]->pp, [[200, 150], [100, 150]], 'clipped line orientation is preserved'; - } - { - # Disabled until Clipper bug #127 is fixed - my $subject = [ - Slic3r::Polyline->new([-90000000, -100000000], [-90000000, 100000000]), # vertical - Slic3r::Polyline->new([-100000000, -10000000], [100000000, -10000000]), # horizontal - Slic3r::Polyline->new([-100000000, 0], [100000000, 0]), # horizontal - Slic3r::Polyline->new([-100000000, 10000000], [100000000, 10000000]), # horizontal - ]; - my $clip = Slic3r::Polygon->new(# a circular, convex, polygon - [99452190, 10452846], [97814760, 20791169], [95105652, 30901699], [91354546, 40673664], [86602540, 50000000], - [80901699, 58778525], [74314483, 66913061], [66913061, 74314483], [58778525, 80901699], [50000000, 86602540], - [40673664, 91354546], [30901699, 95105652], [20791169, 97814760], [10452846, 99452190], [0, 100000000], - [-10452846, 99452190], [-20791169, 97814760], [-30901699, 95105652], [-40673664, 91354546], - [-50000000, 86602540], [-58778525, 80901699], [-66913061, 74314483], [-74314483, 66913061], - [-80901699, 58778525], [-86602540, 50000000], [-91354546, 40673664], [-95105652, 30901699], - [-97814760, 20791169], [-99452190, 10452846], [-100000000, 0], [-99452190, -10452846], - [-97814760, -20791169], [-95105652, -30901699], [-91354546, -40673664], [-86602540, -50000000], - [-80901699, -58778525], [-74314483, -66913061], [-66913061, -74314483], [-58778525, -80901699], - [-50000000, -86602540], [-40673664, -91354546], [-30901699, -95105652], [-20791169, -97814760], - [-10452846, -99452190], [0, -100000000], [10452846, -99452190], [20791169, -97814760], - [30901699, -95105652], [40673664, -91354546], [50000000, -86602540], [58778525, -80901699], - [66913061, -74314483], [74314483, -66913061], [80901699, -58778525], [86602540, -50000000], - [91354546, -40673664], [95105652, -30901699], [97814760, -20791169], [99452190, -10452846], [100000000, 0] - ); - my $result = Slic3r::Geometry::Clipper::intersection_pl($subject, [$clip]); - is scalar(@$result), scalar(@$subject), 'intersection_pl - expected number of polylines'; - is sum(map scalar(@$_), @$result), scalar(@$subject) * 2, 'intersection_pl - expected number of points in polylines'; - } + { + # Clipper does not preserve polyline orientation + my $polyline = Slic3r::Polyline->new([50, 150], [300, 150]); + my $result = Slic3r::Geometry::Clipper::intersection_pl([$polyline], [$square]); + is scalar(@$result), 1, 'intersection_pl - correct number of result lines'; + is_deeply $result->[0]->pp, [[100, 150], [200, 150]], 'clipped line orientation is preserved'; + } + { + # Clipper does not preserve polyline orientation + my $polyline = Slic3r::Polyline->new([300, 150], [50, 150]); + my $result = Slic3r::Geometry::Clipper::intersection_pl([$polyline], [$square]); + is scalar(@$result), 1, 'intersection_pl - correct number of result lines'; + is_deeply $result->[0]->pp, [[200, 150], [100, 150]], 'clipped line orientation is preserved'; + } + { + # Disabled until Clipper bug #127 is fixed + my $subject = [ + Slic3r::Polyline->new([-90000000, -100000000], [-90000000, 100000000]), # vertical + Slic3r::Polyline->new([-100000000, -10000000], [100000000, -10000000]), # horizontal + Slic3r::Polyline->new([-100000000, 0], [100000000, 0]), # horizontal + Slic3r::Polyline->new([-100000000, 10000000], [100000000, 10000000]), # horizontal + ]; + my $clip = Slic3r::Polygon->new(# a circular, convex, polygon + [99452190, 10452846], [97814760, 20791169], [95105652, 30901699], [91354546, 40673664], [86602540, 50000000], + [80901699, 58778525], [74314483, 66913061], [66913061, 74314483], [58778525, 80901699], [50000000, 86602540], + [40673664, 91354546], [30901699, 95105652], [20791169, 97814760], [10452846, 99452190], [0, 100000000], + [-10452846, 99452190], [-20791169, 97814760], [-30901699, 95105652], [-40673664, 91354546], + [-50000000, 86602540], [-58778525, 80901699], [-66913061, 74314483], [-74314483, 66913061], + [-80901699, 58778525], [-86602540, 50000000], [-91354546, 40673664], [-95105652, 30901699], + [-97814760, 20791169], [-99452190, 10452846], [-100000000, 0], [-99452190, -10452846], + [-97814760, -20791169], [-95105652, -30901699], [-91354546, -40673664], [-86602540, -50000000], + [-80901699, -58778525], [-74314483, -66913061], [-66913061, -74314483], [-58778525, -80901699], + [-50000000, -86602540], [-40673664, -91354546], [-30901699, -95105652], [-20791169, -97814760], + [-10452846, -99452190], [0, -100000000], [10452846, -99452190], [20791169, -97814760], + [30901699, -95105652], [40673664, -91354546], [50000000, -86602540], [58778525, -80901699], + [66913061, -74314483], [74314483, -66913061], [80901699, -58778525], [86602540, -50000000], + [91354546, -40673664], [95105652, -30901699], [97814760, -20791169], [99452190, -10452846], [100000000, 0] + ); + my $result = Slic3r::Geometry::Clipper::intersection_pl($subject, [$clip]); + is scalar(@$result), scalar(@$subject), 'intersection_pl - expected number of polylines'; + is sum(map scalar(@$_), @$result), scalar(@$subject) * 2, 'intersection_pl - expected number of points in polylines'; + } #endif } @@ -203,7 +203,7 @@ SCENARIO("Various Clipper operations - t/clipper.t", "[ClipperUtils]") { } } WHEN("diff_ex with another square") { - ExPolygons diff = Slic3r::diff_ex(Polygons{ square, square2 }, Polygons{ hole }); + ExPolygons diff = Slic3r::diff_ex(Polygons{ square, square2 }, Polygons{ hole }); THEN("difference of a cw from two ccw is a contour with one hole") { REQUIRE(diff.size() == 1); REQUIRE(diff.front().area() == Approx(ExPolygon({ {40, 40}, {0, 40}, {0, 0}, {40, 0} }, { {15, 25}, {25, 25}, {25, 15}, {15, 15} }).area())); @@ -225,11 +225,11 @@ SCENARIO("Various Clipper operations - t/clipper.t", "[ClipperUtils]") { } } -template +template double polytree_area(const Tree &tree, std::vector *out) { traverse_pt(tree, out); - + return std::accumulate(out->begin(), out->end(), 0.0, [](double a, const P &p) { return a + p.area(); }); } @@ -238,7 +238,7 @@ size_t count_polys(const ExPolygons& expolys) { size_t c = 0; for (auto &ep : expolys) c += ep.holes.size() + 1; - + return c; } @@ -247,32 +247,32 @@ TEST_CASE("Traversing Clipper PolyTree", "[ClipperUtils]") { Polygon unitbox; const auto UNIT = coord_t(1. / SCALING_FACTOR); unitbox.points = { Vec2crd{0, 0}, Vec2crd{UNIT, 0}, Vec2crd{UNIT, UNIT}, Vec2crd{0, UNIT}}; - + Polygon box_frame = unitbox; box_frame.scale(20, 10); - + Polygon hole_left = unitbox; hole_left.scale(8); hole_left.translate(UNIT, UNIT); hole_left.reverse(); - + Polygon hole_right = hole_left; hole_right.translate(UNIT * 10, 0); - + Polygon inner_left = unitbox; inner_left.scale(4); inner_left.translate(UNIT * 3, UNIT * 3); - + Polygon inner_right = inner_left; inner_right.translate(UNIT * 10, 0); - + Polygons reference = union_({box_frame, hole_left, hole_right, inner_left, inner_right}); - + ClipperLib::PolyTree tree = union_pt(reference); double area_sum = box_frame.area() + hole_left.area() + hole_right.area() + inner_left.area() + inner_right.area(); - + REQUIRE(area_sum > 0); SECTION("Traverse into Polygons WITHOUT spatial ordering") { @@ -280,19 +280,19 @@ TEST_CASE("Traversing Clipper PolyTree", "[ClipperUtils]") { REQUIRE(area_sum == Approx(polytree_area(tree.GetFirst(), &output))); REQUIRE(output.size() == reference.size()); } - + SECTION("Traverse into ExPolygons WITHOUT spatial ordering") { ExPolygons output; REQUIRE(area_sum == Approx(polytree_area(tree.GetFirst(), &output))); REQUIRE(count_polys(output) == reference.size()); } - + SECTION("Traverse into Polygons WITH spatial ordering") { Polygons output; REQUIRE(area_sum == Approx(polytree_area(tree.GetFirst(), &output))); REQUIRE(output.size() == reference.size()); } - + SECTION("Traverse into ExPolygons WITH spatial ordering") { ExPolygons output; REQUIRE(area_sum == Approx(polytree_area(tree.GetFirst(), &output))); diff --git a/tests/libslic3r/test_marchingsquares.cpp b/tests/libslic3r/test_marchingsquares.cpp index fed5cc2d4d..db427a8d70 100644 --- a/tests/libslic3r/test_marchingsquares.cpp +++ b/tests/libslic3r/test_marchingsquares.cpp @@ -17,22 +17,111 @@ #include #include #include +#include using namespace Slic3r; +using namespace Catch::Matchers; -static double area(const sla::RasterBase::PixelDim &pxd) +// Note this tests SLA/RasterToPolygons.hpp, SLA/AGGRaster.hpp, and +// ClipperUtils.hpp at least as much as MarchingSquares.hpp. + +// Get the Point corresponding to a raster column and row. +Point rstPoint(const sla::RasterGrayscaleAA& rst, const size_t c, const size_t r) { - return pxd.w_mm * pxd.h_mm; + size_t rows = rst.resolution().height_px, cols = rst.resolution().width_px; + auto pxd = rst.pixel_dimensions(); + auto tr = rst.trafo(); + coord_t width = scaled(cols * pxd.h_mm), height = scaled(rows * pxd.w_mm); + Point p = Point::new_scale(c * pxd.w_mm, r * pxd.h_mm); + // reverse the raster transformations + if (tr.mirror_y) + p.y() = height - p.y(); + if (tr.mirror_x) + p.x() = width - p.x(); + p.x() -= tr.center_x; + p.y() -= tr.center_y; + if (tr.flipXY) + std::swap(p.x(), p.y()); + return p; } -static Slic3r::sla::RasterGrayscaleAA create_raster( - const sla::RasterBase::Resolution &res, - double disp_w = 100., - double disp_h = 100.) +// Get the size of a raster pixel in coord_t. +static Point rstPixel(const sla::RasterGrayscaleAA& rst) { - sla::RasterBase::PixelDim pixdim{disp_w / res.width_px, disp_h / res.height_px}; - - auto bb = BoundingBox({0, 0}, {scaled(disp_w), scaled(disp_h)}); + auto pxd = rst.pixel_dimensions(); + return Point::new_scale(pxd.w_mm, pxd.h_mm); +} + +// Get the size of a raster in coord_t. +static Point rstSize(const sla::RasterGrayscaleAA& rst) +{ + auto pxd = rst.pixel_dimensions(); + auto res = rst.resolution(); + return Point::new_scale(pxd.w_mm * res.width_px, pxd.h_mm * res.height_px); +} + +// Get the bounding box of a raster. +static BoundingBox rstBBox(const sla::RasterGrayscaleAA& rst) +{ + auto center = rst.trafo().get_center(); + return BoundingBox(Point(0, 0) - center, rstSize(rst) - center); +} + +// Get the ExPolygons directly corresponding to a raster. +static ExPolygons rstGetPolys(sla::RasterGrayscaleAA& rst) +{ + size_t rows = rst.resolution().height_px, cols = rst.resolution().width_px; + Polygons polys; + for (auto r = 0; r < rows; r++) { + // use c0==cols as a sentinel marker for "no start column yet". + size_t c0 = cols; + for (auto c = 0; c <= cols; c++) { + if (c < cols && rst.read_pixel(c, r) > 128) { + // We have set pixels, set the c0 start column if it is not yet set. + if (c0 == cols) + c0 = c; + } else if (c0 < cols) { + // There is no pixel set, but we do have a c0 start column. Output a + // "row-rectangle" poly for this row between the start column c0 and + // the current column. + polys.push_back({rstPoint(rst, c0, r), rstPoint(rst, c0, r + 1), rstPoint(rst, c, r + 1), rstPoint(rst, c, r)}); + // Make sure the poly is anti-clockwise, which it might not be + // depending on how rstPoint() reverses the raster transformations + // from (c,r) to (x,y) coordinates. + if (polys.back().is_clockwise()) + polys.back().reverse(); + // Clear the start column c0 for the next row-rectangle. + c0 = cols; + } + } + } + // Merge all the row-rectangle polys into contiguous raster ExPolygons. + return union_ex(polys); +} + +// Get the length in mm of a "vector" Point. +static double len(const Point& v) { return unscaled(v.norm()); } +// Get the area in mm^2 of a box with corners at the origin and a Point. +static double area(const Point& v) { return unscaled(v.x()) * unscaled(v.y()); } + +// Find the index of the nearest extracted ExPolygon for a reference ExPolygon. +static int find_closest_ext(const ExPolygons& exts, ExPolygon ref) +{ + auto ref_center = ref.contour.bounding_box().center(); + + auto closest = std::min_element(exts.begin(), exts.end(), [&ref_center](auto a, auto b) { + auto a_center = a.contour.bounding_box().center(); + auto b_center = b.contour.bounding_box().center(); + return a_center.distance_to(ref_center) < b_center.distance_to(ref_center); + }); + return std::distance(exts.begin(), closest); +} + +static Slic3r::sla::RasterGrayscaleAA create_raster(const sla::Resolution& res, double disp_w = 100., double disp_h = 100.) +{ + sla::PixelDim pixdim{disp_w / res.width_px, disp_h / res.height_px}; + + auto bb = BoundingBox({0, 0}, {scaled(disp_w), scaled(disp_h)}); sla::RasterBase::Trafo trafo; trafo.center_x = bb.center().x(); trafo.center_y = bb.center().y(); @@ -43,334 +132,439 @@ static Slic3r::sla::RasterGrayscaleAA create_raster( static ExPolygon square(double a, Point center = {0, 0}) { ExPolygon poly; - coord_t V = scaled(a / 2.); - + coord_t V = scaled(a / 2.); + poly.contour.points = {{-V, -V}, {V, -V}, {V, V}, {-V, V}}; poly.translate(center.x(), center.y()); - + return poly; } static ExPolygon square_with_hole(double a, Point center = {0, 0}) { ExPolygon poly = square(a); - + poly.holes.emplace_back(); - coord_t V = scaled(a / 4.); + coord_t V = scaled(a / 4.); poly.holes.front().points = {{-V, V}, {V, V}, {V, -V}, {-V, -V}}; - + poly.translate(center.x(), center.y()); return poly; } -static ExPolygons circle_with_hole(double r, Point center = {0, 0}) { - +static ExPolygons circle_with_hole(double r, Point center = {0, 0}) +{ ExPolygon poly; - + std::vector pis = linspace_vector(0., 2 * PI, 100); - + coord_t rs = scaled(r); for (double phi : pis) { poly.contour.points.emplace_back(rs * std::cos(phi), rs * std::sin(phi)); } - + poly.holes.emplace_back(poly.contour); poly.holes.front().reverse(); - for (auto &p : poly.holes.front().points) p /= 2; - + for (auto& p : poly.holes.front().points) + p /= 2; + poly.translate(center.x(), center.y()); - + return {poly}; } -static const Vec2i32 W4x4 = {4, 4}; static const Vec2i32 W2x2 = {2, 2}; +static const Vec2i32 W1x1 = {1, 1}; template -static void test_expolys(Rst && rst, - const ExPolygons & ref, - Vec2i32 window, - const std::string &name = "test") +static void test_expolys(Rst&& rst, const ExPolygons& ref, Vec2i32 window, const std::string& name = "test", bool strict = true) { - for (const ExPolygon &expoly : ref) rst.draw(expoly); - + auto raster_bb = rstBBox(rst); + Point pixel_size = rstPixel(rst); + Point window_size{coord_t(pixel_size.x() * window.x()), coord_t(pixel_size.y() * window.y())}; + double pixel_area = area(pixel_size); + double pixel_len = len(pixel_size); + double window_area = area(window_size); + double window_len = len(window_size); + + for (const ExPolygon& expoly : ref) + rst.draw(expoly); + std::fstream out(name + ".png", std::ios::out); out << rst.encode(sla::PNGRasterEncoder{}); out.close(); - - ExPolygons extracted = sla::raster_to_polygons(rst, window); - - SVG svg(name + ".svg"); - svg.draw(extracted); - svg.draw(ref, "green"); + + const ExPolygons bmp = rstGetPolys(rst); + const ExPolygons ext = sla::raster_to_polygons(rst, window); + + SVG svg(name + ".svg", raster_bb); + svg.draw(bmp, "green"); + if (pixel_size.x() >= scale_(0.5)) + svg.draw_grid(raster_bb, "grey", scale_(0.05), pixel_size.x()); + if (window_size.x() >= scale_(1.0)) + svg.draw_grid(raster_bb, "grey", scale_(0.10), window_size.x()); + svg.draw_outline(ref, "red", "red", scale_(0.3)); + svg.draw_outline(ext, "blue", "blue"); svg.Close(); - - double max_rel_err = 0.1; - sla::RasterBase::PixelDim pxd = rst.pixel_dimensions(); - double max_abs_err = area(pxd) * scaled(1.) * scaled(1.); - - BoundingBox ref_bb; - for (auto &expoly : ref) ref_bb.merge(expoly.contour.bounding_box()); - - double max_displacement = 4. * (std::pow(pxd.h_mm, 2) + std::pow(pxd.w_mm, 2)); - max_displacement *= scaled(1.) * scaled(1.); - - REQUIRE(extracted.size() == ref.size()); + + // Note all these areas are unscaled back to mm^2. + double raster_area = unscaled(unscaled(area(bmp))); + double reference_area = unscaled(unscaled(area(ref))); + double extracted_area = unscaled(unscaled(area(ext))); + + // Note that errors accumulate with each step going from the reference + // polys to the extracted polys. The rendering of the reference polys to + // the raster does introduce pixelization errors too. This checks for + // acceptable errors going from reference to raster, and raster to + // reference. for (size_t i = 0; i < ref.size(); ++i) { - REQUIRE(extracted[i].contour.is_counter_clockwise()); - REQUIRE(extracted[i].holes.size() == ref[i].holes.size()); - - for (auto &h : extracted[i].holes) REQUIRE(h.is_clockwise()); - - double refa = ref[i].area(); - double abs_err = std::abs(extracted[i].area() - refa); - double rel_err = abs_err / refa; - - REQUIRE((rel_err <= max_rel_err || abs_err <= max_abs_err)); - - BoundingBox bb; - for (auto &expoly : extracted) bb.merge(expoly.contour.bounding_box()); - - Point d = bb.center() - ref_bb.center(); - REQUIRE(double(d.transpose() * d) <= max_displacement); + if (ref[i].contour.size() < 20) + UNSCOPED_INFO("reference ref[" << i << "]: " << ref[i]); + } + CHECK_THAT(raster_area, WithinRel(reference_area, pixel_len * 0.05) || WithinAbs(reference_area, pixel_area)); + for (size_t i = 0; i < ext.size(); ++i) { + if (ext[i].contour.size() < 20) + UNSCOPED_INFO("extracted ext[" << i << "]: " << ext[i]); + } + CHECK_THAT(extracted_area, WithinRel(raster_area, 0.05) || WithinAbs(raster_area, window_area)); + for (auto i = 0; i < ext.size(); ++i) { + CHECK(ext[i].contour.is_counter_clockwise()); + for (auto& h : ext[i].holes) + CHECK(h.is_clockwise()); + } + + BoundingBox ref_bb; + for (auto& expoly : ref) + ref_bb.merge(expoly.contour.bounding_box()); + BoundingBox ext_bb; + for (auto& expoly : ext) + ext_bb.merge(expoly.contour.bounding_box()); + CHECK(len(ext_bb.center() - ref_bb.center()) < pixel_len); + + // In ambigous cases (when polygons just touch) there are multiple equally + // valid interpretations of the raster into polygons. Although + // MarchingSquares currently systematically selects the solution that + // breaks them into separate polygons, that might not always be true. Also, + // SLA/RasterToPolygons.hpp, and in particular union_ex() from + // ClipperUtils.hpp that it uses, can and does sometimes merge them back + // together. This means we cannot reliably make assertions about the + // extracted number of polygons and their shapes in these cases. So we skip + // the individual polygon checks for strict=false. + if (strict) { + CHECK(ext.size() == ref.size()); + for (auto i = 0; i < ext.size(); ++i) { + auto j = find_closest_ext(ref, ext[i]); + INFO("Comparing ext[" << i << "] against closest ref[" << j << "]"); + CHECK(ext[i].holes.size() == ref[j].holes.size()); + double ext_i_area = unscaled(unscaled(ext[i].area())); + double ref_j_area = unscaled(unscaled(ref[j].area())); + CHECK_THAT(ext_i_area, WithinRel(ref_j_area, pixel_len * 0.05) || WithinAbs(ref_j_area, window_area)); + auto ext_i_bb = ext[i].contour.bounding_box(); + auto ref_j_bb = ref[j].contour.bounding_box(); + CHECK(len(ext_i_bb.center() - ref_j_bb.center()) < pixel_len); + } } } -TEST_CASE("Empty raster should result in empty polygons", "[MarchingSquares]") { +TEST_CASE("Empty raster should result in empty polygons", "[MarchingSquares]") +{ sla::RasterGrayscaleAAGammaPower rst{{}, {}, {}}; - ExPolygons extracted = sla::raster_to_polygons(rst); + ExPolygons extracted = sla::raster_to_polygons(rst); REQUIRE(extracted.size() == 0); } -TEST_CASE("Marching squares directions", "[MarchingSquares]") { - marchsq::Coord crd{1, 1}; - - REQUIRE(step(crd, marchsq::__impl::Dir::left).r == 1); - REQUIRE(step(crd, marchsq::__impl::Dir::left).c == 0); - - REQUIRE(step(crd, marchsq::__impl::Dir::down).r == 2); - REQUIRE(step(crd, marchsq::__impl::Dir::down).c == 1); - - REQUIRE(step(crd, marchsq::__impl::Dir::right).r == 1); - REQUIRE(step(crd, marchsq::__impl::Dir::right).c == 2); - - REQUIRE(step(crd, marchsq::__impl::Dir::up).r == 0); - REQUIRE(step(crd, marchsq::__impl::Dir::up).c == 1); +TEST_CASE("Marching squares directions", "[MarchingSquares]") +{ + using namespace marchsq; + Coord crd{0, 0}; + + __impl::step(crd, __impl::Dir::left); + CHECK(crd == Coord(0, -1)); + __impl::step(crd, __impl::Dir::down); + CHECK(crd == Coord(1, -1)); + __impl::step(crd, __impl::Dir::right); + CHECK(crd == Coord(1, 0)); + __impl::step(crd, __impl::Dir::up); + CHECK(crd == Coord(0, 0)); + __impl::step(crd, __impl::Dir::left, 7); + CHECK(crd == Coord(0, -7)); + __impl::step(crd, __impl::Dir::down, 7); + CHECK(crd == Coord(7, -7)); + __impl::step(crd, __impl::Dir::right, 7); + CHECK(crd == Coord(7, 0)); + __impl::step(crd, __impl::Dir::up, 7); + CHECK(crd == Coord(0, 0)); + __impl::step(crd, __impl::Dir::left, -3); + CHECK(crd == Coord(0, 3)); + __impl::step(crd, __impl::Dir::down, -3); + CHECK(crd == Coord(-3, 3)); + __impl::step(crd, __impl::Dir::right, -3); + CHECK(crd == Coord(-3, 0)); + __impl::step(crd, __impl::Dir::up, -3); + CHECK(crd == Coord(0, 0)); } -TEST_CASE("Fully covered raster should result in a rectangle", "[MarchingSquares]") { +TEST_CASE("Fully covered raster should result in a rectangle", "[MarchingSquares]") +{ auto rst = create_raster({4, 4}, 4., 4.); - + ExPolygon rect = square(4); - SECTION("Full accuracy") { - test_expolys(rst, {rect}, W2x2, "fully_covered_full_acc"); - } - - SECTION("Half accuracy") { - test_expolys(rst, {rect}, W4x4, "fully_covered_half_acc"); - } + SECTION("Full accuracy") { test_expolys(rst, {rect}, W1x1, "fully_covered_full_acc"); } + + SECTION("Half accuracy") { test_expolys(rst, {rect}, W2x2, "fully_covered_half_acc"); } } -TEST_CASE("4x4 raster with one ring", "[MarchingSquares]") { - - sla::RasterBase::PixelDim pixdim{1, 1}; - +TEST_CASE("4x4 raster with one ring", "[MarchingSquares]") +{ + sla::PixelDim pixdim{1, 1}; + // We need one additional row and column to detect edges sla::RasterGrayscaleAA rst{{4, 4}, pixdim, {}, agg::gamma_threshold(.5)}; - - // Draw a triangle from individual pixels - rst.draw(square(1., {0500000, 0500000})); - rst.draw(square(1., {1500000, 0500000})); - rst.draw(square(1., {2500000, 0500000})); - - rst.draw(square(1., {1500000, 1500000})); - rst.draw(square(1., {2500000, 1500000})); - - rst.draw(square(1., {2500000, 2500000})); - - std::fstream out("4x4.png", std::ios::out); - out << rst.encode(sla::PNGRasterEncoder{}); - out.close(); - - ExPolygons extracted = sla::raster_to_polygons(rst); - - SVG svg("4x4.svg"); - svg.draw(extracted); - svg.Close(); - - REQUIRE(extracted.size() == 1); + + ExPolygons one = {{{1, 1}, {3, 1}, {3, 3}, {2, 3}, {2, 2}, {1, 2}}}; + for (ExPolygon& p : one) + p.scale(scaled(1.0)); + test_expolys(rst, one, W1x1, "one_4x4"); } -TEST_CASE("4x4 raster with two rings", "[MarchingSquares]") { - - sla::RasterBase::PixelDim pixdim{1, 1}; - +TEST_CASE("10x10 raster with two rings", "[MarchingSquares]") +{ + sla::PixelDim pixdim{1, 1}; + // We need one additional row and column to detect edges - sla::RasterGrayscaleAA rst{{5, 5}, pixdim, {}, agg::gamma_threshold(.5)}; - - SECTION("Ambiguous case with 'ac' square") { - - // Draw a triangle from individual pixels - rst.draw(square(1., {3500000, 2500000})); - rst.draw(square(1., {3500000, 3500000})); - rst.draw(square(1., {2500000, 3500000})); - - rst.draw(square(1., {2500000, 1500000})); - rst.draw(square(1., {1500000, 1500000})); - rst.draw(square(1., {1500000, 2500000})); - - std::fstream out("4x4_ac.png", std::ios::out); - out << rst.encode(sla::PNGRasterEncoder{}); - out.close(); - - ExPolygons extracted = sla::raster_to_polygons(rst); - - SVG svg("4x4_ac.svg"); - svg.draw(extracted); - svg.Close(); - - REQUIRE(extracted.size() == 2); + sla::RasterGrayscaleAA rst{{10, 10}, pixdim, {}, agg::gamma_threshold(.5)}; + + SECTION("Ambiguous case with 'bd' square") + { + ExPolygons ac = {{{1, 1}, {3, 1}, {3, 2}, {2, 2}, {2, 3}, {1, 3}}, {{4, 4}, {2, 4}, {2, 3}, {3, 3}, {3, 2}, {4, 2}}}; + for (ExPolygon& p : ac) + p.scale(scaled(2.0)); + test_expolys(rst, ac, W1x1, "bd_10x10", false); } - - SECTION("Ambiguous case with 'bd' square") { - - // Draw a triangle from individual pixels - rst.draw(square(1., {3500000, 1500000})); - rst.draw(square(1., {3500000, 2500000})); - rst.draw(square(1., {2500000, 1500000})); - - rst.draw(square(1., {1500000, 2500000})); - rst.draw(square(1., {1500000, 3500000})); - rst.draw(square(1., {2500000, 3500000})); - - std::fstream out("4x4_bd.png", std::ios::out); - out << rst.encode(sla::PNGRasterEncoder{}); - out.close(); - - ExPolygons extracted = sla::raster_to_polygons(rst); - - SVG svg("4x4_bd.svg"); - svg.draw(extracted); - svg.Close(); - - REQUIRE(extracted.size() == 2); + + SECTION("Ambiguous case with 'ac' square") + { + ExPolygons bd = {{{1, 4}, {1, 2}, {2, 2}, {2, 3}, {3, 3}, {3, 4}}, {{4, 1}, {4, 3}, {3, 3}, {3, 2}, {2, 2}, {2, 1}}}; + for (ExPolygon& p : bd) + p.scale(scaled(2.0)); + test_expolys(rst, bd, W1x1, "ac_10x10", false); } } -TEST_CASE("Square with hole in the middle", "[MarchingSquares]") { +TEST_CASE("Square with hole in the middle", "[MarchingSquares]") +{ using namespace Slic3r; - + ExPolygons inp = {square_with_hole(50.)}; - - SECTION("Proportional raster, 1x1 mm pixel size, full accuracy") { - test_expolys(create_raster({100, 100}, 100., 100.), inp, W2x2, "square_with_hole_proportional_1x1_mm_px_full"); + + SECTION("Proportional raster, 1x1 mm pixel size, full accuracy") + { + test_expolys(create_raster({100, 100}, 100., 100.), inp, W1x1, "square_with_hole_proportional_1x1_mm_px_full"); } - - SECTION("Proportional raster, 1x1 mm pixel size, half accuracy") { - test_expolys(create_raster({100, 100}, 100., 100.), inp, W4x4, "square_with_hole_proportional_1x1_mm_px_half"); + + SECTION("Proportional raster, 1x1 mm pixel size, half accuracy") + { + test_expolys(create_raster({100, 100}, 100., 100.), inp, W2x2, "square_with_hole_proportional_1x1_mm_px_half"); } - - SECTION("Landscape raster, 1x1 mm pixel size, full accuracy") { - test_expolys(create_raster({150, 100}, 150., 100.), inp, W2x2, "square_with_hole_landsc_1x1_mm_px_full"); + + SECTION("Landscape raster, 1x1 mm pixel size, full accuracy") + { + test_expolys(create_raster({150, 100}, 150., 100.), inp, W1x1, "square_with_hole_landsc_1x1_mm_px_full"); } - - SECTION("Landscape raster, 1x1 mm pixel size, half accuracy") { - test_expolys(create_raster({150, 100}, 150., 100.), inp, W4x4, "square_with_hole_landsc_1x1_mm_px_half"); + + SECTION("Landscape raster, 1x1 mm pixel size, half accuracy") + { + test_expolys(create_raster({150, 100}, 150., 100.), inp, W2x2, "square_with_hole_landsc_1x1_mm_px_half"); } - - SECTION("Portrait raster, 1x1 mm pixel size, full accuracy") { - test_expolys(create_raster({100, 150}, 100., 150.), inp, W2x2, "square_with_hole_portrait_1x1_mm_px_full"); + + SECTION("Portrait raster, 1x1 mm pixel size, full accuracy") + { + test_expolys(create_raster({100, 150}, 100., 150.), inp, W1x1, "square_with_hole_portrait_1x1_mm_px_full"); } - - SECTION("Portrait raster, 1x1 mm pixel size, half accuracy") { - test_expolys(create_raster({100, 150}, 100., 150.), inp, W4x4, "square_with_hole_portrait_1x1_mm_px_half"); + + SECTION("Portrait raster, 1x1 mm pixel size, half accuracy") + { + test_expolys(create_raster({100, 150}, 100., 150.), inp, W2x2, "square_with_hole_portrait_1x1_mm_px_half"); } - - SECTION("Proportional raster, 2x2 mm pixel size, full accuracy") { - test_expolys(create_raster({200, 200}, 100., 100.), inp, W2x2, "square_with_hole_proportional_2x2_mm_px_full"); + + SECTION("Proportional raster, 2x2 mm pixel size, full accuracy") + { + test_expolys(create_raster({50, 50}, 100., 100.), inp, W1x1, "square_with_hole_proportional_2x2_mm_px_full"); } - - SECTION("Proportional raster, 2x2 mm pixel size, half accuracy") { - test_expolys(create_raster({200, 200}, 100., 100.), inp, W4x4, "square_with_hole_proportional_2x2_mm_px_half"); + + SECTION("Proportional raster, 2x2 mm pixel size, half accuracy") + { + test_expolys(create_raster({50, 50}, 100., 100.), inp, W2x2, "square_with_hole_proportional_2x2_mm_px_half"); } - - SECTION("Proportional raster, 0.5x0.5 mm pixel size, full accuracy") { - test_expolys(create_raster({50, 50}, 100., 100.), inp, W2x2, "square_with_hole_proportional_0.5x0.5_mm_px_full"); + + SECTION("Proportional raster, 0.5x0.5 mm pixel size, full accuracy") + { + test_expolys(create_raster({200, 200}, 100., 100.), inp, W1x1, "square_with_hole_proportional_0.5x0.5_mm_px_full"); } - - SECTION("Proportional raster, 0.5x0.5 mm pixel size, half accuracy") { - test_expolys(create_raster({50, 50}, 100., 100.), inp, W4x4, "square_with_hole_proportional_0.5x0.5_mm_px_half"); + + SECTION("Proportional raster, 0.5x0.5 mm pixel size, half accuracy") + { + test_expolys(create_raster({200, 200}, 100., 100.), inp, W2x2, "square_with_hole_proportional_0.5x0.5_mm_px_half"); } } -TEST_CASE("Circle with hole in the middle", "[MarchingSquares]") { +TEST_CASE("Circle with hole in the middle", "[MarchingSquares]") +{ using namespace Slic3r; - - test_expolys(create_raster({1000, 1000}), circle_with_hole(25.), W2x2, "circle_with_hole"); + + test_expolys(create_raster({1000, 1000}), circle_with_hole(25.), W1x1, "circle_with_hole"); } -static void recreate_object_from_rasters(const std::string &objname, float lh) { +static void recreate_object_from_rasters(const std::string& objname, float lh) +{ TriangleMesh mesh = load_model(objname); - - auto bb = mesh.bounding_box(); + + auto bb = mesh.bounding_box(); Vec3f tr = -bb.center().cast(); mesh.translate(tr.x(), tr.y(), tr.z()); bb = mesh.bounding_box(); - + std::vector layers = slice_mesh_ex(mesh.its, grid(float(bb.min.z()) + lh, float(bb.max.z()), lh)); - - sla::RasterBase::Resolution res{2560, 1440}; - double disp_w = 120.96; - double disp_h = 68.04; + + sla::Resolution res{2560, 1440}; + double disp_w = 120.96; + double disp_h = 68.04; #ifndef NDEBUG size_t cntr = 0; #endif - for (ExPolygons &layer : layers) { + for (ExPolygons& layer : layers) { auto rst = create_raster(res, disp_w, disp_h); - - for (ExPolygon &island : layer) { + + for (ExPolygon& island : layer) { rst.draw(island); } - + #ifndef NDEBUG std::fstream out(objname + std::to_string(cntr) + ".png", std::ios::out); out << rst.encode(sla::PNGRasterEncoder{}); out.close(); #endif - + ExPolygons layer_ = sla::raster_to_polygons(rst); -// float delta = scaled(std::min(rst.pixel_dimensions().h_mm, -// rst.pixel_dimensions().w_mm)) / 2; - -// layer_ = expolygons_simplify(layer_, delta); + // float delta = scaled(std::min(rst.pixel_dimensions().h_mm, + // rst.pixel_dimensions().w_mm)) / 2; + + // layer_ = expolygons_simplify(layer_, delta); #ifndef NDEBUG - SVG svg(objname + std::to_string(cntr) + ".svg", BoundingBox(Point{0, 0}, Point{scaled(disp_w), scaled(disp_h)})); + SVG svg(objname + std::to_string(cntr) + ".svg", rstBBox(rst)); svg.draw(layer_); svg.draw(layer, "green"); svg.Close(); #endif - + double layera = 0., layera_ = 0.; - for (auto &p : layer) layera += p.area(); - for (auto &p : layer_) layera_ += p.area(); + for (auto& p : layer) + layera += p.area(); + for (auto& p : layer_) + layera_ += p.area(); #ifndef NDEBUG std::cout << cntr++ << std::endl; #endif double diff = std::abs(layera_ - layera); REQUIRE((diff <= 0.1 * layera || diff < scaled(1.) * scaled(1.))); - + layer = std::move(layer_); } - + indexed_triangle_set out = slices_to_mesh(layers, bb.min.z(), double(lh), double(lh)); its_write_obj(out, "out_from_rasters.obj"); } -TEST_CASE("Recreate object from rasters", "[SL1Import]") { - recreate_object_from_rasters("frog_legs.obj", 0.05f); +TEST_CASE("Recreate object from rasters", "[SL1Import]") { recreate_object_from_rasters("frog_legs.obj", 0.05f); } + +namespace marchsq { + +static constexpr float layerf = 0.20; // layer height in mm (used for z values). +static constexpr float gsizef = 100.0; // grid size in mm (box volume side length). +static constexpr float wsizef = 0.50; // grid window size in mm (roughly line segment length). +static constexpr float psizef = 0.01; // raster pixel size in mm (roughly point accuracy). +static constexpr float isoval = 0.0; // iso value threshold to use. +static constexpr size_t wsize = std::round(wsizef / psizef); + +static float period = 10.0; // gyroid "wavelength" in mm (2x line spacing). +static float freq = 2 * PI / period; // gyroid frequency in waves per mm. + +void set_period(float len = 10.0) +{ + period = len; + freq = 2 * PI / period; } + +static size_t layer_n; +static size_t ring_n; +static size_t point_n; +static size_t get_n; + +void reset_stats() +{ + layer_n = 0; + ring_n = 0; + point_n = 0; + get_n = 0; +} + +using Rings = std::vector; + +template<> struct _RasterTraits +{ + // using Rst = Slic3r::sla::RasterGrayscaleAA; + // The type of pixel cell in the raster + using ValueType = float; + + // Value at a given position + static float get(const size_t& layer, size_t row, size_t col) + { + get_n++; + const float x = col * psizef * freq; + const float y = row * psizef * freq; + const float z = layer * psizef * freq; + + return sinf(x) * cosf(y) + sinf(y) * cosf(z) + sinf(z) * cosf(x); + } + + // Number of rows and cols of the raster + static size_t rows(const size_t& layer) { return std::round(gsizef / psizef); } + static size_t cols(const size_t& layer) { return std::round(gsizef / psizef); } +}; + +Rings get_gyroids(size_t l) +{ + size_t layer = l; + Rings rings = execute(layer, isoval, {wsize, wsize}); + layer_n++; + ring_n += rings.size(); + for (auto r : rings) + point_n += r.size(); + return rings; +} + +}; // namespace marchsq + +void benchmark_gyroid(float period) +{ + marchsq::reset_stats(); + marchsq::set_period(period); + INFO("grid size: " << marchsq::gsizef << "mm\nlayer height: " << marchsq::layerf << "mm\n"); + INFO("window size: " << marchsq::wsizef << "mm\npoint size: " << marchsq::psizef << "mm\n"); + INFO("gyroid period: " << marchsq::period << "mm\n"); + BENCHMARK("indexed", i) { return marchsq::get_gyroids(i); }; + INFO("output avg rings/layer: " << float(marchsq::ring_n) / float(marchsq::layer_n) << "\n"); + INFO("output avg points/layer: " << float(marchsq::point_n) / float(marchsq::layer_n) << "\n"); + INFO("output avg gets/layer: " << float(marchsq::get_n) / float(marchsq::layer_n) << "\n"); + + REQUIRE(marchsq::layer_n > 0); +} + +TEST_CASE("Benchmark gyroid cube period 10.0mm", "[MarchingSquares]") { benchmark_gyroid(10.0); } + +TEST_CASE("Benchmark gyroid cube period 5.0mm", "[MarchingSquares]") { benchmark_gyroid(5.0); }