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

Merge branch 'pm_stable_port_seams' into dev_250

Browse Source
This commit is contained in:
Lukas Matena 2022-06-07 11:13:00 +02:00
commit 06a9114f01
21 changed files with 3413 additions and 1496 deletions
Show Stats 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.
template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
struct IndexedTriangleSetDistancer {
@ -453,74 +504,36 @@ namespace detail {
using IndexedFaceType = AIndexedFaceType;
using TreeType = ATreeType;
using VectorType = AVectorType;
using ScalarType = typename VectorType::Scalar;
const std::vector<VertexType> &vertices;
const std::vector<IndexedFaceType> &faces;
const TreeType &tree;
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 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
};
template<typename IndexedTriangleSetDistancerType, typename Scalar>
static inline Scalar squared_distance_to_indexed_triangle_set_recursive(
IndexedTriangleSetDistancerType &distancer,
template<typename IndexedPrimitivesDistancerType, typename Scalar>
static inline Scalar squared_distance_to_indexed_primitives_recursive(
IndexedPrimitivesDistancerType &distancer,
size_t node_idx,
Scalar low_sqr_d,
Scalar up_sqr_d,
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)
return low_sqr_d;
@ -538,13 +551,9 @@ namespace detail {
assert(node.is_valid());
if (node.is_leaf())
{
const auto &triangle = distancer.faces[node.idx];
Vector c_candidate = closest_point_to_triangle<Vector>(
distancer.origin,
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);
Scalar sqr_dist;
Vector c_candidate = distancer.closest_point_to_origin(node.idx, sqr_dist);
set_min(sqr_dist, node.idx, c_candidate);
}
else
{
@ -561,7 +570,7 @@ namespace detail {
{
size_t i_left;
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);
looked_left = true;
};
@ -569,13 +578,13 @@ namespace detail {
{
size_t i_right;
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);
looked_right = true;
};
// 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>()))
look_left();
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()),
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);
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; });
}
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>
{ vertices, faces, tree, point };
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.
@ -759,22 +773,22 @@ inline bool is_any_triangle_in_radius(
const TreeType &tree,
// Point to which the closest point on the indexed triangle set is searched for.
const VectorType &point,
// Maximum distance in which triangle is search for
typename VectorType::Scalar &max_distance)
//Square of maximum distance in which triangle is searched for
typename VectorType::Scalar &max_distance_squared)
{
using Scalar = typename VectorType::Scalar;
auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
{ vertices, faces, tree, point };
size_t hit_idx;
VectorType hit_point = VectorType::Ones() * (std::nan(""));
size_t hit_idx;
VectorType hit_point = VectorType::Ones() * (NaN<typename VectorType::Scalar>);
if(tree.empty())
{
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();
}

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_ */

10
src/libslic3r/CMakeLists.txt
View File

@ -128,10 +128,12 @@ add_library(libslic3r STATIC
GCodeWriter.hpp
Geometry.cpp
Geometry.hpp
Geometry/Bicubic.hpp
Geometry/Circle.cpp
Geometry/Circle.hpp
Geometry/ConvexHull.cpp
Geometry/ConvexHull.hpp
Geometry/Curves.hpp
Geometry/MedialAxis.cpp
Geometry/MedialAxis.hpp
Geometry/Voronoi.hpp
@ -164,6 +166,8 @@ add_library(libslic3r STATIC
MultiPoint.cpp
MultiPoint.hpp
MutablePriorityQueue.hpp
NormalUtils.cpp
NormalUtils.hpp
ObjectID.cpp
ObjectID.hpp
PerimeterGenerator.cpp
@ -203,6 +207,8 @@ add_library(libslic3r STATIC
QuadricEdgeCollapse.cpp
QuadricEdgeCollapse.hpp
Semver.cpp
ShortEdgeCollapse.cpp
ShortEdgeCollapse.hpp
ShortestPath.cpp
ShortestPath.hpp
SLAPrint.cpp
@ -215,6 +221,8 @@ add_library(libslic3r STATIC
SlicesToTriangleMesh.cpp
SlicingAdaptive.cpp
SlicingAdaptive.hpp
Subdivide.cpp
Subdivide.hpp
SupportMaterial.cpp
SupportMaterial.hpp
Surface.cpp
@ -241,6 +249,8 @@ add_library(libslic3r STATIC
Thread.hpp
TriangleSelector.cpp
TriangleSelector.hpp
TriangleSetSampling.cpp
TriangleSetSampling.hpp
MTUtils.hpp
Zipper.hpp
Zipper.cpp

65
src/libslic3r/GCode.cpp
View File

@ -3,7 +3,6 @@
#include "GCode.hpp"
#include "Exception.hpp"
#include "ExtrusionEntity.hpp"
#include "EdgeGrid.hpp"
#include "Geometry/ConvexHull.hpp"
#include "GCode/PrintExtents.hpp"
#include "GCode/WipeTower.hpp"
@ -1348,7 +1347,8 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
print.throw_if_canceled();
// 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);
if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming)) {
// Set initial extruder only after custom start G-code.
@ -2395,7 +2395,6 @@ GCode::LayerResult GCode::process_layer(
} // for objects
// 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)
{
gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ?
@ -2489,9 +2488,9 @@ GCode::LayerResult GCode::process_layer(
//FIXME the following code prints regions in the order they are defined, the path is not optimized in any way.
if (print.config().infill_first) {
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 {
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);
}
// ironing
@ -2626,51 +2625,22 @@ std::string GCode::change_layer(coordf_t print_z)
return gcode;
}
static std::unique_ptr<EdgeGrid::Grid> calculate_layer_edge_grid(const Layer& layer)
{
auto out = make_unique<EdgeGrid::Grid>();
// Create the distance field for a layer below.
const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5);
out->create(layer.lslices, distance_field_resolution);
out->calculate_sdf();
#if 0
{
static int iRun = 0;
BoundingBox bbox = (*lower_layer_edge_grid)->bbox();
bbox.min(0) -= scale_(5.f);
bbox.min(1) -= scale_(5.f);
bbox.max(0) += scale_(5.f);
bbox.max(1) += scale_(5.f);
EdgeGrid::save_png(*(*lower_layer_edge_grid), bbox, scale_(0.1f), debug_out_path("GCode_extrude_loop_edge_grid-%d.png", iRun++));
}
#endif
return out;
}
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
// 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);
// extrude all loops ccw
bool was_clockwise = loop.make_counter_clockwise();
// find the point of the loop that is closest to the current extruder position
// or randomize if requested
Point last_pos = this->last_pos();
if (m_config.spiral_vase) {
if (! m_config.spiral_vase && description == "perimeter") {
assert(m_layer != nullptr);
m_seam_placer.place_seam(m_layer, loop, m_config.external_perimeters_first, this->last_pos());
} else
loop.split_at(last_pos, false);
}
else
m_seam_placer.place_seam(loop, this->last_pos(), m_config.external_perimeters_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;
// if polyline was shorter than the clipping distance we'd get a null polyline, so
@ -2759,14 +2729,14 @@ std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string
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))
return this->extrude_path(*path, description, speed);
else if (const ExtrusionMultiPath* multipath = dynamic_cast<const ExtrusionMultiPath*>(&entity))
return this->extrude_multi_path(*multipath, description, speed);
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
throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()");
return "";
@ -2787,24 +2757,15 @@ std::string GCode::extrude_path(ExtrusionPath path, std::string description, dou
}
// 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;
for (const ObjectByExtruder::Island::Region &region : by_region)
if (! region.perimeters.empty()) {
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)
gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid);
gcode += this->extrude_entity(*ee, "perimeter", -1.);
}
return gcode;
}

6
src/libslic3r/GCode.hpp
View File

@ -271,8 +271,8 @@ private:
void set_extruders(const std::vector<unsigned int> &extruder_ids);
std::string preamble();
std::string change_layer(coordf_t print_z);
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_loop(ExtrusionLoop loop, 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::string extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.);
std::string extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.);
@ -339,7 +339,7 @@ private:
// For sequential print, the instance of the object to be printing has to be defined.
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_support(const ExtrusionEntityCollection &support_fills);

