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

Merge pull request #8242 from prusa3d/pm_occlusion_guided_seams

Browse Source 
Implementation of occlusion guided seam placement - balance between visibility and angle sharpness.
Improvements in seam drawing and alignment, via cubic splines.
Spline/Polynomial fitting functions.
Refactoring of Point structure - angle functions.
Triangular mesh subdivision function.
Various small additions/refactoring.
! Macro which replaces std::deque with boost version on Windows builds.
This commit is contained in:
Pavel Mikuš 2022-04-25 17:53:30 +02:00 committed by GitHub
commit 1b1f523e3b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
31 changed files with 2603 additions and 1482 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/seam_placement/corner_penalty_function.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 34 KiB

15
src/libslic3r/AABBTreeIndirect.hpp
View File

@ -709,10 +709,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();
@ -759,8 +764,8 @@ 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>
@ -774,7 +779,7 @@ inline bool is_any_triangle_in_radius(
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_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
return hit_point.allFinite();
}

4
src/libslic3r/CMakeLists.txt
View File

@ -138,10 +138,12 @@ set(SLIC3R_SOURCES
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
@ -225,6 +227,8 @@ set(SLIC3R_SOURCES
SlicesToTriangleMesh.cpp
SlicingAdaptive.cpp
SlicingAdaptive.hpp
Subdivide.cpp
Subdivide.hpp
SupportMaterial.cpp
SupportMaterial.hpp
Surface.cpp

64
src/libslic3r/Color.hpp
View File

@ -4,6 +4,8 @@
#include <array>
#include <algorithm>
#include "Point.hpp"
namespace Slic3r {
class ColorRGB
@ -133,41 +135,57 @@ public:
static const ColorRGBA Z() { return { 0.0f, 0.0f, 0.75f, 1.0f }; }
};
extern ColorRGB operator * (float value, const ColorRGB& other);
extern ColorRGBA operator * (float value, const ColorRGBA& other);
ColorRGB operator * (float value, const ColorRGB& other);
ColorRGBA operator * (float value, const ColorRGBA& other);
extern ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t);
extern ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t);
ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t);
ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t);
extern ColorRGB complementary(const ColorRGB& color);
extern ColorRGBA complementary(const ColorRGBA& color);
ColorRGB complementary(const ColorRGB& color);
ColorRGBA complementary(const ColorRGBA& color);
extern ColorRGB saturate(const ColorRGB& color, float factor);
extern ColorRGBA saturate(const ColorRGBA& color, float factor);
ColorRGB saturate(const ColorRGB& color, float factor);
ColorRGBA saturate(const ColorRGBA& color, float factor);
extern ColorRGB opposite(const ColorRGB& color);
extern ColorRGB opposite(const ColorRGB& a, const ColorRGB& b);
ColorRGB opposite(const ColorRGB& color);
ColorRGB opposite(const ColorRGB& a, const ColorRGB& b);
extern bool can_decode_color(const std::string& color);
bool can_decode_color(const std::string& color);
extern bool decode_color(const std::string& color_in, ColorRGB& color_out);
extern bool decode_color(const std::string& color_in, ColorRGBA& color_out);
bool decode_color(const std::string& color_in, ColorRGB& color_out);
bool decode_color(const std::string& color_in, ColorRGBA& color_out);
extern bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGB>& colors_out);
extern bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGBA>& colors_out);
bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGB>& colors_out);
bool decode_colors(const std::vector<std::string>& colors_in, std::vector<ColorRGBA>& colors_out);
extern std::string encode_color(const ColorRGB& color);
extern std::string encode_color(const ColorRGBA& color);
std::string encode_color(const ColorRGB& color);
std::string encode_color(const ColorRGBA& color);
extern ColorRGB to_rgb(const ColorRGBA& other_rgba);
extern ColorRGBA to_rgba(const ColorRGB& other_rgb);
extern ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha);
ColorRGB to_rgb(const ColorRGBA& other_rgba);
ColorRGBA to_rgba(const ColorRGB& other_rgb);
ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha);
extern ColorRGBA picking_decode(unsigned int id);
extern unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b);
// Color mapping of a value into RGB false colors.
inline Vec3f value_to_rgbf(float minimum, float maximum, float value)
{
float ratio = 2.0f * (value - minimum) / (maximum - minimum);
float b = std::max(0.0f, (1.0f - ratio));
float r = std::max(0.0f, (ratio - 1.0f));
float g = 1.0f - b - r;
return Vec3f { r, g, b };
}
// Color mapping of a value into RGB false colors.
inline Vec3i value_to_rgbi(float minimum, float maximum, float value)
{
return (value_to_rgbf(minimum, maximum, value) * 255).cast<int>();
}
ColorRGBA picking_decode(unsigned int id);
unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b);
// Produce an alpha channel checksum for the red green blue components. The alpha channel may then be used to verify, whether the rgb components
// were not interpolated by alpha blending or multi sampling.
extern unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue);
unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue);
} // namespace Slic3r

