3D-printing/OrcaSlicer
1
0
Fork 0
mirror of https://github.com/SoftFever/OrcaSlicer.git synced 2025-10-27 18:51:11 -06:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'master' into tm_colldetection_upgr

Browse source
This commit is contained in:
tamasmeszaros 2019-02-25 13:24:43 +01:00
commit 03079381e1
119 changed files with 11942 additions and 8619 deletions
Download patch file Download diff file Expand all files Collapse all files

739
src/libslic3r/SLA/SLAAutoSupports.cpp
View file

@ -1,47 +1,23 @@
#include "igl/random_points_on_mesh.h"
#include "igl/AABB.h"
#include <tbb/parallel_for.h>
#include "SLAAutoSupports.hpp"
#include "Model.hpp"
#include "ExPolygon.hpp"
#include "SVG.hpp"
#include "Point.hpp"
#include "ClipperUtils.hpp"
#include "Tesselate.hpp"
#include "libslic3r.h"
#include <iostream>
#include <random>
namespace Slic3r {
SLAAutoSupports::SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices, const std::vector<float>& heights,
const Config& config, std::function<void(void)> throw_on_cancel)
: m_config(config), m_V(emesh.V()), m_F(emesh.F()), m_throw_on_cancel(throw_on_cancel)
{
// FIXME: It might be safer to get rid of the rand() calls altogether, because it is probably
// not always thread-safe and can be slow if it is.
srand(time(NULL)); // rand() is used by igl::random_point_on_mesh
// Find all separate islands that will need support. The coord_t number denotes height
// of a point just below the mesh (so that we can later project the point precisely
// on the mesh by raycasting (done by igl) and not risking we will place the point inside).
std::vector<std::pair<ExPolygon, coord_t>> islands = find_islands(slices, heights);
// Uniformly cover each of the islands with support points.
for (const auto& island : islands) {
std::vector<Vec3d> points = uniformly_cover(island);
m_throw_on_cancel();
project_upward_onto_mesh(points);
m_output.insert(m_output.end(), points.begin(), points.end());
m_throw_on_cancel();
}
// We are done with the islands. Let's sprinkle the rest of the mesh.
// The function appends to m_output.
sprinkle_mesh(mesh);
}
float SLAAutoSupports::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2)
/*float SLAAutoSupports::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2)
{
n1.normalize();
n2.normalize();
@ -59,115 +35,6 @@ float SLAAutoSupports::approximate_geodesic_distance(const Vec3d& p1, const Vec3
}
void SLAAutoSupports::sprinkle_mesh(const TriangleMesh& mesh)
{
std::vector<Vec3d> points;
// Loads the ModelObject raw_mesh and transforms it by first instance's transformation matrix (disregarding translation).
// Instances only differ in z-rotation, so it does not matter which of them will be used for the calculation.
// The supports point will be calculated on this mesh (so scaling ang vertical direction is correctly accounted for).
// Results will be inverse-transformed to raw_mesh coordinates.
//TriangleMesh mesh = m_model_object.raw_mesh();
//Transform3d transformation_matrix = m_model_object.instances[0]->get_matrix(true/*dont_translate*/);
//mesh.transform(transformation_matrix);
// Check that the object is thick enough to produce any support points
BoundingBoxf3 bb = mesh.bounding_box();
if (bb.size()(2) < m_config.minimal_z)
return;
// All points that we curretly have must be transformed too, so distance to them is correcly calculated.
//for (Vec3f& point : m_model_object.sla_support_points)
// point = transformation_matrix.cast<float>() * point;
// In order to calculate distance to already placed points, we must keep know which facet the point lies on.
std::vector<Vec3d> facets_normals;
// Only points belonging to islands were added so far - they all lie on horizontal surfaces:
for (unsigned int i=0; i<m_output.size(); ++i)
facets_normals.push_back(Vec3d(0,0,-1));
// The AABB hierarchy will be used to find normals of already placed points.
// The points added automatically will just push_back the new normal on the fly.
/*igl::AABB<Eigen::MatrixXf,3> aabb;
aabb.init(V, F);
for (unsigned int i=0; i<m_model_object.sla_support_points.size(); ++i) {
int facet_idx = 0;
Eigen::Matrix<float, 1, 3> dump;
Eigen::MatrixXf query_point = m_model_object.sla_support_points[i];
aabb.squared_distance(V, F, query_point, facet_idx, dump);
Vec3f a1 = V.row(F(facet_idx,1)) - V.row(F(facet_idx,0));
Vec3f a2 = V.row(F(facet_idx,2)) - V.row(F(facet_idx,0));
Vec3f normal = a1.cross(a2);
normal.normalize();
facets_normals.push_back(normal);
}*/
// New potential support point is randomly generated on the mesh and distance to all already placed points is calculated.
// In case it is never smaller than certain limit (depends on the new point's facet normal), the point is accepted.
// The process stops after certain number of points is refused in a row.
Vec3d point;
Vec3d normal;
int added_points = 0;
int refused_points = 0;
const int refused_limit = 30;
// Angle at which the density reaches zero:
const float threshold_angle = std::min(M_PI_2, M_PI_4 * acos(0.f/m_config.density_at_horizontal) / acos(m_config.density_at_45/m_config.density_at_horizontal));
size_t cancel_test_cntr = 0;
while (refused_points < refused_limit) {
if (++ cancel_test_cntr == 500) {
// Don't call the cancellation routine too often as the multi-core cache synchronization
// may be pretty expensive.
m_throw_on_cancel();
cancel_test_cntr = 0;
}
// Place a random point on the mesh and calculate corresponding facet's normal:
Eigen::VectorXi FI;
Eigen::MatrixXd B;
igl::random_points_on_mesh(1, m_V, m_F, B, FI);
point = B(0,0)*m_V.row(m_F(FI(0),0)) +
B(0,1)*m_V.row(m_F(FI(0),1)) +
B(0,2)*m_V.row(m_F(FI(0),2));
if (point(2) - bb.min(2) < m_config.minimal_z)
continue;
Vec3d a1 = m_V.row(m_F(FI(0),1)) - m_V.row(m_F(FI(0),0));
Vec3d a2 = m_V.row(m_F(FI(0),2)) - m_V.row(m_F(FI(0),0));
normal = a1.cross(a2);
normal.normalize();
// calculate angle between the normal and vertical:
float angle = angle_from_normal(normal.cast<float>());
if (angle > threshold_angle)
continue;
const float limit = distance_limit(angle);
bool add_it = true;
for (unsigned int i=0; i<points.size(); ++i) {
if (approximate_geodesic_distance(points[i], point, facets_normals[i], normal) < limit) {
add_it = false;
++refused_points;
break;
}
}
if (add_it) {
points.push_back(point.cast<double>());
facets_normals.push_back(normal);
++added_points;
refused_points = 0;
}
}
m_output.insert(m_output.end(), points.begin(), points.end());
// Now transform all support points to mesh coordinates:
//for (Vec3f& point : m_model_object.sla_support_points)
// point = transformation_matrix.inverse().cast<float>() * point;
}
float SLAAutoSupports::get_required_density(float angle) const
{
// calculation would be density_0 * cos(angle). To provide one more degree of freedom, we will scale the angle
@ -179,10 +46,470 @@ float SLAAutoSupports::get_required_density(float angle) const
float SLAAutoSupports::distance_limit(float angle) const
{
return 1./(2.4*get_required_density(angle));
}*/
SLAAutoSupports::SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices, const std::vector<float>& heights,
const Config& config, std::function<void(void)> throw_on_cancel)
: m_config(config), m_emesh(emesh), m_throw_on_cancel(throw_on_cancel)
{
process(slices, heights);
project_onto_mesh(m_output);
}
void SLAAutoSupports::project_onto_mesh(std::vector<sla::SupportPoint>& points) const
{
// The function makes sure that all the points are really exactly placed on the mesh.
igl::Hit hit_up{0, 0, 0.f, 0.f, 0.f};
igl::Hit hit_down{0, 0, 0.f, 0.f, 0.f};
// Use a reasonable granularity to account for the worker thread synchronization cost.
tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size(), 64),
[this, &points](const tbb::blocked_range<size_t>& range) {
for (size_t point_id = range.begin(); point_id < range.end(); ++ point_id) {
if ((point_id % 16) == 0)
// Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
m_throw_on_cancel();
Vec3f& p = points[point_id].pos;
// Project the point upward and downward and choose the closer intersection with the mesh.
//bool up = igl::ray_mesh_intersect(p.cast<float>(), Vec3f(0., 0., 1.), m_V, m_F, hit_up);
//bool down = igl::ray_mesh_intersect(p.cast<float>(), Vec3f(0., 0., -1.), m_V, m_F, hit_down);
sla::EigenMesh3D::hit_result hit_up = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., 1.));
sla::EigenMesh3D::hit_result hit_down = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., -1.));
bool up = hit_up.face() != -1;
bool down = hit_down.face() != -1;
if (!up && !down)
continue;
sla::EigenMesh3D::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down;
//int fid = hit.face();
//Vec3f bc(1-hit.u-hit.v, hit.u, hit.v);
//p = (bc(0) * m_V.row(m_F(fid, 0)) + bc(1) * m_V.row(m_F(fid, 1)) + bc(2)*m_V.row(m_F(fid, 2))).cast<float>();
p = p + (hit.distance() * hit.direction()).cast<float>();
}
});
}
static std::vector<SLAAutoSupports::MyLayer> make_layers(
const std::vector<ExPolygons>& slices, const std::vector<float>& heights,
std::function<void(void)> throw_on_cancel)
{
assert(slices.size() == heights.size());
// Allocate empty layers.
std::vector<SLAAutoSupports::MyLayer> layers;
layers.reserve(slices.size());
for (size_t i = 0; i < slices.size(); ++ i)
layers.emplace_back(i, heights[i]);
// FIXME: calculate actual pixel area from printer config:
//const float pixel_area = pow(wxGetApp().preset_bundle->project_config.option<ConfigOptionFloat>("display_width") / wxGetApp().preset_bundle->project_config.option<ConfigOptionInt>("display_pixels_x"), 2.f); //
const float pixel_area = pow(0.047f, 2.f);
// Use a reasonable granularity to account for the worker thread synchronization cost.
tbb::parallel_for(tbb::blocked_range<size_t>(0, layers.size(), 32),
[&layers, &slices, &heights, pixel_area, throw_on_cancel](const tbb::blocked_range<size_t>& range) {
for (size_t layer_id = range.begin(); layer_id < range.end(); ++ layer_id) {
if ((layer_id % 8) == 0)
// Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
throw_on_cancel();
SLAAutoSupports::MyLayer &layer = layers[layer_id];
const ExPolygons &islands = slices[layer_id];
//FIXME WTF?
const float height = (layer_id>2 ? heights[layer_id-3] : heights[0]-(heights[1]-heights[0]));
layer.islands.reserve(islands.size());
for (const ExPolygon &island : islands) {
float area = float(island.area() * SCALING_FACTOR * SCALING_FACTOR);
if (area >= pixel_area)
//FIXME this is not a correct centroid of a polygon with holes.
layer.islands.emplace_back(layer, island, get_extents(island.contour), Slic3r::unscale(island.contour.centroid()).cast<float>(), area, height);
}
}
});
// Calculate overlap of successive layers. Link overlapping islands.
tbb::parallel_for(tbb::blocked_range<size_t>(1, layers.size(), 8),
[&layers, &heights, throw_on_cancel](const tbb::blocked_range<size_t>& range) {
for (size_t layer_id = range.begin(); layer_id < range.end(); ++layer_id) {
if ((layer_id % 2) == 0)
// Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
throw_on_cancel();
SLAAutoSupports::MyLayer &layer_above = layers[layer_id];
SLAAutoSupports::MyLayer &layer_below = layers[layer_id - 1];
//FIXME WTF?
const float height = (layer_id>2 ? heights[layer_id-3] : heights[0]-(heights[1]-heights[0]));
const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]);
const float safe_angle = 5.f * (float(M_PI)/180.f); // smaller number - less supports
const float between_layers_offset = float(scale_(layer_height / std::tan(safe_angle)));
//FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands.
for (SLAAutoSupports::Structure &top : layer_above.islands) {
for (SLAAutoSupports::Structure &bottom : layer_below.islands) {
float overlap_area = top.overlap_area(bottom);
if (overlap_area > 0) {
top.islands_below.emplace_back(&bottom, overlap_area);
bottom.islands_above.emplace_back(&top, overlap_area);
}
}
if (! top.islands_below.empty()) {
Polygons top_polygons = to_polygons(*top.polygon);
Polygons bottom_polygons = top.polygons_below();
top.overhangs = diff_ex(top_polygons, bottom_polygons);
if (! top.overhangs.empty()) {
top.overhangs_area = 0.f;
std::vector<std::pair<ExPolygon*, float>> expolys_with_areas;
for (ExPolygon &ex : top.overhangs) {
float area = float(ex.area());
expolys_with_areas.emplace_back(&ex, area);
top.overhangs_area += area;
}
std::sort(expolys_with_areas.begin(), expolys_with_areas.end(),
[](const std::pair<ExPolygon*, float> &p1, const std::pair<ExPolygon*, float> &p2)
{ return p1.second > p2.second; });
ExPolygons overhangs_sorted;
for (auto &p : expolys_with_areas)
overhangs_sorted.emplace_back(std::move(*p.first));
top.overhangs = std::move(overhangs_sorted);
top.overhangs_area *= float(SCALING_FACTOR * SCALING_FACTOR);
top.dangling_areas = diff_ex(top_polygons, offset(bottom_polygons, between_layers_offset));
}
}
}
}
});
return layers;
}
void SLAAutoSupports::process(const std::vector<ExPolygons>& slices, const std::vector<float>& heights)
{
#ifdef SLA_AUTOSUPPORTS_DEBUG
std::vector<std::pair<ExPolygon, coord_t>> islands;
#endif /* SLA_AUTOSUPPORTS_DEBUG */
std::vector<SLAAutoSupports::MyLayer> layers = make_layers(slices, heights, m_throw_on_cancel);
PointGrid3D point_grid;
point_grid.cell_size = Vec3f(10.f, 10.f, 10.f);
for (unsigned int layer_id = 0; layer_id < layers.size(); ++ layer_id) {
SLAAutoSupports::MyLayer *layer_top = &layers[layer_id];
SLAAutoSupports::MyLayer *layer_bottom = (layer_id > 0) ? &layers[layer_id - 1] : nullptr;
std::vector<float> support_force_bottom;
if (layer_bottom != nullptr) {
support_force_bottom.assign(layer_bottom->islands.size(), 0.f);
for (size_t i = 0; i < layer_bottom->islands.size(); ++ i)
support_force_bottom[i] = layer_bottom->islands[i].supports_force_total();
}
for (Structure &top : layer_top->islands)
for (Structure::Link &bottom_link : top.islands_below) {
Structure &bottom = *bottom_link.island;
float centroids_dist = (bottom.centroid - top.centroid).norm();
// Penalization resulting from centroid offset:
// bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, (1600.f * layer_height) * centroids_dist * centroids_dist / bottom.area));
float &support_force = support_force_bottom[&bottom - layer_bottom->islands.data()];
//FIXME this condition does not reflect a bifurcation into a one large island and one tiny island well, it incorrectly resets the support force to zero.
// One should rather work with the overlap area vs overhang area.
// support_force *= std::min(1.f, 1.f - std::min(1.f, 0.1f * centroids_dist * centroids_dist / bottom.area));
// Penalization resulting from increasing polygon area:
support_force *= std::min(1.f, 20.f * bottom.area / top.area);
}
// Let's assign proper support force to each of them:
if (layer_id > 0) {
for (Structure &below : layer_bottom->islands) {
float below_support_force = support_force_bottom[&below - layer_bottom->islands.data()];
float above_overlap_area = 0.f;
for (Structure::Link &above_link : below.islands_above)
above_overlap_area += above_link.overlap_area;
for (Structure::Link &above_link : below.islands_above)
above_link.island->supports_force_inherited += below_support_force * above_link.overlap_area / above_overlap_area;
}
}
// Now iterate over all polygons and append new points if needed.
for (Structure &s : layer_top->islands) {
// Penalization resulting from large diff from the last layer:
// s.supports_force_inherited /= std::max(1.f, (layer_height / 0.3f) * e_area / s.area);
s.supports_force_inherited /= std::max(1.f, 0.17f * (s.overhangs_area) / s.area);
float force_deficit = s.support_force_deficit(m_config.tear_pressure());
if (s.islands_below.empty()) { // completely new island - needs support no doubt
uniformly_cover({ *s.polygon }, s, point_grid, true);
} else if (! s.dangling_areas.empty()) {
// Let's see if there's anything that overlaps enough to need supports:
// What we now have in polygons needs support, regardless of what the forces are, so we can add them.
//FIXME is it an island point or not? Vojtech thinks it is.
uniformly_cover(s.dangling_areas, s, point_grid);
} else if (! s.overhangs.empty()) {
//FIXME add the support force deficit as a parameter, only cover until the defficiency is covered.
uniformly_cover(s.overhangs, s, point_grid);
}
}
m_throw_on_cancel();
#ifdef SLA_AUTOSUPPORTS_DEBUG
/*std::string layer_num_str = std::string((i<10 ? "0" : "")) + std::string((i<100 ? "0" : "")) + std::to_string(i);
output_expolygons(expolys_top, "top" + layer_num_str + ".svg");
output_expolygons(diff, "diff" + layer_num_str + ".svg");
if (!islands.empty())
output_expolygons(islands, "islands" + layer_num_str + ".svg");*/
#endif /* SLA_AUTOSUPPORTS_DEBUG */
}
}
std::vector<Vec2f> sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng)
{
// Triangulate the polygon with holes into triplets of 3D points.
std::vector<Vec2f> triangles = Slic3r::triangulate_expolygon_2f(expoly);
std::vector<Vec2f> out;
if (! triangles.empty())
{
// Calculate area of each triangle.
std::vector<float> areas;
areas.reserve(triangles.size() / 3);
for (size_t i = 0; i < triangles.size(); ) {
const Vec2f &a = triangles[i ++];
const Vec2f v1 = triangles[i ++] - a;
const Vec2f v2 = triangles[i ++] - a;
areas.emplace_back(0.5f * std::abs(cross2(v1, v2)));
if (i != 3)
// Prefix sum of the areas.
areas.back() += areas[areas.size() - 2];
}
size_t num_samples = size_t(ceil(areas.back() * samples_per_mm2));
std::uniform_real_distribution<> random_triangle(0., double(areas.back()));
std::uniform_real_distribution<> random_float(0., 1.);
for (size_t i = 0; i < num_samples; ++ i) {
double r = random_triangle(rng);
size_t idx_triangle = std::min<size_t>(std::upper_bound(areas.begin(), areas.end(), (float)r) - areas.begin(), areas.size() - 1) * 3;
// Select a random point on the triangle.
double u = float(sqrt(random_float(rng)));
double v = float(random_float(rng));
const Vec2f &a = triangles[idx_triangle ++];
const Vec2f &b = triangles[idx_triangle++];
const Vec2f &c = triangles[idx_triangle];
const Vec2f x = a * (1.f - u) + b * (u * (1.f - v)) + c * (v * u);
out.emplace_back(x);
}
}
return out;
}
std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygon &expoly, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
{
std::vector<Vec2f> out = sample_expolygon(expoly, samples_per_mm2, rng);
double point_stepping_scaled = scale_(1.f) / samples_per_mm_boundary;
for (size_t i_contour = 0; i_contour <= expoly.holes.size(); ++ i_contour) {
const Polygon &contour = (i_contour == 0) ? expoly.contour : expoly.holes[i_contour - 1];
const Points pts = contour.equally_spaced_points(point_stepping_scaled);
for (size_t i = 0; i < pts.size(); ++ i)
out.emplace_back(unscale<float>(pts[i].x()), unscale<float>(pts[i].y()));
}
return out;
}
std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygons &expolys, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
{
std::vector<Vec2f> out;
for (const ExPolygon &expoly : expolys)
append(out, sample_expolygon_with_boundary(expoly, samples_per_mm2, samples_per_mm_boundary, rng));
return out;
}
template<typename REFUSE_FUNCTION>
static inline std::vector<Vec2f> poisson_disk_from_samples(const std::vector<Vec2f> &raw_samples, float radius, REFUSE_FUNCTION refuse_function)
{
Vec2f corner_min(std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
for (const Vec2f &pt : raw_samples) {
corner_min.x() = std::min(corner_min.x(), pt.x());
corner_min.y() = std::min(corner_min.y(), pt.y());
}
// Assign the raw samples to grid cells, sort the grid cells lexicographically.
struct RawSample {
Vec2f coord;
Vec2i cell_id;
};
std::vector<RawSample> raw_samples_sorted;
RawSample sample;
for (const Vec2f &pt : raw_samples) {
sample.coord = pt;
sample.cell_id = ((pt - corner_min) / radius).cast<int>();
raw_samples_sorted.emplace_back(sample);
}
std::sort(raw_samples_sorted.begin(), raw_samples_sorted.end(), [](const RawSample &lhs, const RawSample &rhs)
{ return lhs.cell_id.x() < rhs.cell_id.x() || (lhs.cell_id.x() == rhs.cell_id.x() && lhs.cell_id.y() < rhs.cell_id.y()); });
struct PoissonDiskGridEntry {
// Resulting output sample points for this cell:
enum {
max_positions = 4
};
Vec2f poisson_samples[max_positions];
int num_poisson_samples = 0;
// Index into raw_samples:
int first_sample_idx;
int sample_cnt;
};
struct CellIDHash {
std::size_t operator()(const Vec2i &cell_id) const {
return std::hash<int>()(cell_id.x()) ^ std::hash<int>()(cell_id.y() * 593);
}
};
// Map from cell IDs to hash_data. Each hash_data points to the range in raw_samples corresponding to that cell.
// (We could just store the samples in hash_data. This implementation is an artifact of the reference paper, which
// is optimizing for GPU acceleration that we haven't implemented currently.)
typedef std::unordered_map<Vec2i, PoissonDiskGridEntry, CellIDHash> Cells;
Cells cells;
{
typename Cells::iterator last_cell_id_it;
Vec2i last_cell_id(-1, -1);
for (int i = 0; i < raw_samples_sorted.size(); ++ i) {
const RawSample &sample = raw_samples_sorted[i];
if (sample.cell_id == last_cell_id) {
// This sample is in the same cell as the previous, so just increase the count. Cells are
// always contiguous, since we've sorted raw_samples_sorted by cell ID.
++ last_cell_id_it->second.sample_cnt;
} else {
// This is a new cell.
PoissonDiskGridEntry data;
data.first_sample_idx = i;
data.sample_cnt = 1;
auto result = cells.insert({sample.cell_id, data});
last_cell_id = sample.cell_id;
last_cell_id_it = result.first;
}
}
}
const int max_trials = 5;
const float radius_squared = radius * radius;
for (int trial = 0; trial < max_trials; ++ trial) {
// Create sample points for each entry in cells.
for (auto &it : cells) {
const Vec2i &cell_id = it.first;
PoissonDiskGridEntry &cell_data = it.second;
// This cell's raw sample points start at first_sample_idx. On trial 0, try the first one. On trial 1, try first_sample_idx + 1.
int next_sample_idx = cell_data.first_sample_idx + trial;
if (trial >= cell_data.sample_cnt)
// There are no more points to try for this cell.
continue;
const RawSample &candidate = raw_samples_sorted[next_sample_idx];
// See if this point conflicts with any other points in this cell, or with any points in
// neighboring cells. Note that it's possible to have more than one point in the same cell.
bool conflict = refuse_function(candidate.coord);
for (int i = -1; i < 2 && ! conflict; ++ i) {
for (int j = -1; j < 2; ++ j) {
const auto &it_neighbor = cells.find(cell_id + Vec2i(i, j));
if (it_neighbor != cells.end()) {
const PoissonDiskGridEntry &neighbor = it_neighbor->second;
for (int i_sample = 0; i_sample < neighbor.num_poisson_samples; ++ i_sample)
if ((neighbor.poisson_samples[i_sample] - candidate.coord).squaredNorm() < radius_squared) {
conflict = true;
break;
}
}
}
}
if (! conflict) {
// Store the new sample.
assert(cell_data.num_poisson_samples < cell_data.max_positions);
if (cell_data.num_poisson_samples < cell_data.max_positions)
cell_data.poisson_samples[cell_data.num_poisson_samples ++] = candidate.coord;
}
}
}
// Copy the results to the output.
std::vector<Vec2f> out;
for (const auto& it : cells)
for (int i = 0; i < it.second.num_poisson_samples; ++ i)
out.emplace_back(it.second.poisson_samples[i]);
return out;
}
void SLAAutoSupports::uniformly_cover(const ExPolygons& islands, Structure& structure, PointGrid3D &grid3d, bool is_new_island, bool just_one)
{
//int num_of_points = std::max(1, (int)((island.area()*pow(SCALING_FACTOR, 2) * m_config.tear_pressure)/m_config.support_force));
const float support_force_deficit = structure.support_force_deficit(m_config.tear_pressure());
if (support_force_deficit < 0)
return;
// Number of newly added points.
const size_t poisson_samples_target = size_t(ceil(support_force_deficit / m_config.support_force()));
const float density_horizontal = m_config.tear_pressure() / m_config.support_force();
//FIXME why?
float poisson_radius = std::max(m_config.minimal_distance, 1.f / (5.f * density_horizontal));
// const float poisson_radius = 1.f / (15.f * density_horizontal);
const float samples_per_mm2 = 30.f / (float(M_PI) * poisson_radius * poisson_radius);
// Minimum distance between samples, in 3D space.
// float min_spacing = poisson_radius / 3.f;
float min_spacing = poisson_radius;
//FIXME share the random generator. The random generator may be not so cheap to initialize, also we don't want the random generator to be restarted for each polygon.
std::random_device rd;
std::mt19937 rng(rd());
std::vector<Vec2f> raw_samples = sample_expolygon_with_boundary(islands, samples_per_mm2, 5.f / poisson_radius, rng);
std::vector<Vec2f> poisson_samples;
for (size_t iter = 0; iter < 4; ++ iter) {
poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius,
[&structure, &grid3d, min_spacing](const Vec2f &pos) {
return grid3d.collides_with(pos, &structure, min_spacing);
});
if (poisson_samples.size() >= poisson_samples_target || m_config.minimal_distance > poisson_radius-EPSILON)
break;
float coeff = 0.5f;
if (poisson_samples.size() * 2 > poisson_samples_target)
coeff = float(poisson_samples.size()) / float(poisson_samples_target);
poisson_radius = std::max(m_config.minimal_distance, poisson_radius * coeff);
min_spacing = std::max(m_config.minimal_distance, min_spacing * coeff);
}
#ifdef SLA_AUTOSUPPORTS_DEBUG
{
static int irun = 0;
Slic3r::SVG svg(debug_out_path("SLA_supports-uniformly_cover-%d.svg", irun ++), get_extents(islands));
for (const ExPolygon &island : islands)
svg.draw(island);
for (const Vec2f &pt : raw_samples)
svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "red");
for (const Vec2f &pt : poisson_samples)
svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "blue");
}
#endif /* NDEBUG */
// assert(! poisson_samples.empty());
if (poisson_samples_target < poisson_samples.size()) {
std::shuffle(poisson_samples.begin(), poisson_samples.end(), rng);
poisson_samples.erase(poisson_samples.begin() + poisson_samples_target, poisson_samples.end());
}
for (const Vec2f &pt : poisson_samples) {
m_output.emplace_back(float(pt(0)), float(pt(1)), structure.height, 0.2f, is_new_island);
structure.supports_force_this_layer += m_config.support_force();
grid3d.insert(pt, &structure);
}
}
#ifdef SLA_AUTOSUPPORTS_DEBUG
void SLAAutoSupports::output_expolygons(const ExPolygons& expolys, std::string filename) const
void SLAAutoSupports::output_structures(const std::vector<Structure>& structures)
{
for (unsigned int i=0 ; i<structures.size(); ++i) {
std::stringstream ss;
ss << structures[i].unique_id.count() << "_" << std::setw(10) << std::setfill('0') << 1000 + (int)structures[i].height/1000 << ".png";
output_expolygons(std::vector<ExPolygon>{*structures[i].polygon}, ss.str());
}
}
void SLAAutoSupports::output_expolygons(const ExPolygons& expolys, const std::string &filename)
{
BoundingBox bb(Point(-30000000, -30000000), Point(30000000, 30000000));
Slic3r::SVG svg_cummulative(filename, bb);
@ -198,138 +525,6 @@ void SLAAutoSupports::output_expolygons(const ExPolygons& expolys, std::string f
svg_cummulative.draw_outline(expolys[i].holes, "blue", scale_(0.05));
}
}
#endif /* SLA_AUTOSUPPORTS_DEBUG */
std::vector<std::pair<ExPolygon, coord_t>> SLAAutoSupports::find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights) const
{
std::vector<std::pair<ExPolygon, coord_t>> islands;
struct PointAccessor {
const Point* operator()(const Point &pt) const { return &pt; }
};
typedef ClosestPointInRadiusLookup<Point, PointAccessor> ClosestPointLookupType;
for (unsigned int i = 0; i<slices.size(); ++i) {
const ExPolygons& expolys_top = slices[i];
const ExPolygons& expolys_bottom = (i == 0 ? ExPolygons() : slices[i-1]);
std::string layer_num_str = std::string((i<10 ? "0" : "")) + std::string((i<100 ? "0" : "")) + std::to_string(i);
#ifdef SLA_AUTOSUPPORTS_DEBUG
output_expolygons(expolys_top, "top" + layer_num_str + ".svg");
#endif /* SLA_AUTOSUPPORTS_DEBUG */
ExPolygons diff = diff_ex(expolys_top, expolys_bottom);
#ifdef SLA_AUTOSUPPORTS_DEBUG
output_expolygons(diff, "diff" + layer_num_str + ".svg");
#endif /* SLA_AUTOSUPPORTS_DEBUG */
ClosestPointLookupType cpl(SCALED_EPSILON);
for (const ExPolygon& expol : expolys_top) {
for (const Point& p : expol.contour.points)
cpl.insert(p);
for (const Polygon& hole : expol.holes)
for (const Point& p : hole.points)
cpl.insert(p);
// the lookup structure now contains all points from the top slice
}
for (const ExPolygon& polygon : diff) {
// we want to check all boundary points of the diff polygon
bool island = true;
for (const Point& p : polygon.contour.points) {
if (cpl.find(p).second != 0) { // the point belongs to the bottom slice - this cannot be an island
island = false;
goto NO_ISLAND;
}
}
for (const Polygon& hole : polygon.holes)
for (const Point& p : hole.points)
if (cpl.find(p).second != 0) {
island = false;
goto NO_ISLAND;
}
if (island) { // all points of the diff polygon are from the top slice
islands.push_back(std::make_pair(polygon, scale_(i!=0 ? heights[i-1] : heights[0]-(heights[1]-heights[0]))));
}
NO_ISLAND: ;// continue with next ExPolygon
}
#ifdef SLA_AUTOSUPPORTS_DEBUG
//if (!islands.empty())
// output_expolygons(islands, "islands" + layer_num_str + ".svg");
#endif /* SLA_AUTOSUPPORTS_DEBUG */
m_throw_on_cancel();
}
return islands;
}
std::vector<Vec3d> SLAAutoSupports::uniformly_cover(const std::pair<ExPolygon, coord_t>& island)
{
int num_of_points = std::max(1, (int)(island.first.area()*pow(SCALING_FACTOR, 2) * get_required_density(0)));
// In case there is just one point to place, we'll place it into the polygon's centroid (unless it lies in a hole).
if (num_of_points == 1) {
Point out(island.first.contour.centroid());
for (const auto& hole : island.first.holes)
if (hole.contains(out))
goto HOLE_HIT;
return std::vector<Vec3d>{unscale(out(0), out(1), island.second)};
}
HOLE_HIT:
// In this case either the centroid lies in a hole, or there are multiple points
// to place. We will cover the island another way.
// For now we'll just place the points randomly not too close to the others.
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0., 1.);
std::vector<Vec3d> island_new_points;
const BoundingBox& bb = get_extents(island.first);
const int refused_limit = 30;
int refused_points = 0;
while (refused_points < refused_limit) {
Point out(bb.min(0) + bb.size()(0) * dis(gen),
bb.min(1) + bb.size()(1) * dis(gen)) ;
Vec3d unscaled_out = unscale(out(0), out(1), island.second);
bool add_it = true;
if (!island.first.contour.contains(out))
add_it = false;
else
for (const Polygon& hole : island.first.holes)
if (hole.contains(out))
add_it = false;
if (add_it) {
for (const Vec3d& p : island_new_points) {
if ((p - unscaled_out).squaredNorm() < distance_limit(0)) {
add_it = false;
++refused_points;
break;
}
}
}
if (add_it)
island_new_points.emplace_back(unscaled_out);
}
return island_new_points;
}
void SLAAutoSupports::project_upward_onto_mesh(std::vector<Vec3d>& points) const
{
Vec3f dir(0., 0., 1.);
igl::Hit hit{0, 0, 0.f, 0.f, 0.f};
for (Vec3d& p : points) {
igl::ray_mesh_intersect(p.cast<float>(), dir, m_V, m_F, hit);
int fid = hit.id;
Vec3f bc(1-hit.u-hit.v, hit.u, hit.v);
p = (bc(0) * m_V.row(m_F(fid, 0)) + bc(1) * m_V.row(m_F(fid, 1)) + bc(2)*m_V.row(m_F(fid, 2))).cast<double>();
}
}
#endif
} // namespace Slic3r

