3D-printing/OrcaSlicer
1
0
Fork 0
mirror of https://github.com/SoftFever/OrcaSlicer.git synced 2025-07-08 23:46:24 -06:00
Code Issues Projects Releases Packages Wiki Activity Actions 3

ENH: new seam strategy from prusa2.5

Browse source 
As title. Thanks @Prusa

Signed-off-by: salt.wei <salt.wei@bambulab.com>
Change-Id: I2fa177e27ac53211952ea9b6c62e98182b8f05ce
This commit is contained in:
salt.wei 2022-08-18 14:17:02 +08:00 • committed by Lane.Wei
parent ce082f6e2a
commit d73142c2f9
23 changed files with 3105 additions and 1323 deletions
Download patch file Download diff file Expand all files Collapse all files

160
src/libslic3r/AABBTreeIndirect.hpp
View file

@ -446,6 +446,57 @@ namespace detail {
} }
} }
// Real-time collision detection, Ericson, Chapter 5
template<typename Vector>
static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c)
{
using Scalar = typename Vector::Scalar;
// Check if P in vertex region outside A
Vector ab = b - a;
Vector ac = c - a;
Vector ap = p - a;
Scalar d1 = ab.dot(ap);
Scalar d2 = ac.dot(ap);
if (d1 <= 0 && d2 <= 0)
return a;
// Check if P in vertex region outside B
Vector bp = p - b;
Scalar d3 = ab.dot(bp);
Scalar d4 = ac.dot(bp);
if (d3 >= 0 && d4 <= d3)
return b;
// Check if P in edge region of AB, if so return projection of P onto AB
Scalar vc = d1*d4 - d3*d2;
if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) {
Scalar v = d1 / (d1 - d3);
return a + v * ab;
}
// Check if P in vertex region outside C
Vector cp = p - c;
Scalar d5 = ab.dot(cp);
Scalar d6 = ac.dot(cp);
if (d6 >= 0 && d5 <= d6)
return c;
// Check if P in edge region of AC, if so return projection of P onto AC
Scalar vb = d5*d2 - d1*d6;
if (vb <= 0 && d2 >= 0 && d6 <= 0) {
Scalar w = d2 / (d2 - d6);
return a + w * ac;
}
// Check if P in edge region of BC, if so return projection of P onto BC
Scalar va = d3*d6 - d5*d4;
if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
return b + w * (c - b);
}
// P inside face region. Compute Q through its barycentric coordinates (u,v,w)
Scalar denom = Scalar(1.0) / (va + vb + vc);
Scalar v = vb * denom;
Scalar w = vc * denom;
return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
};
// Nothing to do with COVID-19 social distancing. // Nothing to do with COVID-19 social distancing.
template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType> template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
struct IndexedTriangleSetDistancer { struct IndexedTriangleSetDistancer {
@ -453,74 +504,36 @@ namespace detail {
using IndexedFaceType = AIndexedFaceType; using IndexedFaceType = AIndexedFaceType;
using TreeType = ATreeType; using TreeType = ATreeType;
using VectorType = AVectorType; using VectorType = AVectorType;
using ScalarType = typename VectorType::Scalar;
const std::vector<VertexType> &vertices; const std::vector<VertexType> &vertices;
const std::vector<IndexedFaceType> &faces; const std::vector<IndexedFaceType> &faces;
const TreeType &tree; const TreeType &tree;
const VectorType origin; const VectorType origin;
inline VectorType closest_point_to_origin(size_t primitive_index,
ScalarType& squared_distance){
const auto &triangle = this->faces[primitive_index];
VectorType closest_point = closest_point_to_triangle<VectorType>(origin,
this->vertices[triangle(0)].template cast<ScalarType>(),
this->vertices[triangle(1)].template cast<ScalarType>(),
this->vertices[triangle(2)].template cast<ScalarType>());
squared_distance = (origin - closest_point).squaredNorm();
return closest_point;
}
}; };
// Real-time collision detection, Ericson, Chapter 5 template<typename IndexedPrimitivesDistancerType, typename Scalar>
template<typename Vector> static inline Scalar squared_distance_to_indexed_primitives_recursive(
static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c) IndexedPrimitivesDistancerType &distancer,
{
using Scalar = typename Vector::Scalar;
// Check if P in vertex region outside A
Vector ab = b - a;
Vector ac = c - a;
Vector ap = p - a;
Scalar d1 = ab.dot(ap);
Scalar d2 = ac.dot(ap);
if (d1 <= 0 && d2 <= 0)
return a;
// Check if P in vertex region outside B
Vector bp = p - b;
Scalar d3 = ab.dot(bp);
Scalar d4 = ac.dot(bp);
if (d3 >= 0 && d4 <= d3)
return b;
// Check if P in edge region of AB, if so return projection of P onto AB
Scalar vc = d1*d4 - d3*d2;
if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) {
Scalar v = d1 / (d1 - d3);
return a + v * ab;
}
// Check if P in vertex region outside C
Vector cp = p - c;
Scalar d5 = ab.dot(cp);
Scalar d6 = ac.dot(cp);
if (d6 >= 0 && d5 <= d6)
return c;
// Check if P in edge region of AC, if so return projection of P onto AC
Scalar vb = d5*d2 - d1*d6;
if (vb <= 0 && d2 >= 0 && d6 <= 0) {
Scalar w = d2 / (d2 - d6);
return a + w * ac;
}
// Check if P in edge region of BC, if so return projection of P onto BC
Scalar va = d3*d6 - d5*d4;
if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
return b + w * (c - b);
}
// P inside face region. Compute Q through its barycentric coordinates (u,v,w)
Scalar denom = Scalar(1.0) / (va + vb + vc);
Scalar v = vb * denom;
Scalar w = vc * denom;
return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
};
template<typename IndexedTriangleSetDistancerType, typename Scalar>
static inline Scalar squared_distance_to_indexed_triangle_set_recursive(
IndexedTriangleSetDistancerType &distancer,
size_t node_idx, size_t node_idx,
Scalar low_sqr_d, Scalar low_sqr_d,
Scalar up_sqr_d, Scalar up_sqr_d,
size_t &i, size_t &i,
Eigen::PlainObjectBase<typename IndexedTriangleSetDistancerType::VectorType> &c) Eigen::PlainObjectBase<typename IndexedPrimitivesDistancerType::VectorType> &c)
{ {
using Vector = typename IndexedTriangleSetDistancerType::VectorType; using Vector = typename IndexedPrimitivesDistancerType::VectorType;
if (low_sqr_d > up_sqr_d) if (low_sqr_d > up_sqr_d)
return low_sqr_d; return low_sqr_d;
@ -538,13 +551,9 @@ namespace detail {
assert(node.is_valid()); assert(node.is_valid());
if (node.is_leaf()) if (node.is_leaf())
{ {
const auto &triangle = distancer.faces[node.idx]; Scalar sqr_dist;
Vector c_candidate = closest_point_to_triangle<Vector>( Vector c_candidate = distancer.closest_point_to_origin(node.idx, sqr_dist);
distancer.origin, set_min(sqr_dist, node.idx, c_candidate);
distancer.vertices[triangle(0)].template cast<Scalar>(),
distancer.vertices[triangle(1)].template cast<Scalar>(),
distancer.vertices[triangle(2)].template cast<Scalar>());
set_min((c_candidate - distancer.origin).squaredNorm(), node.idx, c_candidate);
} }
else else
{ {
@ -561,7 +570,7 @@ namespace detail {
{ {
size_t i_left; size_t i_left;
Vector c_left = c; Vector c_left = c;
Scalar sqr_d_left = squared_distance_to_indexed_triangle_set_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left); Scalar sqr_d_left = squared_distance_to_indexed_primitives_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left);
set_min(sqr_d_left, i_left, c_left); set_min(sqr_d_left, i_left, c_left);
looked_left = true; looked_left = true;
}; };
@ -569,13 +578,13 @@ namespace detail {
{ {
size_t i_right; size_t i_right;
Vector c_right = c; Vector c_right = c;
Scalar sqr_d_right = squared_distance_to_indexed_triangle_set_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right); Scalar sqr_d_right = squared_distance_to_indexed_primitives_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right);
set_min(sqr_d_right, i_right, c_right); set_min(sqr_d_right, i_right, c_right);
looked_right = true; looked_right = true;
}; };
// must look left or right if in box // must look left or right if in box
using BBoxScalar = typename IndexedTriangleSetDistancerType::TreeType::BoundingBox::Scalar; using BBoxScalar = typename IndexedPrimitivesDistancerType::TreeType::BoundingBox::Scalar;
if (node_left.bbox.contains(distancer.origin.template cast<BBoxScalar>())) if (node_left.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
look_left(); look_left();
if (node_right.bbox.contains(distancer.origin.template cast<BBoxScalar>())) if (node_right.bbox.contains(distancer.origin.template cast<BBoxScalar>()))
@ -709,10 +718,15 @@ inline bool intersect_ray_all_hits(
origin, dir, VectorType(dir.cwiseInverse()), origin, dir, VectorType(dir.cwiseInverse()),
eps } eps }
}; };
if (! tree.empty()) { if (tree.empty()) {
hits.clear();
} else {
// Reusing the output memory if there is some memory already pre-allocated.
ray_intersector.hits = std::move(hits);
ray_intersector.hits.clear();
ray_intersector.hits.reserve(8); ray_intersector.hits.reserve(8);
detail::intersect_ray_recursive_all_hits(ray_intersector, 0); detail::intersect_ray_recursive_all_hits(ray_intersector, 0);
std::swap(hits, ray_intersector.hits); hits = std::move(ray_intersector.hits);
std::sort(hits.begin(), hits.end(), [](const auto &l, const auto &r) { return l.t < r.t; }); std::sort(hits.begin(), hits.end(), [](const auto &l, const auto &r) { return l.t < r.t; });
} }
return ! hits.empty(); return ! hits.empty();
@ -742,7 +756,7 @@ inline typename VectorType::Scalar squared_distance_to_indexed_triangle_set(
auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType> auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
{ vertices, faces, tree, point }; { vertices, faces, tree, point };
return tree.empty() ? Scalar(-1) : return tree.empty() ? Scalar(-1) :
detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out); detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
} }
// Decides if exists some triangle in defined radius on a 3D indexed triangle set using a pre-built AABBTreeIndirect::Tree. // Decides if exists some triangle in defined radius on a 3D indexed triangle set using a pre-built AABBTreeIndirect::Tree.
@ -759,22 +773,22 @@ inline bool is_any_triangle_in_radius(
const TreeType &tree, const TreeType &tree,
// Point to which the closest point on the indexed triangle set is searched for. // Point to which the closest point on the indexed triangle set is searched for.
const VectorType &point, const VectorType &point,
// Maximum distance in which triangle is search for //Square of maximum distance in which triangle is searched for
typename VectorType::Scalar &max_distance) typename VectorType::Scalar &max_distance_squared)
{ {
using Scalar = typename VectorType::Scalar; using Scalar = typename VectorType::Scalar;
auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType> auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
{ vertices, faces, tree, point }; { vertices, faces, tree, point };
size_t hit_idx; size_t hit_idx;
VectorType hit_point = VectorType::Ones() * (std::nan("")); VectorType hit_point = VectorType::Ones() * (NaN<typename VectorType::Scalar>);
if(tree.empty()) if(tree.empty())
{ {
return false; return false;
} }
detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance, hit_idx, hit_point); detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
return hit_point.allFinite(); return hit_point.allFinite();
} }

112
src/libslic3r/AABBTreeLines.hpp Normal file
View file

