Improved robustness of the cut algorithm, follow-up to 96ab500a13

1) Fixed crack between the trimmed object triangles and the triangles
   closing the cut (newly triangulated) by snapping the existing vertices
   on cutting plane to coord_t scaled coordinates. Such snapping is
   needed for vertices close to zero where float has a higher absolute
   accuracy than coord_t.
2) Improved accuracy of triangle cutting by calculating the cut contour with doubles.
3) Improved accuracy of triangle cutting by implementing rounding to coord_t instead of floor.
This commit is contained in:
Vojtech Bubnik 2022-12-12 10:18:17 +01:00
parent 2bcb62d447
commit cda29fa4ac
2 changed files with 76 additions and 47 deletions

View File

@ -112,7 +112,7 @@ inline double angle(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBas
template<typename Derived>
Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign> to_2d(const Eigen::MatrixBase<Derived> &ptN) {
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) >= 3, "to_2d(): first parameter is not a 3D or higher dimensional vector");
return { ptN.x(), ptN.y() };
return ptN.head<2>();
}
template<typename Derived>

View File

@ -137,18 +137,41 @@ enum class FacetSliceType {
Cutting = 2
};
// Convert an int32_t scaled coordinate into an unscaled 3D floating point coordinate (mesh vertex).
template<typename T>
inline Vec3f contour_point_to_v3f(const Point &pt, const T z)
{
return to_3d(
// unscale using doubles for higher accuracy
unscaled<double>(pt).
// then convert to floats
cast<float>(),
float(z));
}
// Convert 2D projection of an int32_t scaled coordinate into an unscaled 3D floating point coordinate (mesh vertex).
template<typename Derived>
inline Point v3f_scaled_to_contour_point(const Eigen::MatrixBase<Derived> &v)
{
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) >= 2, "v3f_scaled_to_contour_point(): Not a 2D or 3D vector.");
using T = typename Derived::Scalar;
return { coord_t(std::floor(v.x() + T(0.5))), coord_t(std::floor(v.y() + T(0.5))) };
}
// Return true, if the facet has been sliced and line_out has been filled.
static FacetSliceType slice_facet(
template<typename T>
inline FacetSliceType slice_facet(
// Z height of the slice in XY plane. Scaled or unscaled (same as vertices[].z()).
float slice_z,
T slice_z,
// 3 vertices of the triangle, XY scaled. Z scaled or unscaled (same as slice_z).
const stl_vertex *vertices,
const Eigen::Matrix<T, 3, 1, Eigen::DontAlign> *vertices,
const stl_triangle_vertex_indices &indices,
const Vec3i &edge_ids,
const int idx_vertex_lowest,
const bool horizontal,
IntersectionLine &line_out)
{
using Vector = Eigen::Matrix<T, 3, 1, Eigen::DontAlign>;
IntersectionPoint points[3];
size_t num_points = 0;
auto point_on_layer = size_t(-1);
@ -158,7 +181,7 @@ static FacetSliceType slice_facet(
// (external on the right of the line)
for (int j = 0; j < 3; ++ j) { // loop through facet edges
int edge_id;
const stl_vertex *a, *b;
const Vector *a, *b;
int a_id, b_id;
{
int k = (idx_vertex_lowest + j) % 3;
@ -174,16 +197,16 @@ static FacetSliceType slice_facet(
if (a->z() == slice_z && b->z() == slice_z) {
// Edge is horizontal and belongs to the current layer.
// The following rotation of the three vertices may not be efficient, but this branch happens rarely.
const stl_vertex &v0 = vertices[0];
const stl_vertex &v1 = vertices[1];
const stl_vertex &v2 = vertices[2];
const Vector &v0 = vertices[0];
const Vector &v1 = vertices[1];
const Vector &v2 = vertices[2];
// We may ignore this edge for slicing purposes, but we may still use it for object cutting.
FacetSliceType result = FacetSliceType::Slicing;
if (horizontal) {
// All three vertices are aligned with slice_z.
line_out.edge_type = IntersectionLine::FacetEdgeType::Horizontal;
result = FacetSliceType::Cutting;
double normal = (v1.x() - v0.x()) * (v2.y() - v1.y()) - (v1.y() - v0.y()) * (v2.x() - v1.x());
double normal = cross2((to_2d(v1) - to_2d(v0)).cast<double>(), (to_2d(v2) - to_2d(v1)).cast<double>());
if (normal < 0) {
// If normal points downwards this is a bottom horizontal facet so we reverse its point order.
std::swap(a, b);
@ -205,10 +228,8 @@ static FacetSliceType slice_facet(
} else
line_out.edge_type = IntersectionLine::FacetEdgeType::Bottom;
}
line_out.a.x() = a->x();
line_out.a.y() = a->y();
line_out.b.x() = b->x();
line_out.b.y() = b->y();
line_out.a = v3f_scaled_to_contour_point(*a);
line_out.b = v3f_scaled_to_contour_point(*b);
line_out.a_id = a_id;
line_out.b_id = b_id;
assert(line_out.a != line_out.b);
@ -220,8 +241,7 @@ static FacetSliceType slice_facet(
if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != a_id) {
point_on_layer = num_points;
IntersectionPoint &point = points[num_points ++];
point.x() = a->x();
point.y() = a->y();
static_cast<Point&>(point) = v3f_scaled_to_contour_point(*a);
point.point_id = a_id;
}
} else if (b->z() == slice_z) {
@ -229,8 +249,7 @@ static FacetSliceType slice_facet(
if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != b_id) {
point_on_layer = num_points;
IntersectionPoint &point = points[num_points ++];
point.x() = b->x();
point.y() = b->y();
static_cast<Point&>(point) = v3f_scaled_to_contour_point(*b);
point.point_id = b_id;
}
} else if ((a->z() < slice_z && b->z() > slice_z) || (b->z() < slice_z && a->z() > slice_z)) {
@ -270,16 +289,10 @@ static FacetSliceType slice_facet(
}
#else
// Just clamp the intersection point to source triangle edge.
if (t <= 0.) {
point.x() = a->x();
point.y() = a->y();
} else if (t >= 1.) {
point.x() = b->x();
point.y() = b->y();
} else {
point.x() = coord_t(floor(double(a->x()) + (double(b->x()) - double(a->x())) * t + 0.5));
point.y() = coord_t(floor(double(a->y()) + (double(b->y()) - double(a->y())) * t + 0.5));
}
static_cast<Point&>(point) =
t <= 0. ? v3f_scaled_to_contour_point(*a) :
t >= 1. ? v3f_scaled_to_contour_point(*b) :
v3f_scaled_to_contour_point(a->head<2>().cast<double>() * (1. - t) + b->head<2>().cast<double>() * t + Vec2d(0.5, 0.5));
point.edge_id = edge_id;
++ num_points;
#endif
@ -2182,28 +2195,44 @@ void cut_mesh(const indexed_triangle_set &mesh, float z, indexed_triangle_set *u
int idx_vertex_lowest = (vertices[1].z() == min_z) ? 1 : ((vertices[2].z() == min_z) ? 2 : 0);
FacetSliceType slice_type = FacetSliceType::NoSlice;
if (z > min_z - EPSILON && z < max_z + EPSILON) {
Vec3f vertices_scaled[3];
Vec3d vertices_scaled[3];
for (int i = 0; i < 3; ++ i) {
const Vec3f &src = vertices[i];
Vec3f &dst = vertices_scaled[i];
dst.x() = scale_(src.x());
dst.y() = scale_(src.y());
Vec3d &dst = vertices_scaled[i];
dst.x() = scaled<double>(src.x());
dst.y() = scaled<double>(src.y());
dst.z() = src.z();
}
slice_type = slice_facet(z, vertices_scaled, mesh.indices[facet_idx], facets_edge_ids[facet_idx], idx_vertex_lowest, min_z == max_z, line);
slice_type = slice_facet(double(z), vertices_scaled, mesh.indices[facet_idx], facets_edge_ids[facet_idx], idx_vertex_lowest, min_z == max_z, line);
}
if (slice_type != FacetSliceType::NoSlice) {
// Save intersection lines for generating correct triangulations.
if (line.edge_type == IntersectionLine::FacetEdgeType::Top) {
lower_lines.emplace_back(line);
lower_slice_vertices.emplace_back(line.a_id);
lower_slice_vertices.emplace_back(line.b_id);
if (lower) {
lower_lines.emplace_back(line);
if (triangulate_caps) {
// Snap these vertices to coord_t grid, so that they will be matched with the vertices produced
// by triangulating opening on the cut.
lower->vertices[line.a_id] = contour_point_to_v3f(line.a, z);
lower->vertices[line.b_id] = contour_point_to_v3f(line.b, z);
}
}
} else if (line.edge_type == IntersectionLine::FacetEdgeType::Bottom) {
upper_lines.emplace_back(line);
upper_slice_vertices.emplace_back(line.a_id);
upper_slice_vertices.emplace_back(line.b_id);
} else if (line.edge_type == IntersectionLine::FacetEdgeType::General) {
if (upper) {
upper_lines.emplace_back(line);
if (triangulate_caps) {
// Snap these vertices to coord_t grid, so that they will be matched with the vertices produced
// by triangulating opening on the cut.
upper->vertices[line.a_id] = contour_point_to_v3f(line.a, z);
upper->vertices[line.b_id] = contour_point_to_v3f(line.b, z);
}
}
} else if (line.edge_type == IntersectionLine::FacetEdgeType::General && triangulate_caps) {
lower_lines.emplace_back(line);
upper_lines.emplace_back(line);
}
@ -2235,11 +2264,11 @@ void cut_mesh(const indexed_triangle_set &mesh, float z, indexed_triangle_set *u
assert(facets_edge_ids[facet_idx](iv) == line.edge_a_id || facets_edge_ids[facet_idx](iv) == line.edge_b_id);
if (facets_edge_ids[facet_idx](iv) == line.edge_a_id) {
// Unscale to doubles first, then to floats to reach the same accuracy as triangulate_expolygons_2d().
v0v1 = to_3d(unscaled<double>(line.a).cast<float>().eval(), z);
v2v0 = to_3d(unscaled<double>(line.b).cast<float>().eval(), z);
v0v1 = contour_point_to_v3f(line.a, z);
v2v0 = contour_point_to_v3f(line.b, z);
} else {
v0v1 = to_3d(unscaled<double>(line.b).cast<float>().eval(), z);
v2v0 = to_3d(unscaled<double>(line.a).cast<float>().eval(), z);
v0v1 = contour_point_to_v3f(line.b, z);
v2v0 = contour_point_to_v3f(line.a, z);
}
const stl_vertex &v0 = vertices[iv];
const int iv0 = facet[iv];