diff --git a/src/libslic3r/SLA/SupportTreeUtils.hpp b/src/libslic3r/SLA/SupportTreeUtils.hpp index 1a771e028b..6da28f1fdd 100644 --- a/src/libslic3r/SLA/SupportTreeUtils.hpp +++ b/src/libslic3r/SLA/SupportTreeUtils.hpp @@ -473,12 +473,20 @@ inline long build_ground_connection(SupportTreeBuilder &builder, return ret; } -template +template +constexpr bool IsWideningFn = std::is_invocable_r_v; + +template> > GroundConnection deepsearch_ground_connection( Ex policy, const SupportableMesh &sm, const Junction &j, - double end_radius, + WideningFn &&wideningfn, const Vec3d &init_dir = DOWN) { // Score is the total lenght of the route. Feasible routes will have @@ -488,9 +496,6 @@ GroundConnection deepsearch_ground_connection( const auto sd = sm.cfg.safety_distance(j.r); const auto gndlvl = ground_level(sm); - const double widening = end_radius - j.r; - const double base_r = std::max(sm.cfg.base_radius_mm, end_radius); - const double zelev_gap = sm.cfg.pillar_base_safety_distance_mm + base_r; auto criteria = get_criteria(sm.cfg).stop_score(StopScore); @@ -507,7 +512,7 @@ GroundConnection deepsearch_ground_connection( Vec3d bridge_end = j.pos + bridge_len * n; double full_len = bridge_len + bridge_end.z() - gndlvl; - double bridge_r = j.r + widening * bridge_len / full_len; + double bridge_r = wideningfn(Ball{j.pos, j.r}, n, bridge_len); double brhit_dist = 0.; if (bridge_len > EPSILON) { @@ -522,7 +527,8 @@ GroundConnection deepsearch_ground_connection( ret = brhit_dist; } else { // check if pillar can be placed below - auto gp = Vec3d{bridge_end.x(), bridge_end.y(), gndlvl}; + auto gp = Vec3d{bridge_end.x(), bridge_end.y(), gndlvl}; + double end_radius = wideningfn(Ball{bridge_end, bridge_r}, DOWN, bridge_end.z() - gndlvl); Beam gndbeam {{bridge_end, bridge_r}, {gp, end_radius}}; auto gndhit = beam_mesh_hit(policy, sm.emesh, gndbeam, sd); @@ -533,9 +539,11 @@ GroundConnection deepsearch_ground_connection( if (sm.cfg.object_elevation_mm < EPSILON) { // Dealing with zero elevation mode, to not route pillars // into the gap between the optional pad and the model - double gap = std::sqrt(sm.emesh.squared_distance(gp)); - if (gap < zelev_gap) - ret = full_len - zelev_gap + gap; + double gap = std::sqrt(sm.emesh.squared_distance(gp)); + double base_r = std::max(sm.cfg.base_radius_mm, end_radius); + double max_gap = sm.cfg.pillar_base_safety_distance_mm + base_r; + if (gap < max_gap) + ret = full_len - max_gap + gap; else // success ret = StopScore + EPSILON; } else { @@ -560,8 +568,8 @@ GroundConnection deepsearch_ground_connection( optfn, initvals({plr_init, azm_init, 0.}), // start with a zero bridge bounds({ {PI - sm.cfg.bridge_slope, PI}, // bounds for polar angle - {-PI, PI}, // bounds for azimuth - {0., sm.cfg.max_bridge_length_mm} }) // bounds bridge length + {-PI, PI}, // bounds for azimuth + {0., sm.cfg.max_bridge_length_mm} }) // bounds bridge length ); GroundConnection conn; @@ -572,22 +580,153 @@ GroundConnection deepsearch_ground_connection( Vec3d n = spheric_to_dir(plr, azm); Vec3d bridge_end = j.pos + bridge_len * n; + Vec3d gp{bridge_end.x(), bridge_end.y(), gndlvl}; - double full_len = bridge_len + bridge_end.z() - gndlvl; - double bridge_r = j.r + widening * bridge_len / full_len; - Vec3d gp{bridge_end.x(), bridge_end.y(), gndlvl}; + double bridge_r = wideningfn(Ball{j.pos, j.r}, n, bridge_len); + double down_l = bridge_end.z() - gndlvl; + double end_radius = wideningfn(Ball{bridge_end, bridge_r}, DOWN, down_l); + double base_r = std::max(sm.cfg.base_radius_mm, end_radius); conn.path.emplace_back(j); if (bridge_len > EPSILON) conn.path.emplace_back(Junction{bridge_end, bridge_r}); conn.pillar_base = - Pedestal{gp, sm.cfg.base_height_mm, base_r, end_radius}; + Pedestal{gp, sm.cfg.base_height_mm, base_r, end_radius}; } return conn; } +template +GroundConnection deepsearch_ground_connection(Ex policy, + const SupportableMesh &sm, + const Junction &j, + double end_radius, + const Vec3d &init_dir = DOWN) +{ + double gndlvl = ground_level(sm); + auto wfn = [end_radius, gndlvl](Ball src, Vec3d dir, double len) { + Vec3d dst = src.p + len * dir; + double widening = end_radius - src.R; + double zlen = dst.z() - gndlvl; + double full_len = len + zlen; + double r = src.R + widening * len / full_len; + + return r; + }; + + static_assert(IsWideningFn, "Not a widening function"); + + return deepsearch_ground_connection(policy, sm, j, wfn, init_dir); + +// // Score is the total lenght of the route. Feasible routes will have +// // infinite length (rays not colliding with model), thus the stop score +// // should be a reasonably big number. +// constexpr double StopScore = 1e6; + +// const auto sd = sm.cfg.safety_distance(j.r); +// const auto gndlvl = ground_level(sm); +// const double widening = end_radius - j.r; +// const double base_r = std::max(sm.cfg.base_radius_mm, end_radius); +// const double zelev_gap = sm.cfg.pillar_base_safety_distance_mm + base_r; + +// auto criteria = get_criteria(sm.cfg).stop_score(StopScore); + +// Optimizer solver(criteria); +// solver.seed(0); // enforce deterministic behavior + +// auto optfn = [&](const opt::Input<3> &input) { +// double ret = NaNd; + +// // solver suggests polar, azimuth and bridge length values: +// auto &[plr, azm, bridge_len] = input; + +// Vec3d n = spheric_to_dir(plr, azm); +// Vec3d bridge_end = j.pos + bridge_len * n; + +// double full_len = bridge_len + bridge_end.z() - gndlvl; +// double bridge_r = j.r + widening * bridge_len / full_len; +// double brhit_dist = 0.; + +// if (bridge_len > EPSILON) { +// // beam_mesh_hit with a zero lenght bridge is invalid + +// Beam bridgebeam{Ball{j.pos, j.r}, Ball{bridge_end, bridge_r}}; +// auto brhit = beam_mesh_hit(policy, sm.emesh, bridgebeam, sd); +// brhit_dist = brhit.distance(); +// } + +// if (brhit_dist < bridge_len) { +// ret = brhit_dist; +// } else { +// // check if pillar can be placed below +// auto gp = Vec3d{bridge_end.x(), bridge_end.y(), gndlvl}; + +// Beam gndbeam {{bridge_end, bridge_r}, {gp, end_radius}}; +// auto gndhit = beam_mesh_hit(policy, sm.emesh, gndbeam, sd); + +// if (std::isinf(gndhit.distance())) { +// // Ground route is free with this bridge + +// if (sm.cfg.object_elevation_mm < EPSILON) { +// // Dealing with zero elevation mode, to not route pillars +// // into the gap between the optional pad and the model +// double gap = std::sqrt(sm.emesh.squared_distance(gp)); +// if (gap < zelev_gap) +// ret = full_len - zelev_gap + gap; +// else // success +// ret = StopScore + EPSILON; +// } else { +// // No zero elevation, return success +// ret = StopScore + EPSILON; +// } +// } else { +// // Ground route is not free +// ret = bridge_len + gndhit.distance(); +// } +// } + +// return ret; +// }; + +// auto [plr_init, azm_init] = dir_to_spheric(init_dir); + +// // Saturate the polar angle to max tilt defined in config +// plr_init = std::max(plr_init, PI - sm.cfg.bridge_slope); + +// auto oresult = solver.to_max().optimize( +// optfn, +// initvals({plr_init, azm_init, 0.}), // start with a zero bridge +// bounds({ {PI - sm.cfg.bridge_slope, PI}, // bounds for polar angle +// {-PI, PI}, // bounds for azimuth +// {0., sm.cfg.max_bridge_length_mm} }) // bounds bridge length +// ); + +// GroundConnection conn; + +// if (oresult.score >= StopScore) { +// // search was successful, extract and apply the result +// auto &[plr, azm, bridge_len] = oresult.optimum; + +// Vec3d n = spheric_to_dir(plr, azm); +// Vec3d bridge_end = j.pos + bridge_len * n; + +// double full_len = bridge_len + bridge_end.z() - gndlvl; +// double bridge_r = j.r + widening * bridge_len / full_len; +// Vec3d gp{bridge_end.x(), bridge_end.y(), gndlvl}; + +// conn.path.emplace_back(j); +// if (bridge_len > EPSILON) +// conn.path.emplace_back(Junction{bridge_end, bridge_r}); + +// conn.pillar_base = +// Pedestal{gp, sm.cfg.base_height_mm, base_r, end_radius}; +// } + +// return conn; +} + template bool optimize_anchor_placement(Ex policy, const SupportableMesh &sm, @@ -671,6 +810,92 @@ std::optional calculate_anchor_placement(Ex policy, return anchor; } +inline bool is_outside_support_cone(const Vec3f &supp, + const Vec3f &pt, + float angle) +{ + using namespace Slic3r; + + Vec3d D = (pt - supp).cast(); + double dot_sq = -D.z() * std::abs(-D.z()); + + return dot_sq < + D.squaredNorm() * std::cos(angle) * std::abs(std::cos(angle)); +} + +inline // TODO: should be in a cpp +std::optional find_merge_pt(const Vec3f &A, + const Vec3f &B, + float critical_angle) +{ + // The idea is that A and B both have their support cones. But searching + // for the intersection of these support cones is difficult and its enough + // to reduce this problem to 2D and search for the intersection of two + // rays that merge somewhere between A and B. The 2D plane is a vertical + // slice of the 3D scene where the 2D Y axis is equal to the 3D Z axis and + // the 2D X axis is determined by the XY direction of the AB vector. + // + // Z^ + // | A * + // | . . B * + // | . . . . + // | . . . . + // | . x . + // -------------------> XY + + // Determine the transformation matrix for the 2D projection: + Vec3f diff = {B.x() - A.x(), B.y() - A.y(), 0.f}; + Vec3f dir = diff.normalized(); // TODO: avoid normalization + + Eigen::Matrix tr2D; + tr2D.row(0) = Vec3f{dir.x(), dir.y(), dir.z()}; + tr2D.row(1) = Vec3f{0.f, 0.f, 1.f}; + + // Transform the 2 vectors A and B into 2D vector 'a' and 'b'. Here we can + // omit 'a', pretend that its the origin and use BA as the vector b. + Vec2f b = tr2D * (B - A); + + // Get the square sine of the ray emanating from 'a' towards 'b'. This ray might + // exceed the allowed angle but that is corrected subsequently. + // The sign of the original sine is also needed, hence b.y is multiplied by + // abs(b.y) + float b_sqn = b.squaredNorm(); + float sin2sig_a = b_sqn > EPSILON ? (b.y() * std::abs(b.y())) / b_sqn : 0.f; + + // sine2 from 'b' to 'a' is the opposite of sine2 from a to b + float sin2sig_b = -sin2sig_a; + + // Derive the allowed angles from the given critical angle. + // critical_angle is measured from the horizontal X axis. + // The rays need to go downwards which corresponds to negative angles + + float sincrit = std::sin(critical_angle); // sine of the critical angle + float sin2crit = -sincrit * sincrit; // signed sine squared + sin2sig_a = std::min(sin2sig_a, sin2crit); // Do the angle saturation of both rays + sin2sig_b = std::min(sin2sig_b, sin2crit); // + float sin2_a = std::abs(sin2sig_a); // Get cosine squared values + float sin2_b = std::abs(sin2sig_b); + float cos2_a = 1.f - sin2_a; + float cos2_b = 1.f - sin2_b; + + // Derive the new direction vectors. This is by square rooting the sin2 + // and cos2 values and restoring the original signs + Vec2f Da = {std::copysign(std::sqrt(cos2_a), b.x()), std::copysign(std::sqrt(sin2_a), sin2sig_a)}; + Vec2f Db = {-std::copysign(std::sqrt(cos2_b), b.x()), std::copysign(std::sqrt(sin2_b), sin2sig_b)}; + + // Determine where two rays ([0, 0], Da), (b, Db) intersect. + // Based on + // https://stackoverflow.com/questions/27459080/given-two-points-and-two-direction-vectors-find-the-point-where-they-intersect + // One ray is emanating from (0, 0) so the formula is simplified + double t1 = (Db.y() * b.x() - b.y() * Db.x()) / + (Da.x() * Db.y() - Da.y() * Db.x()); + + Vec2f mp = t1 * Da; + Vec3f Mp = A + tr2D.transpose() * mp; + + return t1 >= 0.f ? Mp : Vec3f{}; +} + }} // namespace Slic3r::sla #endif // SLASUPPORTTREEUTILS_H