diff --git a/editor/editorAlignment.cpp b/editor/editorAlignment.cpp index b0b501c1..2907d3af 100644 --- a/editor/editorAlignment.cpp +++ b/editor/editorAlignment.cpp @@ -299,6 +299,51 @@ recommendation recommend(double const Speed, double const Radius, limits const & namespace { +// the tangents of a curve of the radius with transitions of the lengths, from the vertex back along the straight in and on along the +// one out +void curve_tangents(transition_shape const Shape, double const Lengthin, double const Lengthout, double const Radius, glm::dvec2 const &Directionin, glm::dvec2 const &Directionout, int const Side, + double &Tangentin, double &Tangentout) +{ + auto const endin{transition_point(Shape, Lengthin, Radius, 1.0)}; + auto const endout{transition_point(Shape, Lengthout, Radius, 1.0)}; + auto const a{Radius + endin.y - Radius * (1.0 - std::cos(endin.angle))}; + auto const b{Radius + endout.y - Radius * (1.0 - std::cos(endout.angle))}; + auto const normalin{perpendicular(Directionin) * static_cast(Side)}; + auto const normalout{perpendicular(Directionout) * static_cast(Side)}; + auto const determinant{cross(normalin, normalout)}; + glm::dvec2 const offset{(a * normalout.y - normalin.y * b) / determinant, (normalin.x * b - a * normalout.x) / determinant}; + auto const curvestart{offset - normalin * a - Directionin * (endin.x - Radius * std::sin(endin.angle))}; + auto const curveend{offset - normalout * b + Directionout * (endout.x - Radius * std::sin(endout.angle))}; + Tangentin = std::max(0.0, -glm::dot(curvestart, Directionin)); + Tangentout = std::max(0.0, glm::dot(curveend, Directionout)); +} + +} + +// each of the two transitions turns half of the deflection, ending at the curvature of the radius +double spirals_length(transition_shape const Shape, double const Deflection, double const Radius) +{ + auto const half{std::clamp(Deflection * 0.5, 0.0, 0.5 * kPi - 1e-6)}; + return Shape == transition_shape::clothoid ? 2.0 * Radius * half : 2.0 * Radius * std::tan(half); +} + +bool meets_norms(curve_report const &Curve, design const &Design) +{ + auto const &norms{Design.norms}; + auto sharpest{Curve.radius}; + for (auto const radius : Curve.radii) + sharpest = std::min(sharpest, radius); + if (Curve.unbalanced > norms.unbalanced + 1e-6 || Curve.cant > norms.cant_max + 1e-6) + return false; + if (false == Curve.spirals && sharpest < Curve.recommended.radius_min - 1e-6) + return false; + auto const required{transition_length(Design.speed, sharpest, Curve.cant, norms)}; + return Curve.transition_in + 1e-6 >= required && Curve.transition_out + 1e-6 >= required; +} + +namespace +{ + result fit_between(design const &Design) { result r; @@ -403,6 +448,7 @@ result fit_between(design const &Design) double deflection{0.0}; double cant{0.0}; bool compound{false}; + bool spirals{false}; std::vector arcs; std::function tangents; }; @@ -456,6 +502,19 @@ result fit_between(design const &Design) fit.tangents(fit.radius, fit.tangent_in, fit.tangent_out); continue; } + if (vertex.spirals) + { + // the lengths follow from the radius, so a smaller one makes the whole curve smaller in proportion + fit.spirals = true; + fit.deflection = deflection; + fit.transition_in = fit.transition_out = spirals_length(Design.shape, deflection, fit.radius); + fit.tangents = [&, directionin, directionout, side](double const Radius, double &Tangentin, double &Tangentout) { + auto const length{spirals_length(Design.shape, deflection, Radius)}; + curve_tangents(Design.shape, length, length, Radius, directionin, directionout, side, Tangentin, Tangentout); + }; + fit.tangents(fit.radius, fit.tangent_in, fit.tangent_out); + continue; + } auto const angles = [&](double const Factor) { return transition_point(Design.shape, fit.transition_in * Factor, fit.radius, 1.0).angle + transition_point(Design.shape, fit.transition_out * Factor, fit.radius, 1.0).angle; }; @@ -474,18 +533,7 @@ result fit_between(design const &Design) } fit.deflection = deflection; fit.tangents = [&, k, directionin, directionout, side](double const Radius, double &Tangentin, double &Tangentout) { - auto const endin{transition_point(Design.shape, fits[k].transition_in, Radius, 1.0)}; - auto const endout{transition_point(Design.shape, fits[k].transition_out, Radius, 1.0)}; - auto const a{Radius + endin.y - Radius * (1.0 - std::cos(endin.angle))}; - auto const b{Radius + endout.y - Radius * (1.0 - std::cos(endout.angle))}; - auto const normalin{perpendicular(directionin) * static_cast(side)}; - auto const normalout{perpendicular(directionout) * static_cast(side)}; - auto const determinant{cross(normalin, normalout)}; - glm::dvec2 const offset{(a * normalout.y - normalin.y * b) / determinant, (normalin.x * b - a * normalout.x) / determinant}; - auto const curvestart{offset - normalin * a - directionin * (endin.x - Radius * std::sin(endin.angle))}; - auto const curveend{offset - normalout * b + directionout * (endout.x - Radius * std::sin(endout.angle))}; - Tangentin = std::max(0.0, -glm::dot(curvestart, directionin)); - Tangentout = std::max(0.0, glm::dot(curveend, directionout)); + curve_tangents(Design.shape, fits[k].transition_in, fits[k].transition_out, Radius, directionin, directionout, side, Tangentin, Tangentout); }; fit.tangents(fit.radius, fit.tangent_in, fit.tangent_out); } @@ -526,7 +574,7 @@ result fit_between(design const &Design) } low = high * (1.0 + 1e-9); high = fit.radius; - if (false == fit.compound && fits_budget(low)) + if (false == fit.compound && false == fit.spirals && fits_budget(low)) { for (int i = 0; i < 60; ++i) { @@ -549,8 +597,8 @@ result fit_between(design const &Design) { auto &fit{fits[k]}; auto const &vertex{Design.vertices[k]}; - fit.cant = std::max(0.0, vertex.cant); - if (fit.compound || fit.deflection < 1e-7 || fit.transition_in + fit.transition_out <= 0.0) + fit.cant = fit.spirals ? 0.0 : std::max(0.0, vertex.cant); + if (fit.compound || fit.spirals || fit.deflection < 1e-7 || fit.transition_in + fit.transition_out <= 0.0) continue; auto const angles = [&](double const Factor) { return transition_point(Design.shape, fit.transition_in * Factor, fit.radius, 1.0).angle + transition_point(Design.shape, fit.transition_out * Factor, fit.radius, 1.0).angle; @@ -672,6 +720,7 @@ result fit_between(design const &Design) report.tangent_in = glm::distance(position, curvestart); report.tangent_out = glm::distance(position, cursor); report.unbalanced = unbalanced_acceleration(Design.speed, sharpest, std::max(0.0, vertex.cant), Design.norms); + report.cant = std::max(0.0, vertex.cant); report.recommended = recommend(Design.speed, sharpest, Design.norms); r.curves.push_back(report); if (Design.norms.checked && sharpest < report.recommended.radius_min - 1e-6) @@ -699,7 +748,7 @@ result fit_between(design const &Design) auto const curvestart{centre - normalin * a - directionin * middlein}; auto const curveend{centre - normalout * b + directionout * middleout}; auto const arcangle{deflection - endin.angle - endout.angle}; - if (arcangle < -1e-9) + if (arcangle < -1e-9 && false == fit.spirals) { r.errors.emplace_back(format(STR_C("Vertex %zu: transition curves are too long for the deflection angle"), k + 1)); return r; @@ -759,6 +808,7 @@ result fit_between(design const &Design) report.tangent_in = glm::distance(position, curvestart); report.tangent_out = glm::distance(position, curveend); report.unbalanced = unbalanced_acceleration(Design.speed, radius, curve.cant, Design.norms); + report.cant = curve.cant; auto const speed{Design.speed / 3.6}; if (lengthin > 0.0) { @@ -771,17 +821,27 @@ result fit_between(design const &Design) report.jerk_out = std::max(0.0, report.unbalanced) * speed / lengthout; } report.recommended = recommend(Design.speed, radius, Design.norms); + report.spirals = fit.spirals; r.curves.push_back(report); auto const &norms{Design.norms}; if (false == norms.checked) continue; - if (radius < report.recommended.radius_min - 1e-6) - r.warnings.emplace_back(format(STR_C("Vertex %zu: R %.0f m is below the minimum %.0f m for %.0f km/h"), k + 1, radius, report.recommended.radius_min, Design.speed)); + if (fit.spirals) + { + // with no cant all the acceleration across is unbalanced + if (report.unbalanced > norms.unbalanced + 1e-6) + r.warnings.emplace_back(format(STR_C("Vertex %zu: the transition curves meet at R %.0f m, with no cant at least %.0f m for %.0f km/h"), k + 1, radius, Design.speed * Design.speed / (12.96 * norms.unbalanced), Design.speed)); + } + else + { + if (radius < report.recommended.radius_min - 1e-6) + r.warnings.emplace_back(format(STR_C("Vertex %zu: R %.0f m is below the minimum %.0f m for %.0f km/h"), k + 1, radius, report.recommended.radius_min, Design.speed)); + if (report.unbalanced > norms.unbalanced + 1e-6) + r.warnings.emplace_back(format(STR_C("Vertex %zu: unbalanced acceleration %.2f m/s2 exceeds %.2f m/s2 (cant at least %.0f mm)"), k + 1, report.unbalanced, norms.unbalanced, report.recommended.cant_min)); + } if (curve.cant > norms.cant_max + 1e-6) r.warnings.emplace_back(format(STR_C("Vertex %zu: cant %.0f mm exceeds %.0f mm"), k + 1, curve.cant, norms.cant_max)); - if (report.unbalanced > norms.unbalanced + 1e-6) - r.warnings.emplace_back(format(STR_C("Vertex %zu: unbalanced acceleration %.2f m/s2 exceeds %.2f m/s2 (cant at least %.0f mm)"), k + 1, report.unbalanced, norms.unbalanced, report.recommended.cant_min)); auto const required{transition_length(Design.speed, radius, curve.cant, norms)}; if (lengthin + 1e-6 < required || lengthout + 1e-6 < required) r.warnings.emplace_back(format(STR_C("Vertex %zu: transition curves shorter than the required %.0f m"), k + 1, required)); @@ -906,7 +966,7 @@ result compute_recommended(design &Design, bool const Transitions) auto slowest{speed}; for (auto const &curve : r.curves) { - if (curve.vertex < 0 || curve.vertex >= static_cast(attempt.vertices.size())) + if (curve.vertex < 0 || curve.vertex >= static_cast(attempt.vertices.size()) || attempt.vertices[curve.vertex].spirals) continue; auto const allowed{std::max(1.0, std::min(speed, speed_limit(curve.radius, attempt.norms)))}; auto const recommended{recommend(allowed, curve.radius, attempt.norms)}; @@ -2401,4 +2461,598 @@ std::vector pieces(result const &Result, design const &Design, std return result; } +namespace +{ + +// the vertices where fit_between puts them, the start before them and the end after; empty where it can't +std::vector polygon_of(design const &Design) +{ + if (glm::length(Design.start_direction) < 1e-9 || glm::length(Design.end_direction) < 1e-9) + return {}; + auto const count{Design.vertices.size()}; + auto const start{plan_of(Design.start)}; + auto const end{plan_of(Design.end)}; + std::vector polygon{start}; + if (count == 1) + { + glm::dvec2 corner; + if (false == tangent_intersection(Design, corner)) + return {}; + polygon.push_back(corner); + } + else if (count >= 2) + { + polygon.push_back(start + glm::normalize(Design.start_direction) * Design.vertices.front().offset); + for (std::size_t i = 1; i + 1 < count; ++i) + polygon.push_back(Design.vertices[i].position); + polygon.push_back(end - glm::normalize(Design.end_direction) * Design.vertices.back().offset); + } + polygon.push_back(end); + return polygon; +} + +// the cant the transitions of the length can ramp at the speed, in steps of 5 mm +double ramped_cant(double const Cant, double const Length, double const Speed, limits const &Norms) +{ + auto const cant{std::min({Cant, Length * 1000.0 / (Norms.ramp_factor * Speed), Norms.cant_rate * 3.6 * Length / Speed})}; + return cant < 5.0 ? 0.0 : std::floor(cant / 5.0) * 5.0; +} + +// the transitions and the cant of the norms for the radii of the vertex, at the speed its sharpest arc allows +void recommended_curve(vertex &Vertex, double const Speed, limits const &Norms) +{ + auto sharpest{std::max(1.0, Vertex.radius)}; + if (Vertex.compound) + for (auto const &arc : Vertex.arcs) + sharpest = std::min(sharpest, std::max(1.0, arc.radius)); + auto const allowed{std::max(1.0, std::min(Speed, speed_limit(sharpest, Norms)))}; + auto const recommended{recommend(allowed, sharpest, Norms)}; + Vertex.transition_in = Vertex.transition_out = recommended.transition; + Vertex.cant = ramped_cant(recommended.cant, recommended.transition, allowed, Norms); + recommend_between(Vertex, Speed, Norms); +} + +// a curve a bend can be eased with, and how far it reaches along the straights either side of it +struct corner_shape +{ + bool spirals{false}; + double radius{0.0}; + double transition{0.0}; + double cant{0.0}; + double reach{0.0}; +}; + +} // namespace + +result ease_corners(design &Design, std::vector const &Corners, corner_radii const &Radii) +{ + auto const count{Design.vertices.size()}; + auto const corner = [&](std::size_t const K) { return K < Corners.size() && Corners[K] != 0; }; + // the curves of the vertices left as they are reach along the straights as far as they do with the bends kept out of their way + std::vector> kept(count, {0.0, 0.0}); + { + auto bare{Design}; + for (std::size_t k = 0; k < count; ++k) + if (corner(k)) + bare.vertices[k] = {bare.vertices[k].position, bare.vertices[k].offset, 1.0}; + for (auto const &curve : compute(bare).curves) + if (curve.vertex >= 0 && curve.vertex < static_cast(count)) + kept[curve.vertex] = {curve.tangent_in, curve.tangent_out}; + } + auto const polygon{polygon_of(Design)}; + if (polygon.size() != count + 2) + return compute(Design); + auto const &norms{Design.norms}; + auto const speed{std::max(1.0, Design.speed)}; + auto const low{std::max(10.0, std::min(Radii.low, Radii.high))}; + auto const high{std::max(low, Radii.high)}; + + // the shapes each bend can take, the best first + std::vector> ladders(count); + for (std::size_t k = 0; k < count; ++k) + { + if (false == corner(k)) + continue; + auto const in{polygon[k + 1] - polygon[k]}; + auto const out{polygon[k + 2] - polygon[k + 1]}; + if (glm::length(in) < 1e-6 || glm::length(out) < 1e-6) + continue; + auto const directionin{glm::normalize(in)}; + auto const directionout{glm::normalize(out)}; + auto const turn{cross(directionin, directionout)}; + auto const deflection{std::atan2(std::abs(turn), glm::dot(directionin, directionout))}; + int const side{turn > 0.0 ? 1 : -1}; + if (deflection < 1e-7 || deflection > kPi - 1e-3) + continue; + auto const reach = [&](double const Length, double const Radius) { + double tangentin; + double tangentout; + curve_tangents(Design.shape, Length, Length, Radius, directionin, directionout, side, tangentin, tangentout); + return std::max(tangentin, tangentout); + }; + auto &ladder{ladders[k]}; + // two transitions meeting with no cant, long enough for the acceleration across to change gently at the speed: a bend too + // slight for them at the highest radius takes a larger one + auto const long_enough = [&](double const Radius) { return spirals_length(Design.shape, deflection, Radius) >= transition_length(speed, Radius, 0.0, norms) - 1e-6; }; + auto top{high}; + if (false == long_enough(top)) + { + auto below{top}; + auto above{top * 2.0}; + while (false == long_enough(above) && above < 1e7) + { + below = above; + above *= 2.0; + } + for (int i = 0; i < 40; ++i) + { + auto const middle{(below + above) * 0.5}; + (long_enough(middle) ? above : below) = middle; + } + top = std::ceil(above / 10.0) * 10.0; + } + constexpr int kSpirals{30}; + for (int i = 0; i <= kSpirals; ++i) + { + auto const radius{std::round(top * std::pow(low / top, static_cast(i) / kSpirals))}; + if (false == long_enough(radius)) + break; + auto const length{spirals_length(Design.shape, deflection, radius)}; + ladder.push_back({true, radius, length, 0.0, reach(length, radius)}); + if (top <= low) + break; + } + // then arcs with the transitions and the cant of the norms, as large as fit, to any radius: the transitions shortened the way + // fit_between does it, to leave the arc its least length + auto const keep_arc{std::max(norms.arc_min, norms.arc_min_factor * speed)}; + constexpr int kArcs{60}; + for (int i = 0; i <= kArcs; ++i) + { + auto const radius{std::round(high * std::pow(10.0 / high, static_cast(i) / kArcs))}; + auto const allowed{std::max(1.0, std::min(speed, speed_limit(radius, norms)))}; + auto const recommended{recommend(allowed, radius, norms)}; + auto const turn_of = [&](double const Length) { return 2.0 * transition_point(Design.shape, Length, radius, 1.0).angle; }; + auto const least{std::min(radius * deflection, keep_arc)}; + auto length{recommended.transition}; + if (length > 0.0 && radius * (deflection - turn_of(length)) < least - 1e-6) + { + double shorter{0.0}; + double longer{length}; + for (int j = 0; j < 50; ++j) + { + auto const middle{(shorter + longer) * 0.5}; + (radius * (deflection - turn_of(middle)) < least ? longer : shorter) = middle; + } + length = shorter < 1.0 ? 0.0 : shorter; + } + auto const cant{length > 0.0 ? ramped_cant(recommended.cant, length, allowed, norms) : 0.0}; + ladder.push_back({false, radius, length, cant, reach(length, radius)}); + } + } + + // the straight between two curves is shared: each of them gets what it wants of it where they both fit, else the one wanting less + // than a half gets it, the other the rest, or each a half. Then each bend takes the best shape the room its neighbours leave fits + auto const between{norms.tangent_min_time * speed / 3.6}; + auto const leg = [&](std::size_t const Leg) { + return glm::distance(polygon[Leg], polygon[Leg + 1]) + (Leg == 0 ? std::max(0.0, Design.start_reserve) : 0.0) + (Leg == count ? std::max(0.0, Design.end_reserve) : 0.0) - + (Leg > 0 && Leg < count ? between : 0.0); + }; + std::vector picks(count, 0); + auto const reach_of = [&](std::size_t const K, bool const Out) { + if (false == ladders[K].empty()) + return ladders[K][picks[K]].reach; + return corner(K) ? 0.0 : kept[K][Out ? 1 : 0]; + }; + auto const best_within = [&](std::size_t const K, double const Room) { + auto const &ladder{ladders[K]}; + for (std::size_t i = 0; i < ladder.size(); ++i) + if (ladder[i].reach <= Room + 1e-6) + return i; + return ladder.size() - 1; + }; + std::vector roomin(count, std::numeric_limits::max()); + std::vector roomout(count, std::numeric_limits::max()); + for (std::size_t j = 0; j <= count; ++j) + { + auto const room{std::max(0.0, leg(j))}; + auto const before{j > 0 ? reach_of(j - 1, true) : 0.0}; + auto const after{j < count ? reach_of(j, false) : 0.0}; + auto give_before{room - after}; + auto give_after{room - before}; + if (before + after > room) + { + auto const half{room * 0.5}; + give_before = after <= half ? room - after : before <= half ? before : half; + give_after = before <= half ? room - before : after <= half ? after : half; + } + if (j > 0) + roomout[j - 1] = give_before; + if (j < count) + roomin[j] = give_after; + } + for (std::size_t k = 0; k < count; ++k) + if (false == ladders[k].empty()) + picks[k] = best_within(k, std::min(roomin[k], roomout[k])); + for (int pass = 0; pass < 3; ++pass) + for (std::size_t k = 0; k < count; ++k) + { + if (ladders[k].empty()) + continue; + auto const in{leg(k) - (k > 0 ? reach_of(k - 1, true) : 0.0)}; + auto const out{leg(k + 1) - (k + 1 < count ? reach_of(k + 1, false) : 0.0)}; + picks[k] = best_within(k, std::min(in, out)); + } + + for (std::size_t k = 0; k < count; ++k) + { + if (ladders[k].empty()) + continue; + auto const &shape{ladders[k][picks[k]]}; + auto &vertex{Design.vertices[k]}; + vertex.spirals = shape.spirals; + vertex.radius = shape.radius; + vertex.transition_in = vertex.transition_out = shape.transition; + vertex.cant = shape.cant; + vertex.reverse_turn = false; + vertex.compound = false; + vertex.share = 1.0; + vertex.arcs.clear(); + } + return compute(Design); +} + +namespace +{ + +// a curve a single radius follows the points about as closely as this needs no more of them +constexpr double kStretchClose{0.08}; + +// the stretch as the design lays it, how far the points are from it +std::optional stretch_of(design const &Design, std::vector const &Points) +{ + auto const laid{compute(Design)}; + if (false == laid.valid) + return std::nullopt; + trace_stretch stretch; + stretch.valid = true; + stretch.vertices = Design.vertices; + for (std::size_t k = 0; k < stretch.vertices.size() && k < laid.vertices.size(); ++k) + stretch.vertices[k].position = laid.vertices[k]; + double sum{0.0}; + for (auto const &projection : project(laid, Points)) + { + stretch.offsets.push_back(projection.offset); + sum += projection.offset * projection.offset; + stretch.worst = std::max(stretch.worst, std::abs(projection.offset)); + } + stretch.rms = Points.empty() ? 0.0 : std::sqrt(sum / static_cast(Points.size())); + return stretch; +} + +// how far the points are from the line the design lays, a point clicked by mistake far off pulling it less than the square of its +// departure would +bool stretch_residuals(design Trial, std::vector const &Points, std::vector &Residuals) +{ + Trial.norms.checked = false; + auto const laid{compute(Trial)}; + if (false == laid.valid) + return false; + auto const line{sample_plan(laid, 0.0, laid.length, std::clamp(laid.length / 300.0, 0.5, 4.0))}; + Residuals.clear(); + for (auto const &point : Points) + { + auto const offset{project_on(line, point).offset}; + auto const departure{std::abs(offset)}; + Residuals.push_back(departure <= kTraceRobust ? offset : std::copysign(std::sqrt(2.0 * kTraceRobust * departure - kTraceRobust * kTraceRobust), offset)); + } + return true; +} + +// a curve of one radius with transitions of the length: the cant of the norms for it, as much as the transitions ramp +void plain_curve(vertex &Vertex, double const Radius, double const Transition, double const Speed, limits const &Norms) +{ + Vertex.radius = Radius; + Vertex.compound = false; + Vertex.share = 1.0; + Vertex.arcs.clear(); + Vertex.spirals = false; + Vertex.transition_in = Vertex.transition_out = Transition; + auto const allowed{std::max(1.0, std::min(Speed, speed_limit(Radius, Norms)))}; + Vertex.cant = Transition > 0.0 ? ramped_cant(recommend(allowed, Radius, Norms).cant, Transition, allowed, Norms) : 0.0; +} + +// the radius and the length of the transitions of the curve of the vertex following the points best, the rest of the design kept: the +// transitions of a line on the photo are what they are, the norms may have asked for others +void fit_curve(design &Design, std::size_t const Vertex, std::vector const &Points) +{ + auto const speed{std::max(1.0, Design.speed)}; + auto const shaped = [&](std::vector const &Parameters, design &Trial) { + Trial = Design; + auto const radius{std::exp(std::clamp(Parameters[0], std::log(kTraceRadiusMin), std::log(kTraceRadiusMax)))}; + plain_curve(Trial.vertices[Vertex], radius, std::min(Parameters[1] * Parameters[1], 2.0 * radius), speed, Trial.norms); + }; + residuals_of const residuals{[&](std::vector const &Parameters, std::vector &Residuals) { + design trial; + shaped(Parameters, trial); + return stretch_residuals(trial, Points, Residuals); + }}; + // from the transitions it has and from those of the norms; the square root of the length, as a length of nought would hold the + // step at nought + auto const radius{std::max(kTraceRadiusMin, Design.vertices[Vertex].radius)}; + auto const allowed{std::max(1.0, std::min(speed, speed_limit(radius, Design.norms)))}; + std::vector best; + auto lowest{std::numeric_limits::max()}; + for (auto const length : {std::max(Design.vertices[Vertex].transition_in, 10.0), std::max(recommend(allowed, radius, Design.norms).transition, 10.0)}) + { + std::vector parameters{std::log(radius), std::sqrt(length)}; + if (auto const cost{least_squares(parameters, residuals)}; cost < lowest) + { + lowest = cost; + best = parameters; + } + } + if (best.empty()) + return; + design fitted; + shaped(best, fitted); + auto &curve{fitted.vertices[Vertex]}; + plain_curve(curve, std::max(kTraceRadiusMin, std::round(curve.radius)), std::round(curve.transition_in), speed, fitted.norms); + Design = fitted; +} + +// the radii and the shares of the compound curve of the vertex together with the length of its transitions, following the points +// best, from the shape it has +void fit_compound(design &Design, std::size_t const Vertex, std::vector const &Points) +{ + auto const base{Design.vertices[Vertex]}; + if (false == base.compound || base.arcs.empty()) + return; + auto const speed{std::max(1.0, Design.speed)}; + int const radii{static_cast(base.arcs.size()) + 1}; + // as traced_shape takes them, the square root of the length of the transitions last + std::vector start{std::log(std::max(kTraceRadiusMin, base.radius))}; + for (auto const &arc : base.arcs) + start.push_back(std::log(std::max(kTraceRadiusMin, arc.radius))); + for (auto const &arc : base.arcs) + start.push_back(std::log(std::max(1e-3, arc.share) / std::max(1e-3, base.share))); + auto const shaped = [&](std::vector const &Parameters, design &Trial) { + Trial = Design; + auto &curve{Trial.vertices[Vertex]}; + curve = traced_shape(Design, base, Parameters, radii); + curve.transition_in = curve.transition_out = std::min(Parameters.back() * Parameters.back(), 2.0 * curve.radius); + }; + residuals_of const residuals{[&](std::vector const &Parameters, std::vector &Residuals) { + design trial; + shaped(Parameters, trial); + return stretch_residuals(trial, Points, Residuals); + }}; + std::vector best; + auto lowest{std::numeric_limits::max()}; + for (auto const length : {base.transition_in, 0.5 * base.transition_in}) + { + auto parameters{start}; + parameters.push_back(std::sqrt(std::max(length, 10.0))); + if (auto const cost{least_squares(parameters, residuals)}; cost < lowest) + { + lowest = cost; + best = parameters; + } + } + if (best.empty()) + return; + design fitted; + shaped(best, fitted); + auto &curve{fitted.vertices[Vertex]}; + curve.radius = std::max(kTraceRadiusMin, std::round(curve.radius)); + auto sharpest{curve.radius}; + for (auto &arc : curve.arcs) + sharpest = std::min(sharpest, arc.radius = std::max(kTraceRadiusMin, std::round(arc.radius))); + curve.transition_in = curve.transition_out = std::round(curve.transition_in); + auto const allowed{std::max(1.0, std::min(speed, speed_limit(sharpest, fitted.norms)))}; + curve.cant = curve.transition_in > 0.0 ? ramped_cant(recommend(allowed, sharpest, fitted.norms).cant, curve.transition_in, allowed, fitted.norms) : 0.0; + recommend_between(curve, speed, fitted.norms); + Design = fitted; +} + +} // namespace + +trace_stretch fit_trace_stretch(design const &Base, std::vector const &Points, int const Radii, bool const Straights) +{ + trace_stretch chosen; + if (Points.empty() || glm::length(Base.start_direction) < 1e-9 || glm::length(Base.end_direction) < 1e-9) + return chosen; + auto const speed{std::max(1.0, Base.speed)}; + auto const start{plan_of(Base.start)}; + auto const startdirection{glm::normalize(Base.start_direction)}; + auto const end{plan_of(Base.end)}; + auto const enddirection{glm::normalize(Base.end_direction)}; + // a shape with more curves is taken where it follows the points much closer than the simpler one, which doesn't keep to them + auto const take = [&](std::optional const &Candidate) { + if (false == Candidate.has_value()) + return; + if (false == chosen.valid || (chosen.rms > kStretchClose && Candidate->rms < 0.6 * chosen.rms)) + chosen = *Candidate; + }; + + // one straight, where the two are one + if (auto const offset{end - start}; glm::dot(startdirection, enddirection) > std::cos(2e-3) && std::abs(cross(startdirection, offset)) < 0.25 && glm::dot(startdirection, offset) > 0.0) + { + auto straight{Base}; + straight.vertices.clear(); + auto const along{glm::dot(startdirection, offset)}; + straight.end = {Base.start.x + startdirection.x * along, Base.end.y, Base.start.z + startdirection.y * along}; + straight.end_direction = startdirection; + take(stretch_of(straight, Points)); + } + + // one curve where the straights meet: the radius of a plain arc first, then the radius and the transitions together; the radii of + // a compound curve where one radius doesn't keep to the points, or where they're asked for + std::optional one; + if (glm::dvec2 meeting; tangent_intersection(Base, meeting)) + { + auto single{Base}; + single.vertices.assign(1, vertex{}); + auto &curve{single.vertices.front()}; + curve.transition_in = curve.transition_out = curve.cant = 0.0; + if (auto const arc{fit_trace(single, 0, Points, 1)}; arc.valid) + { + curve = arc.fitted; + fit_curve(single, 0, Points); + one = stretch_of(single, Points); + if (one.has_value() && Radii != 1 && (Radii > 1 || one->rms > kStretchClose)) + if (auto const compound{fit_trace(single, 0, Points, Radii)}; compound.valid && compound.fitted.compound) + { + auto several{single}; + several.vertices.front() = compound.fitted; + fit_compound(several, 0, Points); + if (auto const candidate{stretch_of(several, Points)}; candidate.has_value() && (Radii > 1 || candidate->rms < 0.6 * one->rms)) + one = candidate; + } + } + take(one); + } + + // two curves with a straight between them, turning the same way or not: the straight is where the points keep to a line between + // the curves. The parameters: the direction of the straight, its distance from the start across it, the logarithms of the radii, + // first with the transitions of the single curve, then with the square roots of the lengths of the transitions of each curve too + constexpr std::size_t kParameters{4}; + if (Straights && Points.size() > kParameters && (false == chosen.valid || chosen.rms > kStretchClose)) + { + auto const transition{one.has_value() && false == one->vertices.empty() ? one->vertices.front().transition_in : 0.0}; + auto const shaped = [&](std::vector const &Parameters, design &Trial) { + glm::dvec2 const direction{std::cos(Parameters[0]), std::sin(Parameters[0])}; + glm::dvec2 const across{-direction.y, direction.x}; + auto const meet = [&](glm::dvec2 const &From, glm::dvec2 const &Along, double &Distance) { + auto const denominator{glm::dot(across, Along)}; + if (std::abs(denominator) < 1e-9) + return false; + Distance = (Parameters[1] - glm::dot(across, From - start)) / denominator; + return Distance > 0.0; + }; + double first; + double second; + if (false == meet(start, startdirection, first) || false == meet(end, -enddirection, second)) + return false; + if (glm::dot((end - enddirection * second) - (start + startdirection * first), direction) <= 0.0) + return false; + Trial = Base; + Trial.vertices.assign(2, vertex{}); + Trial.vertices[0].offset = first; + Trial.vertices[1].offset = second; + for (std::size_t i = 0; i < 2; ++i) + { + auto const radius{std::exp(std::clamp(Parameters[2 + i], std::log(kTraceRadiusMin), std::log(kTraceRadiusMax)))}; + auto const length{Parameters.size() > kParameters ? Parameters[kParameters + i] * Parameters[kParameters + i] : transition}; + plain_curve(Trial.vertices[i], radius, std::min(length, 2.0 * radius), speed, Trial.norms); + } + return true; + }; + residuals_of const residuals{[&](std::vector const &Parameters, std::vector &Residuals) { + design trial; + return shaped(Parameters, trial) && stretch_residuals(trial, Points, Residuals); + }}; + // started from the lines through neighbouring points spread over the stretch, short of the two on the straight at its end + auto const guess{std::log(one.has_value() && false == one->vertices.empty() ? std::max(kTraceRadiusMin, one->vertices.front().radius) : 1000.0)}; + auto const travel{end - start}; + std::vector best; + auto lowest{std::numeric_limits::max()}; + auto const pairs{Points.size() - 2}; + auto const starts{std::min(4, pairs)}; + for (std::size_t s = 0; s < starts; ++s) + { + auto const i{starts > 1 ? s * (pairs - 1) / (starts - 1) : 0}; + auto along{Points[i + 1] - Points[i]}; + if (glm::length(along) < 1e-3) + continue; + if (glm::dot(along, travel) < 0.0) + along = -along; + auto const angle{std::atan2(along.y, along.x)}; + glm::dvec2 const across{-std::sin(angle), std::cos(angle)}; + std::vector parameters{angle, glm::dot(across, Points[i] - start), guess, guess}; + if (auto const cost{least_squares(parameters, residuals)}; cost < lowest) + { + lowest = cost; + best = parameters; + } + } + if (false == best.empty()) + { + auto longer{best}; + longer.push_back(std::sqrt(std::max(transition, 10.0))); + longer.push_back(longer.back()); + if (least_squares(longer, residuals) < lowest) + best = longer; + } + if (design two; false == best.empty() && shaped(best, two)) + { + for (auto &curve : two.vertices) + plain_curve(curve, std::max(kTraceRadiusMin, std::round(curve.radius)), std::round(curve.transition_in), speed, two.norms); + take(stretch_of(two, Points)); + } + } + return chosen; +} + +stretch_fitter::~stretch_fitter() +{ + { + std::lock_guard lock{m_mutex}; + m_stop = true; + } + m_wake.notify_all(); + if (m_thread.joinable()) + m_thread.join(); +} + +void stretch_fitter::ask(request Request) +{ + { + std::lock_guard lock{m_mutex}; + if ((m_asked.has_value() && m_asked->key == Request.key) || (m_fitting.has_value() && m_fitting->key == Request.key) || (m_result.has_value() && m_result->first == Request.key)) + return; + m_asked = std::move(Request); + if (false == m_thread.joinable()) + m_thread = std::thread{[this]() { work(); }}; + } + m_wake.notify_one(); +} + +bool stretch_fitter::done(std::vector const &Key, trace_stretch &Fit) const +{ + std::lock_guard lock{m_mutex}; + if (false == m_result.has_value() || m_result->first != Key) + return false; + Fit = m_result->second; + return true; +} + +bool stretch_fitter::latest(std::vector &Key, trace_stretch &Fit) const +{ + std::lock_guard lock{m_mutex}; + if (false == m_result.has_value()) + return false; + Key = m_result->first; + Fit = m_result->second; + return true; +} + +void stretch_fitter::work() +{ + std::unique_lock lock{m_mutex}; + for (;;) + { + m_wake.wait(lock, [this]() { return m_stop || m_asked.has_value(); }); + if (m_stop) + return; + m_fitting = std::move(m_asked); + m_asked.reset(); + auto const asked{*m_fitting}; + lock.unlock(); + auto fit{fit_trace_stretch(asked.base, asked.points, asked.radii, asked.straights)}; + lock.lock(); + m_result.emplace(asked.key, std::move(fit)); + m_fitting.reset(); + } +} + } diff --git a/editor/editorAlignment.hpp b/editor/editorAlignment.hpp index c91801af..7b411533 100644 --- a/editor/editorAlignment.hpp +++ b/editor/editorAlignment.hpp @@ -9,8 +9,12 @@ http: #pragma once +#include +#include #include #include +#include +#include #include #include @@ -62,6 +66,9 @@ struct vertex bool compound{false}; double share{1.0}; std::vector arcs; + // two transition curves meeting at the curvature of the radius, with no arc between them and no cant: the bend of two straights + // eased the way the lines are drawn. Their lengths follow from the radius and the deflection, the ones given are left out + bool spirals{false}; }; struct design @@ -144,6 +151,8 @@ struct curve_report double jerk_in{0.0}; double jerk_out{0.0}; recommendation recommended; + bool spirals{false}; // two transition curves with no arc between them + double cant{0.0}; }; struct profile_point @@ -317,4 +326,66 @@ std::optional lay_on(lay_request const &Request); // the transitions between the arcs of a compound curve which the norms ask for at the speed, the cant of the vertex on its sharpest arc void recommend_between(vertex &Vertex, double const Speed, limits const &Limits); +// the length of each of the two transition curves of a vertex with spirals, meeting at the curvature of the radius +double spirals_length(transition_shape const Shape, double const Deflection, double const Radius); +// the curve keeps to the norms at the speed of the design: the radius, the acceleration across, the transitions long enough +bool meets_norms(curve_report const &Curve, design const &Design); + +// the radii the bends of the straights are eased with: the two transition curves meet at a radius from High down to Low, or above +// High where they'd be too short for the speed; where they don't fit, an arc with the transitions and the cant of the norms +struct corner_radii +{ + double low{3000.0}; + double high{6000.0}; +}; +// the vertices marked are given the curves of the bends of straights clicked on the photo, as large as fit between the neighbouring +// ones, these sharing the straights between them; the vertices stay where they are, a curve below the norms is left with its warnings +result ease_corners(design &Design, std::vector const &Corners, corner_radii const &Radii); + +// a stretch of the line traced on the photo: the points clicked along the middle of the track, from the straight Base starts on to +// the one it ends on, both fixed by the points clicked on them +struct trace_stretch +{ + bool valid{false}; + // the curves of the stretch with their places: one, of up to three radii, or two with a straight between them where the points + // follow those better; none where the stretch is straight + std::vector vertices; + std::vector offsets; // of the points from the line laid, positive to the left + double rms{0.0}; + double worst{0.0}; +}; +trace_stretch fit_trace_stretch(design const &Base, std::vector const &Points, int const Radii, bool const Straights); + +// the stretches traced fitted on a thread of their own while the cursor moves: fitting one takes tens of milliseconds. The latest one +// asked for is fitted, the ones asked for meanwhile are left out +class stretch_fitter +{ + public: + struct request + { + std::vector key; // tells the requests apart + design base; + std::vector points; + int radii{0}; + bool straights{true}; + }; + ~stretch_fitter(); + void ask(request Request); + // the fit of the request of the key, once it's done + bool done(std::vector const &Key, trace_stretch &Fit) const; + // the fit done last, of whatever request + bool latest(std::vector &Key, trace_stretch &Fit) const; + + private: + void work(); + + mutable std::mutex m_mutex; + std::condition_variable m_wake; + std::optional m_asked; + std::optional m_fitting; + std::optional, trace_stretch>> m_result; + bool m_stop{false}; + std::thread m_thread; +}; + } diff --git a/editor/editorAlignmentStore.cpp b/editor/editorAlignmentStore.cpp index 37b70795..0b9a6d8d 100644 --- a/editor/editorAlignmentStore.cpp +++ b/editor/editorAlignmentStore.cpp @@ -45,7 +45,8 @@ bool g_loaded{false}; // ends // norms // pieces -// vertex [ ]... +// vertex [ ]... [spirals 1] +// spirals: two transition curves with no arc between them, their lengths following from R // element // radii signed with the turn, 0 for no curvature // end @@ -136,6 +137,8 @@ void mark(scene::layer_handle const Layer) vertex.share, vertex.arcs.size())}; for (auto const &arc : vertex.arcs) text += format(" %.3f %.3f %.9f", arc.radius, arc.transition, arc.share); + if (vertex.spirals) + text += " spirals 1"; lines.push_back(text); } for (auto const &e : entry.result.elements) @@ -307,6 +310,10 @@ std::optional read_vertex(std::istringstream &Line, scene::la } if (Line.fail()) return std::nullopt; + // written since the bends of the straights are eased with two transition curves, the lines of before end with the arcs + for (std::string word; Line >> word;) + if (int value{0}; word == "spirals" && Line >> value) + vertex.spirals = (value != 0); position = Context.to_world(position); vertex.position = {position.x, position.z}; vertex.reverse_turn = (reverse != 0); diff --git a/editor/editorGroundFollow.cpp b/editor/editorGroundFollow.cpp index 926a82e5..19220bb4 100644 --- a/editor/editorGroundFollow.cpp +++ b/editor/editorGroundFollow.cpp @@ -371,4 +371,86 @@ void follow_meanwhile(std::vector &Pieces, ends const &Ends, bed_q lay_on(Pieces, lengths, total, Ends, stretched ? &graded : nullptr, stretched ? last.total / total : 1.0); } +void follow_points(std::vector &Pieces, std::vector const &Points, ends const &Ends) +{ + std::vector lengths; + double total{0.0}; + for (auto const &piece : Pieces) + total += lengths.emplace_back(bezier{piece}.plan_length(256)); + if (total < 1e-3) + return; + // the plan of the pieces as a line of short chords, each point of it with its chainage + std::vector line; + std::vector chainages; + { + double from{0.0}; + for (std::size_t i = 0; i < Pieces.size(); ++i) + { + bezier const curve{Pieces[i]}; + auto const steps{std::max(4, static_cast(std::ceil(lengths[i] / kStep)))}; + for (int j = (i == 0 ? 0 : 1); j <= steps; ++j) + { + line.push_back(plan_of(curve.point(curve.parameter(lengths[i] * j / steps)))); + chainages.push_back(from + lengths[i] * j / steps); + } + from += lengths[i]; + } + } + auto const chainage_of = [&](glm::dvec2 const &Point) { + auto best{std::numeric_limits::max()}; + double chainage{0.0}; + for (std::size_t i = 0; i + 1 < line.size(); ++i) + { + auto const chord{line[i + 1] - line[i]}; + auto const length2{glm::dot(chord, chord)}; + auto const t{length2 > 1e-18 ? std::clamp(glm::dot(Point - line[i], chord) / length2, 0.0, 1.0) : 0.0}; + auto const offset{Point - (line[i] + chord * t)}; + if (auto const distance2{glm::dot(offset, offset)}; distance2 < best) + { + best = distance2; + chainage = chainages[i] + (chainages[i + 1] - chainages[i]) * t; + } + } + return chainage; + }; + // the knots: the ends, and between them the points in their order along the track, none too near the one before + std::vector knots{0.0}; + std::vector heights{Ends.height0}; + for (auto const &point : Points) + { + auto const chainage{chainage_of(plan_of(point))}; + if (chainage < knots.back() + 1.0 || chainage > total - 1.0) + continue; + knots.push_back(chainage); + heights.push_back(point.y); + } + knots.push_back(total); + heights.push_back(Ends.height1); + auto const last{knots.size() - 1}; + // the height at the chainage on the stretch Middle is on, an end of it taken at a knot + auto const at = [&](double const Chainage, double const Middle, double &Grade) { + std::size_t stretch{0}; + while (stretch + 1 < last && Middle > knots[stretch + 1]) + ++stretch; + auto const span{knots[stretch + 1] - knots[stretch]}; + auto const even{(heights[stretch + 1] - heights[stretch]) / span}; + auto const from{stretch == 0 && false == std::isnan(Ends.grade0) ? Ends.grade0 : even}; + auto const to{stretch + 1 == last && false == std::isnan(Ends.grade1) ? Ends.grade1 : even}; + // a cubic of Hermite, a straight grade where both ends keep the even one + auto const t{std::clamp((Chainage - knots[stretch]) / span, 0.0, 1.0)}; + auto const t2{t * t}; + auto const t3{t2 * t}; + Grade = ((6.0 * t2 - 6.0 * t) * heights[stretch] + (3.0 * t2 - 4.0 * t + 1.0) * from * span + (6.0 * t - 6.0 * t2) * heights[stretch + 1] + (3.0 * t2 - 2.0 * t) * to * span) / span; + return (2.0 * t3 - 3.0 * t2 + 1.0) * heights[stretch] + (t3 - 2.0 * t2 + t) * from * span + (3.0 * t2 - 2.0 * t3) * heights[stretch + 1] + (t3 - t2) * to * span; + }; + std::vector> spans; + auto result{cut(Pieces, lengths, std::vector(knots.begin() + 1, knots.end() - 1), spans)}; + for (std::size_t i = 0; i < result.size(); ++i) + { + auto const middle{(spans[i].first + spans[i].second) * 0.5}; + lift(result[i], spans[i], [&](double const Chainage, double &Grade) { return at(Chainage, middle, Grade); }); + } + Pieces = std::move(result); +} + } // namespace ground_follow diff --git a/editor/editorGroundFollow.hpp b/editor/editorGroundFollow.hpp index aa45216f..d8e29943 100644 --- a/editor/editorGroundFollow.hpp +++ b/editor/editorGroundFollow.hpp @@ -102,5 +102,8 @@ void follow(std::vector &Pieces, ends const &Ends, bed_query const // the same, though the line is fitted by the fitter: until it has it, the track gets the line it fitted last stretched over its // length, if it went on from the same end, or the heights from end to end void follow_meanwhile(std::vector &Pieces, ends const &Ends, bed_query const &Bed, grading const &Grading, fitter &Fitter); +// from point to point: the points set square onto the pieces, an even grade between each two of them, the grade breaking only at +// them; from the ends of the heights given, a grade given there led on into the first and the last stretch +void follow_points(std::vector &Pieces, std::vector const &Points, ends const &Ends); } // namespace ground_follow diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 7bb62410..38bf6cb2 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -14,6 +14,9 @@ add_executable(${UNIT_TESTS_NAME} main.cpp vehicle.cpp engine_stubs.cpp + editor_stubs.cpp + # the geometry of the lines the track editor designs, standing on its own + "${CMAKE_SOURCE_DIR}/editor/editorAlignment.cpp" ${UNIT_TEST_CASES}) # the core library brings the game's build settings along diff --git a/tests/unit/editor_stubs.cpp b/tests/unit/editor_stubs.cpp new file mode 100644 index 00000000..54c1f4a9 --- /dev/null +++ b/tests/unit/editor_stubs.cpp @@ -0,0 +1,13 @@ +// stand-ins for the parts of the engine which the code of the editor under test reaches for: the translations give the English +// messages back as they are + +#include "stdafx.h" +#include "utilities/translation.h" + +locale Translations; + +const char *locale::lookup_c(const char *msg, bool) { return msg; } + +const std::string &locale::lookup_s(const std::string &msg, bool) { return msg; } + +const char *locale::lookup_n(const char *singular, const char *plural, long long n) { return n == 1 ? singular : plural; } diff --git a/tests/unit/test_alignment.cpp b/tests/unit/test_alignment.cpp new file mode 100644 index 00000000..8bbb8437 --- /dev/null +++ b/tests/unit/test_alignment.cpp @@ -0,0 +1,247 @@ +#include "stdafx.h" +#include "doctest.h" +#include "editor/editorAlignment.hpp" + +#include + +namespace { + +double constexpr pi{3.14159265358979323846}; + +glm::dvec2 heading(double const Angle) { return {std::cos(Angle), std::sin(Angle)}; } + +// straights from the origin along x, bending at the vertices by the deflections given (radians, to the left positive), each leg the +// length given; the vertices clicked as bends of straights +alignment::design straights(std::vector const &Legs, std::vector const &Deflections, double const Speed) { + alignment::design design; + design.speed = Speed; + glm::dvec2 point{0.0}; + double direction{0.0}; + std::vector corners; + for (std::size_t i = 0; i < Legs.size(); ++i) { + point += heading(direction) * Legs[i]; + if (i < Deflections.size()) { + corners.push_back(point); + direction += Deflections[i]; + } + } + design.start = {0.0, 0.0, 0.0}; + design.end = {point.x, 0.0, point.y}; + design.start_direction = heading(0.0); + design.end_direction = heading(direction); + for (auto const &corner : corners) { + alignment::vertex vertex; + vertex.position = corner; + design.vertices.push_back(vertex); + } + if (design.vertices.size() >= 2) { + design.vertices.front().offset = Legs.front(); + design.vertices.back().offset = Legs.back(); + } + return design; +} + +// points along the curves of the line every Spacing metres, from 30 m before the first curve to its last one, and two on the straight +// after it, shifted across by up to Noise +std::vector traced(alignment::result const &Line, double const Spacing, double const Noise) { + auto from{Line.length}; + auto to{0.0}; + for (auto const &e : Line.elements) + if (e.kind != alignment::element_kind::straight) { + from = std::min(from, e.chainage); + to = std::max(to, e.chainage + e.length); + } + std::vector points; + auto const along = [&](double const Chainage, double const Shift) { + auto const station{alignment::evaluate(Line, Chainage)}; + glm::dvec2 const across{-station.direction.z, station.direction.x}; + return glm::dvec2{station.position.x, station.position.z} + glm::normalize(across) * Shift; + }; + int count{0}; + for (auto chainage{std::max(0.0, from - 30.0)}; chainage <= to; chainage += Spacing, ++count) + points.push_back(along(chainage, Noise * (count % 3 == 0 ? 1.0 : count % 3 == 1 ? -0.6 : 0.2))); + points.push_back(along(to + 30.0, Noise * 0.5)); + points.push_back(along(to + 60.0, -Noise * 0.5)); + return points; +} + +// the stretch traced as the lay tool asks for it: from the start of the line along its first straight, to the last point along the +// straight through the last two +alignment::design stretch_base(alignment::design const &Line, std::vector const &Points) { + alignment::design base{Line}; + base.vertices.clear(); + base.end = {Points.back().x, 0.0, Points.back().y}; + base.end_direction = glm::normalize(Points.back() - Points[Points.size() - 2]); + return base; +} + +} // namespace + +TEST_CASE("alignment: a slight bend of straights is eased with two transition curves meeting, with no cant") { + auto design{straights({500.0, 500.0}, {3.0 * pi / 180.0}, 160.0)}; + auto const line{alignment::ease_corners(design, {1}, {3000.0, 6000.0})}; + + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 1); + auto const &curve{line.curves.front()}; + CHECK(curve.spirals); + CHECK(curve.cant == doctest::Approx(0.0)); + CHECK(curve.arc_length < 0.01); + CHECK(curve.radius >= 3000.0 - 1e-6); + CHECK(curve.radius <= 6000.0 + 1e-6); + CHECK(alignment::meets_norms(curve, design)); + // the bend stays where it was clicked + CHECK(line.vertices.front().x == doctest::Approx(500.0)); + CHECK(line.vertices.front().y == doctest::Approx(0.0)); +} + +TEST_CASE("alignment: a bend too slight for transitions long enough at the radii given takes a larger radius") { + auto design{straights({500.0, 500.0}, {0.2 * pi / 180.0}, 160.0)}; + auto const line{alignment::ease_corners(design, {1}, {3000.0, 6000.0})}; + + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 1); + auto const &curve{line.curves.front()}; + CHECK(curve.spirals); + CHECK(curve.radius > 6000.0); + CHECK(curve.transition_in >= alignment::transition_length(160.0, curve.radius, 0.0, design.norms) - 1e-6); + CHECK(alignment::meets_norms(curve, design)); +} + +TEST_CASE("alignment: a sharp bend between short straights is laid where clicked, an arc with cant below the norms") { + auto design{straights({120.0, 120.0}, {60.0 * pi / 180.0}, 160.0)}; + auto const line{alignment::ease_corners(design, {1}, {3000.0, 6000.0})}; + + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 1); + auto const &curve{line.curves.front()}; + CHECK_FALSE(curve.spirals); + CHECK(curve.cant > 0.0); + CHECK(curve.radius < alignment::recommend(160.0, curve.radius, design.norms).radius_min); + CHECK_FALSE(alignment::meets_norms(curve, design)); + CHECK_FALSE(line.warnings.empty()); + CHECK(curve.tangent_in <= 120.0 + 1e-6); + CHECK(line.vertices.front().x == doctest::Approx(120.0)); +} + +TEST_CASE("alignment: two bends share the straight between them") { + auto design{straights({600.0, 250.0, 600.0}, {4.0 * pi / 180.0, -5.0 * pi / 180.0}, 120.0)}; + auto const line{alignment::ease_corners(design, {1, 1}, {3000.0, 6000.0})}; + + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 2); + CHECK(line.curves[0].tangent_out + line.curves[1].tangent_in <= 250.0 + 1e-3); + CHECK(line.vertices[1].x == doctest::Approx(design.vertices[1].position.x)); + CHECK(line.vertices[1].y == doctest::Approx(design.vertices[1].position.y)); +} + +TEST_CASE("alignment: two transition curves meet with no arc in a design computed again") { + auto design{straights({500.0, 500.0}, {3.0 * pi / 180.0}, 120.0)}; + design.vertices.front().spirals = true; + design.vertices.front().radius = 4000.0; + design.vertices.front().cant = 60.0; + auto const line{alignment::compute(design)}; + + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 1); + auto const &curve{line.curves.front()}; + CHECK(curve.spirals); + CHECK(curve.cant == doctest::Approx(0.0)); + CHECK(curve.radius == doctest::Approx(4000.0)); + CHECK(curve.transition_in == doctest::Approx(alignment::spirals_length(design.shape, 3.0 * pi / 180.0, 4000.0))); + // the heading of the line out is that of the straight after the bend + auto const station{alignment::evaluate(line, line.length)}; + CHECK(std::atan2(station.direction.z, station.direction.x) == doctest::Approx(3.0 * pi / 180.0).epsilon(1e-6)); +} + +TEST_CASE("alignment: the radius and the transitions of a curve traced along a known one") { + auto truth{straights({500.0, 500.0}, {25.0 * pi / 180.0}, 120.0)}; + truth.vertices.front().radius = 1000.0; + truth.vertices.front().transition_in = truth.vertices.front().transition_out = 120.0; + auto const line{alignment::compute(truth)}; + REQUIRE(line.valid); + auto const points{traced(line, 20.0, 0.03)}; + + auto const fit{alignment::fit_trace_stretch(stretch_base(truth, points), points, 0, true)}; + + REQUIRE(fit.valid); + REQUIRE(fit.vertices.size() == 1); + CHECK(fit.vertices.front().radius == doctest::Approx(1000.0).epsilon(0.02)); + CHECK(fit.vertices.front().transition_in == doctest::Approx(120.0).epsilon(0.1)); + CHECK(fit.rms < 0.04); +} + +TEST_CASE("alignment: a plain arc traced keeps short transitions, not the ones of the norms") { + auto truth{straights({500.0, 500.0}, {20.0 * pi / 180.0}, 120.0)}; + truth.vertices.front().radius = 800.0; + auto const line{alignment::compute(truth)}; + REQUIRE(line.valid); + auto const points{traced(line, 15.0, 0.0)}; + + auto const fit{alignment::fit_trace_stretch(stretch_base(truth, points), points, 0, true)}; + + REQUIRE(fit.valid); + REQUIRE(fit.vertices.size() == 1); + CHECK(fit.vertices.front().radius == doctest::Approx(800.0).epsilon(0.01)); + CHECK(fit.vertices.front().transition_in < 20.0); + CHECK(fit.rms < 0.02); +} + +TEST_CASE("alignment: a curve, a straight and a curve the other way traced as two curves") { + auto truth{straights({400.0, 700.0, 400.0}, {20.0 * pi / 180.0, -25.0 * pi / 180.0}, 120.0)}; + truth.vertices[0].radius = 1500.0; + truth.vertices[1].radius = 1000.0; + for (auto &vertex : truth.vertices) + vertex.transition_in = vertex.transition_out = 100.0; + auto const line{alignment::compute(truth)}; + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 2); + auto const points{traced(line, 20.0, 0.02)}; + + auto const fit{alignment::fit_trace_stretch(stretch_base(truth, points), points, 0, true)}; + + REQUIRE(fit.valid); + REQUIRE(fit.vertices.size() == 2); + CHECK(fit.vertices[0].radius == doctest::Approx(line.curves[0].radius).epsilon(0.05)); + CHECK(fit.vertices[1].radius == doctest::Approx(line.curves[1].radius).epsilon(0.05)); + CHECK(fit.rms < 0.05); +} + +TEST_CASE("alignment: two curves turning the same way with a straight between them traced") { + auto truth{straights({400.0, 600.0, 400.0}, {15.0 * pi / 180.0, 15.0 * pi / 180.0}, 120.0)}; + truth.vertices[0].radius = 1200.0; + truth.vertices[1].radius = 900.0; + for (auto &vertex : truth.vertices) + vertex.transition_in = vertex.transition_out = 90.0; + auto const line{alignment::compute(truth)}; + REQUIRE(line.valid); + REQUIRE(line.curves.size() == 2); + auto const points{traced(line, 20.0, 0.02)}; + + auto const fit{alignment::fit_trace_stretch(stretch_base(truth, points), points, 0, true)}; + + REQUIRE(fit.valid); + REQUIRE(fit.vertices.size() == 2); + CHECK(fit.vertices[0].radius == doctest::Approx(1200.0).epsilon(0.05)); + CHECK(fit.vertices[1].radius == doctest::Approx(900.0).epsilon(0.05)); + CHECK(fit.rms < 0.04); +} + +TEST_CASE("alignment: a compound curve traced is followed closely") { + auto truth{straights({500.0, 500.0}, {35.0 * pi / 180.0}, 100.0)}; + auto &vertex{truth.vertices.front()}; + vertex.radius = 1200.0; + vertex.compound = true; + vertex.share = 0.5; + vertex.arcs = {{700.0, 60.0, 0.5}}; + vertex.transition_in = vertex.transition_out = 80.0; + auto const line{alignment::compute(truth)}; + REQUIRE(line.valid); + auto const points{traced(line, 15.0, 0.02)}; + + auto const fit{alignment::fit_trace_stretch(stretch_base(truth, points), points, 0, true)}; + + REQUIRE(fit.valid); + CHECK(fit.rms < 0.05); + CHECK(fit.worst < 0.15); +}