@ -0,0 +1,112 @@
#ifndef SRC_LIBSLIC3R_AABBTREELINES_HPP_
#define SRC_LIBSLIC3R_AABBTREELINES_HPP_
#include "libslic3r/Point.hpp"
#include "libslic3r/EdgeGrid.hpp"
#include "libslic3r/AABBTreeIndirect.hpp"
#include "libslic3r/Line.hpp"
namespace Slic3r {
namespace AABBTreeLines {
namespace detail {
template<typename ALineType, typename ATreeType, typename AVectorType>
struct IndexedLinesDistancer {
using LineType = ALineType;
using TreeType = ATreeType;
using VectorType = AVectorType;
using ScalarType = typename VectorType::Scalar;
const std::vector<LineType> &lines;
const TreeType &tree;
const VectorType origin;
inline VectorType closest_point_to_origin(size_t primitive_index,
ScalarType &squared_distance) {
VectorType nearest_point;
const LineType &line = lines[primitive_index];
squared_distance = line_alg::distance_to_squared(line, origin, &nearest_point);
return nearest_point;
}
};
}
// Build a balanced AABB Tree over a vector of float lines, balancing the tree
// on centroids of the lines.
// Epsilon is applied to the bounding boxes of the AABB Tree to cope with numeric inaccuracies
// during tree traversal.
template<typename LineType>
inline AABBTreeIndirect::Tree<2, typename LineType::Scalar> build_aabb_tree_over_indexed_lines(
const std::vector<LineType> &lines,
//FIXME do we want to apply an epsilon?
const float eps = 0)
{
using TreeType = AABBTreeIndirect::Tree<2, typename LineType::Scalar>;
// using CoordType = typename TreeType::CoordType;
using VectorType = typename TreeType::VectorType;
using BoundingBox = typename TreeType::BoundingBox;
struct InputType {
size_t idx() const {
return m_idx;
}
const BoundingBox& bbox() const {
return m_bbox;
}
const VectorType& centroid() const {
return m_centroid;
}
size_t m_idx;
BoundingBox m_bbox;
VectorType m_centroid;
};
std::vector<InputType> input;
input.reserve(lines.size());
const VectorType veps(eps, eps);
for (size_t i = 0; i < lines.size(); ++i) {
const LineType &line = lines[i];
InputType n;
n.m_idx = i;
n.m_centroid = (line.a + line.b) * 0.5;
n.m_bbox = BoundingBox(line.a, line.a);
n.m_bbox.extend(line.b);
n.m_bbox.min() -= veps;
n.m_bbox.max() += veps;
input.emplace_back(n);
}
TreeType out;
out.build(std::move(input));
return out;
}
// Finding a closest line, its closest point and squared distance to the closest point
// Returns squared distance to the closest point or -1 if the input is empty.
template<typename LineType, typename TreeType, typename VectorType>
inline typename VectorType::Scalar squared_distance_to_indexed_lines(
const std::vector<LineType> &lines,
const TreeType &tree,
const VectorType &point,
size_t &hit_idx_out,
Eigen::PlainObjectBase<VectorType> &hit_point_out)
{
using Scalar = typename VectorType::Scalar;
auto distancer = detail::IndexedLinesDistancer<LineType, TreeType, VectorType>
{ lines, tree, point };
return tree.empty() ?
Scalar(-1) :
AABBTreeIndirect::detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0),
std::numeric_limits<Scalar>::infinity(), hit_idx_out, hit_point_out);
}
}
}
#endif /* SRC_LIBSLIC3R_AABBTREELINES_HPP_ */

8
src/libslic3r/CMakeLists.txt
View file

@ -142,10 +142,12 @@ set(lisbslic3r_sources
GCodeWriter.hpp GCodeWriter.hpp
Geometry.cpp Geometry.cpp
Geometry.hpp Geometry.hpp
Geometry/Bicubic.hpp
Geometry/Circle.cpp Geometry/Circle.cpp
Geometry/Circle.hpp Geometry/Circle.hpp
Geometry/ConvexHull.cpp Geometry/ConvexHull.cpp
Geometry/ConvexHull.hpp Geometry/ConvexHull.hpp
Geometry/Curves.hpp
Geometry/MedialAxis.cpp Geometry/MedialAxis.cpp
Geometry/MedialAxis.hpp Geometry/MedialAxis.hpp
Geometry/Voronoi.hpp Geometry/Voronoi.hpp
@ -176,6 +178,8 @@ set(lisbslic3r_sources
CustomGCode.hpp CustomGCode.hpp
Arrange.hpp Arrange.hpp
Arrange.cpp Arrange.cpp
NormalUtils.cpp
NormalUtils.hpp
Orient.hpp Orient.hpp
Orient.cpp Orient.cpp
MultiPoint.cpp MultiPoint.cpp
@ -222,6 +226,8 @@ set(lisbslic3r_sources
QuadricEdgeCollapse.cpp QuadricEdgeCollapse.cpp
QuadricEdgeCollapse.hpp QuadricEdgeCollapse.hpp
Semver.cpp Semver.cpp
ShortEdgeCollapse.cpp
ShortEdgeCollapse.hpp
ShortestPath.cpp ShortestPath.cpp
ShortestPath.hpp ShortestPath.hpp
SLAPrint.cpp SLAPrint.cpp
@ -264,6 +270,8 @@ set(lisbslic3r_sources
Thread.hpp Thread.hpp
TriangleSelector.cpp TriangleSelector.cpp
TriangleSelector.hpp TriangleSelector.hpp
TriangleSetSampling.cpp
TriangleSetSampling.hpp
MTUtils.hpp MTUtils.hpp
VariableWidth.cpp VariableWidth.cpp
VariableWidth.hpp VariableWidth.hpp

72
src/libslic3r/ExtrusionEntity.cpp
View file

@ -158,11 +158,10 @@ double ExtrusionLoop::length() const
return len; return len;
} }
bool ExtrusionLoop::split_at_vertex(const Point &point) bool ExtrusionLoop::split_at_vertex(const Point &point, const double scaled_epsilon)
{ {
for (ExtrusionPaths::iterator path = this->paths.begin(); path != this->paths.end(); ++path) { for (ExtrusionPaths::iterator path = this->paths.begin(); path != this->paths.end(); ++path) {
int idx = path->polyline.find_point(point); if (int idx = path->polyline.find_point(point, scaled_epsilon); idx != -1) {
if (idx != -1) {
if (this->paths.size() == 1) { if (this->paths.size() == 1) {
// just change the order of points // just change the order of points
Polyline p1, p2; Polyline p1, p2;
@ -207,46 +206,57 @@ bool ExtrusionLoop::split_at_vertex(const Point &point)
return false; return false;
} }
std::pair<size_t, Point> ExtrusionLoop::get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const ExtrusionLoop::ClosestPathPoint ExtrusionLoop::get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const
{ {
// Find the closest path and closest point belonging to that path. Avoid overhangs, if asked for. // Find the closest path and closest point belonging to that path. Avoid overhangs, if asked for.
size_t path_idx = 0; ClosestPathPoint out{0, 0};
Point p; double min2 = std::numeric_limits<double>::max();
{ ClosestPathPoint best_non_overhang{0, 0};
double min = std::numeric_limits<double>::max(); double min2_non_overhang = std::numeric_limits<double>::max();
Point p_non_overhang; for (const ExtrusionPath &path : this->paths) {
size_t path_idx_non_overhang = 0; std::pair<int, Point> foot_pt_ = foot_pt(path.polyline.points, point);
double min_non_overhang = std::numeric_limits<double>::max(); double d2 = (foot_pt_.second - point).cast<double>().squaredNorm();
for (const ExtrusionPath& path : this->paths) { if (d2 < min2) {
Point p_tmp = point.projection_onto(path.polyline); out.foot_pt = foot_pt_.second;
double dist = (p_tmp - point).cast<double>().norm(); out.path_idx = &path - &this->paths.front();
if (dist < min) { out.segment_idx = foot_pt_.first;
p = p_tmp; min2 = d2;
min = dist;
path_idx = &path - &this->paths.front();
}
if (prefer_non_overhang && !is_bridge(path.role()) && dist < min_non_overhang) {
p_non_overhang = p_tmp;
min_non_overhang = dist;
path_idx_non_overhang = &path - &this->paths.front();
}
} }
if (prefer_non_overhang && min_non_overhang != std::numeric_limits<double>::max()) { if (prefer_non_overhang && !is_bridge(path.role()) && d2 < min2_non_overhang) {
// Only apply the non-overhang point if there is one. best_non_overhang.foot_pt = foot_pt_.second;
path_idx = path_idx_non_overhang; best_non_overhang.path_idx = &path - &this->paths.front();
p = p_non_overhang; best_non_overhang.segment_idx = foot_pt_.first;
min2_non_overhang = d2;
} }
} }
return std::make_pair(path_idx, p); if (prefer_non_overhang && min2_non_overhang != std::numeric_limits<double>::max())
// Only apply the non-overhang point if there is one.
out = best_non_overhang;
return out;
} }
// Splitting an extrusion loop, possibly made of multiple segments, some of the segments may be bridging. // Splitting an extrusion loop, possibly made of multiple segments, some of the segments may be bridging.
void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang) void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon)
{ {
if (this->paths.empty()) if (this->paths.empty())
return; return;
auto [path_idx, p] = get_closest_path_and_point(point, prefer_non_overhang); auto [path_idx, segment_idx, p] = get_closest_path_and_point(point, prefer_non_overhang);
// Snap p to start or end of segment_idx if closer than scaled_epsilon.
{
const Point *p1 = this->paths[path_idx].polyline.points.data() + segment_idx;
const Point *p2 = p1;
++p2;
double d2_1 = (point - *p1).cast<double>().squaredNorm();
double d2_2 = (point - *p2).cast<double>().squaredNorm();
const double thr2 = scaled_epsilon * scaled_epsilon;
if (d2_1 < d2_2) {
if (d2_1 < thr2) p = *p1;
} else {
if (d2_2 < thr2) p = *p2;
}
}
// now split path_idx in two parts // now split path_idx in two parts
const ExtrusionPath &path = this->paths[path_idx]; const ExtrusionPath &path = this->paths[path_idx];

12
src/libslic3r/ExtrusionEntity.hpp
View file

@ -314,9 +314,15 @@ public:
const Point& last_point() const override { assert(this->first_point() == this->paths.back().polyline.points.back()); return this->first_point(); } const Point& last_point() const override { assert(this->first_point() == this->paths.back().polyline.points.back()); return this->first_point(); }
Polygon polygon() const; Polygon polygon() const;
double length() const override; double length() const override;
bool split_at_vertex(const Point &point); bool split_at_vertex(const Point &point, const double scaled_epsilon = scaled<double>(0.001));
void split_at(const Point &point, bool prefer_non_overhang); void split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon = scaled<double>(0.001));
std::pair<size_t, Point> get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const; struct ClosestPathPoint
{
size_t path_idx;
size_t segment_idx;
Point foot_pt;
};
ClosestPathPoint get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const;
void clip_end(double distance, ExtrusionPaths* paths) const; void clip_end(double distance, ExtrusionPaths* paths) const;
// Test, whether the point is extruded by a bridging flow. // Test, whether the point is extruded by a bridging flow.
// This used to be used to avoid placing seams on overhangs, but now the EdgeGrid is used instead. // This used to be used to avoid placing seams on overhangs, but now the EdgeGrid is used instead.

46
src/libslic3r/GCode.cpp
View file

@ -1567,7 +1567,8 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
print.throw_if_canceled(); print.throw_if_canceled();
// Collect custom seam data from all objects. // Collect custom seam data from all objects.
m_seam_placer.init(print); std::function<void(void)> throw_if_canceled_func = [&print]() { print.throw_if_canceled(); };
m_seam_placer.init(print, throw_if_canceled_func);
// BBS: priming logic is removed, always set first extruer here. // BBS: priming logic is removed, always set first extruer here.
//if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming)) //if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming))
@ -2792,7 +2793,6 @@ GCode::LayerResult GCode::process_layer(
m_wipe_tower->set_is_first_print(true); m_wipe_tower->set_is_first_print(true);
// Extrude the skirt, brim, support, perimeters, infill ordered by the extruders. // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders.
std::vector<std::unique_ptr<EdgeGrid::Grid>> lower_layer_edge_grids(layers.size());
for (unsigned int extruder_id : layer_tools.extruders) for (unsigned int extruder_id : layer_tools.extruders)
{ {
gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ? gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ?
@ -2981,9 +2981,9 @@ GCode::LayerResult GCode::process_layer(
//This behaviour is same with cura //This behaviour is same with cura
if (is_infill_first && !first_layer) { if (is_infill_first && !first_layer) {
gcode += this->extrude_infill(print, by_region_specific, false); gcode += this->extrude_infill(print, by_region_specific, false);
gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); gcode += this->extrude_perimeters(print, by_region_specific);
} else { } else {
gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); gcode += this->extrude_perimeters(print, by_region_specific);
gcode += this->extrude_infill(print,by_region_specific, false); gcode += this->extrude_infill(print,by_region_specific, false);
} }
// ironing // ironing
@ -3161,14 +3161,11 @@ static std::unique_ptr<EdgeGrid::Grid> calculate_layer_edge_grid(const Layer& la
} }
std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid) std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed)
{ {
// get a copy; don't modify the orientation of the original loop object otherwise // get a copy; don't modify the orientation of the original loop object otherwise
// next copies (if any) would not detect the correct orientation // next copies (if any) would not detect the correct orientation
if (m_layer->lower_layer && lower_layer_edge_grid != nullptr && ! *lower_layer_edge_grid)
*lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
//BBS: extrude contour of wall ccw, hole of wall cw, except spiral mode //BBS: extrude contour of wall ccw, hole of wall cw, except spiral mode
bool was_clockwise = loop.is_clockwise(); bool was_clockwise = loop.is_clockwise();
if (m_config.spiral_mode || !is_perimeter(loop.role())) if (m_config.spiral_mode || !is_perimeter(loop.role()))
@ -3178,17 +3175,13 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
// find the point of the loop that is closest to the current extruder position // find the point of the loop that is closest to the current extruder position
// or randomize if requested // or randomize if requested
Point last_pos = this->last_pos(); Point last_pos = this->last_pos();
if (m_config.spiral_mode) { if (!m_config.spiral_mode && description == "perimeter") {
assert(m_layer != nullptr);
bool is_outer_wall_first = m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill
|| m_config.wall_infill_order == WallInfillOrder::InfillOuterInner;
m_seam_placer.place_seam(m_layer, loop, is_outer_wall_first, this->last_pos());
} else
loop.split_at(last_pos, false); loop.split_at(last_pos, false);
}
else {
//BBS
bool is_outer_wall_first =
m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill ||
m_config.wall_infill_order == WallInfillOrder::InfillOuterInner;
m_seam_placer.place_seam(loop, this->last_pos(), is_outer_wall_first,
EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr);
}
// clip the path to avoid the extruder to get exactly on the first point of the loop; // clip the path to avoid the extruder to get exactly on the first point of the loop;
// if polyline was shorter than the clipping distance we'd get a null polyline, so // if polyline was shorter than the clipping distance we'd get a null polyline, so
@ -3298,14 +3291,14 @@ std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string
return gcode; return gcode;
} }
std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed, std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid) std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed)
{ {
if (const ExtrusionPath* path = dynamic_cast<const ExtrusionPath*>(&entity)) if (const ExtrusionPath* path = dynamic_cast<const ExtrusionPath*>(&entity))
return this->extrude_path(*path, description, speed); return this->extrude_path(*path, description, speed);
else if (const ExtrusionMultiPath* multipath = dynamic_cast<const ExtrusionMultiPath*>(&entity)) else if (const ExtrusionMultiPath* multipath = dynamic_cast<const ExtrusionMultiPath*>(&entity))
return this->extrude_multi_path(*multipath, description, speed); return this->extrude_multi_path(*multipath, description, speed);
else if (const ExtrusionLoop* loop = dynamic_cast<const ExtrusionLoop*>(&entity)) else if (const ExtrusionLoop* loop = dynamic_cast<const ExtrusionLoop*>(&entity))
return this->extrude_loop(*loop, description, speed, lower_layer_edge_grid); return this->extrude_loop(*loop, description, speed);
else else
throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()"); throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()");
return ""; return "";
@ -3327,24 +3320,15 @@ std::string GCode::extrude_path(ExtrusionPath path, std::string description, dou
} }
// Extrude perimeters: Decide where to put seams (hide or align seams). // Extrude perimeters: Decide where to put seams (hide or align seams).
std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid) std::string GCode::extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region)
{ {
std::string gcode; std::string gcode;
for (const ObjectByExtruder::Island::Region &region : by_region) for (const ObjectByExtruder::Island::Region &region : by_region)
if (! region.perimeters.empty()) { if (! region.perimeters.empty()) {
m_config.apply(print.get_print_region(&region - &by_region.front()).config()); m_config.apply(print.get_print_region(&region - &by_region.front()).config());
// plan_perimeters tries to place seams, it needs to have the lower_layer_edge_grid calculated already.
if (m_layer->lower_layer && ! lower_layer_edge_grid)
lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
m_seam_placer.plan_perimeters(std::vector<const ExtrusionEntity*>(region.perimeters.begin(), region.perimeters.end()),
*m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter),
(m_layer == NULL ? nullptr : m_layer->object()),
(lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr));
for (const ExtrusionEntity* ee : region.perimeters) for (const ExtrusionEntity* ee : region.perimeters)
gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid); gcode += this->extrude_entity(*ee, "perimeter", -1.);
} }
return gcode; return gcode;
} }

