TYAcousticModel9613Solver2024.cpp

Go to the documentation of this file.

/*
 * Copyright (C) <2012-2024> <EDF-DTG> <FRANCE>
 * This file is part of Code_TYMPAN (R).
 * Code_TYMPAN (R) is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * Code_TYMPAN (R) is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License along
 * with Code_TYMPAN (R). If not, see <https://www.gnu.org/licenses/>.
 */
 
#include "TYAcousticModel9613Solver2024.h"
#include <fstream>
#include <string>
 
namespace
{
size_t readReflectionOrderFromConfig(const std::string& filePath, size_t defaultValue)
{
    std::ifstream file(filePath);
    if (!file.is_open())
    {
        return defaultValue;
    }
 
    const std::string key = "REFLECTION_ORDER=";
    std::string line;
 
    while (std::getline(file, line))
    {
        if (line.compare(0, key.size(), key) == 0)
        {
            try
            {
                return static_cast<size_t>(std::stoul(line.substr(key.size())));
            }
            catch (...)
            {
                return defaultValue;
            }
        }
    }
 
    return defaultValue;
}
 
bool intersectSegmentsFastInTheirPlane(const OSegment3D& segA, const OSegment3D& segB, OPoint3D& pt,
                                       double seuilConfondus)
{
    const OPoint3D& A = segA._ptA;
    const OPoint3D& B = segA._ptB;
    const OPoint3D& C = segB._ptA;
    const OPoint3D& D = segB._ptB;
 
    const double ux = B._x - A._x;
    const double uy = B._y - A._y;
    const double uz = B._z - A._z;
 
    const double vx = D._x - C._x;
    const double vy = D._y - C._y;
    const double vz = D._z - C._z;
 
    // Normal of the plane spanned by the 2 segment directions.
    const double nx = uy * vz - uz * vy;
    const double ny = uz * vx - ux * vz;
    const double nz = ux * vy - uy * vx;
 
    const double anx = std::abs(nx);
    const double any = std::abs(ny);
    const double anz = std::abs(nz);
 
    // Degenerate or quasi-parallel case: keep legacy robust behavior.
    if ((anx <= TYSEUILCONFONDUS) && (any <= TYSEUILCONFONDUS) && (anz <= TYSEUILCONFONDUS))
    {
        return segA.intersects(segB, pt, seuilConfondus) != INTERS_NULLE;
    }
 
    // Project to the coordinate plane where the 2D determinant is the most stable.
    double a1, a2, u1, u2, c1, c2, v1, v2;
 
    if ((anx >= any) && (anx >= anz))
    {
        // Drop X => work in YZ
        a1 = A._y;
        a2 = A._z;
        u1 = uy;
        u2 = uz;
        c1 = C._y;
        c2 = C._z;
        v1 = vy;
        v2 = vz;
    }
    else if ((any >= anx) && (any >= anz))
    {
        // Drop Y => work in XZ
        a1 = A._x;
        a2 = A._z;
        u1 = ux;
        u2 = uz;
        c1 = C._x;
        c2 = C._z;
        v1 = vx;
        v2 = vz;
    }
    else
    {
        // Drop Z => work in XY
        a1 = A._x;
        a2 = A._y;
        u1 = ux;
        u2 = uy;
        c1 = C._x;
        c2 = C._y;
        v1 = vx;
        v2 = vy;
    }
 
    const double w1 = c1 - a1;
    const double w2 = c2 - a2;
 
    const double det = u1 * v2 - u2 * v1;
 
    // Near-parallel in projected plane: keep legacy behavior.
    if (std::abs(det) <= TYSEUILCONFONDUS)
    {
        return segA.intersects(segB, pt, seuilConfondus) != INTERS_NULLE;
    }
 
    const double t = (w1 * v2 - w2 * v1) / det;
    const double s = (w1 * u2 - w2 * u1) / det;
 
    const double lenU = std::sqrt(ux * ux + uy * uy + uz * uz);
    const double lenV = std::sqrt(vx * vx + vy * vy + vz * vz);
    const double paramTol = seuilConfondus / std::max(1.0, std::max(lenU, lenV));
 
    if ((t < -paramTol) || (t > 1.0 + paramTol) || (s < -paramTol) || (s > 1.0 + paramTol))
    {
        return false;
    }
 
    pt._x = A._x + t * ux;
    pt._y = A._y + t * uy;
    pt._z = A._z + t * uz;
 
    return true;
}
} // namespace
 
