3DPrinting/PrusaSlicer-NonPlainar
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Refactoring tests

Browse source
This commit is contained in:
Filip Sykala 2021-10-19 12:39:16 +02:00
parent c11245d82e
commit e8f7f3245f
9 changed files with 1612 additions and 717 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

15
src/libslic3r/Triangulation.cpp
View file

@ -26,7 +26,9 @@ inline void insert_edge(Triangulation::HalfEdges &edges, uint32_t &offset, const
} // namespace Private
Triangulation::Indices Triangulation::triangulate(const Points &points, const HalfEdges &half_edges)
Triangulation::Indices Triangulation::triangulate(const Points & points,
const HalfEdges &half_edges,
bool allow_opposit_edge)
{
// IMPROVE use int point insted of float !!!
@ -67,12 +69,11 @@ Triangulation::Indices Triangulation::triangulate(const Points &points, const Ha
for (size_t i = 0; i < 3; ++i) pi[i] = map[face->vertex(i)];
// Do not use triangles with opposit edges
if (half_edges.find(std::make_pair(pi[1], pi[0])) != half_edges.end())
continue;
if (half_edges.find(std::make_pair(pi[2], pi[1])) != half_edges.end())
continue;
if (half_edges.find(std::make_pair(pi[0], pi[2])) != half_edges.end())
continue;
if (!allow_opposit_edge) {
if (half_edges.find(std::make_pair(pi[1], pi[0])) != half_edges.end()) continue;
if (half_edges.find(std::make_pair(pi[2], pi[1])) != half_edges.end()) continue;
if (half_edges.find(std::make_pair(pi[0], pi[2])) != half_edges.end()) continue;
}
indices.emplace_back(pi[0], pi[1], pi[2]);
}

7
src/libslic3r/Triangulation.hpp
View file

@ -26,9 +26,12 @@ public:
/// </summary>
/// <param name="points">Points to connect</param>
/// <param name="edges">Constraint for edges, pair is from point(first) to
/// point(second)</param> <returns>Triangles</returns>
/// point(second)</param>
/// <param name="allow_opposit_edge">Flag for filtration result indices by opposit half edge</param>
/// <returns>Triangles</returns>
static Indices triangulate(const Points & points,
const HalfEdges &half_edges);
const HalfEdges &half_edges,
bool allow_opposit_edge = false);
static Indices triangulate(const Polygon &polygon);
static Indices triangulate(const Polygons &polygons);
static Indices triangulate(const ExPolygons &expolygons);

3
tests/libslic3r/CMakeLists.txt
View file

@ -22,6 +22,9 @@ add_executable(${_TEST_NAME}_tests
test_optimizers.cpp
test_png_io.cpp
test_timeutils.cpp
test_quadric_edge_collapse.cpp
test_triangulation.cpp
test_emboss.cpp
test_indexed_triangle_set.cpp
../libnest2d/printer_parts.cpp
)

1036
tests/libslic3r/Simplify.h Normal file
View file

File diff suppressed because it is too large Load diff

181
tests/libslic3r/test_emboss.cpp Normal file
View file

@ -0,0 +1,181 @@
#include <catch2/catch.hpp>
#include <libslic3r/Emboss.hpp>
#include <libslic3r/SVG.hpp> // only debug visualization
#include <optional>
#include <libslic3r/AABBTreeIndirect.hpp>
using namespace Slic3r;
namespace Private{
// calculate multiplication of ray dir to intersect - inspired by
// segment_segment_intersection when ray dir is normalized retur distance from
// ray point to intersection No value mean no intersection
std::optional<double> ray_segment_intersection(const Vec2d &r_point,
const Vec2d &r_dir,
const Vec2d &s0,
const Vec2d &s1)
{
auto denominate = [](const Vec2d &v0, const Vec2d &v1) -> double {
return v0.x() * v1.y() - v1.x() * v0.y();
};
Vec2d segment_dir = s1 - s0;
double d = denominate(segment_dir, r_dir);
if (std::abs(d) < std::numeric_limits<double>::epsilon())
// Line and ray are collinear.
return {};
Vec2d s12 = s0 - r_point;
double s_number = denominate(r_dir, s12);
bool change_sign = false;
if (d < 0.) {
change_sign = true;
d = -d;
s_number = -s_number;
}
if (s_number < 0. || s_number > d)
// Intersection outside of segment.
return {};
double r_number = denominate(segment_dir, s12);
if (change_sign) r_number = -r_number;
if (r_number < 0.)
// Intersection before ray start.
return {};
return r_number / d;
}
Vec2d get_intersection(const Vec2d & point,
const Vec2d & dir,
const std::array<Vec2d, 3> &triangle)
{
std::optional<double> t;
for (size_t i = 0; i < 3; ++i) {
size_t i2 = i + 1;
if (i2 == 3) i2 = 0;
if (!t.has_value()) {
t = ray_segment_intersection(point, dir, triangle[i],
triangle[i2]);
continue;
}
// small distance could be preccission inconsistance
std::optional<double> t2 = ray_segment_intersection(point, dir,
triangle[i],
triangle[i2]);
if (t2.has_value() && *t2 > *t) t = t2;
}
assert(t.has_value()); // Not found intersection.
return point + dir * (*t);
}
Vec3d calc_hit_point(const igl::Hit & h,
const Vec3i & triangle,
const std::vector<Vec3f> &vertices)
{
double c1 = h.u;
double c2 = h.v;
double c0 = 1.0 - c1 - c2;
Vec3d v0 = vertices[triangle[0]].cast<double>();
Vec3d v1 = vertices[triangle[1]].cast<double>();
Vec3d v2 = vertices[triangle[2]].cast<double>();
return v0 * c0 + v1 * c1 + v2 * c2;
}
Vec3d calc_hit_point(const igl::Hit &h, indexed_triangle_set &its)
{
return calc_hit_point(h, its.indices[h.id], its.vertices);
}
} // namespace Private
TEST_CASE("Emboss text", "[Emboss]")
{
const char *font_name = "C:/windows/fonts/arialbd.ttf";
char letter = '%';
float flatness = 2.;
std::optional<Emboss::Font> font = Emboss::load_font(font_name);
REQUIRE(font.has_value());
std::optional<Emboss::Glyph> glyph = Emboss::letter2glyph(*font, letter, flatness);
REQUIRE(glyph.has_value());
ExPolygons shape = glyph->shape;
REQUIRE(!shape.empty());
float z_depth = 1.f;
Emboss::ProjectZ projection(z_depth);
indexed_triangle_set its = Emboss::polygons2model(shape, projection);
CHECK(!its.indices.empty());
}
TEST_CASE("Test hit point", "[AABBTreeIndirect]")
{
indexed_triangle_set its;
its.vertices = {
Vec3f(1, 1, 1),
Vec3f(2, 10, 2),
Vec3f(10, 0, 2),
};
its.indices = {Vec3i(0, 2, 1)};
auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
its.vertices, its.indices);
Vec3d ray_point(8, 1, 0);
Vec3d ray_dir(0, 0, 1);
igl::Hit hit;
AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices, tree,
ray_point, ray_dir, hit);
Vec3d hp = Private::calc_hit_point(hit, its);
CHECK(abs(hp.x() - ray_point.x()) < .1);
CHECK(abs(hp.y() - ray_point.y()) < .1);
}
TEST_CASE("ray segment intersection", "[MeshBoolean]")
{
Vec2d r_point(1, 1);
Vec2d r_dir(1, 0);
// colinear
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 0), Vec2d(2, 0)).has_value());
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(0, 0)).has_value());
// before ray
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 0), Vec2d(0, 2)).has_value());
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(0, 2), Vec2d(0, 0)).has_value());
// above ray
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 2), Vec2d(2, 3)).has_value());
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 3), Vec2d(2, 2)).has_value());
// belove ray
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(2, -1)).has_value());
CHECK(!Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, -1), Vec2d(2, 0)).has_value());
// intersection at [2,1] distance 1
auto t1 = Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 0), Vec2d(2, 2));
REQUIRE(t1.has_value());
auto t2 = Private::ray_segment_intersection(r_point, r_dir, Vec2d(2, 2), Vec2d(2, 0));
REQUIRE(t2.has_value());
CHECK(abs(*t1 - *t2) < std::numeric_limits<double>::epsilon());
}
TEST_CASE("triangle intersection", "[]")
{
Vec2d point(1, 1);
Vec2d dir(-1, 0);
std::array<Vec2d, 3> triangle = {Vec2d(0, 0), Vec2d(5, 0), Vec2d(0, 5)};
Vec2d i = Private::get_intersection(point, dir, triangle);
CHECK(abs(i.x()) < std::numeric_limits<double>::epsilon());
CHECK(abs(i.y() - 1.) < std::numeric_limits<double>::epsilon());
}