6
src/libslic3r/GCode.hpp
View file

@ -325,8 +325,8 @@ private:
std::string preamble(); std::string preamble();
// BBS // BBS
std::string change_layer(coordf_t print_z, bool lazy_raise = false); std::string change_layer(coordf_t print_z, bool lazy_raise = false);
std::string extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr); std::string extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1.);
std::string extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1., std::unique_ptr<EdgeGrid::Grid> *lower_layer_edge_grid = nullptr); std::string extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1.);
std::string extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.); std::string extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.);
std::string extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.); std::string extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.);
@ -393,7 +393,7 @@ private:
// For sequential print, the instance of the object to be printing has to be defined. // For sequential print, the instance of the object to be printing has to be defined.
const size_t single_object_instance_idx); const size_t single_object_instance_idx);
std::string extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, std::unique_ptr<EdgeGrid::Grid> &lower_layer_edge_grid); std::string extrude_perimeters(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region);
std::string extrude_infill(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, bool ironing); std::string extrude_infill(const Print &print, const std::vector<ObjectByExtruder::Island::Region> &by_region, bool ironing);
std::string extrude_support(const ExtrusionEntityCollection &support_fills); std::string extrude_support(const ExtrusionEntityCollection &support_fills);

2177
src/libslic3r/GCode/SeamPlacer.cpp
View file

File diff suppressed because it is too large Load diff

257
src/libslic3r/GCode/SeamPlacer.hpp
View file

