mirror of
https://github.com/OrcaSlicer/OrcaSlicer.git
synced 2026-04-06 00:32:05 +02:00
Optimize and simplify MarchingSquares.hpp. (#10747)
* Optimize and simplify MarchingSquares.hpp, and fix it's test. This changes the implementation to get the possible next directions for a cell when building the tags and clearing them as the cells are visited during the march, instead of adding the visited previous direction to the tags during the march. The Dir enum has been turned into bit flags that for the possible next directions with boolean operators for testing/setting/clearing them. This simplifies and optimizes many operations during the march and building the polygons. The complicated/broken and unused partial support for cell overlap has been removed, simplifying the overly confusing grid iteration logic. The broken test has been fixed by removing the now gone `RasterBase` namespace from `sla::RasterBase::Pixeldim` and `sla:RasterBase:Resolution`, and the CMakeLists.txt entry uncommented. make Dir into flags * Further optimize MarchingSquares.hpp and improve comments. * Switch from a single byte-vector containing tags and dirs for each cell to a m_tags vector of bit-packed tags for each grid corner and an m_dirs vector of packed 4bit dirs for each cell. Since each grid corner tag is shared by the 4 adjacent cells this significantly reduces storage space and avoids redundantly calculating each tag 4x. It also significantly improves memory locality with each phase of calculating tags, calculating dirs, calculating rings operating only on the tags or dirs data required without them being interleaved with the data they don't need. * Change NEXT_CCW to be initialized with a static constexpr lambda instead of a manually entered table. This avoids typo errors manually building the table. * Optimize search_start_cell() so it can efficiently skip over cleared blocks of 8 dirs in the packed m_dirs vector. * Change the tags logical labeling to better suit the packed tags vector data. This makes it a tiny bit more efficient to extract from the m_tags bitmap. * Remove the now unused SquareTag enum class. * Add comments explaining the algorithm, including corner-cases in cell iteration. * Remove unused Dir operators and get_dirs() argument, and clang-format. * Fix some bugs and add stream output operators for debugging. * Fix a bug building tags where `step(gcrd, Dir::right)` was not assigned to update the gcrd grid point. Perhaps this should be a mutating method, or even a += operator? Also when wrapping at the end of a row it was updating the gcrd grid point by mutating the p raster point instead of itself. Perhaps Grid and Raster points should be different types? Maybe even templated? * Fix a bug in get_tags() when the second row tags are packed into any of the 2 LSB's of the uint32_t blocks. In hind-sight obviously `>>(o - 2)` will not shift left when `o < 2`. * Move interpolation of the edge-crossings into a `interpolate()` method, and make it shift bottom and right side points "out" by one to account for raster pixel width. This makes the results track the raster shapes much more accurately for very small windows. * Make `interpolate_rings()` check for and remove duplicated points. It turns out it's pretty common that two edge-crossing-points at a corner interpolate to the same point. This can also happen for the first and last points. * For Coord add `==` and `!=` operators, and use them wherever Coord's are compared. * Add `<<` stream output operators for Coord, Ring, and Dir classes. Add `streamtags(<stream>)` and `streamdirs(<stream>)` methods for dumping the tags and dirs data in an easy to understand text format. These make print-debugging much easier. * Add `assert(idx < m_gridlen)` in a bunch of places where grid-indexes are used. * For test_clipper_utils.cpp fix three "ambiguous overloading" compiler errors. This just adds three `Polygons` qualifications to fix compiler errors about ambiguous overloaded methods. Note this file was formated with a mixture of tabs and spaces and had lots of trailing whitespace. My editor cleaned these up resulting in a large looking diff, but if you use `git diff -w` to ignore the whitespace changes you will see it is actually tiny. errros * Update SLA/RasterToPolygons.* for MarchingSquares.hpp improvements. Change the minimum and default window size from 2x2 to 1x1. Also remove the strange pixel size re-scaling by (resolution/resolution-1). The old MarchingSquares implementation had complications around a default minimum 1 pixel "overlap" between cells which messed with the scaling a tiny bit and meant when you requested a 2x2 window size it actually used a 1x1 window. Both of these meant you had to specify a window 1 pixel larger than you really wanted, and you needed to undo the strange scaling artifact for accurate dimensions of your results. This has been fixed/removed in the new implementation, so the window is the window, there is no overlap, and no strange miss-scaling. * Fix test_marchingsquares.cpp and add StreamUtils.hpp. This fixes the MarchingSquares unittests to both pass and be more strict than they were before. It also adds libslic3r/StreamUtils.hpp which includes some handy streaming operators for standard libslic3r classes used to show extracted polys in the unittests. * Change Format/SL1.cpp to support the min 1x1 window for MarchingSquares. * Fix the ring-walk termination condition. Terminate the ring-walk when we return to the starting cell instead of when we reach a cell with no remaining directions. This ensures we don't merge two polygons if we started on an ambiguous cell. * Revert the removal of duplicate points in interpolate_rings(). It turns out that duplicate points are only relatively common when using a 1x1 window. These happen when the line passes through the corner pixel on a top-left corner in the raster, and the probability of this rapidly declines as the window increases, so in many cases this filtering is just overhead. It can also be potentially useful to see the points for every edge crossing even if they are duplicates. This kind of filtering is already done and done better in the polygon post-processing. * rename `interpolate()` to `interpolate_edge()`, make it update the point in-place, and add asserts to ensure the input point is a valid edge interpolation point. * Remove the duplicate point filtering from `interpolate_rings()` and simplify it. * Optimize directions building. This optimizes `get_dirs_block8()` to rapidly skip over blocks where the tags produce no directions (all tags are 1's or 0's), and also to build the directions faster when it has to by fetching the whole blocks worth of tags at once instead of cell-by-cell. * Rename `get_tags()` to `get_tags9()` and make it fetch a row of nine tags instead of the tags for a single cell. * Optimize `get_dirs_block8()` to use `get_tags9()` to get the next nine tags for the current and next rows and then shift through them to generate the tags and directions for each cell in the block. Also abort early and just return an empty block if the tags are all 0's or all 1's. * Tiny optimization for `get_tags_block32()`. This avoids using the `step()` method for a simple step-right that can be done with a simple increment of the column. It also avoids re-calculating the raster-coodinates for every corner, instead incrementing the column by `m_window.c` until the end of a row. * Fix svg output in test_marchingsquares.cpp for recreate_object_from_rasters. These SVG's were not properly centered... * Fix 2 static_casts for compiling on Windows. Thanks to RF47 for pointing this out on the #10747 pull request. * Make edge iteration use O(ln(N)) binary search instead of linear. This should be much faster when the window size is large. * Make `CellIt` into a `std::random_access_iterator_tag` so that `std::lower_bound()` can use a binary search to find the point on the edge instead of a linear search. * Change `step()` to support an optional distance argument and make it modify the `Coord` in-place instead of return a new one. * Update tests for the `step()` change. * Add Catch2 BENCHMARK tests for MarchingSquares. This required enabling the benchmarks in the tests/CMakeLists.txt config. * Add a _Loop<> specialization for parallel execution using ExecutionTBB. This is something that could be added wherever you are going to use this, but I intend on using this in multiple places so we might as add this once in one place where it can be reused. * Fix whitespace in messed up by tab-replacements. My editor renders, and replaces, tabs as 8 spaces. This messed up the indenting in tests/libslic3r/CMakeLists.txt and tests/libslic3r/test_clipper_utils.cpp when I made tiny changes in them. This fixes the indenting using 4 chars. Note it will still show as a diff because it is replacing tabs with 4 spaces, and removing trailing whitespace. But at least it's now indented correctly... --------- Co-authored-by: Donovan Baarda <dbaarda@google.com> Co-authored-by: SoftFever <softfeverever@gmail.com>
This commit is contained in:
Donovan Baarda
GitHub
@@ -305,8 +305,8 @@ ConfigSubstitutions import_sla_archive(
|
||||
std::function<bool(int)> 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);
|
||||
|
||||
@@ -1,49 +1,160 @@
|
||||
#ifndef MARCHINGSQUARES_HPP
|
||||
#define MARCHINGSQUARES_HPP
|
||||
|
||||
#include "Execution/ExecutionTBB.hpp"
|
||||
#include <type_traits>
|
||||
#include <cstdint>
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
|
||||
// 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 <cbda>;
|
||||
//
|
||||
// ^ 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 <cbda> 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
|
||||
// <c_b_> 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
|
||||
// <cb__> 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 <urdl>
|
||||
// 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<Coord>;
|
||||
|
||||
// Specialize this struct to register a raster type for the Marching squares alg
|
||||
template<class T, class Enable = void> 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<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);
|
||||
|
||||
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<class ExecutionPolicy, class Enable = void> struct _Loop {
|
||||
template<class It, class Fn> static void for_each(It from, It to, Fn &&fn)
|
||||
template<class ExecutionPolicy, class Enable = void> struct _Loop
|
||||
{
|
||||
template<class It, class Fn> 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<ExecutionTBB>
|
||||
{
|
||||
template<class It, class Fn> 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<class T> using RasterTraits = _RasterTraits<std::decay_t<T>>;
|
||||
template<class T> using TRasterValue = typename RasterTraits<T>::ValueType;
|
||||
|
||||
template<class T> size_t rows(const T &raster)
|
||||
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 T> TRasterValue<T> isoval(const T& rst, const Coord& crd) { return RasterTraits<T>::get(rst, crd.r, crd.c); }
|
||||
|
||||
template<class ExecutionPolicy, class It, class Fn> void for_each_idx(ExecutionPolicy&& policy, It from, It to, Fn&& fn)
|
||||
{
|
||||
return RasterTraits<T>::rows(raster);
|
||||
_Loop<ExecutionPolicy>::for_each_idx(from, to, fn);
|
||||
}
|
||||
|
||||
template<class T> size_t cols(const T &raster)
|
||||
{
|
||||
return RasterTraits<T>::cols(raster);
|
||||
}
|
||||
template<class E> constexpr std::underlying_type_t<E> _t(E e) noexcept { return static_cast<std::underlying_type_t<E>>(e); }
|
||||
|
||||
template<class T> TRasterValue<T> isoval(const T &rst, const Coord &crd)
|
||||
{
|
||||
return RasterTraits<T>::get(rst, crd.r, crd.c);
|
||||
}
|
||||
|
||||
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
|
||||
// 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<class E> constexpr std::underlying_type_t<E> _t(E e) noexcept
|
||||
inline std::ostream& operator<<(std::ostream& os, const Dir& d) { return os << ".<v#>x##^#X#####"[_t(d)]; }
|
||||
|
||||
// This maps square tag column/row order <cbda> 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<uint32_t, 16> 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<std::underlying_type_t<E>>(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 Rst> 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<uint32_t> m_tags; // bit-packed squaretags for each corner.
|
||||
std::vector<uint32_t> m_dirs; // bit-packed next directions for each cell.
|
||||
|
||||
template<class Rst> class Grid {
|
||||
const Rst * m_rst = nullptr;
|
||||
Coord m_cellsize, m_res_1, m_window, m_gridsize, m_grid_1;
|
||||
std::vector<uint8_t> 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<Rst> 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<Rst> 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<Rst> 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 <urdl> 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<uint32_t>(d) << (o * 4));
|
||||
}
|
||||
|
||||
// Get directions in a cell's <uldr> 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<Rst>;
|
||||
using pointer = TRasterValue<Rst> *;
|
||||
using reference = TRasterValue<Rst> &;
|
||||
using pointer = TRasterValue<Rst>*;
|
||||
using reference = TRasterValue<Rst>&;
|
||||
using difference_type = long;
|
||||
using iterator_category = std::forward_iterator_tag;
|
||||
|
||||
Coord crd;
|
||||
Dir dir = Dir::none;
|
||||
const Rst* grid = nullptr;
|
||||
|
||||
inline TRasterValue<Rst> operator*() const { return isoval(*grid, crd); }
|
||||
inline TRasterValue<Rst>& 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<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);
|
||||
});
|
||||
|
||||
void interpolate_edge(Coord& ecrd, TRasterValue<Rst> 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<long>(rows(rst)), static_cast<long>(cols(rst))}
|
||||
, m_gridsize{2 + m_rastsize.r / m_window.r, 2 + m_rastsize.c / m_window.c}
|
||||
, m_gridlen{static_cast<size_t>(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<class ExecutionPolicy> void tag_grid(ExecutionPolicy&& policy, TRasterValue<Rst> isoval)
|
||||
{
|
||||
// Get all the tags. parallel?
|
||||
for_each_idx(std::forward<ExecutionPolicy>(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<ExecutionPolicy>(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<Ring> scan_rings()
|
||||
{
|
||||
std::vector<Ring> 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<class ExecutionPolicy>
|
||||
void interpolate_rings(ExecutionPolicy && policy,
|
||||
std::vector<Ring> &rings,
|
||||
TRasterValue<Rst> isov)
|
||||
// sequential 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;
|
||||
}
|
||||
for_each_idx(std::forward<ExecutionPolicy>(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<class Raster, class ExecutionPolicy>
|
||||
std::vector<marchsq::Ring> execute_with_policy(ExecutionPolicy && policy,
|
||||
const Raster & raster,
|
||||
std::vector<marchsq::Ring> execute_with_policy(ExecutionPolicy&& policy,
|
||||
const Raster& raster,
|
||||
TRasterValue<Raster> 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<Raster> grid{raster, windowsize, overlap};
|
||||
|
||||
|
||||
Grid<Raster> grid{raster, windowsize};
|
||||
|
||||
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 = {})
|
||||
template<class Raster> std::vector<marchsq::Ring> execute(const Raster& raster, TRasterValue<Raster> isoval, Coord windowsize = {})
|
||||
{
|
||||
return execute_with_policy(nullptr, raster, isoval, windowsize);
|
||||
}
|
||||
|
||||
@@ -10,13 +10,13 @@ namespace marchsq {
|
||||
// Specialize this struct to register a raster type for the Marching squares alg
|
||||
template<> struct _RasterTraits<Slic3r::sla::RasterGrayscaleAA> {
|
||||
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<class Fn> 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<marchsq::Ring> 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;
|
||||
}
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
60
src/libslic3r/StreamUtils.hpp
Normal file
60
src/libslic3r/StreamUtils.hpp
Normal file
@@ -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 <sstream>
|
||||
#include <vector>
|
||||
|
||||
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_
|
||||
@@ -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
|
||||
|
||||
@@ -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}: ")
|
||||
|
||||
@@ -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<e_ordering o = e_ordering::OFF, class P, class Tree, class Alloc>
|
||||
template<e_ordering o = e_ordering::OFF, class P, class Tree, class Alloc>
|
||||
double polytree_area(const Tree &tree, std::vector<P, Alloc> *out)
|
||||
{
|
||||
traverse_pt<o>(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<e_ordering::ON>(tree.GetFirst(), &output)));
|
||||
REQUIRE(output.size() == reference.size());
|
||||
}
|
||||
|
||||
|
||||
SECTION("Traverse into ExPolygons WITH spatial ordering") {
|
||||
ExPolygons output;
|
||||
REQUIRE(area_sum == Approx(polytree_area<e_ordering::ON>(tree.GetFirst(), &output)));
|
||||
|
||||
@@ -17,22 +17,111 @@
|
||||
#include <libslic3r/TriangulateWall.hpp>
|
||||
#include <libslic3r/Tesselate.hpp>
|
||||
#include <libslic3r/SlicesToTriangleMesh.hpp>
|
||||
#include <libslic3r/StreamUtils.hpp>
|
||||
|
||||
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<double> 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<class Rst>
|
||||
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<double>(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<float>();
|
||||
mesh.translate(tr.x(), tr.y(), tr.z());
|
||||
bb = mesh.bounding_box();
|
||||
|
||||
|
||||
std::vector<ExPolygons> 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<double>(1.) * scaled<double>(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<Ring>;
|
||||
|
||||
template<> struct _RasterTraits<size_t>
|
||||
{
|
||||
// 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); }
|
||||
|
||||
Reference in New Issue
Block a user