134
src/libslic3r/GCode.cpp
View File

@ -2312,7 +2312,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) ?
@ -2337,7 +2336,7 @@ GCode::LayerResult GCode::process_layer(
path.mm3_per_mm = mm3_per_mm;
}
//FIXME using the support_material_speed of the 1st object printed.
gcode += this->extrude_loop(loop, "skirt", m_config.support_material_speed.value);
gcode += this->extrude_loop(loop, "skirt"sv, m_config.support_material_speed.value);
}
m_avoid_crossing_perimeters.use_external_mp(false);
// Allow a straight travel move to the first object point if this is the first layer (but don't in next layers).
@ -2350,7 +2349,7 @@ GCode::LayerResult GCode::process_layer(
this->set_origin(0., 0.);
m_avoid_crossing_perimeters.use_external_mp();
for (const ExtrusionEntity *ee : print.brim().entities) {
gcode += this->extrude_entity(*ee, "brim", m_config.support_material_speed.value);
gcode += this->extrude_entity(*ee, "brim"sv, m_config.support_material_speed.value);
}
m_brim_done = true;
m_avoid_crossing_perimeters.use_external_mp(false);
@ -2406,9 +2405,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
@ -2543,51 +2542,28 @@ 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;
static const auto comment_perimeter = "perimeter"sv;
// Comparing string_view pointer & length for speed.
static inline bool comment_is_perimeter(const std::string_view comment) {
return comment.data() == comment_perimeter.data() && comment.size() == comment_perimeter.size();
}
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, const std::string_view 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 && comment_is_perimeter(description)) {
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
@ -2607,11 +2583,9 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
// extrude along the path
std::string gcode;
for (ExtrusionPaths::iterator path = paths.begin(); path != paths.end(); ++path) {
// description += ExtrusionLoop::role_to_string(loop.loop_role());
// description += ExtrusionEntity::role_to_string(path->role);
path->simplify(m_scaled_resolution);
gcode += this->_extrude(*path, description, speed);
for (ExtrusionPath &path : paths) {
path.simplify(m_scaled_resolution);
gcode += this->_extrude(path, description, speed);
}
// reset acceleration
@ -2626,18 +2600,19 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
// the side depends on the original winding order of the polygon (left for contours, right for holes)
//FIXME improve the algorithm in case the loop is tiny.
//FIXME improve the algorithm in case the loop is split into segments with a low number of points (see the Point b query).
Point a = paths.front().polyline.points[1]; // second point
Point b = *(paths.back().polyline.points.end()-3); // second to last point
// Angle from the 2nd point to the last point.
double angle_inside = angle(paths.front().polyline.points[1] - paths.front().first_point(),
*(paths.back().polyline.points.end()-3) - paths.front().first_point());
assert(angle_inside >= -M_PI && angle_inside <= M_PI);
// 3rd of this angle will be taken, thus make the angle monotonic before interpolation.
if (was_clockwise) {
// swap points
Point c = a; a = b; b = c;
if (angle_inside > 0)
angle_inside -= 2.0 * M_PI;
} else {
if (angle_inside < 0)
angle_inside += 2.0 * M_PI;
}
double angle = paths.front().first_point().ccw_angle(a, b) / 3;
// turn left if contour, turn right if hole
if (was_clockwise) angle *= -1;
// create the destination point along the first segment and rotate it
// we make sure we don't exceed the segment length because we don't know
// the rotation of the second segment so we might cross the object boundary
@ -2649,7 +2624,8 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
// Shift by no more than a nozzle diameter.
//FIXME Hiding the seams will not work nicely for very densely discretized contours!
Point pt = ((nd * nd >= l2) ? p2 : (p1 + v * (nd / sqrt(l2)))).cast<coord_t>();
pt.rotate(angle, paths.front().polyline.points.front());
// Rotate pt inside around the seam point.
pt.rotate(angle_inside / 3., paths.front().polyline.points.front());
// generate the travel move
gcode += m_writer.travel_to_xy(this->point_to_gcode(pt), "move inwards before travel");
}
@ -2657,13 +2633,11 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
return gcode;
}
std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string description, double speed)
std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed)
{
// extrude along the path
std::string gcode;
for (ExtrusionPath path : multipath.paths) {
// description += ExtrusionLoop::role_to_string(loop.loop_role());
// description += ExtrusionEntity::role_to_string(path->role);
path.simplify(m_scaled_resolution);
gcode += this->_extrude(path, description, speed);
}
@ -2676,22 +2650,21 @@ 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, const std::string_view 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 "";
}
std::string GCode::extrude_path(ExtrusionPath path, std::string description, double speed)
std::string GCode::extrude_path(ExtrusionPath path, std::string_view description, double speed)
{
// description += ExtrusionEntity::role_to_string(path.role());
path.simplify(m_scaled_resolution);
std::string gcode = this->_extrude(path, description, speed);
if (m_wipe.enable) {
@ -2704,24 +2677,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, comment_perimeter, -1.);
}
return gcode;
}
@ -2731,7 +2695,7 @@ std::string GCode::extrude_infill(const Print &print, const std::vector<ObjectBy
{
std::string gcode;
ExtrusionEntitiesPtr extrusions;
const char* extrusion_name = ironing ? "ironing" : "infill";
const auto extrusion_name = ironing ? "ironing"sv : "infill"sv;
for (const ObjectByExtruder::Island::Region &region : by_region)
if (! region.infills.empty()) {
extrusions.clear();
@ -2757,8 +2721,8 @@ std::string GCode::extrude_infill(const Print &print, const std::vector<ObjectBy
std::string GCode::extrude_support(const ExtrusionEntityCollection &support_fills)
{
static constexpr const char *support_label = "support material";
static constexpr const char *support_interface_label = "support material interface";
static constexpr const auto support_label = "support material"sv;
static constexpr const auto support_interface_label = "support material interface"sv;
std::string gcode;
if (! support_fills.entities.empty()) {
@ -2767,7 +2731,7 @@ std::string GCode::extrude_support(const ExtrusionEntityCollection &support_fill
for (const ExtrusionEntity *ee : support_fills.entities) {
ExtrusionRole role = ee->role();
assert(role == erSupportMaterial || role == erSupportMaterialInterface);
const char *label = (role == erSupportMaterial) ? support_label : support_interface_label;
const auto label = (role == erSupportMaterial) ? support_label : support_interface_label;
const double speed = (role == erSupportMaterial) ? support_speed : support_interface_speed;
const ExtrusionPath *path = dynamic_cast<const ExtrusionPath*>(ee);
if (path)
@ -2855,20 +2819,18 @@ void GCode::GCodeOutputStream::write_format(const char* format, ...)
va_end(args);
}
std::string GCode::_extrude(const ExtrusionPath &path, std::string description, double speed)
std::string GCode::_extrude(const ExtrusionPath &path, const std::string_view description, double speed)
{
std::string gcode;
if (is_bridge(path.role()))
description += " (bridge)";
const std::string_view description_bridge = is_bridge(path.role()) ? " (bridge)"sv : ""sv;
// go to first point of extrusion path
if (!m_last_pos_defined || m_last_pos != path.first_point()) {
gcode += this->travel_to(
path.first_point(),
path.role(),
"move to first " + description + " point"
);
std::string comment = "move to first ";
comment += description;
comment += description_bridge;
comment += " point";
gcode += this->travel_to(path.first_point(), path.role(), comment);
}
// compensate retraction
@ -3006,7 +2968,11 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
gcode += m_writer.set_speed(F, "", comment);
double path_length = 0.;
{
std::string comment = m_config.gcode_comments ? description : "";
std::string comment;
if (m_config.gcode_comments) {
comment = description;
comment += description_bridge;
}
Vec2d prev = this->point_to_gcode_quantized(path.polyline.points.front());
auto it = path.polyline.points.begin();
auto end = path.polyline.points.end();

12
src/libslic3r/GCode.hpp
View File

@ -274,10 +274,10 @@ 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_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.);
std::string extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.);
std::string extrude_entity(const ExtrusionEntity &entity, const std::string_view description, double speed = -1.);
std::string extrude_loop(ExtrusionLoop loop, const std::string_view description, double speed = -1.);
std::string extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed = -1.);
std::string extrude_path(ExtrusionPath path, const std::string_view description, double speed = -1.);
// Extruding multiple objects with soluble / non-soluble / combined supports
// on a multi-material printer, trying to minimize tool switches.
@ -342,7 +342,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);
@ -428,7 +428,7 @@ private:
// Processor
GCodeProcessor m_processor;
std::string _extrude(const ExtrusionPath &path, std::string description = "", double speed = -1);
std::string _extrude(const ExtrusionPath &path, const std::string_view description, double speed = -1);
void print_machine_envelope(GCodeOutputStream &file, Print &print);
void _print_first_layer_bed_temperature(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait);
void _print_first_layer_extruder_temperatures(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait);

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

File diff suppressed because it is too large Load Diff

221
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,148 @@ class PrintObject;
class ExtrusionLoop;
class Print;
class Layer;
namespace EdgeGrid { class Grid; }
namespace EdgeGrid {
class Grid;
}
class SeamHistory {
public:
SeamHistory() { clear(); }
std::optional<Point> get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb);
void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb);
void clear();
namespace SeamPlacerImpl {
private:
struct SeamPoint {
Point m_pos;
BoundingBox m_island_bb;
};
struct GlobalModelInfo;
struct SeamComparator;
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;
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 FaceVisibilityInfo {
float visibility;
};
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);
static constexpr size_t raycasting_decimation_target_triangle_count = 10000;
static constexpr float raycasting_subdivision_target_length = 2.0f;
//square of number of rays per triangle
static constexpr size_t sqr_rays_per_triangle = 7;
// 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.5f;
// increases angle importance at the cost of deacreasing visibility info importance. must be > 0
static constexpr float additional_angle_importance = 0.6f;
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);
// 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;
void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter,
const EdgeGrid::Grid* lower_layer_edge_grid);
// 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.5f;
// 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;
using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
void init(const Print &print);
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);
bool find_next_seam_in_layer(
const std::vector<PrintObjectSeamData::LayerSeams> &layers,
std::pair<size_t, size_t> &last_point_indexes,
const size_t layer_idx, const float slice_z,
const SeamPlacerImpl::SeamComparator &comparator,
std::vector<std::pair<size_t, size_t>> &seam_string) const;
};
}
} // namespace Slic3r
#endif // libslic3r_SeamPlacer_hpp_