@ -3,12 +3,16 @@
#include <optional> #include <optional>
#include <vector> #include <vector>
#include <memory>
#include <atomic>
#include "libslic3r/libslic3r.h"
#include "libslic3r/ExtrusionEntity.hpp" #include "libslic3r/ExtrusionEntity.hpp"
#include "libslic3r/Polygon.hpp" #include "libslic3r/Polygon.hpp"
#include "libslic3r/PrintConfig.hpp" #include "libslic3r/PrintConfig.hpp"
#include "libslic3r/BoundingBox.hpp" #include "libslic3r/BoundingBox.hpp"
#include "libslic3r/AABBTreeIndirect.hpp" #include "libslic3r/AABBTreeIndirect.hpp"
#include "libslic3r/KDTreeIndirect.hpp"
namespace Slic3r { namespace Slic3r {
@ -16,119 +20,148 @@ class PrintObject;
class ExtrusionLoop; class ExtrusionLoop;
class Print; class Print;
class Layer; class Layer;
namespace EdgeGrid { class Grid; }
class SeamHistory {
public:
SeamHistory() { clear(); }
std::optional<Point> get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb);
void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb);
void clear();
private:
struct SeamPoint {
Point m_pos;
BoundingBox m_island_bb;
};
std::map<const PrintObject*, std::vector<SeamPoint>> m_data_last_layer;
std::map<const PrintObject*, std::vector<SeamPoint>> m_data_this_layer;
size_t m_layer_id;
};
class SeamPlacer {
public:
void init(const Print& print);
// When perimeters are printed, first call this function with the respective
// external perimeter. SeamPlacer will find a location for its seam and remember it.
// Subsequent calls to get_seam will return this position.
void plan_perimeters(const std::vector<const ExtrusionEntity*> perimeters,
const Layer& layer, SeamPosition seam_position,
Point last_pos, coordf_t nozzle_dmr, const PrintObject* po,
const EdgeGrid::Grid* lower_layer_edge_grid);
void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
const EdgeGrid::Grid* lower_layer_edge_grid);
using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
private:
// When given an external perimeter (!), returns the seam.
Point calculate_seam(const Layer& layer, const SeamPosition seam_position,
const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po,
const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos);
struct CustomTrianglesPerLayer {
Polygons polys;
TreeType tree;
};
// Just a cache to save some lookups.
const Layer* m_last_layer_po = nullptr;
coordf_t m_last_print_z = -1.;
const PrintObject* m_last_po = nullptr;
struct SeamPoint {
Point pt;
bool precalculated = false;
bool external = false;
const Layer* layer = nullptr;
SeamPosition seam_position;
const PrintObject* po = nullptr;
};
std::vector<SeamPoint> m_plan;
size_t m_plan_idx;
std::vector<std::vector<CustomTrianglesPerLayer>> m_enforcers;
std::vector<std::vector<CustomTrianglesPerLayer>> m_blockers;
std::vector<const PrintObject*> m_po_list;
//std::map<const PrintObject*, Point> m_last_seam_position;
SeamHistory m_seam_history;
// Get indices of points inside enforcers and blockers.
void get_enforcers_and_blockers(size_t layer_id,
const Polygon& polygon,
size_t po_id,
std::vector<size_t>& enforcers_idxs,
std::vector<size_t>& blockers_idxs) const;
// Apply penalties to points inside enforcers/blockers.
void apply_custom_seam(const Polygon& polygon, size_t po_id,
std::vector<float>& penalties,
const std::vector<float>& lengths,
int layer_id, SeamPosition seam_position) const;
// Return random point of a polygon. The distribution will be uniform
// along the contour and account for enforcers and blockers.
Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id,
bool* saw_custom = nullptr) const;
// Is there any enforcer/blocker on this layer?
bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const {
return is_custom_enforcer_on_layer(layer_id, po_idx)
|| is_custom_blocker_on_layer(layer_id, po_idx);
}
bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const {
return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty());
}
bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const {
return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty());
}
};
namespace EdgeGrid {
class Grid;
} }
namespace SeamPlacerImpl {
// ************ FOR BACKPORT COMPATIBILITY ONLY ***************
// Angle from v1 to v2, returning double atan2(y, x) normalized to <-PI, PI>.
template<typename Derived, typename Derived2> inline double angle(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2)
{
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "angle(): first parameter is not a 2D vector");
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "angle(): second parameter is not a 2D vector");
auto v1d = v1.template cast<double>();
auto v2d = v2.template cast<double>();
return atan2(cross2(v1d, v2d), v1d.dot(v2d));
}
// ***************************
struct GlobalModelInfo;
struct SeamComparator;
enum class EnforcedBlockedSeamPoint {
Blocked = 0,
Neutral = 1,
Enforced = 2,
};
// struct representing single perimeter loop
struct Perimeter
{
size_t start_index{};
size_t end_index{}; // inclusive!
size_t seam_index{};
float flow_width{};
// During alignment, a final position may be stored here. In that case, finalized is set to true.
// Note that final seam position is not limited to points of the perimeter loop. In theory it can be any position
// Random position also uses this flexibility to set final seam point position
bool finalized = false;
Vec3f final_seam_position = Vec3f::Zero();
};
// Struct over which all processing of perimeters is done. For each perimeter point, its respective candidate is created,
// then all the needed attributes are computed and finally, for each perimeter one point is chosen as seam.
// This seam position can be then further aligned
struct SeamCandidate
{
SeamCandidate(const Vec3f &pos, Perimeter &perimeter, float local_ccw_angle, EnforcedBlockedSeamPoint type)
: position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), embedded_distance(0.0f), local_ccw_angle(local_ccw_angle), type(type), central_enforcer(false)
{}
const Vec3f position;
// pointer to Perimeter loop of this point. It is shared across all points of the loop
Perimeter &perimeter;
float visibility;
float overhang;
// distance inside the merged layer regions, for detecting perimeter points which are hidden indside the print (e.g. multimaterial join)
// Negative sign means inside the print, comes from EdgeGrid structure
float embedded_distance;
float local_ccw_angle;
EnforcedBlockedSeamPoint type;
bool central_enforcer; // marks this candidate as central point of enforced segment on the perimeter - important for alignment
};
struct SeamCandidateCoordinateFunctor
{
SeamCandidateCoordinateFunctor(const std::vector<SeamCandidate> &seam_candidates) : seam_candidates(seam_candidates) {}
const std::vector<SeamCandidate> &seam_candidates;
float operator()(size_t index, size_t dim) const { return seam_candidates[index].position[dim]; }
};
} // namespace SeamPlacerImpl
struct PrintObjectSeamData
{
using SeamCandidatesTree = KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>;
struct LayerSeams
{
Slic3r::deque<SeamPlacerImpl::Perimeter> perimeters;
std::vector<SeamPlacerImpl::SeamCandidate> points;
std::unique_ptr<SeamCandidatesTree> points_tree;
};
// Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter
std::vector<LayerSeams> layers;
// Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD
// tree of all points of the given layer
void clear() { layers.clear(); }
};
class SeamPlacer
{
public:
// Number of samples generated on the mesh. There are sqr_rays_per_sample_point*sqr_rays_per_sample_point rays casted from each samples
static constexpr size_t raycasting_visibility_samples_count = 30000;
// square of number of rays per sample point
static constexpr size_t sqr_rays_per_sample_point = 5;
// arm length used during angles computation
static constexpr float polygon_local_angles_arm_distance = 0.3f;
// snapping angle - angles larger than this value will be snapped to during seam painting
static constexpr float sharp_angle_snapping_threshold = 55.0f * float(PI) / 180.0f;
// overhang angle for seam placement that still yields good results, in degrees, measured from vertical direction
static constexpr float overhang_angle_threshold = 45.0f * float(PI) / 180.0f;
// determines angle importance compared to visibility ( neutral value is 1.0f. )
static constexpr float angle_importance_aligned = 0.6f;
static constexpr float angle_importance_nearest = 1.0f; // use much higher angle importance for nearest mode, to combat the visibility info noise
// For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size
static constexpr float enforcer_oversampling_distance = 0.2f;
// When searching for seam clusters for alignment:
// following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer
static constexpr float seam_align_score_tolerance = 0.3f;
// seam_align_tolerable_dist_factor - how far to search for seam from current position, final dist is seam_align_tolerable_dist_factor * flow_width
static constexpr float seam_align_tolerable_dist_factor = 4.0f;
// minimum number of seams needed in cluster to make alignment happen
static constexpr size_t seam_align_minimum_string_seams = 6;
// millimeters covered by spline; determines number of splines for the given string
static constexpr size_t seam_align_mm_per_segment = 4.0f;
// The following data structures hold all perimeter points for all PrintObject.
std::unordered_map<const PrintObject *, PrintObjectSeamData> m_seam_per_object;
void init(const Print &print, std::function<void(void)> throw_if_canceled_func);
void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const;
private:
void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info, const SeamPosition configured_seam_preference);
void calculate_candidates_visibility(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info);
void calculate_overhangs_and_layer_embedding(const PrintObject *po);
void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator);
std::vector<std::pair<size_t, size_t>> find_seam_string(const PrintObject *po, std::pair<size_t, size_t> start_seam, const SeamPlacerImpl::SeamComparator &comparator) const;
std::optional<std::pair<size_t, size_t>> find_next_seam_in_layer(const std::vector<PrintObjectSeamData::LayerSeams> &layers,
const Vec3f & projected_position,
const size_t layer_idx,
const float max_distance,
const SeamPlacerImpl::SeamComparator & comparator) const;
};
} // namespace Slic3r
#endif // libslic3r_SeamPlacer_hpp_ #endif // libslic3r_SeamPlacer_hpp_

291
src/libslic3r/Geometry/Bicubic.hpp Normal file
View file

@ -0,0 +1,291 @@
#ifndef BICUBIC_HPP
#define BICUBIC_HPP
#include <algorithm>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
namespace Slic3r {
namespace Geometry {
namespace BicubicInternal {
// Linear kernel, to be able to test cubic methods with hat kernels.
template<typename T>
struct LinearKernel
{
typedef T FloatType;
static T a00() {
return T(0.);
}
static T a01() {
return T(0.);
}
static T a02() {
return T(0.);
}
static T a03() {
return T(0.);
}
static T a10() {
return T(1.);
}
static T a11() {
return T(-1.);
}
static T a12() {
return T(0.);
}
static T a13() {
return T(0.);
}
static T a20() {
return T(0.);
}
static T a21() {
return T(1.);
}
static T a22() {
return T(0.);
}
static T a23() {
return T(0.);
}
static T a30() {
return T(0.);
}
static T a31() {
return T(0.);
}
static T a32() {
return T(0.);
}
static T a33() {
return T(0.);
}
};
// Interpolation kernel aka Catmul-Rom aka Keyes kernel.
template<typename T>
struct CubicCatmulRomKernel
{
typedef T FloatType;
static T a00() {
return 0;
}
static T a01() {
return T( -0.5);
}
static T a02() {
return T( 1.);
}
static T a03() {
return T( -0.5);
}
static T a10() {
return T( 1.);
}
static T a11() {
return 0;
}
static T a12() {
return T( -5. / 2.);
}
static T a13() {
return T( 3. / 2.);
}
static T a20() {
return 0;
}
static T a21() {
return T( 0.5);
}
static T a22() {
return T( 2.);
}
static T a23() {
return T( -3. / 2.);
}
static T a30() {
return 0;
}
static T a31() {
return 0;
}
static T a32() {
return T( -0.5);
}
static T a33() {
return T( 0.5);
}
};
// B-spline kernel
template<typename T>
struct CubicBSplineKernel
{
typedef T FloatType;
static T a00() {
return T( 1. / 6.);
}
static T a01() {
return T( -3. / 6.);
}
static T a02() {
return T( 3. / 6.);
}
static T a03() {
return T( -1. / 6.);
}
static T a10() {
return T( 4. / 6.);
}
static T a11() {
return 0;
}
static T a12() {
return T( -6. / 6.);
}
static T a13() {
return T( 3. / 6.);
}
static T a20() {
return T( 1. / 6.);
}
static T a21() {
return T( 3. / 6.);
}
static T a22() {
return T( 3. / 6.);
}
static T a23() {
return T( -3. / 6.);
}
static T a30() {
return 0;
}
static T a31() {
return 0;
}
static T a32() {
return 0;
}
static T a33() {
return T( 1. / 6.);
}
};
template<class T>
inline T clamp(T a, T lower, T upper)
{
return (a < lower) ? lower :
(a > upper) ? upper : a;
}
}
template<typename Kernel>
struct CubicKernelWrapper
{
typedef typename Kernel::FloatType FloatType;
static constexpr size_t kernel_span = 4;
static FloatType kernel(FloatType x)
{
x = fabs(x);
if (x >= (FloatType) 2.)
return 0.0f;
if (x <= (FloatType) 1.) {
FloatType x2 = x * x;
FloatType x3 = x2 * x;
return Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3;
}
assert(x > (FloatType )1. && x < (FloatType )2.);
x -= (FloatType) 1.;
FloatType x2 = x * x;
FloatType x3 = x2 * x;
return Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3;
}
static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x)
{
const FloatType x2 = x * x;
const FloatType x3 = x * x * x;
return f0 * (Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3) +
f1 * (Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3) +
f2 * (Kernel::a20() + Kernel::a21() * x + Kernel::a22() * x2 + Kernel::a23() * x3) +
f3 * (Kernel::a30() + Kernel::a31() * x + Kernel::a32() * x2 + Kernel::a33() * x3);
}
};
// Linear splines
template<typename NumberType>
using LinearKernel = CubicKernelWrapper<BicubicInternal::LinearKernel<NumberType>>;
// Catmul-Rom splines
template<typename NumberType>
using CubicCatmulRomKernel = CubicKernelWrapper<BicubicInternal::CubicCatmulRomKernel<NumberType>>;
// Cubic B-splines
template<typename NumberType>
using CubicBSplineKernel = CubicKernelWrapper<BicubicInternal::CubicBSplineKernel<NumberType>>;
template<typename KernelWrapper>
static typename KernelWrapper::FloatType cubic_interpolate(const Eigen::ArrayBase<typename KernelWrapper::FloatType> &F,
const typename KernelWrapper::FloatType pt) {
typedef typename KernelWrapper::FloatType T;
const int w = int(F.size());
const int ix = (int) floor(pt);
const T s = pt - T( ix);
if (ix > 1 && ix + 2 < w) {
// Inside the fully interpolated region.
return KernelWrapper::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
}
// Transition region. Extend with a constant function.
auto f = [&F, w](T x) {
return F[BicubicInternal::clamp(x, 0, w - 1)];
};
return KernelWrapper::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s);
}
template<typename Kernel, typename Derived>
static float bicubic_interpolate(const Eigen::MatrixBase<Derived> &F,
const Eigen::Matrix<typename Kernel::FloatType, 2, 1, Eigen::DontAlign> &pt) {
typedef typename Kernel::FloatType T;
const int w = F.cols();
const int h = F.rows();
const int ix = (int) floor(pt[0]);
const int iy = (int) floor(pt[1]);
const T s = pt[0] - T( ix);
const T t = pt[1] - T( iy);
if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) {
// Inside the fully interpolated region.
return Kernel::interpolate(
Kernel::interpolate(F(ix - 1, iy - 1), F(ix, iy - 1), F(ix + 1, iy - 1), F(ix + 2, iy - 1), s),
Kernel::interpolate(F(ix - 1, iy), F(ix, iy), F(ix + 1, iy), F(ix + 2, iy), s),
Kernel::interpolate(F(ix - 1, iy + 1), F(ix, iy + 1), F(ix + 1, iy + 1), F(ix + 2, iy + 1), s),
Kernel::interpolate(F(ix - 1, iy + 2), F(ix, iy + 2), F(ix + 1, iy + 2), F(ix + 2, iy + 2), s), t);
}
// Transition region. Extend with a constant function.
auto f = [&F, w, h](int x, int y) {
return F(BicubicInternal::clamp(x, 0, w - 1), BicubicInternal::clamp(y, 0, h - 1));
};
return Kernel::interpolate(
Kernel::interpolate(f(ix - 1, iy - 1), f(ix, iy - 1), f(ix + 1, iy - 1), f(ix + 2, iy - 1), s),
Kernel::interpolate(f(ix - 1, iy), f(ix, iy), f(ix + 1, iy), f(ix + 2, iy), s),
Kernel::interpolate(f(ix - 1, iy + 1), f(ix, iy + 1), f(ix + 1, iy + 1), f(ix + 2, iy + 1), s),
Kernel::interpolate(f(ix - 1, iy + 2), f(ix, iy + 2), f(ix + 1, iy + 2), f(ix + 2, iy + 2), s), t);
}
} //namespace Geometry
} // namespace Slic3r
#endif /* BICUBIC_HPP */