217
tests/libslic3r/test_indexed_triangle_set.cpp
View file

@ -101,220 +101,3 @@ TEST_CASE("Split two watertight meshes", "[its_split][its]") {
debug_write_obj(res, "parts_watertight");
}
#include <libslic3r/QuadricEdgeCollapse.hpp>
static float triangle_area(const Vec3f &v0, const Vec3f &v1, const Vec3f &v2)
{
Vec3f ab = v1 - v0;
Vec3f ac = v2 - v0;
return ab.cross(ac).norm() / 2.f;
}
static float triangle_area(const Vec3crd &triangle_inices, const std::vector<Vec3f> &vertices)
{
return triangle_area(vertices[triangle_inices[0]],
vertices[triangle_inices[1]],
vertices[triangle_inices[2]]);
}
static std::mt19937 create_random_generator() {
std::random_device rd;
std::mt19937 gen(rd());
return gen;
}
std::vector<Vec3f> its_sample_surface(const indexed_triangle_set &its,
double sample_per_mm2,
std::mt19937 random_generator = create_random_generator())
{
std::vector<Vec3f> samples;
std::uniform_real_distribution<float> rand01(0.f, 1.f);
for (const auto &triangle_indices : its.indices) {
float area = triangle_area(triangle_indices, its.vertices);
float countf;
float fractional = std::modf(area * sample_per_mm2, &countf);
int count = static_cast<int>(countf);
float generate = rand01(random_generator);
if (generate < fractional) ++count;
if (count == 0) continue;
const Vec3f &v0 = its.vertices[triangle_indices[0]];
const Vec3f &v1 = its.vertices[triangle_indices[1]];
const Vec3f &v2 = its.vertices[triangle_indices[2]];
for (int c = 0; c < count; c++) {
// barycentric coordinate
Vec3f b;
b[0] = rand01(random_generator);
b[1] = rand01(random_generator);
if ((b[0] + b[1]) > 1.f) {
b[0] = 1.f - b[0];
b[1] = 1.f - b[1];
}
b[2] = 1.f - b[0] - b[1];
Vec3f pos;
for (int i = 0; i < 3; i++) {
pos[i] = b[0] * v0[i] + b[1] * v1[i] + b[2] * v2[i];
}
samples.push_back(pos);
}
}
return samples;
}
#include "libslic3r/AABBTreeIndirect.hpp"
struct CompareConfig
{
float max_distance = 3.f;
float max_average_distance = 2.f;
};
bool is_similar(const indexed_triangle_set &from,
const indexed_triangle_set &to,
const CompareConfig &cfg)
{
// create ABBTree
auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
from.vertices, from.indices);
float sum_distance = 0.f;
float max_distance = 0.f;
auto collect_distances = [&](const Vec3f &surface_point) {
size_t hit_idx;
Vec3f hit_point;
float distance2 =
AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
from.vertices, from.indices, tree, surface_point, hit_idx, hit_point);
float distance = sqrt(distance2);
if (max_distance < distance) max_distance = distance;
sum_distance += distance;
};
for (const Vec3f &vertex : to.vertices) {
collect_distances(vertex);
}
for (const Vec3i &t : to.indices) {
Vec3f center(0,0,0);
for (size_t i = 0; i < 3; ++i) {
center += to.vertices[t[i]] / 3;
}
collect_distances(center);
}
size_t count = to.vertices.size() + to.indices.size();
float avg_distance = sum_distance / count;
if (avg_distance > cfg.max_average_distance ||
max_distance > cfg.max_distance)
return false;
return true;
}
TEST_CASE("Reduce one edge by Quadric Edge Collapse", "[its]")
{
indexed_triangle_set its;
its.vertices = {Vec3f(-1.f, 0.f, 0.f), Vec3f(0.f, 1.f, 0.f),
Vec3f(1.f, 0.f, 0.f), Vec3f(0.f, 0.f, 1.f),
// vertex to be removed
Vec3f(0.9f, .1f, -.1f)};
its.indices = {Vec3i(1, 0, 3), Vec3i(2, 1, 3), Vec3i(0, 2, 3),
Vec3i(0, 1, 4), Vec3i(1, 2, 4), Vec3i(2, 0, 4)};
// edge to remove is between vertices 2 and 4 on trinagles 4 and 5
indexed_triangle_set its_ = its; // copy
// its_write_obj(its, "tetrhedron_in.obj");
uint32_t wanted_count = its.indices.size() - 1;
its_quadric_edge_collapse(its, wanted_count);
// its_write_obj(its, "tetrhedron_out.obj");
CHECK(its.indices.size() == 4);
CHECK(its.vertices.size() == 4);
for (size_t i = 0; i < 3; i++) {
CHECK(its.indices[i] == its_.indices[i]);
}
for (size_t i = 0; i < 4; i++) {
if (i == 2) continue;
CHECK(its.vertices[i] == its_.vertices[i]);
}
const Vec3f &v = its.vertices[2]; // new vertex
const Vec3f &v2 = its_.vertices[2]; // moved vertex
const Vec3f &v4 = its_.vertices[4]; // removed vertex
for (size_t i = 0; i < 3; i++) {
bool is_between = (v[i] < v4[i] && v[i] > v2[i]) ||
(v[i] > v4[i] && v[i] < v2[i]);
CHECK(is_between);
}
CompareConfig cfg;
cfg.max_average_distance = 0.014f;
cfg.max_distance = 0.75f;
CHECK(is_similar(its, its_, cfg));
CHECK(is_similar(its_, its, cfg));
}
#include "test_utils.hpp"
TEST_CASE("Simplify mesh by Quadric edge collapse to 5%", "[its]")
{
TriangleMesh mesh = load_model("frog_legs.obj");
double original_volume = its_volume(mesh.its);
uint32_t wanted_count = mesh.its.indices.size() * 0.05;
REQUIRE_FALSE(mesh.empty());
indexed_triangle_set its = mesh.its; // copy
float max_error = std::numeric_limits<float>::max();
its_quadric_edge_collapse(its, wanted_count, &max_error);
//its_write_obj(its, "frog_legs_qec.obj");
CHECK(its.indices.size() <= wanted_count);
double volume = its_volume(its);
CHECK(fabs(original_volume - volume) < 33.);
CompareConfig cfg;
cfg.max_average_distance = 0.043f;
cfg.max_distance = 0.32f;
CHECK(is_similar(mesh.its, its, cfg));
CHECK(is_similar(its, mesh.its, cfg));
}
bool exist_triangle_with_twice_vertices(const std::vector<stl_triangle_vertex_indices>& indices)
{
for (const auto &face : indices)
if (face[0] == face[1] ||
face[0] == face[2] ||
face[1] == face[2]) return true;
return false;
}
TEST_CASE("Simplify trouble case", "[its]")
{
TriangleMesh tm = load_model("simplification.obj");
REQUIRE_FALSE(tm.empty());
float max_error = std::numeric_limits<float>::max();
uint32_t wanted_count = 0;
its_quadric_edge_collapse(tm.its, wanted_count, &max_error);
CHECK(!exist_triangle_with_twice_vertices(tm.its.indices));
}
TEST_CASE("Simplified cube should not be empty.", "[its]")
{
auto its = its_make_cube(1, 2, 3);
float max_error = std::numeric_limits<float>::max();
uint32_t wanted_count = 0;
its_quadric_edge_collapse(its, wanted_count, &max_error);
CHECK(!its.indices.empty());
}
TEST_CASE("Neighbors in cube", "[its]")
{
auto its = its_make_cube(1, 2, 3);
auto neighbors = its_face_neighbors(its);
int no_value = -1;
for (auto &neighbor : neighbors) {
for (size_t i = 0; i < 3; i++) {
CHECK(neighbor[i] != no_value);
}
}
}

