mirror of
https://github.com/MaSzyna-EU07/maszyna.git
synced 2026-10-20 18:34:53 +02:00
ff1710d9ea
Labels: - every button, field, tooltip and message of the editor reread for a beginner: short nouns for the fields with the sentence in the tooltip, the less used options under More settings, one glossary (field of work, template, track, grade), sides told from directions, degrees and per mille as signs, messages built from one format string - plural forms through STR_P; the last entry of a po file is read too, and the Polish plural rule is in the header - lang/editor_pl.po: the code pass entries merged, 391 entries no code uses and 3 duplicates dropped - roads: the notes and errors naming a road translated; the build options and the layout fold the texture length, kerb sizes, traffic share, roundabout part, height over the ground and junction speed limit under More settings, the lanes start folded Grade lines fitted to the terrain model (NMT), a line or the whole scenery at once: - fit_ground: a robust fit on nodes 5 m apart; the ground parting from the line (a bridge, a viaduct, a wall) is gone past rather than followed, the bends keep off the switches as far as their curves reach, a run of bends the same way is one point of intersection, a curve with no room is as sharp as the room lets it be - the ends of the line and the other paths of its switches leading into lines fitted already bound the line by what the tracks between can climb; ends meeting in the plan, laid apart in height by hand, are joined before the run - apply_profile: the paths are split where split_path takes the parameter (a straight was cut in the wrong place, into pieces of a few centimetres the game takes for zero length), no piece shorter than the margin; the paths off the ends of the route are brought to its height - Bytom: 1600 lines, grade breaks at the joints 89 -> 20; Tarnowskie Gory: 848 lines, 223 -> 3; no step over 2 cm
896 lines
35 KiB
C++
896 lines
35 KiB
C++
/*
|
|
This Source Code Form is subject to the
|
|
terms of the Mozilla Public License, v.
|
|
2.0. If a copy of the MPL was not
|
|
distributed with this file, You can
|
|
obtain one at
|
|
http://mozilla.org/MPL/2.0/.
|
|
*/
|
|
|
|
#include "stdafx.h"
|
|
#include "editor/editorProfile.hpp"
|
|
#include "editor/editorFormat.hpp"
|
|
#include "utilities/translation.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <cstdio>
|
|
#include <optional>
|
|
#include <tuple>
|
|
|
|
namespace profile
|
|
{
|
|
|
|
namespace
|
|
{
|
|
|
|
double polyline(line const &Line, double const Chainage)
|
|
{
|
|
auto const &points{Line.points};
|
|
if (points.empty())
|
|
return 0.0;
|
|
if (points.size() == 1 || Chainage <= points.front().chainage)
|
|
return points.front().elevation + (points.size() > 1 ? grade_after(Line, 0) * (Chainage - points.front().chainage) : 0.0);
|
|
if (Chainage >= points.back().chainage)
|
|
return points.back().elevation + grade_after(Line, points.size() - 2) * (Chainage - points.back().chainage);
|
|
auto const next{std::upper_bound(points.begin(), points.end(), Chainage, [](double const Value, pvi const &Point) { return Value < Point.chainage; })};
|
|
auto const index{static_cast<std::size_t>(std::distance(points.begin(), next)) - 1};
|
|
return points[index].elevation + grade_after(Line, index) * (Chainage - points[index].chainage);
|
|
}
|
|
|
|
bool overlaps(double const From, double const To, context::zone const &Zone)
|
|
{
|
|
return From < Zone.to - 1e-6 && To > Zone.from + 1e-6;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
double grade_after(line const &Line, std::size_t const Index)
|
|
{
|
|
auto const &points{Line.points};
|
|
if (Index + 1 >= points.size())
|
|
return 0.0;
|
|
auto const run{points[Index + 1].chainage - points[Index].chainage};
|
|
return run > 1e-9 ? (points[Index + 1].elevation - points[Index].elevation) / run : 0.0;
|
|
}
|
|
|
|
double required_radius(line const &Line)
|
|
{
|
|
return std::max(Line.norms.radius_min, Line.norms.radius_factor * Line.speed * Line.speed);
|
|
}
|
|
|
|
void assign_radii(line &Line)
|
|
{
|
|
auto const required{std::ceil(required_radius(Line) / 10.0) * 10.0};
|
|
for (std::size_t k = 1; k + 1 < Line.points.size(); ++k)
|
|
{
|
|
auto &point{Line.points[k]};
|
|
if (false == point.automatic || point.joint)
|
|
continue;
|
|
auto const change{std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) * 1000.0};
|
|
point.radius = change > Line.norms.curve_threshold + 1e-6 ? required : 0.0;
|
|
}
|
|
}
|
|
|
|
void complete(line &Line, context const &Context)
|
|
{
|
|
auto &points{Line.points};
|
|
points.erase(std::remove_if(points.begin(), points.end(), [](pvi const &Point) { return Point.joint; }), points.end());
|
|
if (points.size() < 2)
|
|
return;
|
|
bool const start{Line.join_start && Context.start_joined};
|
|
bool const end{Line.join_end && Context.end_joined};
|
|
if (start)
|
|
points.front().elevation = Context.start_elevation;
|
|
if (end)
|
|
points.back().elevation = Context.end_elevation;
|
|
|
|
// the vertical curve of the joint starts at the end of the route, off the adjoining track, or past the switches and the shared
|
|
// timbers the end lies on: over them the grade of the adjoining track goes on unchanged
|
|
auto const joint = [&](pvi const &End, pvi const &Neighbour, double const Grade, double const Direction) -> std::optional<pvi> {
|
|
auto lead{0.0};
|
|
for (bool grown{true}; grown;)
|
|
{
|
|
grown = false;
|
|
for (auto const &zone : Context.switches)
|
|
{
|
|
auto const to_start{Direction > 0.0 ? zone.from - End.chainage : End.chainage - zone.to};
|
|
auto const to_end{Direction > 0.0 ? zone.to - End.chainage : End.chainage - zone.from};
|
|
if (to_start <= lead + 0.5 && to_end > lead)
|
|
{
|
|
lead = to_end;
|
|
grown = true;
|
|
}
|
|
}
|
|
}
|
|
auto const run{(Neighbour.chainage - End.chainage) * Direction - lead};
|
|
if (run < 1.0)
|
|
return std::nullopt;
|
|
auto const offset{std::abs(Neighbour.elevation - (End.elevation + Grade * (lead + run) * Direction))};
|
|
if (offset / run * 1000.0 <= Line.norms.curve_threshold + 1e-6)
|
|
{
|
|
// no curve needed, but the grade breaks only past the timbers the end lies on
|
|
if (lead < 0.5)
|
|
return std::nullopt;
|
|
pvi result;
|
|
result.chainage = End.chainage + lead * Direction;
|
|
result.elevation = End.elevation + Grade * lead * Direction;
|
|
result.automatic = false;
|
|
result.joint = true;
|
|
return result;
|
|
}
|
|
auto radius{std::ceil(required_radius(Line) / 10.0) * 10.0};
|
|
double tangent;
|
|
if (run * run >= 2.0 * radius * offset)
|
|
tangent = (run - std::sqrt(run * run - 2.0 * radius * offset)) * 0.5;
|
|
else
|
|
{
|
|
tangent = run * 0.5;
|
|
radius = run * run / (2.0 * offset);
|
|
}
|
|
pvi result;
|
|
result.chainage = End.chainage + (lead + tangent) * Direction;
|
|
result.elevation = End.elevation + Grade * (lead + tangent) * Direction;
|
|
result.radius = radius;
|
|
result.automatic = false;
|
|
result.joint = true;
|
|
return result;
|
|
};
|
|
if (start)
|
|
if (auto const point{joint(points.front(), points[1], Context.start_grade, 1.0)})
|
|
points.insert(points.begin() + 1, *point);
|
|
if (end)
|
|
if (auto const point{joint(points.back(), points[points.size() - 2], Context.end_grade, -1.0)})
|
|
points.insert(points.end() - 1, *point);
|
|
assign_radii(Line);
|
|
}
|
|
|
|
std::vector<curve> curves(line const &Line)
|
|
{
|
|
std::vector<curve> result;
|
|
auto const &points{Line.points};
|
|
for (std::size_t k = 1; k + 1 < points.size(); ++k)
|
|
{
|
|
auto const radius{points[k].radius};
|
|
if (radius <= 0.0)
|
|
continue;
|
|
curve c;
|
|
c.point = k;
|
|
c.change = grade_after(Line, k) - grade_after(Line, k - 1);
|
|
c.tangent = radius * std::abs(c.change) * 0.5;
|
|
if (c.tangent < 1e-6)
|
|
continue;
|
|
c.sagitta = c.tangent * c.tangent / (2.0 * radius);
|
|
c.start = points[k].chainage - c.tangent;
|
|
c.end = points[k].chainage + c.tangent;
|
|
result.push_back(c);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
double elevation(line const &Line, double const Chainage)
|
|
{
|
|
auto height{polyline(Line, Chainage)};
|
|
for (auto const &c : curves(Line))
|
|
{
|
|
auto const distance{std::abs(Chainage - Line.points[c.point].chainage)};
|
|
if (distance >= c.tangent)
|
|
continue;
|
|
height += c.change * (c.tangent - distance) * (c.tangent - distance) / (4.0 * c.tangent);
|
|
}
|
|
return height;
|
|
}
|
|
|
|
double grade(line const &Line, double const Chainage)
|
|
{
|
|
auto const &points{Line.points};
|
|
if (points.size() < 2)
|
|
return 0.0;
|
|
std::size_t index{0};
|
|
while (index + 2 < points.size() && Chainage >= points[index + 1].chainage)
|
|
++index;
|
|
auto slope{grade_after(Line, index)};
|
|
for (auto const &c : curves(Line))
|
|
{
|
|
auto const offset{Chainage - points[c.point].chainage};
|
|
if (std::abs(offset) >= c.tangent)
|
|
continue;
|
|
slope -= c.change * (c.tangent - std::abs(offset)) / (2.0 * c.tangent) * (offset < 0.0 ? -1.0 : 1.0);
|
|
}
|
|
return slope;
|
|
}
|
|
|
|
double radius_for_sagitta(line const &Line, std::size_t const Point, double const Sagitta)
|
|
{
|
|
if (Point == 0 || Point + 1 >= Line.points.size())
|
|
return 0.0;
|
|
auto const change{std::abs(grade_after(Line, Point) - grade_after(Line, Point - 1))};
|
|
if (change < 1e-9)
|
|
return Line.points[Point].radius;
|
|
return std::max(0.0, 8.0 * std::abs(Sagitta) / (change * change));
|
|
}
|
|
|
|
std::vector<issue> check(line const &Line, context const &Context)
|
|
{
|
|
std::vector<issue> result;
|
|
auto const &points{Line.points};
|
|
auto const &norms{Line.norms};
|
|
auto const add = [&](bool const Error, int const Point, double const Chainage, std::string Text) {
|
|
result.push_back({Error, Point, Chainage, std::move(Text)});
|
|
};
|
|
if (points.size() < 2)
|
|
{
|
|
add(true, -1, 0.0, STR_C("The grade line needs at least the two ends"));
|
|
return result;
|
|
}
|
|
for (std::size_t k = 0; k + 1 < points.size(); ++k)
|
|
{
|
|
if (points[k + 1].chainage - points[k].chainage < 0.5)
|
|
add(true, static_cast<int>(k + 1), points[k + 1].chainage, format(STR_C("Points %zu and %zu are closer than 0.5 m"), k + 1, k + 2));
|
|
auto const slope{grade_after(Line, k) * 1000.0};
|
|
if (std::abs(slope) > norms.grade_max + 1e-6)
|
|
add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Grade %.2f per mille between points %zu and %zu exceeds %.2f"), slope, k + 1, k + 2, norms.grade_max));
|
|
}
|
|
|
|
auto const list{curves(Line)};
|
|
auto const required{required_radius(Line)};
|
|
for (std::size_t k = 1; k + 1 < points.size(); ++k)
|
|
{
|
|
auto const change{std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) * 1000.0};
|
|
auto const radius{points[k].radius};
|
|
if (radius <= 0.0 && change > norms.curve_threshold + 1e-6)
|
|
add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: change of grade %.2f per mille needs a vertical curve"), k + 1, change));
|
|
else if (radius > 0.0 && radius < required - 1e-6)
|
|
add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: R %.0f m is below the required %.0f m for %.0f km/h"), k + 1, radius, required, Line.speed));
|
|
if (radius <= 0.0 && change > 1e-3)
|
|
for (auto const &zone : Context.switches)
|
|
// a break on the bound of a switch is let be, within what the chainages of routes laid over each other differ by
|
|
if (points[k].chainage > zone.from + 0.05 && points[k].chainage < zone.to - 0.05)
|
|
add(true, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: the grade changes inside switch %s"), k + 1, zone.name.c_str()));
|
|
}
|
|
for (std::size_t i = 0; i < list.size(); ++i)
|
|
{
|
|
auto const &c{list[i]};
|
|
auto const point{static_cast<int>(c.point)};
|
|
if (c.start < points.front().chainage - 1e-6 || c.end > points.back().chainage + 1e-6)
|
|
add(true, point, points[c.point].chainage, format(STR_C("Point %zu: the vertical curve reaches past the end of the route"), c.point + 1));
|
|
if (i + 1 < list.size() && c.end > list[i + 1].start + 1e-6)
|
|
add(true, point, c.end, format(STR_C("Vertical curves of points %zu and %zu overlap by %.2f m"), c.point + 1, list[i + 1].point + 1, c.end - list[i + 1].start));
|
|
for (auto const &zone : Context.switches)
|
|
if (overlaps(c.start, c.end, zone))
|
|
add(true, point, points[c.point].chainage, format(STR_C("Point %zu: the vertical curve lies on switch %s"), c.point + 1, zone.name.c_str()));
|
|
}
|
|
if (norms.element_min > 0.0)
|
|
{
|
|
for (std::size_t k = 1; k + 2 < points.size(); ++k)
|
|
{
|
|
auto const tangent = [&](std::size_t const Point) {
|
|
for (auto const &c : list)
|
|
if (c.point == Point)
|
|
return c.tangent;
|
|
return 0.0;
|
|
};
|
|
auto const free{(points[k + 1].chainage - tangent(k + 1)) - (points[k].chainage + tangent(k))};
|
|
if (free < norms.element_min - 1e-6)
|
|
add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Constant grade between points %zu and %zu is %.1f m long, less than %.0f m"), k + 1, k + 2, std::max(0.0, free), norms.element_min));
|
|
}
|
|
}
|
|
for (auto const &zone : Context.switches)
|
|
{
|
|
auto const slope{std::abs(grade(Line, (zone.from + zone.to) * 0.5)) * 1000.0};
|
|
if (norms.grade_switch_max > 0.0 && slope > norms.grade_switch_max + 1e-6)
|
|
add(false, -1, zone.from, format(STR_C("Grade %.2f per mille over switch %s exceeds %.2f"), slope, zone.name.c_str(), norms.grade_switch_max));
|
|
}
|
|
for (auto const &fixed : Context.fixed)
|
|
{
|
|
auto const departure{elevation(Line, fixed.first) - fixed.second};
|
|
if (std::abs(departure) > norms.fixed_tolerance)
|
|
{
|
|
add(true, -1, fixed.first, format(STR_C("At %s the grade line departs %.3f m from a fixed point"), format_chainage(fixed.first).c_str(), departure));
|
|
break;
|
|
}
|
|
}
|
|
auto const ends = [&](bool const Atend) {
|
|
auto const joined{Atend ? Context.end_joined : Context.start_joined};
|
|
if (false == joined)
|
|
return;
|
|
auto const chainage{Atend ? points.back().chainage : points.front().chainage};
|
|
auto const departure{elevation(Line, chainage) - (Atend ? Context.end_elevation : Context.start_elevation)};
|
|
if (std::abs(departure) > norms.fixed_tolerance)
|
|
add(false, Atend ? static_cast<int>(points.size()) - 1 : 0, chainage, format(Atend ? STR_C("The end of the grade line is %.3f m off the adjoining track") : STR_C("The start of the grade line is %.3f m off the adjoining track"), departure));
|
|
auto const slope{grade(Line, chainage)};
|
|
auto const change{std::abs(slope - (Atend ? Context.end_grade : Context.start_grade)) * 1000.0};
|
|
if (change > norms.curve_threshold + 1e-6)
|
|
add(false, Atend ? static_cast<int>(points.size()) - 1 : 0, chainage, format(Atend ? STR_C("Change of grade %.2f per mille against the adjoining track at the end") : STR_C("Change of grade %.2f per mille against the adjoining track at the start"), change));
|
|
};
|
|
ends(false);
|
|
ends(true);
|
|
std::stable_sort(result.begin(), result.end(), [](issue const &A, issue const &B) { return A.error && false == B.error; });
|
|
return result;
|
|
}
|
|
|
|
line recognize(std::vector<double> const &Chainages, std::vector<double> const &Elevations, std::vector<double> const &Grades, recognition const &Options)
|
|
{
|
|
line result;
|
|
auto const count{std::min({Chainages.size(), Elevations.size(), Grades.size()})};
|
|
if (count < 2)
|
|
return result;
|
|
|
|
struct part
|
|
{
|
|
std::size_t first, last;
|
|
double slope{0.0}, offset{0.0};
|
|
};
|
|
auto const fit = [&](part &Part) {
|
|
double sx{0.0}, sy{0.0}, sxx{0.0}, sxy{0.0};
|
|
auto const origin{Chainages[Part.first]};
|
|
auto const n{static_cast<double>(Part.last - Part.first + 1)};
|
|
for (auto i = Part.first; i <= Part.last; ++i)
|
|
{
|
|
auto const x{Chainages[i] - origin};
|
|
sx += x;
|
|
sy += Elevations[i];
|
|
sxx += x * x;
|
|
sxy += x * Elevations[i];
|
|
}
|
|
auto const denominator{n * sxx - sx * sx};
|
|
Part.slope = std::abs(denominator) > 1e-12 ? (n * sxy - sx * sy) / denominator : 0.0;
|
|
Part.offset = (sy - Part.slope * sx) / n - Part.slope * origin;
|
|
};
|
|
|
|
std::vector<part> parts;
|
|
for (std::size_t i = 0; i < count;)
|
|
{
|
|
auto j{i};
|
|
auto low{Grades[i]}, high{Grades[i]};
|
|
while (j + 1 < count && std::max(high, Grades[j + 1]) - std::min(low, Grades[j + 1]) <= Options.grade_tolerance)
|
|
{
|
|
++j;
|
|
low = std::min(low, Grades[j]);
|
|
high = std::max(high, Grades[j]);
|
|
}
|
|
if (Chainages[j] - Chainages[i] >= Options.minimum_length)
|
|
{
|
|
part p{i, j};
|
|
fit(p);
|
|
parts.push_back(p);
|
|
i = j + 1;
|
|
}
|
|
else
|
|
++i;
|
|
}
|
|
auto const apart = [&](part const &A, part const &B, double const Chainage) { return std::abs((A.slope - B.slope) * Chainage + A.offset - B.offset); };
|
|
for (std::size_t i = 0; i + 1 < parts.size();)
|
|
{
|
|
if (std::abs(parts[i].slope - parts[i + 1].slope) < Options.grade_tolerance * 0.5 && apart(parts[i], parts[i + 1], (Chainages[parts[i].last] + Chainages[parts[i + 1].first]) * 0.5) <= Options.step)
|
|
{
|
|
parts[i].last = parts[i + 1].last;
|
|
fit(parts[i]);
|
|
parts.erase(parts.begin() + i + 1);
|
|
}
|
|
else
|
|
++i;
|
|
}
|
|
|
|
result.points.push_back({Chainages.front(), Elevations.front(), 0.0});
|
|
for (std::size_t i = 0; i + 1 < parts.size(); ++i)
|
|
{
|
|
auto const &a{parts[i]};
|
|
auto const &b{parts[i + 1]};
|
|
auto const from{Chainages[a.last]};
|
|
auto const to{Chainages[b.first]};
|
|
auto const change{b.slope - a.slope};
|
|
auto chainage{std::abs(change) > 1e-9 ? (a.offset - b.offset) / change : (from + to) * 0.5};
|
|
chainage = std::clamp(chainage, from, std::max(from, to));
|
|
// the grades don't meet between the parts, as at a stretch lifted with the grade it had: a break where each ends
|
|
if (apart(a, b, chainage) > Options.step)
|
|
{
|
|
for (auto const &[at, slope, offset] : {std::tuple{from, a.slope, a.offset}, std::tuple{to, b.slope, b.offset}})
|
|
if (at - result.points.back().chainage >= 0.5)
|
|
{
|
|
pvi point;
|
|
point.automatic = false;
|
|
point.chainage = at;
|
|
point.elevation = slope * at + offset;
|
|
result.points.push_back(point);
|
|
}
|
|
continue;
|
|
}
|
|
pvi point;
|
|
point.automatic = false;
|
|
point.chainage = chainage;
|
|
point.elevation = a.slope * chainage + a.offset;
|
|
auto const gap{to - from};
|
|
if (gap >= Options.sharp_gap && std::abs(change) > 1e-9)
|
|
point.radius = std::round(gap / std::abs(change));
|
|
if (point.chainage - result.points.back().chainage >= 0.5)
|
|
result.points.push_back(point);
|
|
}
|
|
if (Chainages[count - 1] - result.points.back().chainage < 0.5 && result.points.size() > 1)
|
|
result.points.pop_back();
|
|
result.points.push_back({Chainages[count - 1], Elevations[count - 1], 0.0});
|
|
|
|
// the gap between constant grades is shorter than the curve, the sagitta gives the radius
|
|
for (std::size_t i = 1; i + 1 < result.points.size(); ++i)
|
|
{
|
|
auto &point{result.points[i]};
|
|
if (point.radius <= 0.0)
|
|
continue;
|
|
auto const next{std::lower_bound(Chainages.begin(), Chainages.begin() + count, point.chainage)};
|
|
if (next == Chainages.begin() || next == Chainages.begin() + count)
|
|
continue;
|
|
auto const index{static_cast<std::size_t>(std::distance(Chainages.begin(), next))};
|
|
auto const run{Chainages[index] - Chainages[index - 1]};
|
|
auto const measured{run > 1e-9 ? Elevations[index - 1] + (Elevations[index] - Elevations[index - 1]) * (point.chainage - Chainages[index - 1]) / run : Elevations[index]};
|
|
auto const radius{radius_for_sagitta(result, i, measured - point.elevation)};
|
|
if (radius > 0.0)
|
|
point.radius = std::round(radius);
|
|
}
|
|
if (Options.detail <= 0.0)
|
|
return result;
|
|
std::vector<char> off(count, 0);
|
|
for (std::size_t i = 0; i < count; ++i)
|
|
off[i] = std::abs(elevation(result, Chainages[i]) - Elevations[i]) > Options.detail ? 1 : 0;
|
|
std::vector<pvi> bends;
|
|
for (std::size_t i = 0; i < count;)
|
|
{
|
|
if (off[i] == 0)
|
|
{
|
|
++i;
|
|
continue;
|
|
}
|
|
auto const first{i > 0 ? i - 1 : i};
|
|
while (i < count && off[i] != 0)
|
|
++i;
|
|
auto const last{std::min(i, count - 1)};
|
|
std::vector<char> keep(count, 0);
|
|
keep[first] = keep[last] = 1;
|
|
std::vector<std::pair<std::size_t, std::size_t>> stack{{first, last}};
|
|
while (false == stack.empty())
|
|
{
|
|
auto const [a, b]{stack.back()};
|
|
stack.pop_back();
|
|
if (b <= a + 1)
|
|
continue;
|
|
auto const run{std::max(1e-9, Chainages[b] - Chainages[a])};
|
|
double worst{0.0};
|
|
std::size_t index{a};
|
|
for (auto k = a + 1; k < b; ++k)
|
|
{
|
|
auto const departure{std::abs(Elevations[k] - (Elevations[a] + (Elevations[b] - Elevations[a]) * (Chainages[k] - Chainages[a]) / run))};
|
|
if (departure > worst)
|
|
{
|
|
worst = departure;
|
|
index = k;
|
|
}
|
|
}
|
|
if (worst > Options.detail * 0.5)
|
|
{
|
|
keep[index] = 1;
|
|
stack.push_back({a, index});
|
|
stack.push_back({index, b});
|
|
}
|
|
}
|
|
auto const low{Chainages[first]}, high{Chainages[last]};
|
|
result.points.erase(std::remove_if(result.points.begin() + 1, result.points.end() - 1, [&](pvi const &Point) { return Point.chainage >= low - 0.5 && Point.chainage <= high + 0.5; }), result.points.end() - 1);
|
|
for (auto k = first; k <= last; ++k)
|
|
if (keep[k] != 0 && k > 0 && k + 1 < count)
|
|
{
|
|
pvi point;
|
|
point.automatic = false;
|
|
point.chainage = Chainages[k];
|
|
point.elevation = Elevations[k];
|
|
bends.push_back(point);
|
|
}
|
|
}
|
|
if (bends.empty())
|
|
return result;
|
|
result.points.insert(result.points.end() - 1, bends.begin(), bends.end());
|
|
std::sort(result.points.begin(), result.points.end(), [](pvi const &Left, pvi const &Right) { return Left.chainage < Right.chainage; });
|
|
std::vector<pvi> points;
|
|
for (auto const &point : result.points)
|
|
if (points.empty() || point.chainage - points.back().chainage >= 0.5 || &point == &result.points.back())
|
|
{
|
|
if (false == points.empty() && point.chainage - points.back().chainage < 0.5)
|
|
points.pop_back();
|
|
points.push_back(point);
|
|
}
|
|
result.points = std::move(points);
|
|
return result;
|
|
}
|
|
|
|
line fit_ground(std::vector<double> const &Chainages, std::vector<double> const &Ground, context const &Context, line const &Settings, ground_fit const &Options, std::vector<double> const &Tolerances,
|
|
std::vector<double> const &Reaches)
|
|
{
|
|
line result;
|
|
result.speed = Settings.speed;
|
|
result.norms = Settings.norms;
|
|
result.join_start = Settings.join_start;
|
|
result.join_end = Settings.join_end;
|
|
auto const count{std::min(Chainages.size(), Ground.size())};
|
|
if (count < 2 || std::none_of(Ground.begin(), Ground.begin() + count, [](double const Height) { return std::isfinite(Height); }))
|
|
return result;
|
|
auto const origin{Chainages.front()};
|
|
auto const length{Chainages[count - 1] - origin};
|
|
if (length < 1.0)
|
|
return result;
|
|
|
|
// the line is sought on nodes a few metres apart, straight between them: the grade can change only at a node
|
|
auto const nodes{std::max<std::size_t>(3, static_cast<std::size_t>(std::ceil(length / 5.0)) + 1)};
|
|
auto const step{length / static_cast<double>(nodes - 1)};
|
|
std::vector<std::size_t> cell(count);
|
|
std::vector<double> part(count);
|
|
for (std::size_t i = 0; i < count; ++i)
|
|
{
|
|
auto const at{std::clamp((Chainages[i] - origin) / step, 0.0, static_cast<double>(nodes - 1))};
|
|
cell[i] = std::min(nodes - 2, static_cast<std::size_t>(at));
|
|
part[i] = at - static_cast<double>(cell[i]);
|
|
}
|
|
auto const node = [&](std::size_t const K) { return origin + step * static_cast<double>(K); };
|
|
// a change of grade at a node is the curve of the required radius around it, which keeps off the switches: none on a
|
|
// switch, and beside one no more than a curve reaching up to it
|
|
auto const required{std::ceil(required_radius(result) / 10.0) * 10.0};
|
|
auto const sharpest{step / (required * 1.001)};
|
|
std::vector<double> bend_max(nodes, sharpest);
|
|
for (std::size_t k = 1; k + 1 < nodes; ++k)
|
|
for (auto const &zone : Context.switches)
|
|
bend_max[k] = std::min(bend_max[k], std::max(0.0, 2.0 * std::max(zone.from - node(k), node(k) - zone.to) / (required * 1.001)));
|
|
struct held
|
|
{
|
|
std::size_t cell;
|
|
double part;
|
|
double low;
|
|
double high;
|
|
};
|
|
std::vector<held> bounds;
|
|
for (auto const &bound : Options.bounds)
|
|
{
|
|
auto const at{std::clamp((bound.chainage - origin) / step, 0.0, static_cast<double>(nodes - 1))};
|
|
auto const cell{std::min(nodes - 2, static_cast<std::size_t>(at))};
|
|
bounds.push_back({cell, at - static_cast<double>(cell), std::min(bound.low, bound.high), std::max(bound.low, bound.high)});
|
|
}
|
|
bool const start_fixed{Settings.join_start && Context.start_joined};
|
|
bool const end_fixed{Settings.join_end && Context.end_joined};
|
|
|
|
// a change of grade costs the price times its size, so the line breaks where the ground keeps away from it long enough to
|
|
// be worth it. The ground costs the square of its distance over its tolerance; past its reach ever less, up to a bound:
|
|
// what parts from the line, a bridge, a viaduct, a wall, a track on another level, is gone past rather than followed. The
|
|
// bends are no sharper than the radius of the norms and the grades no steeper than the limit
|
|
auto const tolerance{std::max(0.01, Options.tolerance)};
|
|
auto const own = [&](std::size_t const I) { return I < Tolerances.size() && Tolerances[I] > 0.0 ? std::max(0.005, Tolerances[I]) : tolerance; };
|
|
auto const reach = [&](std::size_t const I) { return I < Reaches.size() && Reaches[I] > 0.0 ? std::max(0.005, Reaches[I]) : own(I); };
|
|
auto const price{24000.0 * tolerance * tolerance};
|
|
auto const grade_max{(Options.grade_max > 0.0 ? Options.grade_max : result.norms.grade_max) / 1000.0};
|
|
// what has to hold, the heights of the joined ends and the switches in one plane, holds before the grades of the ends, and
|
|
// these before the limits of the norms
|
|
double constexpr stiff{1e11};
|
|
double constexpr firm{1e9};
|
|
|
|
std::vector<double> x(nodes);
|
|
{
|
|
std::size_t i{0};
|
|
for (std::size_t k = 0; k < nodes; ++k)
|
|
{
|
|
while (i + 1 < count && Chainages[i] < node(k))
|
|
++i;
|
|
auto j{i};
|
|
while (j > 0 && false == std::isfinite(Ground[j]))
|
|
--j;
|
|
while (j + 1 < count && false == std::isfinite(Ground[j]))
|
|
++j;
|
|
x[k] = Ground[j];
|
|
}
|
|
}
|
|
// the cost of a line, with its slope and the band of its curvature when asked. The distance of the ground counts up to
|
|
// its tolerance times the root of the spread, the size of a change of grade is rounded off within the flat
|
|
auto spread{256.0};
|
|
auto flat{sharpest};
|
|
std::vector<double> gradient(nodes), diagonal(nodes), first(nodes), second(nodes);
|
|
auto const cost = [&](std::vector<double> const &Line, bool const Derivatives) {
|
|
if (Derivatives)
|
|
{
|
|
std::fill(gradient.begin(), gradient.end(), 0.0);
|
|
std::fill(diagonal.begin(), diagonal.end(), 1e-9);
|
|
std::fill(first.begin(), first.end(), 0.0);
|
|
std::fill(second.begin(), second.end(), 0.0);
|
|
}
|
|
auto total{0.0};
|
|
auto const add = [&](std::initializer_list<std::pair<std::size_t, double>> Terms, double const Value, double const Slope, double const Curvature) {
|
|
total += Value;
|
|
if (false == Derivatives)
|
|
return;
|
|
for (auto const &[a, ca] : Terms)
|
|
{
|
|
gradient[a] += Slope * ca;
|
|
for (auto const &[b, cb] : Terms)
|
|
{
|
|
if (b == a)
|
|
diagonal[a] += Curvature * ca * cb;
|
|
else if (b == a + 1)
|
|
first[a] += Curvature * ca * cb;
|
|
else if (b == a + 2)
|
|
second[a] += Curvature * ca * cb;
|
|
}
|
|
}
|
|
};
|
|
auto const square = [&](std::initializer_list<std::pair<std::size_t, double>> Terms, double const Off, double const Weight) { add(Terms, Weight * Off * Off, 2.0 * Weight * Off, 2.0 * Weight); };
|
|
for (std::size_t i = 0; i < count; ++i)
|
|
{
|
|
if (false == std::isfinite(Ground[i]))
|
|
continue;
|
|
auto const off{Line[cell[i]] * (1.0 - part[i]) + Line[cell[i] + 1] * part[i] - Ground[i]};
|
|
auto const agreeing{spread * reach(i) * reach(i)};
|
|
auto const share{agreeing / (agreeing + off * off)};
|
|
auto const weight{tolerance * tolerance / (own(i) * own(i))};
|
|
add({{cell[i], 1.0 - part[i]}, {cell[i] + 1, part[i]}}, weight * off * off * share, weight * 2.0 * off * share * share, weight * 2.0 * share * share);
|
|
}
|
|
for (std::size_t k = 1; k + 1 < nodes; ++k)
|
|
{
|
|
auto const terms = {std::pair{k - 1, 1.0 / step}, std::pair{k, -2.0 / step}, std::pair{k + 1, 1.0 / step}};
|
|
auto const bend{(Line[k - 1] - 2.0 * Line[k] + Line[k + 1]) / step};
|
|
if (bend_max[k] <= 0.0)
|
|
{
|
|
square(terms, bend, stiff);
|
|
continue;
|
|
}
|
|
auto const root{std::sqrt(bend * bend + flat * flat)};
|
|
add(terms, price * root, price * bend / root, price * flat * flat / (root * root * root));
|
|
if (std::abs(bend) > bend_max[k])
|
|
square(terms, bend - std::copysign(bend_max[k], bend), firm);
|
|
}
|
|
for (std::size_t k = 0; k + 1 < nodes; ++k)
|
|
{
|
|
auto const grade{(Line[k + 1] - Line[k]) / step};
|
|
if (std::abs(grade) > grade_max)
|
|
square({{k, -1.0 / step}, {k + 1, 1.0 / step}}, grade - std::copysign(grade_max, grade), firm);
|
|
}
|
|
if (start_fixed)
|
|
{
|
|
square({{0, 1.0}}, Line[0] - Context.start_elevation, stiff);
|
|
square({{0, -1.0 / step}, {1, 1.0 / step}}, (Line[1] - Line[0]) / step - Context.start_grade, firm * 10.0);
|
|
}
|
|
if (end_fixed)
|
|
{
|
|
square({{nodes - 1, 1.0}}, Line[nodes - 1] - Context.end_elevation, stiff);
|
|
square({{nodes - 2, -1.0 / step}, {nodes - 1, 1.0 / step}}, (Line[nodes - 1] - Line[nodes - 2]) / step - Context.end_grade, firm * 10.0);
|
|
}
|
|
for (auto const &bound : bounds)
|
|
{
|
|
auto const height{Line[bound.cell] * (1.0 - bound.part) + Line[bound.cell + 1] * bound.part};
|
|
if (height < bound.low || height > bound.high)
|
|
square({{bound.cell, 1.0 - bound.part}, {bound.cell + 1, bound.part}}, height - std::clamp(height, bound.low, bound.high), firm * 10.0);
|
|
}
|
|
return total;
|
|
};
|
|
// the step of Newton, the band factored as L D L^T with two diagonals below the main one
|
|
auto const newton = [&]() {
|
|
std::vector<double> d(nodes), l1(nodes, 0.0), l2(nodes, 0.0), z(nodes);
|
|
for (std::size_t i = 0; i < nodes; ++i)
|
|
{
|
|
d[i] = diagonal[i] - (i >= 1 ? l1[i - 1] * l1[i - 1] * d[i - 1] : 0.0) - (i >= 2 ? l2[i - 2] * l2[i - 2] * d[i - 2] : 0.0);
|
|
if (i + 1 < nodes)
|
|
l1[i] = (first[i] - (i >= 1 ? l2[i - 1] * l1[i - 1] * d[i - 1] : 0.0)) / d[i];
|
|
if (i + 2 < nodes)
|
|
l2[i] = second[i] / d[i];
|
|
}
|
|
for (std::size_t i = 0; i < nodes; ++i)
|
|
z[i] = -gradient[i] - (i >= 1 ? l1[i - 1] * z[i - 1] : 0.0) - (i >= 2 ? l2[i - 2] * z[i - 2] : 0.0);
|
|
for (std::size_t i = 0; i < nodes; ++i)
|
|
z[i] /= d[i];
|
|
for (auto i = nodes - 1; i-- > 0;)
|
|
z[i] -= l1[i] * z[i + 1] + (i + 2 < nodes ? l2[i] * z[i + 2] : 0.0);
|
|
return z;
|
|
};
|
|
// the spread and the flat narrow down level by level, each level going on from the line of the one before: from nearly
|
|
// all the ground agreeing and changes of grade costing their square, to the tolerance and their size
|
|
for (;;)
|
|
{
|
|
auto current{cost(x, true)};
|
|
for (int pass = 0; pass < 100; ++pass)
|
|
{
|
|
auto const direction{newton()};
|
|
auto descent{0.0};
|
|
for (std::size_t k = 0; k < nodes; ++k)
|
|
descent += gradient[k] * direction[k];
|
|
std::vector<double> trial(nodes);
|
|
auto length{1.0};
|
|
for (;; length *= 0.5)
|
|
{
|
|
for (std::size_t k = 0; k < nodes; ++k)
|
|
trial[k] = x[k] + length * direction[k];
|
|
if (cost(trial, false) <= current + 1e-4 * length * descent || length < 1e-6)
|
|
break;
|
|
}
|
|
auto moved{0.0};
|
|
for (std::size_t k = 0; k < nodes; ++k)
|
|
moved = std::max(moved, std::abs(length * direction[k]));
|
|
x = std::move(trial);
|
|
current = cost(x, true);
|
|
if (moved < 1e-5)
|
|
break;
|
|
}
|
|
if (spread <= 1.0 && flat <= 1e-6)
|
|
break;
|
|
spread = std::max(1.0, spread * 0.5);
|
|
flat = std::max(1e-6, flat * 0.5);
|
|
}
|
|
|
|
// the design: the changes of grade along a bend, one after another the same way, are one point of intersection where
|
|
// the grades before and after them meet. No change at a node is sharper than the required radius, so the curve of that
|
|
// radius stays within the nodes of its bend; where the ends held the line sharper, the curve is as sharp as the line
|
|
auto const change = [&](std::size_t const K) { return (x[K - 1] - 2.0 * x[K] + x[K + 1]) / step; };
|
|
auto const significant{sharpest * 0.01};
|
|
auto const between = [&](std::size_t const K, std::size_t const L) {
|
|
return std::any_of(Context.switches.begin(), Context.switches.end(), [&](auto const &Zone) { return Zone.to > node(K) && Zone.from < node(L); });
|
|
};
|
|
auto const bends = [&](std::size_t const K) { return bend_max[K] > 0.0 && std::abs(change(K)) >= significant; };
|
|
// an end whose grade the line couldn't take over keeps only its height, the grade breaks at the joint
|
|
if (start_fixed && std::abs((x[1] - x[0]) / step - Context.start_grade) > significant)
|
|
result.join_start = false;
|
|
if (end_fixed && std::abs((x[nodes - 1] - x[nodes - 2]) / step - Context.end_grade) > significant)
|
|
result.join_end = false;
|
|
result.points.push_back({origin, x[0]});
|
|
std::vector<double> rooms{0.0};
|
|
for (std::size_t k = 1; k + 1 < nodes;)
|
|
{
|
|
if (false == bends(k))
|
|
{
|
|
++k;
|
|
continue;
|
|
}
|
|
auto const first{k};
|
|
auto const sign{change(k) > 0.0};
|
|
double sum{0.0}, moment{0.0};
|
|
for (; k + 1 < nodes && bends(k) && (change(k) > 0.0) == sign && (k == first || false == between(k - 1, k)); ++k)
|
|
{
|
|
sum += change(k);
|
|
moment += change(k) * node(k);
|
|
}
|
|
auto const at{moment / sum};
|
|
pvi point{at, x[first] + (x[first] - x[first - 1]) / step * (at - node(first))};
|
|
auto room{std::min(at - node(first) + step * 0.5, node(k - 1) + step * 0.5 - at)};
|
|
for (auto const &zone : Context.switches)
|
|
room = std::min(room, std::max(0.0, std::max(zone.from - at, at - zone.to)));
|
|
result.points.push_back(point);
|
|
rooms.push_back(room);
|
|
}
|
|
result.points.push_back({node(nodes - 1), x[nodes - 1]});
|
|
rooms.push_back(0.0);
|
|
// the curve of the required radius where it has room, off the switches and the curves beside; where it hasn't, as sharp as
|
|
// the room lets it be
|
|
for (std::size_t i = 1; i + 1 < result.points.size(); ++i)
|
|
{
|
|
auto const turn{std::abs(grade_after(result, i) - grade_after(result, i - 1))};
|
|
if (turn * 1000.0 > result.norms.curve_threshold + 1e-6 && required * turn * 0.5 > rooms[i] - 1e-4)
|
|
{
|
|
result.points[i].radius = 2.0 * std::max(0.0, rooms[i] - 1e-4) / turn;
|
|
result.points[i].automatic = false;
|
|
}
|
|
}
|
|
assign_radii(result);
|
|
|
|
// then the simplest design which keeps within half the tolerance of the line found, the one of the ground agreeing with
|
|
// it: a point is dropped, or two bending the same way become one, as long as that departs the least and adds no fault
|
|
std::vector<double> limit(nodes, tolerance);
|
|
for (std::size_t i = 0; i < count; ++i)
|
|
if (std::isfinite(Ground[i]) && std::abs(Ground[i] - (x[cell[i]] * (1.0 - part[i]) + x[cell[i] + 1] * part[i])) <= reach(i))
|
|
for (auto const k : {cell[i], cell[i] + 1})
|
|
limit[k] = std::min(limit[k], own(i));
|
|
for (auto const &bound : bounds)
|
|
for (auto const k : {bound.cell, bound.cell + 1})
|
|
limit[k] = std::min(limit[k], std::max(0.001, (bound.high - bound.low) * 0.5));
|
|
auto const errors = [&](line const &Line) {
|
|
auto completed{Line};
|
|
complete(completed, Context);
|
|
auto const found{check(completed, Context)};
|
|
return std::count_if(found.begin(), found.end(), [](issue const &Issue) { return Issue.error; });
|
|
};
|
|
auto const departure = [&](line const &Line, double const From, double const To) {
|
|
auto worst{0.0};
|
|
for (std::size_t k = 0; k < nodes; ++k)
|
|
if (node(k) >= From && node(k) <= To)
|
|
worst = std::max(worst, std::abs(elevation(Line, node(k)) - x[k]) / limit[k]);
|
|
return worst;
|
|
};
|
|
auto faults{errors(result)};
|
|
// an end whose grade the curve of the joint can't take over without a fault breaks at the joint
|
|
for (auto *join : {&result.join_start, &result.join_end})
|
|
if (faults > 0 && *join)
|
|
{
|
|
*join = false;
|
|
if (auto const fewer{errors(result)}; fewer < faults)
|
|
faults = fewer;
|
|
else
|
|
*join = true;
|
|
}
|
|
for (;;)
|
|
{
|
|
auto const &points{result.points};
|
|
auto const last{points.size() - 1};
|
|
std::vector<std::pair<double, line>> candidates;
|
|
for (std::size_t i = 1; i < last; ++i)
|
|
{
|
|
auto const from{points[i >= 2 ? i - 2 : 0].chainage};
|
|
auto dropped{result};
|
|
dropped.points.erase(dropped.points.begin() + i);
|
|
assign_radii(dropped);
|
|
candidates.emplace_back(departure(dropped, from, points[std::min(last, i + 2)].chainage), std::move(dropped));
|
|
if (i + 1 >= last)
|
|
continue;
|
|
auto const before{grade_after(result, i - 1)};
|
|
auto const between{grade_after(result, i)};
|
|
auto const after{grade_after(result, i + 1)};
|
|
if ((between - before) * (after - between) <= 0.0)
|
|
continue;
|
|
auto const at{(points[i + 1].elevation - points[i].elevation + before * points[i].chainage - after * points[i + 1].chainage) / (before - after)};
|
|
auto merged{result};
|
|
merged.points[i] = {at, points[i].elevation + before * (at - points[i].chainage)};
|
|
merged.points.erase(merged.points.begin() + i + 1);
|
|
assign_radii(merged);
|
|
candidates.emplace_back(departure(merged, from, points[std::min(last, i + 3)].chainage), std::move(merged));
|
|
}
|
|
std::sort(candidates.begin(), candidates.end(), [](auto const &A, auto const &B) { return A.first < B.first; });
|
|
auto const chosen{std::find_if(candidates.begin(), candidates.end(), [&](auto const &Candidate) { return Candidate.first <= 0.5 && errors(Candidate.second) <= faults; })};
|
|
if (chosen == candidates.end())
|
|
break;
|
|
result = std::move(chosen->second);
|
|
faults = errors(result);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
std::vector<departure> departures(line const &Line, std::vector<double> const &Chainages, std::vector<double> const &Ground, double const Tolerance)
|
|
{
|
|
std::vector<departure> result;
|
|
auto const count{std::min(Chainages.size(), Ground.size())};
|
|
auto const off = [&](std::size_t const I) { return std::isfinite(Ground[I]) ? Ground[I] - elevation(Line, Chainages[I]) : 0.0; };
|
|
for (std::size_t i = 0; i < count;)
|
|
{
|
|
if (std::abs(off(i)) <= Tolerance)
|
|
{
|
|
++i;
|
|
continue;
|
|
}
|
|
departure stretch;
|
|
stretch.from = Chainages[i];
|
|
auto const below{off(i) < 0.0};
|
|
for (; i < count && std::abs(off(i)) > Tolerance && (off(i) < 0.0) == below; ++i)
|
|
{
|
|
if (std::abs(off(i)) > std::abs(stretch.depth))
|
|
stretch.depth = off(i);
|
|
stretch.to = Chainages[i];
|
|
}
|
|
result.push_back(stretch);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
std::vector<double> breaks(line const &Line)
|
|
{
|
|
std::vector<double> result;
|
|
auto const &points{Line.points};
|
|
for (std::size_t k = 1; k + 1 < points.size(); ++k)
|
|
{
|
|
if (points[k].radius > 0.0)
|
|
continue;
|
|
if (std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) > 1e-7)
|
|
result.push_back(points[k].chainage);
|
|
}
|
|
for (auto const &c : curves(Line))
|
|
{
|
|
result.push_back(c.start);
|
|
result.push_back(c.end);
|
|
}
|
|
std::sort(result.begin(), result.end());
|
|
return result;
|
|
}
|
|
|
|
std::string format_chainage(double const Chainage)
|
|
{
|
|
auto const value{std::abs(Chainage)};
|
|
auto const kilometres{static_cast<long>(value / 1000.0)};
|
|
return format("%s%ld+%06.2f", Chainage < 0.0 ? "-" : "", kilometres, value - kilometres * 1000.0);
|
|
}
|
|
|
|
}
|