218
src/libslic3r/Geometry/Curves.hpp Normal file
View file

@ -0,0 +1,218 @@
#ifndef SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
#define SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_
#include "libslic3r/Point.hpp"
#include "Bicubic.hpp"
#include <iostream>
//#define LSQR_DEBUG
namespace Slic3r {
namespace Geometry {
template<int Dimension, typename NumberType>
struct PolynomialCurve {
Eigen::MatrixXf coefficients;
Vec<Dimension, NumberType> get_fitted_value(const NumberType& value) const {
Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
size_t order = this->coefficients.rows() - 1;
auto x = NumberType(1.);
for (size_t index = 0; index < order + 1; ++index, x *= value)
result += x * this->coefficients.col(index);
return result;
}
};
//https://towardsdatascience.com/least-square-polynomial-CURVES-using-c-eigen-package-c0673728bd01
template<int Dimension, typename NumberType>
PolynomialCurve<Dimension, NumberType> fit_polynomial(const std::vector<Vec<Dimension, NumberType>> &observations,
const std::vector<NumberType> &observation_points,
const std::vector<NumberType> &weights, size_t order) {
// check to make sure inputs are correct
size_t cols = order + 1;
assert(observation_points.size() >= cols);
assert(observation_points.size() == weights.size());
assert(observations.size() == weights.size());
Eigen::MatrixXf data_points(Dimension, observations.size());
Eigen::MatrixXf T(observations.size(), cols);
for (size_t i = 0; i < weights.size(); ++i) {
auto squared_weight = sqrt(weights[i]);
data_points.col(i) = observations[i] * squared_weight;
// Populate the matrix
auto x = squared_weight;
auto c = observation_points[i];
for (size_t j = 0; j < cols; ++j, x *= c)
T(i, j) = x;
}
const auto QR = T.householderQr();
Eigen::MatrixXf coefficients(Dimension, cols);
// Solve for linear least square fit
for (size_t dim = 0; dim < Dimension; ++dim) {
coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
}
return {std::move(coefficients)};
}
template<size_t Dimension, typename NumberType, typename KernelType>
struct PiecewiseFittedCurve {
using Kernel = KernelType;
Eigen::MatrixXf coefficients;
NumberType start;
NumberType segment_size;
size_t endpoints_level_of_freedom;
Vec<Dimension, NumberType> get_fitted_value(const NumberType &observation_point) const {
Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
//find corresponding segment index; expects kernels to be centered
int middle_right_segment_index = floor((observation_point - start) / segment_size);
//find index of first segment that is affected by the point i; this can be deduced from kernel_span
int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1;
for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
segment_index++) {
NumberType segment_start = start + segment_index * segment_size;
NumberType normalized_segment_distance = (segment_start - observation_point) / segment_size;
int parameter_index = segment_index + endpoints_level_of_freedom;
parameter_index = std::clamp(parameter_index, 0, int(coefficients.cols()) - 1);
result += Kernel::kernel(normalized_segment_distance) * coefficients.col(parameter_index);
}
return result;
}
};
// observations: data to be fitted by the curve
// observation points: growing sequence of points where the observations were made.
// In other words, for function f(x) = y, observations are y0...yn, and observation points are x0...xn
// weights: how important the observation is
// segments_count: number of segments inside the valid length of the curve
// endpoints_level_of_freedom: number of additional parameters at each end; reasonable values depend on the kernel span
template<typename Kernel, int Dimension, typename NumberType>
PiecewiseFittedCurve<Dimension, NumberType, Kernel> fit_curve(
const std::vector<Vec<Dimension, NumberType>> &observations,
const std::vector<NumberType> &observation_points,
const std::vector<NumberType> &weights,
size_t segments_count,
size_t endpoints_level_of_freedom) {
// check to make sure inputs are correct
assert(segments_count > 0);
assert(observations.size() > 0);
assert(observation_points.size() == observations.size());
assert(observation_points.size() == weights.size());
assert(segments_count <= observations.size());
//prepare sqrt of weights, which will then be applied to both matrix T and observed data: https://en.wikipedia.org/wiki/Weighted_least_squares
std::vector<NumberType> sqrt_weights(weights.size());
for (size_t index = 0; index < weights.size(); ++index) {
assert(weights[index] > 0);
sqrt_weights[index] = sqrt(weights[index]);
}
// prepare result and compute metadata
PiecewiseFittedCurve<Dimension, NumberType, Kernel> result { };
NumberType valid_length = observation_points.back() - observation_points.front();
NumberType segment_size = valid_length / NumberType(segments_count);
result.start = observation_points.front();
result.segment_size = segment_size;
result.endpoints_level_of_freedom = endpoints_level_of_freedom;
// prepare observed data
// Eigen defaults to column major memory layout.
Eigen::MatrixXf data_points(Dimension, observations.size());
for (size_t index = 0; index < observations.size(); ++index) {
data_points.col(index) = observations[index] * sqrt_weights[index];
}
// parameters count is always increased by one to make the parametric space of the curve symmetric.
// without this fix, the end of the curve is less flexible than the beginning
size_t parameters_count = segments_count + 1 + 2 * endpoints_level_of_freedom;
//Create weight matrix T for each point and each segment;
Eigen::MatrixXf T(observation_points.size(), parameters_count);
T.setZero();
//Fill the weight matrix
for (size_t i = 0; i < observation_points.size(); ++i) {
NumberType observation_point = observation_points[i];
//find corresponding segment index; expects kernels to be centered
int middle_right_segment_index = floor((observation_point - result.start) / result.segment_size);
//find index of first segment that is affected by the point i; this can be deduced from kernel_span
int start_segment_idx = middle_right_segment_index - int(Kernel::kernel_span / 2) + 1;
for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span);
segment_index++) {
NumberType segment_start = result.start + segment_index * result.segment_size;
NumberType normalized_segment_distance = (segment_start - observation_point) / result.segment_size;
int parameter_index = segment_index + endpoints_level_of_freedom;
parameter_index = std::clamp(parameter_index, 0, int(parameters_count) - 1);
T(i, parameter_index) += Kernel::kernel(normalized_segment_distance) * sqrt_weights[i];
}
}
#ifdef LSQR_DEBUG
std::cout << "weight matrix: " << std::endl;
for (int obs = 0; obs < observation_points.size(); ++obs) {
std::cout << std::endl;
for (int segment = 0; segment < parameters_count; ++segment) {
std::cout << T(obs, segment) << " ";
}
}
std::cout << std::endl;
#endif
// Solve for linear least square fit
result.coefficients.resize(Dimension, parameters_count);
const auto QR = T.fullPivHouseholderQr();
for (size_t dim = 0; dim < Dimension; ++dim) {
result.coefficients.row(dim) = QR.solve(data_points.row(dim).transpose());
}
return result;
}
template<int Dimension, typename NumberType>
PiecewiseFittedCurve<Dimension, NumberType, LinearKernel<NumberType>>
fit_linear_spline(
const std::vector<Vec<Dimension, NumberType>> &observations,
std::vector<NumberType> observation_points,
std::vector<NumberType> weights,
size_t segments_count,
size_t endpoints_level_of_freedom = 0) {
return fit_curve<LinearKernel<NumberType>>(observations, observation_points, weights, segments_count,
endpoints_level_of_freedom);
}
template<int Dimension, typename NumberType>
PiecewiseFittedCurve<Dimension, NumberType, CubicBSplineKernel<NumberType>>
fit_cubic_bspline(
const std::vector<Vec<Dimension, NumberType>> &observations,
std::vector<NumberType> observation_points,
std::vector<NumberType> weights,
size_t segments_count,
size_t endpoints_level_of_freedom = 0) {
return fit_curve<CubicBSplineKernel<NumberType>>(observations, observation_points, weights, segments_count,
endpoints_level_of_freedom);
}
template<int Dimension, typename NumberType>
PiecewiseFittedCurve<Dimension, NumberType, CubicCatmulRomKernel<NumberType>>
fit_catmul_rom_spline(
const std::vector<Vec<Dimension, NumberType>> &observations,
std::vector<NumberType> observation_points,
std::vector<NumberType> weights,
size_t segments_count,
size_t endpoints_level_of_freedom = 0) {
return fit_curve<CubicCatmulRomKernel<NumberType>>(observations, observation_points, weights, segments_count,
endpoints_level_of_freedom);
}
}
}
#endif /* SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ */

498
src/libslic3r/KDTreeIndirect.hpp
View file