492
tests/libslic3r/test_meshboolean.cpp
View file

@ -39,7 +39,7 @@ Vec3d calc_normal(const Vec3i &triangle, const std::vector<Vec3f> &vertices)
TEST_CASE("Add TriangleMeshes", "[MeshBoolean]")
{
TriangleMesh tm1 = make_sphere(1.6, 1.6);
Vec3f move(5, -3, 7);
Vec3f move(5, -3, 7);
move.normalize();
tm1.translate(0.3 * move);
its_write_obj(tm1.its, "tm1.obj");
@ -48,493 +48,3 @@ TEST_CASE("Add TriangleMeshes", "[MeshBoolean]")
MeshBoolean::cgal::plus(tm1, tm2);
its_write_obj(tm1.its, "test_add.obj");
}
#include "libslic3r/Emboss.hpp"
ExPolygons ttf2polygons(const char * font_name, char letter, float flatness = 1.f) {
auto font = Emboss::load_font(font_name);
if (!font.has_value()) return ExPolygons();
return Emboss::letter2glyph(*font, letter, flatness)->shape;
}
#include "libslic3r/SVG.hpp"
void store_to_svg(Polygons polygons,std::string file_name = "letter.svg")
{
double scale = 1e6;
BoundingBox bb;
for (auto& p : polygons) {
p.scale(scale);
bb.merge(p.points);
}
SVG svg(file_name, bb);
svg.draw(polygons);
}
struct Plane
{
Vec3d point = Vec3d(0., 0., 0.); // lay on plane - define zero position
Vec3d normal = Vec3d(0., 0., 1.);// [unit vector] - define orientation
};
struct OrientedPlane: public Plane
{
// Must be perpendiculat to normal
Vec3d up = Vec3d(0., 1., 0.); // [unit vector]
};
struct EmbossConfig
{
// emboss plane must be above surface
// point define zero for 2d polygon and must be out of model
// normal is direction to model
// up define orientation of polygon
OrientedPlane projection;
// Move surface distance
// Positive value out of model (Raised)
// Negative value into model (Engraved)
float height = 1.; // [in milimeters]
};
#include <optional>
#include <libslic3r/AABBTreeIndirect.hpp>
Vec3d calc_hit_point(const igl::Hit & h,
const Vec3i & triangle,
const std::vector<Vec3f> &vertices)
{
double c1 = h.u;
double c2 = h.v;
double c0 = 1.0 - c1 - c2;
Vec3d v0 = vertices[triangle[0]].cast<double>();
Vec3d v1 = vertices[triangle[1]].cast<double>();
Vec3d v2 = vertices[triangle[2]].cast<double>();
return v0 * c0 + v1 * c1 + v2 * c2;
}
Vec3d calc_hit_point(const igl::Hit &h, indexed_triangle_set &its)
{
return calc_hit_point(h, its.indices[h.id], its.vertices);
}
TEST_CASE("Test hit point", "[AABBTreeIndirect]")
{
indexed_triangle_set its;
its.vertices = {
Vec3f(1,1,1),
Vec3f(2, 10, 2),
Vec3f(10, 0, 2),
};
its.indices = {Vec3i(0, 2, 1)};
auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
its.vertices, its.indices);
Vec3d ray_point(8, 1, 0);
Vec3d ray_dir(0,0,1);
igl::Hit hit;
AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices, tree,
ray_point, ray_dir, hit);
Vec3d hp = calc_hit_point(hit, its);
CHECK(abs(hp.x() - ray_point.x()) < .1);
CHECK(abs(hp.y() - ray_point.y()) < .1);
}
// represents triangle extend by seam
struct TrianglePath
{
// input edge, output edge, when cross triangle border
std::optional<std::pair<char, char>> edges;
// when edges has value than first and last point lay on triangle edge
std::vector<Vec3f> points;
// first point has index offset_id in result vertices
uint32_t offset_id;
};
using TrianglePaths = std::vector<TrianglePath>;
// create transformation matrix to convert direction vectors
// do not care about up vector
// Directions are normalized
Eigen::Matrix3d create_transformation(const Vec3d &from_dir, const Vec3d &to_dir)
{
Vec3d axis = from_dir.cross(to_dir);
axis.normalize();
double angle = acos(from_dir.dot(to_dir));
auto rotation = Eigen::AngleAxisd(angle, axis);
return rotation.matrix();
}
TEST_CASE("Transformation matrix", "[]") {
Vec3d d1(3, -7, 13);
Vec3d d2(-9, 5, 1);
d1.normalize();
d2.normalize();
auto tr_mat = create_transformation(d1, d2);
Vec3d d1_tr = tr_mat * d1;
Vec3d diff = d1_tr - d2;
double eps = 1e-12;
for (double d : diff)
CHECK(abs(d) < std::numeric_limits<double>::epsilon());
}
// calculate multiplication of ray dir to intersect - inspired by segment_segment_intersection
// when ray dir is normalized retur distance from ray point to intersection
// No value mean no intersection
std::optional<double> ray_segment_intersection(
const Vec2d& r_point, const Vec2d& r_dir, const Vec2d& s0, const Vec2d& s1){
auto denominate = [](const Vec2d& v0,const Vec2d& v1)->double{
return v0.x() * v1.y() - v1.x() * v0.y();
};
Vec2d segment_dir = s1 - s0;
double d = denominate(segment_dir, r_dir);
if (std::abs(d) < std::numeric_limits<double>::epsilon())
// Line and ray are collinear.
return {};
Vec2d s12 = s0 - r_point;
double s_number = denominate(r_dir, s12);
bool change_sign = false;
if (d < 0.) {
change_sign = true;
d = -d;
s_number = -s_number;
}
if (s_number < 0.|| s_number > d)
// Intersection outside of segment.
return {};
double r_number = denominate(segment_dir, s12);
if (change_sign)
r_number = - r_number;
if(r_number < 0.)
// Intersection before ray start.
return {};
return r_number / d;
}
TEST_CASE("ray segment intersection", "[MeshBoolean]")
{
Vec2d r_point(1,1);
Vec2d r_dir(1,0);
// colinear
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,0), Vec2d(2,0)).has_value());
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(0,0)).has_value());
// before ray
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,0), Vec2d(0,2)).has_value());
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(0,2), Vec2d(0,0)).has_value());
// above ray
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,2), Vec2d(2,3)).has_value());
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,3), Vec2d(2,2)).has_value());
// belove ray
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(2,-1)).has_value());
CHECK(!ray_segment_intersection(r_point, r_dir, Vec2d(2,-1), Vec2d(2,0)).has_value());
// intersection at [2,1] distance 1
auto t1 = ray_segment_intersection(r_point, r_dir, Vec2d(2,0), Vec2d(2,2));
REQUIRE(t1.has_value());
auto t2 = ray_segment_intersection(r_point, r_dir, Vec2d(2,2), Vec2d(2,0));
REQUIRE(t2.has_value());
CHECK(abs(*t1 - *t2) < std::numeric_limits<double>::epsilon());
}
Vec2d get_intersection(const Vec2d& point, const Vec2d& dir, const std::array<Vec2d, 3>& triangle)
{
std::optional<double> t;
for (size_t i = 0; i < 3; ++i) {
size_t i2 = i+1;
if(i2 == 3) i2 = 0;
if (!t.has_value()) {
t = ray_segment_intersection(point, dir, triangle[i], triangle[i2]);
continue;
}
// small distance could be preccission inconsistance
std::optional<double> t2 = ray_segment_intersection(
point, dir, triangle[i], triangle[i2]);
if (t2.has_value() && *t2 > *t) t = t2;
}
assert(t.has_value()); //Not found intersection.
return point + dir * (*t);
}
TEST_CASE("triangle intersection", "[]")
{
Vec2d point(1,1);
Vec2d dir(-1, 0);
std::array<Vec2d, 3> triangle = {
Vec2d(0, 0),
Vec2d(5, 0),
Vec2d(0, 5)
};
Vec2d i = get_intersection(point, dir, triangle);
CHECK(abs(i.x()) < std::numeric_limits<double>::epsilon());
CHECK(abs(i.y() - 1.) < std::numeric_limits<double>::epsilon());
}
#include <libslic3r/Geometry.hpp>
#include <libslic3r/Triangulation.hpp>
indexed_triangle_set emboss3d(const Polygon &shape, float height) {
// CW order of triangle indices
std::vector<Vec3i> shape_triangles = Triangulation::triangulate(shape);
indexed_triangle_set result;
const Points &pts = shape.points;
size_t count_point = pts.size();
result.vertices.reserve(2 * count_point);
// top points
std::transform(pts.begin(), pts.end(), std::back_inserter(result.vertices),
[](const Point &p) { return Vec3f(p.x(), p.y(), 0); });
// bottom points
std::transform(pts.begin(), pts.end(), std::back_inserter(result.vertices),
[height](const Point &p) { return Vec3f(p.x(), p.y(), height); });
result.indices.reserve(shape_triangles.size() * 2 + count_point*2);
// top triangles - change to CCW
for (const Vec3i &t : shape_triangles)
result.indices.emplace_back(t.x(), t.z(), t.y());
// bottom triangles - use CW
for (const Vec3i &t : shape_triangles)
result.indices.emplace_back(
t.x() + count_point,
t.y() + count_point,
t.z() + count_point
);
// quads around - zig zag by triangles
for (uint32_t i = 0; i < count_point; ++i) {
// previous index
uint32_t ip = (i == 0) ? (count_point - 1) : (i - 1);
// bottom indices
uint32_t i2 = i + count_point;
uint32_t ip2 = ip + count_point;
result.indices.emplace_back(i, i2, ip);
result.indices.emplace_back(ip2, ip, i2);
}
return result;
}
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_exact_constructions_kernel EK;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = PMP::parameters;
// is it realy neccessary to copy all?
Mesh its_to_mesh(const indexed_triangle_set &its)
{
Mesh mesh;
size_t vertices_count = its.vertices.size();
size_t edges_count = (its.indices.size() * 3) / 2;
size_t faces_count = its.indices.size();
mesh.reserve(vertices_count, edges_count, faces_count);
for (auto &v : its.vertices)
mesh.add_vertex(typename Mesh::Point{v.x(), v.y(), v.z()});
using VI = typename Mesh::Vertex_index;
for (auto &f : its.indices)
mesh.add_face(VI(f(0)), VI(f(1)), VI(f(2)));
return mesh;
}
indexed_triangle_set mesh_to_its(const Mesh &mesh)
{
indexed_triangle_set res;
res.vertices.reserve(mesh.num_vertices());
res.indices.reserve(mesh.num_faces());
for (auto &vi : mesh.vertices()) {
auto &v = mesh.point(vi); // index addresing, to compare same float point
res.vertices.emplace_back(v.x(),v.y(), v.z());
}
for (auto &face : mesh.faces()) {
auto vtc = mesh.vertices_around_face(mesh.halfedge(face));
int i = 0;
Vec3i facet;
for (auto v : vtc) {
if (i > 2) {
i = 0;
break;
}
facet(i++) = v;
}
if (i == 3) res.indices.emplace_back(facet);
}
return res;
}
void emboss3d_(const Polygon& shape, const EmbossConfig &cfg, indexed_triangle_set &its)
{
indexed_triangle_set shape3d = emboss3d(shape, cfg.height);
//its_write_obj(shape3d,"shape3d.obj");
Mesh m1 = its_to_mesh(its);
Mesh m2 = its_to_mesh(shape3d);
size_t in_vertices_count = its.vertices.size();
size_t in_indices_count = its.indices.size();
auto tt = m1.property_map<Mesh::Edge_index, bool>("e:removed");
Mesh::Property_map<Mesh::Edge_index, bool> ecm1 = tt.first; // hope in copy
PMP::corefine(m1, m2
, PMP::parameters::edge_is_constrained_map(ecm1)
//, PMP::parameters::do_not_modify(true)
);
//PMP::Corefinement::Intersection_of_triangle_meshes()
size_t count_true = 0;
for (const Mesh::Edge_index &e : edges(m1))
if (ecm1[e]) {
++count_true;
const Mesh::Halfedge_index &h = e.halfedge();
const Mesh::Halfedge_index &h2 = m1.opposite(h);
const Mesh::Vertex_index &v = m1.target(h);
std::cout <<
"edge " << e <<
" halfedge " << h <<
" target " << v <<
" index0 " << v.idx() <<
" index1 " << m1.target(h2).idx() <<
std::endl;
}
its = mesh_to_its(m1);
// twice and move additional vertices
Vec3f move = (cfg.projection.normal * cfg.height).cast<float>();
size_t new_add_vertice = its.vertices.size() - in_vertices_count;
its.vertices.reserve(its.vertices.size() + new_add_vertice);
for (size_t i = in_vertices_count; i < its.vertices.size(); ++i)
its.vertices.emplace_back(its.vertices[i] + move);
// zig zag edge collected edge
for (size_t i = in_indices_count; i < its.indices.size(); ++i) {
// new added triangle
Vec3i &t = its.indices[i];
std::array<bool, 3> is_new;
bool is_inner = true;
for (size_t i = 0; i < 3; ++i) {
bool is_n = t[0] >= in_vertices_count;
is_inner |= is_n;
is_new[i] = is_n;
}
}
its_write_obj(mesh_to_its(m1), "corefine1.obj");
its_write_obj(mesh_to_its(m2), "corefine2.obj");
/*MeshBoolean::cgal::CGALMesh cm1 = m1->m;
My_visitor<MeshBoolean::cgal::CGALMesh> sm_v;
CGAL::Polygon_mesh_processing::corefine(m1->m, m2->m,
CGAL::Polygon_mesh_processing::parameters::visitor(sm_v)
, CGAL::parameters::do_not_modify(true)
);*/
//TriangleMesh tm1 = cgal_to_triangle_mesh(*m1);
//store_obj("tm1.obj", &tm1);
//TriangleMesh tm2 = cgal_to_triangle_mesh(*m2);
//store_obj("tm2.obj", &tm1);
its = mesh_to_its(m1);
return;
}
void emboss3d(const Polygon& shape, const EmbossConfig &cfg, indexed_triangle_set &its)
{
indexed_triangle_set shape3d = emboss3d(shape, abs(cfg.height));
//its_write_obj(shape3d,"shape3d.obj");
auto m1 = MeshBoolean::cgal::triangle_mesh_to_cgal(its);
auto m2 = MeshBoolean::cgal::triangle_mesh_to_cgal(shape3d);
bool is_plus = shape.is_counter_clockwise() == (cfg.height > 0.f);
if (is_plus) {
MeshBoolean::cgal::plus(*m1, *m2);
} else {
MeshBoolean::cgal::minus(*m1, *m2);
}
TriangleMesh tm1 = cgal_to_triangle_mesh(*m1);
store_obj("tm1.obj", &tm1);
//TriangleMesh tm2 = cgal_to_triangle_mesh(*m2);
//store_obj("tm2.obj", &tm1);
its = tm1.its; // copy
}
TEST_CASE("Emboss polygon example", "[MeshBoolean]")
{
const char *font_name = "C:/windows/fonts/arialbd.ttf";
char letter = '%';
float flatness = 2.;
ExPolygons espolygons = ttf2polygons(font_name, letter, flatness);
//TriangleMesh tm = make_sphere(1., 1.); tm.scale(10.f);
TriangleMesh tm = make_cube(10., 5., 2.);
tm.translate(Vec3f(0, 0, 1.7));
Polygons polygons;
polygons = {Polygon({{1, 1}, {1, 2}, {2, 2}, {2, 1}})}; // rectangle CW
polygons = {Polygon({{1, 1}, {2, 1}, {2, 2}, {1, 2}})}; // rectangle CCW
EmbossConfig ec;
ec.height = 3;
indexed_triangle_set its = tm.its; // copy
for (const auto& polygon : polygons) {
// TODO: differ CW and CCW
emboss3d_(polygon, ec, its);
}
}
TEST_CASE("Triangulate by cgal", "[Triangulation]")
{
Points points = {Point(1, 1), Point(2, 1), Point(2, 2), Point(1, 2)};
Triangulation::HalfEdges edges1 = {{1, 3}};
std::vector<Vec3i> indices1 = Triangulation::triangulate(points, edges1);
auto check = [](int i1, int i2, Vec3i t)->bool { return true;
return (t[0] == i1 || t[1] == i1 || t[2] == i1) &&
(t[0] == i2 || t[1] == i2 || t[2] == i2);
};
REQUIRE(indices1.size() == 2);
int i1 = edges1.begin()->first,
i2 = edges1.begin()->second;
CHECK(check(i1, i2, indices1[0]));
CHECK(check(i1, i2, indices1[1]));
Triangulation::HalfEdges edges2 = {{0, 2}};
std::vector<Vec3i> indices2 = Triangulation::triangulate(points, edges2);
REQUIRE(indices2.size() == 2);
i1 = edges2.begin()->first;
i2 = edges2.begin()->second;
CHECK(check(i1, i2, indices2[0]));
CHECK(check(i1, i2, indices2[1]));
}