2367
src/libslic3r/GCode/SeamPlacer.cpp
View File

File diff suppressed because it is too large Load Diff

241
src/libslic3r/GCode/SeamPlacer.hpp
View File

@ -3,12 +3,16 @@
#include <optional>
#include <vector>
#include <memory>
#include <atomic>
#include "libslic3r/libslic3r.h"
#include "libslic3r/ExtrusionEntity.hpp"
#include "libslic3r/Polygon.hpp"
#include "libslic3r/PrintConfig.hpp"
#include "libslic3r/BoundingBox.hpp"
#include "libslic3r/AABBTreeIndirect.hpp"
#include "libslic3r/KDTreeIndirect.hpp"
namespace Slic3r {
@ -16,119 +20,166 @@ class PrintObject;
class ExtrusionLoop;
class Print;
class Layer;
namespace EdgeGrid { class Grid; }
namespace EdgeGrid {
class Grid;
}
namespace SeamPlacerImpl {
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();
// ************ 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));
}
// ***************************
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;
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;
};
//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:
void init(const Print& print);
// 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;
// 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.
// arm length used during angles computation
static constexpr float polygon_local_angles_arm_distance = 0.3f;
// max tolerable distance from the previous layer is overhang_distance_tolerance_factor * flow_width
static constexpr float overhang_distance_tolerance_factor = 0.5f;
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);
// 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 = 2.0f; // use much higher angle importance for nearest mode, to combat the visiblity info noise
void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
const EdgeGrid::Grid* lower_layer_edge_grid);
// If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer
static constexpr float enforcer_blocker_distance_tolerance = 0.35f;
// 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;
using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
// 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.27f;
// seam_align_tolerable_dist - if next layer closes point is too far away, break string
static constexpr float seam_align_tolerable_dist = 1.0f;
// if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped.
// this param limits the number of allowed skips
static constexpr size_t seam_align_tolerable_skips = 4;
// minimum number of seams needed in cluster to make alignment happen
static constexpr size_t seam_align_minimum_string_seams = 6;
// points covered by spline; determines number of splines for the given string
static constexpr size_t seam_align_seams_per_segment = 8;
//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:
// 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, bool prefer_nearest);
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());
}
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,
std::optional<std::pair<size_t, size_t>> &out_best_moved_seam,
size_t& out_moved_seams_count) const;
std::optional<std::pair<size_t, size_t>> find_next_seam_in_layer(
const std::vector<PrintObjectSeamData::LayerSeams> &layers,
const std::pair<size_t, size_t> &prev_point_index,
const size_t layer_idx, const float slice_z,
const SeamPlacerImpl::SeamComparator &comparator) const;
};
}
} // namespace Slic3r
#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 */