6
src/libslic3r/Geometry.hpp
View File

@ -36,9 +36,9 @@ enum Orientation
static inline Orientation orient(const Point &a, const Point &b, const Point &c)
{
static_assert(sizeof(coord_t) * 2 == sizeof(int64_t), "orient works with 32 bit coordinates");
int64_t u = int64_t(b(0)) * int64_t(c(1)) - int64_t(b(1)) * int64_t(c(0));
int64_t v = int64_t(a(0)) * int64_t(c(1)) - int64_t(a(1)) * int64_t(c(0));
int64_t w = int64_t(a(0)) * int64_t(b(1)) - int64_t(a(1)) * int64_t(b(0));
int64_t u = int64_t(b.x()) * int64_t(c.y()) - int64_t(b.y()) * int64_t(c.x());
int64_t v = int64_t(a.x()) * int64_t(c.y()) - int64_t(a.y()) * int64_t(c.x());
int64_t w = int64_t(a.x()) * int64_t(b.y()) - int64_t(a.y()) * int64_t(b.x());
int64_t d = u - v + w;
return (d > 0) ? ORIENTATION_CCW : ((d == 0) ? ORIENTATION_COLINEAR : ORIENTATION_CW);
}

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

9
src/libslic3r/Geometry/ConvexHull.cpp
View File

