From 4cec347511938f58632cf04a45bded2ecced182e Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Thu, 25 May 2023 10:02:26 +0200 Subject: [PATCH] WIP Organic Supports: New smoothing of organic trees. --- src/libslic3r/Support/OrganicSupport.cpp | 428 +++++++++++++++++++++++ 1 file changed, 428 insertions(+) diff --git a/src/libslic3r/Support/OrganicSupport.cpp b/src/libslic3r/Support/OrganicSupport.cpp index b26a3c730..ca92629ff 100644 --- a/src/libslic3r/Support/OrganicSupport.cpp +++ b/src/libslic3r/Support/OrganicSupport.cpp @@ -26,6 +26,434 @@ namespace Slic3r namespace FFFTreeSupport { +// Single slice through a single branch or trough a number of branches. +struct Slice +{ + // All polygons collected for this slice. + Polygons polygons; + // All bottom contacts collected for this slice. + Polygons bottom_contacts; + // How many branches were merged in this slice? Used to decide whether ClipperLib union is needed. + size_t num_branches{ 0 }; +}; + +struct Element +{ + // Current position of the centerline including the Z coordinate, unscaled. + Vec3f position; + float radius; + + // Index of this layer, including the raft layers. + LayerIndex layer_idx; + + // Limits where the centerline could be placed at the current layer Z. + Polygons influence_area; + + // Locked node should not be moved. Locked nodes are at the top of an object or at the tips of branches. + bool locked; + + // Previous position, for Laplacian smoothing, unscaled. + Vec3f prev_position; + + // For sphere tracing and other collision detection optimizations. + Vec3f last_collision; + double last_collision_depth; + + struct CollisionSphere { + // 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() }; + // Span of layers to test collision of this sphere against. + uint32_t layer_begin; + uint32_t layer_end; + }; + + CollisionSphere collision_sphere; +}; + +struct Branch; + +struct Bifurcation +{ + Branch *branch; + double area; +}; + +// Single branch of a tree. +struct Branch +{ + std::vector path; + + using Bifurcations = +#ifdef NDEBUG + // To reduce memory allocation in release mode. + boost::container::small_vector; +#else // NDEBUG + // To ease debugging. + std::vector; +#endif // NDEBUG + + Bifurcations up; + Bifurcation down; + + // How many of the thick up branches are considered continuation of the trunk? + // These will be smoothed out together. + size_t num_up_trunk; + + bool has_root() const { return this->down.branch == nullptr; } + bool has_tip() const { return this->up.empty(); } +}; + +struct Tree +{ + // Branches: Store of all branches. + // The first branch is the root of the tree. + Slic3r::deque branches; + + Branch& root() { return branches.front(); } + const Branch& root() const { return branches.front(); } + + // Result of slicing the branches. + std::vector slices; + // First layer index of the first slice in the vector above. + LayerIndex first_layer_id{ -1 }; +}; + +using Forest = std::vector; +using Trees = std::vector; + +Element to_tree_element(const TreeSupportSettings &config, const SlicingParameters &slicing_params, SupportElement &element, bool is_root) +{ + Element out; + out.position = to_3d(unscaled(element.state.result_on_layer), float(layer_z(slicing_params, config, element.state.layer_idx))); + out.radius = support_element_radius(config, element); + out.layer_idx = element.state.layer_idx; + out.influence_area = std::move(element.influence_area); + out.locked = (is_root && element.state.layer_idx > 0) || element.state.locked(); + return out; +} + +// Convert move bounds into a forest of trees, each tree made of a graph of branches and bifurcation points. +// Destroys move_bounds. +Forest make_forest(const TreeSupportSettings &config, const SlicingParameters &slicing_params, std::vector &&move_bounds) +{ + struct TreeVisitor { + void visit_recursive(std::vector &move_bounds, SupportElement &start_element, Branch *parent_branch, Tree &out) const { + assert(! start_element.state.marked && ! start_element.parents.empty()); + // Collect elements up to a bifurcation above. + start_element.state.marked = true; + // For each branch bifurcating from this point: + SupportElements &layer = move_bounds[start_element.state.layer_idx]; + SupportElements &layer_above = move_bounds[start_element.state.layer_idx + 1]; + for (size_t parent_idx = 0; parent_idx < start_element.parents.size(); ++ parent_idx) { + Branch branch; + if (parent_branch) + // Duplicate the last element of the trunk below. + // If this branch has a smaller diameter than the trunk below, its centerline will not be aligned with the centerline of the trunk. + branch.path.emplace_back(parent_branch->path.back()); + branch.path.emplace_back(to_tree_element(config, slicing_params, start_element, parent_branch == nullptr)); + // Traverse each branch until it branches again. + SupportElement &first_parent = layer_above[start_element.parents[parent_idx]]; + assert(! first_parent.state.marked); + assert(branch.path.back().layer_idx + 1 == first_parent.state.layer_idx); + branch.path.emplace_back(to_tree_element(config, slicing_params, first_parent, false)); + if (first_parent.parents.size() < 2) + first_parent.state.marked = true; + SupportElement *next_branch = nullptr; + if (first_parent.parents.size() == 1) { + for (SupportElement *parent = &first_parent;;) { + assert(parent->state.marked); + SupportElement &next_parent = move_bounds[parent->state.layer_idx + 1][parent->parents.front()]; + assert(! next_parent.state.marked); + assert(branch.path.back().layer_idx + 1 == next_parent.state.layer_idx); + branch.path.emplace_back(to_tree_element(config, slicing_params, next_parent, false)); + if (next_parent.parents.size() > 1) { + // Branching point was reached. + next_branch = &next_parent; + break; + } + next_parent.state.marked = true; + if (next_parent.parents.size() == 0) + // Tip is reached. + break; + parent = &next_parent; + } + } else if (first_parent.parents.size() > 1) + // Branching point was reached. + next_branch = &first_parent; + assert(branch.path.size() >= 2); + assert(next_branch == nullptr || ! next_branch->state.marked); + out.branches.emplace_back(std::move(branch)); + Branch *pbranch = &out.branches.back(); + if (parent_branch) { + parent_branch->up.push_back({ pbranch }); + pbranch->down = { parent_branch }; + } + if (next_branch) + this->visit_recursive(move_bounds, *next_branch, pbranch, out); + } + + if (parent_branch) { + // Update initial radii of thin branches merging with a trunk. + auto it_up_max_r = std::max_element(parent_branch->up.begin(), parent_branch->up.end(), + [](const Bifurcation &l, const Bifurcation &r){ return l.branch->path[1].radius < r.branch->path[1].radius; }); + const float r1 = it_up_max_r->branch->path[1].radius; + const float radius_increment = unscaled(config.branch_radius_increase_per_layer); + for (auto it = parent_branch->up.begin(); it != parent_branch->up.end(); ++ it) + if (it != it_up_max_r) { + Element &el = it->branch->path.front(); + Element &el2 = it->branch->path[1]; + if (! is_approx(r1, el2.radius)) + el.radius = std::min(el.radius, el2.radius + radius_increment); + } + // Sort children of parent_branch by decreasing radius. + std::sort(parent_branch->up.begin(), parent_branch->up.end(), + [](const Bifurcation &l, const Bifurcation &r){ return l.branch->path.front().radius > r.branch->path.front().radius; }); + // Update number of branches to be considered a continuation of the trunk during smoothing. + { + const float r_trunk = 0.75 * it_up_max_r->branch->path.front().radius; + parent_branch->num_up_trunk = 0; + for (const Bifurcation& up : parent_branch->up) + if (up.branch->path.front().radius < r_trunk) + break; + else + ++ parent_branch->num_up_trunk; + } + } + } + + const TreeSupportSettings &config; + const SlicingParameters &slicing_params; + }; + + TreeVisitor visitor{ config, slicing_params }; + + for (SupportElements &elements : move_bounds) + for (SupportElement &el : elements) + el.state.marked = false; + + Trees trees; + for (LayerIndex layer_idx = 0; layer_idx + 1 < LayerIndex(move_bounds.size()); ++ layer_idx) { + for (SupportElement &start_element : move_bounds[layer_idx]) { + if (! start_element.state.marked && ! start_element.parents.empty()) { +#if 0 + { + // Verify that this node is a root, such that there is no element in the layer below + // that points to it. + int ielement = &start_element - move_bounds.data(); + int found = 0; + if (layer_idx > 0) { + for (auto &el : move_bounds[layer_idx - 1]) { + for (auto iparent : el.parents) + if (iparent == ielement) + ++ found; + } + if (found != 0) + printf("Found: %d\n", found); + } + } +#endif + trees.push_back({}); + visitor.visit_recursive(move_bounds, start_element, nullptr, trees.back()); + assert(! trees.back().branches.empty()); + assert(! trees.back().branches.front().path.empty()); +#if 0 + // Debugging: Only build trees with specific properties. + if (start_element.state.lost) { + } + else if (start_element.state.verylost) { + } + else + trees.pop_back(); +#endif + } + } + } + +#if 1 + move_bounds.clear(); +#else + for (SupportElements &elements : move_bounds) + for (SupportElement &el : elements) + el.state.marked = false; +#endif + + return trees; +} + +// Move bounds were propagated top to bottom. At each joint of branches the move bounds were reduced significantly. +// Now reflect the reduction of tree space by propagating the reduction of tree centerline space +// bottom-up starting with the bottom-most joint. +void trim_influence_areas_bottom_up(Forest &forest, const float dxy_dlayer) +{ + struct Trimmer { + static void trim_recursive(Branch &branch, const float delta_r, const float dxy_dlayer) { + assert(delta_r >= 0); + if (delta_r > 0) + branch.path.front().influence_area = offset(branch.path.front().influence_area, delta_r); + for (size_t i = 1; i < branch.path.size(); ++ i) + branch.path[i].influence_area = intersection(branch.path[i].influence_area, offset(branch.path[i - 1].influence_area, dxy_dlayer)); + const float r0 = branch.path.back().radius; + for (Bifurcation &up : branch.up) { + up.branch->path.front().influence_area = branch.path.back().influence_area; + trim_recursive(*up.branch, r0 - up.branch->path.front().radius, dxy_dlayer); + } + } + }; + + for (Tree &tree : forest) { + Branch &root = tree.root(); + const float r0 = root.path.back().radius; + for (Bifurcation &up : root.up) + Trimmer::trim_recursive(*up.branch, r0 - up.branch->path.front().radius, dxy_dlayer); + } +} + +// Straighten up and smooth centerlines inside their influence areas. +void smooth_trees_inside_influence_areas(Branch &root, bool is_root) +{ + // Smooth the subtree: + // + // Apply laplacian and bilaplacian smoothing inside a branch, + // apply laplacian smoothing only at a bifurcation point. + // + // Applying a bilaplacian smoothing inside a branch should ensure curvature of the brach to be lower + // than the radius at each particular point of the centerline, + // while omitting bilaplacian smoothing at bifurcation points will create sharp bifurcations. + // Sharp bifurcations have a smaller volume, but just a tiny bit larger surfaces than smooth bifurcations + // where each continuation of the trunk satifies the path radius > centerline element radius. + const size_t num_iterations = 100; + struct StackElement { + Branch &branch; + size_t idx_up; + }; + std::vector stack; + + auto adjust_position = [](Element &el, Vec2f new_pos) { + Point new_pos_scaled = scaled(new_pos); + if (! contains(el.influence_area, new_pos_scaled)) { + int64_t min_dist = std::numeric_limits::max(); + Point min_proj_scaled; + for (const Polygon& polygon : el.influence_area) { + Point proj_scaled = polygon.point_projection(new_pos_scaled); + if (int64_t dist = (proj_scaled - new_pos_scaled).cast().squaredNorm(); dist < min_dist) { + min_dist = dist; + min_proj_scaled = proj_scaled; + } + } + new_pos = unscaled(min_proj_scaled); + } + el.position.head<2>() = new_pos; + }; + + for (size_t iter = 0; iter < num_iterations; ++ iter) { + // 1) Back-up the current positions. + stack.push_back({ root, 0 }); + while (! stack.empty()) { + StackElement &state = stack.back(); + if (state.idx_up == state.branch.num_up_trunk) { + // Process this path. + for (auto &el : state.branch.path) + el.prev_position = el.position; + stack.pop_back(); + } else { + // Open another up node of this branch. + stack.push_back({ *state.branch.up[state.idx_up].branch, 0 }); + ++ state.idx_up; + } + } + // 2) Calculate new position. + stack.push_back({ root, 0 }); + while (! stack.empty()) { + StackElement &state = stack.back(); + if (state.idx_up == state.branch.num_up_trunk) { + // Process this path. + for (size_t i = 1; i + 1 < state.branch.path.size(); ++ i) + if (auto &el = state.branch.path[i]; ! el.locked) { + // Laplacian smoothing with 0.5 weight. + const Vec3f &p0 = state.branch.path[i - 1].prev_position; + const Vec3f &p1 = el.prev_position; + const Vec3f &p2 = state.branch.path[i + 1].prev_position; + adjust_position(el, 0.5 * p1.head<2>() + 0.25 * (p0.head<2>() + p2.head<2>())); +#if 0 + // Only apply bilaplacian smoothing if the current curvature is smaller than el.radius. + // Interpolate p0, p1, p2 with a circle. + // First project p0, p1, p2 into a common plane. + const Vec3f n = (p1 - p0).cross(p2 - p1); + const Vec3f y = Vec3f(n.y(), n.x(), 0).normalized(); + const Vec2f q0{ p0.z(), p0.dot(y) }; + const Vec2f q1{ p1.z(), p1.dot(y) }; + const Vec2f q2{ p2.z(), p2.dot(y) }; + // Interpolate q0, q1, q2 with a circle, calculate its radius. + Vec2f b = q1 - q0; + Vec2f c = q2 - q0; + float lb = b.squaredNorm(); + float lc = c.squaredNorm(); + if (float d = b.x() * c.y() - b.y() * c.x(); std::abs(d) > EPSILON) { + Vec2f v = lc * b - lb * c; + float r2 = 0.25f * v.squaredNorm() / sqr(d); + if (r2 ) + } +#endif + } + { + // Laplacian smoothing with 0.5 weight, branching point. + const Vec3f &p0 = state.branch.path[state.branch.path.size() - 2].prev_position; + const Vec3f &p1 = state.branch.path.back().prev_position; + float weight = 0; + Vec2f new_pos = Vec2f::Zero(); + for (size_t i = 0; i < state.branch.num_up_trunk; ++i) { + const Element &el = state.branch.up[i].branch->path.front(); + new_pos += el.prev_position.head<2>(); + weight += el.radius; + } + { + const Element &el = state.branch.path[state.branch.path.size() - 2]; + new_pos += el.prev_position.head<2>(); + weight *= 2.f; + } + adjust_position(state.branch.path.back(), 0.5f * state.branch.path.back().prev_position.head<2>() + 0.5f * weight * new_pos); + } + stack.pop_back(); + } else { + // Open another up node of this branch. + stack.push_back({ *state.branch.up[state.idx_up].branch, 0 }); + ++ state.idx_up; + } + } + } + // Also smoothen start of the path. + if (Element &first = root.path.front(); ! first.locked) { + Element &second = root.path[1]; + Vec2f new_pos = 0.75f * first.prev_position.head<2>() + 0.25f * second.prev_position.head<2>(); + if (is_root) + // Let the root of the tree float inside its influence area. + adjust_position(first, new_pos); + else { + // Keep the start of a thin branch inside the trunk. + const Element &trunk = root.down.branch->path.back(); + const float rdif = trunk.radius - root.path.front().radius; + assert(rdif >= 0); + Vec2f vdif = new_pos - trunk.prev_position.head<2>(); + float ldif = vdif.squaredNorm(); + if (ldif > sqr(rdif)) + // Clamp new position. + new_pos = trunk.prev_position.head<2>() + vdif * rdif / sqrt(ldif); + first.position.head<2>() = new_pos; + } + } +} + +void smooth_trees_inside_influence_areas(Forest &forest) +{ + // Parallel for! + for (Tree &tree : forest) + smooth_trees_inside_influence_areas(tree.root(), true); +} + // Test whether two circles, each on its own plane in 3D intersect. // Circles are considered intersecting, if the lowest point on one circle is below the other circle's plane. // Assumption: The two planes are oriented the same way.