meshprofileinterp.h

Go to the documentation of this file.

 
#ifndef MESH_PROFILE_INTERP_H
#define MESH_PROFILE_INTERP_H
 
#include "plot/meshprofilesampler.h"
 
#include <QVector>
 
#include <cmath>
#include <functional>
#include <limits>
 
namespace MeshProfileInterp
{
 
inline constexpr double kFilm = 1e-6;
 
inline QVector<double> bridgedTops(
    const QVector<MeshProfileSampler::Sample> &s,
    const std::function<double(const MeshProfileSampler::Sample &)> &topElev)
{
    const int n = s.size();
    const double kNaN = std::numeric_limits<double>::quiet_NaN();
    QVector<double> top(n, kNaN);
 
    // 1) Raw per-sample surface — films included via kFilm.
    for (int i = 0; i < n; ++i) {
        if (!std::isfinite(s[i].ground)) continue;
        const double t = topElev(s[i]);
        if (t > s[i].ground + kFilm) top[i] = t;
    }
 
    // 2) Bridge dry gaps that lie between two wet runs.
    int i = 0;
    while (i < n) {
        if (std::isnan(top[i])) { ++i; continue; }
 
        // End of the current wet run.
        int runEnd = i;
        while (runEnd + 1 < n && !std::isnan(top[runEnd + 1])) ++runEnd;
 
        // Walk the gap to the next wet sample. An off-mesh (NaN ground) sample
        // blocks the bridge — a genuine mesh hole must stay a true gap. A gap
        // whose samples ALL carry a valid free surface (cellHasSurface) is
        // genuinely dry ground, not missing data — never bridge it.
        const int gapStart = runEnd + 1;
        int next = gapStart;
        bool blocked = false;
        bool sawNoData = false;
        while (next < n && std::isnan(top[next])) {
            if (!std::isfinite(s[next].ground)) { blocked = true; break; }
            if (!s[next].cellHasSurface) sawNoData = true;
            ++next;
        }
 
        if (!blocked && sawNoData && next < n && gapStart < next) {
            const double cL = s[runEnd].chainage, wseL = top[runEnd];
            const double cR = s[next].chainage,    wseR = top[next];
            const double dC = cR - cL;
            const double wseUp = std::max(wseL, wseR);   // upstream (higher) level
            for (int k = gapStart; k < next; ++k) {
                if (!std::isfinite(s[k].ground)) continue;   // safety
                double wse = (dC > 1e-12)
                    ? wseL + (s[k].chainage - cL) / dC * (wseR - wseL)
                    : wseL;
                if (wse > wseUp) wse = wseUp;                // no uphill water
                if (wse > s[k].ground + kFilm)               // above the bed → wet
                    top[k] = wse;
                // else: dry crest — leave NaN so the run splits at the hill.
            }
        }
 
        i = blocked ? (runEnd + 1) : std::max(next, runEnd + 1);
    }
 
    return top;
}
 
inline bool shorelineIntercept(const QVector<MeshProfileSampler::Sample> &s,
                               const QVector<double> &top,
                               int runFirst, int runLast, bool trailing,
                               double *outChainage, double *outElev)
{
    const int n = s.size();
    if (runFirst < 0 || runLast >= n || runFirst > runLast || top.size() != n)
        return false;
    const int i = trailing ? runLast : runFirst;          // boundary wet sample
    const int j = trailing ? runLast - 1 : runFirst + 1;  // inner wet neighbour
    const int k = trailing ? runLast + 1 : runFirst - 1;  // adjacent dry sample
    if (k < 0 || k >= n) return false;                    // run ends at the data edge
    if (!std::isfinite(s[k].ground)) return false;        // off-mesh — keep the hard edge
    if (!std::isnan(top[k])) return false;                // neighbour not dry (defensive)
    if (std::isnan(top[i]) || !std::isfinite(s[i].ground)) return false;
 
    const double ci = s[i].chainage, ck = s[k].chainage;
    const double dSeg = ck - ci;                          // signed, toward the dry side
    if (std::abs(dSeg) < 1e-12) return false;
 
    // Wet-side surface slope from the two boundary-most wet samples; a
    // single-sample run extends flat.
    double mTop = 0.0;
    if (j >= runFirst && j <= runLast) {
        const double dcw = ci - s[j].chainage;
        if (std::abs(dcw) > 1e-12 && !std::isnan(top[j]))
            mTop = (top[i] - top[j]) / dcw;
    }
 
    // Along t ∈ [0,1] from sample i to sample k:
    //   surface(t) = top_i + mTop·dSeg·t
    //   ground(t)  = g_i + (g_k − g_i)·t
    // They cross where surface == ground.
    const double gi = s[i].ground, gk = s[k].ground;
    const double denom = (gk - gi) - mTop * dSeg;
    if (denom <= 1e-12) return false;      // ground never rises through the surface
    const double t = (top[i] - gi) / denom;
    if (t <= 0.0 || t > 1.0) return false; // no crossing inside the segment
 
    if (outChainage) *outChainage = ci + t * dSeg;
    if (outElev)     *outElev     = gi + t * (gk - gi);
    return true;
}
 
} // namespace MeshProfileInterp
 
#endif // MESH_PROFILE_INTERP_H