eu07/maszyna
17
0
Fork 0
mirror of https://github.com/MaSzyna-EU07/maszyna.git synced 2026-10-23 17:15:15 +02:00
Code Issues 152 Wiki Activity

feat(editor): labels for a layman, light at the start; whole sceneries fitted to the terrain model

Browse Source 
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
This commit is contained in:
maj00r
2026-10-19 15:27:50 +02:00
parent 5a00d88be9
commit ff1710d9ea
78 changed files with 15011 additions and 4386 deletions
Download Patch File Download Diff File Expand all files Collapse all files

681
editor/editorProfile.cpp
View File

@@ -10,11 +10,13 @@ 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
{
@@ -83,14 +85,40 @@ void complete(line &Line, context const &Context)
if (end)
points.back().elevation = Context.end_elevation;
// the vertical curve of the joint starts right at the end of the route, off the adjoining track
// 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 const run{(Neighbour.chainage - End.chainage) * Direction};
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 * run * Direction))};
auto const offset{std::abs(Neighbour.elevation - (End.elevation + Grade * (lead + run) * Direction))};
if (offset / run * 1000.0 <= Line.norms.curve_threshold + 1e-6)
return std::nullopt;
{
// 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)
@@ -101,8 +129,8 @@ void complete(line &Line, context const &Context)
radius = run * run / (2.0 * offset);
}
pvi result;
result.chainage = End.chainage + tangent * Direction;
result.elevation = End.elevation + Grade * tangent * Direction;
result.chainage = End.chainage + (lead + tangent) * Direction;
result.elevation = End.elevation + Grade * (lead + tangent) * Direction;
result.radius = radius;
result.automatic = false;
result.joint = true;
@@ -192,16 +220,16 @@ std::vector<issue> check(line const &Line, context const &Context)
};
if (points.size() < 2)
{
add(true, -1, 0.0, "The grade line needs at least the two ends");
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("Points %zu and %zu are closer than 0.5 m", k + 1, k + 2));
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("Grade %.2f per mille between points %zu and %zu exceeds %.2f", slope, k + 1, k + 2, norms.grade_max));
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)};
@@ -211,25 +239,26 @@ std::vector<issue> check(line const &Line, context const &Context)
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("Point %zu: change of grade %.2f per mille needs a vertical curve", k + 1, change));
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("Point %zu: R %.0f m is below the required %.0f m for %.0f km/h", k + 1, radius, required, Line.speed));
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)
if (points[k].chainage > zone.from + 1e-3 && points[k].chainage < zone.to - 1e-3)
add(true, static_cast<int>(k), points[k].chainage, format("Point %zu: the grade changes inside switch %s", k + 1, zone.name.c_str()));
// 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("Point %zu: the vertical curve reaches past the end of the route", c.point + 1));
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("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));
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("Point %zu: the vertical curve lies on switch %s", c.point + 1, zone.name.c_str()));
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)
{
@@ -243,21 +272,21 @@ std::vector<issue> check(line const &Line, context const &Context)
};
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("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));
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 (slope > norms.grade_switch_max + 1e-6)
add(false, -1, zone.from, format("Grade %.2f per mille over switch %s exceeds %.2f", slope, zone.name.c_str(), norms.grade_switch_max));
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("At %s the grade line departs %.3f m from a fixed point", format_chainage(fixed.first).c_str(), departure));
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;
}
}
@@ -268,11 +297,11 @@ std::vector<issue> check(line const &Line, context const &Context)
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("The %s of the grade line is %.3f m off the adjoining track", Atend ? "end" : "start", departure));
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("Change of grade %.2f per mille against the adjoining track at the %s", change, Atend ? "end" : "start"));
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);
@@ -330,9 +359,10 @@ line recognize(std::vector<double> const &Chainages, std::vector<double> const &
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)
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]);
@@ -352,6 +382,20 @@ line recognize(std::vector<double> const &Chainages, std::vector<double> const &
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;
@@ -382,72 +426,26 @@ line recognize(std::vector<double> const &Chainages, std::vector<double> const &
if (radius > 0.0)
point.radius = std::round(radius);
}
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)
{
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)
if (Options.detail <= 0.0)
return result;
// gaps in the ground bridged straight, past the outermost samples it stays level
std::vector<double> ground(Ground.begin(), Ground.begin() + count);
std::vector<std::size_t> known;
std::vector<char> off(count, 0);
for (std::size_t i = 0; i < count; ++i)
if (std::isfinite(ground[i]))
known.push_back(i);
if (known.empty())
return result;
for (std::size_t i = 0, k = 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 (std::isfinite(ground[i]))
if (off[i] == 0)
{
++i;
continue;
while (k + 1 < known.size() && known[k + 1] < i)
++k;
if (i < known.front())
ground[i] = Ground[known.front()];
else if (k + 1 >= known.size())
ground[i] = Ground[known.back()];
else
{
auto const a{known[k]}, b{known[k + 1]};
ground[i] = Ground[a] + (Ground[b] - Ground[a]) * (Chainages[i] - Chainages[a]) / std::max(1e-9, Chainages[b] - Chainages[a]);
}
}
if (Options.smoothing > 1.0)
{
std::vector<double> sums(count + 1, 0.0);
for (std::size_t i = 0; i < count; ++i)
sums[i + 1] = sums[i] + ground[i];
std::vector<double> smooth(count);
std::size_t low{0}, high{0};
for (std::size_t i = 0; i < count; ++i)
{
while (Chainages[i] - Chainages[low] > Options.smoothing * 0.5)
++low;
while (high + 1 < count && Chainages[high + 1] - Chainages[i] <= Options.smoothing * 0.5)
++high;
high = std::max(high, i);
smooth[i] = (sums[high + 1] - sums[low]) / static_cast<double>(high - low + 1);
}
ground.swap(smooth);
}
bool const start_fixed{Settings.join_start && Context.start_joined};
bool const end_fixed{Settings.join_end && Context.end_joined};
// first guess of the points: the ground simplified to within the tolerance
std::vector<std::size_t> knots{0, count - 1};
{
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.front() = keep.back() = 1;
std::vector<std::pair<std::size_t, std::size_t>> stack{{0, count - 1}};
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()};
@@ -457,179 +455,412 @@ line fit_ground(std::vector<double> const &Chainages, std::vector<double> const
auto const run{std::max(1e-9, Chainages[b] - Chainages[a])};
double worst{0.0};
std::size_t index{a};
for (auto i = a + 1; i < b; ++i)
for (auto k = a + 1; k < b; ++k)
{
auto const chord{ground[a] + (ground[b] - ground[a]) * (Chainages[i] - Chainages[a]) / run};
auto const departure{std::abs(ground[i] - chord)};
auto const departure{std::abs(Elevations[k] - (Elevations[a] + (Elevations[b] - Elevations[a]) * (Chainages[k] - Chainages[a]) / run))};
if (departure > worst)
{
worst = departure;
index = i;
index = k;
}
}
if (worst > Options.tolerance)
if (worst > Options.detail * 0.5)
{
keep[index] = 1;
stack.push_back({a, index});
stack.push_back({index, b});
}
}
knots.clear();
for (std::size_t i = 0; i < count; ++i)
if (keep[i] != 0)
knots.push_back(i);
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;
}
auto const segment_limit = [&](double const From, double const To) {
auto limit{result.norms.grade_max};
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)
if (overlaps(From, To, zone))
limit = std::min(limit, result.norms.grade_switch_max);
return limit / 1000.0;
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;
};
// elevations of the points with the least squares of the departures from the ground
auto const solve = [&](std::vector<std::size_t> const &Knots) {
auto const n{Knots.size()};
std::vector<double> lower(n, 0.0), diagonal(n, 1e-9), upper(n, 0.0), right(n, 0.0);
auto const add = [&](std::size_t const Interval, double const T, double const Value, double const Weight) {
diagonal[Interval] += Weight * (1.0 - T) * (1.0 - T);
diagonal[Interval + 1] += Weight * T * T;
upper[Interval] += Weight * (1.0 - T) * T;
lower[Interval + 1] += Weight * (1.0 - T) * T;
right[Interval] += Weight * (1.0 - T) * Value;
right[Interval + 1] += Weight * T * Value;
};
auto const interval_of = [&](double const Chainage) {
std::size_t k{0};
while (k + 2 < n && Chainage >= Chainages[Knots[k + 1]])
++k;
return k;
};
for (std::size_t k = 0; k + 1 < n; ++k)
{
auto const from{Chainages[Knots[k]]};
auto const run{std::max(1e-9, Chainages[Knots[k + 1]] - from)};
auto const last{k + 2 == n ? Knots[k + 1] : Knots[k + 1] - 1};
for (auto i = Knots[k]; i <= last; ++i)
add(k, (Chainages[i] - from) / run, ground[i], 1.0);
}
for (auto const &fixed : Context.fixed)
{
auto const k{interval_of(fixed.first)};
auto const from{Chainages[Knots[k]]};
auto const t{std::clamp((fixed.first - from) / std::max(1e-9, Chainages[Knots[k + 1]] - from), 0.0, 1.0)};
add(k, t, fixed.second, 1e6);
}
auto const pin = [&](std::size_t const Index, double const Value) {
lower[Index] = upper[Index] = 0.0;
diagonal[Index] = 1.0;
right[Index] = Value;
};
if (start_fixed)
pin(0, Context.start_elevation);
if (end_fixed)
pin(n - 1, Context.end_elevation);
for (std::size_t k = 1; k < n; ++k)
{
auto const factor{lower[k] / diagonal[k - 1]};
diagonal[k] -= factor * upper[k - 1];
right[k] -= factor * right[k - 1];
}
std::vector<double> z(n);
z[n - 1] = right[n - 1] / diagonal[n - 1];
for (auto k = n - 1; k-- > 0;)
z[k] = (right[k] - upper[k] * z[k + 1]) / diagonal[k];
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};
// grades kept within the limits, the fixed ends don't move
for (int pass = 0; pass < 4; ++pass)
// 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)
{
for (std::size_t k = 0; k + 1 < n; ++k)
{
if (k + 2 == n && end_fixed)
continue;
auto const run{Chainages[Knots[k + 1]] - Chainages[Knots[k]]};
auto const step{segment_limit(Chainages[Knots[k]], Chainages[Knots[k + 1]]) * run};
z[k + 1] = std::clamp(z[k + 1], z[k] - step, z[k] + step);
}
for (auto k = n - 1; k-- > 0;)
{
if (k == 0 && start_fixed)
continue;
auto const run{Chainages[Knots[k + 1]] - Chainages[Knots[k]]};
auto const step{segment_limit(Chainages[Knots[k]], Chainages[Knots[k + 1]]) * run};
z[k] = std::clamp(z[k], z[k + 1] - step, z[k + 1] + step);
}
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;
};
auto const radius{std::ceil(required_radius(result) / 10.0) * 10.0};
auto const threshold{result.norms.curve_threshold / 1000.0};
std::vector<double> z;
for (int iteration = 0; iteration < 1000; ++iteration)
// 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 (;;)
{
z = solve(knots);
auto const n{knots.size()};
if (n <= 2)
break;
auto const chainage = [&](std::size_t const K) { return Chainages[knots[K]]; };
auto const slope = [&](std::size_t const K) { return (z[K + 1] - z[K]) / std::max(1e-9, chainage(K + 1) - chainage(K)); };
std::vector<double> change(n, 0.0), tangent(n, 0.0);
for (std::size_t k = 1; k + 1 < n; ++k)
auto current{cost(x, true)};
for (int pass = 0; pass < 100; ++pass)
{
change[k] = std::abs(slope(k) - slope(k - 1));
tangent[k] = change[k] > threshold + 1e-9 ? radius * change[k] * 0.5 : 0.0;
}
// the curve of a joint starts at the end of the route and takes its whole length
if (start_fixed)
{
change[0] = std::abs(slope(0) - Context.start_grade);
tangent[0] = change[0] > threshold + 1e-9 ? radius * change[0] : 0.0;
}
if (end_fixed)
{
change[n - 1] = std::abs(slope(n - 2) - Context.end_grade);
tangent[n - 1] = change[n - 1] > threshold + 1e-9 ? radius * change[n - 1] : 0.0;
}
std::size_t removed{0};
for (std::size_t k = 1; k + 1 < n && removed == 0; ++k)
for (auto const &zone : Context.switches)
if (overlaps(chainage(k) - tangent[k] - 1.0, chainage(k) + tangent[k] + 1.0, zone))
{
removed = k;
break;
}
if (removed == 0)
{
double worst{0.0};
for (std::size_t k = 0; k + 1 < n; ++k)
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)
{
bool const inner{(k > 0 || tangent[0] > 0.0) && (k + 2 < n || tangent[n - 1] > 0.0)};
auto const needed{tangent[k] + tangent[k + 1] + (inner ? result.norms.element_min : 0.0)};
auto const shortage{needed - (chainage(k + 1) - chainage(k))};
if (shortage <= 1e-6 || shortage <= worst)
continue;
worst = shortage;
if (k == 0)
removed = 1;
else if (k + 2 == n)
removed = n - 2;
else
removed = change[k] <= change[k + 1] ? k : k + 1;
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 (removed == 0)
if (spread <= 1.0 && flat <= 1e-6)
break;
knots.erase(knots.begin() + removed);
spread = std::max(1.0, spread * 0.5);
flat = std::max(1e-6, flat * 0.5);
}
for (std::size_t k = 0; k < knots.size(); ++k)
// 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;)
{
pvi point;
point.chainage = Chainages[knots[k]];
point.elevation = z[k];
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;
}