vertexdepthreconstruct.h

Go to the documentation of this file.

 
#ifndef VERTEX_DEPTH_RECONSTRUCT_H
#define VERTEX_DEPTH_RECONSTRUCT_H
 
#include <algorithm>
#include <array>
#include <cmath>
#include <vector>
 
namespace VertexDepthReconstruct
{
 
inline double cellEtaFromMeanDepth(double h, double za, double zb, double zc)
{
    double z1 = za, z2 = zb, z3 = zc;
    if (z1 > z2) std::swap(z1, z2);
    if (z2 > z3) std::swap(z2, z3);
    if (z1 > z2) std::swap(z1, z2);
 
    if (!(h > 0.0)) return z1;
 
    const double zbar   = (z1 + z2 + z3) / 3.0;
    const double relief = z3 - z1;
    if (relief < 1.0e-9 || h >= z3 - zbar)
        return zbar + h;                              // flat / fully wet
 
    const double h_at_z2 = (z2 - z1) * (z2 - z1) / (3.0 * relief);
    if (h <= h_at_z2)                                 // waterline below z2
        return z1 + std::cbrt(3.0 * h * (z2 - z1) * relief);
 
    // Waterline between z2 and z3: safeguarded Newton on the bracket.
    const double denom = 3.0 * relief * (z3 - z2);
    double lo = z2, hi = z3;
    double eta = zbar + h;
    if (eta <= lo || eta >= hi) eta = 0.5 * (lo + hi);
    for (int it = 0; it < 64; ++it) {
        const double dz3 = z3 - eta;
        const double f  = (eta - zbar) + dz3 * dz3 * dz3 / denom - h;
        if (f > 0.0) hi = eta; else lo = eta;
        const double df = 1.0 - dz3 * dz3 / (relief * (z3 - z2));
        double next = (df > 1.0e-12) ? eta - f / df : 0.5 * (lo + hi);
        if (next <= lo || next >= hi) next = 0.5 * (lo + hi);
        if (std::abs(next - eta) < 1.0e-12 * (1.0 + relief)) return next;
        eta = next;
    }
    return eta;
}
 
inline double triMeanDepthFromEta(double eta, double za, double zb, double zc)
{
    double z1 = za, z2 = zb, z3 = zc;
    if (z1 > z2) std::swap(z1, z2);
    if (z2 > z3) std::swap(z2, z3);
    if (z1 > z2) std::swap(z1, z2);
    if (eta <= z1) return 0.0;
    const double relief = z3 - z1;
    const double zbar   = (z1 + z2 + z3) / 3.0;
    if (relief < 1.0e-9 || eta >= z3) return eta - zbar;
    if (eta <= z2) {
        const double d = eta - z1;
        return d * d * d / (3.0 * (z2 - z1) * relief);
    }
    const double dz3 = z3 - eta;
    return (eta - zbar) + dz3 * dz3 * dz3 / (3.0 * relief * (z3 - z2));
}
 
struct CellSplit
{
    std::array<int, 4>                 v{{-1, -1, -1, -1}}; 
    int                                nSub = 1;            
    std::array<std::array<int, 3>, 2>  sub{};               
    std::array<double, 2>              area{{0.0, 0.0}};    
 