@ -11,224 +11,362 @@
namespace Slic3r { namespace Slic3r {
enum class VisitorReturnMask : unsigned int {
CONTINUE_LEFT = 1,
CONTINUE_RIGHT = 2,
STOP = 4,
};
// KD tree for N-dimensional closest point search. // KD tree for N-dimensional closest point search.
template<size_t ANumDimensions, typename ACoordType, typename ACoordinateFn> template<size_t ANumDimensions, typename ACoordType, typename ACoordinateFn>
class KDTreeIndirect class KDTreeIndirect
{ {
public: public:
static constexpr size_t NumDimensions = ANumDimensions; static constexpr size_t NumDimensions = ANumDimensions;
using CoordinateFn = ACoordinateFn; using CoordinateFn = ACoordinateFn;
using CoordType = ACoordType; using CoordType = ACoordType;
// Following could be static constexpr size_t, but that would not link in C++11 // Following could be static constexpr size_t, but that would not link in C++11
enum : size_t { enum : size_t {
npos = size_t(-1) npos = size_t(-1)
}; };
KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {} KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {}
KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> indices) : coordinate(coordinate) { this->build(std::move(indices)); } KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> indices) : coordinate(coordinate) { this->build(indices); }
KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> &&indices) : coordinate(coordinate) { this->build(std::move(indices)); } KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); }
KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); } KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {}
KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {} KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; }
KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; } void clear() { m_nodes.clear(); }
void clear() { m_nodes.clear(); }
void build(size_t num_indices) void build(size_t num_indices)
{ {
std::vector<size_t> indices; std::vector<size_t> indices;
indices.reserve(num_indices); indices.reserve(num_indices);
for (size_t i = 0; i < num_indices; ++ i) for (size_t i = 0; i < num_indices; ++ i)
indices.emplace_back(i); indices.emplace_back(i);
this->build(std::move(indices)); this->build(indices);
} }
void build(std::vector<size_t> &&indices) void build(std::vector<size_t> &indices)
{ {
if (indices.empty()) if (indices.empty())
clear(); clear();
else { else {
// Allocate enough memory for a full binary tree. // Allocate enough memory for a full binary tree.
m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos); m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos);
build_recursive(indices, 0, 0, 0, indices.size() - 1); build_recursive(indices, 0, 0, 0, indices.size() - 1);
} }
indices.clear(); indices.clear();
} }
enum class VisitorReturnMask : unsigned int template<typename CoordType>
{ unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
CONTINUE_LEFT = 1, {
CONTINUE_RIGHT = 2, CoordType dist = point_coord - this->coordinate(idx, dimension);
STOP = 4, return (dist * dist < search_radius + CoordType(EPSILON)) ?
}; // The plane intersects a hypersphere centered at point_coord of search_radius.
template<typename CoordType> ((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) :
unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const // The plane does not intersect the hypersphere.
{ (dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT);
CoordType dist = point_coord - this->coordinate(idx, dimension); }
return (dist * dist < search_radius + CoordType(EPSILON)) ?
// The plane intersects a hypersphere centered at point_coord of search_radius.
((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) :
// The plane does not intersect the hypersphere.
(dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT);
}
// Visitor is supposed to return a bit mask of VisitorReturnMask. // Visitor is supposed to return a bit mask of VisitorReturnMask.
template<typename Visitor> template<typename Visitor>
void visit(Visitor &visitor) const void visit(Visitor &visitor) const
{ {
visit_recursive(0, 0, visitor); visit_recursive(0, 0, visitor);
} }
CoordinateFn coordinate; CoordinateFn coordinate;
private: private:
// Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension. // Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension.
void build_recursive(std::vector<size_t> &input, size_t node, const size_t dimension, const size_t left, const size_t right) void build_recursive(std::vector<size_t> &input, size_t node, const size_t dimension, const size_t left, const size_t right)
{ {
if (left > right) if (left > right)
return; return;
assert(node < m_nodes.size()); assert(node < m_nodes.size());
if (left == right) { if (left == right) {
// Insert a node into the balanced tree. // Insert a node into the balanced tree.
m_nodes[node] = input[left]; m_nodes[node] = input[left];
return; return;
} }
// Partition the input to left / right pieces of the same length to produce a balanced tree. // Partition the input to left / right pieces of the same length to produce a balanced tree.
size_t center = (left + right) / 2; size_t center = (left + right) / 2;
partition_input(input, dimension, left, right, center); partition_input(input, dimension, left, right, center);
// Insert a node into the tree. // Insert a node into the tree.
m_nodes[node] = input[center]; m_nodes[node] = input[center];
// Build up the left / right subtrees. // Build up the left / right subtrees.
size_t next_dimension = dimension; size_t next_dimension = dimension;
if (++ next_dimension == NumDimensions) if (++ next_dimension == NumDimensions)
next_dimension = 0; next_dimension = 0;
if (center > left) if (center > left)
build_recursive(input, node * 2 + 1, next_dimension, left, center - 1); build_recursive(input, node * 2 + 1, next_dimension, left, center - 1);
build_recursive(input, node * 2 + 2, next_dimension, center + 1, right); build_recursive(input, node * 2 + 2, next_dimension, center + 1, right);
} }
// Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method: // Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
// https://en.wikipedia.org/wiki/Quickselect // https://en.wikipedia.org/wiki/Quickselect
// Items left of the k'th item are lower than the k'th item in the "dimension", // Items left of the k'th item are lower than the k'th item in the "dimension",
// items right of the k'th item are higher than the k'th item in the "dimension", // items right of the k'th item are higher than the k'th item in the "dimension",
void partition_input(std::vector<size_t> &input, const size_t dimension, size_t left, size_t right, const size_t k) const void partition_input(std::vector<size_t> &input, const size_t dimension, size_t left, size_t right, const size_t k) const
{ {
while (left < right) { while (left < right) {
size_t center = (left + right) / 2; size_t center = (left + right) / 2;
CoordType pivot; CoordType pivot;
{ {
// Bubble sort the input[left], input[center], input[right], so that a median of the three values // Bubble sort the input[left], input[center], input[right], so that a median of the three values
// will end up in input[center]. // will end up in input[center].
CoordType left_value = this->coordinate(input[left], dimension); CoordType left_value = this->coordinate(input[left], dimension);
CoordType center_value = this->coordinate(input[center], dimension); CoordType center_value = this->coordinate(input[center], dimension);
CoordType right_value = this->coordinate(input[right], dimension); CoordType right_value = this->coordinate(input[right], dimension);
if (left_value > center_value) { if (left_value > center_value) {
std::swap(input[left], input[center]); std::swap(input[left], input[center]);
std::swap(left_value, center_value); std::swap(left_value, center_value);
} }
if (left_value > right_value) { if (left_value > right_value) {
std::swap(input[left], input[right]); std::swap(input[left], input[right]);
right_value = left_value; right_value = left_value;
} }
if (center_value > right_value) { if (center_value > right_value) {
std::swap(input[center], input[right]); std::swap(input[center], input[right]);
center_value = right_value; center_value = right_value;
} }
pivot = center_value; pivot = center_value;
} }
if (right <= left + 2) if (right <= left + 2)
// The <left, right> interval is already sorted. // The <left, right> interval is already sorted.
break; break;
size_t i = left; size_t i = left;
size_t j = right - 1; size_t j = right - 1;
std::swap(input[center], input[j]); std::swap(input[center], input[j]);
// Partition the set based on the pivot. // Partition the set based on the pivot.
for (;;) { for (;;) {
// Skip left points that are already at correct positions. // Skip left points that are already at correct positions.
// Search will certainly stop at position (right - 1), which stores the pivot. // Search will certainly stop at position (right - 1), which stores the pivot.
while (this->coordinate(input[++ i], dimension) < pivot) ; while (this->coordinate(input[++ i], dimension) < pivot) ;
// Skip right points that are already at correct positions. // Skip right points that are already at correct positions.
while (this->coordinate(input[-- j], dimension) > pivot && i < j) ; while (this->coordinate(input[-- j], dimension) > pivot && i < j) ;
if (i >= j) if (i >= j)
break; break;
std::swap(input[i], input[j]); std::swap(input[i], input[j]);
} }
// Restore pivot to the center of the sequence. // Restore pivot to the center of the sequence.
std::swap(input[i], input[right - 1]); std::swap(input[i], input[right - 1]);
// Which side the kth element is in? // Which side the kth element is in?
if (k < i) if (k < i)
right = i - 1; right = i - 1;
else if (k == i) else if (k == i)
// Sequence is partitioned, kth element is at its place. // Sequence is partitioned, kth element is at its place.
break; break;
else else
left = i + 1; left = i + 1;
} }
} }
template<typename Visitor> template<typename Visitor>
void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const
{ {
assert(! m_nodes.empty()); assert(! m_nodes.empty());
if (node >= m_nodes.size() || m_nodes[node] == npos) if (node >= m_nodes.size() || m_nodes[node] == npos)
return; return;
// Left / right child node index. // Left / right child node index.
size_t left = node * 2 + 1; size_t left = node * 2 + 1;
size_t right = left + 1; size_t right = left + 1;
unsigned int mask = visitor(m_nodes[node], dimension); unsigned int mask = visitor(m_nodes[node], dimension);
if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) { if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) {
size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension; size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension;
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT) if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT)
visit_recursive(left, next_dimension, visitor); visit_recursive(left, next_dimension, visitor);
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT) if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT)
visit_recursive(right, next_dimension, visitor); visit_recursive(right, next_dimension, visitor);
} }
} }
std::vector<size_t> m_nodes; std::vector<size_t> m_nodes;
}; };
// Find a closest point using Euclidian metrics. // Find a closest point using Euclidian metrics.
// Returns npos if not found. // Returns npos if not found.
template<typename KDTreeIndirectType, typename PointType, typename FilterFn> template<size_t K,
size_t find_closest_point(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) typename PointType,
typename FilterFn,
size_t D,
typename CoordT,
typename CoordFn>
std::array<size_t, K> find_closest_points(
const KDTreeIndirect<D, CoordT, CoordFn> &kdtree,
const PointType &point,
FilterFn filter)
{ {
using CoordType = typename KDTreeIndirectType::CoordType; using Tree = KDTreeIndirect<D, CoordT, CoordFn>;
struct Visitor { struct Visitor
const KDTreeIndirectType &kdtree; {
const PointType &point; const Tree &kdtree;
const FilterFn filter; const PointType &point;
size_t min_idx = KDTreeIndirectType::npos; const FilterFn filter;
CoordType min_dist = std::numeric_limits<CoordType>::max();
Visitor(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) : kdtree(kdtree), point(point), filter(filter) {} std::array<std::pair<size_t, CoordT>, K> results;
unsigned int operator()(size_t idx, size_t dimension) {
if (this->filter(idx)) {
auto dist = CoordType(0);
for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++ i) {
CoordType d = point[i] - kdtree.coordinate(idx, i);
dist += d * d;
}
if (dist < min_dist) {
min_dist = dist;
min_idx = idx;
}
}
return kdtree.descent_mask(point[dimension], min_dist, idx, dimension);
}
} visitor(kdtree, point, filter);
kdtree.visit(visitor); Visitor(const Tree &kdtree, const PointType &point, FilterFn filter)
return visitor.min_idx; : kdtree(kdtree), point(point), filter(filter)
{
results.fill(std::make_pair(Tree::npos,
std::numeric_limits<CoordT>::max()));
}
unsigned int operator()(size_t idx, size_t dimension)
{
if (this->filter(idx)) {
auto dist = CoordT(0);
for (size_t i = 0; i < D; ++i) {
CoordT d = point[i] - kdtree.coordinate(idx, i);
dist += d * d;
}
auto res = std::make_pair(idx, dist);
auto it = std::lower_bound(results.begin(), results.end(),
res, [](auto &r1, auto &r2) {
return r1.second < r2.second;
});
if (it != results.end()) {
std::rotate(it, std::prev(results.end()), results.end());
*it = res;
}
}
return kdtree.descent_mask(point[dimension],
results.front().second, idx,
dimension);
}
} visitor(kdtree, point, filter);
kdtree.visit(visitor);
std::array<size_t, K> ret;
for (size_t i = 0; i < K; i++) ret[i] = visitor.results[i].first;
return ret;
}
template<size_t K, typename PointType, size_t D, typename CoordT, typename CoordFn>
std::array<size_t, K> find_closest_points(
const KDTreeIndirect<D, CoordT, CoordFn> &kdtree, const PointType &point)
{
return find_closest_points<K>(kdtree, point, [](size_t) { return true; });
}
template<typename PointType,
typename FilterFn,
size_t D,
typename CoordT,
typename CoordFn>
size_t find_closest_point(const KDTreeIndirect<D, CoordT, CoordFn> &kdtree,
const PointType &point,
FilterFn filter)
{
return find_closest_points<1>(kdtree, point, filter)[0];
} }
template<typename KDTreeIndirectType, typename PointType> template<typename KDTreeIndirectType, typename PointType>
size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& point) size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& point)
{ {
return find_closest_point(kdtree, point, [](size_t) { return true; }); return find_closest_point(kdtree, point, [](size_t) { return true; });
}
// Find nearby points (spherical neighbourhood) using Euclidian metrics.
template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
const typename KDTreeIndirectType::CoordType& max_distance, FilterFn filter)
{
using CoordType = typename KDTreeIndirectType::CoordType;
struct Visitor {
const KDTreeIndirectType &kdtree;
const PointType center;
const CoordType max_distance_squared;
const FilterFn filter;
std::vector<size_t> result;
Visitor(const KDTreeIndirectType &kdtree, const PointType& center, const CoordType &max_distance,
FilterFn filter) :
kdtree(kdtree), center(center), max_distance_squared(max_distance*max_distance), filter(filter) {
}
unsigned int operator()(size_t idx, size_t dimension) {
if (this->filter(idx)) {
auto dist = CoordType(0);
for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
CoordType d = center[i] - kdtree.coordinate(idx, i);
dist += d * d;
}
if (dist < max_distance_squared) {
result.push_back(idx);
}
}
return kdtree.descent_mask(center[dimension], max_distance_squared, idx, dimension);
}
} visitor(kdtree, center, max_distance, filter);
kdtree.visit(visitor);
return visitor.result;
}
template<typename KDTreeIndirectType, typename PointType>
std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
const typename KDTreeIndirectType::CoordType& max_distance)
{
return find_nearby_points(kdtree, center, max_distance, [](size_t) {
return true;
});
}
// Find nearby points (spherical neighbourhood) using Euclidian metrics.
template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree,
const PointType &bb_min,
const PointType &bb_max,
FilterFn filter)
{
struct Visitor {
const KDTreeIndirectType &kdtree;
const PointType &bb_min, &bb_max;
const FilterFn filter;
std::vector<size_t> result;
Visitor(const KDTreeIndirectType &kdtree, const PointType& bbmin, const PointType& bbmax,
FilterFn filter) :
kdtree(kdtree), bb_min{bbmin}, bb_max{bbmax}, filter(filter) {
}
unsigned int operator()(size_t idx, size_t dimension) {
unsigned int ret =
static_cast<unsigned int>(VisitorReturnMask::CONTINUE_LEFT) |
static_cast<unsigned int>(VisitorReturnMask::CONTINUE_RIGHT);
if (this->filter(idx)) {
PointType p;
bool contains = true;
for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
p(i) = kdtree.coordinate(idx, i);
contains = contains && bb_min(i) <= p(i) && p(i) <= bb_max(i);
}
if (p(dimension) < bb_min(dimension))
ret = static_cast<unsigned int>(VisitorReturnMask::CONTINUE_RIGHT);
if (p(dimension) > bb_max(dimension))
ret = static_cast<unsigned int>(VisitorReturnMask::CONTINUE_LEFT);
if (contains)
result.emplace_back(idx);
}
return ret;
}
} visitor(kdtree, bb_min, bb_max, filter);
kdtree.visit(visitor);
return visitor.result;
} }
} // namespace Slic3r } // namespace Slic3r