242
tests/libslic3r/test_quadric_edge_collapse.cpp Normal file
View file

@ -0,0 +1,242 @@
#include <catch2/catch.hpp>
#include <test_utils.hpp>
#include <libslic3r/QuadricEdgeCollapse.hpp>
#include <libslic3r/TriangleMesh.hpp> // its - indexed_triangle_set
#include <libslic3r/SimplifyMesh.hpp> // no priority queue
#include "libslic3r/AABBTreeIndirect.hpp" // is similar
using namespace Slic3r;
namespace Private {
struct Similarity
{
float max_distance = 0.f;
float average_distance = 0.f;
Similarity() = default;
Similarity(float max_distance, float average_distance)
: max_distance(max_distance), average_distance(average_distance)
{}
};
// border for our algorithm with frog_leg model and decimation to 5%
Similarity frog_leg_5(0.32f, 0.043f);
Similarity get_similarity(const indexed_triangle_set &from,
const indexed_triangle_set &to)
{
// create ABBTree
auto tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
from.vertices, from.indices);
float sum_distance = 0.f;
float max_distance = 0.f;
auto collect_distances = [&](const Vec3f &surface_point) {
size_t hit_idx;
Vec3f hit_point;
float distance2 =
AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
from.vertices, from.indices, tree, surface_point, hit_idx,
hit_point);
float distance = sqrt(distance2);
if (max_distance < distance) max_distance = distance;
sum_distance += distance;
};
for (const Vec3f &vertex : to.vertices) { collect_distances(vertex); }
for (const Vec3i &t : to.indices) {
Vec3f center(0, 0, 0);
for (size_t i = 0; i < 3; ++i) { center += to.vertices[t[i]] / 3; }
collect_distances(center);
}
size_t count = to.vertices.size() + to.indices.size();
float average_distance = sum_distance / count;
std::cout << "max_distance = " << max_distance << ", average_distance = " << average_distance << std::endl;
return Similarity(max_distance, average_distance);
}
void is_better_similarity(const indexed_triangle_set &its_first,
const indexed_triangle_set &its_second,
const Similarity & compare)
{
Similarity s1 = get_similarity(its_first, its_second);
Similarity s2 = get_similarity(its_second, its_first);
CHECK(s1.average_distance < compare.average_distance);
CHECK(s1.max_distance < compare.max_distance);
CHECK(s2.average_distance < compare.average_distance);
CHECK(s2.max_distance < compare.max_distance);
}
void is_worse_similarity(const indexed_triangle_set &its_first,
const indexed_triangle_set &its_second,
const Similarity & compare)
{
Similarity s1 = get_similarity(its_first, its_second);
Similarity s2 = get_similarity(its_second, its_first);
if (s1.max_distance < compare.max_distance &&
s2.max_distance < compare.max_distance)
CHECK(false);
}
bool exist_triangle_with_twice_vertices(const std::vector<stl_triangle_vertex_indices> &indices)
{
for (const auto &face : indices)
if (face[0] == face[1] || face[0] == face[2] || face[1] == face[2])
return true;
return false;
}
} // namespace Private
TEST_CASE("Reduce one edge by Quadric Edge Collapse", "[its]")
{
indexed_triangle_set its;
its.vertices = {Vec3f(-1.f, 0.f, 0.f), Vec3f(0.f, 1.f, 0.f),
Vec3f(1.f, 0.f, 0.f), Vec3f(0.f, 0.f, 1.f),
// vertex to be removed
Vec3f(0.9f, .1f, -.1f)};
its.indices = {Vec3i(1, 0, 3), Vec3i(2, 1, 3), Vec3i(0, 2, 3),
Vec3i(0, 1, 4), Vec3i(1, 2, 4), Vec3i(2, 0, 4)};
// edge to remove is between vertices 2 and 4 on trinagles 4 and 5
indexed_triangle_set its_ = its; // copy
// its_write_obj(its, "tetrhedron_in.obj");
uint32_t wanted_count = its.indices.size() - 1;
its_quadric_edge_collapse(its, wanted_count);
// its_write_obj(its, "tetrhedron_out.obj");
CHECK(its.indices.size() == 4);
CHECK(its.vertices.size() == 4);
for (size_t i = 0; i < 3; i++) {
CHECK(its.indices[i] == its_.indices[i]);
}
for (size_t i = 0; i < 4; i++) {
if (i == 2) continue;
CHECK(its.vertices[i] == its_.vertices[i]);
}
const Vec3f &v = its.vertices[2]; // new vertex
const Vec3f &v2 = its_.vertices[2]; // moved vertex
const Vec3f &v4 = its_.vertices[4]; // removed vertex
for (size_t i = 0; i < 3; i++) {
bool is_between = (v[i] < v4[i] && v[i] > v2[i]) ||
(v[i] > v4[i] && v[i] < v2[i]);
CHECK(is_between);
}
Private::Similarity max_similarity(0.75f, 0.014f);
Private::is_better_similarity(its, its_, max_similarity);
}
TEST_CASE("Simplify frog_legs.obj to 5% by Quadric edge collapse", "[its][quadric_edge_collapse]")
{
TriangleMesh mesh = load_model("frog_legs.obj");
double original_volume = its_volume(mesh.its);
uint32_t wanted_count = mesh.its.indices.size() * 0.05;
REQUIRE_FALSE(mesh.empty());
indexed_triangle_set its = mesh.its; // copy
float max_error = std::numeric_limits<float>::max();
its_quadric_edge_collapse(its, wanted_count, &max_error);
// its_write_obj(its, "frog_legs_qec.obj");
CHECK(its.indices.size() <= wanted_count);
double volume = its_volume(its);
CHECK(fabs(original_volume - volume) < 33.);
Private::is_better_similarity(mesh.its, its, Private::frog_leg_5);
}
#include <libigl/igl/qslim.h>
#include "Simplify.h"
TEST_CASE("Simplify frog_legs.obj to 5% by IGL/qslim", "[]")
{
std::string obj_filename = "frog_legs.obj";
TriangleMesh mesh = load_model(obj_filename);
REQUIRE_FALSE(mesh.empty());
indexed_triangle_set &its = mesh.its;
double original_volume = its_volume(its);
uint32_t wanted_count = its.indices.size() * 0.05;
Eigen::MatrixXd V(its.vertices.size(), 3);
Eigen::MatrixXi F(its.indices.size(), 3);
for (size_t j = 0; j < its.vertices.size(); ++j) {
Vec3d vd = its.vertices[j].cast<double>();
for (int i = 0; i < 3; ++i) V(j, i) = vd(i);
}
for (size_t j = 0; j < its.indices.size(); ++j) {
const auto &f = its.indices[j];
for (int i = 0; i < 3; ++i) F(j, i) = f(i);
}
size_t max_m = wanted_count;
Eigen::MatrixXd U;
Eigen::MatrixXi G;
Eigen::VectorXi J, I;
CHECK(igl::qslim(V, F, max_m, U, G, J, I));
// convert to its
indexed_triangle_set its_out;
its_out.vertices.reserve(U.size()/3);
its_out.indices.reserve(G.size()/3);
for (size_t i = 0; i < U.size()/3; i++)
its_out.vertices.emplace_back(U(i, 0), U(i, 1), U(i, 2));
for (size_t i = 0; i < G.size()/3; i++)
its_out.indices.emplace_back(G(i, 0), G(i, 1), G(i, 2));
// check if algorithm is still worse than our
Private::is_worse_similarity(its_out, its, Private::frog_leg_5);
// its_out, its --> avg_distance: 0.0351217, max_distance 0.364316
// its, its_out --> avg_distance: 0.0412358, max_distance 0.238913
}
TEST_CASE("Simplify frog_legs.obj to 5% by simplify", "[]") {
std::string obj_filename = "frog_legs.obj";
TriangleMesh mesh = load_model(obj_filename);
uint32_t wanted_count = mesh.its.indices.size() * 0.05;
Simplify::load_obj((TEST_DATA_DIR PATH_SEPARATOR + obj_filename).c_str());
Simplify::simplify_mesh(wanted_count, 5, true);
// convert to its
indexed_triangle_set its_out;
its_out.vertices.reserve(Simplify::vertices.size());
its_out.indices.reserve(Simplify::triangles.size());
for (size_t i = 0; i < Simplify::vertices.size(); i++) {
const Simplify::Vertex &v = Simplify::vertices[i];
its_out.vertices.emplace_back(v.p.x, v.p.y, v.p.z);
}
for (size_t i = 0; i < Simplify::triangles.size(); i++) {
const Simplify::Triangle &t = Simplify::triangles[i];
its_out.indices.emplace_back(t.v[0], t.v[1], t.v[2]);
}
// check if algorithm is still worse than our
Private::is_worse_similarity(its_out, mesh.its, Private::frog_leg_5);
// its_out, mesh.its --> max_distance = 0.700494, average_distance = 0.0902524
// mesh.its, its_out --> max_distance = 0.393184, average_distance = 0.0537392
}
TEST_CASE("Simplify trouble case", "[its]")
{
TriangleMesh tm = load_model("simplification.obj");
REQUIRE_FALSE(tm.empty());
float max_error = std::numeric_limits<float>::max();
uint32_t wanted_count = 0;
its_quadric_edge_collapse(tm.its, wanted_count, &max_error);
CHECK(!Private::exist_triangle_with_twice_vertices(tm.its.indices));
}
TEST_CASE("Simplified cube should not be empty.", "[its]")
{
auto its = its_make_cube(1, 2, 3);
float max_error = std::numeric_limits<float>::max();
uint32_t wanted_count = 0;
its_quadric_edge_collapse(its, wanted_count, &max_error);
CHECK(!its.indices.empty());
}