    int vertexCount() const noexcept { return v[3] >= 0 ? 4 : 3; }
};
 
inline double quadEtaFromMeanDepth(const double zs[6], double a1, double a2, double h)
{
    const double A = a1 + a2;
    if (!(A > 0.0)) {
        // Degenerate quad — fall back to the flat closure over its mean bed.
        double zm = 0.0;
        for (int k = 0; k < 6; ++k) zm += zs[k];
        return zm / 6.0 + ((h > 0.0) ? h : 0.0);
    }
    const double zbar1 = (zs[0] + zs[1] + zs[2]) / 3.0;
    const double zbar2 = (zs[3] + zs[4] + zs[5]) / 3.0;
    const double zw    = (a1 * zbar1 + a2 * zbar2) / A;
    double zlow = zs[0], ztop = zs[0];
    for (int k = 1; k < 6; ++k) {
        if (zs[k] < zlow) zlow = zs[k];
        if (zs[k] > ztop) ztop = zs[k];
    }
    if (!(h > 0.0)) return zlow;
    const double relief = ztop - zlow;
    if (relief < 1.0e-9 || h >= ztop - zw) return zw + h;      // flat / fully wet
 
    auto meanDepth = [&](double eta) {
        return (a1 * triMeanDepthFromEta(eta, zs[0], zs[1], zs[2])
              + a2 * triMeanDepthFromEta(eta, zs[3], zs[4], zs[5])) / A;
    };
    double lo = zlow, hi = ztop;
    for (int it = 0; it < 64; ++it) {
        const double mid = 0.5 * (lo + hi);
        if (meanDepth(mid) < h) lo = mid; else hi = mid;
        if (hi - lo < 1.0e-12 * (1.0 + relief)) break;
    }
    return 0.5 * (lo + hi);
}
 
inline void reconstructVertexSignedDepths(
    const std::vector<CellSplit>& cells,
    const std::vector<float>&  cellDepths,
    const std::vector<float>&  cellZc,
    const std::vector<double>& vz,
    float dryF,
    std::vector<float>& vsum,
    std::vector<float>& wsum,
    std::vector<float>& outVertexDepth)
{
    const int nVert = static_cast<int>(vz.size());
    vsum.assign(static_cast<size_t>(nVert), 0.0f);
    wsum.assign(static_cast<size_t>(nVert), 0.0f);
    const int nCell = std::min<int>(static_cast<int>(cells.size()),
                                    static_cast<int>(cellDepths.size()));
    for (int i = 0; i < nCell; ++i) {
        const float h = cellDepths[i];
        // NaN-robust dry skip: `h < dryF` is false for NaN, so a non-finite
        // depth would NOT be skipped and would poison vsum/wsum at all of the
        // cell's vertices (→ streaked triangle fans in the Gouraud fill).
        if (!(h >= dryF)) continue;                // only wetted cells contribute
        const CellSplit& c = cells[i];
        const int nv = c.vertexCount();
        bool zOk = true;
        for (int k = 0; k < nv; ++k) {
            const int vi = c.v[k];
            if (vi < 0 || vi >= nVert || !std::isfinite(vz[vi])) { zOk = false; break; }
        }
        // Cell free surface via the planar-bed stage–storage inversion when
        // the vertex elevations are usable; flat closure z_c + h as the
        // fallback (out-of-range index / nodata z).
        double eta;
        if (!zOk) {
            eta = double(cellZc[i]) + double(h);
        } else if (nv == 3) {
            eta = cellEtaFromMeanDepth(double(h), vz[c.v[0]], vz[c.v[1]], vz[c.v[2]]);
        } else {
            const double zs[6] = { vz[c.sub[0][0]], vz[c.sub[0][1]], vz[c.sub[0][2]],
                                   vz[c.sub[1][0]], vz[c.sub[1][1]], vz[c.sub[1][2]] };
            eta = quadEtaFromMeanDepth(zs, c.area[0], c.area[1], double(h));
        }
        const float w  = (nv == 3) ? h : h * 0.75f;  // depth weight × 3/nv
        const float we = w * float(eta);            // weighted η contribution
        if (!std::isfinite(we)) continue;           // non-finite z_c must not spread
        for (int k = 0; k < nv; ++k) {
            const int vi = c.v[k];
            if (vi < 0 || vi >= nVert) continue;
            // Wetted-contact gate (mirror of the engine): this cell's water
            // votes at corner vi only if its surface reaches the corner.
            // NaN vz[vi] compares false → skipped, consistent with the
            // non-finite handling at output.
            if (!(eta > vz[vi])) continue;
            vsum[vi] += we;
            wsum[vi] += w;
        }
    }
    outVertexDepth.assign(static_cast<size_t>(nVert), 0.0f);
    for (int v = 0; v < nVert; ++v)
        if (wsum[v] > 0.0f) {
            const double d = double(vsum[v]) / double(wsum[v]) - vz[v];
            // Non-finite vertex elevation (e.g. DTM nodata) must yield a dry
            // vertex, not a NaN that the colour ramp turns into garbage.
            outVertexDepth[v] = std::isfinite(d) ? float(d) : 0.0f;
        }
}
 
inline void reconstructVertexSignedDepths(
    const std::vector<std::array<int, 3>>& tris,
    const std::vector<float>&  cellDepths,
    const std::vector<float>&  cellZc,
    const std::vector<double>& vz,
    float dryF,
    std::vector<float>& vsum,
    std::vector<float>& wsum,
    std::vector<float>& outVertexDepth)
{
    std::vector<CellSplit> cells(tris.size());
    for (size_t i = 0; i < tris.size(); ++i) {
        cells[i].v      = {tris[i][0], tris[i][1], tris[i][2], -1};
        cells[i].sub[0] = tris[i];
    }
    reconstructVertexSignedDepths(cells, cellDepths, cellZc, vz, dryF,
                                  vsum, wsum, outVertexDepth);
}
 
inline void extrapolateDryCorners(double z0, double z1, double z2,
                                  float& sd0, float& sd1, float& sd2)
{
    float* const sd[3] = { &sd0, &sd1, &sd2 };
    const double z[3]  = { z0, z1, z2 };
 
    bool   any    = false;
    double maxEta = 0.0;
    for (int k = 0; k < 3; ++k) {
        if (!(*sd[k] > 0.0f)) continue;
        const double e = z[k] + double(*sd[k]);
        if (!std::isfinite(e)) continue;
        if (!any || e > maxEta) { maxEta = e; any = true; }
    }
    if (!any) return;                       // fully dry — nothing to extend
 
    for (int k = 0; k < 3; ++k) {
        if (*sd[k] > 0.0f) continue;        // wet corner keeps its own η
        if (!std::isfinite(z[k])) continue; // nodata bed stays the sentinel
        const double d = maxEta - z[k];
        // ADVERSE SLOPE ONLY. A dry corner standing above the driving head is
        // the pooling case: the surface runs level into the cell and meets the
        // rising bed at the intercept, so d <= 0 carries that geometry. A dry
        // corner BELOW the driving head is the opposite — the bed falls away
        // from the water, which is where the solver's dryness is meaningful
        // (the water drained, or never arrived). Extrapolating there would
        // inject d metres of standing water into a dry cell and spread the
        // pool downhill one cell in every direction, which is the flood-fill
        // behaviour this feature explicitly does not do. Leave the sentinel.
        if (d > 0.0) continue;
        *sd[k] = float(d);
    }
}
 
} // namespace VertexDepthReconstruct
 
#endif // VERTEX_DEPTH_RECONSTRUCT_H