diff --git a/application/editormode.h b/application/editormode.h index 4453aa5e..72f02054 100644 --- a/application/editormode.h +++ b/application/editormode.h @@ -1464,8 +1464,14 @@ class editor_mode : public application_mode, private editor_track::observer, pri lay_line_plan *Plan = nullptr) const; bool lay_line(std::vector const &Points, std::vector const &Modes, editor_track::snap_target const &End, lay_line_plan &Plan, std::string &Error) const; // the curves of a stretch traced, the ones fitted lately kept for the same points - // none while the stretch the cursor moves along is still being fitted - std::optional traced_stretch(alignment::design const &Base, std::vector const &Points) const; + // none while the stretch the cursor moves along is still being fitted; Meanwhile: the fit done last of the same points but the one + // at the cursor until then + std::optional traced_stretch(alignment::design const &Base, std::vector const &Points, int const Radii, bool const Straights, bool const Meanwhile = true) const; + // the curve through a row of bends clicked close along one, by traced_stretch + alignment::stretch_fit bends_fit() const + { + return [this](alignment::design const &Base, std::vector const &Points) { return traced_stretch(Base, Points, 0, false, false); }; + } mutable std::vector, alignment::trace_stretch>> m_stretch_cache; // the stretch the cursor moves along is fitted on a thread, the frames drawn meanwhile with the fit done last mutable alignment::stretch_fitter m_stretch_fitter; diff --git a/application/editormodelay.cpp b/application/editormodelay.cpp index 1c550abe..0161a1a5 100644 --- a/application/editormodelay.cpp +++ b/application/editormodelay.cpp @@ -1005,7 +1005,7 @@ std::vector editor_mode::lay_pieces(std::vector const easeddesign.vertices.front().offset = std::max(0.1, glm::dot(shape.vertices.front() - plan_of(easeddesign.start), glm::normalize(easeddesign.start_direction))); easeddesign.vertices.back().offset = std::max(0.1, glm::dot(plan_of(easeddesign.end) - shape.vertices.back(), glm::normalize(easeddesign.end_direction))); } - auto const easedshape{alignment::ease_corners(easeddesign, plan.corners, {EditorSettings.corner_low(), EditorSettings.corner_high()})}; + auto const easedshape{alignment::ease_corners(easeddesign, plan.corners, {EditorSettings.corner_low(), EditorSettings.corner_high()}, bends_fit())}; if (easedshape.valid) { Length = editor_track::eased_length(in, easedshape.length, out); @@ -1181,7 +1181,7 @@ bool editor_mode::lay_line(std::vector const &Points, std::vector const &Points, std::vector const &Points, std::vector editor_mode::traced_stretch(alignment::design const &Base, std::vector const &Points) const +std::optional editor_mode::traced_stretch(alignment::design const &Base, std::vector const &Points, int const Radii, bool const Straights, bool const Meanwhile) const { - auto const &tool{m_lay}; // the stretches traced before the one the cursor moves along are the same from frame to frame: their fits are kept std::vector key{Base.start.x, Base.start.z, Base.start_direction.x, Base.start_direction.y, Base.end.x, Base.end.z, Base.end_direction.x, Base.end_direction.y, Base.speed, - static_cast(Base.shape), static_cast(tool.radii), tool.straights_between ? 1.0 : 0.0, Base.start_reserve, Base.end_reserve}; + static_cast(Base.shape), static_cast(Radii), Straights ? 1.0 : 0.0, Base.start_reserve, Base.end_reserve}; for (auto const &point : Points) { key.push_back(point.x); @@ -1308,7 +1307,7 @@ std::optional editor_mode::traced_stretch(alignment::d alignment::trace_stretch fit; if (false == m_follow_preview) { - fit = alignment::fit_trace_stretch(Base, Points, tool.radii, tool.straights_between); + fit = alignment::fit_trace_stretch(Base, Points, Radii, Straights); keep(fit); return fit; } @@ -1319,10 +1318,10 @@ std::optional editor_mode::traced_stretch(alignment::d keep(fit); return fit; } - m_stretch_fitter.ask({key, Base, Points, tool.radii, tool.straights_between}); + m_stretch_fitter.ask({key, Base, Points, Radii, Straights}); std::vector last; auto const header{std::size_t{14}}; - if (m_stretch_fitter.latest(last, fit) && last.size() >= header + 2 && key.size() >= header + 2 && std::equal(key.begin(), key.begin() + 4, last.begin()) && + if (Meanwhile && m_stretch_fitter.latest(last, fit) && last.size() >= header + 2 && key.size() >= header + 2 && std::equal(key.begin(), key.begin() + 4, last.begin()) && std::equal(key.begin() + 8, key.begin() + header, last.begin() + 8)) { auto const common{std::min(last.size(), key.size()) - 2}; diff --git a/editor/editorAlignment.cpp b/editor/editorAlignment.cpp index bbbe3ab3..6bdeed29 100644 --- a/editor/editorAlignment.cpp +++ b/editor/editorAlignment.cpp @@ -62,7 +62,8 @@ local_point transition_point(transition_shape const Shape, double const Length, } auto const distance{Length * Parameter}; auto const angle = [&](double const S) { return S * S / (2.0 * Radius * Length); }; - int const steps{std::max(8, static_cast(std::ceil(distance / 0.5)) * 2)}; + // Simpson over 64 steps at most keeps a transition of a radian within a millimetre, however long + int const steps{std::clamp(static_cast(std::ceil(distance / 0.5)) * 2, 8, 64)}; double const h{distance / steps}; double x{0.0}; double y{0.0}; @@ -80,8 +81,16 @@ local_point spiral_point(double const Start, double const End, double const Leng { if (Distance <= 0.0) return {}; + // an arc, or a straight, has its point at once + if (Start == End) + { + auto const angle{Start * Distance}; + if (std::abs(angle) < 1e-9) + return {Distance, 0.0, 0.0}; + return {std::sin(angle) / Start, (1.0 - std::cos(angle)) / Start, angle}; + } auto const heading = [&](double const S) { return Start * S + (End - Start) * S * S / (2.0 * std::max(Length, 1e-9)); }; - int const steps{std::max(8, static_cast(std::ceil(Distance / 0.5)) * 2)}; + int const steps{std::clamp(static_cast(std::ceil(Distance / 0.5)) * 2, 8, 64)}; double const h{Distance / steps}; double x{0.0}; double y{0.0}; @@ -1973,9 +1982,16 @@ double least_squares(std::vector &Parameters, residuals_of const &Residu return cost; } +bool compound_residuals(design const &Trial, std::vector const &Points, std::vector &Residuals); + constexpr double kTraceRadiusMin{15.0}; constexpr double kTraceRadiusMax{50000.0}; constexpr double kTraceRobust{0.25}; +// m, how far the points clicked along a curve should be from the curve laid through them: a shape keeping to them no closer is +// given more radii +constexpr double kPointsWithin{0.3}; +// m, how far they may be from it, clicked by hand far out: a row of bends making a curve no closer is cut in parts +constexpr double kPointsLoose{1.0}; // the parameters: the logarithms of the radii, then the weights of the shares of all but the first radius, as logarithms too vertex traced_shape(design const &Design, vertex const &Base, std::vector const &Parameters, int const Radii) @@ -2035,6 +2051,9 @@ trace_fit fit_trace(design const &Design, std::size_t const Vertex, std::vector< auto trial{Design}; trial.norms.checked = false; trial.vertices[Vertex] = traced_shape(Design, base, Parameters, Count); + // a compound curve alone made straight away + if (Count >= 2 && trial.vertices.size() == 1) + return compound_residuals(trial, Points, Residuals); auto const result{compute(trial)}; if (false == result.valid) return false; @@ -2491,6 +2510,19 @@ std::vector polygon_of(design const &Design) return polygon; } +// the angle the line turns by at the vertex, nought where it can't be told +double vertex_deflection(design const &Design, std::size_t const Vertex) +{ + auto const polygon{polygon_of(Design)}; + if (Vertex + 2 >= polygon.size()) + return 0.0; + auto const in{polygon[Vertex + 1] - polygon[Vertex]}; + auto const out{polygon[Vertex + 2] - polygon[Vertex + 1]}; + if (glm::length(in) < 1e-9 || glm::length(out) < 1e-9) + return 0.0; + return std::atan2(std::abs(cross(in, out)), glm::dot(in, out)); +} + // 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) { @@ -2522,9 +2554,8 @@ struct corner_shape double reach{0.0}; }; -} // namespace - -result ease_corners(design &Design, std::vector const &Corners, corner_radii const &Radii) +// the bends eased one by one, each with the curve as large as fits between its neighbours +result ease_bends(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; }; @@ -2701,6 +2732,159 @@ result ease_corners(design &Design, std::vector const &Corners, corner_rad return compute(Design); } +// the curves the line lays at its vertices: the side each turns to, where each begins and ends +struct laid_curve +{ + int vertex{-1}; + int side{0}; + double from{0.0}; + double to{0.0}; +}; + +std::vector laid_curves(result const &Line) +{ + std::vector curves; + for (auto const &e : Line.elements) + { + if (e.vertex < 0 || e.kind == element_kind::straight) + continue; + if (curves.empty() || curves.back().vertex != e.vertex) + curves.push_back({e.vertex, e.turn, e.chainage, e.chainage + e.length}); + else + curves.back().to = e.chainage + e.length; + } + return curves; +} + +// the shortest straight the line keeps between two curves turning the same way, as many metres as km/h of the speed, 30 m at least: +// two curves with a shorter one between them are one compound curve +double shortest_between(design const &Design) +{ + return std::max(30.0, Design.speed); +} + +} // namespace + +result ease_corners(design &Design, std::vector &Corners, corner_radii const &Radii, stretch_fit const &Fit) +{ + // a point following a slight bend clicked closer than the curves leave a straight of the shortest length was clicked on a curve + // seen on the photo: the bends of a row of such, turning the same way, are the points of a curve traced through them, between + // the straights before and after them, one compound curve rather than short straights between curves, or several where one + // doesn't keep to them + auto line{ease_bends(Design, Corners, Radii)}; + auto const shortest{shortest_between(Design)}; + auto const bend = [&](int const K) { return K >= 0 && K < static_cast(Corners.size()) && Corners[K] == 1; }; + auto pending{false}; + for (int pass = 0; pass < 32 && line.valid && false == pending; ++pass) + { + auto const curves{laid_curves(line)}; + auto const polygon{polygon_of(Design)}; + if (polygon.size() != Design.vertices.size() + 2) + break; + // the rows of bends: each as a pair of the first and the last curve, the longest first + std::vector> rows; + for (std::size_t i = 0; i < curves.size();) + { + auto j{i}; + while (j + 1 < curves.size() && bend(curves[i].vertex) && bend(curves[j + 1].vertex) && curves[j + 1].vertex == curves[j].vertex + 1 && curves[j + 1].side == curves[i].side && + curves[j + 1].from - curves[j].to < shortest) + ++j; + if (j > i) + rows.emplace_back(i, j); + i = j + 1; + } + auto merged{false}; + for (auto row{rows.begin()}; row != rows.end() && false == merged && false == pending; ++row) + { + // a row which doesn't make one curve is tried in parts, cut at its longest straight + std::vector> parts{*row}; + while (false == parts.empty() && false == merged && false == pending) + { + auto const [first, last]{parts.back()}; + parts.pop_back(); + if (first == last) + continue; + auto const k{curves[first].vertex}; + auto const m{curves[last].vertex}; + design base{Design}; + base.vertices.clear(); + base.start = {polygon[k].x, Design.start.y, polygon[k].y}; + base.start_direction = glm::normalize(polygon[k + 1] - polygon[k]); + base.end = {polygon[m + 2].x, Design.end.y, polygon[m + 2].y}; + base.end_direction = glm::normalize(polygon[m + 2] - polygon[m + 1]); + base.start_reserve = 0.0; + base.end_reserve = 0.0; + std::vector points(polygon.begin() + k + 1, polygon.begin() + m + 2); + // two bends don't show the curve between them: the curve goes through the line their two curves laid, with no + // straight between them, from where the first begins to where the second ends, the straights either side its room + auto within{kPointsLoose}; + if (last == first + 1) + { + auto const from{curves[first].from}; + auto const to{curves[last].to}; + points.clear(); + auto const count{std::max(4, static_cast(std::ceil((to - from) / 5.0)))}; + for (int i = 0; i <= count; ++i) + { + auto const station{evaluate(line, from + (to - from) * i / count)}; + points.push_back(plan_of(station.position)); + } + auto const before{evaluate(line, from)}; + auto const after{evaluate(line, to)}; + base.start = before.position; + base.start_direction = glm::normalize(plan_of(before.direction)); + base.end = after.position; + base.end_direction = glm::normalize(plan_of(after.direction)); + base.start_reserve = from - (first > 0 ? curves[first - 1].to : 0.0); + base.end_reserve = (last + 1 < curves.size() ? curves[last + 1].from : line.length) - to; + within = kPointsWithin; + } + auto const fit{Fit ? Fit(base, points) : std::optional{fit_trace_stretch(base, points, 0, false)}}; + // still being fitted on the thread: the bends stay apart until it's done + if (false == fit.has_value()) + { + pending = true; + break; + } + if (fit->valid && fit->vertices.size() == 1 && fit->worst <= within) + { + // the row of bends one vertex, kept as it is + std::vector vertices(Design.vertices.begin(), Design.vertices.begin() + k); + vertices.push_back(fit->vertices.front()); + vertices.insert(vertices.end(), Design.vertices.begin() + m + 1, Design.vertices.end()); + std::vector corners(Corners.begin(), Corners.begin() + k); + corners.push_back(0); + corners.insert(corners.end(), Corners.begin() + m + 1, Corners.end()); + for (std::size_t i = 0; i < vertices.size(); ++i) + if (i < static_cast(k) || i > static_cast(k)) + vertices[i].position = polygon[(i < static_cast(k) ? i : i + (m - k)) + 1]; + if (vertices.size() >= 2) + { + vertices.front().offset = glm::dot(vertices.front().position - plan_of(Design.start), glm::normalize(Design.start_direction)); + vertices.back().offset = glm::dot(plan_of(Design.end) - vertices.back().position, glm::normalize(Design.end_direction)); + } + Design.vertices = std::move(vertices); + Corners = std::move(corners); + merged = true; + break; + } + if (last - first < 2) + continue; + auto cut{first}; + for (auto i{first}; i < last; ++i) + if (curves[i + 1].from - curves[i].to > curves[cut + 1].from - curves[cut].to) + cut = i; + parts.emplace_back(first, cut); + parts.emplace_back(cut + 1, last); + } + } + if (false == merged || pending) + break; + line = ease_bends(Design, Corners, Radii); + } + return line; +} + namespace { @@ -2748,6 +2932,70 @@ bool stretch_residuals(design Trial, std::vector const &Points, std: return true; } +// the same for a design of one compound curve, its shape made straight away as fit_between makes it, many times quicker than the whole +// design laid: the fits of compound curves try it some thousand times +bool compound_residuals(design const &Trial, std::vector const &Points, std::vector &Residuals) +{ + if (Trial.vertices.size() != 1 || false == Trial.vertices.front().compound || Trial.vertices.front().arcs.empty()) + return stretch_residuals(Trial, Points, Residuals); + glm::dvec2 corner; + if (false == tangent_intersection(Trial, corner)) + return false; + auto const start{plan_of(Trial.start)}; + auto const end{plan_of(Trial.end)}; + auto const in{glm::normalize(Trial.start_direction)}; + auto const out{glm::normalize(Trial.end_direction)}; + auto const turn{cross(in, out)}; + auto const deflection{std::atan2(std::abs(turn), glm::dot(in, out))}; + if (deflection < 1e-7 || deflection > kPi - 1e-3) + return false; + int const side{turn > 0.0 ? 1 : -1}; + auto const &vertex{Trial.vertices.front()}; + std::vector arcs{{std::max(1.0, vertex.radius), 0.0, vertex.share}}; + for (auto const &arc : vertex.arcs) + arcs.push_back({std::max(1.0, arc.radius), std::max(0.0, arc.transition), arc.share}); + auto const shape{build_compound(deflection, arcs, std::max(0.0, vertex.transition_in), std::max(0.0, vertex.transition_out))}; + auto const normalin{perpendicular(in) * static_cast(side)}; + auto const offset{in * shape.end.x + normalin * shape.end.y}; + auto const along{cross(out, offset) / cross(out, in)}; + auto const after{glm::dot(offset - in * along, out)}; + if (along <= 0.0 || after <= 0.0) + return false; + // a curve too large for the straights is made smaller as a whole, its radii and its transitions alike, as fit_between does it + auto const scale{std::min({1.0, (glm::distance(corner, start) + std::max(0.0, Trial.start_reserve)) / along, (glm::distance(end, corner) + std::max(0.0, Trial.end_reserve)) / after})}; + // the curve in chords of 4 m at most, the straights either side carried on far + constexpr double kFar{5000.0}; + auto const begin{corner - in * along * scale}; + auto origin{glm::dvec2{0.0}}; + auto heading{in}; + std::vector line{begin - in * kFar, begin}; + for (auto const &segment : shape.segments) + { + auto const normal{perpendicular(heading) * static_cast(side)}; + auto const steps{std::max(1, static_cast(std::ceil(segment[2] * scale / 4.0)))}; + local_point local; + for (int i = 1; i <= steps; ++i) + { + local = spiral_point(segment[0], segment[1], segment[2], segment[2] * i / steps); + line.push_back(begin + (origin + heading * local.x + normal * local.y) * scale); + } + origin += heading * local.x + normal * local.y; + heading = turned(heading, side * local.angle); + } + line.push_back(line.back() + out * kFar); + plan_line plan; + plan.points = std::move(line); + plan.chainages.assign(plan.points.size(), 0.0); + Residuals.clear(); + for (auto const &point : Points) + { + auto const offset{project_on(plan, 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) { @@ -2766,10 +3014,12 @@ void plain_curve(vertex &Vertex, double const Radius, double const Transition, d void fit_curve(design &Design, std::size_t const Vertex, std::vector const &Points) { auto const speed{std::max(1.0, Design.speed)}; + // the transitions together turn no more than the curve: longer ones would only be cut to it + auto const deflection{vertex_deflection(Design, Vertex)}; 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); + plain_curve(Trial.vertices[Vertex], radius, std::min(Parameters[1] * Parameters[1], radius * deflection), speed, Trial.norms); }; residuals_of const residuals{[&](std::vector const &Parameters, std::vector &Residuals) { design trial; @@ -2801,14 +3051,24 @@ void fit_curve(design &Design, std::size_t const Vertex, std::vector } // 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) +// best, from the shape it has; Between: the transitions of the norms between the arcs, else the arcs meet with none +void fit_compound(design &Design, std::size_t const Vertex, std::vector const &Points, bool const Between = true) { 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}; + // the transitions into the first arc and out of the last one turn no more than those arcs + auto const deflection{vertex_deflection(Design, Vertex)}; + auto const reach = [&](vertex const &Curve) { + auto total{std::max(0.0, Curve.share)}; + for (auto const &arc : Curve.arcs) + total += std::max(0.0, arc.share); + auto const first{total > 0.0 ? std::max(0.0, Curve.share) / total : 1.0}; + auto const last{total > 0.0 && false == Curve.arcs.empty() ? std::max(0.0, Curve.arcs.back().share) / total : 1.0}; + return std::min(Curve.radius * first, Curve.arcs.empty() ? Curve.radius : Curve.arcs.back().radius * last) * deflection; + }; // 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) @@ -2818,13 +3078,18 @@ void fit_compound(design &Design, std::size_t const Vertex, 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); + auto shape{Design}; + shape.norms.checked = Between; + curve = traced_shape(shape, base, Parameters, radii); + if (false == Between) + for (auto &arc : curve.arcs) + arc.transition = 0.0; + curve.transition_in = curve.transition_out = std::min(Parameters.back() * Parameters.back(), reach(curve)); }; residuals_of const residuals{[&](std::vector const &Parameters, std::vector &Residuals) { design trial; shaped(Parameters, trial); - return stretch_residuals(trial, Points, Residuals); + return compound_residuals(trial, Points, Residuals); }}; std::vector best; auto lowest{std::numeric_limits::max()}; @@ -2850,10 +3115,89 @@ void fit_compound(design &Design, std::size_t const Vertex, std::vector 0.0 ? ramped_cant(recommend(allowed, sharpest, fitted.norms).cant, curve.transition_in, allowed, fitted.norms) : 0.0; - recommend_between(curve, speed, fitted.norms); + if (Between) + recommend_between(curve, speed, fitted.norms); Design = fitted; } +// the one curve of the design, where three radii don't keep within kPointsWithin of the points: a compound curve of more of them, +// up to Most, as few as keep to the points; each arc started from the curvature the points show along its share of the curve +std::optional fit_arcs(design const &Single, std::vector const &Points, int const Most) +{ + auto const laid{compute(Single)}; + if (false == laid.valid || Points.size() < 3) + return std::nullopt; + auto from{laid.length}; + auto to{0.0}; + for (auto const &e : laid.elements) + if (e.vertex == 0) + { + from = std::min(from, e.chainage); + to = std::max(to, e.chainage + e.length); + } + if (from >= to) + return std::nullopt; + // the curvature at each point of the circle through it and the points either side, at the chainage of the point on the curve + auto const projections{project(laid, Points)}; + std::vector> curvatures; + for (std::size_t i = 1; i + 1 < Points.size(); ++i) + { + auto const a{glm::distance(Points[i - 1], Points[i])}; + auto const b{glm::distance(Points[i], Points[i + 1])}; + auto const c{glm::distance(Points[i - 1], Points[i + 1])}; + if (a * b * c < 1e-9) + continue; + curvatures.emplace_back(projections[i].chainage, 2.0 * std::abs(cross(Points[i] - Points[i - 1], Points[i + 1] - Points[i - 1])) / (a * b * c)); + } + auto const &single{Single.vertices.front()}; + std::optional best; + for (int arcs = 4; arcs <= Most; ++arcs) + { + auto shaped{Single}; + auto &curve{shaped.vertices.front()}; + curve.compound = true; + curve.spirals = false; + curve.arcs.clear(); + auto const span{(to - from) / arcs}; + for (int j = 0; j < arcs; ++j) + { + double sum{0.0}; + int count{0}; + for (auto const &[chainage, curvature] : curvatures) + if (chainage >= from + span * j && chainage < from + span * (j + 1)) + { + sum += curvature; + ++count; + } + auto const curvature{count > 0 ? sum / count : 1.0 / std::max(kTraceRadiusMin, single.radius)}; + auto const radius{std::clamp(1.0 / std::max(curvature, 1.0 / kTraceRadiusMax), kTraceRadiusMin, kTraceRadiusMax)}; + // the share of the angle: of the curvature along the span + auto const share{span * std::max(curvature, 1e-6)}; + if (j == 0) + { + curve.radius = radius; + curve.share = share; + } + else + curve.arcs.push_back({radius, 0.0, share}); + } + // arcs meeting with no transitions between them: the points of a curve clicked by hand ask for the radius to change + // quicker than the transitions of the norms let it + fit_compound(shaped, 0, Points, false); + // none of an arc too short to be one: a kink the points made + auto const fitted{compute(shaped)}; + auto const kinked{std::any_of(fitted.elements.begin(), fitted.elements.end(), [](element const &Element) { + return Element.vertex == 0 && Element.curvature_start == Element.curvature_end && Element.curvature_start > 0.0 && Element.length < 10.0; + })}; + auto const candidate{kinked ? std::nullopt : stretch_of(shaped, Points)}; + if (candidate.has_value() && (false == best.has_value() || candidate->worst < best->worst)) + best = candidate; + if (best.has_value() && best->worst <= kPointsWithin) + break; + } + return best; +} + } // namespace trace_stretch fit_trace_stretch(design const &Base, std::vector const &Points, int const Radii, bool const Straights) @@ -2898,16 +3242,44 @@ trace_stretch fit_trace_stretch(design const &Base, std::vector cons { curve = arc.fitted; fit_curve(single, 0, Points); - one = stretch_of(single, Points); - if (one.has_value() && Radii != 1 && (Radii > 1 || one->rms > kStretchClose)) + // a fit running off to the largest radius, its transitions eating the arc, is two transitions meeting + auto const &fitted{single.vertices.front()}; + auto const runaway{fitted.radius > 0.95 * kTraceRadiusMax}; + if (false == runaway) + one = stretch_of(single, Points); + auto spirals{single}; + spirals.vertices.front() = {}; + spirals.vertices.front().spirals = true; + residuals_of const peak{[&](std::vector const &Parameters, std::vector &Residuals) { + auto trial{spirals}; + trial.vertices.front().radius = std::exp(std::clamp(Parameters[0], std::log(kTraceRadiusMin), std::log(kTraceRadiusMax))); + return stretch_residuals(trial, Points, Residuals); + }}; + std::vector radius{std::log(std::clamp(arc.fitted.radius, kTraceRadiusMin, kTraceRadiusMax))}; + least_squares(radius, peak); + spirals.vertices.front().radius = std::round(std::exp(std::clamp(radius[0], std::log(kTraceRadiusMin), std::log(kTraceRadiusMax)))); + if (auto const candidate{stretch_of(spirals, Points)}; candidate.has_value() && (false == one.has_value() || candidate->rms < one->rms)) + one = candidate; + // a compound curve is the one curve which keeps to the points the closest: of up to three radii, of more where those leave + // a point further than kPointsWithin off + auto const loose = [&]() { return one->rms > kStretchClose || one->worst > kPointsWithin; }; + if (one.has_value() && Radii != 1 && (Radii > 1 || loose())) 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)) + if (auto const candidate{stretch_of(several, Points)}; + candidate.has_value() && (Radii > 1 || candidate->rms < 0.75 * one->rms || (one->worst > kPointsWithin && candidate->worst < one->worst))) one = candidate; } + if (one.has_value() && Radii == 0 && one->worst > kPointsWithin) + { + auto best{single}; + best.vertices.front() = one->vertices.front(); + if (auto const candidate{fit_arcs(best, Points, 8)}; candidate.has_value() && candidate->worst < one->worst) + one = candidate; + } } take(one); } @@ -2943,7 +3315,7 @@ trace_stretch fit_trace_stretch(design const &Base, std::vector cons { 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); + plain_curve(Trial.vertices[i], radius, std::min(length, radius * vertex_deflection(Trial, i)), speed, Trial.norms); } return true; }; @@ -2993,8 +3365,13 @@ trace_stretch fit_trace_stretch(design const &Base, std::vector cons auto const sound{laid.valid && laid.radius_warnings.empty() && laid.curves.size() == 2 && std::all_of(laid.curves.begin(), laid.curves.end(), [](curve_report const &Curve) { return Curve.radius > 1.05 * kTraceRadiusMin && Curve.radius < 0.95 * kTraceRadiusMax && Curve.deflection > 0.25 * kPi / 180.0; })}; - if (sound) - take(stretch_of(two, Points)); + // one compound curve rather than two turning the same way with a short straight between them; two such only where the + // straight is long and they keep to the points much closer + auto const candidate{sound ? stretch_of(two, Points) : std::nullopt}; + auto const curves{laid_curves(laid)}; + auto const same{curves.size() == 2 && curves[0].side == curves[1].side}; + if (candidate.has_value() && (false == same || (curves[1].from - curves[0].to >= shortest_between(Base) && (false == chosen.valid || candidate->rms < 0.5 * chosen.rms)))) + take(candidate); } } return chosen; @@ -3015,9 +3392,15 @@ 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)) + auto const same = [&](std::vector const &Key) { return Key == Request.key; }; + if (std::any_of(m_asked.begin(), m_asked.end(), [&](request const &Asked) { return same(Asked.key); }) || (m_fitting.has_value() && same(m_fitting->key)) || + std::any_of(m_results.begin(), m_results.end(), [&](auto const &Result) { return same(Result.first); })) return; - m_asked = std::move(Request); + m_asked.push_back(std::move(Request)); + // the cursor moves on: what was asked for long ago waits no more + constexpr std::size_t kWaiting{8}; + if (m_asked.size() > kWaiting) + m_asked.erase(m_asked.begin()); if (false == m_thread.joinable()) m_thread = std::thread{[this]() { work(); }}; } @@ -3027,19 +3410,22 @@ void stretch_fitter::ask(request Request) 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; + for (auto const &result : m_results) + if (result.first == Key) + { + Fit = result.second; + return true; + } + return false; } bool stretch_fitter::latest(std::vector &Key, trace_stretch &Fit) const { std::lock_guard lock{m_mutex}; - if (false == m_result.has_value()) + if (m_results.empty()) return false; - Key = m_result->first; - Fit = m_result->second; + Key = m_results.front().first; + Fit = m_results.front().second; return true; } @@ -3048,16 +3434,20 @@ void stretch_fitter::work() std::unique_lock lock{m_mutex}; for (;;) { - m_wake.wait(lock, [this]() { return m_stop || m_asked.has_value(); }); + m_wake.wait(lock, [this]() { return m_stop || false == m_asked.empty(); }); if (m_stop) return; - m_fitting = std::move(m_asked); - m_asked.reset(); + // the one asked for last first: it's what the cursor shows now + m_fitting = std::move(m_asked.back()); + m_asked.pop_back(); 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)); + constexpr std::size_t kKept{16}; + m_results.insert(m_results.begin(), {asked.key, std::move(fit)}); + if (m_results.size() > kKept) + m_results.pop_back(); m_fitting.reset(); } } diff --git a/editor/editorAlignment.hpp b/editor/editorAlignment.hpp index 7b411533..0dc9ef6d 100644 --- a/editor/editorAlignment.hpp +++ b/editor/editorAlignment.hpp @@ -10,6 +10,7 @@ http: #pragma once #include +#include #include #include #include @@ -338,26 +339,32 @@ 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 + // follow those much 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 fit of a stretch traced, with up to three radii and no straight in it; none while it's being fitted elsewhere +using stretch_fit = std::function(design const &, std::vector const &)>; -// 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 +// the vertices marked 1 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. A row of bends turning the same way, clicked so close that the straights between their curves would be shorter than +// the shortest element, is taken for points along a curve: one curve, compound where it has to, is fitted through them, or one for +// each part of the row where the whole doesn't keep to them; such a vertex is marked 0 then. Fit fits it, by fit_trace_stretch +// where none is given +result ease_corners(design &Design, std::vector &Corners, corner_radii const &Radii, stretch_fit const &Fit = {}); + +// the stretches traced fitted on a thread of their own while the cursor moves: fitting one takes tens of milliseconds. The one asked +// for last is fitted first, the oldest of too many waiting left out, the fits of the last ones kept class stretch_fitter { public: @@ -381,9 +388,9 @@ class stretch_fitter mutable std::mutex m_mutex; std::condition_variable m_wake; - std::optional m_asked; + std::vector m_asked; std::optional m_fitting; - std::optional, trace_stretch>> m_result; + std::vector, trace_stretch>> m_results; // the last done first bool m_stop{false}; std::thread m_thread; }; diff --git a/tests/unit/test_alignment.cpp b/tests/unit/test_alignment.cpp index 8bbb8437..4651eb55 100644 --- a/tests/unit/test_alignment.cpp +++ b/tests/unit/test_alignment.cpp @@ -79,7 +79,8 @@ alignment::design stretch_base(alignment::design const &Line, std::vector corners{1}; + auto const line{alignment::ease_corners(design, corners, {3000.0, 6000.0})}; REQUIRE(line.valid); REQUIRE(line.curves.size() == 1); @@ -97,7 +98,8 @@ TEST_CASE("alignment: a slight bend of straights is eased with two transition cu 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})}; + std::vector corners{1}; + auto const line{alignment::ease_corners(design, corners, {3000.0, 6000.0})}; REQUIRE(line.valid); REQUIRE(line.curves.size() == 1); @@ -110,7 +112,8 @@ TEST_CASE("alignment: a bend too slight for transitions long enough at the radii 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})}; + std::vector corners{1}; + auto const line{alignment::ease_corners(design, corners, {3000.0, 6000.0})}; REQUIRE(line.valid); REQUIRE(line.curves.size() == 1); @@ -126,7 +129,8 @@ TEST_CASE("alignment: a sharp bend between short straights is laid where clicked 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})}; + std::vector corners{1, 1}; + auto const line{alignment::ease_corners(design, corners, {3000.0, 6000.0})}; REQUIRE(line.valid); REQUIRE(line.curves.size() == 2); @@ -245,3 +249,42 @@ TEST_CASE("alignment: a compound curve traced is followed closely") { CHECK(fit.rms < 0.05); CHECK(fit.worst < 0.15); } + +TEST_CASE("alignment: bends clicked close along a curve make one compound curve, not curves with short straights between them") { + // the points clicked in the scenery abc, at 60 km/h: from a straight along z into a curve to the left + alignment::design design; + design.speed = 60.0; + design.shape = alignment::transition_shape::clothoid; + design.start = {-399.8169, 0.0, 419.0340}; + design.end = {-665.1078, 0.0, 1138.5287}; + design.start_direction = glm::normalize(glm::dvec2{0.024027152, 0.999711306}); + design.end_direction = glm::normalize(glm::dvec2{-0.831198797, 0.555975323}); + for (auto const &point : std::vector{{-394.5230, 639.2978}, {-395.0216, 680.5362}, {-400.4562, 753.8415}, {-405.1477, 807.6113}, {-410.4658, 847.8606}, {-419.7362, 885.8522}, + {-435.3280, 924.7706}, {-458.2043, 961.0863}, {-485.8050, 994.2479}, {-531.7847, 1039.0865}, {-565.5341, 1067.6467}, {-637.3993, 1119.9949}}) { + alignment::vertex vertex; + vertex.position = point; + design.vertices.push_back(vertex); + } + design.vertices.front().offset = glm::dot(design.vertices.front().position - glm::dvec2{design.start.x, design.start.z}, design.start_direction); + design.vertices.back().offset = glm::dot(glm::dvec2{design.end.x, design.end.z} - design.vertices.back().position, design.end_direction); + std::vector corners(design.vertices.size(), 1); + + auto const line{alignment::ease_corners(design, corners, {3000.0, 6000.0})}; + + REQUIRE(line.valid); + auto const shortest{std::max(30.0, design.speed)}; + int straights{0}; + for (std::size_t i = 0; i < line.elements.size(); ++i) { + auto const &e{line.elements[i]}; + if (e.kind != alignment::element_kind::straight || e.length < 1e-3) + continue; + ++straights; + // none short between two curves turning the same way + if (i > 0 && i + 1 < line.elements.size() && line.elements[i - 1].turn == line.elements[i + 1].turn) + CHECK(e.length >= shortest - 1e-6); + } + // the straight before the curve and the one after it + CHECK(straights <= 2); + CHECK(line.curves.size() <= 2); + CHECK(design.vertices.size() < 12); +}