@ -1,6 +1,7 @@
#include "libslic3r.h"
#include "ConvexHull.hpp"
#include "BoundingBox.hpp"
#include "../Geometry.hpp"
#include <boost/multiprecision/integer.hpp>
@ -19,13 +20,13 @@ Polygon convex_hull(Points pts)
hull.points.resize(2 * n);
// Build lower hull
for (int i = 0; i < n; ++ i) {
while (k >= 2 && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
while (k >= 2 && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
-- k;
hull[k ++] = pts[i];
}
// Build upper hull
for (int i = n-2, t = k+1; i >= 0; i--) {
while (k >= t && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
while (k >= t && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
-- k;
hull[k ++] = pts[i];
}
@ -58,7 +59,7 @@ Pointf3s convex_hull(Pointf3s points)
Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
if (p.ccw(k2, k1) <= 0)
if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
--k;
else
break;
@ -76,7 +77,7 @@ Pointf3s convex_hull(Pointf3s points)
Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
if (p.ccw(k2, k1) <= 0)
if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
--k;
else
break;

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

47
src/libslic3r/KDTreeIndirect.hpp
View File

@ -231,6 +231,53 @@ size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& poi
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;
});
}
} // namespace Slic3r
#endif /* slic3r_KDTreeIndirect_hpp_ */

1
src/libslic3r/Line.hpp
View File

@ -113,7 +113,6 @@ public:
Vector vector() const { return this->b - this->a; }
Vector normal() const { return Vector((this->b(1) - this->a(1)), -(this->b(0) - this->a(0))); }
bool intersection(const Line& line, Point* intersection) const;
double ccw(const Point& point) const { return point.ccw(*this); }
// Clip a line with a bounding box. Returns false if the line is completely outside of the bounding box.
bool clip_with_bbox(const BoundingBox &bbox);
// Extend the line from both sides by an offset.

30
src/libslic3r/Point.cpp
View File