136
tests/libslic3r/test_triangulation.cpp Normal file
View file

@ -0,0 +1,136 @@
#include <catch2/catch.hpp>
#include <libslic3r/Triangulation.hpp>
#include <libslic3r/SVG.hpp> // only debug visualization
using namespace Slic3r;
namespace Private{
void store_trinagulation(const ExPolygons & shape,
const std::vector<Vec3i> &triangles,
const char* file_name = "Triangulation_visualization.svg",
double scale = 1e5)
{
BoundingBox bb;
for (const auto &expoly : shape) bb.merge(expoly.contour.points);
bb.scale(scale);
SVG svg_vis(file_name, bb);
svg_vis.draw(shape, "gray", .7f);
size_t count = count_points(shape);
Points points;
points.reserve(count);
auto insert_point = [](const Slic3r::Polygon &polygon, Points &points) {
for (const Point &p : polygon.points) points.emplace_back(p);
};
for (const ExPolygon &expolygon : shape) {
insert_point(expolygon.contour, points);
for (const Slic3r::Polygon &hole : expolygon.holes)
insert_point(hole, points);
}
for (const auto &t : triangles) {
Polygon triangle({points[t[0]], points[t[1]], points[t[2]]});
triangle.scale(scale);
svg_vis.draw(triangle, "green");
}
// prevent visualization in test
CHECK(false);
}
} // namespace
TEST_CASE("Triangulate rectangle with restriction on edge", "[Triangulation]")
{
// 0 1 2 3
Points points = {Point(1, 1), Point(2, 1), Point(2, 2), Point(1, 2)};
Triangulation::HalfEdges edges1 = {{1, 3}};
std::vector<Vec3i> indices1 = Triangulation::triangulate(points, edges1, true);
auto check = [](int i1, int i2, Vec3i t) -> bool {
return true;
return (t[0] == i1 || t[1] == i1 || t[2] == i1) &&
(t[0] == i2 || t[1] == i2 || t[2] == i2);
};
REQUIRE(indices1.size() == 2);
int i1 = edges1.begin()->first, i2 = edges1.begin()->second;
CHECK(check(i1, i2, indices1[0]));
CHECK(check(i1, i2, indices1[1]));
Triangulation::HalfEdges edges2 = {{0, 2}};
std::vector<Vec3i> indices2 = Triangulation::triangulate(points, edges2, true);
REQUIRE(indices2.size() == 2);
i1 = edges2.begin()->first;
i2 = edges2.begin()->second;
CHECK(check(i1, i2, indices2[0]));
CHECK(check(i1, i2, indices2[1]));
}
TEST_CASE("Triangulation polygon", "[triangulation]")
{
Points points = {Point(416, 346), Point(445, 362), Point(463, 389),
Point(469, 427), Point(445, 491)};
Polygon polygon(points);
Polygons polygons({polygon});
ExPolygon expolygon(points);
ExPolygons expolygons({expolygon});
std::vector<Vec3i> tp = Triangulation::triangulate(polygon);
std::vector<Vec3i> tps = Triangulation::triangulate(polygons);
std::vector<Vec3i> tep = Triangulation::triangulate(expolygon);
std::vector<Vec3i> teps = Triangulation::triangulate(expolygons);
//Private::store_trinagulation(expolygons, teps);
CHECK(tp.size() == tps.size());
CHECK(tep.size() == teps.size());
CHECK(tp.size() == tep.size());
CHECK(tp.size() == 3);
}
TEST_CASE("Triangulation M shape polygon", "[triangulation]")
{
// 0 1 2 3 4
Polygon shape_M = {Point(0, 0), Point(2, 0), Point(2, 2), Point(1, 1), Point(0, 2)};
std::vector<Vec3i> triangles = Triangulation::triangulate(shape_M);
// Check outer triangle is not contain
std::set<int> outer_triangle = {2, 3, 4};
bool is_in = false;
for (const Vec3i &t : triangles) {
for (size_t i = 0; i < 3; i++) {
int index = t[i];
if (outer_triangle.find(index) == outer_triangle.end()) {
is_in = false;
break;
} else {
is_in = true;
}
}
if (is_in) break;
}
//Private::store_trinagulation({ExPolygon(shape_M)}, triangles);
CHECK(triangles.size() == 3);
CHECK(!is_in);
}
// same point in triangulation are not Supported
//TEST_CASE("Triangulation 2 polygons with same point", "[triangulation][not supported]")
//{
// Slic3r::Polygon polygon1 = {
// Point(416, 346), Point(445, 362),
// Point(463, 389), Point(469, 427) /* This point */,
// Point(445, 491)
// };
// Slic3r::Polygon polygon2 = {
// Point(495, 488), Point(469, 427) /* This point */,
// Point(495, 364)
// };
// ExPolygons shape2d = {ExPolygon(polygon1), ExPolygon(polygon2)};
// std::vector<Vec3i> shape_triangles = Triangulation::triangulate(shape2d);
// store_trinagulation(shape2d, shape_triangles);
//}