17
src/libslic3r/MultiPoint.cpp
View file

@ -63,6 +63,23 @@ int MultiPoint::find_point(const Point &point) const
return -1; // not found return -1; // not found
} }
int MultiPoint::find_point(const Point &point, double scaled_epsilon) const
{
if (scaled_epsilon == 0) return this->find_point(point);
auto dist2_min = std::numeric_limits<double>::max();
auto eps2 = scaled_epsilon * scaled_epsilon;
int idx_min = -1;
for (const Point &pt : this->points) {
double d2 = (pt - point).cast<double>().squaredNorm();
if (d2 < dist2_min) {
idx_min = int(&pt - &this->points.front());
dist2_min = d2;
}
}
return dist2_min < eps2 ? idx_min : -1;
}
bool MultiPoint::has_boundary_point(const Point &point) const bool MultiPoint::has_boundary_point(const Point &point) const
{ {
double dist = (point.projection_onto(*this) - point).cast<double>().norm(); double dist = (point.projection_onto(*this) - point).cast<double>().norm();

5
src/libslic3r/MultiPoint.hpp
View file

@ -43,7 +43,12 @@ public:
double length() const; double length() const;
bool is_valid() const { return this->points.size() >= 2; } bool is_valid() const { return this->points.size() >= 2; }
// Return index of a polygon point exactly equal to point.
// Return -1 if no such point exists.
int find_point(const Point &point) const; int find_point(const Point &point) const;
// Return index of the closest point to point closer than scaled_epsilon.
// Return -1 if no such point exists.
int find_point(const Point &point, const double scaled_epsilon) const;
bool has_boundary_point(const Point &point) const; bool has_boundary_point(const Point &point) const;
int closest_point_index(const Point &point) const { int closest_point_index(const Point &point) const {
int idx = -1; int idx = -1;

142
src/libslic3r/NormalUtils.cpp Normal file
View file

@ -0,0 +1,142 @@
#include "NormalUtils.hpp"
using namespace Slic3r;
Vec3f NormalUtils::create_triangle_normal(
const stl_triangle_vertex_indices &indices,
const std::vector<stl_vertex> & vertices)
{
const stl_vertex &v0 = vertices[indices[0]];
const stl_vertex &v1 = vertices[indices[1]];
const stl_vertex &v2 = vertices[indices[2]];
Vec3f direction = (v1 - v0).cross(v2 - v0);
direction.normalize();
return direction;
}
std::vector<Vec3f> NormalUtils::create_triangle_normals(
const indexed_triangle_set &its)
{
std::vector<Vec3f> normals;
normals.reserve(its.indices.size());
for (const Vec3crd &index : its.indices) {
normals.push_back(create_triangle_normal(index, its.vertices));
}
return normals;
}
NormalUtils::Normals NormalUtils::create_normals_average_neighbor(
const indexed_triangle_set &its)
{
size_t count_vertices = its.vertices.size();
std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
std::vector<unsigned int> count(count_vertices, 0);
for (const Vec3crd &indice : its.indices) {
Vec3f normal = create_triangle_normal(indice, its.vertices);
for (int i = 0; i < 3; ++i) {
normals[indice[i]] += normal;
++count[indice[i]];
}
}
// normalize to size 1
for (auto &normal : normals) {
size_t index = &normal - &normals.front();
normal /= static_cast<float>(count[index]);
}
return normals;
}
// calc triangle angle of vertex defined by index to triangle indices
float NormalUtils::indice_angle(int i,
const Vec3crd & indice,
const std::vector<stl_vertex> &vertices)
{
int i1 = (i == 0) ? 2 : (i - 1);
int i2 = (i == 2) ? 0 : (i + 1);
Vec3f v1 = vertices[i1] - vertices[i];
Vec3f v2 = vertices[i2] - vertices[i];
v1.normalize();
v2.normalize();
float w = v1.dot(v2);
if (w > 1.f)
w = 1.f;
else if (w < -1.f)
w = -1.f;
return acos(w);
}
NormalUtils::Normals NormalUtils::create_normals_angle_weighted(
const indexed_triangle_set &its)
{
size_t count_vertices = its.vertices.size();
std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
std::vector<float> count(count_vertices, 0.f);
for (const Vec3crd &indice : its.indices) {
Vec3f normal = create_triangle_normal(indice, its.vertices);
Vec3f angles(indice_angle(0, indice, its.vertices),
indice_angle(1, indice, its.vertices), 0.f);
angles[2] = (M_PI - angles[0] - angles[1]);
for (int i = 0; i < 3; ++i) {
const float &weight = angles[i];
normals[indice[i]] += normal * weight;
count[indice[i]] += weight;
}
}
// normalize to size 1
for (auto &normal : normals) {
size_t index = &normal - &normals.front();
normal /= count[index];
}
return normals;
}
NormalUtils::Normals NormalUtils::create_normals_nelson_weighted(
const indexed_triangle_set &its)
{
size_t count_vertices = its.vertices.size();
std::vector<Vec3f> normals(count_vertices, Vec3f(.0f, .0f, .0f));
std::vector<float> count(count_vertices, 0.f);
const std::vector<stl_vertex> &vertices = its.vertices;
for (const Vec3crd &indice : its.indices) {
Vec3f normal = create_triangle_normal(indice, vertices);
const stl_vertex &v0 = vertices[indice[0]];
const stl_vertex &v1 = vertices[indice[1]];
const stl_vertex &v2 = vertices[indice[2]];
float e0 = (v0 - v1).norm();
float e1 = (v1 - v2).norm();
float e2 = (v2 - v0).norm();
Vec3f coefs(e0 * e2, e0 * e1, e1 * e2);
for (int i = 0; i < 3; ++i) {
const float &weight = coefs[i];
normals[indice[i]] += normal * weight;
count[indice[i]] += weight;
}
}
// normalize to size 1
for (auto &normal : normals) {
size_t index = &normal - &normals.front();
normal /= count[index];
}
return normals;
}
// calculate normals by averaging normals of neghbor triangles
std::vector<Vec3f> NormalUtils::create_normals(
const indexed_triangle_set &its, VertexNormalType type)
{
switch (type) {
case VertexNormalType::AverageNeighbor:
return create_normals_average_neighbor(its);
case VertexNormalType::AngleWeighted:
return create_normals_angle_weighted(its);
case VertexNormalType::NelsonMaxWeighted:
default:
return create_normals_nelson_weighted(its);
}
}

69
src/libslic3r/NormalUtils.hpp Normal file
View file

@ -0,0 +1,69 @@
#ifndef slic3r_NormalUtils_hpp_
#define slic3r_NormalUtils_hpp_
#include "Point.hpp"
#include "Model.hpp"
namespace Slic3r {
/// <summary>
/// Collection of static function
/// to create normals
/// </summary>
class NormalUtils
{
public:
using Normal = Vec3f;
using Normals = std::vector<Normal>;
NormalUtils() = delete; // only static functions
enum class VertexNormalType {
AverageNeighbor,
AngleWeighted,
NelsonMaxWeighted
};
/// <summary>
/// Create normal for triangle defined by indices from vertices
/// </summary>
/// <param name="indices">index into vertices</param>
/// <param name="vertices">vector of vertices</param>
/// <returns>normal to triangle(normalized to size 1)</returns>
static Normal create_triangle_normal(
const stl_triangle_vertex_indices &indices,
const std::vector<stl_vertex> & vertices);
/// <summary>
/// Create normals for each vertices
/// </summary>
/// <param name="its">indices and vertices</param>
/// <returns>Vector of normals</returns>
static Normals create_triangle_normals(const indexed_triangle_set &its);
/// <summary>
/// Create normals for each vertex by averaging neighbor triangles normal
/// </summary>
/// <param name="its">Triangle indices and vertices</param>
/// <param name="type">Type of calculation normals</param>
/// <returns>Normal for each vertex</returns>
static Normals create_normals(
const indexed_triangle_set &its,
VertexNormalType type = VertexNormalType::NelsonMaxWeighted);
static Normals create_normals_average_neighbor(const indexed_triangle_set &its);
static Normals create_normals_angle_weighted(const indexed_triangle_set &its);
static Normals create_normals_nelson_weighted(const indexed_triangle_set &its);
/// <summary>
/// Calculate angle of trinagle side.
/// </summary>
/// <param name="i">index to indices, define angle point</param>
/// <param name="indice">address to vertices</param>
/// <param name="vertices">vertices data</param>
/// <returns>Angle [in radian]</returns>
static float indice_angle(int i,
const Vec3crd & indice,
const std::vector<stl_vertex> &vertices);
};
} // namespace Slic3r
#endif // slic3r_NormalUtils_hpp_

22
src/libslic3r/Polyline.cpp
View file

@ -517,6 +517,28 @@ bool remove_degenerate(Polylines &polylines)
return modified; return modified;
} }
std::pair<int, Point> foot_pt(const Points &polyline, const Point &pt)
{
if (polyline.size() < 2) return std::make_pair(-1, Point(0, 0));
auto d2_min = std::numeric_limits<double>::max();
Point foot_pt_min;
Point prev = polyline.front();
auto it = polyline.begin();
auto it_proj = polyline.begin();
for (++it; it != polyline.end(); ++it) {
Point foot_pt = pt.projection_onto(Line(prev, *it));
double d2 = (foot_pt - pt).cast<double>().squaredNorm();
if (d2 < d2_min) {
d2_min = d2;
foot_pt_min = foot_pt;
it_proj = it;
}
prev = *it;
}
return std::make_pair(int(it_proj - polyline.begin()) - 1, foot_pt_min);
}
ThickLines ThickPolyline::thicklines() const ThickLines ThickPolyline::thicklines() const
{ {
ThickLines lines; ThickLines lines;

3
src/libslic3r/Polyline.hpp
View file

@ -222,6 +222,9 @@ const Point& leftmost_point(const Polylines &polylines);
bool remove_degenerate(Polylines &polylines); bool remove_degenerate(Polylines &polylines);
// Returns index of a segment of a polyline and foot point of pt on polyline.
std::pair<int, Point> foot_pt(const Points &polyline, const Point &pt);
class ThickPolyline : public Polyline { class ThickPolyline : public Polyline {
public: public:
ThickPolyline() : endpoints(std::make_pair(false, false)) {} ThickPolyline() : endpoints(std::make_pair(false, false)) {}

183
src/libslic3r/ShortEdgeCollapse.cpp Normal file
View file

@ -0,0 +1,183 @@
#include "ShortEdgeCollapse.hpp"
#include "libslic3r/NormalUtils.hpp"
#include <unordered_map>
#include <unordered_set>
#include <random>
#include <algorithm>
namespace Slic3r {
void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) {
// whenever vertex is removed, its mapping is update to the index of vertex with wich it merged
std::vector<size_t> vertices_index_mapping(mesh.vertices.size());
for (size_t idx = 0; idx < vertices_index_mapping.size(); ++idx) {
vertices_index_mapping[idx] = idx;
}
// Algorithm uses get_final_index query to get the actual vertex index. The query also updates all mappings on the way, essentially flattening the mapping
std::vector<size_t> flatten_queue;
auto get_final_index = [&vertices_index_mapping, &flatten_queue](const size_t &orig_index) {
flatten_queue.clear();
size_t idx = orig_index;
while (vertices_index_mapping[idx] != idx) {
flatten_queue.push_back(idx);
idx = vertices_index_mapping[idx];
}
for (size_t i : flatten_queue) {
vertices_index_mapping[i] = idx;
}
return idx;
};
// if face is removed, mark it here
std::vector<bool> face_removal_flags(mesh.indices.size(), false);
std::vector<Vec3i> triangles_neighbors = its_face_neighbors_par(mesh);
// now compute vertices dot product - this is used during edge collapse,
// to determine which vertex to remove and which to keep; We try to keep the one with larger angle, because it defines the shape "more".
// The min vertex dot product is lowest dot product of its normal with the normals of faces around it.
// the lower the dot product, the more we want to keep the vertex
// NOTE: This score is not updated, even though the decimation does change the mesh. It saves computation time, and there are no strong reasons to update.
std::vector<float> min_vertex_dot_product(mesh.vertices.size(), 1);
{
std::vector<Vec3f> face_normals = its_face_normals(mesh);
std::vector<Vec3f> vertex_normals = NormalUtils::create_normals(mesh);
for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) {
Vec3i t = mesh.indices[face_idx];
Vec3f n = face_normals[face_idx];
min_vertex_dot_product[t[0]] = std::min(min_vertex_dot_product[t[0]], n.dot(vertex_normals[t[0]]));
min_vertex_dot_product[t[1]] = std::min(min_vertex_dot_product[t[1]], n.dot(vertex_normals[t[1]]));
min_vertex_dot_product[t[2]] = std::min(min_vertex_dot_product[t[2]], n.dot(vertex_normals[t[2]]));
}
}
// lambda to remove face. It flags the face as removed, and updates neighbourhood info
auto remove_face = [&triangles_neighbors, &face_removal_flags](int face_idx, int other_face_idx) {
if (face_idx < 0) {
return;
}
face_removal_flags[face_idx] = true;
Vec3i neighbors = triangles_neighbors[face_idx];
int n_a = neighbors[0] != other_face_idx ? neighbors[0] : neighbors[1];
int n_b = neighbors[2] != other_face_idx ? neighbors[2] : neighbors[1];
if (n_a > 0)
for (int &n : triangles_neighbors[n_a]) {
if (n == face_idx) {
n = n_b;
break;
}
}
if (n_b > 0)
for (int &n : triangles_neighbors[n_b]) {
if (n == face_idx) {
n = n_a;
break;
}
}
};
std::mt19937_64 generator { 27644437 };// default constant seed! so that results are deterministic
std::vector<size_t> face_indices(mesh.indices.size());
for (size_t idx = 0; idx < face_indices.size(); ++idx) {
face_indices[idx] = idx;
}
//tmp face indices used only for swapping
std::vector<size_t> tmp_face_indices(mesh.indices.size());
float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all
float edge_len = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased
while (face_indices.size() > target_triangle_count) {
// simpple func to increase the edge len - if decimation ratio is low, it increases the len up to twice, if decimation ratio is high, increments are low
edge_len = edge_len * (1.0f + 1.0 - decimation_ratio);
float max_edge_len_squared = edge_len * edge_len;
//shuffle the faces and traverse in random order, this MASSIVELY improves the quality of the result
std::shuffle(face_indices.begin(), face_indices.end(), generator);
for (const size_t &face_idx : face_indices) {
if (face_removal_flags[face_idx]) {
// if face already removed from previous collapses, skip (each collapse removes two triangles [at least] )
continue;
}
// look at each edge if it is good candidate for collapse
for (size_t edge_idx = 0; edge_idx < 3; ++edge_idx) {
size_t vertex_index_keep = get_final_index(mesh.indices[face_idx][edge_idx]);
size_t vertex_index_remove = get_final_index(mesh.indices[face_idx][(edge_idx + 1) % 3]);
//check distance, skip long edges
if ((mesh.vertices[vertex_index_keep] - mesh.vertices[vertex_index_remove]).squaredNorm()
> max_edge_len_squared) {
continue;
}
// swap indexes if vertex_index_keep has higher dot product (we want to keep low dot product vertices)
if (min_vertex_dot_product[vertex_index_remove] < min_vertex_dot_product[vertex_index_keep]) {
size_t tmp = vertex_index_keep;
vertex_index_keep = vertex_index_remove;
vertex_index_remove = tmp;
}
//remove vertex
{
// map its index to the index of the kept vertex
vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep];
}
int neighbor_to_remove_face_idx = triangles_neighbors[face_idx][edge_idx];
// remove faces
remove_face(face_idx, neighbor_to_remove_face_idx);
remove_face(neighbor_to_remove_face_idx, face_idx);
// break. this triangle is done
break;
}
}
// filter face_indices, remove those that have been collapsed
size_t prev_size = face_indices.size();
tmp_face_indices.clear();
for (size_t face_idx : face_indices) {
if (!face_removal_flags[face_idx]){
tmp_face_indices.push_back(face_idx);
}
}
face_indices.swap(tmp_face_indices);
decimation_ratio = float(prev_size - face_indices.size()) / float(prev_size);
//std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl;
}
//Extract the result mesh
std::unordered_map<size_t, size_t> final_vertices_mapping;
std::vector<Vec3f> final_vertices;
std::vector<Vec3i> final_indices;
final_indices.reserve(face_indices.size());
for (size_t idx : face_indices) {
Vec3i final_face;
for (size_t i = 0; i < 3; ++i) {
final_face[i] = get_final_index(mesh.indices[idx][i]);
}
if (final_face[0] == final_face[1] || final_face[1] == final_face[2] || final_face[2] == final_face[0]) {
continue; // discard degenerate triangles
}
for (size_t i = 0; i < 3; ++i) {
if (final_vertices_mapping.find(final_face[i]) == final_vertices_mapping.end()) {
final_vertices_mapping[final_face[i]] = final_vertices.size();
final_vertices.push_back(mesh.vertices[final_face[i]]);
}
final_face[i] = final_vertices_mapping[final_face[i]];
}
final_indices.push_back(final_face);
}
mesh.vertices = final_vertices;
mesh.indices = final_indices;
}
} //namespace Slic3r

16
src/libslic3r/ShortEdgeCollapse.hpp Normal file
View file

@ -0,0 +1,16 @@
#ifndef SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_
#define SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_
#include "libslic3r/TriangleMesh.hpp"
namespace Slic3r{
// Decimates the model by collapsing short edges. It starts with very small edges and gradually increases the collapsible length,
// until the target triangle count is reached (the algorithm will certainly undershoot the target count, result will have less triangles than target count)
// The algorithm does not check for triangle flipping, disconnections, self intersections or any other degeneration that can appear during mesh processing.
void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count);
}
#endif /* SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_ */

71
src/libslic3r/TriangleSetSampling.cpp Normal file
View file

@ -0,0 +1,71 @@
#include "TriangleSetSampling.hpp"
#include <map>
#include <random>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
namespace Slic3r {
TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set) {
std::vector<double> triangles_area(triangle_set.indices.size());
tbb::parallel_for(tbb::blocked_range<size_t>(0, triangle_set.indices.size()),
[&triangle_set, &triangles_area](
tbb::blocked_range<size_t> r) {
for (size_t t_idx = r.begin(); t_idx < r.end(); ++t_idx) {
const Vec3f &a = triangle_set.vertices[triangle_set.indices[t_idx].x()];
const Vec3f &b = triangle_set.vertices[triangle_set.indices[t_idx].y()];
const Vec3f &c = triangle_set.vertices[triangle_set.indices[t_idx].z()];
double area = double(0.5 * (b - a).cross(c - a).norm());
triangles_area[t_idx] = area;
}
});
std::map<double, size_t> area_sum_to_triangle_idx;
float area_sum = 0;
for (size_t t_idx = 0; t_idx < triangles_area.size(); ++t_idx) {
area_sum += triangles_area[t_idx];
area_sum_to_triangle_idx[area_sum] = t_idx;
}
std::mt19937_64 mersenne_engine { 27644437 };
// random numbers on interval [0, 1)
std::uniform_real_distribution<double> fdistribution;
auto get_random = [&fdistribution, &mersenne_engine]() {
return Vec3d { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) };
};
std::vector<Vec3d> random_samples(samples_count);
std::generate(random_samples.begin(), random_samples.end(), get_random);
TriangleSetSamples result;
result.total_area = area_sum;
result.positions.resize(samples_count);
result.normals.resize(samples_count);
result.triangle_indices.resize(samples_count);
tbb::parallel_for(tbb::blocked_range<size_t>(0, samples_count),
[&triangle_set, &area_sum_to_triangle_idx, &area_sum, &random_samples, &result](
tbb::blocked_range<size_t> r) {
for (size_t s_idx = r.begin(); s_idx < r.end(); ++s_idx) {
double t_sample = random_samples[s_idx].x() * area_sum;
size_t t_idx = area_sum_to_triangle_idx.upper_bound(t_sample)->second;
double sq_u = std::sqrt(random_samples[s_idx].y());
double v = random_samples[s_idx].z();
Vec3f A = triangle_set.vertices[triangle_set.indices[t_idx].x()];
Vec3f B = triangle_set.vertices[triangle_set.indices[t_idx].y()];
Vec3f C = triangle_set.vertices[triangle_set.indices[t_idx].z()];
result.positions[s_idx] = A * (1 - sq_u) + B * (sq_u * (1 - v)) + C * (v * sq_u);
result.normals[s_idx] = ((B - A).cross(C - B)).normalized();
result.triangle_indices[s_idx] = t_idx;
}
});
return result;
}
}