@ -108,36 +108,6 @@ bool Point::nearest_point(const Points &points, Point* point) const
return true;
}
/* Three points are a counter-clockwise turn if ccw > 0, clockwise if
* ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
* gives the signed area of the triangle formed by p1, p2 and this point.
* In other words it is the 2D cross product of p1-p2 and p1-this, i.e.
* z-component of their 3D cross product.
* We return double because it must be big enough to hold 2*max(|coordinate|)^2
*/
double Point::ccw(const Point &p1, const Point &p2) const
{
static_assert(sizeof(coord_t) == 4, "Point::ccw() requires a 32 bit coord_t");
return cross2((p2 - p1).cast<int64_t>(), (*this - p1).cast<int64_t>());
// return cross2((p2 - p1).cast<double>(), (*this - p1).cast<double>());
}
double Point::ccw(const Line &line) const
{
return this->ccw(line.a, line.b);
}
// returns the CCW angle between this-p1 and this-p2
// i.e. this assumes a CCW rotation from p1 to p2 around this
double Point::ccw_angle(const Point &p1, const Point &p2) const
{
//FIXME this calculates an atan2 twice! Project one vector into the other!
double angle = atan2(p1.x() - (*this).x(), p1.y() - (*this).y())
- atan2(p2.x() - (*this).x(), p2.y() - (*this).y());
// we only want to return only positive angles
return angle <= 0 ? angle + 2*PI : angle;
}
Point Point::projection_onto(const MultiPoint &poly) const
{
Point running_projection = poly.first_point();

39
src/libslic3r/Point.hpp
View File

@ -28,6 +28,9 @@ using Mat = Eigen::Matrix<T, N, M, Eigen::DontAlign, N, M>;
template<int N, class T> using Vec = Mat<N, 1, T>;
template<typename NumberType>
using DynVec = Eigen::Matrix<NumberType, Eigen::Dynamic, 1>;
// Eigen types, to replace the Slic3r's own types in the future.
// Vector types with a fixed point coordinate base type.
using Vec2crd = Eigen::Matrix<coord_t, 2, 1, Eigen::DontAlign>;
@ -76,24 +79,35 @@ inline const auto &identity3d = identity<3, double>;
inline bool operator<(const Vec2d &lhs, const Vec2d &rhs) { return lhs.x() < rhs.x() || (lhs.x() == rhs.x() && lhs.y() < rhs.y()); }
template<int Options>
int32_t cross2(const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v2) = delete;
template<typename T, int Options>
inline T cross2(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v2)
{
return v1.x() * v2.y() - v1.y() * v2.x();
}
// Cross product of two 2D vectors.
// None of the vectors may be of int32_t type as the result would overflow.
template<typename Derived, typename Derived2>
inline typename Derived::Scalar cross2(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2)
{
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
static_assert(! std::is_same<typename Derived::Scalar, int32_t>::value, "cross2(): Scalar type must not be int32_t, otherwise the cross product would overflow.");
static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value, "cross2(): Scalar types of 1st and 2nd operand must be equal.");
return v1.x() * v2.y() - v1.y() * v2.x();
}
template<typename T, int Options>
inline Eigen::Matrix<T, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v) { return Eigen::Matrix<T, 2, 1, Eigen::DontAlign>(- v.y(), v.x()); }
// 2D vector perpendicular to the argument.
template<typename Derived>
inline Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Derived> &v)
{
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "perp(): parameter is not a 2D vector");
return { - v.y(), v.x() };
}
// 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));
}
template<class T, int N, int Options>
Eigen::Matrix<T, 2, 1, Eigen::DontAlign> to_2d(const Eigen::MatrixBase<Eigen::Matrix<T, N, 1, Options>> &ptN) { return { ptN.x(), ptN.y() }; }
@ -165,9 +179,6 @@ public:
int nearest_point_index(const PointConstPtrs &points) const;
int nearest_point_index(const PointPtrs &points) const;
bool nearest_point(const Points &points, Point* point) const;
double ccw(const Point &p1, const Point &p2) const;
double ccw(const Line &line) const;
double ccw_angle(const Point &p1, const Point &p2) const;
Point projection_onto(const MultiPoint &poly) const;
Point projection_onto(const Line &line) const;
};

84
src/libslic3r/Polygon.cpp
View File

@ -171,50 +171,56 @@ Point Polygon::centroid() const
return Point(Vec2d(c / (3. * area_sum)));
}
// find all concave vertices (i.e. having an internal angle greater than the supplied angle)
// (external = right side, thus we consider ccw orientation)
Points Polygon::concave_points(double angle) const
// Filter points from poly to the output with the help of FilterFn.
// filter function receives two vectors:
// v1: this_point - previous_point
// v2: next_point - this_point
// and returns true if the point is to be copied to the output.
template<typename FilterFn>
Points filter_points_by_vectors(const Points &poly, FilterFn filter)
{
Points points;
angle = 2. * PI - angle + EPSILON;
// Last point is the first point visited.
Point p1 = poly.back();
// Previous vector to p1.
Vec2d v1 = (p1 - *(poly.end() - 2)).cast<double>();
// check whether first point forms a concave angle
if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) <= angle)
points.push_back(this->points.front());
// check whether points 1..(n-1) form concave angles
for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++ p)
if (p->ccw_angle(*(p-1), *(p+1)) <= angle)
points.push_back(*p);
// check whether last point forms a concave angle
if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) <= angle)
points.push_back(this->points.back());
return points;
}
// find all convex vertices (i.e. having an internal angle smaller than the supplied angle)
// (external = right side, thus we consider ccw orientation)
Points Polygon::convex_points(double angle) const
{
Points points;
angle = 2*PI - angle - EPSILON;
// check whether first point forms a convex angle
if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) >= angle)
points.push_back(this->points.front());
// check whether points 1..(n-1) form convex angles
for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++p) {
if (p->ccw_angle(*(p-1), *(p+1)) >= angle) points.push_back(*p);
Points out;
for (Point p2 : poly) {
// p2 is next point to the currently visited point p1.
Vec2d v2 = (p2 - p1).cast<double>();
if (filter(v1, v2))
out.emplace_back(p2);
v1 = v2;
p1 = p2;
}
// check whether last point forms a convex angle
if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) >= angle)
points.push_back(this->points.back());
return out;
}
return points;
template<typename ConvexConcaveFilterFn>
Points filter_convex_concave_points_by_angle_threshold(const Points &poly, double angle_threshold, ConvexConcaveFilterFn convex_concave_filter)
{
assert(angle_threshold >= 0.);
if (angle_threshold < EPSILON) {
double cos_angle = cos(angle_threshold);
return filter_points_by_vectors(poly, [convex_concave_filter, cos_angle](const Vec2d &v1, const Vec2d &v2){
return convex_concave_filter(v1, v2) && v1.normalized().dot(v2.normalized()) < cos_angle;
});
} else {
return filter_points_by_vectors(poly, [convex_concave_filter](const Vec2d &v1, const Vec2d &v2){
return convex_concave_filter(v1, v2);
});
}
}
Points Polygon::convex_points(double angle_threshold) const
{
return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) > 0.; });
}
Points Polygon::concave_points(double angle_threshold) const
{
return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) < 0.; });
}
// Projection of a point onto the polygon.

