diff --git a/src/libslic3r/TreeModelVolumes.cpp b/src/libslic3r/TreeModelVolumes.cpp index 3c7cd8af3..63b0ab7a3 100644 --- a/src/libslic3r/TreeModelVolumes.cpp +++ b/src/libslic3r/TreeModelVolumes.cpp @@ -355,6 +355,14 @@ const Polygons& TreeModelVolumes::getCollision(const coord_t orig_radius, LayerI return getCollision(orig_radius, layer_idx, min_xy_dist); } +// Get a collision area at a given layer for a radius that is a lower or equial to the key radius. +// It is expected that the collision area is precalculated for a given layer at least for the radius zero. +// Used for pushing tree supports away from object during the final Organic optimization step. +std::optional>> TreeModelVolumes::get_collision_lower_bound_area(LayerIndex layer_id, coord_t max_radius) const +{ + return m_collision_cache.get_lower_bound_area({ max_radius, layer_id }); +} + // Private. Only called internally by calculateAvoidance() and calculateAvoidanceToModel(), radius is already snapped to grid. const Polygons& TreeModelVolumes::getCollisionHolefree(coord_t radius, LayerIndex layer_idx) const { diff --git a/src/libslic3r/TreeModelVolumes.hpp b/src/libslic3r/TreeModelVolumes.hpp index 407025e1f..0c3c4cc8a 100644 --- a/src/libslic3r/TreeModelVolumes.hpp +++ b/src/libslic3r/TreeModelVolumes.hpp @@ -237,6 +237,11 @@ public: */ const Polygons& getCollision(const coord_t radius, LayerIndex layer_idx, bool min_xy_dist) const; + // Get a collision area at a given layer for a radius that is a lower or equial to the key radius. + // It is expected that the collision area is precalculated for a given layer at least for the radius zero. + // Used for pushing tree supports away from object during the final Organic optimization step. + std::optional>> get_collision_lower_bound_area(LayerIndex layer_id, coord_t max_radius) const; + /*! * \brief Provides the areas that have to be avoided by the tree's branches * in order to reach the build plate. diff --git a/src/libslic3r/TreeSupport.cpp b/src/libslic3r/TreeSupport.cpp index 7f5f497e5..9b373312f 100644 --- a/src/libslic3r/TreeSupport.cpp +++ b/src/libslic3r/TreeSupport.cpp @@ -8,6 +8,7 @@ #include "TreeSupport.hpp" #include "AABBTreeIndirect.hpp" +#include "AABBTreeLines.hpp" #include "BuildVolume.hpp" #include "ClipperUtils.hpp" #include "EdgeGrid.hpp" @@ -801,6 +802,14 @@ static double layer_z(const SlicingParameters &slicing_params, const size_t laye { return slicing_params.object_print_z_min + slicing_params.first_object_layer_height + layer_idx * slicing_params.layer_height; } +static LayerIndex layer_idx_ceil(const SlicingParameters &slicing_params, const double z) +{ + return LayerIndex(ceil((z - slicing_params.object_print_z_min - slicing_params.first_object_layer_height) / slicing_params.layer_height)); +} +static LayerIndex layer_idx_floor(const SlicingParameters &slicing_params, const double z) +{ + return LayerIndex(floor((z - slicing_params.object_print_z_min - slicing_params.first_object_layer_height) / slicing_params.layer_height)); +} static inline SupportGeneratorLayer& layer_initialize( SupportGeneratorLayer &layer_new, @@ -2533,6 +2542,7 @@ static void create_nodes_from_area( double radius_increase = config.getRadius(elem.state) - config.getRadius(parent.state); assert(radius_increase >= 0); double shift = (elem.state.result_on_layer - parent.state.result_on_layer).cast().norm(); + //FIXME this assert fails a lot. Is it correct? assert(shift < radius_increase + 2. * config.maximum_move_distance_slow); } } @@ -2557,6 +2567,7 @@ static void create_nodes_from_area( double radius_increase = config.getRadius(elem.state) - config.getRadius(parent.state); assert(radius_increase >= 0); double shift = (elem.state.result_on_layer - parent.state.result_on_layer).cast().norm(); + //FIXME this assert fails a lot. Is it correct? assert(shift < radius_increase + 2. * config.maximum_move_distance_slow); } } @@ -3378,6 +3389,306 @@ static void extrude_branch( } #endif +// #define TREE_SUPPORT_ORGANIC_NUDGE_NEW 1 + +#ifdef TREE_SUPPORT_ORGANIC_NUDGE_NEW +// New version using per layer AABB trees of lines for nudging spheres away from an object. +static void organic_smooth_branches_avoid_collisions( + const PrintObject &print_object, + const TreeModelVolumes &volumes, + const TreeSupportSettings &config, + std::vector &move_bounds, + const std::vector> &elements_with_link_down, + const std::vector &linear_data_layers) +{ + struct LayerCollisionCache { + coord_t min_element_radius{ std::numeric_limits::max() }; + bool min_element_radius_known() const { return this->min_element_radius != std::numeric_limits::max(); } + coord_t collision_radius{ 0 }; + std::vector lines; + AABBTreeIndirect::Tree<2, double> aabbtree_lines; + bool empty() const { return this->lines.empty(); } + }; + std::vector layer_collision_cache; + layer_collision_cache.reserve(1024); + const SlicingParameters &slicing_params = print_object.slicing_parameters(); + for (const std::pair& element : elements_with_link_down) { + LayerIndex layer_idx = element.first->state.layer_idx; + if (size_t num_layers = layer_idx + 1; num_layers > layer_collision_cache.size()) { + if (num_layers > layer_collision_cache.capacity()) + reserve_power_of_2(layer_collision_cache, num_layers); + layer_collision_cache.resize(num_layers, {}); + } + auto& l = layer_collision_cache[layer_idx]; + l.min_element_radius = std::min(l.min_element_radius, config.getRadius(element.first->state)); + } + for (LayerIndex layer_idx = 0; layer_idx < LayerIndex(layer_collision_cache.size()); ++layer_idx) + if (LayerCollisionCache& l = layer_collision_cache[layer_idx]; !l.min_element_radius_known()) + l.min_element_radius = 0; + else { + //FIXME + l.min_element_radius = 0; + std::optional>> res = volumes.get_collision_lower_bound_area(layer_idx, l.min_element_radius); + assert(res.has_value()); + l.collision_radius = res->first; + Lines alines = to_lines(res->second.get()); + l.lines.reserve(alines.size()); + for (const Line &line : alines) + l.lines.push_back({ unscaled(line.a), unscaled(line.b) }); + l.aabbtree_lines = AABBTreeLines::build_aabb_tree_over_indexed_lines(l.lines); + } + + struct CollisionSphere { + const SupportElement& element; + int element_below_id; + const bool locked; + float radius; + // Current position, when nudged away from the collision. + Vec3f position; + // Previous position, for Laplacian smoothing. + Vec3f prev_position; + // + Vec3f last_collision; + double last_collision_depth; + // Minimum Z for which the sphere collision will be evaluated. + // Limited by the minimum sloping angle and by the bottom of the tree. + float min_z{ -std::numeric_limits::max() }; + // Maximum Z for which the sphere collision will be evaluated. + // Limited by the minimum sloping angle and by the tip of the current branch. + float max_z{ std::numeric_limits::max() }; + uint32_t layer_begin; + uint32_t layer_end; + }; + + std::vector collision_spheres; + collision_spheres.reserve(elements_with_link_down.size()); + for (const std::pair &element_with_link : elements_with_link_down) { + const SupportElement &element = *element_with_link.first; + const int link_down = element_with_link.second; + collision_spheres.push_back({ + element, + link_down, + // locked + element.parents.empty() || (link_down == -1 && element.state.layer_idx > 0), + unscaled(config.getRadius(element.state)), + // 3D position + to_3d(unscaled(element.state.result_on_layer), float(layer_z(slicing_params, element.state.layer_idx))) + }); + // Update min_z coordinate to min_z of the tree below. + CollisionSphere &collision_sphere = collision_spheres.back(); + if (link_down != -1) { + const size_t offset_below = linear_data_layers[element.state.layer_idx - 1]; + collision_sphere.min_z = collision_spheres[offset_below + link_down].min_z; + } else + collision_sphere.min_z = collision_sphere.position.z(); + } + // Update max_z by propagating max_z from the tips of the branches. + for (int collision_sphere_id = int(collision_spheres.size()) - 1; collision_sphere_id >= 0; -- collision_sphere_id) { + CollisionSphere &collision_sphere = collision_spheres[collision_sphere_id]; + if (collision_sphere.element.parents.empty()) + // Tip + collision_sphere.max_z = collision_sphere.position.z(); + else { + // Below tip + const size_t offset_above = linear_data_layers[collision_sphere.element.state.layer_idx + 1]; + for (auto iparent : collision_sphere.element.parents) { + float parent_z = collision_spheres[offset_above + iparent].max_z; +// collision_sphere.max_z = collision_sphere.max_z == std::numeric_limits::max() ? parent_z : std::max(collision_sphere.max_z, parent_z); + collision_sphere.max_z = std::min(collision_sphere.max_z, parent_z); + } + } + } + // Update min_z / max_z to limit the search Z span of a given sphere for collision detection. + for (CollisionSphere &collision_sphere : collision_spheres) { + //FIXME limit the collision span by the tree slope. + collision_sphere.min_z = std::max(collision_sphere.min_z, collision_sphere.position.z() - collision_sphere.radius); + collision_sphere.max_z = std::min(collision_sphere.max_z, collision_sphere.position.z() + collision_sphere.radius); + collision_sphere.layer_begin = std::min(collision_sphere.element.state.layer_idx, layer_idx_ceil(slicing_params, collision_sphere.min_z)); + collision_sphere.layer_end = std::max(collision_sphere.element.state.layer_idx, layer_idx_floor(slicing_params, collision_sphere.max_z)) + 1; + } + + static constexpr const double collision_extra_gap = 0.1; + static constexpr const double max_nudge_collision_avoidance = 0.2; + static constexpr const double max_nudge_smoothing = 0.2; + static constexpr const size_t num_iter = 100; // 1000; + for (size_t iter = 0; iter < num_iter; ++ iter) { + // Back up prev position before Laplacian smoothing. + for (CollisionSphere &collision_sphere : collision_spheres) + collision_sphere.prev_position = collision_sphere.position; + std::atomic num_moved{ 0 }; + tbb::parallel_for(tbb::blocked_range(0, collision_spheres.size()), + [&collision_spheres, &layer_collision_cache, &slicing_params, &move_bounds, &linear_data_layers, &num_moved](const tbb::blocked_range range) { + for (size_t collision_sphere_id = range.begin(); collision_sphere_id < range.end(); ++ collision_sphere_id) + if (CollisionSphere &collision_sphere = collision_spheres[collision_sphere_id]; ! collision_sphere.locked) { + // Calculate collision of multiple 2D layers against a collision sphere. + collision_sphere.last_collision_depth = - std::numeric_limits::max(); + for (uint32_t layer_id = collision_sphere.layer_begin; layer_id != collision_sphere.layer_end; ++ layer_id) { + double dz = (layer_id - collision_sphere.element.state.layer_idx) * slicing_params.layer_height; + if (double r2 = sqr(collision_sphere.radius) - sqr(dz); r2 > 0) { + if (const LayerCollisionCache &layer_collision_cache_item = layer_collision_cache[layer_id]; ! layer_collision_cache.empty()) { + size_t hit_idx_out; + Vec2d hit_point_out; + double dist = sqrt(AABBTreeLines::squared_distance_to_indexed_lines( + layer_collision_cache_item.lines, layer_collision_cache_item.aabbtree_lines, Vec2d(to_2d(collision_sphere.position).cast()), + hit_idx_out, hit_point_out, r2)); + double collision_depth = sqrt(r2) - dist; + if (collision_depth > collision_sphere.last_collision_depth) { + collision_sphere.last_collision_depth = collision_depth; + collision_sphere.last_collision = to_3d(hit_point_out.cast(), float(layer_z(slicing_params, layer_id))); + } + } + } + } + if (collision_sphere.last_collision_depth > 0) { + // Collision detected to be removed. + // Nudge the circle center away from the collision. + if (collision_sphere.last_collision_depth > EPSILON) + // a little bit of hysteresis to detect end of + ++ num_moved; + // Shift by maximum 2mm. + double nudge_dist = std::min(std::max(0., collision_sphere.last_collision_depth + collision_extra_gap), max_nudge_collision_avoidance); + Vec2d nudge_vector = (to_2d(collision_sphere.position) - to_2d(collision_sphere.last_collision)).cast().normalized() * nudge_dist; + collision_sphere.position.head<2>() += (nudge_vector * nudge_dist).cast(); + } + // Laplacian smoothing + Vec2d avg{ 0, 0 }; + const SupportElements &above = move_bounds[collision_sphere.element.state.layer_idx + 1]; + const size_t offset_above = linear_data_layers[collision_sphere.element.state.layer_idx + 1]; + double weight = 0.; + for (auto iparent : collision_sphere.element.parents) { + double w = collision_sphere.radius; + avg += w * to_2d(collision_spheres[offset_above + iparent].prev_position.cast()); + weight += w; + } + if (collision_sphere.element_below_id != -1) { + const size_t offset_below = linear_data_layers[collision_sphere.element.state.layer_idx - 1]; + const double w = weight; // config.getRadius(move_bounds[element.state.layer_idx - 1][below].state); + avg += w * to_2d(collision_spheres[offset_below + collision_sphere.element_below_id].prev_position.cast()); + weight += w; + } + avg /= weight; + static constexpr const double smoothing_factor = 0.5; + Vec2d old_pos = to_2d(collision_sphere.position).cast(); + Vec2d new_pos = (1. - smoothing_factor) * old_pos + smoothing_factor * avg; + Vec2d shift = new_pos - old_pos; + double nudge_dist_max = shift.norm(); + // Shift by maximum 1mm, less than the collision avoidance factor. + double nudge_dist = std::min(std::max(0., nudge_dist_max), max_nudge_smoothing); + collision_sphere.position.head<2>() += (shift.normalized() * nudge_dist).cast(); + } + }); + // printf("iteration: %d, moved: %d\n", int(iter), int(num_moved)); + if (num_moved == 0) + break; + } + + for (size_t i = 0; i < collision_spheres.size(); ++ i) + elements_with_link_down[i].first->state.result_on_layer = scaled(to_2d(collision_spheres[i].position)); +} +#else // TREE_SUPPORT_ORGANIC_NUDGE_NEW +// Old version using OpenVDB, works but it is extremely slow for complex meshes. +static void organic_smooth_branches_avoid_collisions( + const PrintObject &print_object, + const TreeModelVolumes &volumes, + const TreeSupportSettings &config, + std::vector &move_bounds, + const std::vector> &elements_with_link_down, + const std::vector &linear_data_layers) +{ + TriangleMesh mesh = print_object.model_object()->raw_mesh(); + mesh.transform(print_object.trafo_centered()); + double scale = 10.; + openvdb::FloatGrid::Ptr grid = mesh_to_grid(mesh.its, openvdb::math::Transform{}, scale, 0., 0.); + std::unique_ptr> closest_surface_point = openvdb::tools::ClosestSurfacePoint::create(*grid); + std::vector pts, prev, projections; + std::vector distances; + for (const std::pair& element : elements_with_link_down) { + Vec3d pt = to_3d(unscaled(element.first->state.result_on_layer), layer_z(print_object.slicing_parameters(), element.first->state.layer_idx)) * scale; + pts.push_back({ pt.x(), pt.y(), pt.z() }); + } + + const double collision_extra_gap = 1. * scale; + const double max_nudge_collision_avoidance = 2. * scale; + const double max_nudge_smoothing = 1. * scale; + + static constexpr const size_t num_iter = 100; // 1000; + for (size_t iter = 0; iter < num_iter; ++ iter) { + prev = pts; + projections = pts; + distances.assign(pts.size(), std::numeric_limits::max()); + closest_surface_point->searchAndReplace(projections, distances); + size_t num_moved = 0; + for (size_t i = 0; i < projections.size(); ++ i) { + const SupportElement &element = *elements_with_link_down[i].first; + const int below = elements_with_link_down[i].second; + const bool locked = below == -1 && element.state.layer_idx > 0; + if (! locked && pts[i] != projections[i]) { + // Nudge the circle center away from the collision. + Vec3d v{ projections[i].x() - pts[i].x(), projections[i].y() - pts[i].y(), projections[i].z() - pts[i].z() }; + double depth = v.norm(); + assert(std::abs(distances[i] - depth) < EPSILON); + double radius = unscaled(config.getRadius(element.state)) * scale; + if (depth < radius) { + // Collision detected to be removed. + ++ num_moved; + double dxy = sqrt(sqr(radius) - sqr(v.z())); + double nudge_dist_max = dxy - std::hypot(v.x(), v.y()) + //FIXME 1mm gap + + collision_extra_gap; + // Shift by maximum 2mm. + double nudge_dist = std::min(std::max(0., nudge_dist_max), max_nudge_collision_avoidance); + Vec2d nudge_v = to_2d(v).normalized() * (- nudge_dist); + pts[i].x() += nudge_v.x(); + pts[i].y() += nudge_v.y(); + } + } + // Laplacian smoothing + if (! locked && ! element.parents.empty()) { + Vec2d avg{ 0, 0 }; + const SupportElements &above = move_bounds[element.state.layer_idx + 1]; + const size_t offset_above = linear_data_layers[element.state.layer_idx + 1]; + double weight = 0.; + for (auto iparent : element.parents) { + double w = config.getRadius(above[iparent].state); + avg.x() += w * prev[offset_above + iparent].x(); + avg.y() += w * prev[offset_above + iparent].y(); + weight += w; + } + size_t cnt = element.parents.size(); + if (below != -1) { + const size_t offset_below = linear_data_layers[element.state.layer_idx - 1]; + const double w = weight; // config.getRadius(move_bounds[element.state.layer_idx - 1][below].state); + avg.x() += w * prev[offset_below + below].x(); + avg.y() += w * prev[offset_below + below].y(); + ++ cnt; + weight += w; + } + //avg /= double(cnt); + avg /= weight; + static constexpr const double smoothing_factor = 0.5; + Vec2d old_pos{ pts[i].x(), pts[i].y() }; + Vec2d new_pos = (1. - smoothing_factor) * old_pos + smoothing_factor * avg; + Vec2d shift = new_pos - old_pos; + double nudge_dist_max = shift.norm(); + // Shift by maximum 1mm, less than the collision avoidance factor. + double nudge_dist = std::min(std::max(0., nudge_dist_max), max_nudge_smoothing); + Vec2d nudge_v = shift.normalized() * nudge_dist; + pts[i].x() += nudge_v.x(); + pts[i].y() += nudge_v.y(); + } + } +// printf("iteration: %d, moved: %d\n", int(iter), int(num_moved)); + if (num_moved == 0) + break; + } + + for (size_t i = 0; i < projections.size(); ++ i) { + elements_with_link_down[i].first->state.result_on_layer.x() = scaled(pts[i].x()) / scale; + elements_with_link_down[i].first->state.result_on_layer.y() = scaled(pts[i].y()) / scale; + } +} +#endif // TREE_SUPPORT_ORGANIC_NUDGE_NEW + static void draw_branches( PrintObject &print_object, const TreeModelVolumes &volumes, @@ -3392,8 +3703,6 @@ static void draw_branches( { static int irun = 0; - const SlicingParameters& slicing_params = print_object.slicing_parameters(); - // All SupportElements are put into a layer independent storage to improve parallelization. std::vector> elements_with_link_down; std::vector linear_data_layers; @@ -3432,102 +3741,7 @@ static void draw_branches( } } - std::unique_ptr> closest_surface_point; - { - TriangleMesh mesh = print_object.model_object()->raw_mesh(); - mesh.transform(print_object.trafo_centered()); - double scale = 10.; - openvdb::FloatGrid::Ptr grid = mesh_to_grid(mesh.its, openvdb::math::Transform{}, scale, 0., 0.); - closest_surface_point = openvdb::tools::ClosestSurfacePoint::create(*grid); - std::vector pts, prev, projections; - std::vector distances; - for (const std::pair &element : elements_with_link_down) { - Vec3d pt = to_3d(unscaled(element.first->state.result_on_layer), layer_z(slicing_params, element.first->state.layer_idx)) * scale; - pts.push_back({ pt.x(), pt.y(), pt.z() }); - } - - const double collision_extra_gap = 1. * scale; - const double max_nudge_collision_avoidance = 2. * scale; - const double max_nudge_smoothing = 1. * scale; - - static constexpr const size_t num_iter = 100; // 1000; - for (size_t iter = 0; iter < num_iter; ++ iter) { - prev = pts; - projections = pts; - distances.assign(pts.size(), std::numeric_limits::max()); - closest_surface_point->searchAndReplace(projections, distances); - size_t num_moved = 0; - for (size_t i = 0; i < projections.size(); ++ i) { - const SupportElement &element = *elements_with_link_down[i].first; - const int below = elements_with_link_down[i].second; - const bool locked = below == -1 && element.state.layer_idx > 0; - if (! locked && pts[i] != projections[i]) { - // Nudge the circle center away from the collision. - Vec3d v{ projections[i].x() - pts[i].x(), projections[i].y() - pts[i].y(), projections[i].z() - pts[i].z() }; - double depth = v.norm(); - assert(std::abs(distances[i] - depth) < EPSILON); - double radius = unscaled(config.getRadius(element.state)) * scale; - if (depth < radius) { - // Collision detected to be removed. - ++ num_moved; - double dxy = sqrt(sqr(radius) - sqr(v.z())); - double nudge_dist_max = dxy - std::hypot(v.x(), v.y()) - //FIXME 1mm gap - + collision_extra_gap; - // Shift by maximum 2mm. - double nudge_dist = std::min(std::max(0., nudge_dist_max), max_nudge_collision_avoidance); - Vec2d nudge_v = to_2d(v).normalized() * (- nudge_dist); - pts[i].x() += nudge_v.x(); - pts[i].y() += nudge_v.y(); - } - } - // Laplacian smoothing - if (! locked && ! element.parents.empty()) { - Vec2d avg{ 0, 0 }; - const SupportElements &above = move_bounds[element.state.layer_idx + 1]; - const size_t offset_above = linear_data_layers[element.state.layer_idx + 1]; - double weight = 0.; - for (auto iparent : element.parents) { - double w = config.getRadius(above[iparent].state); - avg.x() += w * prev[offset_above + iparent].x(); - avg.y() += w * prev[offset_above + iparent].y(); - weight += w; - } - size_t cnt = element.parents.size(); - if (below != -1) { - const size_t offset_below = linear_data_layers[element.state.layer_idx - 1]; - const double w = weight; // config.getRadius(move_bounds[element.state.layer_idx - 1][below].state); - avg.x() += w * prev[offset_below + below].x(); - avg.y() += w * prev[offset_below + below].y(); - ++ cnt; - weight += w; - } - //avg /= double(cnt); - avg /= weight; - static constexpr const double smoothing_factor = 0.5; - Vec2d old_pos{ pts[i].x(), pts[i].y() }; - Vec2d new_pos = (1. - smoothing_factor) * old_pos + smoothing_factor * avg; - Vec2d shift = new_pos - old_pos; - double nudge_dist_max = shift.norm(); - // Shift by maximum 1mm, less than the collision avoidance factor. - double nudge_dist = std::min(std::max(0., nudge_dist_max), max_nudge_smoothing); - Vec2d nudge_v = shift.normalized() * nudge_dist; - pts[i].x() += nudge_v.x(); - pts[i].y() += nudge_v.y(); - } - } -// printf("iteration: %d, moved: %d\n", int(iter), int(num_moved)); - if (num_moved == 0) - break; - } - -#if 1 - for (size_t i = 0; i < projections.size(); ++ i) { - elements_with_link_down[i].first->state.result_on_layer.x() = scaled(pts[i].x()) / scale; - elements_with_link_down[i].first->state.result_on_layer.y() = scaled(pts[i].y()) / scale; - } -#endif - } + organic_smooth_branches_avoid_collisions(print_object, volumes, config, move_bounds, elements_with_link_down, linear_data_layers); std::vector support_layer_storage(move_bounds.size()); std::vector support_roof_storage(move_bounds.size()); @@ -3539,6 +3753,7 @@ static void draw_branches( // Traverse all nodes, generate tubes. // Traversal stack with nodes and thier current parent + const SlicingParameters &slicing_params = print_object.slicing_parameters(); std::vector path; indexed_triangle_set cummulative_mesh; indexed_triangle_set partial_mesh;