20
src/libslic3r/TriangleSetSampling.hpp Normal file
View file

@ -0,0 +1,20 @@
#ifndef SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_
#define SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_
#include <admesh/stl.h>
#include "libslic3r/Point.hpp"
namespace Slic3r {
struct TriangleSetSamples {
float total_area;
std::vector<Vec3f> positions;
std::vector<Vec3f> normals;
std::vector<size_t> triangle_indices;
};
TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set);
}
#endif /* SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_ */

23
src/libslic3r/libslic3r.h
View file

@ -24,6 +24,13 @@
#include <cmath> #include <cmath>
#include <type_traits> #include <type_traits>
#ifdef _WIN32
// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740
// Thus it is recommended to use boost::container::deque instead.
#include <boost/container/deque.hpp>
#endif // _WIN32
#include "Technologies.hpp" #include "Technologies.hpp"
#include "Semver.hpp" #include "Semver.hpp"
@ -96,6 +103,16 @@ namespace Slic3r {
extern Semver SEMVER; extern Semver SEMVER;
// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation.
template<class T, class Allocator = std::allocator<T>>
using deque =
#ifdef _WIN32
// Use boost implementation, which allocates blocks of 512 bytes instead of blocks of 8 bytes.
boost::container::deque<T, Allocator>;
#else // _WIN32
std::deque<T, Allocator>;
#endif // _WIN32
template<typename T, typename Q> template<typename T, typename Q>
inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); } inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); }
@ -341,6 +358,12 @@ public:
inline bool empty() const { return size() == 0; } inline bool empty() const { return size() == 0; }
}; };
template<class T, class = FloatingOnly<T>>
constexpr T NaN = std::numeric_limits<T>::quiet_NaN();
constexpr float NaNf = NaN<float>;
constexpr double NaNd = NaN<double>;
} // namespace Slic3r } // namespace Slic3r
#endif #endif