7
src/libslic3r/Polygon.hpp
View File

@ -65,8 +65,11 @@ public:
void densify(float min_length, std::vector<float>* lengths = nullptr);
void triangulate_convex(Polygons* polygons) const;
Point centroid() const;
Points concave_points(double angle = PI) const;
Points convex_points(double angle = PI) const;
// Considering CCW orientation of this polygon, find all convex resp. concave points
// with the angle at the vertex larger than a threshold.
// Zero angle_threshold means to accept all convex resp. concave points.
Points convex_points(double angle_threshold = 0.) const;
Points concave_points(double angle_threshold = 0.) const;
// Projection of a point onto the polygon.
Point point_projection(const Point &point) const;
std::vector<float> parameter_by_length() const;

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

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_

17
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); }

6
src/libslic3r/pchheader.hpp
View File

@ -64,6 +64,12 @@
#include <boost/config.hpp>
#include <boost/config/warning_disable.hpp>
#include <boost/container/small_vector.hpp>
#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 <boost/date_time/local_time/local_time.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/filesystem.hpp>

59
t/geometry.t
View File

@ -2,7 +2,7 @@ use Test::More;
use strict;
use warnings;
plan tests => 42;
plan tests => 30;
BEGIN {
use FindBin;
@ -189,63 +189,6 @@ my $polygons = [
is +Slic3r::Point->new(10, 0)->distance_to_line($line), 0, 'distance_to';
}
#==========================================================
{
my $square = Slic3r::Polygon->new_scale(
[100,100],
[200,100],
[200,200],
[100,200],
);
is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in ccw square';
is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in ccw square';
$square->make_clockwise;
is scalar(@{$square->concave_points(PI*4/3)}), 4, 'fuor concave vertices detected in cw square';
is scalar(@{$square->convex_points(PI*2/3)}), 0, 'no convex vertices detected in cw square';
}
{
my $square = Slic3r::Polygon->new_scale(
[150,100],
[200,100],
[200,200],
[100,200],
[100,100],
);
is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
}
{
my $square = Slic3r::Polygon->new_scale(
[200,200],
[100,200],
[100,100],
[150,100],
[200,100],
);
is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
}
{
my $triangle = Slic3r::Polygon->new(
[16000170,26257364], [714223,461012], [31286371,461008],
);
is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle';
is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle';
}
{
my $triangle = Slic3r::Polygon->new(
[16000170,26257364], [714223,461012], [20000000,461012], [31286371,461012],
);
is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle having collinear point';
is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle having collinear point';
}
{
my $triangle = Slic3r::Polygon->new(
[16000170,26257364], [714223,461012], [31286371,461008],

2
t/perimeters.t
View File

@ -223,7 +223,7 @@ use Slic3r::Test;
if ($model eq 'cube_with_concave_hole') {
# check that loop starts at a concave vertex
my $ccw_angle = $loop->first_point->ccw_angle(@$loop[-2,1]);
my $ccw_angle = $loop->[-2]->ccw($loop->first_point, $loop->[1]);
my $convex = ($ccw_angle > PI); # whether the angle on the *right* side is convex
$starts_on_convex_point = 1
if ($convex && $is_contour) || (!$convex && $is_hole);

1
tests/libslic3r/CMakeLists.txt
View File

@ -8,6 +8,7 @@ add_executable(${_TEST_NAME}_tests
test_clipper_utils.cpp
test_color.cpp
test_config.cpp
test_curve_fitting.cpp
test_elephant_foot_compensation.cpp
test_geometry.cpp
test_placeholder_parser.cpp

118
tests/libslic3r/test_curve_fitting.cpp Normal file
View File

@ -0,0 +1,118 @@
#include <catch2/catch.hpp>
#include <test_utils.hpp>
#include <libslic3r/Geometry/Curves.hpp>
#include <libslic3r/Utils.hpp>
#include <libslic3r/SVG.hpp>
TEST_CASE("Curves: cubic b spline fit test", "[Curves]") {
using namespace Slic3r;
using namespace Slic3r::Geometry;
auto fx = [&](size_t index) {
return float(index) / 200.0f;
};
auto fy = [&](size_t index) {
return 1.0f;
};
std::vector<Vec<1, float>> observations { };
std::vector<float> observation_points { };
std::vector<float> weights { };
for (size_t index = 0; index < 200; ++index) {
observations.push_back(Vec<1, float> { fy(index) });
observation_points.push_back(fx(index));
weights.push_back(1);
}
Vec2f fmin { fx(0), fy(0) };
Vec2f fmax { fx(200), fy(200) };
auto bspline = fit_cubic_bspline(observations, observation_points, weights, 1);
Approx ap(1.0f);
ap.epsilon(0.1f);
for (int p = 0; p < 200; ++p) {
float fitted_val = bspline.get_fitted_value(fx(p))(0);
float expected = fy(p);
REQUIRE(fitted_val == ap(expected));
}
}
TEST_CASE("Curves: quadratic f cubic b spline fit test", "[Curves]") {
using namespace Slic3r;
using namespace Slic3r::Geometry;
auto fx = [&](size_t index) {
return float(index) / 100.0f;
};
auto fy = [&](size_t index) {
return (fx(index) - 1) * (fx(index) - 1);
};
std::vector<Vec<1, float>> observations { };
std::vector<float> observation_points { };
std::vector<float> weights { };
for (size_t index = 0; index < 200; ++index) {
observations.push_back(Vec<1, float> { fy(index) });
observation_points.push_back(fx(index));
weights.push_back(1);
}
Vec2f fmin { fx(0), fy(0) };
Vec2f fmax { fx(200), fy(200) };
auto bspline = fit_cubic_bspline(observations, observation_points, weights, 10);
for (int p = 0; p < 200; ++p) {
float fitted_val = bspline.get_fitted_value(fx(p))(0);
float expected = fy(p);
auto check = [](float a, float b) {
return abs(a - b) < 0.2f;
};
//Note: checking is problematic, splines will not perfectly align
REQUIRE(check(fitted_val, expected));
}
}
TEST_CASE("Curves: polynomial fit test", "[Curves]") {
using namespace Slic3r;
using namespace Slic3r::Geometry;
auto fx = [&](size_t index) {
return float(index) / 100.0f;
};
auto fy = [&](size_t index) {
return (fx(index) - 1) * (fx(index) - 1);
};
std::vector<Vec<1, float>> observations { };
std::vector<float> observation_points { };
std::vector<float> weights { };
for (size_t index = 0; index < 200; ++index) {
observations.push_back(Vec<1, float> { fy(index) });
observation_points.push_back(fx(index));
weights.push_back(1);
}
Vec2f fmin { fx(0), fy(0) };
Vec2f fmax { fx(200), fy(200) };
Approx ap(1.0f);
ap.epsilon(0.1f);
auto poly = fit_polynomial(observations, observation_points, weights, 2);
REQUIRE(poly.coefficients(0, 0) == ap(1));
REQUIRE(poly.coefficients(0, 1) == ap(-2));
REQUIRE(poly.coefficients(0, 2) == ap(1));
}

64
tests/libslic3r/test_geometry.cpp
View File

@ -373,7 +373,43 @@ SCENARIO("Line distances", "[Geometry]"){
}
}
SCENARIO("Calculating angles", "[Geometry]")
{
GIVEN(("Vectors 30 degrees apart"))
{
std::vector<std::pair<Point, Point>> pts {
{ {1000, 0}, { 866, 500 } },
{ { 866, 500 }, { 500, 866 } },
{ { 500, 866 }, { 0, 1000 } },
{ { -500, 866 }, { -866, 500 } }
};
THEN("Angle detected is 30 degrees")
{
for (auto &p : pts)
REQUIRE(is_approx(angle(p.first, p.second), M_PI / 6.));
}
}
GIVEN(("Vectors 30 degrees apart"))
{
std::vector<std::pair<Point, Point>> pts {
{ { 866, 500 }, {1000, 0} },
{ { 500, 866 }, { 866, 500 } },
{ { 0, 1000 }, { 500, 866 } },
{ { -866, 500 }, { -500, 866 } }
};
THEN("Angle detected is -30 degrees")
{
for (auto &p : pts)
REQUIRE(is_approx(angle(p.first, p.second), - M_PI / 6.));
}
}
}
SCENARIO("Polygon convex/concave detection", "[Geometry]"){
static constexpr const double angle_threshold = M_PI / 3.;
GIVEN(("A Square with dimension 100")){
auto square = Slic3r::Polygon /*new_scale*/(std::vector<Point>({
Point(100,100),
@ -381,13 +417,13 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(200,200),
Point(100,200)}));
THEN("It has 4 convex points counterclockwise"){
REQUIRE(square.concave_points(PI*4/3).size() == 0);
REQUIRE(square.convex_points(PI*2/3).size() == 4);
REQUIRE(square.concave_points(angle_threshold).size() == 0);
REQUIRE(square.convex_points(angle_threshold).size() == 4);
}
THEN("It has 4 concave points clockwise"){
square.make_clockwise();
REQUIRE(square.concave_points(PI*4/3).size() == 4);
REQUIRE(square.convex_points(PI*2/3).size() == 0);
REQUIRE(square.concave_points(angle_threshold).size() == 4);
REQUIRE(square.convex_points(angle_threshold).size() == 0);
}
}
GIVEN("A Square with an extra colinearvertex"){
@ -398,8 +434,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(100,200),
Point(100,100)}));
THEN("It has 4 convex points counterclockwise"){
REQUIRE(square.concave_points(PI*4/3).size() == 0);
REQUIRE(square.convex_points(PI*2/3).size() == 4);
REQUIRE(square.concave_points(angle_threshold).size() == 0);
REQUIRE(square.convex_points(angle_threshold).size() == 4);
}
}
GIVEN("A Square with an extra collinear vertex in different order"){
@ -410,8 +446,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(150,100),
Point(200,100)}));
THEN("It has 4 convex points counterclockwise"){
REQUIRE(square.concave_points(PI*4/3).size() == 0);
REQUIRE(square.convex_points(PI*2/3).size() == 4);
REQUIRE(square.concave_points(angle_threshold).size() == 0);
REQUIRE(square.convex_points(angle_threshold).size() == 4);
}
}
@ -422,8 +458,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(31286371,461008)
}));
THEN("it has three convex vertices"){
REQUIRE(triangle.concave_points(PI*4/3).size() == 0);
REQUIRE(triangle.convex_points(PI*2/3).size() == 3);
REQUIRE(triangle.concave_points(angle_threshold).size() == 0);
REQUIRE(triangle.convex_points(angle_threshold).size() == 3);
}
}
@ -435,8 +471,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(31286371,461012)
}));
THEN("it has three convex vertices"){
REQUIRE(triangle.concave_points(PI*4/3).size() == 0);
REQUIRE(triangle.convex_points(PI*2/3).size() == 3);
REQUIRE(triangle.concave_points(angle_threshold).size() == 0);
REQUIRE(triangle.convex_points(angle_threshold).size() == 3);
}
}
GIVEN("A polygon with concave vertices with angles of specifically 4/3pi"){
@ -453,8 +489,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){
Point(38092663,692699),Point(52100125,692699)
}));
THEN("the correct number of points are detected"){
REQUIRE(polygon.concave_points(PI*4/3).size() == 6);
REQUIRE(polygon.convex_points(PI*2/3).size() == 10);
REQUIRE(polygon.concave_points(angle_threshold).size() == 6);
REQUIRE(polygon.convex_points(angle_threshold).size() == 10);
}
}
}

