From 5764b328bec3b5f3f2a68596d0078d1194045c65 Mon Sep 17 00:00:00 2001 From: maj00r Date: Wed, 21 Oct 2026 18:54:14 +0200 Subject: [PATCH] feat(editor): a curve of the plan traced on the orthophoto, the points clicked along the track give its radius or the radii and the shares of a compound curve --- application/editormode.cpp | 9 + application/editormode.h | 18 +- application/editormoderoute.cpp | 158 +++++++++ application/editormodeselftesttrack.cpp | 88 +++++ application/editormodetrackux.cpp | 22 +- editor/editorAlignment.cpp | 429 ++++++++++++++++++++++++ editor/editorAlignment.hpp | 19 ++ lang/editor_pl.po | 54 +++ 8 files changed, 795 insertions(+), 2 deletions(-) diff --git a/application/editormode.cpp b/application/editormode.cpp index 886546bb..7ed6c292 100644 --- a/application/editormode.cpp +++ b/application/editormode.cpp @@ -6465,6 +6465,8 @@ void editor_mode::on_key(int const Key, int const Scancode, int const Action, in { if (m_tape.active) m_tape = {}; + else if (m_route.trace.active && route_active()) + route_trace_end(false); else if (m_lay.active && false == m_lay.points.empty()) lay_cancel(); else if (m_track_style.picking) @@ -6493,12 +6495,19 @@ void editor_mode::on_key(int const Key, int const Scancode, int const Action, in m_fill_points.pop_back(); if (is_press(Action) && ui()->mode() == nodebank_panel::TRACK && m_lay.active && !m_lay.points.empty()) lay_back(); + if (is_press(Action) && m_route.trace.active && route_active() && !m_route.trace.points.empty()) + { + m_route.trace.points.pop_back(); + route_trace_fit(); + } break; case GLFW_KEY_ENTER: case GLFW_KEY_KP_ENTER: if (is_press(Action) && ui()->mode() == nodebank_panel::TRACK && m_lay.active) lay_enter(); + if (is_press(Action) && m_route.trace.active && route_active()) + route_trace_end(true); break; case GLFW_KEY_K: diff --git a/application/editormode.h b/application/editormode.h index c9053ffd..2b6bcdfa 100644 --- a/application/editormode.h +++ b/application/editormode.h @@ -926,6 +926,16 @@ class editor_mode : public application_mode, private editor_track::observer, pri float metres{1.0f}; // pixels of a metre and of a unit of the curvature, kept as the drag goes float curvature{1.0f}; } band; + // points clicked along the middle of the track seen on the orthophoto, the radii of the curves and their shares fitted to them after each click + struct trace_state + { + bool active{false}; + int radii{0}; // 0: as many as the points need, up to three + std::vector points; + alignment::design before; // taken back by Esc + double rms{0.0}; + double worst{0.0}; + } trace; }; void render_route_ui(); bool render_route_band(); @@ -933,6 +943,12 @@ class editor_mode : public application_mode, private editor_track::observer, pri bool render_route_vertices(); bool render_route_vertex(int const Index); void render_route_result(); + void render_route_trace(); + void route_trace_start(); + void route_trace_fit(); + void route_trace_end(bool const Keep); + int route_trace_hit() const; + glm::dvec3 route_trace_cursor() const; void route_reset(); void route_bind_ends(); bool route_from_curve(TTrack &Track); @@ -1586,7 +1602,7 @@ class editor_mode : public application_mode, private editor_track::observer, pri void draw_track_hover() const; struct track_intent { - enum class kind { none, gizmo, select, point, box, straight_set, straight_break, detour, straight_tool, straight_handle, route_vertex, route_grip, switch_at, lay, lay_end, sweep, turntable_place, turntable_exit, turntable_fit, spread, extend_end, signal_place, signal_move, switch_slide, pick_style, fouling, measure }; + enum class kind { none, gizmo, select, point, box, straight_set, straight_break, detour, straight_tool, straight_handle, route_vertex, route_grip, switch_at, lay, lay_end, sweep, turntable_place, turntable_exit, turntable_fit, spread, extend_end, signal_place, signal_move, switch_slide, pick_style, fouling, measure, route_trace }; kind what{kind::none}; TTrack *track{nullptr}; int path{0}; diff --git a/application/editormoderoute.cpp b/application/editormoderoute.cpp index 9fd4f5b5..d4c2419b 100644 --- a/application/editormoderoute.cpp +++ b/application/editormoderoute.cpp @@ -83,6 +83,7 @@ void editor_mode::route_reset() route.vertex = -1; route.grip = -1; route.result = {}; + route.trace = {}; if (route.from == nullptr || route.to == nullptr) return; if (false == editor_track::find_chain(route.from, route.to, route.chain, route.error, true)) @@ -1086,6 +1087,162 @@ void editor_mode::draw_route_overlay() const } drawlist->AddText(ImVec2(screen.x + 9.0f, screen.y - 18.0f), overlay_color::marked, format("W%d", i + 1).c_str()); } + if (route.trace.points.empty()) + return; + std::vector plan; + for (auto const &point : route.trace.points) + plan.push_back(plan_of(point)); + auto const projections{alignment::project(result, plan)}; + auto const hit{route_trace_hit()}; + for (std::size_t i = 0; i < route.trace.points.size(); ++i) + { + ImVec2 screen; + if (false == projection.project(route.trace.points[i], screen)) + continue; + auto const departure{i < projections.size() ? std::abs(projections[i].offset) : 0.0}; + auto const colour{departure <= 0.1 ? IM_COL32(90, 220, 110, 255) : departure <= 0.3 ? IM_COL32(255, 200, 60, 255) : IM_COL32(255, 90, 70, 255)}; + drawlist->AddCircle(screen, 5.0f, IM_COL32(0, 0, 0, 200), 12, 3.5f); + drawlist->AddCircle(screen, 5.0f, colour, 12, 2.0f); + if (static_cast(i) == hit) + drawlist->AddCircle(screen, 10.0f, overlay_color::highlight, 16, 2.0f); + if (i < projections.size()) + drawlist->AddText(ImVec2(screen.x + 8.0f, screen.y - 7.0f), colour, format("%.0f cm", departure * 100.0).c_str()); + } +} + +void editor_mode::render_route_trace() +{ + auto &trace{m_route.trace}; + if (false == trace.active) + { + if (ImGui::Button(STR_C("Trace the curve on the photo"))) + route_trace_start(); + item_tooltip("Click points along the middle of the track on the orthophoto: after each click the curve takes the radius, or the radii of a compound curve and their shares, which follow the points best. The straights and the transition curves stay"); + return; + } + ImGui::SeparatorText(STR_C("Tracing the curve")); + if (segmented("traceradii", trace.radii, + {{0, STR_C("Auto"), "As few radii as follow the points, up to three"}, {1, STR_C("1 R"), "One radius"}, {2, STR_C("2 R"), "A compound curve of two radii"}, {3, STR_C("3 R"), "A compound curve of three radii"}})) + route_trace_fit(); + if (trace.points.empty()) + { + ImGui::PushTextWrapPos(0.0f); + ImGui::TextDisabled("%s", STR_C("Click points along the middle of the track on the photo, spread over the whole curve; a few on the straights at its ends help too")); + ImGui::PopTextWrapPos(); + } + else + { + auto const colour = [](double const Departure) { return Departure <= 0.1 ? IM_COL32(90, 220, 110, 255) : Departure <= 0.3 ? IM_COL32(255, 200, 60, 255) : IM_COL32(255, 90, 70, 255); }; + metric_tiles({{format("%zu", trace.points.size()), STR_C("points"), ImGui::GetColorU32(ImGuiCol_Text)}, + {format("%.0f cm", trace.rms * 100.0), STR_C("departure"), colour(trace.rms)}, + {format("%.0f cm", trace.worst * 100.0), STR_C("at most"), colour(trace.worst)}}); + } + if (main_button(STR_C("Apply Enter"), ImVec2(-1.0f, 0.0f))) + route_trace_end(true); + if (ImGui::Button(STR_C("Drop Esc"))) + route_trace_end(false); + ImGui::SameLine(); + ImGui::BeginDisabled(trace.points.empty()); + if (ImGui::Button(STR_C("Take out the last point Backspace"))) + { + trace.points.pop_back(); + route_trace_fit(); + } + ImGui::EndDisabled(); +} + +void editor_mode::route_trace_start() +{ + auto &trace{m_route.trace}; + trace.active = true; + trace.points.clear(); + trace.before = m_route.design; + trace.rms = 0.0; + trace.worst = 0.0; +} + +// every point goes to the curve nearest to it along the line, the curves without points keep their shape +void editor_mode::route_trace_fit() +{ + auto &route{m_route}; + auto &trace{route.trace}; + auto design{route.design}; + trace.rms = 0.0; + trace.worst = 0.0; + if (trace.points.empty() || design.vertices.empty() || false == route.result.valid) + return; + std::vector plan; + for (auto const &point : trace.points) + plan.push_back(plan_of(point)); + auto const projections{alignment::project(route.result, plan)}; + std::vector> owned(design.vertices.size()); + for (std::size_t i = 0; i < plan.size(); ++i) + { + std::size_t nearest{0}; + for (std::size_t k = 1; k < route.result.vertex_chainages.size() && k < owned.size(); ++k) + if (std::abs(route.result.vertex_chainages[k] - projections[i].chainage) < std::abs(route.result.vertex_chainages[nearest] - projections[i].chainage)) + nearest = k; + owned[nearest].push_back(plan[i]); + } + for (std::size_t k = 0; k < owned.size(); ++k) + { + if (owned[k].empty()) + continue; + auto const fit{alignment::fit_trace(design, k, owned[k], trace.radii)}; + if (fit.valid) + design.vertices[k] = fit.fitted; + } + route.design = design; + route_update(); + double sum{0.0}; + for (auto const &projection : alignment::project(route.result, plan)) + { + sum += projection.offset * projection.offset; + trace.worst = std::max(trace.worst, std::abs(projection.offset)); + } + trace.rms = std::sqrt(sum / static_cast(plan.size())); +} + +void editor_mode::route_trace_end(bool const Keep) +{ + auto &trace{m_route.trace}; + if (false == trace.active) + return; + auto const traced{false == trace.points.empty()}; + if (false == Keep) + { + m_route.design = trace.before; + route_update(); + } + trace = {}; + if (Keep && traced) + route_apply(); +} + +// the photo under the cursor where the depth finds it, else the level of the design there +glm::dvec3 editor_mode::route_trace_cursor() const +{ + if (GfxRenderer->Mouse_Hit()) + return cursor_ground(); + auto const &result{m_route.result}; + auto point{cursor_level(alignment::elevation(result, result.length * 0.5))}; + if (auto const projections{alignment::project(result, {plan_of(point)})}; false == projections.empty()) + point = cursor_level(alignment::elevation(result, projections.front().chainage)); + return point; +} + +int editor_mode::route_trace_hit() const +{ + screen_projection const projection; + float best{kHandleRadius * kHandleRadius}; + int hit{-1}; + for (int i = 0; i < static_cast(m_route.trace.points.size()); ++i) + if (auto const distance{projection.mouse_distance2(m_route.trace.points[i])}; distance < best) + { + best = distance; + hit = i; + } + return hit; } void editor_mode::render_route_gizmo() @@ -1166,6 +1323,7 @@ void editor_mode::render_route_ui() } item_tooltip("Reads the curves of the line from the scenery again, dropping the changes not applied"); auto changed{render_route_band()}; + render_route_trace(); render_route_stored(); changed |= render_route_vertices(); if (ImGui::CollapsingHeader(STR_C("Design speed and limits"))) diff --git a/application/editormodeselftesttrack.cpp b/application/editormodeselftesttrack.cpp index 845e9fd2..05ca9361 100644 --- a/application/editormodeselftesttrack.cpp +++ b/application/editormodeselftesttrack.cpp @@ -263,6 +263,94 @@ editor_mode::selftest_result editor_mode::selftest_track(std::string const &comm WriteLog(format("SELFTEST %s: valid %d, %zu elements, extension %.3f / %.3f", command.c_str(), m_route.result.valid ? 1 : 0, m_route.result.elements.size(), m_route.result.start_extension, m_route.result.end_extension) + errors); return selftest_result::handled; } + if (command == "routeplain") + { + int index{0}; + double radius{1000.0}; + words >> index >> radius; + if (index >= 0 && index < static_cast(m_route.design.vertices.size())) + { + auto &vertex{m_route.design.vertices[index]}; + vertex.radius = radius; + vertex.compound = false; + vertex.share = 1.0; + vertex.arcs.clear(); + route_update(); + } + WriteLog(format("SELFTEST routeplain: valid %d", m_route.result.valid ? 1 : 0)); + return selftest_result::handled; + } + // points along the curves of the design as it is now, every Spacing metres from 20 m before the first curve to 20 m past the last, + // shifted across by up to Noise to and fro, as the trace takes them + if (command == "tracesample") + { + double spacing{20.0}; + double noise{0.0}; + words >> spacing >> noise; + route_trace_start(); + auto const &result{m_route.result}; + auto from{result.length}; + auto to{0.0}; + for (auto const &e : result.elements) + if (e.kind != alignment::element_kind::straight) + { + from = std::min(from, e.chainage); + to = std::max(to, e.chainage + e.length); + } + int count{0}; + for (auto chainage{std::max(0.0, from - 20.0)}; chainage <= std::min(result.length, to + 20.0); chainage += spacing, ++count) + { + auto const station{alignment::evaluate(result, chainage)}; + glm::dvec3 const across{-station.direction.z, 0.0, station.direction.x}; + auto const shift{noise * (count % 3 == 0 ? 1.0 : count % 3 == 1 ? -0.6 : 0.2)}; + m_route.trace.points.push_back(station.position + glm::normalize(across) * shift); + } + auto const &middle{m_route.trace.points.empty() ? glm::dvec3{0.0} : m_route.trace.points[m_route.trace.points.size() / 2]}; + WriteLog(format("SELFTEST tracesample: %zu points over %.1f m, the middle one at %.3f %.3f", m_route.trace.points.size(), to - from, middle.x, middle.z)); + return selftest_result::handled; + } + if (command == "traceradii" || command == "traceclick" || command == "traceundo") + { + auto const started{std::chrono::steady_clock::now()}; + if (command == "traceradii") + { + words >> m_route.trace.radii; + route_trace_fit(); + } + else if (command == "traceundo") + { + if (false == m_route.trace.points.empty()) + m_route.trace.points.pop_back(); + route_trace_fit(); + } + else + { + track_intent intent; + words >> intent.position.x >> intent.position.z; + intent.what = track_intent::kind::route_trace; + intent.index = -1; + track_press(intent); + } + auto const took{std::chrono::duration(std::chrono::steady_clock::now() - started).count()}; + WriteLog(format("SELFTEST %s: %zu points, departure %.3f m, at most %.3f m, valid %d, %.0f ms", command.c_str(), m_route.trace.points.size(), m_route.trace.rms, m_route.trace.worst, m_route.result.valid ? 1 : 0, took)); + for (std::size_t i = 0; i < m_route.design.vertices.size(); ++i) + { + auto const &v{m_route.design.vertices[i]}; + auto text{format("SELFTEST %s vertex %zu: R %.0f L %.1f / %.1f compound %d share %.3f", command.c_str(), i, v.radius, v.transition_in, v.transition_out, v.compound ? 1 : 0, v.share)}; + for (auto const &arc : v.arcs) + text += format(" [R %.0f L %.1f share %.3f]", arc.radius, arc.transition, arc.share); + WriteLog(text); + } + return selftest_result::handled; + } + if (command == "traceend") + { + int keep{1}; + words >> keep; + route_trace_end(keep != 0); + WriteLog(format("SELFTEST traceend: active %d, valid %d, ", m_route.trace.active ? 1 : 0, m_route.result.valid ? 1 : 0) + m_route.status + " " + m_route.error); + return selftest_result::handled; + } if (command == "routekm") { double km{0.0}; diff --git a/application/editormodetrackux.cpp b/application/editormodetrackux.cpp index 0869faba..f40348d3 100644 --- a/application/editormodetrackux.cpp +++ b/application/editormodetrackux.cpp @@ -549,6 +549,18 @@ editor_mode::track_intent editor_mode::track_intent_at(int const Mods) return make(kind::straight_handle, hit == 2 ? "Drag: shift the straight sideways" : "Drag: move this end, the curves at it follow"); } } + if (m_route.trace.active && route_active()) + { + intent.index = route_trace_hit(); + if (control && intent.index >= 0) + { + intent.position = m_route.trace.points[intent.index]; + return make(kind::route_trace, "Click: take this point out"); + } + intent.index = -1; + intent.position = route_trace_cursor(); + return make(kind::route_trace, m_route.trace.points.empty() ? "Click: a point in the middle of the track on the photo. Enter: apply, Esc: drop" : "Click: one more point, the curve follows. Ctrl+click: take a point out"); + } if (int vertex, grip; route_hit(vertex, grip)) { intent.index = vertex >= 0 ? vertex : grip; @@ -711,6 +723,13 @@ void editor_mode::track_press(track_intent const &Intent) m_tape.points.clear(); m_tape.points.push_back(m_tape.mouse); return; + case kind::route_trace: + if (Intent.index >= 0) + m_route.trace.points.erase(m_route.trace.points.begin() + Intent.index); + else + m_route.trace.points.push_back(Intent.position); + route_trace_fit(); + return; case kind::lay_end: if (m_lay.points.empty() && start_extend(Intent.snap)) return; @@ -849,7 +868,8 @@ void editor_mode::draw_track_intent() const case kind::extend_end: case kind::turntable_fit: drawlist->AddCircleFilled(at, 5.0f, IM_COL32(255, 60, 255, 255)); break; case kind::turntable_exit: drawlist->AddCircle(at, 11.0f, overlay_color::grip, 20, 2.0f); break; - case kind::lay: drawlist->AddCircle(at, 6.0f, IM_COL32(255, 255, 255, 200), 12, 1.5f); break; + case kind::lay: + case kind::route_trace: drawlist->AddCircle(at, 6.0f, IM_COL32(255, 255, 255, 200), 12, 1.5f); break; default: break; } } diff --git a/editor/editorAlignment.cpp b/editor/editorAlignment.cpp index 232bb9f4..dfc86eb2 100644 --- a/editor/editorAlignment.cpp +++ b/editor/editorAlignment.cpp @@ -1237,6 +1237,435 @@ station evaluate(result const &Result, double const Chainage) return result; } +namespace +{ + +// the plan from From to To in chords short enough to keep to the curves, the directions of its ends carried on far, so that a point +// past them is measured across the track +struct plan_line +{ + std::vector points; + std::vector chainages; +}; + +plan_line sample_plan(result const &Result, double const From, double const To, double const Step) +{ + constexpr double kFar{5000.0}; + plan_line line; + if (Result.elements.empty()) + return line; + auto const from{std::clamp(From, 0.0, Result.length)}; + auto const to{std::clamp(To, from, Result.length)}; + auto const push = [&](glm::dvec2 const &Position, double const Chainage) { + line.points.push_back(Position); + line.chainages.push_back(Chainage); + }; + for (auto const &e : Result.elements) + { + auto const low{std::max(from, e.chainage)}; + auto const high{std::min(to, e.chainage + e.length)}; + if (high < low) + continue; + auto const count{std::max(1, static_cast(std::ceil((high - low) / Step)))}; + auto const at = [&](int const Index) { return low + (high - low) * Index / count; }; + if (e.kind != element_kind::spiral) + { + glm::dvec2 position; + glm::dvec2 direction; + for (int i = 0; i <= count; ++i) + { + element_point(e, Result.shape, at(i) - e.chainage, position, direction); + push(position, at(i)); + } + continue; + } + // spiral_point integrates from the start of the element for every point, here it goes on from the previous one + auto const heading = [&](double const S) { return e.curvature_start * S + (e.curvature_end - e.curvature_start) * S * S / (2.0 * std::max(e.length, 1e-9)); }; + double x{0.0}; + double y{0.0}; + double s{0.0}; + auto const advance = [&](double const Distance) { + auto const pieces{std::max(1, static_cast(std::ceil((Distance - s) / 2.0)))}; + auto const h{(Distance - s) / pieces}; + for (int i = 0; i < pieces; ++i) + { + auto const a{heading(s + h * i)}; + auto const m{heading(s + h * (i + 0.5))}; + auto const b{heading(s + h * (i + 1))}; + x += h / 6.0 * (std::cos(a) + 4.0 * std::cos(m) + std::cos(b)); + y += h / 6.0 * (std::sin(a) + 4.0 * std::sin(m) + std::sin(b)); + } + s = Distance; + }; + for (int i = 0; i <= count; ++i) + { + advance(at(i) - e.chainage); + push(e.origin + e.direction * x + e.normal * y, at(i)); + } + } + if (line.points.size() < 2) + return line; + auto const away = [](glm::dvec2 const &End, glm::dvec2 const &Inner) { + auto const offset{End - Inner}; + return glm::length(offset) > 1e-9 ? offset / glm::length(offset) : glm::dvec2{0.0}; + }; + std::size_t first{1}; + while (first + 1 < line.points.size() && glm::distance(line.points[first], line.points.front()) < 1e-6) + ++first; + auto last{line.points.size() - 2}; + while (last > 0 && glm::distance(line.points[last], line.points.back()) < 1e-6) + --last; + auto const before{away(line.points.front(), line.points[first])}; + auto const after{away(line.points.back(), line.points[last])}; + line.points.insert(line.points.begin(), line.points.front() + before * kFar); + line.chainages.insert(line.chainages.begin(), line.chainages.front() - kFar); + push(line.points.back() + after * kFar, line.chainages.back() + kFar); + return line; +} + +plan_projection project_on(plan_line const &Line, glm::dvec2 const &Point) +{ + plan_projection best; + auto nearest{std::numeric_limits::max()}; + for (std::size_t i = 0; i + 1 < Line.points.size(); ++i) + { + auto const segment{Line.points[i + 1] - Line.points[i]}; + auto const length2{glm::dot(segment, segment)}; + if (length2 < 1e-18) + continue; + auto const t{std::clamp(glm::dot(Point - Line.points[i], segment) / length2, 0.0, 1.0)}; + auto const offset{Point - (Line.points[i] + segment * t)}; + auto const distance2{glm::dot(offset, offset)}; + if (distance2 >= nearest) + continue; + nearest = distance2; + best.chainage = Line.chainages[i] + (Line.chainages[i + 1] - Line.chainages[i]) * t; + best.offset = std::copysign(std::sqrt(distance2), cross(segment, offset)); + } + return best; +} + +double squares(std::vector const &Values) +{ + double sum{0.0}; + for (auto const value : Values) + sum += value * value; + return sum; +} + +using residuals_of = std::function const &, std::vector &)>; + +// Levenberg-Marquardt: the parameters of the least sum of the squares of the residuals; a step to where Residuals fails isn't taken +double least_squares(std::vector &Parameters, residuals_of const &Residuals) +{ + constexpr double kDifference{1e-4}; + std::vector residuals; + if (false == Residuals(Parameters, residuals)) + return std::numeric_limits::max(); + auto cost{squares(residuals)}; + auto const n{Parameters.size()}; + auto const m{residuals.size()}; + double damping{1e-3}; + std::vector trial; + for (int iteration = 0; iteration < 40 && cost > 1e-12; ++iteration) + { + std::vector> jacobian(n, std::vector(m, 0.0)); + for (std::size_t j = 0; j < n; ++j) + { + auto parameters{Parameters}; + parameters[j] += kDifference; + auto sign{1.0}; + if (false == Residuals(parameters, trial) || trial.size() != m) + { + parameters[j] = Parameters[j] - kDifference; + sign = -1.0; + if (false == Residuals(parameters, trial) || trial.size() != m) + continue; + } + for (std::size_t i = 0; i < m; ++i) + jacobian[j][i] = sign * (trial[i] - residuals[i]) / kDifference; + } + std::vector normal(n * n, 0.0); + std::vector gradient(n, 0.0); + for (std::size_t j = 0; j < n; ++j) + { + for (std::size_t i = 0; i < m; ++i) + gradient[j] += jacobian[j][i] * residuals[i]; + for (std::size_t k = 0; k < n; ++k) + for (std::size_t i = 0; i < m; ++i) + normal[j * n + k] += jacobian[j][i] * jacobian[k][i]; + } + auto better{false}; + auto settled{false}; + for (int attempt = 0; attempt < 12 && false == better; ++attempt) + { + auto system{normal}; + std::vector step(n); + for (std::size_t j = 0; j < n; ++j) + { + system[j * n + j] += damping * std::max(normal[j * n + j], 1e-9); + step[j] = -gradient[j]; + } + for (std::size_t c = 0; c < n; ++c) + { + auto pivot{c}; + for (std::size_t r = c + 1; r < n; ++r) + if (std::abs(system[r * n + c]) > std::abs(system[pivot * n + c])) + pivot = r; + if (std::abs(system[pivot * n + c]) < 1e-300) + break; + for (std::size_t k = 0; k < n; ++k) + std::swap(system[c * n + k], system[pivot * n + k]); + std::swap(step[c], step[pivot]); + for (std::size_t r = c + 1; r < n; ++r) + { + auto const factor{system[r * n + c] / system[c * n + c]}; + for (std::size_t k = c; k < n; ++k) + system[r * n + k] -= factor * system[c * n + k]; + step[r] -= factor * step[c]; + } + } + for (std::size_t c = n; c-- > 0;) + { + for (std::size_t k = c + 1; k < n; ++k) + step[c] -= system[c * n + k] * step[k]; + step[c] = std::abs(system[c * n + c]) > 1e-300 ? step[c] / system[c * n + c] : 0.0; + } + auto candidate{Parameters}; + for (std::size_t j = 0; j < n; ++j) + candidate[j] += step[j]; + if (Residuals(candidate, trial) && trial.size() == m && squares(trial) < cost) + { + auto const lower{squares(trial)}; + settled = cost - lower < 1e-6 * cost; + Parameters = candidate; + residuals = trial; + cost = lower; + damping = std::max(1e-9, damping * 0.3); + better = true; + } + else + { + damping *= 8.0; + } + } + if (false == better || settled) + break; + } + return cost; +} + +constexpr double kTraceRadiusMin{15.0}; +constexpr double kTraceRadiusMax{50000.0}; +constexpr double kTraceRobust{0.25}; + +// 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(vertex const &Base, std::vector const &Parameters, int const Radii) +{ + auto shape{Base}; + auto const radius = [&](int const Index) { return std::exp(std::clamp(Parameters[Index], std::log(kTraceRadiusMin), std::log(kTraceRadiusMax))); }; + shape.radius = radius(0); + shape.share = 1.0; + if (Radii == 1) + { + shape.compound = false; + return shape; + } + std::vector weights{1.0}; + for (int i = 1; i < Radii; ++i) + weights.push_back(std::exp(std::clamp(Parameters[Radii + i - 1], -5.0, 5.0))); + double total{0.0}; + for (auto const weight : weights) + total += weight; + shape.compound = true; + shape.share = weights[0] / total; + shape.arcs.clear(); + for (int i = 1; i < Radii; ++i) + { + compound_arc arc; + arc.radius = radius(i); + arc.transition = i - 1 < static_cast(Base.arcs.size()) ? Base.arcs[i - 1].transition : 0.0; + arc.share = weights[i] / total; + shape.arcs.push_back(arc); + } + return shape; +} + +} // namespace + +std::vector project(result const &Result, std::vector const &Points) +{ + std::vector projections; + if (Result.elements.empty()) + return projections; + auto const line{sample_plan(Result, 0.0, Result.length, std::max(0.5, Result.length / 4000.0))}; + for (auto const &point : Points) + projections.push_back(project_on(line, point)); + return projections; +} + +trace_fit fit_trace(design const &Design, std::size_t const Vertex, std::vector const &Points, int const Radii) +{ + trace_fit best; + if (Vertex >= Design.vertices.size() || Points.empty()) + return best; + auto const &base{Design.vertices[Vertex]}; + auto const residuals = [&](int const Count) -> residuals_of { + return [&, Count](std::vector const &Parameters, std::vector &Residuals) { + auto trial{Design}; + trial.norms.checked = false; + trial.vertices[Vertex] = traced_shape(base, Parameters, Count); + auto const result{compute(trial)}; + if (false == result.valid) + return false; + auto from{std::numeric_limits::max()}; + auto to{std::numeric_limits::lowest()}; + for (auto const &e : result.elements) + { + if (e.vertex != static_cast(Vertex)) + continue; + from = std::min(from, e.chainage); + to = std::max(to, e.chainage + e.length); + } + if (from > to) + return false; + auto const line{sample_plan(result, from - 1.0, to + 1.0, std::clamp((to - from) / 300.0, 0.25, 2.0))}; + Residuals.clear(); + for (auto const &point : Points) + { + // a point clicked by mistake far off pulls the curve less than the square of its departure would + 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; + }; + }; + auto const cost_of = [](residuals_of const &Residuals, std::vector const &Parameters) { + std::vector values; + return Residuals(Parameters, values) ? squares(values) : std::numeric_limits::max(); + }; + + // one radius: a scan over the whole range, then the golden section around the best of it + auto const one{residuals(1)}; + auto const low{std::log(kTraceRadiusMin)}; + auto const high{std::log(kTraceRadiusMax)}; + constexpr int kScan{60}; + auto const spacing{(high - low) / kScan}; + auto bestlog{std::log(std::max(kTraceRadiusMin, base.radius))}; + auto bestcost{cost_of(one, {bestlog})}; + for (int i = 0; i <= kScan; ++i) + if (auto const cost{cost_of(one, {low + spacing * i})}; cost < bestcost) + { + bestcost = cost; + bestlog = low + spacing * i; + } + if (bestcost == std::numeric_limits::max()) + return best; + { + constexpr double kGolden{0.6180339887498949}; + auto a{bestlog - spacing}; + auto b{bestlog + spacing}; + auto c{b - (b - a) * kGolden}; + auto d{a + (b - a) * kGolden}; + auto fc{cost_of(one, {c})}; + auto fd{cost_of(one, {d})}; + for (int i = 0; i < 50; ++i) + { + if (fc < fd) + { + b = d; + d = c; + fd = fc; + c = b - (b - a) * kGolden; + fc = cost_of(one, {c}); + } + else + { + a = c; + c = d; + fc = fd; + d = a + (b - a) * kGolden; + fd = cost_of(one, {d}); + } + } + if (auto const middle{(a + b) * 0.5}; cost_of(one, {middle}) < bestcost) + bestlog = middle; + } + + struct candidate + { + std::vector parameters; + double cost{std::numeric_limits::max()}; + }; + std::array fits; + fits[1] = {{bestlog}, cost_of(one, {bestlog})}; + auto const most{std::clamp(Radii > 0 ? Radii : 3, 1, std::max(1, static_cast(Points.size()) / 2))}; + auto const shift{std::log(1.3)}; + auto const share = [](double const Second, double const First) { return std::log(Second / First); }; + auto const rms = [&](int const Count) { return std::sqrt(fits[Count].cost / static_cast(Points.size())); }; + // the next radius has to halve the departure about, a curve which keeps to the points already gets none + constexpr double kClose{0.08}; + auto chosen{1}; + auto const wanted = [&](int const Count) { return Count <= most && (Radii > 0 || rms(chosen) > kClose); }; + auto const take = [&](int const Count) { + if (fits[Count].cost == std::numeric_limits::max()) + return; + if (Radii > 0 || rms(Count) < 0.6 * rms(chosen)) + chosen = Count; + }; + if (wanted(2)) + { + auto const two{residuals(2)}; + for (auto const fraction : {0.3, 0.5, 0.7}) + for (auto const sign : {-1.0, 1.0}) + { + std::vector parameters{bestlog + sign * shift, bestlog - sign * shift, share(1.0 - fraction, fraction)}; + if (auto const cost{least_squares(parameters, two)}; cost < fits[2].cost) + fits[2] = {parameters, cost}; + } + take(2); + } + if (wanted(3)) + { + auto const three{residuals(3)}; + std::vector> starts; + for (auto const sign : {-1.0, 1.0}) + starts.push_back({bestlog + sign * shift, bestlog - sign * shift, bestlog + sign * shift, 0.0, 0.0}); + if (false == fits[2].parameters.empty()) + { + auto const &two{fits[2].parameters}; + auto const second{std::exp(two[2])}; + starts.push_back({two[0], two[0], two[1], 0.0, std::log(2.0 * second)}); + starts.push_back({two[0], two[1], two[1], std::log(0.5 * second), std::log(0.5 * second)}); + } + for (auto parameters : starts) + if (auto const cost{least_squares(parameters, three)}; cost < fits[3].cost) + fits[3] = {parameters, cost}; + take(3); + } + + auto shape{traced_shape(base, fits[chosen].parameters, chosen)}; + shape.radius = std::max(kTraceRadiusMin, std::round(shape.radius)); + for (auto &arc : shape.arcs) + arc.radius = std::max(kTraceRadiusMin, std::round(arc.radius)); + auto trial{Design}; + trial.vertices[Vertex] = shape; + auto const result{compute(trial)}; + if (false == result.valid) + return best; + best.valid = true; + best.fitted = shape; + double sum{0.0}; + for (auto const &projection : project(result, Points)) + { + sum += projection.offset * projection.offset; + best.worst = std::max(best.worst, std::abs(projection.offset)); + } + best.rms = std::sqrt(sum / static_cast(Points.size())); + return best; +} + std::size_t minimum_pieces(result const &Result, design const &Design) { std::size_t count{0}; diff --git a/editor/editorAlignment.hpp b/editor/editorAlignment.hpp index 48cefcc5..eaba18a0 100644 --- a/editor/editorAlignment.hpp +++ b/editor/editorAlignment.hpp @@ -187,6 +187,25 @@ station evaluate(result const &Result, double const Chainage); double elevation(result const &Result, double const Chainage); double grade(result const &Result, double const Chainage); +// a point of the plan projected on the alignment, its offset positive to the left; past the ends it's measured from the straights carried on +struct plan_projection +{ + double chainage{0.0}; + double offset{0.0}; +}; +std::vector project(result const &Result, std::vector const &Points); + +struct trace_fit +{ + bool valid{false}; + vertex fitted; + double rms{0.0}; + double worst{0.0}; +}; +// the radii of the curve of the vertex and their shares which follow the points traced along the track best; the tangents, the +// transitions and the cant stay. Radii 0 takes the fewest radii, up to three, which follow the points about as well as more would +trace_fit fit_trace(design const &Design, std::size_t const Vertex, std::vector const &Points, int const Radii); + std::size_t minimum_pieces(result const &Result, design const &Design); // Breaks: ascending chainages at which pieces have to end. Counts: least number of the pieces between the consecutive breaks, // one more than the breaks; Intervals receives the number of the interval of each piece diff --git a/lang/editor_pl.po b/lang/editor_pl.po index d1d66e7f..77d020a8 100644 --- a/lang/editor_pl.po +++ b/lang/editor_pl.po @@ -10168,3 +10168,57 @@ msgstr "Kładzenie nowych torów" msgid "Ask as the laying starts" msgstr "Pytaj przy rozpoczęciu kładzenia" + +msgid "Trace the curve on the photo" +msgstr "Odrysuj łuk ze zdjęcia" + +msgid "Click points along the middle of the track on the orthophoto: after each click the curve takes the radius, or the radii of a compound curve and their shares, which follow the points best. The straights and the transition curves stay" +msgstr "Klikaj punkty po środku toru na ortofotomapie: po każdym kliknięciu łuk dobiera promień, albo promienie łuku koszowego i ich udziały, tak żeby szedł po punktach. Proste i krzywe przejściowe zostają" + +msgid "Tracing the curve" +msgstr "Odrysowywanie łuku" + +msgid "Auto" +msgstr "Auto" + +msgid "As few radii as follow the points, up to three" +msgstr "Najmniej promieni, przy których łuk trzyma się punktów, najwyżej trzy" + +msgid "One radius" +msgstr "Jeden promień" + +msgid "A compound curve of two radii" +msgstr "Łuk koszowy z dwóch promieni" + +msgid "A compound curve of three radii" +msgstr "Łuk koszowy z trzech promieni" + +msgid "Click points along the middle of the track on the photo, spread over the whole curve; a few on the straights at its ends help too" +msgstr "Klikaj punkty po środku toru na zdjęciu, rozłożone na cały łuk; kilka na prostych przy jego końcach też pomaga" + +msgid "points" +msgstr "punktów" + +msgid "departure" +msgstr "odchyłka" + +msgid "at most" +msgstr "największa" + +msgid "Apply Enter" +msgstr "Zastosuj Enter" + +msgid "Drop Esc" +msgstr "Porzuć Esc" + +msgid "Take out the last point Backspace" +msgstr "Usuń ostatni punkt Backspace" + +msgid "Click: take this point out" +msgstr "LPM: usuń ten punkt" + +msgid "Click: a point in the middle of the track on the photo. Enter: apply, Esc: drop" +msgstr "LPM: punkt na środku toru na zdjęciu. Enter: zastosuj, Esc: porzuć" + +msgid "Click: one more point, the curve follows. Ctrl+click: take a point out" +msgstr "LPM: kolejny punkt, łuk idzie za nimi. Ctrl+LPM: usuń punkt"