TYReflectionPathFinder.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 "TYReflectionPathFinder.h"
#include "Tympan/geometric_methods/AcousticRaytracer/Geometry/UniformCircularSampler.h"
 
TYReflectionPathFinderART::TYReflectionPathFinderART(unsigned int nbRays, double rcptRadius)
    : TYAbstractReflectionPathFinder(), _nbRays(nbRays), _rcptRadius(rcptRadius)
{
    _extrusionLength = 10.; // hard-coded value, should not have any impact on results
 
    // AcousticRaytracer configuration
    constexpr double maxDouble{std::numeric_limits<double>::max()};
    _ARTSimulation.getConfiguration()->Accelerator = 2;        // BVH data-structure
    _ARTSimulation.getConfiguration()->MaxTreeDepth = 12;      // BVH parameter
    _ARTSimulation.getConfiguration()->MaxDiffraction = 0;     // should disable diffractions
    _ARTSimulation.getConfiguration()->MaxLength = maxDouble;  // disable length selector
    _ARTSimulation.getConfiguration()->UseSol = false;         // disable reflections on ground triangles
    _ARTSimulation.getConfiguration()->UsePostFilters = false; // disable non-default selectors
    _ARTSimulation.getConfiguration()->KeepDebugRay = false;   // do not store invalidated rays
}
 
std::unique_ptr<TYAbstractReflectionPathFinder> TYReflectionPathFinderART::instanciate() const
{
    return std::make_unique<TYReflectionPathFinderART>(this->_nbRays, this->_rcptRadius);
}
 
bool TYReflectionPathFinderART::setup(size_t maxOrder, const std::deque<TYSIntersection>& tabIntersect,
                                      const OPoint3D& source, const OPoint3D& receptor)
{
    setupConfiguration(maxOrder);
    if (!setupScene(tabIntersect, source, receptor))
    {
        return false;
    }
    setupAcceleratingDataStructure();
    setupSolver();
    computeRaytracing();
    buildReflectionPaths(source, receptor);
 
    return true;
}
 
bool TYReflectionPathFinderART::remainingPaths() const
{
    return _paths.size() > 0;
}
 
TYReflectionPath TYReflectionPathFinderART::popPath()
{
    if (_paths.size() > 0)
    {
        TYReflectionPath res = _paths.back();
        _paths.pop_back();
        return res;
    }
    else
    {
        return TYReflectionPath();
    }
}
 
void TYReflectionPathFinderART::setupConfiguration(size_t maxOrder)
{
    _ARTSimulation.getConfiguration()->MaxProfondeur = maxOrder;
    _ARTSimulation.getConfiguration()->MaxReflexion = maxOrder;
}
 
bool TYReflectionPathFinderART::setupScene(const std::deque<TYSIntersection>& tabIntersect,
                                           const OPoint3D& sourcePoint, const OPoint3D& receptorPoint)
{
    _ELnormal = computeELPlaneNormal(sourcePoint, receptorPoint);
 
    // add faces
    for (const TYSIntersection& inter : tabIntersect)
    {
        if (!processIntersectingSegment(inter))
        {
            return false;
        }
    }
 
    // create initial rays sampler
    UniformCircularSampler sampler{_nbRays, OVector3Dtovec3(_ELnormal)};
 
    // add source, using the sampler
    Source source;
    source.setPosition(OPoint3Dtovec3(sourcePoint));
    source.setSampler(&sampler); // no transfer of ownership
                                 // it may segfault further once sampler is popped from the stack
    source.setInitialRayCount(_nbRays);
    _ARTSimulation.addSource(source); // a copy of the sampler is in fact hidden here using
                                      // `UniformCircularSampler(UniformCircularSampler*)`
                                      // in practice it does not segfault (I assume that only the copy is
                                      // used when setupScene is finished)
 
    // add receptor
    Recepteur receptor;
    receptor.setPosition(OPoint3Dtovec3(receptorPoint));
    receptor.setRadius(_rcptRadius);
    _ARTSimulation.addRecepteur(receptor);
 
    return true;
}
 
void TYReflectionPathFinderART::setupAcceleratingDataStructure()
{
    const unsigned int acceleratorID{_ARTSimulation.getConfiguration()->Accelerator};
 
    // this instruction builds the accelerating data-structure of the scene
    // the second argument sets the behavior when intersections between a ray and scene elements are found:
    // here only the first intersection is considered
    _ARTSimulation.getScene()->finish(acceleratorID, leafTreatment::FIRST);
 
    // this instruction builds the accelerating data-structure of receptors
    // the second argument sets the behavior when intersections between a ray and scene elements are found:
    // here all intersections considered
    _ARTSimulation.get_receptors_landscape()->finish(acceleratorID, leafTreatment::ALL);
}
 