2
xs/xsp/Line.xsp
View File

@ -41,7 +41,7 @@
Clone<Point> normal();
Clone<Point> vector();
double ccw(Point* point)
%code{% RETVAL = THIS->ccw(*point); %};
%code{% RETVAL = cross2((THIS->a - *point).cast<double>(), (THIS->b - THIS->a).cast<double>()); %};
%{
Line*

8
xs/xsp/Point.xsp
View File

@ -39,13 +39,7 @@
double perp_distance_to_line(Line* line)
%code{% RETVAL = line->perp_distance_to(*THIS); %};
double ccw(Point* p1, Point* p2)
%code{% RETVAL = THIS->ccw(*p1, *p2); %};
double ccw_angle(Point* p1, Point* p2)
%code{% RETVAL = THIS->ccw_angle(*p1, *p2); %};
Point* projection_onto_polygon(Polygon* polygon)
%code{% RETVAL = new Point(THIS->projection_onto(*polygon)); %};
Point* projection_onto_polyline(Polyline* polyline)
%code{% RETVAL = new Point(THIS->projection_onto(*polyline)); %};
%code{% RETVAL = cross2((*p1 - *THIS).cast<double>(), (*p2 - *p1).cast<double>()); %};
Point* projection_onto_line(Line* line)
%code{% RETVAL = new Point(THIS->projection_onto(*line)); %};
Point* negative()

2
xs/xsp/Polygon.xsp
View File

@ -39,8 +39,6 @@
%code{% THIS->triangulate_convex(&RETVAL); %};
Clone<Point> centroid();
Clone<BoundingBox> bounding_box();
Points concave_points(double angle);
Points convex_points(double angle);
Clone<Point> point_projection(Point* point)
%code{% RETVAL = THIS->point_projection(*point); %};
Clone<Point> intersection(Line* line)