Widening improvements in SupportTreeUtils

Fix failing tests after introducing wideningfn into ground route search
This commit is contained in:
tamasmeszaros 2022-11-28 11:01:49 +01:00
parent 2cd6a20254
commit 9cdd6738ae

View file

@ -473,12 +473,20 @@ inline long build_ground_connection(SupportTreeBuilder &builder,
return ret;
}
template<class Ex>
template<class Fn>
constexpr bool IsWideningFn = std::is_invocable_r_v</*retval*/ double,
Fn,
Ball /*source*/,
Vec3d /*dir*/,
double /*length*/>;
template<class Ex, class WideningFn,
class = std::enable_if_t<IsWideningFn<WideningFn>> >
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<class Ex>
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<decltype(wfn)>, "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<opt::AlgNLoptMLSL> 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<class Ex>
bool optimize_anchor_placement(Ex policy,
const SupportableMesh &sm,
@ -671,6 +810,92 @@ std::optional<Anchor> 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>();
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<Vec3f> 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<float, 2, 3> 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