river_tools.h

Go to the documentation of this file.

 
#ifndef UG4_RIVER_TOOLS_H
#define UG4_RIVER_TOOLS_H
 
#include <vector>
#include <queue>
#include <map>
#include <set>
#include "common/common.h"
#include "lib_grid/algorithms/trees/kd_tree_static.h"
#include "lib_grid/file_io/file_io_ugx.h"
 
namespace ug {
namespace d3f {
 
template <typename TDomain>
void project_rivers_on_terrain(SmartPtr<TDomain> river_domain, SmartPtr<TDomain> soil_domain, const char* top_subset) {
    typedef KDTreeStatic<typename TDomain::position_attachment_type, TDomain::dim, typename TDomain::position_type> tree_type;
    typedef typename TDomain::grid_type TGrid;
    static const int dim = TDomain::dim;
 
    auto river_posAcc = river_domain->position_accessor();
    auto soil_posAcc = soil_domain->position_accessor();
    auto soil_ssHandler = soil_domain->subset_handler();
    auto river_grid = river_domain->grid();
    auto soil_grid = soil_domain->grid();
    
    const int top_subset_index = soil_ssHandler->get_subset_index(top_subset);
    UG_ASSERT(top_subset_index > -1, "project_rivers_on_terrain: Top subset '" << top_subset << "' not found.");
 
    UG_LOG("Projecting river network onto terrain subset: " << top_subset << "\n");
 
    // 1. Build a spatial index for the terrain surface to speed up vertex lookup
    auto tree = tree_type();
    tree.create_from_grid(*soil_grid,
                          soil_ssHandler->template begin<Vertex>(top_subset_index, soil_grid->top_level()),
                          soil_ssHandler->template end<Vertex>(top_subset_index, soil_grid->top_level()),
                          soil_posAcc,
                          15, 20);
 
    MathVector<dim> vDir = MathVector<dim>(0.0);
    vDir[dim - 1] = 1.0; // Raycast direction is +z (assuming terrain is above river nodes)
 
    std::set<Vertex*> skipped_vertices;
    
    // 2. Initial projection via raycasting
    for (size_t lvl = 0; lvl <= river_grid->top_level(); ++lvl) {
        auto vIt = river_grid->template begin<Vertex>(lvl);
        auto end = river_grid->template end<Vertex>(lvl);
 
        for (; vIt != end; ++vIt) {
            Vertex* vRiver = *vIt;
            std::vector<Vertex*> neighborhood;
            tree.get_neighbourhood(neighborhood, river_posAcc[vRiver], 1);
            
            if (neighborhood.empty()) {
                skipped_vertices.insert(vRiver);
                continue;
            }
 
            Vertex* closest_soil_vertex = neighborhood[0];
            typename TGrid::template traits<Face>::secure_container faces;
            soil_grid->template associated_elements<Vertex>(faces, closest_soil_vertex);
 
            bool hit = false;
            for (size_t i = 0; i < faces.size(); ++i) {
                if (soil_ssHandler->get_subset_index(faces[i]) != top_subset_index) continue;
 
                number smin, smax;
                if (RayElementIntersection(smin, smax, river_posAcc[vRiver], vDir, faces[i], 
                                           *soil_grid, soil_posAcc, SMALL)) {
                    hit = true;
                    river_posAcc[vRiver][dim - 1] = PointOnRay(river_posAcc[vRiver], vDir, smin)[dim - 1];
                    break;
                }
            }
 
            if (!hit) {
                skipped_vertices.insert(vRiver);
            }
        }
    }
 
    // 3. Handle skipped vertices by inheriting elevation from the nearest projected neighbor
    if (!skipped_vertices.empty()) {
        UG_LOG("Handling " << skipped_vertices.size() << " vertices outside terrain extent...\n");
        
        // We use a copy of the set to iterate, as we will remove elements from the original
        std::vector<Vertex*> to_process(skipped_vertices.begin(), skipped_vertices.end());
        
        for (Vertex* v : to_process) {
            if (skipped_vertices.find(v) == skipped_vertices.end()) continue;
 
            std::queue<Vertex*> q;
            std::set<Vertex*> visited;
            q.push(v);
            visited.insert(v);
 
            bool success = false;
            while (!q.empty()) {
                Vertex* curr = q.front();
                q.pop();
 
                // If curr is not in skipped_vertices, it means it has a valid elevation
                if (skipped_vertices.find(curr) == skipped_vertices.end()) {
                    river_posAcc[v][dim - 1] = river_posAcc[curr][dim - 1];
                    skipped_vertices.erase(v);
                    success = true;
                    break;
                }
 
                typename TGrid::template traits<Edge>::secure_container edges;
                river_grid->template associated_elements<Vertex>(edges, curr);
                for (size_t i = 0; i < edges.size(); ++i) {
                    for (int j = 0; j < 2; ++j) {
                        Vertex* nb = edges[i]->vertex(j);
                        if (visited.find(nb) == visited.end()) {
                            visited.insert(nb);
                            q.push(nb);
                        }
                    }
                }
            }
            if (!success) {
                UG_THROW("Could not find a projected neighbor for vertex at " << river_posAcc[v]);
            }
        }
    }
 
    // 4. Debug output: Save the projected grid
    std::stringstream ss;
    ss << "projected_river";
    #ifdef UG_PARALLEL
    ss << "_p" << pcl::ProcRank();
    #endif
    ss << ".ugx";
    
    UG_LOG("Saving projected river grid to: " << ss.str() << "\n");
    SaveGridToUGX(*river_grid, *(river_domain->subset_handler()), ss.str().c_str());
}
 
} // namespace d3f
} // namespace ug
 
#endif // UG4_RIVER_TOOLS_H