TYAcousticModel9613Solver2024::TYAcousticModel9613Solver2024(TYSolver9613Solver2024& solver)
    : TYAcousticModel9613Solver(solver)
{
}
 
OSpectreOctave TYAcousticModel9613Solver2024::calculZMin(const double C2, const OSpectreOctave& C3) const
{
    // zmin = -2 * lambda / (C2 * C3)
    OSpectreOctave zmin = _lambda * (-2.0);
    zmin = zmin.div(C3 * C2);
    return zmin;
}
 
OSpectreOctave TYAcousticModel9613Solver2024::calculKmeteo(const bool vertical, const double d_SS,
                                                           const double d_SR, const double d, const double z,
                                                           const double e, const OSpectreOctave& z_min) const
{
    // K_meteo = exp{-(1 / 2000) * sqrt(([max(d_SS, d_SR) + e] * min(d_SS, d_SR) * d) / [2 * (z - z_min)])}
    // for vertical diffraction K_meteo = 1 for lateral diffraction
    //
    // This is not specified explicitly in the standard, but we also set:
    // K_meteo = 1 when z < z_min
    // which avoids a division by zero and has no effect on the computation of Dz
 
    OSpectreOctave K_meteo{1.0};
 
    if (!vertical)
    {
        return K_meteo;
    }
 
    double* kMeteo = K_meteo.getTabValReel();
    const double* zMin = z_min.getTabValReel();
 
    const double dMax = std::max(d_SS, d_SR);
    const double dMin = std::min(d_SS, d_SR);
    const double numerator = (dMax + e) * dMin * d;
 
    for (std::size_t i = 0; i < TY_SPECTRE_OCT_NB_ELMT; ++i)
    {
        if (z > zMin[i])
        {
            kMeteo[i] = std::exp(-(1.0 / 2000.0) * std::sqrt(numerator / (2.0 * (z - zMin[i]))));
        }
    }
 
    return K_meteo;
}
 
OSpectreOctave TYAcousticModel9613Solver2024::calculDz(const double z, const double C2,
                                                       const OSpectreOctave& C3, const OSpectreOctave& Kmeteo,
                                                       const OSpectreOctave& zmin) const
{
    // Dz = 10 * log10[1 + (2 + (C2 / lambda) * C3 * z) * Kmeteo] dB for z > zmin
    // Dz = 0 dB for z <= zmin
 
    OSpectreOctave Dz{0.0};
 
    double* dz = Dz.getTabValReel();
    const double* lambda = _lambda.getTabValReel();
    const double* c3 = C3.getTabValReel();
    const double* kmeteo = Kmeteo.getTabValReel();
    const double* zMin = zmin.getTabValReel();
 
    for (std::size_t i = 0; i < TY_SPECTRE_OCT_NB_ELMT; ++i)
    {
        if (z > zMin[i])
        {
            const double term = 2.0 + (C2 / lambda[i]) * c3[i] * z;
            dz[i] = 10.0 * std::log10(1.0 + term * kmeteo[i]);
        }
    }
 
    return Dz;
}
 
void TYAcousticModel9613Solver2024::computeCheminReflexion(const std::deque<TYSIntersection>& tabIntersect,
                                                           const OSegment3D& ray,
                                                           const tympan::AcousticSource& source,
                                                           TYTabChemin9613Solver& TabChemins,
                                                           double distance) const
{
    // Multiple reflection order
    // TODO : Externalize in computation parameters
    const std::string solverConfigPath = "C:\\projects\\tympan_tools\\17534-3\\solver_conf.txt";
    const size_t defaultReflectionOrder = 1;
    const size_t multipleReflectionOrder =
        readReflectionOrderFromConfig(solverConfigPath, defaultReflectionOrder);
 
    if (multipleReflectionOrder == 0)
    {
        return;
    }
 
    std::vector<ReflectingSegmentCache> reflectingSegments;
    reflectingSegments.reserve(tabIntersect.size());
 
    // Keep only barriers that are both infrastructure and intersecting the EL plane.
    // Their support-line geometry is cached once here, then reused during the DFS.
    for (const auto& inter : tabIntersect)
    {
        if (inter.isInfra && inter.bIntersect[1])
        {
            reflectingSegments.push_back(makeReflectingSegmentCache(inter));
        }
    }
 
    if (reflectingSegments.empty())
    {
        return;
    }
 
    std::vector<const TYSIntersection*> currentCombination;
    currentCombination.reserve(multipleReflectionOrder);
 
    std::deque<OPoint3D> currentImageSources;
    currentImageSources.push_back(ray._ptA);
 
    // Streaming exploration of ordered reflection combinations.
    // Each valid prefix is turned directly into an acoustic path, avoiding the
    // construction of an intermediate global set of combinations/candidates.
    buildReflectionPathsStreaming(multipleReflectionOrder, reflectingSegments, currentCombination,
                                  currentImageSources, tabIntersect, ray, source, TabChemins, distance);
}
 