195
src/libslic3r/SLA/SLAAutoSupports.hpp
View file

@ -1,9 +1,12 @@
#ifndef SLAAUTOSUPPORTS_HPP_
#define SLAAUTOSUPPORTS_HPP_
#include <libslic3r/ClipperUtils.hpp>
#include <libslic3r/Point.hpp>
#include <libslic3r/TriangleMesh.hpp>
#include <libslic3r/SLA/SLASupportTree.hpp>
#include <libslic3r/SLA/SLACommon.hpp>
#include <boost/container/small_vector.hpp>
// #define SLA_AUTOSUPPORTS_DEBUG
@ -12,36 +15,184 @@ namespace Slic3r {
class SLAAutoSupports {
public:
struct Config {
float density_at_horizontal;
float density_at_45;
float minimal_z;
float density_relative;
float minimal_distance;
///////////////
inline float support_force() const { return 10.f / density_relative; } // a force one point can support (arbitrary force unit)
inline float tear_pressure() const { return 1.f; } // pressure that the display exerts (the force unit per mm2)
};
SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices,
const std::vector<float>& heights, const Config& config, std::function<void(void)> throw_on_cancel);
const std::vector<Vec3d>& output() { return m_output; }
SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3D& emesh, const std::vector<ExPolygons>& slices,
const std::vector<float>& heights, const Config& config, std::function<void(void)> throw_on_cancel);
const std::vector<sla::SupportPoint>& output() { return m_output; }
private:
std::vector<Vec3d> m_output;
std::vector<Vec3d> m_normals;
TriangleMesh mesh;
static float angle_from_normal(const stl_normal& normal) { return acos((-normal.normalized())(2)); }
float get_required_density(float angle) const;
float distance_limit(float angle) const;
static float approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2);
std::vector<std::pair<ExPolygon, coord_t>> find_islands(const std::vector<ExPolygons>& slices, const std::vector<float>& heights) const;
void sprinkle_mesh(const TriangleMesh& mesh);
std::vector<Vec3d> uniformly_cover(const std::pair<ExPolygon, coord_t>& island);
void project_upward_onto_mesh(std::vector<Vec3d>& points) const;
struct MyLayer;
struct Structure {
Structure(MyLayer &layer, const ExPolygon& poly, const BoundingBox &bbox, const Vec2f &centroid, float area, float h) :
layer(&layer), polygon(&poly), bbox(bbox), centroid(centroid), area(area), height(h)
#ifdef SLA_AUTOSUPPORTS_DEBUG
void output_expolygons(const ExPolygons& expolys, std::string filename) const;
, unique_id(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()))
#endif /* SLA_AUTOSUPPORTS_DEBUG */
{}
MyLayer *layer;
const ExPolygon* polygon = nullptr;
const BoundingBox bbox;
const Vec2f centroid = Vec2f::Zero();
const float area = 0.f;
float height = 0;
// How well is this ExPolygon held to the print base?
// Positive number, the higher the better.
float supports_force_this_layer = 0.f;
float supports_force_inherited = 0.f;
float supports_force_total() const { return this->supports_force_this_layer + this->supports_force_inherited; }
#ifdef SLA_AUTOSUPPORTS_DEBUG
std::chrono::milliseconds unique_id;
#endif /* SLA_AUTOSUPPORTS_DEBUG */
struct Link {
Link(Structure *island, float overlap_area) : island(island), overlap_area(overlap_area) {}
Structure *island;
float overlap_area;
};
#ifdef NDEBUG
// In release mode, use the optimized container.
boost::container::small_vector<Link, 4> islands_above;
boost::container::small_vector<Link, 4> islands_below;
#else
// In debug mode, use the standard vector, which is well handled by debugger visualizer.
std::vector<Link> islands_above;
std::vector<Link> islands_below;
#endif
ExPolygons dangling_areas;
ExPolygons overhangs;
float overhangs_area;
bool overlaps(const Structure &rhs) const {
return this->bbox.overlap(rhs.bbox) && (this->polygon->overlaps(*rhs.polygon) || rhs.polygon->overlaps(*this->polygon));
}
float overlap_area(const Structure &rhs) const {
double out = 0.;
if (this->bbox.overlap(rhs.bbox)) {
Polygons polys = intersection(to_polygons(*this->polygon), to_polygons(*rhs.polygon), false);
for (const Polygon &poly : polys)
out += poly.area();
}
return float(out);
}
float area_below() const {
float area = 0.f;
for (const Link &below : this->islands_below)
area += below.island->area;
return area;
}
Polygons polygons_below() const {
size_t cnt = 0;
for (const Link &below : this->islands_below)
cnt += 1 + below.island->polygon->holes.size();
Polygons out;
out.reserve(cnt);
for (const Link &below : this->islands_below) {
out.emplace_back(below.island->polygon->contour);
append(out, below.island->polygon->holes);
}
return out;
}
ExPolygons expolygons_below() const {
ExPolygons out;
out.reserve(this->islands_below.size());
for (const Link &below : this->islands_below)
out.emplace_back(*below.island->polygon);
return out;
}
// Positive deficit of the supports. If negative, this area is well supported. If positive, more supports need to be added.
float support_force_deficit(const float tear_pressure) const { return this->area * tear_pressure - this->supports_force_total(); }
};
struct MyLayer {
MyLayer(const size_t layer_id, coordf_t print_z) : layer_id(layer_id), print_z(print_z) {}
size_t layer_id;
coordf_t print_z;
std::vector<Structure> islands;
};
struct RichSupportPoint {
Vec3f position;
Structure *island;
};
struct PointGrid3D {
struct GridHash {
std::size_t operator()(const Vec3i &cell_id) const {
return std::hash<int>()(cell_id.x()) ^ std::hash<int>()(cell_id.y() * 593) ^ std::hash<int>()(cell_id.z() * 7919);
}
};
typedef std::unordered_multimap<Vec3i, RichSupportPoint, GridHash> Grid;
Vec3f cell_size;
Grid grid;
Vec3i cell_id(const Vec3f &pos) {
return Vec3i(int(floor(pos.x() / cell_size.x())),
int(floor(pos.y() / cell_size.y())),
int(floor(pos.z() / cell_size.z())));
}
void insert(const Vec2f &pos, Structure *island) {
RichSupportPoint pt;
pt.position = Vec3f(pos.x(), pos.y(), float(island->layer->print_z));
pt.island = island;
grid.emplace(cell_id(pt.position), pt);
}
bool collides_with(const Vec2f &pos, Structure *island, float radius) {
Vec3f pos3d(pos.x(), pos.y(), float(island->layer->print_z));
Vec3i cell = cell_id(pos3d);
std::pair<Grid::const_iterator, Grid::const_iterator> it_pair = grid.equal_range(cell);
if (collides_with(pos3d, radius, it_pair.first, it_pair.second))
return true;
for (int i = -1; i < 2; ++ i)
for (int j = -1; j < 2; ++ j)
for (int k = -1; k < 1; ++ k) {
if (i == 0 && j == 0 && k == 0)
continue;
it_pair = grid.equal_range(cell + Vec3i(i, j, k));
if (collides_with(pos3d, radius, it_pair.first, it_pair.second))
return true;
}
return false;
}
private:
bool collides_with(const Vec3f &pos, float radius, Grid::const_iterator it_begin, Grid::const_iterator it_end) {
for (Grid::const_iterator it = it_begin; it != it_end; ++ it) {
float dist2 = (it->second.position - pos).squaredNorm();
if (dist2 < radius * radius)
return true;
}
return false;
}
};
private:
std::vector<sla::SupportPoint> m_output;
SLAAutoSupports::Config m_config;
float m_supports_force_total = 0.f;
void process(const std::vector<ExPolygons>& slices, const std::vector<float>& heights);
void uniformly_cover(const ExPolygons& islands, Structure& structure, PointGrid3D &grid3d, bool is_new_island = false, bool just_one = false);
void project_onto_mesh(std::vector<sla::SupportPoint>& points) const;
#ifdef SLA_AUTOSUPPORTS_DEBUG
static void output_expolygons(const ExPolygons& expolys, const std::string &filename);
static void output_structures(const std::vector<Structure> &structures);
#endif // SLA_AUTOSUPPORTS_DEBUG
std::function<void(void)> m_throw_on_cancel;
const Eigen::MatrixXd& m_V;
const Eigen::MatrixXi& m_F;
const sla::EigenMesh3D& m_emesh;
};