205
src/libslic3r/Geometry/Curves.hpp Normal file
View File

@ -0,0 +1,205 @@
#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;
Vec3f get_fitted_value(const NumberType value) const {
auto 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 - 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, 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 {
enum class VisitorReturnMask : unsigned int {
CONTINUE_LEFT = 1,
CONTINUE_RIGHT = 2,
STOP = 4,
};
// KD tree for N-dimensional closest point search.
template<size_t ANumDimensions, typename ACoordType, typename ACoordinateFn>
class KDTreeIndirect
{
public:
static constexpr size_t NumDimensions = ANumDimensions;
using CoordinateFn = ACoordinateFn;
using CoordType = ACoordType;
static constexpr size_t NumDimensions = ANumDimensions;
using CoordinateFn = ACoordinateFn;
using CoordType = ACoordType;
// Following could be static constexpr size_t, but that would not link in C++11
enum : size_t {
npos = size_t(-1)
};
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(std::move(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& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; }
void clear() { m_nodes.clear(); }
KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {}
KDTreeIndirect(CoordinateFn coordinate, std::vector<size_t> indices) : coordinate(coordinate) { this->build(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& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; }
void clear() { m_nodes.clear(); }
void build(size_t num_indices)
{
std::vector<size_t> indices;
indices.reserve(num_indices);
for (size_t i = 0; i < num_indices; ++ i)
indices.emplace_back(i);
this->build(std::move(indices));
}
void build(size_t num_indices)
{
std::vector<size_t> indices;
indices.reserve(num_indices);
for (size_t i = 0; i < num_indices; ++ i)
indices.emplace_back(i);
this->build(indices);
}
void build(std::vector<size_t> &&indices)
{
if (indices.empty())
clear();
else {
// Allocate enough memory for a full binary tree.
m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos);
build_recursive(indices, 0, 0, 0, indices.size() - 1);
}
indices.clear();
}
void build(std::vector<size_t> &indices)
{
if (indices.empty())
clear();
else {
// Allocate enough memory for a full binary tree.
m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos);
build_recursive(indices, 0, 0, 0, indices.size() - 1);
}
indices.clear();
}
enum class VisitorReturnMask : unsigned int
{
CONTINUE_LEFT = 1,
CONTINUE_RIGHT = 2,
STOP = 4,
};
template<typename CoordType>
unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
{
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);
}
template<typename CoordType>
unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const
{
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.
template<typename Visitor>
void visit(Visitor &visitor) const
{
// Visitor is supposed to return a bit mask of VisitorReturnMask.
template<typename Visitor>
void visit(Visitor &visitor) const
{
visit_recursive(0, 0, visitor);
}
}
CoordinateFn coordinate;
CoordinateFn coordinate;
private:
// 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)
{
if (left > right)
return;
// 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)
{
if (left > right)
return;
assert(node < m_nodes.size());
assert(node < m_nodes.size());
if (left == right) {
// Insert a node into the balanced tree.
m_nodes[node] = input[left];
return;
}
if (left == right) {
// Insert a node into the balanced tree.
m_nodes[node] = input[left];
return;
}
// Partition the input to left / right pieces of the same length to produce a balanced tree.
size_t center = (left + right) / 2;
partition_input(input, dimension, left, right, center);
// Insert a node into the tree.
m_nodes[node] = input[center];
// Build up the left / right subtrees.
size_t next_dimension = dimension;
if (++ next_dimension == NumDimensions)
next_dimension = 0;
if (center > left)
build_recursive(input, node * 2 + 1, next_dimension, left, center - 1);
build_recursive(input, node * 2 + 2, next_dimension, center + 1, right);
}
// Partition the input to left / right pieces of the same length to produce a balanced tree.
size_t center = (left + right) / 2;
partition_input(input, dimension, left, right, center);
// Insert a node into the tree.
m_nodes[node] = input[center];
// Build up the left / right subtrees.
size_t next_dimension = dimension;
if (++ next_dimension == NumDimensions)
next_dimension = 0;
if (center > left)
build_recursive(input, node * 2 + 1, next_dimension, left, center - 1);
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:
// https://en.wikipedia.org/wiki/Quickselect
// 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",
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) {
size_t center = (left + right) / 2;
CoordType pivot;
{
// Bubble sort the input[left], input[center], input[right], so that a median of the three values
// will end up in input[center].
CoordType left_value = this->coordinate(input[left], dimension);
CoordType center_value = this->coordinate(input[center], dimension);
CoordType right_value = this->coordinate(input[right], dimension);
if (left_value > center_value) {
std::swap(input[left], input[center]);
std::swap(left_value, center_value);
}
if (left_value > right_value) {
std::swap(input[left], input[right]);
right_value = left_value;
}
if (center_value > right_value) {
std::swap(input[center], input[right]);
center_value = right_value;
}
pivot = center_value;
}
if (right <= left + 2)
// The <left, right> interval is already sorted.
break;
size_t i = left;
size_t j = right - 1;
std::swap(input[center], input[j]);
// Partition the set based on the pivot.
for (;;) {
// Skip left points that are already at correct positions.
// Search will certainly stop at position (right - 1), which stores the pivot.
while (this->coordinate(input[++ i], dimension) < pivot) ;
// Skip right points that are already at correct positions.
while (this->coordinate(input[-- j], dimension) > pivot && i < j) ;
if (i >= j)
break;
std::swap(input[i], input[j]);
}
// Restore pivot to the center of the sequence.
std::swap(input[i], input[right - 1]);
// Which side the kth element is in?
if (k < i)
right = i - 1;
else if (k == i)
// Sequence is partitioned, kth element is at its place.
break;
else
left = i + 1;
}
}
// Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
// https://en.wikipedia.org/wiki/Quickselect
// 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",
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) {
size_t center = (left + right) / 2;
CoordType pivot;
{
// Bubble sort the input[left], input[center], input[right], so that a median of the three values
// will end up in input[center].
CoordType left_value = this->coordinate(input[left], dimension);
CoordType center_value = this->coordinate(input[center], dimension);
CoordType right_value = this->coordinate(input[right], dimension);
if (left_value > center_value) {
std::swap(input[left], input[center]);
std::swap(left_value, center_value);
}
if (left_value > right_value) {
std::swap(input[left], input[right]);
right_value = left_value;
}
if (center_value > right_value) {
std::swap(input[center], input[right]);
center_value = right_value;
}
pivot = center_value;
}
if (right <= left + 2)
// The <left, right> interval is already sorted.
break;
size_t i = left;
size_t j = right - 1;
std::swap(input[center], input[j]);
// Partition the set based on the pivot.
for (;;) {
// Skip left points that are already at correct positions.
// Search will certainly stop at position (right - 1), which stores the pivot.
while (this->coordinate(input[++ i], dimension) < pivot) ;
// Skip right points that are already at correct positions.
while (this->coordinate(input[-- j], dimension) > pivot && i < j) ;
if (i >= j)
break;
std::swap(input[i], input[j]);
}
// Restore pivot to the center of the sequence.
std::swap(input[i], input[right - 1]);
// Which side the kth element is in?
if (k < i)
right = i - 1;
else if (k == i)
// Sequence is partitioned, kth element is at its place.
break;
else
left = i + 1;
}
}
template<typename Visitor>
void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const
{
assert(! m_nodes.empty());
if (node >= m_nodes.size() || m_nodes[node] == npos)
return;
template<typename Visitor>
void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const
{
assert(! m_nodes.empty());
if (node >= m_nodes.size() || m_nodes[node] == npos)
return;
// Left / right child node index.
size_t left = node * 2 + 1;
size_t right = left + 1;
unsigned int mask = visitor(m_nodes[node], dimension);
if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) {
size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension;
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT)
visit_recursive(left, next_dimension, visitor);
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT)
visit_recursive(right, next_dimension, visitor);
}
}
// Left / right child node index.
size_t left = node * 2 + 1;
size_t right = left + 1;
unsigned int mask = visitor(m_nodes[node], dimension);
if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) {
size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension;
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT)
visit_recursive(left, next_dimension, visitor);
if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT)
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.
// Returns npos if not found.
template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
size_t find_closest_point(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter)
template<size_t K,
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 {
const KDTreeIndirectType &kdtree;
const PointType &point;
const FilterFn filter;
size_t min_idx = KDTreeIndirectType::npos;
CoordType min_dist = std::numeric_limits<CoordType>::max();
struct Visitor
{
const Tree &kdtree;
const PointType &point;
const FilterFn filter;
Visitor(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) : kdtree(kdtree), point(point), 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 = 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);
std::array<std::pair<size_t, CoordT>, K> results;
kdtree.visit(visitor);
return visitor.min_idx;
Visitor(const Tree &kdtree, const PointType &point, FilterFn filter)
: 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>
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

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_

14
src/libslic3r/QuadricEdgeCollapse.cpp
View File

@ -595,8 +595,7 @@ bool QuadricEdgeCollapse::is_flipped(const Vec3f & new_vertex,
const EdgeInfos & e_infos,
const indexed_triangle_set &its)
{
static const float thr_pos = 1.0f - std::numeric_limits<float>::epsilon();
static const float thr_neg = -thr_pos;
static const float triangle_beauty_threshold = 1.0f - std::numeric_limits<float>::epsilon();
static const float dot_thr = 0.2f; // Value from simplify mesh cca 80 DEG
// for each vertex triangles
@ -615,7 +614,16 @@ bool QuadricEdgeCollapse::is_flipped(const Vec3f & new_vertex,
d2.normalize();
float dot = d1.dot(d2);
if (dot > thr_pos || dot < thr_neg) return true;
if (dot > triangle_beauty_threshold || dot < -triangle_beauty_threshold) { // OK, the new triangle is suspiciously ugly, but it can still be better than the original
const Vec3f &v_orig = its.vertices[t[(e_info.edge) % 3]];
Vec3f d1_orig = vf - v_orig;
d1_orig.normalize();
Vec3f d2_orig = vs - v_orig;
d2_orig.normalize();
if (std::fabs(d1_orig.dot(d2_orig)) < std::fabs(dot)) { // original was not that ugly, so return flipped
return true;
} // else original triangle was worse than the new, so don't discard the new yet
}
// IMPROVE: propagate new normal
Vec3f n = d1.cross(d2);
n.normalize();

186
src/libslic3r/SLA/bicubic.h
View File

@ -1,186 +0,0 @@
#ifndef BICUBIC_HPP
#define BICUBIC_HPP
#include <algorithm>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
namespace Slic3r {
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 CubicKernel
{
typedef typename KERNEL KernelInternal;
typedef typename KERNEL::FloatType FloatType;
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
typedef CubicKernel<BicubicInternal::LinearKernel<float>> LinearKernelf;
typedef CubicKernel<BicubicInternal::LinearKernel<double>> LinearKerneld;
// Catmul-Rom splines
typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>> CubicCatmulRomKernelf;
typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>> CubicCatmulRomKerneld;
typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>> CubicInterpolationKernelf;
typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>> CubicInterpolationKerneld;
// Cubic B-splines
typedef CubicKernel<BicubicInternal::CubicBSplineKernel<float>> CubicBSplineKernelf;
typedef CubicKernel<BicubicInternal::CubicBSplineKernel<double>> CubicBSplineKerneld;
template<typename KERNEL, typename Derived>
static float cubic_interpolate(const Eigen::ArrayBase<Derived> &F, const typename KERNEL::FloatType pt, const typename KERNEL::FloatType dx)
{
typedef typename KERNEL::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 KERNEL::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
}
// Transition region. Extend with a constant function.
auto f = [&F, w](x) { return F[BicubicInternal::clamp(x, 0, w - 1)]; }
return KERNEL::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, const typename KERNEL::FloatType dx)
{
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 Slic3r
#endif /* BICUBIC_HPP */

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_ */

218
src/libslic3r/Subdivide.cpp Normal file
View File

@ -0,0 +1,218 @@
#include "Subdivide.hpp"
#include "Point.hpp"
namespace Slic3r{
indexed_triangle_set its_subdivide(
const indexed_triangle_set &its, float max_length)
{
// same order as key order in Edge Divides
struct VerticesSequence
{
size_t start_index;
bool positive_order;
VerticesSequence(size_t start_index, bool positive_order = true)
: start_index(start_index), positive_order(positive_order){}
};
// vertex index small, big vertex index from key.first to key.second
using EdgeDivides = std::map<std::pair<size_t, size_t>, VerticesSequence>;
struct Edges
{
Vec3f data[3];
Vec3f lengths;
Edges(const Vec3crd &indices, const std::vector<Vec3f> &vertices)
: lengths(-1.f,-1.f,-1.f)
{
const Vec3f &v0 = vertices[indices[0]];
const Vec3f &v1 = vertices[indices[1]];
const Vec3f &v2 = vertices[indices[2]];
data[0] = v0 - v1;
data[1] = v1 - v2;
data[2] = v2 - v0;
}
float abs_sum(const Vec3f &v)
{
return abs(v[0]) + abs(v[1]) + abs(v[2]);
}
bool is_dividable(const float& max_length) {
Vec3f sum(abs_sum(data[0]), abs_sum(data[1]), abs_sum(data[2]));
Vec3i biggest_index = (sum[0] > sum[1]) ?
((sum[0] > sum[2]) ?
((sum[2] > sum[1]) ?
Vec3i(0, 2, 1) :
Vec3i(0, 1, 2)) :
Vec3i(2, 0, 1)) :
((sum[1] > sum[2]) ?
((sum[2] > sum[0]) ?
Vec3i(1, 2, 0) :
Vec3i(1, 0, 2)) :
Vec3i(2, 1, 0));
for (int i = 0; i < 3; i++) {
int index = biggest_index[i];
if (sum[index] <= max_length) return false;
lengths[index] = data[index].norm();
if (lengths[index] <= max_length) continue;
// calculate rest of lengths
for (int j = i + 1; j < 3; j++) {
index = biggest_index[j];
lengths[index] = data[index].norm();
}
return true;
}
return false;
}
};
struct TriangleLengths
{
Vec3crd indices;
Vec3f l; // lengths
TriangleLengths(const Vec3crd &indices, const Vec3f &lengths)
: indices(indices), l(lengths)
{}
int get_divide_index(float max_length) {
if (l[0] > l[1] && l[0] > l[2]) {
if (l[0] > max_length) return 0;
} else if (l[1] > l[2]) {
if (l[1] > max_length) return 1;
} else {
if (l[2] > max_length) return 2;
}
return -1;
}
// divide triangle add new vertex to vertices
std::pair<TriangleLengths, TriangleLengths> divide(
int divide_index, float max_length,
std::vector<Vec3f> &vertices,
EdgeDivides &edge_divides)
{
// index to lengths and indices
size_t i0 = divide_index;
size_t i1 = (divide_index + 1) % 3;
size_t vi0 = indices[i0];
size_t vi1 = indices[i1];
std::pair<size_t, size_t> key(vi0, vi1);
bool key_swap = false;
if (key.first > key.second) {
std::swap(key.first, key.second);
key_swap = true;
}
float length = l[divide_index];
size_t count_edge_vertices = static_cast<size_t>(floor(length / max_length));
float count_edge_segments = static_cast<float>(count_edge_vertices + 1);
auto it = edge_divides.find(key);
if (it == edge_divides.end()) {
// Create new vertices
VerticesSequence new_vs(vertices.size());
Vec3f vf = vertices[key.first]; // copy
const Vec3f &vs = vertices[key.second];
Vec3f dir = vs - vf;
for (size_t i = 1; i <= count_edge_vertices; ++i) {
float ratio = i / count_edge_segments;
vertices.push_back(vf + dir * ratio);
}
bool success;
std::tie(it,success) = edge_divides.insert({key, new_vs});
assert(success);
}
const VerticesSequence &vs = it->second;
int index_offset = count_edge_vertices/2;
size_t i2 = (divide_index + 2) % 3;
if (count_edge_vertices % 2 == 0 && key_swap == (l[i1] < l[i2])) {
--index_offset;
}
int sign = (vs.positive_order) ? 1 : -1;
size_t new_index = vs.start_index + sign*index_offset;
size_t vi2 = indices[i2];
const Vec3f &v2 = vertices[vi2];
Vec3f new_edge = v2 - vertices[new_index];
float new_len = new_edge.norm();
float ratio = (1 + index_offset) / count_edge_segments;
float len1 = l[i0] * ratio;
float len2 = l[i0] - len1;
if (key_swap) std::swap(len1, len2);
Vec3crd indices1(vi0, new_index, vi2);
Vec3f lengths1(len1, new_len, l[i2]);
Vec3crd indices2(new_index, vi1, vi2);
Vec3f lengths2(len2, l[i1], new_len);
// append key for divided edge when neccesary
if (index_offset > 0) {
std::pair<size_t, size_t> new_key(key.first, new_index);
bool new_key_swap = false;
if (new_key.first > new_key.second) {
std::swap(new_key.first, new_key.second);
new_key_swap = true;
}
if (edge_divides.find(new_key) == edge_divides.end()) {
// insert new
edge_divides.insert({new_key, (new_key_swap) ?
VerticesSequence(new_index - sign, !vs.positive_order)
: vs});
}
}
if (index_offset < int(count_edge_vertices)-1) {
std::pair<size_t, size_t> new_key(new_index, key.second);
bool new_key_swap = false;
if (new_key.first > new_key.second) {
std::swap(new_key.first, new_key.second);
new_key_swap = true;
}
// bad order
if (edge_divides.find(new_key) == edge_divides.end()) {
edge_divides.insert({new_key, (new_key_swap) ?
VerticesSequence(vs.start_index + sign*(count_edge_vertices-1), !vs.positive_order)
: VerticesSequence(new_index + sign, vs.positive_order)});
}
}
return {TriangleLengths(indices1, lengths1),
TriangleLengths(indices2, lengths2)};
}
};
indexed_triangle_set result;
result.indices.reserve(its.indices.size());
const std::vector<Vec3f> &vertices = its.vertices;
result.vertices = vertices; // copy
std::queue<TriangleLengths> tls;
EdgeDivides edge_divides;
for (const Vec3crd &indices : its.indices) {
Edges edges(indices, vertices);
// speed up only sum not sqrt is apply
if (!edges.is_dividable(max_length)) {
// small triangle
result.indices.push_back(indices);
continue;
}
TriangleLengths tl(indices, edges.lengths);
do {
int divide_index = tl.get_divide_index(max_length);
if (divide_index < 0) {
// no dividing
result.indices.push_back(tl.indices);
if (tls.empty()) break;
tl = tls.front(); // copy
tls.pop();
} else {
auto [tl1, tl2] = tl.divide(divide_index, max_length,
result.vertices, edge_divides);
tl = tl1;
tls.push(tl2);
}
} while (true);
}
return result;
}
}

12
src/libslic3r/Subdivide.hpp Normal file
View File

@ -0,0 +1,12 @@
#ifndef libslic3r_Subdivide_hpp_
#define libslic3r_Subdivide_hpp_
#include "TriangleMesh.hpp"
namespace Slic3r {
indexed_triangle_set its_subdivide(const indexed_triangle_set &its, float max_length);
}
#endif //libslic3r_Subdivide_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

@ -23,6 +23,13 @@
#include <cmath>
#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 "Semver.hpp"
@ -73,6 +80,16 @@ namespace Slic3r {
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>
inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); }
@ -314,6 +331,12 @@ public:
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
#endif