bool TYAcousticModel9613Solver2024::getGroundfactors(
    const std::deque<std::deque<TYSIntersection>>& tabIntersectSegments,
    const std::deque<OPoint3D>& pathPoints2D, double hs, double hr, double& Gs, double& Gm, double& Gr) const
{
    Gs = 0.5;
    Gm = 0.5;
    Gr = 0.5;
 
    if (!computeGroundFactorSourceZone(tabIntersectSegments, pathPoints2D, hs, Gs))
    {
        return false;
    }
 
    if (!computeGroundFactorMiddleZone(tabIntersectSegments, pathPoints2D, hs, hr, Gm))
    {
        return false;
    }
 
    if (!computeGroundFactorReceiverZone(tabIntersectSegments, pathPoints2D, hr, Gr))
    {
        return false;
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::computeGroundFactorSourceZone(
    const std::deque<std::deque<TYSIntersection>>& tabIntersectSegments,
    const std::deque<OPoint3D>& pathPoints2D, double hs, double& Gs) const
{
    double heightRatio = 30.0;
    Gs = 0.5;
 
    if (pathPoints2D.size() < 2)
    {
        return false;
    }
 
    if (tabIntersectSegments.size() != pathPoints2D.size() - 1)
    {
        return false;
    }
 
    const double srcZoneLength = heightRatio * hs;
 
    // Same fallback behavior as legacy code.
    if (srcZoneLength <= TYSEUILCONFONDUS)
    {
        Gs = 0.5;
        return true;
    }
 
    double traversedLength = 0.0; // curvilinear length along the broken path
    double coveredLength = 0.0;   // sum of dp returned by computeGZone
    double weightedGroundSum = 0.0;
 
    for (size_t i = 0; i + 1 < pathPoints2D.size(); ++i)
    {
        const OPoint3D& ptA = pathPoints2D[i];
        const OPoint3D& ptB = pathPoints2D[i + 1];
 
        const OSegment3D currentSegment(ptA, ptB);
        const double currentSegmentLength = currentSegment.longueur();
 
        if (currentSegmentLength <= TYSEUILCONFONDUS)
        {
            continue;
        }
 
        const double remainingLength = srcZoneLength - traversedLength;
 
        if (remainingLength <= TYSEUILCONFONDUS)
        {
            break;
        }
 
        OPoint3D ptEndZone = ptB;
 
        // The source-zone limit lies inside the current segment.
        if (remainingLength < currentSegmentLength)
        {
            const double alpha = remainingLength / currentSegmentLength;
 
            ptEndZone._x = ptA._x + alpha * (ptB._x - ptA._x);
            ptEndZone._y = ptA._y + alpha * (ptB._y - ptA._y);
            ptEndZone._z = ptA._z + alpha * (ptB._z - ptA._z);
 
            double GCurSeg = 0.0;
            double dpCurSeg = 0.0;
            computeGZone(ptA, ptEndZone, GCurSeg, dpCurSeg, tabIntersectSegments[i]);
 
            weightedGroundSum += GCurSeg;
            coveredLength += dpCurSeg;
            traversedLength = srcZoneLength;
 
            break;
        }
 
        // The whole current segment belongs to the source zone.
        double GCurSeg = 0.0;
        double dpCurSeg = 0.0;
        computeGZone(ptA, ptB, GCurSeg, dpCurSeg, tabIntersectSegments[i]);
 
        weightedGroundSum += GCurSeg;
        coveredLength += dpCurSeg;
        traversedLength += currentSegmentLength;
    }
 
    if (coveredLength > TYSEUILCONFONDUS)
    {
        Gs = weightedGroundSum / coveredLength;
    }
    else
    {
        Gs = 0.5;
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::computeGroundFactorReceiverZone(
    const std::deque<std::deque<TYSIntersection>>& tabIntersectSegments,
    const std::deque<OPoint3D>& pathPoints2D, double hr, double& Gr) const
{
    double heightRatio = 30.0;
    Gr = 0.5;
 
    if (pathPoints2D.size() < 2)
    {
        return false;
    }
 
    if (tabIntersectSegments.size() != pathPoints2D.size() - 1)
    {
        return false;
    }
 
    const double rcpZoneLength = heightRatio * hr;
 
    // Same fallback behavior as legacy code.
    if (rcpZoneLength <= TYSEUILCONFONDUS)
    {
        Gr = 0.5;
        return true;
    }
 
    double traversedLength = 0.0; // curvilinear length traversed backward from receiver
    double coveredLength = 0.0;   // sum of dp returned by computeGZone
    double weightedGroundSum = 0.0;
 
    for (size_t i = pathPoints2D.size() - 1; i > 0; --i)
    {
        const OPoint3D& ptA = pathPoints2D[i - 1];
        const OPoint3D& ptB = pathPoints2D[i];
 
        const OSegment3D currentSegment(ptA, ptB);
        const double currentSegmentLength = currentSegment.longueur();
 
        if (currentSegmentLength <= TYSEUILCONFONDUS)
        {
            continue;
        }
 
        const double remainingLength = rcpZoneLength - traversedLength;
 
        if (remainingLength <= TYSEUILCONFONDUS)
        {
            break;
        }
 
        // The receiver-zone limit lies inside the current segment.
        if (remainingLength < currentSegmentLength)
        {
            const double alpha = remainingLength / currentSegmentLength;
 
            // Start point of the receiver zone on the current segment,
            // located at distance "remainingLength" from ptB toward ptA.
            OPoint3D ptStartZone;
            ptStartZone._x = ptB._x + alpha * (ptA._x - ptB._x);
            ptStartZone._y = ptB._y + alpha * (ptA._y - ptB._y);
            ptStartZone._z = ptB._z + alpha * (ptA._z - ptB._z);
 
            double GCurSeg = 0.0;
            double dpCurSeg = 0.0;
            computeGZone(ptStartZone, ptB, GCurSeg, dpCurSeg, tabIntersectSegments[i - 1]);
 
            weightedGroundSum += GCurSeg;
            coveredLength += dpCurSeg;
            traversedLength = rcpZoneLength;
 
            break;
        }
 
        // The whole current segment belongs to the receiver zone.
        double GCurSeg = 0.0;
        double dpCurSeg = 0.0;
        computeGZone(ptA, ptB, GCurSeg, dpCurSeg, tabIntersectSegments[i - 1]);
 
        weightedGroundSum += GCurSeg;
        coveredLength += dpCurSeg;
        traversedLength += currentSegmentLength;
    }
 
    if (coveredLength > TYSEUILCONFONDUS)
    {
        Gr = weightedGroundSum / coveredLength;
    }
    else
    {
        Gr = 0.5;
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::computeGroundFactorMiddleZone(
    const std::deque<std::deque<TYSIntersection>>& tabIntersectSegments,
    const std::deque<OPoint3D>& pathPoints2D, double hs, double hr, double& Gm) const
{
    double heightRatio = 30.0;
    Gm = 0.5;
 
    if (pathPoints2D.size() < 2)
    {
        return false;
    }
 
    if (tabIntersectSegments.size() != pathPoints2D.size() - 1)
    {
        return false;
    }
 
    const double srcZoneLength = heightRatio * hs;
    const double rcpZoneLength = heightRatio * hr;
 
    // Compute total curvilinear length of the broken path.
    double dp = 0.0;
    for (size_t i = 0; i + 1 < pathPoints2D.size(); ++i)
    {
        dp += OSegment3D(pathPoints2D[i], pathPoints2D[i + 1]).longueur();
    }
 
    // No middle zone when source and receiver zones touch or overlap.
    if ((srcZoneLength + rcpZoneLength) >= (dp - TYSEUILCONFONDUS))
    {
        Gm = 0.5;
        return true;
    }
 
    const double middleZoneBegin = srcZoneLength;
    const double middleZoneEnd = dp - rcpZoneLength;
 
    double traversedLength = 0.0; // curvilinear abscissa at segment start
    double coveredLength = 0.0;   // sum of dp returned by computeGZone
    double weightedGroundSum = 0.0;
 
    for (size_t i = 0; i + 1 < pathPoints2D.size(); ++i)
    {
        const OPoint3D& ptA = pathPoints2D[i];
        const OPoint3D& ptB = pathPoints2D[i + 1];
 
        const OSegment3D currentSegment(ptA, ptB);
        const double currentSegmentLength = currentSegment.longueur();
 
        if (currentSegmentLength <= TYSEUILCONFONDUS)
        {
            continue;
        }
 
        const double segmentBegin = traversedLength;
        const double segmentEnd = traversedLength + currentSegmentLength;
 
        // Overlap between current segment and middle zone in curvilinear abscissa.
        const double overlapBegin = std::max(segmentBegin, middleZoneBegin);
        const double overlapEnd = std::min(segmentEnd, middleZoneEnd);
 
        if ((overlapEnd - overlapBegin) > TYSEUILCONFONDUS)
        {
            const double alphaBegin = (overlapBegin - segmentBegin) / currentSegmentLength;
            const double alphaEnd = (overlapEnd - segmentBegin) / currentSegmentLength;
 
            OPoint3D ptZoneBegin;
            ptZoneBegin._x = ptA._x + alphaBegin * (ptB._x - ptA._x);
            ptZoneBegin._y = ptA._y + alphaBegin * (ptB._y - ptA._y);
            ptZoneBegin._z = ptA._z + alphaBegin * (ptB._z - ptA._z);
 
            OPoint3D ptZoneEnd;
            ptZoneEnd._x = ptA._x + alphaEnd * (ptB._x - ptA._x);
            ptZoneEnd._y = ptA._y + alphaEnd * (ptB._y - ptA._y);
            ptZoneEnd._z = ptA._z + alphaEnd * (ptB._z - ptA._z);
 
            double GCurSeg = 0.0;
            double dpCurSeg = 0.0;
            computeGZone(ptZoneBegin, ptZoneEnd, GCurSeg, dpCurSeg, tabIntersectSegments[i]);
 
            weightedGroundSum += GCurSeg;
            coveredLength += dpCurSeg;
        }
 
        traversedLength = segmentEnd;
    }
 
    if (coveredLength > TYSEUILCONFONDUS)
    {
        Gm = weightedGroundSum / coveredLength;
    }
    else
    {
        Gm = 0.5;
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::sameReflectingSegment(const TYSIntersection& lhs,
                                                          const TYSIntersection& rhs) const
{
    const OSegment3D& seg1 = lhs.segInter[1];
    const OSegment3D& seg2 = rhs.segInter[1];
 
    const bool sameOrder = (seg1._ptA == seg2._ptA) && (seg1._ptB == seg2._ptB);
 
    const bool reverseOrder = (seg1._ptA == seg2._ptB) && (seg1._ptB == seg2._ptA);
 
    return sameOrder || reverseOrder;
}
 
bool TYAcousticModel9613Solver2024::buildReflectionPoints(
    const std::vector<const TYSIntersection*>& currentCombination,
    const std::deque<OPoint3D>& imageSourcesList, const OPoint3D& receptorPoint,
    std::deque<OPoint3D>& reflectionPointsList) const
{
    reflectionPointsList.clear();
 
    if (currentCombination.empty())
    {
        return true;
    }
 
    // By convention:
    //   imageSourcesList[0] = S0 = source point
    //   imageSourcesList[i] = image source after i reflections
    if (imageSourcesList.size() != currentCombination.size() + 1)
    {
        return false;
    }
 
    OPoint3D nextPoint = receptorPoint; // P(n+1) = R
 
    // Reflection points are reconstructed backward:
    // starting from the receptor, intersect the current image ray with the
    // corresponding reflecting segment, then continue with the previous order.
    for (int i = static_cast<int>(currentCombination.size()) - 1; i >= 0; --i)
    {
        const OPoint3D& currentImageSource = imageSourcesList[i + 1];          // Si
        const OSegment3D& currentBarrier = currentCombination[i]->segInter[1]; // Bi
 
        OPoint3D currentReflectionPoint;
 
        // Intersection must be computed on segments, not on infinite lines.
        if (!intersectSegmentsFastInTheirPlane(currentBarrier, OSegment3D(currentImageSource, nextPoint),
                                               currentReflectionPoint, TYSEUILCONFONDUS))
        {
            reflectionPointsList.clear();
            return false;
        }
 
        reflectionPointsList.push_front(currentReflectionPoint);
        nextPoint = currentReflectionPoint;
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::validateReflectionCandidate(
    const std::vector<const TYSIntersection*>& barrierCombination,
    const std::deque<OPoint3D>& reflectionPointsList, const std::deque<TYSIntersection>& tabIntersect,
    const OPoint3D& sourcePoint, const OPoint3D& receptorPoint) const
{
    const size_t reflectionOrder = barrierCombination.size();
 
    if (reflectionPointsList.size() != reflectionOrder)
    {
        return false;
    }
 
    std::deque<OPoint3D> pathPoints;
    pathPoints.push_back(sourcePoint);
    for (const auto& point : reflectionPointsList)
    {
        pathPoints.push_back(point);
    }
    pathPoints.push_back(receptorPoint);
 
    // For each leg [Pi, P(i+1)], only the reflecting barriers on which belong Pi or P(i+1) are allowed to
    // touch the path. Any other barrier intersection invalidates the candidate.
    for (size_t i = 0; i <= reflectionOrder; ++i)
    {
        const OSegment3D pathSegment(pathPoints[i], pathPoints[i + 1]);
 
        for (const auto& sceneBarrier : tabIntersect)
        {
            if ((i >= 1) && sameReflectingSegment(sceneBarrier, *barrierCombination[i - 1]))
            {
                continue;
            }
 
            if ((i < reflectionOrder) && sameReflectingSegment(sceneBarrier, *barrierCombination[i]))
            {
                continue;
            }
 
            OPoint3D intersectionPoint;
            if (sceneBarrier.segInter[1].intersects(pathSegment, intersectionPoint, TYSEUILCONFONDUS))
            {
                return false;
            }
        }
    }
 
    return true;
}
 
bool TYAcousticModel9613Solver2024::buildReflectionPath(
    const std::vector<const TYSIntersection*>& barrierCombination,
    const std::deque<OPoint3D>& imageSourcesList, const std::deque<OPoint3D>& reflectionPointsList,
    const OSegment3D& directRay, const tympan::AcousticSource& source, TYTabChemin9613Solver& TabChemins,
    double distance) const
{
    const size_t reflectionOrder = barrierCombination.size();
 
    if (reflectionOrder == 0)
    {
        return false;
    }
 
    if (imageSourcesList.size() != reflectionOrder + 1)
    {
        return false;
    }
 
    if (reflectionPointsList.size() != reflectionOrder)
    {
        return false;
    }
 
    // Build path points:
    // P0 = S
    // P1 ... Pn = reflection points
    // P(n+1) = R
    std::deque<OPoint3D> pathPoints;
    std::deque<OPoint3D> pathPoints2D;
    pathPoints.push_back(directRay._ptA);
 
    for (const auto& point : reflectionPointsList)
    {
        pathPoints.push_back(point);
    }
 
    pathPoints.push_back(directRay._ptB);
 
    TYTabEtape9613Solver tabEtapes;
    double pathLength = 0.0;
    double dp = 0.0;
 
    OPoint3D A2D;
    OPoint3D B2D;
    for (size_t i = 0; i + 1 < pathPoints.size(); ++i)
    {
        OSegment3D segment(pathPoints[i], pathPoints[i + 1]);
        pathLength += segment.longueur();
 
        A2D = OPoint3D{segment._ptA._x, segment._ptA._y, 0.0};
        B2D = OPoint3D{segment._ptB._x, segment._ptB._y, 0.0};
 
        // dp is the broken-path length projected in the horizontal plane.
        dp += OSegment3D{A2D, B2D}.longueur();
        pathPoints2D.push_back(A2D);
    }
    pathPoints2D.push_back(B2D);
 
    {
        TYEtape9613Solver etape;
        etape._pt = pathPoints[0];
        etape._type = TYSOURCE;
 
        const OSegment3D firstLeg(pathPoints[0], pathPoints[1]);
        etape._Absorption =
            source.directivity->lwAdjustment(OVector3D(firstLeg._ptA, firstLeg._ptB), firstLeg.longueur());
 
        tabEtapes.push_back(etape);
    }
 
    for (size_t i = 0; i < reflectionOrder; ++i)
    {
        tympan::AcousticBuildingMaterial* material =
            dynamic_cast<tympan::AcousticBuildingMaterial*>(barrierCombination[i]->material);
 
        if (material == nullptr)
        {
            return false;
        }
 
        if (!barrierCombination[i]->pFaceGeomData)
        {
            return false;
        }
 
        TYEtape9613Solver etape;
        etape._pt = reflectionPointsList[i];
        etape._type = TYREFLEXION;
        etape._Absorption = material->spectrum;
 
        tabEtapes.push_back(etape);
    }
 
    OSegment3D segMeanSlope;
    meanSlope(directRay, segMeanSlope);
 
    double hs{0.0}, hr{0.0};
    computeSegmentEdgesHeights(hs, hr, segMeanSlope, directRay);
 
    double Gs{0.0}, Gm{0.0}, Gr{0.0};
 
    // Recompute the scene intersections for each broken-path leg [Pi, P(i+1)].
    // Ground factors must be evaluated on the reflected path, not on the direct ray.
    std::deque<std::deque<TYSIntersection>> tabIntersectSegments;
    OSegment3D seg;
    tabIntersectSegments.resize(reflectionOrder + 1);
    for (size_t i = 0; i < reflectionOrder + 1; ++i)
    {
        seg = OSegment3D{pathPoints[i], pathPoints[i + 1]};
        _solver.selectFaces(tabIntersectSegments[i], seg, source.volume_id);
    }
 
    getGroundfactors(tabIntersectSegments, pathPoints2D, hs, hr, Gs, Gm, Gr);
 
    std::unique_ptr<TYChemin9613Solver> chemin = createChemin();
    chemin->setType(TYTypeChemin::CHEMIN_REFLEX);
    chemin->setLongueur(pathLength);
    chemin->setDistance(distance);
    chemin->calcAttenuation(tabEtapes, *_pSolverAtmos, dp, hs, hr, Gs, Gm, Gr);
 
    // The minimal extension condition is evaluated reflection by reflection on
    // the broken reflected path.
    for (size_t i = 0; i < reflectionOrder; ++i)
    {
        const OPoint3D& S = pathPoints[i];
        const OPoint3D& O = pathPoints[i + 1];
        const OPoint3D& R = pathPoints[i + 2];
 
        const double a = barrierCombination[i]->pFaceGeomData->a;
        const double h = barrierCombination[i]->pFaceGeomData->h;
        const OVector3D& normal = barrierCombination[i]->pFaceGeomData->n;
 
        chemin->calcMinimalExtensionCondition(S, O, R, a, h, normal);
    }
 
    TabChemins.push_back(*chemin);
    return true;
}
 
TYAcousticModel9613Solver2024::ReflectingSegmentCache
TYAcousticModel9613Solver2024::makeReflectingSegmentCache(const TYSIntersection& inter) const
{
    ReflectingSegmentCache cache;
    cache.barrier = &inter;
 
    const OSegment3D& seg = inter.segInter[1];
 
    cache.ax = seg._ptA._x;
    cache.ay = seg._ptA._y;
    cache.az = seg._ptA._z;
 
    cache.ux = seg._ptB._x - seg._ptA._x;
    cache.uy = seg._ptB._y - seg._ptA._y;
    cache.uz = seg._ptB._z - seg._ptA._z;
 
    // The cache stores the support line of the reflecting segment.
    // invNorm2 is used to project points on this line without recomputing the
    // segment norm at each recursive step.
    const double norm2 = cache.ux * cache.ux + cache.uy * cache.uy + cache.uz * cache.uz;
    if (norm2 > TYSEUILCONFONDUS * TYSEUILCONFONDUS)
    {
        cache.invNorm2 = 1.0 / norm2;
    }
    else
    {
        cache.invNorm2 = 0.0;
    }
 
    return cache;
}
 
void TYAcousticModel9613Solver2024::reflectPointAboutCachedLine(const ReflectingSegmentCache& cache,
                                                                const OPoint3D& inputPoint,
                                                                OPoint3D& reflectedPoint) const
{
    // Keep legacy behavior for degenerate cases.
    if ((cache.barrier == nullptr) || (cache.invNorm2 <= 0.0))
    {
        if (cache.barrier != nullptr)
        {
            cache.barrier->segInter[1].symetrieOf(inputPoint, reflectedPoint);
        }
        return;
    }
 
    const double apx = inputPoint._x - cache.ax;
    const double apy = inputPoint._y - cache.ay;
    const double apz = inputPoint._z - cache.az;
 
    const double t = (apx * cache.ux + apy * cache.uy + apz * cache.uz) * cache.invNorm2;
 
    const double hx = cache.ax + t * cache.ux;
    const double hy = cache.ay + t * cache.uy;
    const double hz = cache.az + t * cache.uz;
 
    // Reflection is performed about the support line of the segment in plane EL.
    reflectedPoint._x = 2.0 * hx - inputPoint._x;
    reflectedPoint._y = 2.0 * hy - inputPoint._y;
    reflectedPoint._z = 2.0 * hz - inputPoint._z;
}
 
void TYAcousticModel9613Solver2024::buildReflectionPathsStreaming(
    size_t reflectionOrder, const std::vector<ReflectingSegmentCache>& reflectingSegments,
    std::vector<const TYSIntersection*>& currentCombination, std::deque<OPoint3D>& currentImageSources,
    const std::deque<TYSIntersection>& tabIntersect, const OSegment3D& ray,
    const tympan::AcousticSource& source, TYTabChemin9613Solver& TabChemins, double distance) const
{
    if (currentCombination.size() == reflectionOrder)
    {
        return;
    }
 
    for (const auto& segment : reflectingSegments)
    {
        // Only consecutive reuse of the same reflecting segment is forbidden.
        if (!currentCombination.empty() && currentCombination.back() == segment.barrier)
        {
            continue;
        }
 
        const OPoint3D& previousImageSource = currentImageSources.back();
 
        OPoint3D nextImageSource;
        reflectPointAboutCachedLine(segment, previousImageSource, nextImageSource);
 
        currentCombination.push_back(segment.barrier);
        currentImageSources.push_back(nextImageSource);
 
        // currentImageSources always stores:
        //   S0 = source point, then one image source per current reflection.
        buildReflectionPathForCombination(currentCombination, currentImageSources, tabIntersect, ray, source,
                                          TabChemins, distance);
 
        buildReflectionPathsStreaming(reflectionOrder, reflectingSegments, currentCombination,
                                      currentImageSources, tabIntersect, ray, source, TabChemins, distance);
 
        currentImageSources.pop_back();
        currentCombination.pop_back();
    }
}
 
bool TYAcousticModel9613Solver2024::buildReflectionPathForCombination(
    const std::vector<const TYSIntersection*>& currentCombination,
    const std::deque<OPoint3D>& imageSourcesList, const std::deque<TYSIntersection>& tabIntersect,
    const OSegment3D& ray, const tympan::AcousticSource& source, TYTabChemin9613Solver& TabChemins,
    double distance) const
{
    std::deque<OPoint3D> reflectionPointsList;
 
    // Geometry is built in 3 stages:
    //   1. reconstruct reflection points from current image sources,
    //   2. validate the broken path against the scene,
    //   3. build the acoustic path only for valid combinations.
    if (!buildReflectionPoints(currentCombination, imageSourcesList, ray._ptB, reflectionPointsList))
    {
        return false;
    }
 
    if (!validateReflectionCandidate(currentCombination, reflectionPointsList, tabIntersect, ray._ptA,
                                     ray._ptB))
    {
        return false;
    }
 
    return buildReflectionPath(currentCombination, imageSourcesList, reflectionPointsList, ray, source,
                               TabChemins, distance);
}