363
src/libslic3r/SLA/SLABasePool.cpp
View file

@ -4,78 +4,177 @@
#include "boost/log/trivial.hpp"
#include "SLABoostAdapter.hpp"
#include "ClipperUtils.hpp"
#include "Tesselate.hpp"
// For debugging:
//#include <fstream>
//#include <libnest2d/tools/benchmark.h>
//#include "SVG.hpp"
//#include "benchmark.h"
namespace Slic3r { namespace sla {
/// Convert the triangulation output to an intermediate mesh.
Contour3D convert(const Polygons& triangles, coord_t z, bool dir) {
Pointf3s points;
points.reserve(3*triangles.size());
Indices indices;
indices.reserve(points.size());
for(auto& tr : triangles) {
auto c = coord_t(points.size()), b = c++, a = c++;
if(dir) indices.emplace_back(a, b, c);
else indices.emplace_back(c, b, a);
for(auto& p : tr.points) {
points.emplace_back(unscale(x(p), y(p), z));
}
}
return {points, indices};
}
Contour3D walls(const ExPolygon& floor_plate, const ExPolygon& ceiling,
double floor_z_mm, double ceiling_z_mm,
ThrowOnCancel thr)
/// This function will return a triangulation of a sheet connecting an upper
/// and a lower plate given as input polygons. It will not triangulate the
/// plates themselves only the sheet. The caller has to specify the lower and
/// upper z levels in world coordinates as well as the offset difference
/// between the sheets. If the lower_z_mm is higher than upper_z_mm or the
/// offset difference is negative, the resulting triangle orientation will be
/// reversed.
///
/// IMPORTANT: This is not a universal triangulation algorithm. It assumes
/// that the lower and upper polygons are offsetted versions of the same
/// original polygon. In general, it assumes that one of the polygons is
/// completely inside the other. The offset difference is the reference
/// distance from the inner polygon's perimeter to the outer polygon's
/// perimeter. The real distance will be variable as the clipper offset has
/// different strategies (rounding, etc...). This algorithm should have
/// O(2n + 3m) complexity where n is the number of upper vertices and m is the
/// number of lower vertices.
Contour3D walls(const Polygon& lower, const Polygon& upper,
double lower_z_mm, double upper_z_mm,
double offset_difference_mm, ThrowOnCancel thr)
{
using std::transform; using std::back_inserter;
ExPolygon poly;
poly.contour.points = floor_plate.contour.points;
poly.holes.emplace_back(ceiling.contour);
auto& h = poly.holes.front();
std::reverse(h.points.begin(), h.points.end());
Polygons tri = triangulate(poly);
Contour3D ret;
ret.points.reserve(tri.size() * 3);
double fz = floor_z_mm;
double cz = ceiling_z_mm;
auto& rp = ret.points;
auto& rpi = ret.indices;
ret.indices.reserve(tri.size() * 3);
if(upper.points.size() < 3 || lower.size() < 3) return ret;
coord_t idx = 0;
// The concept of the algorithm is relatively simple. It will try to find
// the closest vertices from the upper and the lower polygon and use those
// as starting points. Then it will create the triangles sequentially using
// an edge from the upper polygon and a vertex from the lower or vice versa,
// depending on the resulting triangle's quality.
// The quality is measured by a scalar value. So far it looks like it is
// enough to derive it from the slope of the triangle's two edges connecting
// the upper and the lower part. A reference slope is calculated from the
// height and the offset difference.
auto hlines = h.lines();
auto is_upper = [&hlines](const Point& p) {
return std::any_of(hlines.begin(), hlines.end(),
[&p](const Line& l) {
return l.distance_to(p) < mm(1e-6);
});
// Offset in the index array for the ceiling
const auto offs = upper.points.size();
// Shorthand for the vertex arrays
auto& upoints = upper.points, &lpoints = lower.points;
auto& rpts = ret.points; auto& rfaces = ret.indices;
// If the Z levels are flipped, or the offset difference is negative, we
// will interpret that as the triangles normals should be inverted.
bool inverted = upper_z_mm < lower_z_mm || offset_difference_mm < 0;
// Copy the points into the mesh, convert them from 2D to 3D
rpts.reserve(upoints.size() + lpoints.size());
rfaces.reserve(2*upoints.size() + 2*lpoints.size());
const double sf = SCALING_FACTOR;
for(auto& p : upoints) rpts.emplace_back(p.x()*sf, p.y()*sf, upper_z_mm);
for(auto& p : lpoints) rpts.emplace_back(p.x()*sf, p.y()*sf, lower_z_mm);
// Create pointing indices into vertex arrays. u-upper, l-lower
size_t uidx = 0, lidx = offs, unextidx = 1, lnextidx = offs + 1;
// Simple squared distance calculation.
auto distfn = [](const Vec3d& p1, const Vec3d& p2) {
auto p = p1 - p2; return p.transpose() * p;
};
std::for_each(tri.begin(), tri.end(),
[&rp, &rpi, thr, &idx, is_upper, fz, cz](const Polygon& pp)
{
thr(); // may throw if cancellation was requested
// We need to find the closest point on lower polygon to the first point on
// the upper polygon. These will be our starting points.
double distmin = std::numeric_limits<double>::max();
for(size_t l = lidx; l < rpts.size(); ++l) {
thr();
double d = distfn(rpts[l], rpts[uidx]);
if(d < distmin) { lidx = l; distmin = d; }
}
for(auto& p : pp.points)
if(is_upper(p))
rp.emplace_back(unscale(x(p), y(p), mm(cz)));
else rp.emplace_back(unscale(x(p), y(p), mm(fz)));
// Set up lnextidx to be ahead of lidx in cyclic mode
lnextidx = lidx + 1;
if(lnextidx == rpts.size()) lnextidx = offs;
coord_t a = idx++, b = idx++, c = idx++;
if(fz > cz) rpi.emplace_back(c, b, a);
else rpi.emplace_back(a, b, c);
});
// This will be the flip switch to toggle between upper and lower triangle
// creation mode
enum class Proceed {
UPPER, // A segment from the upper polygon and one vertex from the lower
LOWER // A segment from the lower polygon and one vertex from the upper
} proceed = Proceed::UPPER;
// Flags to help evaluating loop termination.
bool ustarted = false, lstarted = false;
// The variables for the fitness values, one for the actual and one for the
// previous.
double current_fit = 0, prev_fit = 0;
// Every triangle of the wall has two edges connecting the upper plate with
// the lower plate. From the length of these two edges and the zdiff we
// can calculate the momentary squared offset distance at a particular
// position on the wall. The average of the differences from the reference
// (squared) offset distance will give us the driving fitness value.
const double offsdiff2 = std::pow(offset_difference_mm, 2);
const double zdiff2 = std::pow(upper_z_mm - lower_z_mm, 2);
// Mark the current vertex iterator positions. If the iterators return to
// the same position, the loop can be terminated.
size_t uendidx = uidx, lendidx = lidx;
do { thr(); // check throw if canceled
prev_fit = current_fit;
switch(proceed) { // proceed depending on the current state
case Proceed::UPPER:
if(!ustarted || uidx != uendidx) { // there are vertices remaining
// Get the 3D vertices in order
const Vec3d& p_up1 = rpts[size_t(uidx)];
const Vec3d& p_low = rpts[size_t(lidx)];
const Vec3d& p_up2 = rpts[size_t(unextidx)];
// Calculate fitness: the average of the two connecting edges
double a = offsdiff2 - (distfn(p_up1, p_low) - zdiff2);
double b = offsdiff2 - (distfn(p_up2, p_low) - zdiff2);
current_fit = (std::abs(a) + std::abs(b)) / 2;
if(current_fit > prev_fit) { // fit is worse than previously
proceed = Proceed::LOWER;
} else { // good to go, create the triangle
inverted? rfaces.emplace_back(unextidx, lidx, uidx) :
rfaces.emplace_back(uidx, lidx, unextidx) ;
// Increment the iterators, rotate if necessary
++uidx; ++unextidx;
if(unextidx == offs) unextidx = 0;
if(uidx == offs) uidx = 0;
ustarted = true; // mark the movement of the iterators
// so that the comparison to uendidx can be made correctly
}
} else proceed = Proceed::LOWER;
break;
case Proceed::LOWER:
// Mode with lower segment, upper vertex. Same structure:
if(!lstarted || lidx != lendidx) {
const Vec3d& p_low1 = rpts[size_t(lidx)];
const Vec3d& p_low2 = rpts[size_t(lnextidx)];
const Vec3d& p_up = rpts[size_t(uidx)];
double a = offsdiff2 - (distfn(p_up, p_low1) - zdiff2);
double b = offsdiff2 - (distfn(p_up, p_low2) - zdiff2);
current_fit = (std::abs(a) + std::abs(b)) / 2;
if(current_fit > prev_fit) {
proceed = Proceed::UPPER;
} else {
inverted? rfaces.emplace_back(uidx, lnextidx, lidx) :
rfaces.emplace_back(lidx, lnextidx, uidx);
++lidx; ++lnextidx;
if(lnextidx == rpts.size()) lnextidx = offs;
if(lidx == rpts.size()) lidx = offs;
lstarted = true;
}
} else proceed = Proceed::UPPER;
break;
} // end of switch
} while(!ustarted || !lstarted || uidx != uendidx || lidx != lendidx);
return ret;
}
@ -207,20 +306,31 @@ ExPolygons unify(const ExPolygons& shapes) {
/// Only a debug function to generate top and bottom plates from a 2D shape.
/// It is not used in the algorithm directly.
inline Contour3D roofs(const ExPolygon& poly, coord_t z_distance) {
Polygons triangles = triangulate(poly);
auto lower = convert(triangles, 0, false);
auto upper = convert(triangles, z_distance, true);
lower.merge(upper);
return lower;
auto lower = triangulate_expolygon_3d(poly);
auto upper = triangulate_expolygon_3d(poly, z_distance*SCALING_FACTOR, true);
Contour3D ret;
ret.merge(lower); ret.merge(upper);
return ret;
}
/// This method will create a rounded edge around a flat polygon in 3d space.
/// 'base_plate' parameter is the target plate.
/// 'radius' is the radius of the edges.
/// 'degrees' is tells how much of a circle should be created as the rounding.
/// It should be in degrees, not radians.
/// 'ceilheight_mm' is the Z coordinate of the flat polygon in 3D space.
/// 'dir' Is the direction of the round edges: inward or outward
/// 'thr' Throws if a cancel signal was received
/// 'last_offset' An auxiliary output variable to save the last offsetted
/// version of 'base_plate'
/// 'last_height' An auxiliary output to save the last z coordinate of the
/// offsetted base_plate. In other words, where the rounded edges end.
Contour3D round_edges(const ExPolygon& base_plate,
double radius_mm,
double degrees,
double ceilheight_mm,
bool dir,
ThrowOnCancel throw_on_cancel,
ThrowOnCancel thr,
ExPolygon& last_offset, double& last_height)
{
auto ob = base_plate;
@ -236,10 +346,10 @@ Contour3D round_edges(const ExPolygon& base_plate,
// we use sin for x distance because we interpret the angle starting from
// PI/2
int tos = degrees < 90?
int(radius_mm*std::cos(degrees * PI / 180 - PI/2) / stepx) : steps;
int(radius_mm*std::cos(degrees * PI / 180 - PI/2) / stepx) : steps;
for(int i = 1; i <= tos; ++i) {
throw_on_cancel();
thr();
ob = base_plate;
@ -252,7 +362,8 @@ Contour3D round_edges(const ExPolygon& base_plate,
wh = ceilheight_mm - radius_mm + stepy;
Contour3D pwalls;
pwalls = walls(ob, ob_prev, wh, wh_prev, throw_on_cancel);
double prev_x = xx - (i - 1) * stepx;
pwalls = walls(ob.contour, ob_prev.contour, wh, wh_prev, s*prev_x, thr);
curvedwalls.merge(pwalls);
ob_prev = ob;
@ -264,7 +375,7 @@ Contour3D round_edges(const ExPolygon& base_plate,
int tos = int(tox / stepx);
for(int i = 1; i <= tos; ++i) {
throw_on_cancel();
thr();
ob = base_plate;
double r2 = radius_mm * radius_mm;
@ -275,7 +386,9 @@ Contour3D round_edges(const ExPolygon& base_plate,
wh = ceilheight_mm - radius_mm - stepy;
Contour3D pwalls;
pwalls = walls(ob_prev, ob, wh_prev, wh, throw_on_cancel);
double prev_x = xx - radius_mm + (i - 1)*stepx;
pwalls =
walls(ob_prev.contour, ob.contour, wh_prev, wh, s*prev_x, thr);
curvedwalls.merge(pwalls);
ob_prev = ob;
@ -291,15 +404,17 @@ Contour3D round_edges(const ExPolygon& base_plate,
/// Generating the concave part of the 3D pool with the bottom plate and the
/// side walls.
Contour3D inner_bed(const ExPolygon& poly, double depth_mm,
double begin_h_mm = 0) {
Polygons triangles = triangulate(poly);
Contour3D inner_bed(const ExPolygon& poly,
double depth_mm,
double begin_h_mm = 0)
{
Contour3D bottom;
Pointf3s triangles = triangulate_expolygon_3d(poly, -depth_mm + begin_h_mm);
bottom.merge(triangles);
coord_t depth = mm(depth_mm);
coord_t begin_h = mm(begin_h_mm);
auto bottom = convert(triangles, -depth + begin_h, false);
auto lines = poly.lines();
// Generate outer walls
@ -469,6 +584,9 @@ void base_plate(const TriangleMesh &mesh, ExPolygons &output, float h,
void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
const PoolConfig& cfg)
{
// for debugging:
// Benchmark bench;
// bench.start();
double mergedist = 2*(1.8*cfg.min_wall_thickness_mm + 4*cfg.edge_radius_mm)+
cfg.max_merge_distance_mm;
@ -478,27 +596,28 @@ void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
// serve as the bottom plate of the pad. We will offset this concave hull
// and then offset back the result with clipper with rounding edges ON. This
// trick will create a nice rounded pad shape.
auto concavehs = concave_hull(ground_layer, mergedist, cfg.throw_on_cancel);
ExPolygons concavehs = concave_hull(ground_layer, mergedist, cfg.throw_on_cancel);
const double thickness = cfg.min_wall_thickness_mm;
const double wingheight = cfg.min_wall_height_mm;
const double fullheight = wingheight + thickness;
const double tilt = PI/4;
const double tilt = cfg.wall_tilt;
const double wingdist = wingheight / std::tan(tilt);
// scaled values
const coord_t s_thickness = mm(thickness);
const coord_t s_eradius = mm(cfg.edge_radius_mm);
const coord_t s_safety_dist = 2*s_eradius + coord_t(0.8*s_thickness);
// const coord_t wheight = mm(cfg.min_wall_height_mm);
coord_t s_wingdist = mm(wingdist);
const coord_t s_wingdist = mm(wingdist);
auto& thrcl = cfg.throw_on_cancel;
Contour3D pool;
for(ExPolygon& concaveh : concavehs) {
if(concaveh.contour.points.empty()) return;
// Get rif of any holes in the concave hull output.
// Get rid of any holes in the concave hull output.
concaveh.holes.clear();
// Here lies the trick that does the smooting only with clipper offset
@ -508,15 +627,22 @@ void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
auto outer_base = concaveh;
outer_base.holes.clear();
offset(outer_base, s_safety_dist + s_wingdist + s_thickness);
auto inner_base = outer_base;
offset(inner_base, -(s_thickness + s_wingdist));
ExPolygon bottom_poly = outer_base;
bottom_poly.holes.clear();
if(s_wingdist > 0) offset(bottom_poly, -s_wingdist);
// Punching a hole in the top plate for the cavity
ExPolygon top_poly;
ExPolygon middle_base;
ExPolygon inner_base;
top_poly.contour = outer_base.contour;
if(wingheight > 0) {
inner_base = outer_base;
offset(inner_base, -(s_thickness + s_wingdist + s_eradius));
middle_base = outer_base;
offset(middle_base, -s_thickness);
top_poly.holes.emplace_back(middle_base.contour);
@ -524,8 +650,6 @@ void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
std::reverse(tph.begin(), tph.end());
}
Contour3D pool;
ExPolygon ob = outer_base; double wh = 0;
// now we will calculate the angle or portion of the circle from
@ -557,60 +681,53 @@ void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
// Generate the smoothed edge geometry
auto walledges = round_edges(ob,
r,
phi,
0, // z position of the input plane
true,
thrcl,
ob, wh);
pool.merge(walledges);
pool.merge(round_edges(ob,
r,
phi,
0, // z position of the input plane
true,
thrcl,
ob, wh));
// Now that we have the rounded edge connencting the top plate with
// Now that we have the rounded edge connecting the top plate with
// the outer side walls, we can generate and merge the sidewall geometry
auto pwalls = walls(ob, inner_base, wh, -fullheight, thrcl);
pool.merge(pwalls);
pool.merge(walls(ob.contour, bottom_poly.contour, wh, -fullheight,
wingdist, thrcl));
if(wingheight > 0) {
// Generate the smoothed edge geometry
auto cavityedges = round_edges(middle_base,
r,
phi - 90, // from tangent lines
0,
false,
thrcl,
ob, wh);
pool.merge(cavityedges);
pool.merge(round_edges(middle_base,
r,
phi - 90, // from tangent lines
0, // z position of the input plane
false,
thrcl,
ob, wh));
// Next is the cavity walls connecting to the top plate's
// artificially created hole.
auto cavitywalls = walls(inner_base, ob, -wingheight, wh, thrcl);
pool.merge(cavitywalls);
pool.merge(walls(inner_base.contour, ob.contour, -wingheight,
wh, -wingdist, thrcl));
}
// Now we need to triangulate the top and bottom plates as well as the
// cavity bottom plate which is the same as the bottom plate but it is
// eleveted by the thickness.
Polygons top_triangles, bottom_triangles;
// elevated by the thickness.
pool.merge(triangulate_expolygon_3d(top_poly));
pool.merge(triangulate_expolygon_3d(bottom_poly, -fullheight, true));
triangulate(top_poly, top_triangles);
triangulate(inner_base, bottom_triangles);
if(wingheight > 0)
pool.merge(triangulate_expolygon_3d(inner_base, -wingheight));
auto top_plate = convert(top_triangles, 0, false);
auto bottom_plate = convert(bottom_triangles, -mm(fullheight), true);
pool.merge(top_plate);
pool.merge(bottom_plate);
if(wingheight > 0) {
Polygons middle_triangles;
triangulate(inner_base, middle_triangles);
auto middle_plate = convert(middle_triangles, -mm(wingheight), false);
pool.merge(middle_plate);
}
out.merge(mesh(pool));
}
// For debugging:
// bench.stop();
// std::cout << "Pad creation time: " << bench.getElapsedSec() << std::endl;
// std::fstream fout("pad_debug.obj", std::fstream::out);
// if(fout.good()) pool.to_obj(fout);
out.merge(mesh(pool));
}
}

7
src/libslic3r/SLA/SLABasePool.hpp
View file

@ -3,6 +3,7 @@
#include <vector>
#include <functional>
#include <cmath>
namespace Slic3r {
@ -27,15 +28,17 @@ struct PoolConfig {
double min_wall_height_mm = 5;
double max_merge_distance_mm = 50;
double edge_radius_mm = 1;
double wall_tilt = std::atan(1.0); // Universal constant for Pi/4
ThrowOnCancel throw_on_cancel = [](){};
inline PoolConfig() {}
inline PoolConfig(double wt, double wh, double md, double er):
inline PoolConfig(double wt, double wh, double md, double er, double tilt):
min_wall_thickness_mm(wt),
min_wall_height_mm(wh),
max_merge_distance_mm(md),
edge_radius_mm(er) {}
edge_radius_mm(er),
wall_tilt(tilt) {}
};
/// Calculate the pool for the mesh for SLA printing

21
src/libslic3r/SLA/SLABoilerPlate.hpp
View file

@ -36,14 +36,6 @@ inline coord_t x(const Vec3crd& p) { return p(0); }
inline coord_t y(const Vec3crd& p) { return p(1); }
inline coord_t z(const Vec3crd& p) { return p(2); }
inline void triangulate(const ExPolygon& expoly, Polygons& triangles) {
expoly.triangulate_p2t(&triangles);
}
inline Polygons triangulate(const ExPolygon& expoly) {
Polygons tri; triangulate(expoly, tri); return tri;
}
using Indices = std::vector<Vec3crd>;
/// Intermediate struct for a 3D mesh
@ -63,6 +55,15 @@ struct Contour3D {
}
}
void merge(const Pointf3s& triangles) {
const size_t offs = points.size();
points.insert(points.end(), triangles.begin(), triangles.end());
indices.reserve(indices.size() + points.size() / 3);
for(int i = (int)offs; i < (int)points.size(); i += 3)
indices.emplace_back(i, i + 1, i + 2);
}
// Write the index triangle structure to OBJ file for debugging purposes.
void to_obj(std::ostream& stream) {
for(auto& p : points) {
@ -75,13 +76,9 @@ struct Contour3D {
}
};
//using PointSet = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::DontAlign>; //Eigen::MatrixXd;
using ClusterEl = std::vector<unsigned>;
using ClusteredPoints = std::vector<ClusterEl>;
/// Convert the triangulation output to an intermediate mesh.
Contour3D convert(const Polygons& triangles, coord_t z, bool dir);
/// Mesh from an existing contour.
inline TriangleMesh mesh(const Contour3D& ctour) {
return {ctour.points, ctour.indices};

137
src/libslic3r/SLA/SLACommon.hpp Normal file
View file

@ -0,0 +1,137 @@
#ifndef SLACOMMON_HPP
#define SLACOMMON_HPP
#include <Eigen/Geometry>
// #define SLIC3R_SLA_NEEDS_WINDTREE
namespace Slic3r {
// Typedefs from Point.hpp
typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> Vec3f;
typedef Eigen::Matrix<double, 3, 1, Eigen::DontAlign> Vec3d;
class TriangleMesh;
namespace sla {
struct SupportPoint {
Vec3f pos;
float head_front_radius;
bool is_new_island;
SupportPoint() :
pos(Vec3f::Zero()), head_front_radius(0.f), is_new_island(false) {}
SupportPoint(float pos_x, float pos_y, float pos_z, float head_radius, bool new_island) :
pos(pos_x, pos_y, pos_z), head_front_radius(head_radius), is_new_island(new_island) {}
SupportPoint(Vec3f position, float head_radius, bool new_island) :
pos(position), head_front_radius(head_radius), is_new_island(new_island) {}
SupportPoint(Eigen::Matrix<float, 5, 1, Eigen::DontAlign> data) :
pos(data(0), data(1), data(2)), head_front_radius(data(3)), is_new_island(data(4) != 0.f) {}
bool operator==(const SupportPoint& sp) const { return (pos==sp.pos) && head_front_radius==sp.head_front_radius && is_new_island==sp.is_new_island; }
bool operator!=(const SupportPoint& sp) const { return !(sp == (*this)); }
};
/// An index-triangle structure for libIGL functions. Also serves as an
/// alternative (raw) input format for the SLASupportTree
/*struct EigenMesh3D {
Eigen::MatrixXd V;
Eigen::MatrixXi F;
double ground_level = 0;
};*/
/// An index-triangle structure for libIGL functions. Also serves as an
/// alternative (raw) input format for the SLASupportTree
class EigenMesh3D {
class AABBImpl;
Eigen::MatrixXd m_V;
Eigen::MatrixXi m_F;
double m_ground_level = 0;
std::unique_ptr<AABBImpl> m_aabb;
public:
EigenMesh3D(const TriangleMesh&);
EigenMesh3D(const EigenMesh3D& other);
EigenMesh3D& operator=(const EigenMesh3D&);
~EigenMesh3D();
inline double ground_level() const { return m_ground_level; }
inline const Eigen::MatrixXd& V() const { return m_V; }
inline const Eigen::MatrixXi& F() const { return m_F; }
// Result of a raycast
class hit_result {
double m_t = std::numeric_limits<double>::infinity();
int m_face_id = -1;
const EigenMesh3D& m_mesh;
Vec3d m_dir;
inline hit_result(const EigenMesh3D& em): m_mesh(em) {}
friend class EigenMesh3D;
public:
inline double distance() const { return m_t; }
inline const Vec3d& direction() const { return m_dir; }
inline int face() const { return m_face_id; }
inline Vec3d normal() const {
if(m_face_id < 0) return {};
auto trindex = m_mesh.m_F.row(m_face_id);
const Vec3d& p1 = m_mesh.V().row(trindex(0));
const Vec3d& p2 = m_mesh.V().row(trindex(1));
const Vec3d& p3 = m_mesh.V().row(trindex(2));
Eigen::Vector3d U = p2 - p1;
Eigen::Vector3d V = p3 - p1;
return U.cross(V).normalized();
}
inline bool is_inside() {
return m_face_id >= 0 && normal().dot(m_dir) > 0;
}
};
// Casting a ray on the mesh, returns the distance where the hit occures.
hit_result query_ray_hit(const Vec3d &s, const Vec3d &dir) const;
class si_result {
double m_value;
int m_fidx;
Vec3d m_p;
si_result(double val, int i, const Vec3d& c):
m_value(val), m_fidx(i), m_p(c) {}
friend class EigenMesh3D;
public:
si_result() = delete;
double value() const { return m_value; }
operator double() const { return m_value; }
const Vec3d& point_on_mesh() const { return m_p; }
int F_idx() const { return m_fidx; }
};
#ifdef SLIC3R_SLA_NEEDS_WINDTREE
// The signed distance from a point to the mesh. Outputs the distance,
// the index of the triangle and the closest point in mesh coordinate space.
si_result signed_distance(const Vec3d& p) const;
bool inside(const Vec3d& p) const;
#endif /* SLIC3R_SLA_NEEDS_WINDTREE */
};
} // namespace sla
} // namespace Slic3r
#endif // SLASUPPORTTREE_HPP

11
src/libslic3r/SLA/SLASupportTree.cpp
View file

@ -551,10 +551,16 @@ enum { // For indexing Eigen vectors as v(X), v(Y), v(Z) instead of numbers
X, Y, Z
};
PointSet to_point_set(const std::vector<Vec3d> &v)
PointSet to_point_set(const std::vector<SupportPoint> &v)
{
PointSet ret(v.size(), 3);
{ long i = 0; for(const Vec3d& p : v) ret.row(i++) = p; }
long i = 0;
for(const SupportPoint& support_point : v) {
ret.row(i)(0) = support_point.pos(0);
ret.row(i)(1) = support_point.pos(1);
ret.row(i)(2) = support_point.pos(2);
++i;
}
return ret;
}
@ -679,6 +685,7 @@ double pinhead_mesh_intersect(const Vec3d& s,
return *mit;
}
// Checking bridge (pillar and stick as well) intersection with the model. If
// the function is used for headless sticks, the ins_check parameter have to be
// true as the beginning of the stick might be inside the model geometry.

85
src/libslic3r/SLA/SLASupportTree.hpp
View file

@ -7,6 +7,9 @@
#include <memory>
#include <Eigen/Geometry>
#include "SLACommon.hpp"
namespace Slic3r {
// Needed types from Point.hpp
@ -105,86 +108,6 @@ struct Controller {
std::function<void(void)> cancelfn = [](){};
};
/// An index-triangle structure for libIGL functions. Also serves as an
/// alternative (raw) input format for the SLASupportTree
class EigenMesh3D {
class AABBImpl;
Eigen::MatrixXd m_V;
Eigen::MatrixXi m_F;
double m_ground_level = 0;
std::unique_ptr<AABBImpl> m_aabb;
public:
EigenMesh3D(const TriangleMesh&);
EigenMesh3D(const EigenMesh3D& other);
EigenMesh3D& operator=(const EigenMesh3D&);
~EigenMesh3D();
inline double ground_level() const { return m_ground_level; }
inline const Eigen::MatrixXd& V() const { return m_V; }
inline const Eigen::MatrixXi& F() const { return m_F; }
// Result of a raycast
class hit_result {
double m_t = std::numeric_limits<double>::infinity();
int m_face_id = -1;
const EigenMesh3D& m_mesh;
Vec3d m_dir;
inline hit_result(const EigenMesh3D& em): m_mesh(em) {}
friend class EigenMesh3D;
public:
inline double distance() const { return m_t; }
inline int face() const { return m_face_id; }
inline Vec3d normal() const {
if(m_face_id < 0) return {};
auto trindex = m_mesh.m_F.row(m_face_id);
const Vec3d& p1 = m_mesh.V().row(trindex(0));
const Vec3d& p2 = m_mesh.V().row(trindex(1));
const Vec3d& p3 = m_mesh.V().row(trindex(2));
Eigen::Vector3d U = p2 - p1;
Eigen::Vector3d V = p3 - p1;
return U.cross(V).normalized();
}
inline bool is_inside() {
return m_face_id >= 0 && normal().dot(m_dir) > 0;
}
};
// Casting a ray on the mesh, returns the distance where the hit occures.
hit_result query_ray_hit(const Vec3d &s, const Vec3d &dir) const;
class si_result {
double m_value;
int m_fidx;
Vec3d m_p;
si_result(double val, int i, const Vec3d& c):
m_value(val), m_fidx(i), m_p(c) {}
friend class EigenMesh3D;
public:
si_result() = delete;
double value() const { return m_value; }
operator double() const { return m_value; }
const Vec3d& point_on_mesh() const { return m_p; }
int F_idx() const { return m_fidx; }
};
// The signed distance from a point to the mesh. Outputs the distance,
// the index of the triangle and the closest point in mesh coordinate space.
si_result signed_distance(const Vec3d& p) const;
bool inside(const Vec3d& p) const;
};
using PointSet = Eigen::MatrixXd;
//EigenMesh3D to_eigenmesh(const TriangleMesh& m);
@ -193,7 +116,7 @@ using PointSet = Eigen::MatrixXd;
//EigenMesh3D to_eigenmesh(const ModelObject& model);
// Simple conversion of 'vector of points' to an Eigen matrix
PointSet to_point_set(const std::vector<Vec3d>&);
PointSet to_point_set(const std::vector<sla::SupportPoint>&);
/* ************************************************************************** */

14
src/libslic3r/SLA/SLASupportTreeIGL.cpp
View file

@ -95,7 +95,9 @@ size_t SpatIndex::size() const
class EigenMesh3D::AABBImpl: public igl::AABB<Eigen::MatrixXd, 3> {
public:
#ifdef SLIC3R_SLA_NEEDS_WINDTREE
igl::WindingNumberAABB<Vec3d, Eigen::MatrixXd, Eigen::MatrixXi> windtree;
#endif /* SLIC3R_SLA_NEEDS_WINDTREE */
};
EigenMesh3D::EigenMesh3D(const TriangleMesh& tmesh): m_aabb(new AABBImpl()) {
@ -136,7 +138,9 @@ EigenMesh3D::EigenMesh3D(const TriangleMesh& tmesh): m_aabb(new AABBImpl()) {
// Build the AABB accelaration tree
m_aabb->init(m_V, m_F);
#ifdef SLIC3R_SLA_NEEDS_WINDTREE
m_aabb->windtree.set_mesh(m_V, m_F);
#endif /* SLIC3R_SLA_NEEDS_WINDTREE */
}
EigenMesh3D::~EigenMesh3D() {}
@ -168,6 +172,7 @@ EigenMesh3D::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
return ret;
}
#ifdef SLIC3R_SLA_NEEDS_WINDTREE
EigenMesh3D::si_result EigenMesh3D::signed_distance(const Vec3d &p) const {
double sign = 0; double sqdst = 0; int i = 0; Vec3d c;
igl::signed_distance_winding_number(*m_aabb, m_V, m_F, m_aabb->windtree,
@ -179,6 +184,7 @@ EigenMesh3D::si_result EigenMesh3D::signed_distance(const Vec3d &p) const {
bool EigenMesh3D::inside(const Vec3d &p) const {
return m_aabb->windtree.inside(p);
}
#endif /* SLIC3R_SLA_NEEDS_WINDTREE */
/* ****************************************************************************
* Misc functions
@ -199,9 +205,11 @@ template<class Vec> double distance(const Vec& pp1, const Vec& pp2) {
return std::sqrt(p.transpose() * p);
}
PointSet normals(const PointSet& points, const EigenMesh3D& mesh,
PointSet normals(const PointSet& points,
const EigenMesh3D& mesh,
double eps,
std::function<void()> throw_on_cancel) {
std::function<void()> throw_on_cancel)
{
if(points.rows() == 0 || mesh.V().rows() == 0 || mesh.F().rows() == 0)
return {};
@ -222,7 +230,7 @@ PointSet normals(const PointSet& points, const EigenMesh3D& mesh,
const Vec3d& p3 = mesh.V().row(trindex(2));
// We should check if the point lies on an edge of the hosting triangle.
// If it does than all the other triangles using the same two points
// If it does then all the other triangles using the same two points
// have to be searched and the final normal should be some kind of
// aggregation of the participating triangle normals. We should also
// consider the cases where the support point lies right on a vertex