void TYReflectionPathFinderART::setupSolver()
{
    Scene& scene = *_ARTSimulation.getScene();
    _pSolver = std::make_unique<BasicSolver>();
    _ARTSimulation.setSolver(_pSolver.get()); // no transfer of ownership
 
    // this instruction builds selectors (only a CleanerSelector and a LengthSelector using default Simulation
    // configuration)
    _pSolver->postTreatmentScene(&scene, _ARTSimulation.getSources(), _ARTSimulation.getRecepteurs());
}
 
void TYReflectionPathFinderART::computeRaytracing()
{
    _ARTSimulation.launchSimulation();
}
 
void TYReflectionPathFinderART::buildReflectionPaths(const OPoint3D& source, const OPoint3D& receptor)
{
    for (Ray* pRay : *_pSolver->getValidRays())
    {
        // ignore direct path
        const bool isReflectionPath{pRay->getNbEvents() > 0};
        if (!isReflectionPath)
        {
            continue;
        }
 
        // build the path assuming it is valid (to minimize TYReflectionPath constructions)
        _paths.emplace_back();
        TYReflectionPath& path = _paths.back();
        path.source = source;
        path.receptor = receptor;
        for (boost::shared_ptr<Event> pEvent : *pRay->getEvents())
        {
            // it is assumed that each event is a reflection
            const OPoint3D reflectionPoint = vec3toOPoint3D(pEvent->getPosition());
            path.reflectionPoints.push_back(reflectionPoint);
            const TYSIntersection* pInterReflection = _shapeToTYSIntersection[pEvent->getShape()];
            path.reflectingSegments.push_back(pInterReflection);
        }
 
        // validate the path and remove it if not valid
        const bool validPath = isPathValid(path);
        if (validPath)
        {
            _foundSequences.insert(path.reflectingSegments);
        }
        else
        {
            _paths.pop_back();
        }
    }
 
    _ARTSimulation.clean();
}
 
bool TYReflectionPathFinderART::processIntersectingSegment(const TYSIntersection& inter)
{
    if (inter.bIntersect[1])
    {
        if (inter.isInfra)
        {
            return processInfraIntersectingSegment(inter);
        }
        else
        {
            processTopoIntersectingSegment(inter);
            return true;
        }
    }
    else
    {
        // do nothing with intersecting segment which do not intersect EL-place
        return true;
    }
}
 
bool TYReflectionPathFinderART::processInfraIntersectingSegment(const TYSIntersection& inter)
{
    if (!inter.pFaceGeomData)
    {
        // should never happen
        return false;
    }
 
    const OVector3D& normal{inter.pFaceGeomData->n};
    OPoint3D p1Top, p1Bottom, p2Top, p2Bottom;
    buildQuadranglePointsFromInter(inter, p1Top, p1Bottom, p2Top, p2Bottom);
    addOrientedQuadrangle(p1Top, p1Bottom, p2Top, p2Bottom, normal, inter);
 
    return true;
}
 
bool TYReflectionPathFinderART::processTopoIntersectingSegment(const TYSIntersection& inter)
{
    const OSegment3D& segment{inter.segInter[1]};
 
    const OPoint3D p1{segment._ptA};
    const OPoint3D p2{segment._ptB};
 
    // build a normal of the segment in EL-plane
    // its orientation has no importance
    const OVector3D tangent{p1, p2};
    OVector3D normal{-tangent._y, tangent._x, 0.};
    normal.normalize();
 
    OPoint3D p1Top, p1Bottom, p2Top, p2Bottom;
    buildQuadranglePointsFromInter(inter, p1Top, p1Bottom, p2Top, p2Bottom);
 
    addOrientedQuadrangle(p1Top + normal * EPSILON_6, p1Bottom + normal * EPSILON_6,
                          p2Top + normal * EPSILON_6, p2Bottom + normal * EPSILON_6, normal, inter);
    addOrientedQuadrangle(p1Top + normal * -EPSILON_6, p1Bottom + normal * -EPSILON_6,
                          p2Top + normal * -EPSILON_6, p2Bottom + normal * -EPSILON_6, normal * -1., inter);
 
    return true;
}
 
void TYReflectionPathFinderART::buildQuadranglePointsFromInter(const TYSIntersection& inter, OPoint3D& p1Top,
                                                               OPoint3D& p1Bottom, OPoint3D& p2Top,
                                                               OPoint3D& p2Bottom) const
{
    const OSegment3D& segment{inter.segInter[1]};
    const OPoint3D& p1{segment._ptA};
    const OPoint3D& p2{segment._ptB};
    const OVector3D upExtrusion = _ELnormal * _extrusionLength;
    const OVector3D downExtrusion = _ELnormal * -_extrusionLength;
    p1Top = p1 + upExtrusion;
    p1Bottom = p1 + downExtrusion;
    p2Top = p2 + upExtrusion;
    p2Bottom = p2 + downExtrusion;
}
 
void TYReflectionPathFinderART::addOrientedQuadrangle(const OPoint3D& p1Top, const OPoint3D& p1Bottom,
                                                      const OPoint3D& p2Top, const OPoint3D& p2Bottom,
                                                      const OVector3D vExt,
                                                      const TYSIntersection& sourceInter)
{
    Shape* pT1 = addOrientedTriangle(p1Bottom, p2Top, p1Top, vExt);
    Shape* pT2 = addOrientedTriangle(p1Bottom, p2Bottom, p2Top, vExt);
    _shapeToTYSIntersection[pT1] = &sourceInter;
    _shapeToTYSIntersection[pT2] = &sourceInter;
}
 
Shape* TYReflectionPathFinderART::addOrientedTriangle(const OPoint3D& p1, const OPoint3D& p2,
                                                      const OPoint3D& p3, const OVector3D& vExt)
{
    Scene& scene = *_ARTSimulation.getScene();
    unsigned int i1, i2, i3;
 
    const OVector3D p1p2{p1, p2};
    const OVector3D p1p3{p1, p3};
    const OVector3D nCandidate = p1p2.cross(p1p3);
 
    if (nCandidate.scalar(vExt) > 0)
    {
        // v1 = p1, v2 = p2, v3 = p3
        scene.addVertex(OPoint3Dtovec3(p1), i1);
        scene.addVertex(OPoint3Dtovec3(p2), i2);
        scene.addVertex(OPoint3Dtovec3(p3), i3);
    }
    else
    {
        // v1 = p1, v2 = p3, v3 = p2
        scene.addVertex(OPoint3Dtovec3(p1), i1);
        scene.addVertex(OPoint3Dtovec3(p3), i2);
        scene.addVertex(OPoint3Dtovec3(p2), i3);
    }
 
    return scene.addTriangle(i1, i2, i3, nullptr, false);
}
 
OVector3D TYReflectionPathFinderART::computeELPlaneNormal(const OPoint3D& source, const OPoint3D& receptor)
{
    const OVector3D e3{0., 0., 1.};
    const OVector3D v{source, receptor};
    const OVector3D u = v.cross(e3);
    const OVector3D n = u.cross(v); // normal vector, so that (u, v, n) is a right-handed basis
    const OVector3D nNormalized = n * (1. / n.norme());
    return nNormalized;
}
 
bool TYReflectionPathFinderART::isPathValid(const TYReflectionPath& path)
{
    if (isPathAlreadyFound(path))
    {
        return false;
    }
    if (isLastSegmentIntersectingSceneElement(path))
    {
        return false;
    }
    return true;
}
 
bool TYReflectionPathFinderART::isPathAlreadyFound(const TYReflectionPath& path) const
{
    return _foundSequences.count(path.reflectingSegments) > 0;
}
 
bool TYReflectionPathFinderART::isLastSegmentIntersectingSceneElement(const TYReflectionPath& path)
{
    const OPoint3D& lastReflection = path.reflectionPoints.back();
    const OPoint3D& receptor = path.receptor;
 
    const double distLastReflectionReceptor = lastReflection.distFrom(receptor);
 
    const OVector3D rayDirection{lastReflection, receptor};
    Ray ray{OPoint3Dtovec3(lastReflection), OPoint3Dtovec3(rayDirection)};
    std::list<Intersection> intersections;
    double distToClosestElement =
        static_cast<double>(_ARTSimulation.getScene()->getAccelerator()->traverse(&ray, intersections));
 
    return distToClosestElement > 0. && distToClosestElement < distLastReflectionReceptor;
}