From 70798e1c63fd8c72bacbfa05bdeb992b2533fa89 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 14:52:29 +0300 Subject: [PATCH 01/17] visualizer: adopt backend-neutral cScene3DNode (OSG+VSG coexistence) Companion to the omnetpp-dev change that makes cScene3DNode a renderer-neutral handle rather than an alias of osg::Node. Wrap OSG scene roots with omnetpp::createScene3DNode() when calling cOsgCanvas::setScene(), and recover the typed root via getOsgRoot()/getVsgRoot() instead of casting cScene3DNode directly to a scene-graph type. Updated both visualizer trees: - OSG: SceneOsgVisualizer, SceneOsgEarthVisualizer, SceneOsgVisualizerBase, OsgScene, and GeographicCoordinateSystem (osgEarth map node lookup) - VSG: SceneVsgVisualizer, VsgScene Lets INET build and run against an OMNeT++ that contains both 3-D backends, with the renderer selected per simulation. Co-Authored-By: Claude Opus 4.8 --- .../common/geometry/common/GeographicCoordinateSystem.cc | 3 ++- src/inet/visualizer/osg/base/SceneOsgVisualizerBase.cc | 4 +++- src/inet/visualizer/osg/scene/SceneOsgEarthVisualizer.cc | 3 ++- src/inet/visualizer/osg/scene/SceneOsgVisualizer.cc | 4 +++- src/inet/visualizer/osg/util/OsgScene.cc | 6 ++++-- src/inet/visualizer/vsg/scene/SceneVsgVisualizer.cc | 2 +- src/inet/visualizer/vsg/util/VsgScene.cc | 2 +- 7 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/inet/common/geometry/common/GeographicCoordinateSystem.cc b/src/inet/common/geometry/common/GeographicCoordinateSystem.cc index a251d12ea38..8233f763735 100644 --- a/src/inet/common/geometry/common/GeographicCoordinateSystem.cc +++ b/src/inet/common/geometry/common/GeographicCoordinateSystem.cc @@ -10,6 +10,7 @@ #if defined(WITH_OSGEARTH) && defined(INET_WITH_VISUALIZATIONOSG) #include #include +#include "qtenv/osg/osgscenehandle.h" // omnetpp::getOsgRoot() #endif namespace inet { @@ -49,7 +50,7 @@ Define_Module(OsgGeographicCoordinateSystem); void OsgGeographicCoordinateSystem::initialize(int stage) { if (stage == INITSTAGE_LOCAL) { - auto mapScene = getParentModule()->getOsgCanvas()->getScene(); + auto mapScene = omnetpp::getOsgRoot(getParentModule()->getOsgCanvas()->getScene()); mapNode = osgEarth::MapNode::findMapNode(mapScene); if (mapNode == nullptr) throw cRuntimeError("Count not find map node in the scene"); diff --git a/src/inet/visualizer/osg/base/SceneOsgVisualizerBase.cc b/src/inet/visualizer/osg/base/SceneOsgVisualizerBase.cc index e0567c18289..8be1bb28d94 100644 --- a/src/inet/visualizer/osg/base/SceneOsgVisualizerBase.cc +++ b/src/inet/visualizer/osg/base/SceneOsgVisualizerBase.cc @@ -12,6 +12,8 @@ #include #include +#include "qtenv/osg/osgscenehandle.h" // omnetpp::createScene3DNode(osg::Node*) + #include "inet/common/ModuleAccess.h" #include "inet/visualizer/osg/scene/NetworkNodeOsgVisualizer.h" #include "inet/visualizer/osg/util/OsgScene.h" @@ -28,7 +30,7 @@ void SceneOsgVisualizerBase::initializeScene() throw cRuntimeError("OSG canvas scene at '%s' has been already initialized", visualizationTargetModule->getFullPath().c_str()); else { auto topLevelScene = new inet::osg::TopLevelScene(); - osgCanvas->setScene(topLevelScene); + osgCanvas->setScene(omnetpp::createScene3DNode(topLevelScene)); const char *clearColor = par("clearColor"); if (*clearColor != '\0') osgCanvas->setClearColor(cFigure::Color(clearColor)); diff --git a/src/inet/visualizer/osg/scene/SceneOsgEarthVisualizer.cc b/src/inet/visualizer/osg/scene/SceneOsgEarthVisualizer.cc index 0bf764846c8..3c84bed9b5e 100644 --- a/src/inet/visualizer/osg/scene/SceneOsgEarthVisualizer.cc +++ b/src/inet/visualizer/osg/scene/SceneOsgEarthVisualizer.cc @@ -19,6 +19,7 @@ #include #include #include +#include "qtenv/osg/osgscenehandle.h" // omnetpp::getOsgRoot() #endif // ifdef WITH_OSGEARTH namespace inet { @@ -62,7 +63,7 @@ void SceneOsgEarthVisualizer::initializeScene() throw cRuntimeError("Could not read earth map file '%s'", mapFileString.c_str()); auto osgCanvas = visualizationTargetModule->getOsgCanvas(); osgCanvas->setViewerStyle(cOsgCanvas::STYLE_EARTH); - auto topLevelScene = check_and_cast(osgCanvas->getScene()); + auto topLevelScene = check_and_cast(omnetpp::getOsgRoot(osgCanvas->getScene())); topLevelScene->addChild(mapScene); mapNode = MapNode::findMapNode(mapScene); if (mapNode == nullptr) diff --git a/src/inet/visualizer/osg/scene/SceneOsgVisualizer.cc b/src/inet/visualizer/osg/scene/SceneOsgVisualizer.cc index 89f03d7aec6..c7788dac5d1 100644 --- a/src/inet/visualizer/osg/scene/SceneOsgVisualizer.cc +++ b/src/inet/visualizer/osg/scene/SceneOsgVisualizer.cc @@ -10,6 +10,8 @@ #include #include +#include "qtenv/osg/osgscenehandle.h" // omnetpp::getOsgRoot() + #include "inet/common/ModuleAccess.h" #include "inet/visualizer/osg/util/OsgScene.h" #include "inet/visualizer/osg/util/OsgUtils.h" @@ -39,7 +41,7 @@ void SceneOsgVisualizer::initialize(int stage) void SceneOsgVisualizer::initializeScene() { SceneOsgVisualizerBase::initializeScene(); - auto topLevelScene = check_and_cast(visualizationTargetModule->getOsgCanvas()->getScene()); + auto topLevelScene = check_and_cast(omnetpp::getOsgRoot(visualizationTargetModule->getOsgCanvas()->getScene())); topLevelScene->addChild(new inet::osg::SimulationScene()); } diff --git a/src/inet/visualizer/osg/util/OsgScene.cc b/src/inet/visualizer/osg/util/OsgScene.cc index c267057fb73..1b4bb1e23e7 100644 --- a/src/inet/visualizer/osg/util/OsgScene.cc +++ b/src/inet/visualizer/osg/util/OsgScene.cc @@ -7,6 +7,8 @@ #include "inet/visualizer/osg/util/OsgScene.h" +#include "qtenv/osg/osgscenehandle.h" // omnetpp::createScene3DNode(osg::Node*), getOsgRoot() + namespace inet { namespace osg { @@ -35,7 +37,7 @@ SimulationScene *TopLevelScene::getSimulationScene() SimulationScene *TopLevelScene::getSimulationScene(cModule *module) { auto osgCanvas = module->getOsgCanvas(); - auto topLevelScene = dynamic_cast(osgCanvas->getScene()); + auto topLevelScene = dynamic_cast(omnetpp::getOsgRoot(osgCanvas->getScene())); if (topLevelScene != nullptr) return topLevelScene->getSimulationScene(); else { @@ -43,7 +45,7 @@ SimulationScene *TopLevelScene::getSimulationScene(cModule *module) topLevelScene = new TopLevelScene(); topLevelScene->addChild(simulationScene); // NOTE: these are the default values when there's no SceneOsgVisualizer - osgCanvas->setScene(topLevelScene); + osgCanvas->setScene(omnetpp::createScene3DNode(topLevelScene)); osgCanvas->setClearColor(cFigure::Color("#FFFFFF")); osgCanvas->setZNear(0.1); osgCanvas->setZFar(100000); diff --git a/src/inet/visualizer/vsg/scene/SceneVsgVisualizer.cc b/src/inet/visualizer/vsg/scene/SceneVsgVisualizer.cc index f1dc4d3eaa0..f1391e52f90 100644 --- a/src/inet/visualizer/vsg/scene/SceneVsgVisualizer.cc +++ b/src/inet/visualizer/vsg/scene/SceneVsgVisualizer.cc @@ -39,7 +39,7 @@ void SceneVsgVisualizer::initializeScene() { SceneVsgVisualizerBase::initializeScene(); auto sceneNode = visualizationTargetModule->getOsgCanvas()->getScene(); - auto topLevelScene = sceneNode != nullptr ? dynamic_cast(sceneNode->getRoot().get()) : nullptr; + auto topLevelScene = dynamic_cast(omnetpp::getVsgRoot(sceneNode).get()); if (topLevelScene == nullptr) throw cRuntimeError("Cannot find the VSG top level scene"); topLevelScene->addChild(inet::vsg::SimulationScene::create()); diff --git a/src/inet/visualizer/vsg/util/VsgScene.cc b/src/inet/visualizer/vsg/util/VsgScene.cc index 91378cceaea..8dd23db92f3 100644 --- a/src/inet/visualizer/vsg/util/VsgScene.cc +++ b/src/inet/visualizer/vsg/util/VsgScene.cc @@ -31,7 +31,7 @@ SimulationScene *TopLevelScene::getSimulationScene(cModule *module) // Under WITH_VSG the cOsgCanvas scene is an opaque omnetpp::cScene3DNode that // wraps the VSG scene-root group; recover our TopLevelScene from its root. auto sceneNode = osgCanvas->getScene(); - ::vsg::ref_ptr root = sceneNode != nullptr ? sceneNode->getRoot() : ::vsg::ref_ptr(); + ::vsg::ref_ptr root = omnetpp::getVsgRoot(sceneNode); // null / wrong-backend safe auto topLevelScene = root ? dynamic_cast(root.get()) : nullptr; if (topLevelScene != nullptr) return topLevelScene->getSimulationScene(); From ac141dd081d3238ded28bea3582f8bd545e19e05 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 16:24:42 +0300 Subject: [PATCH 02/17] =?UTF-8?q?docs:=20vsg=20README=20=E2=80=94=20OSG=20?= =?UTF-8?q?and=20VSG=20can=20coexist;=20how=20to=20build=20with=20both?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Corrects the "VisualizationOsg and VisualizationVsg are mutually exclusive" note (true only before the backend-neutral cScene3DNode change) and documents that an OMNeT++ configured with both WITH_OSG=yes and WITH_VSG=yes contains both 3D backends, selected per simulation. Co-Authored-By: Claude Opus 4.8 --- src/inet/visualizer/vsg/README.md | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/inet/visualizer/vsg/README.md b/src/inet/visualizer/vsg/README.md index d2e15d85043..b1f11fc7c8a 100644 --- a/src/inet/visualizer/vsg/README.md +++ b/src/inet/visualizer/vsg/README.md @@ -387,16 +387,23 @@ switching to the 3D view, and pressing Run. ## Building & running - Requires a VSG-enabled OMNeT++ build (in this environment: the - `omnetpp-dev` checkout, configured/built with VSG support, - i.e. `WITH_VSG=yes`). -- INET side: the `VisualizationVsg` feature (`.oppfeatures`, id - `VisualizationVsg`) must be enabled; it compiles with - `-DINET_WITH_VISUALIZATIONVSG`, requires the `PhysicalEnvironment` and - `VisualizationCommon` features, and is **mutually exclusive** with - `VisualizationOsg` — an OMNeT++ build has either OSG or VSG support, never - both, so the two visualizer features cannot be enabled simultaneously. + `omnetpp-dev` checkout, configured with `WITH_VSG=yes`). +- **OSG and VSG can now coexist in one build.** Since OMNeT++'s `cScene3DNode` + became a backend-neutral handle (omnetpp-dev `topic/VSG`, the + "backend-neutral `cScene3DNode`" change), a single OMNeT++ build can contain + **both** 3D backends and choose one **per simulation**. To get this, + configure OMNeT++ with **both** `WITH_OSG=yes` *and* `WITH_VSG=yes`: the + Qtenv loader then builds and loads both `liboppqtenv-osg` and + `liboppqtenv-vsg` and selects the matching one per `cOsgCanvas` (via + `cScene3DNode::getBackendType()`). Building with only one backend still works. +- INET side: enable the `VisualizationVsg` feature (`.oppfeatures`; compiles + with `-DINET_WITH_VISUALIZATIONVSG`; requires `PhysicalEnvironment` and + `VisualizationCommon`). `VisualizationOsg` **may be enabled at the same + time** — a scene renders on whichever backend its visualizer creates. (An + earlier revision of this doc said the two were mutually exclusive; that was + true only before the backend-neutral `cScene3DNode` change.) - Build with `make MODE=release` (from the INET root, after - `opp_featuretool` / configure has enabled `VisualizationVsg`). + `opp_featuretool` / configure has enabled the desired feature(s)). - Run an example with `inet -u Qtenv` from its directory (e.g. `examples/visualizer/vsgsignalwave3d`), then switch to the 3D view in the Qtenv toolbar and press Run. From f3bd13cbe0eabafabe29c871ac3021c2694588a4 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 16:44:56 +0300 Subject: [PATCH 03/17] .oppfeatures: VisualizationVsg no longer mutually exclusive with OSG The feature description said VisualizationOsg and VisualizationVsg are mutually exclusive; that is no longer true after the backend-neutral cScene3DNode change. There is no `conflicts` between them, so both can be enabled together in an OMNeT++ built with WITH_OSG + WITH_VSG, each scene rendering on the backend its visualizer creates. (VisualizationVsg is still initiallyEnabled=false, so enable it explicitly for coexistence.) Co-Authored-By: Claude Opus 4.8 --- .oppfeatures | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/.oppfeatures b/.oppfeatures index b28a9f38bb9..d1f040b2189 100644 --- a/.oppfeatures +++ b/.oppfeatures @@ -1759,10 +1759,12 @@ name = "Visualization VSG (3D)" description = "Provides network-level 3D visualization features for physical layer, data link layer and network layer communications, and more. - This is the VulkanSceneGraph-based replacement for the OpenSceneGraph + This is the VulkanSceneGraph-based counterpart to the OpenSceneGraph (VisualizationOsg) 3D visualizer; it requires OMNeT++ built with VSG - support (WITH_VSG=yes). VisualizationOsg and VisualizationVsg are - mutually exclusive (OMNeT++ is built with either OSG or VSG, not both)." + support (WITH_VSG=yes). Since OMNeT++'s cScene3DNode became + backend-neutral, VisualizationOsg and VisualizationVsg may be enabled + together (in an OMNeT++ built with both WITH_OSG and WITH_VSG); each + scene renders on the backend its visualizer creates." initiallyEnabled = "false" requires = "PhysicalEnvironment VisualizationCommon" recommended = "" From 6d3244880e300ce1ad7713708340123433a5a536 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 17:13:05 +0300 Subject: [PATCH 04/17] =?UTF-8?q?visualizer:=20fix=20OSG-only=20build=20?= =?UTF-8?q?=E2=80=94=20add=20-I$(OMNETPP=5FROOT)/src=20for=20osgscenehandl?= =?UTF-8?q?e.h?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The five OSG source files adopted in this PR #include "qtenv/osg/osgscenehandle.h", which lives under $(OMNETPP_ROOT)/src. That include path was only added by the VisualizationVsg makefrag block (line 81), so an OSG-only build (VisualizationOsg enabled, VisualizationVsg disabled — the common existing configuration) could not find the header and failed to compile. Add -I$(OMNETPP_ROOT)/src to the VisualizationOsg block, mirroring the VSG block. Co-Authored-By: Claude Opus 4.8 --- src/makefrag | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/makefrag b/src/makefrag index b84ee35e483..d691d2f1eae 100644 --- a/src/makefrag +++ b/src/makefrag @@ -60,6 +60,11 @@ endif WITH_VISUALIZATIONOSG := $(shell (cd .. && $(FEATURETOOL) -q isenabled VisualizationOsg && echo enabled) ) ifeq ($(WITH_VISUALIZATIONOSG), enabled) ifeq ($(WITH_OSG), yes) + # -I$(OMNETPP_ROOT)/src gives access to the plugin's qtenv/osg/osgscenehandle.h + # (createScene3DNode / getOsgRoot / cScene3DNode), needed by the OSG scene + # visualizers now that cScene3DNode is backend-neutral. Mirrors the + # VisualizationVsg block below; required for OSG-only builds (VSG disabled). + CFLAGS += -I$(OMNETPP_ROOT)/src OMNETPP_LIBS += -losg -losgText -losgDB -losgGA -losgViewer -losgUtil -lOpenThreads # TODO: use $(OSG_LIBS) from Makefile.inc? Does that include -losgText? Why not? ifeq ($(WITH_OSGEARTH), yes) OMNETPP_LIBS += -losgEarth -losgEarthUtil # TODO: use $(OSGEARTH_LIBS) from Makefile.inc? -lgeos_c might also be needed. Is that included? Should it be? From e12ac52b71717e25ba630c1484722fbd0fff32d3 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 18:52:26 +0300 Subject: [PATCH 05/17] common: PLY point-cloud reader + 2.5D heightfield (LIDAR terrain groundwork) PlyPointCloudReader: the PLY parser extracted from the VSG scene visualizer (ascii + binary_little_endian, x/y/z + optional rgb located by property name, arbitrary scalar types/order, list properties skipped), returning retained point arrays + bounding box so non-visualizer code can consume LIDAR tiles. Heightfield: regular 2.5D grid rasterized from a point cloud with max-z-per- cell (DSM) semantics; O(1) bilinear elevation lookup, surface normals, elevation profiles along a segment (for future line-of-sight/diffraction models), bounded hole filling, isolated-spike clamping (stray LIDAR returns), peak/extent queries, and a maxCells guard. Order-independent rasterization keeps results deterministic across platforms. Both are pure C++ with no visualizer/VSG dependency; covered by tests/unit/{PlyPointCloudReader_1,Heightfield_1}.test. Co-Authored-By: Claude Fable 5 --- .../common/geometry/common/Heightfield.cc | 204 ++++++++++++++++++ src/inet/common/geometry/common/Heightfield.h | 91 ++++++++ .../geometry/common/PlyPointCloudReader.cc | 154 +++++++++++++ .../geometry/common/PlyPointCloudReader.h | 60 ++++++ tests/unit/Heightfield_1.test | 121 +++++++++++ tests/unit/PlyPointCloudReader_1.test | 114 ++++++++++ 6 files changed, 744 insertions(+) create mode 100644 src/inet/common/geometry/common/Heightfield.cc create mode 100644 src/inet/common/geometry/common/Heightfield.h create mode 100644 src/inet/common/geometry/common/PlyPointCloudReader.cc create mode 100644 src/inet/common/geometry/common/PlyPointCloudReader.h create mode 100644 tests/unit/Heightfield_1.test create mode 100644 tests/unit/PlyPointCloudReader_1.test diff --git a/src/inet/common/geometry/common/Heightfield.cc b/src/inet/common/geometry/common/Heightfield.cc new file mode 100644 index 00000000000..80c1da5b5fd --- /dev/null +++ b/src/inet/common/geometry/common/Heightfield.cc @@ -0,0 +1,204 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/common/geometry/common/Heightfield.h" + +#include +#include + +namespace inet { + +void Heightfield::buildFromPoints(const std::vector& xs, const std::vector& ys, const std::vector& zs, + double requestedCellSize, int64_t maxCells, int holeFillPasses, double despikeThreshold) +{ + size_t count = xs.size(); + if (count == 0 || ys.size() != count || zs.size() != count) + throw cRuntimeError("Heightfield: empty or inconsistent point arrays (%zu/%zu/%zu points)", xs.size(), ys.size(), zs.size()); + + double lo_x = xs[0], hi_x = xs[0], lo_y = ys[0], hi_y = ys[0]; + for (size_t i = 1; i < count; i++) { + lo_x = std::min(lo_x, xs[i]); hi_x = std::max(hi_x, xs[i]); + lo_y = std::min(lo_y, ys[i]); hi_y = std::max(hi_y, ys[i]); + } + double width = hi_x - lo_x, height = hi_y - lo_y; + if (width <= 0 && height <= 0) + throw cRuntimeError("Heightfield: degenerate point cloud (all points at the same x/y)"); + + // auto cell size: approximate mean point spacing over the covered area + double cs = requestedCellSize; + if (cs <= 0) + cs = std::sqrt(std::max(width, 1e-9) * std::max(height, 1e-9) / count); + int64_t nx = std::max((int64_t)1, (int64_t)std::ceil(width / cs)); + int64_t ny = std::max((int64_t)1, (int64_t)std::ceil(height / cs)); + if (nx * ny > maxCells) + throw cRuntimeError("Heightfield: grid of %lld x %lld cells exceeds the limit of %lld cells; use a larger cellSize (>= %g m) or raise maxCells", + (long long)nx, (long long)ny, (long long)maxCells, std::sqrt((double)nx * ny / maxCells) * cs); + + minX = lo_x; + minY = lo_y; + cellSize = cs; + numCellsX = (int)nx; + numCellsY = (int)ny; + cells.assign((size_t)(nx * ny), std::numeric_limits::quiet_NaN()); + + // rasterize: max z per cell (DSM semantics; order-independent, deterministic) + for (size_t i = 0; i < count; i++) { + int ix = std::min(numCellsX - 1, (int)((xs[i] - minX) / cellSize)); + int iy = std::min(numCellsY - 1, (int)((ys[i] - minY) / cellSize)); + float& cell = cells[iy * numCellsX + ix]; + float z = (float)zs[i]; + if (std::isnan(cell) || z > cell) + cell = z; + } + + // clamp isolated single-cell spikes (stray LIDAR returns) that exceed ALL + // non-NaN neighbors by more than the threshold; real structures survive + // because their edge cells have same-height neighbors + if (despikeThreshold > 0) { + std::vector> clamps; + for (int iy = 0; iy < numCellsY; iy++) { + for (int ix = 0; ix < numCellsX; ix++) { + float v = getCell(ix, iy); + if (std::isnan(v)) + continue; + float maxNeighbor = -std::numeric_limits::infinity(); + int n = 0; + for (int dy = -1; dy <= 1; dy++) { + for (int dx = -1; dx <= 1; dx++) { + if (dx == 0 && dy == 0) continue; + int jx = ix + dx, jy = iy + dy; + if (jx < 0 || jx >= numCellsX || jy < 0 || jy >= numCellsY) continue; + float w = getCell(jx, jy); + if (!std::isnan(w)) { maxNeighbor = std::max(maxNeighbor, w); n++; } + } + } + if (n >= 3 && v > maxNeighbor + despikeThreshold) + clamps.emplace_back((size_t)(iy * numCellsX + ix), maxNeighbor); + } + } + for (auto& c : clamps) + cells[c.first] = c.second; + } + + // bounded hole filling: NaN cells take the average of their non-NaN 8-neighbors + for (int pass = 0; pass < holeFillPasses; pass++) { + std::vector> fills; + for (int iy = 0; iy < numCellsY; iy++) { + for (int ix = 0; ix < numCellsX; ix++) { + if (!std::isnan(getCell(ix, iy))) + continue; + double sum = 0; + int n = 0; + for (int dy = -1; dy <= 1; dy++) { + for (int dx = -1; dx <= 1; dx++) { + if (dx == 0 && dy == 0) continue; + int jx = ix + dx, jy = iy + dy; + if (jx < 0 || jx >= numCellsX || jy < 0 || jy >= numCellsY) continue; + float v = getCell(jx, jy); + if (!std::isnan(v)) { sum += v; n++; } + } + } + if (n >= 3) // require some support to avoid smearing data into large empty regions + fills.emplace_back((size_t)(iy * numCellsX + ix), (float)(sum / n)); + } + } + if (fills.empty()) + break; + for (auto& f : fills) + cells[f.first] = f.second; + } +} + +double Heightfield::getElevation(double x, double y) const +{ + if (!isValid()) + return NaN; + if (x < minX || x > getMaxX() || y < minY || y > getMaxY()) + return NaN; + // bilinear between cell centers, clamped at the borders + double gx = (x - minX) / cellSize - 0.5; + double gy = (y - minY) / cellSize - 0.5; + int ix0 = (int)std::floor(gx), iy0 = (int)std::floor(gy); + double fx = gx - ix0, fy = gy - iy0; + ix0 = std::max(0, std::min(numCellsX - 1, ix0)); + iy0 = std::max(0, std::min(numCellsY - 1, iy0)); + int ix1 = std::min(numCellsX - 1, ix0 + 1); + int iy1 = std::min(numCellsY - 1, iy0 + 1); + float v00 = getCell(ix0, iy0), v10 = getCell(ix1, iy0); + float v01 = getCell(ix0, iy1), v11 = getCell(ix1, iy1); + if (std::isnan(v00) || std::isnan(v10) || std::isnan(v01) || std::isnan(v11)) { + // fall back to the nearest cell's value (may itself be NaN in an unfilled hole) + int ix = std::min(numCellsX - 1, (int)((x - minX) / cellSize)); + int iy = std::min(numCellsY - 1, (int)((y - minY) / cellSize)); + return getCell(ix, iy); + } + fx = std::max(0.0, std::min(1.0, fx)); + fy = std::max(0.0, std::min(1.0, fy)); + double v0 = v00 + (v10 - v00) * fx; + double v1 = v01 + (v11 - v01) * fx; + return v0 + (v1 - v0) * fy; +} + +Coord Heightfield::getNormal(double x, double y) const +{ + double h = cellSize; + double zxm = getElevation(x - h, y), zxp = getElevation(x + h, y); + double zym = getElevation(x, y - h), zyp = getElevation(x, y + h); + // fall back to one-sided differences at the borders + double zc = NaN; + if (std::isnan(zxm) || std::isnan(zxp) || std::isnan(zym) || std::isnan(zyp)) + zc = getElevation(x, y); + double dx, dy; + if (!std::isnan(zxm) && !std::isnan(zxp)) + dx = (zxp - zxm) / (2 * h); + else if (!std::isnan(zxp) && !std::isnan(zc)) + dx = (zxp - zc) / h; + else if (!std::isnan(zxm) && !std::isnan(zc)) + dx = (zc - zxm) / h; + else + return Coord(NaN, NaN, NaN); + if (!std::isnan(zym) && !std::isnan(zyp)) + dy = (zyp - zym) / (2 * h); + else if (!std::isnan(zyp) && !std::isnan(zc)) + dy = (zyp - zc) / h; + else if (!std::isnan(zym) && !std::isnan(zc)) + dy = (zc - zym) / h; + else + return Coord(NaN, NaN, NaN); + Coord normal(-dx, -dy, 1); + normal.normalize(); + return normal; +} + +Coord Heightfield::getPeakLocation() const +{ + int bestX = 0, bestY = 0; + float best = -std::numeric_limits::infinity(); + for (int iy = 0; iy < numCellsY; iy++) + for (int ix = 0; ix < numCellsX; ix++) { + float v = getCell(ix, iy); + if (!std::isnan(v) && v > best) { best = v; bestX = ix; bestY = iy; } + } + return Coord(minX + (bestX + 0.5) * cellSize, minY + (bestY + 0.5) * cellSize, best); +} + +std::vector Heightfield::computeProfile(const Coord& a, const Coord& b, double step) const +{ + if (step <= 0) + throw cRuntimeError("Heightfield::computeProfile: step must be positive (got %g)", step); + double dx = b.x - a.x, dy = b.y - a.y; + double distance = std::sqrt(dx * dx + dy * dy); + int numSamples = std::max(2, (int)std::floor(distance / step) + 1); + std::vector profile(numSamples); + for (int i = 0; i < numSamples; i++) { + double t = (numSamples == 1) ? 0 : (double)i / (numSamples - 1); + profile[i] = getElevation(a.x + dx * t, a.y + dy * t); + } + return profile; +} + +} // namespace inet diff --git a/src/inet/common/geometry/common/Heightfield.h b/src/inet/common/geometry/common/Heightfield.h new file mode 100644 index 00000000000..71a7465d66e --- /dev/null +++ b/src/inet/common/geometry/common/Heightfield.h @@ -0,0 +1,91 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#ifndef __INET_HEIGHTFIELD_H +#define __INET_HEIGHTFIELD_H + +#include + +#include "inet/common/geometry/common/Coord.h" + +namespace inet { + +/** + * A regular 2.5D elevation grid (digital surface model) rasterized from a + * point cloud: each grid cell stores the maximum z of the points falling into + * it (DSM semantics — building tops and canopy are preserved), with small + * holes filled by bounded neighbor-averaging passes. Provides O(1) bilinear + * elevation lookup, surface normals, and elevation profiles along a segment + * (for line-of-sight / diffraction models). + * + * Elevation queries outside the data extent (or in unfilled holes) return NaN; + * callers decide how to substitute. Rasterization uses a max-reduce, so the + * result is independent of point order (deterministic across platforms). + */ +class INET_API Heightfield +{ + protected: + double minX = 0, minY = 0; + double cellSize = 0; + int numCellsX = 0, numCellsY = 0; + std::vector cells; // row-major [y * numCellsX + x], NaN = no data + + protected: + float getCell(int ix, int iy) const { return cells[iy * numCellsX + ix]; } + + public: + /** + * Rasterizes the points into a grid. A cellSize <= 0 selects an automatic + * cell size approximating the mean point spacing, sqrt(area / count). + * A positive despikeThreshold clamps isolated single-cell spikes that + * exceed ALL their neighbors by more than the threshold (raw LIDAR tiles + * commonly contain stray high returns — birds, atmospheric noise — which + * max-z rasterization would otherwise keep as phantom needles); genuine + * structures survive because their edge cells have same-height neighbors. + * Throws cRuntimeError if the grid would exceed maxCells (the message + * suggests a coarser cellSize) or if the input is empty/degenerate. + */ + void buildFromPoints(const std::vector& xs, const std::vector& ys, const std::vector& zs, + double cellSize = 0, int64_t maxCells = 67108864, int holeFillPasses = 8, double despikeThreshold = 0); + + bool isValid() const { return numCellsX > 0 && numCellsY > 0; } + double getCellSize() const { return cellSize; } + int getNumCellsX() const { return numCellsX; } + int getNumCellsY() const { return numCellsY; } + double getMinX() const { return minX; } + double getMinY() const { return minY; } + double getMaxX() const { return minX + numCellsX * cellSize; } + double getMaxY() const { return minY + numCellsY * cellSize; } + + /** + * Bilinearly interpolated elevation at (x, y) between cell centers; + * NaN outside the extent or where no data survived hole filling. + */ + double getElevation(double x, double y) const; + + /** + * Upward unit surface normal from central differences of the elevation; + * (NaN, NaN, NaN) where the elevation is undetermined. + */ + Coord getNormal(double x, double y) const; + + /** + * Elevations sampled every 'step' meters along the segment from a to b + * (endpoints included; only x/y of the inputs are used). step must be > 0. + */ + std::vector computeProfile(const Coord& a, const Coord& b, double step) const; + + /** + * The center of the highest cell and its elevation (e.g. the tallest + * building top) — useful for placing nodes relative to landmarks. + */ + Coord getPeakLocation() const; +}; + +} // namespace inet + +#endif diff --git a/src/inet/common/geometry/common/PlyPointCloudReader.cc b/src/inet/common/geometry/common/PlyPointCloudReader.cc new file mode 100644 index 00000000000..29d4f6a60c0 --- /dev/null +++ b/src/inet/common/geometry/common/PlyPointCloudReader.cc @@ -0,0 +1,154 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/common/geometry/common/PlyPointCloudReader.h" + +#include +#include +#include + +namespace inet { + +static int plyTypeSize(const std::string& t) +{ + if (t == "double" || t == "float64") return 8; + if (t == "float" || t == "float32" || t == "int" || t == "int32" || t == "uint" || t == "uint32") return 4; + if (t == "short" || t == "int16" || t == "ushort" || t == "uint16") return 2; + if (t == "char" || t == "int8" || t == "uchar" || t == "uint8") return 1; + return 0; +} + +static double plyRead(const char *p, const std::string& t) +{ + if (t == "double" || t == "float64") { double v; std::memcpy(&v, p, 8); return v; } + if (t == "float" || t == "float32") { float v; std::memcpy(&v, p, 4); return (double)v; } + if (t == "int" || t == "int32") { int32_t v; std::memcpy(&v, p, 4); return (double)v; } + if (t == "uint" || t == "uint32") { uint32_t v; std::memcpy(&v, p, 4); return (double)v; } + if (t == "short" || t == "int16") { int16_t v; std::memcpy(&v, p, 2); return (double)v; } + if (t == "ushort" || t == "uint16") { uint16_t v; std::memcpy(&v, p, 2); return (double)v; } + if (t == "char" || t == "int8") { int8_t v; std::memcpy(&v, p, 1); return (double)v; } + if (t == "uchar" || t == "uint8") { uint8_t v; std::memcpy(&v, p, 1); return (double)v; } + return 0.0; +} + +PlyPointCloud PlyPointCloudReader::read(const std::string& path) +{ + std::ifstream in(path, std::ios::binary); + if (!in) + throw cRuntimeError("Cannot open PLY file '%s' (resolved relative to the working directory)", path.c_str()); + + // --- header --- + std::string line; + std::getline(in, line); + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.rfind("ply", 0) != 0) + throw cRuntimeError("'%s' is not a PLY file (missing 'ply' magic)", path.c_str()); + std::string format; + int vertexCount = 0; + std::vector> props; // (type, name) of the vertex element, in order + std::string curElement; + while (std::getline(in, line)) { + if (!line.empty() && line.back() == '\r') line.pop_back(); + std::istringstream ss(line); + std::string tok; ss >> tok; + if (tok == "format") ss >> format; + else if (tok == "element") { ss >> curElement; if (curElement == "vertex") ss >> vertexCount; } + else if (tok == "property" && curElement == "vertex") { + std::string type; ss >> type; + if (type == "list") continue; // e.g. face vertex-index lists — not a vertex scalar + std::string name; ss >> name; + props.emplace_back(type, name); + } + else if (tok == "end_header") break; + } + bool ascii = (format.rfind("ascii", 0) == 0); + bool binaryLE = (format.rfind("binary_little_endian", 0) == 0); + if (!ascii && !binaryLE) + throw cRuntimeError("'%s' has unsupported PLY format '%s' (only ascii and binary_little_endian are supported)", path.c_str(), format.c_str()); + if (vertexCount <= 0 || props.empty()) + throw cRuntimeError("'%s' has no readable vertex element", path.c_str()); + + // locate x/y/z (+ optional r/g/b) by property name; compute byte offsets and column indices + int stride = 0, offX = -1, offY = -1, offZ = -1, offR = -1, offG = -1, offB = -1; + int idxX = -1, idxY = -1, idxZ = -1, idxR = -1, idxG = -1, idxB = -1; + std::string tX, tY, tZ, tR, tG, tB; + for (size_t i = 0; i < props.size(); i++) { + const std::string& ty = props[i].first; + const std::string& n = props[i].second; + if (n == "x") { offX = stride; tX = ty; idxX = (int)i; } + else if (n == "y") { offY = stride; tY = ty; idxY = (int)i; } + else if (n == "z") { offZ = stride; tZ = ty; idxZ = (int)i; } + else if (n == "red" || n == "r") { offR = stride; tR = ty; idxR = (int)i; } + else if (n == "green" || n == "g") { offG = stride; tG = ty; idxG = (int)i; } + else if (n == "blue" || n == "b") { offB = stride; tB = ty; idxB = (int)i; } + stride += plyTypeSize(ty); + } + if (offX < 0 || offY < 0 || offZ < 0) + throw cRuntimeError("'%s' has no x/y/z vertex properties", path.c_str()); + PlyPointCloud cloud; + cloud.hasRGB = (offR >= 0 && offG >= 0 && offB >= 0); + // 8-bit channels are 0..255, float channels 0..1 — normalise each channel by ITS OWN type. + auto colScale = [](const std::string& t) { return (t == "uchar" || t == "uint8") ? 1.0 / 255.0 : 1.0; }; + double scaleR = colScale(tR), scaleG = colScale(tG), scaleB = colScale(tB); + + // --- read the vertices --- + cloud.xs.resize(vertexCount); + cloud.ys.resize(vertexCount); + cloud.zs.resize(vertexCount); + if (cloud.hasRGB) { + cloud.rs.resize(vertexCount); + cloud.gs.resize(vertexCount); + cloud.bs.resize(vertexCount); + } + if (binaryLE) { + std::vector buf(stride); + for (int i = 0; i < vertexCount; i++) { + in.read(buf.data(), stride); + if (!in) + throw cRuntimeError("'%s' is truncated (expected %d vertices)", path.c_str(), vertexCount); + cloud.xs[i] = plyRead(buf.data() + offX, tX); + cloud.ys[i] = plyRead(buf.data() + offY, tY); + cloud.zs[i] = plyRead(buf.data() + offZ, tZ); + if (cloud.hasRGB) { + cloud.rs[i] = plyRead(buf.data() + offR, tR) * scaleR; + cloud.gs[i] = plyRead(buf.data() + offG, tG) * scaleG; + cloud.bs[i] = plyRead(buf.data() + offB, tB) * scaleB; + } + } + } + else { // ascii + for (int i = 0; i < vertexCount; i++) { + if (!std::getline(in, line)) + throw cRuntimeError("'%s' is truncated (expected %d vertices)", path.c_str(), vertexCount); + std::istringstream ss(line); + std::vector v; double d; while (ss >> d) v.push_back(d); + if ((int)v.size() < (int)props.size()) + throw cRuntimeError("'%s' has a short vertex line (%d values, expected %d)", path.c_str(), (int)v.size(), (int)props.size()); + cloud.xs[i] = v[idxX]; + cloud.ys[i] = v[idxY]; + cloud.zs[i] = v[idxZ]; + if (cloud.hasRGB) { + cloud.rs[i] = v[idxR] * scaleR; + cloud.gs[i] = v[idxG] * scaleG; + cloud.bs[i] = v[idxB] * scaleB; + } + } + } + + // --- bounding box --- + cloud.minX = cloud.maxX = cloud.xs[0]; + cloud.minY = cloud.maxY = cloud.ys[0]; + cloud.minZ = cloud.maxZ = cloud.zs[0]; + for (int i = 1; i < vertexCount; i++) { + cloud.minX = std::min(cloud.minX, cloud.xs[i]); cloud.maxX = std::max(cloud.maxX, cloud.xs[i]); + cloud.minY = std::min(cloud.minY, cloud.ys[i]); cloud.maxY = std::max(cloud.maxY, cloud.ys[i]); + cloud.minZ = std::min(cloud.minZ, cloud.zs[i]); cloud.maxZ = std::max(cloud.maxZ, cloud.zs[i]); + } + return cloud; +} + +} // namespace inet diff --git a/src/inet/common/geometry/common/PlyPointCloudReader.h b/src/inet/common/geometry/common/PlyPointCloudReader.h new file mode 100644 index 00000000000..9e88385ef6f --- /dev/null +++ b/src/inet/common/geometry/common/PlyPointCloudReader.h @@ -0,0 +1,60 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#ifndef __INET_PLYPOINTCLOUDREADER_H +#define __INET_PLYPOINTCLOUDREADER_H + +#include +#include + +#include "inet/common/INETDefs.h" + +namespace inet { + +/** + * A point cloud loaded from a PLY file (e.g. a LIDAR scan). Coordinates are + * kept exactly as stored in the file (no recentering or scaling); consumers + * apply their own transform. Color channels, when present, are normalized to + * [0,1] (8-bit channels divided by 255, wider/float channels taken as-is). + */ +class INET_API PlyPointCloud +{ + public: + std::vector xs; + std::vector ys; + std::vector zs; + std::vector rs; // empty unless hasRGB + std::vector gs; // empty unless hasRGB + std::vector bs; // empty unless hasRGB + bool hasRGB = false; + double minX = 0, maxX = 0; + double minY = 0, maxY = 0; + double minZ = 0, maxZ = 0; + + int getNumPoints() const { return (int)xs.size(); } +}; + +/** + * Reads PLY point clouds: ascii and binary_little_endian formats, vertex + * properties located by name (x/y/z required, red/green/blue or r/g/b + * optional) with arbitrary scalar types and property order. List properties + * (e.g. face vertex-index lists) are skipped — only points are read. + * Throws cRuntimeError with a descriptive message on malformed input. + * + * Extracted from the VSG scene visualizer's terrain loader so that both the + * visualization and the physical models (ground/obstacle) read LIDAR data + * through one implementation. + */ +class INET_API PlyPointCloudReader +{ + public: + static PlyPointCloud read(const std::string& path); +}; + +} // namespace inet + +#endif diff --git a/tests/unit/Heightfield_1.test b/tests/unit/Heightfield_1.test new file mode 100644 index 00000000000..85b53e98021 --- /dev/null +++ b/tests/unit/Heightfield_1.test @@ -0,0 +1,121 @@ +%description: +Test Heightfield: bilinear elevation on a flat plane, normals on flat and +ramp surfaces, elevation profile across a step (building wall), out-of-extent +NaN behavior, and the maxCells guard. + +%includes: +#include +#include "inet/common/geometry/common/Heightfield.h" + +using namespace inet; + +%global: + +#define CHECK(label, cond) EV << label << ": " << ((cond) ? "OK" : "FAIL") << "\n"; + +static void makeGrid(std::vector& xs, std::vector& ys, std::vector& zs, + double x0, double x1, double y0, double y1, double spacing, double (*zf)(double, double)) +{ + for (double x = x0; x <= x1 + 1e-9; x += spacing) + for (double y = y0; y <= y1 + 1e-9; y += spacing) { + xs.push_back(x); + ys.push_back(y); + zs.push_back(zf(x, y)); + } +} + +static double zFlat(double, double) { return 5; } +static double zRamp(double x, double) { return x; } +static double zStep(double x, double) { return x < 10 ? 0 : 30; } + +%activity: + +// --- flat plane at z=5: bilinear elevation is exact, normal is straight up --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 10, 0, 10, 0.5, zFlat); + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 1); + CHECK("flat elevation", std::fabs(hf.getElevation(3.3, 4.7) - 5) < 1e-6); + Coord n = hf.getNormal(5, 5); + CHECK("flat normal", n.z > 0.999 && std::fabs(n.x) < 1e-6 && std::fabs(n.y) < 1e-6); + CHECK("out of extent NaN", std::isnan(hf.getElevation(-100, -100))); +} + +// --- ramp z=x: elevation increases along x, normal tilts against the slope --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 20, 0, 10, 0.25, zRamp); + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 1); + CHECK("ramp monotonic", hf.getElevation(15, 5) > hf.getElevation(5, 5) + 5); + Coord n = hf.getNormal(10, 5); + CHECK("ramp normal tilt", n.x < -0.5 && n.z > 0.5 && std::fabs(n.y) < 0.1); +} + +// --- step (building wall) at x=10: profile from low to high side crosses it --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 20, 0, 20, 0.5, zStep); + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 1); + CHECK("step low side", std::fabs(hf.getElevation(3, 10)) < 1e-6); + CHECK("step high side", std::fabs(hf.getElevation(17, 10) - 30) < 1e-6); + std::vector profile = hf.computeProfile(Coord(2, 10, 0), Coord(18, 10, 0), 1); + double lo = profile[0], hi = profile[0]; + for (double e : profile) { lo = std::min(lo, e); hi = std::max(hi, e); } + CHECK("profile spans the wall", profile.size() >= 16 && lo < 1 && hi > 29 && profile.front() < 1 && profile.back() > 29); +} + +// --- maxCells guard rejects oversized grids with a helpful error --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 100, 0, 100, 10, zFlat); + Heightfield hf; + try { + hf.buildFromPoints(xs, ys, zs, 0.01, 10); + EV << "maxCells guard: FAIL\n"; + } catch (cRuntimeError& e) { + EV << "maxCells guard: OK\n"; + } +} + +// --- despike: an isolated 1-cell spike is clamped, a real multi-cell building survives --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 20, 0, 20, 0.5, zStep); // "building" occupying x>=10 at z=30 + xs.push_back(5); ys.push_back(5); zs.push_back(500); // stray LIDAR return + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 1, 67108864, 8, 50); + CHECK("spike clamped", hf.getElevation(5, 5) < 10); + CHECK("building survives despike", std::fabs(hf.getElevation(17, 10) - 30) < 1e-6); +} + +// --- degenerate input rejected --- +{ + std::vector xs(3, 1), ys(3, 1), zs(3, 1); // all points at the same x/y + Heightfield hf; + try { + hf.buildFromPoints(xs, ys, zs, 1); + EV << "degenerate rejected: FAIL\n"; + } catch (cRuntimeError& e) { + EV << "degenerate rejected: OK\n"; + } +} + +EV << ".\n"; + +%contains: stdout +flat elevation: OK +flat normal: OK +out of extent NaN: OK +ramp monotonic: OK +ramp normal tilt: OK +step low side: OK +step high side: OK +profile spans the wall: OK +maxCells guard: OK +spike clamped: OK +building survives despike: OK +degenerate rejected: OK +. diff --git a/tests/unit/PlyPointCloudReader_1.test b/tests/unit/PlyPointCloudReader_1.test new file mode 100644 index 00000000000..193644f82ad --- /dev/null +++ b/tests/unit/PlyPointCloudReader_1.test @@ -0,0 +1,114 @@ +%description: +Test PlyPointCloudReader: ascii and binary_little_endian formats, property +location by name with mixed types and extra columns, per-channel color +normalization, bounding box, and rejection of malformed files. + +%includes: +#include +#include "inet/common/geometry/common/PlyPointCloudReader.h" + +using namespace inet; + +%global: + +static void writeTextFile(const char *name, const char *content) +{ + std::ofstream out(name, std::ios::binary); + out << content; +} + +#define CHECK(label, cond) EV << label << ": " << ((cond) ? "OK" : "FAIL") << "\n"; + +%activity: + +// --- ascii PLY: float/double coordinates, uchar RGB, extra (ignored) column --- +writeTextFile("ascii.ply", + "ply\n" + "format ascii 1.0\n" + "comment synthetic test cloud\n" + "element vertex 3\n" + "property float x\n" + "property float y\n" + "property double z\n" + "property uchar red\n" + "property uchar green\n" + "property uchar blue\n" + "property float intensity\n" + "end_header\n" + "0 0 1 255 0 0 42\n" + "10 0 2 0 255 0 42\n" + "10 20 3.5 0 0 255 42\n"); +PlyPointCloud a = PlyPointCloudReader::read("ascii.ply"); +CHECK("ascii count", a.getNumPoints() == 3); +CHECK("ascii coords", a.xs[2] == 10 && a.ys[2] == 20 && a.zs[2] == 3.5); +CHECK("ascii rgb normalized", a.hasRGB && a.rs[0] == 1 && a.gs[0] == 0 && a.bs[0] == 0 && a.gs[1] == 1 && a.bs[2] == 1); +CHECK("ascii bbox", a.minX == 0 && a.maxX == 10 && a.minY == 0 && a.maxY == 20 && a.minZ == 1 && a.maxZ == 3.5); + +// --- binary_little_endian PLY: double x/y + float z, no color --- +{ + std::ofstream out("binary.ply", std::ios::binary); + out << "ply\n" + "format binary_little_endian 1.0\n" + "element vertex 2\n" + "property double x\n" + "property double y\n" + "property float z\n" + "end_header\n"; + double xy[2]; + float z; + xy[0] = 1.25; xy[1] = -2.5; z = 10.0f; + out.write((char *)xy, 16); out.write((char *)&z, 4); + xy[0] = 3.0; xy[1] = 4.0; z = -1.5f; + out.write((char *)xy, 16); out.write((char *)&z, 4); +} +PlyPointCloud b = PlyPointCloudReader::read("binary.ply"); +CHECK("binary count", b.getNumPoints() == 2); +CHECK("binary coords", b.xs[0] == 1.25 && b.ys[0] == -2.5 && b.zs[0] == 10 && b.zs[1] == -1.5); +CHECK("binary no rgb", !b.hasRGB); +CHECK("binary bbox", b.minX == 1.25 && b.maxX == 3 && b.minZ == -1.5 && b.maxZ == 10); + +// --- malformed inputs are rejected with cRuntimeError --- +writeTextFile("nomagic.ply", "solid nonsense\n"); +try { + PlyPointCloudReader::read("nomagic.ply"); + EV << "missing magic: FAIL\n"; +} catch (cRuntimeError& e) { + EV << "missing magic: OK\n"; +} +writeTextFile("bigendian.ply", "ply\nformat binary_big_endian 1.0\nelement vertex 1\nproperty float x\nproperty float y\nproperty float z\nend_header\n"); +try { + PlyPointCloudReader::read("bigendian.ply"); + EV << "big endian rejected: FAIL\n"; +} catch (cRuntimeError& e) { + EV << "big endian rejected: OK\n"; +} +writeTextFile("truncated.ply", "ply\nformat ascii 1.0\nelement vertex 5\nproperty float x\nproperty float y\nproperty float z\nend_header\n1 2 3\n"); +try { + PlyPointCloudReader::read("truncated.ply"); + EV << "truncation detected: FAIL\n"; +} catch (cRuntimeError& e) { + EV << "truncation detected: OK\n"; +} +try { + PlyPointCloudReader::read("doesnotexist.ply"); + EV << "missing file: FAIL\n"; +} catch (cRuntimeError& e) { + EV << "missing file: OK\n"; +} + +EV << ".\n"; + +%contains: stdout +ascii count: OK +ascii coords: OK +ascii rgb normalized: OK +ascii bbox: OK +binary count: OK +binary coords: OK +binary no rgb: OK +binary bbox: OK +missing magic: OK +big endian rejected: OK +truncation detected: OK +missing file: OK +. From 072e1767f4ba19d918c4c0de4191e22af94497db Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 18:52:26 +0300 Subject: [PATCH 06/17] visualizer: vsg: load sceneModel terrain through the shared PLY reader Behavior-preserving refactor: createTerrainFromPLY now reads the point cloud via PlyPointCloudReader (colors arrive normalized) and keeps the display transform (bbox recenter + aspect-fit into the scene bounds) and elevation coloring unchanged, so the physical terrain models and the visualization consume LIDAR data through one implementation. Co-Authored-By: Claude Fable 5 --- src/inet/visualizer/vsg/util/VsgUtils.cc | 123 +++-------------------- 1 file changed, 13 insertions(+), 110 deletions(-) diff --git a/src/inet/visualizer/vsg/util/VsgUtils.cc b/src/inet/visualizer/vsg/util/VsgUtils.cc index 12723570835..45a40cc3755 100644 --- a/src/inet/visualizer/vsg/util/VsgUtils.cc +++ b/src/inet/visualizer/vsg/util/VsgUtils.cc @@ -15,6 +15,8 @@ #include #include +#include "inet/common/geometry/common/PlyPointCloudReader.h" + namespace inet { namespace vsg { @@ -708,24 +710,6 @@ static ref_ptr getPointCloudPipeline() return pipeline; } -static int plyTypeSize(const std::string& t) { - if (t == "double" || t == "float64") return 8; - if (t == "float" || t == "float32" || t == "int" || t == "int32" || t == "uint" || t == "uint32") return 4; - if (t == "short" || t == "int16" || t == "ushort" || t == "uint16") return 2; - if (t == "char" || t == "int8" || t == "uchar" || t == "uint8") return 1; - return 0; -} -static double plyRead(const char* p, const std::string& t) { - if (t == "double" || t == "float64") { double v; std::memcpy(&v, p, 8); return v; } - if (t == "float" || t == "float32") { float v; std::memcpy(&v, p, 4); return (double)v; } - if (t == "int" || t == "int32") { int32_t v; std::memcpy(&v, p, 4); return (double)v; } - if (t == "uint" || t == "uint32") { uint32_t v; std::memcpy(&v, p, 4); return (double)v; } - if (t == "short" || t == "int16") { int16_t v; std::memcpy(&v, p, 2); return (double)v; } - if (t == "ushort" || t == "uint16") { uint16_t v; std::memcpy(&v, p, 2); return (double)v; } - if (t == "char" || t == "int8") { int8_t v; std::memcpy(&v, p, 1); return (double)v; } - if (t == "uchar" || t == "uint8") { uint8_t v; std::memcpy(&v, p, 1); return (double)v; } - return 0.0; -} // Terrain elevation ramp for a normalised height t in [0,1]: teal (low) -> green -> tan -> brown -> white. static ::vsg::vec4 elevationColor(double t) { t = std::max(0.0, std::min(1.0, t)); @@ -741,99 +725,17 @@ static ::vsg::vec4 elevationColor(double t) { ref_ptr createTerrainFromPLY(const std::string& path, const Coord& sceneMin, const Coord& sceneMax) { - std::ifstream in(path, std::ios::binary); - if (!in) - throw cRuntimeError("sceneModel: cannot open PLY file '%s' (resolved relative to the working directory)", path.c_str()); - - // --- header --- - std::string line; - std::getline(in, line); - if (!line.empty() && line.back() == '\r') line.pop_back(); - if (line.rfind("ply", 0) != 0) - throw cRuntimeError("sceneModel: '%s' is not a PLY file (missing 'ply' magic)", path.c_str()); - std::string format; - int vertexCount = 0; - std::vector> props; // (type, name) of the vertex element, in order - std::string curElement; - while (std::getline(in, line)) { - if (!line.empty() && line.back() == '\r') line.pop_back(); - std::istringstream ss(line); - std::string tok; ss >> tok; - if (tok == "format") ss >> format; - else if (tok == "element") { ss >> curElement; if (curElement == "vertex") ss >> vertexCount; } - else if (tok == "property" && curElement == "vertex") { - std::string type; ss >> type; - if (type == "list") continue; // e.g. face vertex-index lists — not a vertex scalar - std::string name; ss >> name; - props.emplace_back(type, name); - } - else if (tok == "end_header") break; - } - bool ascii = (format.rfind("ascii", 0) == 0); - bool binaryLE = (format.rfind("binary_little_endian", 0) == 0); - if (!ascii && !binaryLE) - throw cRuntimeError("sceneModel: '%s' has unsupported PLY format '%s' (only ascii and binary_little_endian are supported)", path.c_str(), format.c_str()); - if (vertexCount <= 0 || props.empty()) - throw cRuntimeError("sceneModel: '%s' has no readable vertex element", path.c_str()); - - // locate x/y/z (+ optional r/g/b) by property name; compute byte offsets and column indices - int stride = 0, offX = -1, offY = -1, offZ = -1, offR = -1, offG = -1, offB = -1; - int idxX = -1, idxY = -1, idxZ = -1, idxR = -1, idxG = -1, idxB = -1; - std::string tX, tY, tZ, tR, tG, tB; - for (size_t i = 0; i < props.size(); i++) { - const std::string& ty = props[i].first; - const std::string& n = props[i].second; - if (n == "x") { offX = stride; tX = ty; idxX = (int)i; } - else if (n == "y") { offY = stride; tY = ty; idxY = (int)i; } - else if (n == "z") { offZ = stride; tZ = ty; idxZ = (int)i; } - else if (n == "red" || n == "r") { offR = stride; tR = ty; idxR = (int)i; } - else if (n == "green" || n == "g") { offG = stride; tG = ty; idxG = (int)i; } - else if (n == "blue" || n == "b") { offB = stride; tB = ty; idxB = (int)i; } - stride += plyTypeSize(ty); - } - if (offX < 0 || offY < 0 || offZ < 0) - throw cRuntimeError("sceneModel: '%s' has no x/y/z vertex properties", path.c_str()); - bool hasRGB = (offR >= 0 && offG >= 0 && offB >= 0); - // 8-bit channels are 0..255, float channels 0..1 — normalise each channel by ITS OWN type. - auto colScale = [](const std::string& t) { return (t == "uchar" || t == "uint8") ? 1.0 / 255.0 : 1.0; }; - double scaleR = colScale(tR), scaleG = colScale(tG), scaleB = colScale(tB); - - // --- read the vertices --- - std::vector xs(vertexCount), ys(vertexCount), zs(vertexCount); - std::vector<::vsg::vec4> rgbs(hasRGB ? vertexCount : 0); - if (binaryLE) { - std::vector buf(stride); - for (int i = 0; i < vertexCount; i++) { - in.read(buf.data(), stride); - if (!in) - throw cRuntimeError("sceneModel: '%s' is truncated (expected %d vertices)", path.c_str(), vertexCount); - xs[i] = plyRead(buf.data() + offX, tX); ys[i] = plyRead(buf.data() + offY, tY); zs[i] = plyRead(buf.data() + offZ, tZ); - if (hasRGB) - rgbs[i] = ::vsg::vec4((float)(plyRead(buf.data() + offR, tR) * scaleR), - (float)(plyRead(buf.data() + offG, tG) * scaleG), - (float)(plyRead(buf.data() + offB, tB) * scaleB), 1.0f); - } - } - else { // ascii - for (int i = 0; i < vertexCount; i++) { - if (!std::getline(in, line)) - throw cRuntimeError("sceneModel: '%s' is truncated (expected %d vertices)", path.c_str(), vertexCount); - std::istringstream ss(line); - std::vector v; double d; while (ss >> d) v.push_back(d); - if ((int)v.size() < (int)props.size()) - throw cRuntimeError("sceneModel: '%s' has a short vertex line (%d values, expected %d)", path.c_str(), (int)v.size(), (int)props.size()); - xs[i] = v[idxX]; ys[i] = v[idxY]; zs[i] = v[idxZ]; - if (hasRGB) rgbs[i] = ::vsg::vec4((float)(v[idxR] * scaleR), (float)(v[idxG] * scaleG), (float)(v[idxB] * scaleB), 1.0f); - } - } + // shared PLY reader (also used by the physical terrain models); coordinates + // arrive untransformed, colors (if any) already normalized to [0,1] + PlyPointCloud cloud = PlyPointCloudReader::read(path); + int vertexCount = cloud.getNumPoints(); + const std::vector& xs = cloud.xs; + const std::vector& ys = cloud.ys; + const std::vector& zs = cloud.zs; + bool hasRGB = cloud.hasRGB; // --- recentre, fit into the scene box (aspect-preserving), colour by elevation if no RGB --- - double minX = xs[0], maxX = xs[0], minY = ys[0], maxY = ys[0], minZ = zs[0], maxZ = zs[0]; - for (int i = 1; i < vertexCount; i++) { - minX = std::min(minX, xs[i]); maxX = std::max(maxX, xs[i]); - minY = std::min(minY, ys[i]); maxY = std::max(maxY, ys[i]); - minZ = std::min(minZ, zs[i]); maxZ = std::max(maxZ, zs[i]); - } + double minX = cloud.minX, maxX = cloud.maxX, minY = cloud.minY, maxY = cloud.maxY, minZ = cloud.minZ, maxZ = cloud.maxZ; double cx = (minX + maxX) / 2, cy = (minY + maxY) / 2, dz = (maxZ - minZ); double plyW = maxX - minX, plyH = maxY - minY; double sceneW = sceneMax.x - sceneMin.x, sceneH = sceneMax.y - sceneMin.y; @@ -850,7 +752,8 @@ ref_ptr createTerrainFromPLY(const std::string& path, const Coord& sceneMi vertices->set(i, ::vsg::vec3((float)((xs[i] - cx) * scale + centerX), (float)((ys[i] - cy) * scale + centerY), (float)((zs[i] - minZ) * scale + baseZ))); - colors->set(i, hasRGB ? rgbs[i] : elevationColor(dz > 0 ? (zs[i] - minZ) / dz : 0.0)); + colors->set(i, hasRGB ? ::vsg::vec4((float)cloud.rs[i], (float)cloud.gs[i], (float)cloud.bs[i], 1.0f) + : elevationColor(dz > 0 ? (zs[i] - minZ) / dz : 0.0)); } // --- build the node (fit baked into the vertices, so no transform needed) --- From 97e3fdbb29fba891ed3e9a03f774fd2961b05337 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 18:52:26 +0300 Subject: [PATCH 07/17] =?UTF-8?q?environment:=20PointCloudGround=20?= =?UTF-8?q?=E2=80=94=20LIDAR=20PLY=20terrain=20as=20the=20physical=20groun?= =?UTF-8?q?d?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New IGround implementation rasterizing a PLY point cloud (e.g. a LIDAR scan) into a heightfield, so ground-aware radio models use the real surface elevation. Fulfills the environment/ground/__TODO wishes for a renderer- independent, file-loaded ground. Transform modes: "metric" (default; translate-only, 1 PLY unit = 1 sim meter, physically exact), "fit" (reproduces the VSG scene visualizer's recenter + aspect-fit display transform so physics matches the displayed terrain; logs a scaling warning), "manual" (explicit scale/offsets). Out-of-extent queries return outOfBoundsElevation (NaN by default per the IGround contract). Init log reports grid size, extent and peak location for placing nodes. vsglidar example: new TerrainTwoRay config ([General] untouched) — a gcs ground station on an open street near the tallest building (289m tower at (563,277)), five drones flying among the buildings and reporting to it, and TwoRayGroundReflection computing path loss over the real surface. The physicalEnvironment space is pinned to the scene box because the scene visualizer derives its bounds from it once a PhysicalEnvironment module exists (NaN space otherwise shifts the displayed terrain). Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/VsgLidar.ned | 22 ++++ examples/visualizer/vsglidar/omnetpp.ini | 82 ++++++++++++++ .../environment/ground/PointCloudGround.cc | 103 ++++++++++++++++++ .../environment/ground/PointCloudGround.h | 45 ++++++++ .../environment/ground/PointCloudGround.ned | 61 +++++++++++ 5 files changed, 313 insertions(+) create mode 100644 src/inet/environment/ground/PointCloudGround.cc create mode 100644 src/inet/environment/ground/PointCloudGround.h create mode 100644 src/inet/environment/ground/PointCloudGround.ned diff --git a/examples/visualizer/vsglidar/VsgLidar.ned b/examples/visualizer/vsglidar/VsgLidar.ned index 10091aeb386..d77966668ca 100644 --- a/examples/visualizer/vsglidar/VsgLidar.ned +++ b/examples/visualizer/vsglidar/VsgLidar.ned @@ -6,6 +6,7 @@ package inet.examples.visualizer.vsglidar; +import inet.environment.common.PhysicalEnvironment; import inet.networklayer.configurator.ipv4.Ipv4NetworkConfigurator; import inet.node.inet.WirelessHost; import inet.physicallayer.wireless.ieee80211.packetlevel.Ieee80211ScalarRadioMedium; @@ -53,3 +54,24 @@ network VsgLidarShowcase @display("p=150,330;i=misc/drone"); } } + +// +// The same drone swarm, but with a physical environment so the LIDAR terrain +// can also act on the radio physics (see the TerrainTwoRay config): the +// physicalEnvironment.ground submodule loads the SAME terrain.ply through +// ~PointCloudGround, giving ground-aware path-loss models (e.g. +// ~TwoRayGroundReflection) the real surface elevation under each drone. +// A stationary ground station (gcs) sits at street level next to the tallest +// building; the drones report to it, so drone-to-base links pass right by the +// tower — the natural test case for terrain occlusion models. +// +network VsgLidarTerrainShowcase extends VsgLidarShowcase +{ + submodules: + physicalEnvironment: PhysicalEnvironment { + @display("p=60,300"); + } + gcs: WirelessHost { // ground control station at the foot of the tower + @display("p=400,400;i=device/antennatower"); + } +} diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index a6e6d414da5..655967f7eba 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -124,3 +124,85 @@ description = "Drone swarm over a real LIDAR landscape: PLY point-cloud terrain *.visualizer.mobilityVisualizer.displayMovementTrails = true *.visualizer.mobilityVisualizer.animationSpeed = 1 + +# -------------------------------------------------------------------------------------------- +# The LIDAR terrain as PHYSICS, not just scenery: the same terrain.ply is loaded a second time +# by physicalEnvironment.ground (PointCloudGround), and the two-ray ground-reflection path-loss +# model then uses the real surface elevation under each drone instead of a flat ground. +# transformMode "fit" reproduces the scene visualizer's display transform (recenter + aspect-fit +# into the scene bounds), so the physical terrain is exactly the surface you SEE the drones fly +# over; for physically exact studies use the default "metric" mode (1 PLY meter = 1 sim meter). +# -------------------------------------------------------------------------------------------- +[Config TerrainTwoRay] +network = VsgLidarTerrainShowcase +description = "Two-ray ground reflection over the real LIDAR surface (PointCloudGround)" +# Pin the physical environment's space to the scene box. Without this the scene visualizer +# derives its bounds from physicalEnvironment.getSpaceMin/Max (NaN by default with no objects) +# instead of the 800x450 background, which shifts the displayed terrain away from both the +# node positions and the physics heightfield. +*.physicalEnvironment.spaceMinX = 0m +*.physicalEnvironment.spaceMinY = 0m +*.physicalEnvironment.spaceMinZ = 0m +*.physicalEnvironment.spaceMaxX = 800m +*.physicalEnvironment.spaceMaxY = 450m +*.physicalEnvironment.spaceMaxZ = 400m +*.physicalEnvironment.ground.typename = "PointCloudGround" +*.physicalEnvironment.ground.file = "terrain.ply" +*.physicalEnvironment.ground.transformMode = "fit" +*.physicalEnvironment.ground.fitMinX = 0m +*.physicalEnvironment.ground.fitMinY = 0m +*.physicalEnvironment.ground.fitMinZ = 0m +*.physicalEnvironment.ground.fitMaxX = 800m +*.physicalEnvironment.ground.fitMaxY = 450m +*.radioMedium.pathLoss.typename = "TwoRayGroundReflection" + +# Unlike [General] (drones above the tower, permanent line of sight everywhere), fly the swarm +# AMONG the buildings so the terrain matters: the altitude above the local surface — and with it +# the two-ray path loss — changes as drones cross streets, rooftops and the tower. +# The terrain heightfield spans x=[0,801] y=[5,445]; its highest structure (the tower, ~289m +# after despiking stray LIDAR returns) stands at (563,277) — see the PointCloudGround init log. +# The swarm is kept inside that footprint, and drones 3-4 start on the FAR side of the tower +# from the gcs, so their reports pass right by it (the occlusion test case for later phases). +# NOTE: the tile's lowest point (mapped to z=0) is not street level — the streets across this +# scan sit at z ≈ 41-55m in scene coordinates. Altitudes below keep everything above ground. +*.drone*.mobility.constraintAreaMinX = 60m +*.drone*.mobility.constraintAreaMinY = 40m +*.drone*.mobility.constraintAreaMinZ = 70m +*.drone*.mobility.constraintAreaMaxX = 740m +*.drone*.mobility.constraintAreaMaxY = 410m +*.drone*.mobility.constraintAreaMaxZ = 240m +*.drone1.mobility.initialX = 300m +*.drone1.mobility.initialY = 150m +*.drone1.mobility.initialZ = 120m +*.drone2.mobility.initialX = 150m +*.drone2.mobility.initialY = 350m +*.drone2.mobility.initialZ = 90m +*.drone3.mobility.initialX = 680m +*.drone3.mobility.initialY = 240m +*.drone3.mobility.initialZ = 80m +*.drone4.mobility.initialX = 640m +*.drone4.mobility.initialY = 180m +*.drone4.mobility.initialZ = 100m +*.drone5.mobility.initialX = 400m +*.drone5.mobility.initialY = 80m +*.drone5.mobility.initialZ = 160m + +# --- ground control station at the foot of the tower; all drones report to it, so every link +# crosses the built-up center of the map --- +*.gcs.mobility.initFromDisplayString = false +*.gcs.mobility.initialX = 520m +*.gcs.mobility.initialY = 360m +*.gcs.mobility.initialZ = 55m # open street south of the tower (~41m grade) + 14m antenna mast, clear of the ~50m low roofs +*.gcs.numApps = 1 +*.gcs.app[0].typename = "UdpSink" +*.gcs.app[0].localPort = 1000 +*.gcs.wlan[*].radio.transmitter.power = 10mW +*.gcs.wlan[*].bitrate = 54Mbps +*.gcs.wlan[*].mgmt.typename = "Ieee80211MgmtAdhoc" +*.gcs.wlan[*].agent.typename = "" +*.drone1.app[1].destAddresses = "gcs" +*.drone1.app[1].destPort = 1000 +*.drone2.app[1].destAddresses = "gcs" +*.drone3.app[0].destAddresses = "gcs" +*.drone4.app[0].destAddresses = "gcs" +*.drone5.app[0].destAddresses = "gcs" diff --git a/src/inet/environment/ground/PointCloudGround.cc b/src/inet/environment/ground/PointCloudGround.cc new file mode 100644 index 00000000000..6faa1c3b8c3 --- /dev/null +++ b/src/inet/environment/ground/PointCloudGround.cc @@ -0,0 +1,103 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/environment/ground/PointCloudGround.h" + +#include +#include + +#include "inet/common/geometry/common/PlyPointCloudReader.h" + +namespace inet { + +namespace physicalenvironment { + +Define_Module(PointCloudGround); + +void PointCloudGround::initialize(int stage) +{ + if (stage == INITSTAGE_LOCAL) { + const char *file = par("file"); + double cellSize = par("cellSize"); + int64_t maxCells = par("maxCells"); + outOfBoundsElevation = par("outOfBoundsElevation"); + const char *transformMode = par("transformMode"); + + long startTime = clock(); + PlyPointCloud cloud = PlyPointCloudReader::read(file); + + // reduce every mode to point' = point * scale + offset (per axis) + double scale = 1, offsetX = 0, offsetY = 0, offsetZ = 0; + if (!strcmp(transformMode, "metric")) { + // translate only: 1 PLY unit = 1 simulation meter (physics-correct) + scale = 1; + offsetX = par("originX").doubleValue() - cloud.minX; + offsetY = par("originY").doubleValue() - cloud.minY; + offsetZ = par("baseElevation").doubleValue() - cloud.minZ; + } + else if (!strcmp(transformMode, "fit")) { + // reproduce the VSG scene visualizer's display transform: recenter the + // cloud's bbox onto the scene center and aspect-fit it into the scene + // footprint; fitMin/fitMax must match the visualizer's scene bounds + double fitMinX = par("fitMinX"), fitMinY = par("fitMinY"), fitMinZ = par("fitMinZ"); + double fitMaxX = par("fitMaxX"), fitMaxY = par("fitMaxY"); + if (std::isnan(fitMinX) || std::isnan(fitMinY) || std::isnan(fitMinZ) || std::isnan(fitMaxX) || std::isnan(fitMaxY)) + throw cRuntimeError("transformMode=\"fit\" requires the fitMinX/fitMinY/fitMinZ/fitMaxX/fitMaxY parameters (set them to the scene visualizer's scene bounds)"); + double plyW = cloud.maxX - cloud.minX, plyH = cloud.maxY - cloud.minY; + double sceneW = fitMaxX - fitMinX, sceneH = fitMaxY - fitMinY; + scale = (plyW > 0 && plyH > 0 && sceneW > 0 && sceneH > 0) ? std::min(sceneW / plyW, sceneH / plyH) : 1.0; + double cx = (cloud.minX + cloud.maxX) / 2, cy = (cloud.minY + cloud.maxY) / 2; + offsetX = (fitMinX + fitMaxX) / 2 - cx * scale; + offsetY = (fitMinY + fitMaxY) / 2 - cy * scale; + offsetZ = fitMinZ - cloud.minZ * scale; + EV_WARN << "PointCloudGround: transformMode=\"fit\" scales the terrain by " << scale + << " to fit the scene; physical distances over the terrain are scaled accordingly — use \"metric\" for physically exact studies\n"; + } + else if (!strcmp(transformMode, "manual")) { + scale = par("scale"); + offsetX = par("offsetX"); + offsetY = par("offsetY"); + offsetZ = par("offsetZ"); + } + else + throw cRuntimeError("Unknown transformMode '%s' (expected \"metric\", \"fit\" or \"manual\")", transformMode); + + std::vector xs(cloud.getNumPoints()), ys(cloud.getNumPoints()), zs(cloud.getNumPoints()); + for (int i = 0; i < cloud.getNumPoints(); i++) { + xs[i] = cloud.xs[i] * scale + offsetX; + ys[i] = cloud.ys[i] * scale + offsetY; + zs[i] = cloud.zs[i] * scale + offsetZ; + } + heightfield.buildFromPoints(xs, ys, zs, cellSize, maxCells, 8, par("despikeThreshold")); + + double elapsed = (double)(clock() - startTime) / CLOCKS_PER_SEC; + Coord peak = heightfield.getPeakLocation(); + EV_INFO << "PointCloudGround: loaded " << cloud.getNumPoints() << " points from '" << file + << "', rasterized into a " << heightfield.getNumCellsX() << "x" << heightfield.getNumCellsY() + << " heightfield (cell size " << heightfield.getCellSize() << " m, transformMode " << transformMode + << ", scale " << scale << ") in " << elapsed << " s; extent x=[" << heightfield.getMinX() << ", " << heightfield.getMaxX() + << "] y=[" << heightfield.getMinY() << ", " << heightfield.getMaxY() + << "], peak " << peak.z << " m at (" << peak.x << ", " << peak.y << ")\n"; + } +} + +Coord PointCloudGround::computeGroundProjection(const Coord& position) const +{ + double elevation = heightfield.getElevation(position.x, position.y); + if (std::isnan(elevation)) + elevation = outOfBoundsElevation; + return Coord(position.x, position.y, elevation); +} + +Coord PointCloudGround::computeGroundNormal(const Coord& position) const +{ + return heightfield.getNormal(position.x, position.y); +} + +} // namespace physicalenvironment + +} // namespace inet diff --git a/src/inet/environment/ground/PointCloudGround.h b/src/inet/environment/ground/PointCloudGround.h new file mode 100644 index 00000000000..86e855e1c21 --- /dev/null +++ b/src/inet/environment/ground/PointCloudGround.h @@ -0,0 +1,45 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#ifndef __INET_POINTCLOUDGROUND_H +#define __INET_POINTCLOUDGROUND_H + +#include "inet/common/Module.h" +#include "inet/common/geometry/common/Heightfield.h" +#include "inet/environment/contract/IGround.h" + +namespace inet { + +namespace physicalenvironment { + +/** + * Ground model whose surface is a digital surface model rasterized from a PLY + * point cloud (e.g. a LIDAR scan) — the physics-side counterpart of the VSG + * scene visualizer's sceneModel terrain. See PointCloudGround.ned for the + * parameters (file, cell size, and the PLY-to-simulation-coordinate transform). + */ +class INET_API PointCloudGround : public IGround, public Module +{ + protected: + Heightfield heightfield; + double outOfBoundsElevation = NaN; + + protected: + virtual void initialize(int stage) override; + + public: + virtual Coord computeGroundProjection(const Coord& position) const override; + virtual Coord computeGroundNormal(const Coord& position) const override; + + const Heightfield& getHeightfield() const { return heightfield; } +}; + +} // namespace physicalenvironment + +} // namespace inet + +#endif diff --git a/src/inet/environment/ground/PointCloudGround.ned b/src/inet/environment/ground/PointCloudGround.ned new file mode 100644 index 00000000000..9526ddc2ae6 --- /dev/null +++ b/src/inet/environment/ground/PointCloudGround.ned @@ -0,0 +1,61 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +package inet.environment.ground; + +import inet.common.Module; +import inet.environment.contract.IGround; + +// +// Models the ground surface with a digital surface model rasterized from a +// PLY point cloud, e.g. a LIDAR scan (the same kind of file the VSG scene +// visualizer's sceneModel parameter displays). The points are binned into a +// regular grid keeping the maximum elevation per cell (so building tops and +// canopy are preserved), small gaps are filled from neighboring cells, and +// elevation queries interpolate bilinearly between cells. This makes +// ground-aware radio models (e.g. ~TwoRayGroundReflection) use the real +// terrain elevation under each node. +// +// The transformMode parameter controls how PLY coordinates map to simulation +// coordinates: +// - "metric" (default): translate only — 1 PLY unit = 1 simulation meter, +// the cloud's bounding-box minimum is placed at (originX, originY) and its +// lowest point at baseElevation. Physically exact; use this for studies. +// - "fit": reproduces the VSG scene visualizer's display transform (recenter +// + aspect-preserving fit into the scene footprint), so the physical +// terrain matches what the visualizer shows. Requires fitMinX/fitMinY/ +// fitMinZ/fitMaxX/fitMaxY to be set to the visualizer's scene bounds. +// Note: this scales distances — a warning is logged. +// - "manual": explicit scale and offsets, point' = point * scale + offset. +// +// Elevation queries outside the data extent return outOfBoundsElevation +// (NaN by default, following the ~IGround contract); set it to a number to +// clamp instead. +// +module PointCloudGround extends Module like IGround +{ + parameters: + string file; // PLY point-cloud file (resolved relative to the working directory) + double cellSize @unit(m) = default(0m); // heightfield cell size; 0 = automatic (approximates the mean point spacing) + int maxCells = default(67108864); // guard against accidental huge grids; raise it or coarsen cellSize for very large tiles + double despikeThreshold @unit(m) = default(50m); // clamp isolated single-cell spikes exceeding all neighbors by more than this (stray LIDAR returns); 0 disables + string transformMode = default("metric"); // "metric", "fit" or "manual" (see the module description) + double originX @unit(m) = default(0m); // metric mode: simulation x of the cloud's bounding-box minimum + double originY @unit(m) = default(0m); // metric mode: simulation y of the cloud's bounding-box minimum + double baseElevation @unit(m) = default(0m); // metric mode: simulation z of the cloud's lowest point + double fitMinX @unit(m) = default(nan m); // fit mode: the scene visualizer's sceneMinX + double fitMinY @unit(m) = default(nan m); // fit mode: the scene visualizer's sceneMinY + double fitMinZ @unit(m) = default(nan m); // fit mode: the scene visualizer's sceneMinZ + double fitMaxX @unit(m) = default(nan m); // fit mode: the scene visualizer's sceneMaxX + double fitMaxY @unit(m) = default(nan m); // fit mode: the scene visualizer's sceneMaxY + double offsetX @unit(m) = default(0m); // manual mode: x offset + double offsetY @unit(m) = default(0m); // manual mode: y offset + double offsetZ @unit(m) = default(0m); // manual mode: z offset + double scale = default(1); // manual mode: uniform scale factor + double outOfBoundsElevation @unit(m) = default(nan m); // elevation reported outside the data extent + @class(PointCloudGround); +} From 8aae1f4e7b31e9c0bfef096bb5440224decf5992 Mon Sep 17 00:00:00 2001 From: tabgab Date: Fri, 3 Jul 2026 19:29:13 +0300 Subject: [PATCH 08/17] =?UTF-8?q?physicallayer:=20TerrainObstacleLoss=20?= =?UTF-8?q?=E2=80=94=20minimal=20LIDAR=20terrain=20occlusion=20(LOS)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New IObstacleLoss implementation: a transmission is blocked entirely (loss factor 0) when the terrain surface — the PointCloudGround heightfield — rises above the straight line between transmitter and receiver, and is unaffected otherwise. The ray is sampled at the heightfield cell size (capped at 4096 samples); endpoints are excluded so surface-mounted antennas don't occlude themselves; out-of-extent samples don't block. The ground is resolved lazily since the physical environment may initialize after the radio medium. Composes multiplicatively with any path loss via the analog model — no core changes. This is the binary "los" mode of the planned terrain occlusion support; Fresnel-zone attenuation and knife-edge diffraction remain future work pending design review. Observability: logs per-node-pair line-of-sight transitions (established / lost / re-established / cannot be established), with node names resolved from positions (logLinkEvents parameter, on by default). vsglidar: new TerrainOcclusion config on top of TerrainTwoRay, with the traffic reshaped into a star: each drone pings the gcs (now a UdpEchoApp) twice a second, so live drone<->gcs connections hold a data-link arrow pair that fades ~1.5s after terrain occludes the drone — connectivity at a glance. Verified: 60s same-seed run delivers 597 echo round trips without occlusion vs 279 with it. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 29 +++- .../obstacleloss/TerrainObstacleLoss.cc | 129 ++++++++++++++++++ .../common/obstacleloss/TerrainObstacleLoss.h | 60 ++++++++ .../obstacleloss/TerrainObstacleLoss.ned | 34 +++++ 4 files changed, 249 insertions(+), 3 deletions(-) create mode 100644 src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc create mode 100644 src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h create mode 100644 src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index 655967f7eba..6f5dacbd6a1 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -187,22 +187,45 @@ description = "Two-ray ground reflection over the real LIDAR surface (PointCloud *.drone5.mobility.initialY = 80m *.drone5.mobility.initialZ = 160m -# --- ground control station at the foot of the tower; all drones report to it, so every link -# crosses the built-up center of the map --- +# --- ground control station at the foot of the tower: a UDP ECHO server. Every drone pings it +# twice a second and the echo comes straight back, so each live drone<->gcs connection keeps +# refreshing a data-link arrow in BOTH directions; when terrain occludes a drone, its packets +# stop getting through and the arrow pair fades out within ~1.5s — connectivity at a glance --- *.gcs.mobility.initFromDisplayString = false *.gcs.mobility.initialX = 520m *.gcs.mobility.initialY = 360m *.gcs.mobility.initialZ = 55m # open street south of the tower (~41m grade) + 14m antenna mast, clear of the ~50m low roofs *.gcs.numApps = 1 -*.gcs.app[0].typename = "UdpSink" +*.gcs.app[0].typename = "UdpEchoApp" *.gcs.app[0].localPort = 1000 *.gcs.wlan[*].radio.transmitter.power = 10mW *.gcs.wlan[*].bitrate = 54Mbps *.gcs.wlan[*].mgmt.typename = "Ieee80211MgmtAdhoc" *.gcs.wlan[*].agent.typename = "" +# drones talk ONLY to the gcs: frequent small pings instead of the [General] swarm chatter *.drone1.app[1].destAddresses = "gcs" *.drone1.app[1].destPort = 1000 *.drone2.app[1].destAddresses = "gcs" *.drone3.app[0].destAddresses = "gcs" *.drone4.app[0].destAddresses = "gcs" *.drone5.app[0].destAddresses = "gcs" +*.drone1.app[1].messageLength = 500byte +*.drone2.app[1].messageLength = 500byte +*.drone{3..5}.app[0].messageLength = 500byte +*.drone1.app[1].sendInterval = 0.5s +*.drone2.app[1].sendInterval = 0.5s +*.drone{3..5}.app[0].sendInterval = 0.5s +# persistent-looking arrows: refreshed every ping, gone ~1.5s (sim time) after line of sight is lost +*.visualizer.dataLinkVisualizer.fadeOutMode = "simulationTime" +*.visualizer.dataLinkVisualizer.fadeOutTime = 1.5s + +# -------------------------------------------------------------------------------------------- +# Terrain OCCLUSION on top of the two-ray setup: the same heightfield now also blocks any link +# whose direct ray passes through the terrain — a drone dropping behind the tower (or below +# roof level across town) loses its link to the gcs outright, and the data-link lines vanish. +# Minimal binary line-of-sight model (TerrainObstacleLoss); Fresnel/diffraction are future work. +# -------------------------------------------------------------------------------------------- +[Config TerrainOcclusion] +extends = TerrainTwoRay +description = "LIDAR terrain blocks non-line-of-sight links (TerrainObstacleLoss over PointCloudGround)" +*.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc new file mode 100644 index 00000000000..8817f0ff412 --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc @@ -0,0 +1,129 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h" + +#include + +#include "inet/environment/ground/PointCloudGround.h" + +namespace inet { + +namespace physicallayer { + +Define_Module(TerrainObstacleLoss); + +void TerrainObstacleLoss::initialize(int stage) +{ + if (stage == INITSTAGE_LOCAL) { + physicalEnvironment.reference(this, "physicalEnvironmentModule", true); + sampleStep = par("sampleStep"); + logLinkEvents = par("logLinkEvents"); + } +} + +const cModule *TerrainObstacleLoss::findNodeAt(const Coord& position) const +{ + if (nodes.empty()) { + // collect the network nodes that have a mobility submodule (built lazily, after the network is up) + auto networkModule = getSimulation()->getSystemModule(); + for (cModule::SubmoduleIterator it(networkModule); !it.end(); ++it) { + auto mobilityModule = (*it)->getSubmodule("mobility"); + auto mobility = dynamic_cast(mobilityModule); + if (mobility != nullptr) + nodes.emplace_back(*it, mobility); + } + } + const cModule *best = nullptr; + double bestDistance = 10; // best-effort labeling: ignore matches further than 10m from any node + for (auto& node : nodes) { + double distance = const_cast(node.second)->getCurrentPosition().distance(position); + if (distance < bestDistance) { + bestDistance = distance; + best = node.first; + } + } + return best; +} + +void TerrainObstacleLoss::logLineOfSightChange(const Coord& transmissionPosition, const Coord& receptionPosition, bool blocked) const +{ + auto transmitterNode = findNodeAt(transmissionPosition); + auto receiverNode = findNodeAt(receptionPosition); + if (transmitterNode == nullptr || receiverNode == nullptr) + return; // could not attribute the positions to nodes + auto key = std::make_pair(transmitterNode, receiverNode); + auto it = pairBlocked.find(key); + if (it == pairBlocked.end()) { + if (blocked) + EV_WARN << "Terrain occlusion: NO line of sight between " << transmitterNode->getFullName() << " and " << receiverNode->getFullName() + << " — link cannot be established (terrain blocks the direct ray)\n"; + else + EV_INFO << "Terrain occlusion: line of sight established between " << transmitterNode->getFullName() << " and " << receiverNode->getFullName() << "\n"; + pairBlocked[key] = blocked; + } + else if (it->second != blocked) { + if (blocked) + EV_WARN << "Terrain occlusion: line of sight LOST between " << transmitterNode->getFullName() << " and " << receiverNode->getFullName() + << " — link lost (terrain blocks the direct ray)\n"; + else + EV_INFO << "Terrain occlusion: line of sight RE-ESTABLISHED between " << transmitterNode->getFullName() << " and " << receiverNode->getFullName() << "\n"; + it->second = blocked; + } +} + +const Heightfield *TerrainObstacleLoss::getHeightfield() const +{ + if (heightfield == nullptr) { + auto ground = physicalEnvironment->getGround(); + auto pointCloudGround = dynamic_cast(ground); + if (pointCloudGround == nullptr) + throw cRuntimeError("TerrainObstacleLoss requires the physical environment's ground to be a PointCloudGround (found %s)", + ground == nullptr ? "no ground" : "another ground type"); + heightfield = &pointCloudGround->getHeightfield(); + } + return heightfield; +} + +std::ostream& TerrainObstacleLoss::printToStream(std::ostream& stream, int level, int evFlags) const +{ + stream << "TerrainObstacleLoss"; + if (level <= PRINT_LEVEL_TRACE) + stream << EV_FIELD(sampleStep); + return stream; +} + +double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& transmissionPosition, const Coord& receptionPosition) const +{ + const Heightfield *hf = getHeightfield(); + double dx = receptionPosition.x - transmissionPosition.x; + double dy = receptionPosition.y - transmissionPosition.y; + double dz = receptionPosition.z - transmissionPosition.z; + double distance = std::sqrt(dx * dx + dy * dy); + double step = sampleStep > 0 ? sampleStep : hf->getCellSize(); + int numSamples = std::min(4096, std::max(2, (int)(distance / step) + 1)); + // test interior samples only, so an antenna standing on the surface does not occlude itself + bool blocked = false; + for (int i = 1; i < numSamples - 1; i++) { + double t = (double)i / (numSamples - 1); + double groundZ = hf->getElevation(transmissionPosition.x + dx * t, transmissionPosition.y + dy * t); + if (std::isnan(groundZ)) + continue; // outside the data extent: nothing to block there + double rayZ = transmissionPosition.z + dz * t; + if (groundZ > rayZ) { + blocked = true; // terrain or building blocks the direct ray + break; + } + } + if (logLinkEvents) + logLineOfSightChange(transmissionPosition, receptionPosition, blocked); + return blocked ? 0 : 1; +} + +} // namespace physicallayer + +} // namespace inet diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h new file mode 100644 index 00000000000..33315a2706b --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h @@ -0,0 +1,60 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#ifndef __INET_TERRAINOBSTACLELOSS_H +#define __INET_TERRAINOBSTACLELOSS_H + +#include +#include + +#include "inet/common/Module.h" +#include "inet/common/ModuleRefByPar.h" +#include "inet/common/geometry/common/Heightfield.h" +#include "inet/environment/contract/IPhysicalEnvironment.h" +#include "inet/mobility/contract/IMobility.h" +#include "inet/physicallayer/wireless/common/contract/packetlevel/IObstacleLoss.h" + +namespace inet { + +namespace physicallayer { + +/** + * Minimal terrain occlusion: blocks a transmission entirely (loss factor 0) + * when the terrain surface — as rasterized by the PointCloudGround ground + * module — rises above the straight line between the transmitter and the + * receiver; otherwise the transmission is unaffected (loss factor 1). + * Endpoints are excluded from the test so antennas standing on the surface + * do not occlude themselves. See TerrainObstacleLoss.ned. + */ +class INET_API TerrainObstacleLoss : public Module, public IObstacleLoss +{ + protected: + ModuleRefByPar physicalEnvironment; + mutable const Heightfield *heightfield = nullptr; // resolved lazily: the ground module may initialize after this one + double sampleStep = 0; + bool logLinkEvents = false; + + // line-of-sight bookkeeping for logging: last known blocked/clear state per node pair + mutable std::vector> nodes; // network nodes with a mobility submodule (lazy) + mutable std::map, bool> pairBlocked; + + protected: + virtual void initialize(int stage) override; + const Heightfield *getHeightfield() const; + const cModule *findNodeAt(const Coord& position) const; + void logLineOfSightChange(const Coord& transmissionPosition, const Coord& receptionPosition, bool blocked) const; + + public: + virtual std::ostream& printToStream(std::ostream& stream, int level, int evFlags = 0) const override; + virtual double computeObstacleLoss(Hz frequency, const Coord& transmissionPosition, const Coord& receptionPosition) const override; +}; + +} // namespace physicallayer + +} // namespace inet + +#endif diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned new file mode 100644 index 00000000000..c7893d85c5a --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned @@ -0,0 +1,34 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +package inet.physicallayer.wireless.common.obstacleloss; + +import inet.common.Module; +import inet.physicallayer.wireless.common.contract.packetlevel.IObstacleLoss; + +// +// Minimal terrain occlusion model: a transmission is blocked entirely (loss +// factor 0) when the terrain surface rises above the straight line between +// the transmitter and the receiver, and is unaffected (loss factor 1) +// otherwise. The surface comes from the physical environment's ground +// module, which must be a ~PointCloudGround (e.g. a LIDAR scan, so buildings +// present in the scan occlude links passing behind them). The line is sampled +// at the ground heightfield's cell size by default; endpoints are excluded so +// antennas standing on the surface do not occlude themselves. +// +// This is the binary "los" mode of the planned terrain occlusion support; +// graded Fresnel-zone attenuation and knife-edge diffraction are future work. +// +module TerrainObstacleLoss extends Module like IObstacleLoss +{ + parameters: + string physicalEnvironmentModule = default("physicalEnvironment"); // the physical environment whose ground provides the terrain + double sampleStep @unit(m) = default(0m); // spacing of the line-of-sight samples; 0 = the ground heightfield's cell size + bool logLinkEvents = default(true); // log line-of-sight transitions per node pair (established / lost / re-established / cannot be established) + @class(TerrainObstacleLoss); + @display("i=block/control"); +} From 000b6ca02104f88cc18a2f6d8d8d5233fff650e0 Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 13:34:30 +0300 Subject: [PATCH 09/17] visualizer: osg: include where osg::Node is used omnetpp topic/VSG commit fc5a453385 removed the namespace osg { class Node; } forward declaration from the omnetpp headers (it belonged to the dead refOsgNode/unrefOsgNode chain). These three headers hold an osg::ref_ptr but only included , silently relying on that leaked declaration, and stopped compiling against the updated omnetpp. Include directly, like the sibling OSG visualizer headers already do. Same change as PR #1118 (fix/osg-node-includes); duplicated here so this branch builds standalone against current topic/VSG, deduplicates on merge. Co-Authored-By: Claude Fable 5 --- src/inet/visualizer/osg/common/PacketDropOsgVisualizer.h | 1 + src/inet/visualizer/osg/linklayer/LinkBreakOsgVisualizer.h | 1 + src/inet/visualizer/osg/networklayer/RoutingTableOsgVisualizer.h | 1 + 3 files changed, 3 insertions(+) diff --git a/src/inet/visualizer/osg/common/PacketDropOsgVisualizer.h b/src/inet/visualizer/osg/common/PacketDropOsgVisualizer.h index 1026ba97ebc..d6f99b8203b 100644 --- a/src/inet/visualizer/osg/common/PacketDropOsgVisualizer.h +++ b/src/inet/visualizer/osg/common/PacketDropOsgVisualizer.h @@ -8,6 +8,7 @@ #ifndef __INET_PACKETDROPOSGVISUALIZER_H #define __INET_PACKETDROPOSGVISUALIZER_H +#include #include #include "inet/visualizer/base/PacketDropVisualizerBase.h" diff --git a/src/inet/visualizer/osg/linklayer/LinkBreakOsgVisualizer.h b/src/inet/visualizer/osg/linklayer/LinkBreakOsgVisualizer.h index 4c08d14e018..1b14e48b264 100644 --- a/src/inet/visualizer/osg/linklayer/LinkBreakOsgVisualizer.h +++ b/src/inet/visualizer/osg/linklayer/LinkBreakOsgVisualizer.h @@ -8,6 +8,7 @@ #ifndef __INET_LINKBREAKOSGVISUALIZER_H #define __INET_LINKBREAKOSGVISUALIZER_H +#include #include #include "inet/visualizer/base/LinkBreakVisualizerBase.h" diff --git a/src/inet/visualizer/osg/networklayer/RoutingTableOsgVisualizer.h b/src/inet/visualizer/osg/networklayer/RoutingTableOsgVisualizer.h index ed0192d310a..5b1b0bb49be 100644 --- a/src/inet/visualizer/osg/networklayer/RoutingTableOsgVisualizer.h +++ b/src/inet/visualizer/osg/networklayer/RoutingTableOsgVisualizer.h @@ -8,6 +8,7 @@ #ifndef __INET_ROUTINGTABLEOSGVISUALIZER_H #define __INET_ROUTINGTABLEOSGVISUALIZER_H +#include #include #include "inet/visualizer/base/RoutingTableVisualizerBase.h" From 60ea036ea9291c76463d1cafd69cb03829a10239 Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 14:21:32 +0300 Subject: [PATCH 10/17] physicallayer: TwoRayInterference: optional terrain-aware antenna heights MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The model derived both antenna heights from raw z coordinates, silently assuming flat ground at z=0 — over a terrain ground model (e.g. LIDAR tiles whose streets sit at z≈40-55 m) the heights were off by the full terrain elevation. New optional physicalEnvironmentModule parameter: empty by default (byte-identical legacy Veins behavior); when set, heights are measured above the environment's ground model via computeGroundProjection(), mirroring TwoRayGroundReflection, with a fallback to raw z where the ground is undefined (NaN contract). Co-Authored-By: Claude Fable 5 --- .../common/pathloss/TwoRayInterference.cc | 29 +++++++++++++++++-- .../common/pathloss/TwoRayInterference.h | 3 ++ .../common/pathloss/TwoRayInterference.ned | 4 +++ 3 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.cc b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.cc index a2239f05583..d827c219ea2 100644 --- a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.cc +++ b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.cc @@ -1,6 +1,10 @@ #include "inet/physicallayer/wireless/common/pathloss/TwoRayInterference.h" +#include + +#include "inet/common/ModuleAccess.h" +#include "inet/environment/contract/IGround.h" #include "inet/physicallayer/wireless/common/contract/packetlevel/INarrowbandSignalAnalogModel.h" #include "inet/physicallayer/wireless/common/contract/packetlevel/IRadioMedium.h" @@ -26,6 +30,8 @@ TwoRayInterference::TwoRayInterference() : void TwoRayInterference::initialize(int stage) { if (stage == INITSTAGE_LOCAL) { + // optional: when unset (the default), ground is assumed flat at z=0, as in the original Veins model + physicalEnvironment = findModuleFromPar(par("physicalEnvironmentModule"), this); epsilon_r = par("epsilon_r"); const std::string polarization_str = par("polarization"); if (polarization_str == "horizontal") { @@ -59,10 +65,29 @@ double TwoRayInterference::computePathLoss(const ITransmission *transmission, co return computeTwoRayInterference(transmission->getStartPosition(), arrival->getStartPosition(), waveLength); } +double TwoRayInterference::computeHeightAboveGround(const Coord& position) const +{ + // The model needs antenna heights above the reflecting ground. Without a + // physical environment the ground is assumed flat at z=0 (the original + // Veins behavior); with one, heights are measured above its ground model + // (e.g. terrain), falling back to raw z where the ground is undefined. + if (physicalEnvironment != nullptr) { + auto ground = physicalEnvironment->getGround(); + if (ground != nullptr) { + Coord projection = ground->computeGroundProjection(position); + if (std::isfinite(projection.z)) + return position.distance(projection); + } + } + return position.z; +} + double TwoRayInterference::computeTwoRayInterference(const Coord& pos_t, const Coord& pos_r, m lambda) const { - const double h_sum = pos_t.z + pos_r.z; - const double h_diff = pos_t.z - pos_r.z; + const double h_t = computeHeightAboveGround(pos_t); + const double h_r = computeHeightAboveGround(pos_r); + const double h_sum = h_t + h_r; + const double h_diff = h_t - h_r; // direct line of sight between Tx and Rx antenna const double d_los = pos_r.distance(pos_t); diff --git a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.h b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.h index 48ffb3cfdfb..72568c4710f 100644 --- a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.h +++ b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.h @@ -6,6 +6,7 @@ #define __INET_TWORAYINTERFERENCE_H #include "inet/common/Module.h" +#include "inet/environment/contract/IPhysicalEnvironment.h" #include namespace inet { @@ -32,10 +33,12 @@ class INET_API TwoRayInterference : public Module, public IPathLoss protected: double epsilon_r; char polarization; + const physicalenvironment::IPhysicalEnvironment *physicalEnvironment = nullptr; // optional: antenna heights above its ground model instead of raw z protected: virtual double computeTwoRayInterference(const Coord& posTx, const Coord& posRx, m waveLength) const; virtual double reflectionCoefficient(double cos_theta, double sin_theta) const; + virtual double computeHeightAboveGround(const Coord& position) const; }; } // namespace physicallayer diff --git a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.ned b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.ned index 00191559202..0ea9d18ab49 100644 --- a/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.ned +++ b/src/inet/physicallayer/wireless/common/pathloss/TwoRayInterference.ned @@ -30,5 +30,9 @@ module TwoRayInterference extends Module like IPathLoss @display("i=block/control"); double epsilon_r = default(1.02); string polarization = default("horizontal"); + string physicalEnvironmentModule = default(""); // optional path of a physical environment module; when set, antenna + // heights are measured above its ground model (e.g. terrain) instead + // of assuming flat ground at z=0 (the default, which preserves the + // original Veins-model behavior) } From 07640af3c32c42f659e72bfb2bc0e7db79cff85c Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 14:21:48 +0300 Subject: [PATCH 11/17] physicallayer: TerrainObstacleLoss: graded fresnel mode + obstacle-loss tracing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mode="fresnel": instead of the binary line-of-sight cliff, the worst first-Fresnel-zone intrusion along the terrain profile is mapped through the ITU-R P.526 single knife-edge curve J(v) — no loss with >~0.55 r1 clearance, ~6 dB with terrain right at the direct ray, deepening smoothly (and frequency-dependently) into the shadow. J(v) is a public static for unit testing (tests/unit/TerrainObstacleLoss_1.test). mode="los" (the default) is behavior-identical to before. The module now extends TracingObstacleLossBase and emits obstaclePenetrated events at the critical terrain point, so the existing tracing obstacle loss visualizers display where links are obstructed. Terrain has no associated physical object: the event carries nullptr and world-frame coordinates (documented in ITracingObstacleLoss.h), and the canvas/osg/vsg renderers now handle that with identity rotation and no offset. vsglidar demos: TerrainFresnel (graded loss + intersection markers), intersection markers on TerrainOcclusion, and TerrainTwoRayInterference (the terrain-aware TwoRayInterference over the same scene). Headless 60s/seed-1 check: TerrainTwoRay and TerrainOcclusion app-level counts unchanged (597/279); fresnel differs from binary occlusion in 342 radio-level scalars while app counts coincide — at sharp roof edges the graded transition shell is thin, so end-to-end ping outcomes quantize the same way. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 33 +++++++- .../packetlevel/ITracingObstacleLoss.h | 6 ++ .../obstacleloss/TerrainObstacleLoss.cc | 80 +++++++++++++++++-- .../common/obstacleloss/TerrainObstacleLoss.h | 31 ++++--- .../obstacleloss/TerrainObstacleLoss.ned | 40 ++++++---- .../TracingObstacleLossCanvasVisualizer.cc | 6 +- .../TracingObstacleLossOsgVisualizer.cc | 6 +- .../TracingObstacleLossVsgVisualizer.cc | 6 +- tests/unit/TerrainObstacleLoss_1.test | 57 +++++++++++++ 9 files changed, 226 insertions(+), 39 deletions(-) create mode 100644 tests/unit/TerrainObstacleLoss_1.test diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index 6f5dacbd6a1..49ebbc03f67 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -223,9 +223,40 @@ description = "Two-ray ground reflection over the real LIDAR surface (PointCloud # Terrain OCCLUSION on top of the two-ray setup: the same heightfield now also blocks any link # whose direct ray passes through the terrain — a drone dropping behind the tower (or below # roof level across town) loses its link to the gcs outright, and the data-link lines vanish. -# Minimal binary line-of-sight model (TerrainObstacleLoss); Fresnel/diffraction are future work. +# Binary line-of-sight mode of TerrainObstacleLoss; see TerrainFresnel for the graded model. +# The tracing obstacle loss visualizer marks the exact terrain point that blocks each link. # -------------------------------------------------------------------------------------------- [Config TerrainOcclusion] extends = TerrainTwoRay description = "LIDAR terrain blocks non-line-of-sight links (TerrainObstacleLoss over PointCloudGround)" *.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" +*.visualizer.obstacleLossVisualizer.displayIntersections = true +*.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" +*.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s + +# -------------------------------------------------------------------------------------------- +# Graded terrain loss: instead of the binary cliff, the worst first-Fresnel-zone clearance on +# each link is mapped through the ITU-R P.526 single knife-edge curve — links weaken smoothly +# (~6 dB with terrain right at the ray) as a drone sinks behind a building, then fade out in +# the deepening shadow instead of dying at an invisible boundary. Frequency-dependent. +# -------------------------------------------------------------------------------------------- +[Config TerrainFresnel] +extends = TerrainTwoRay +description = "graded Fresnel/knife-edge terrain attenuation (TerrainObstacleLoss mode=fresnel)" +*.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" +*.radioMedium.obstacleLoss.mode = "fresnel" +*.visualizer.obstacleLossVisualizer.displayIntersections = true +*.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" +*.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s + +# -------------------------------------------------------------------------------------------- +# Phase-correct two-ray path loss over the real terrain: TwoRayInterference normally assumes +# flat ground at z=0 (its antenna "heights" would be ~340 m here, since the Sandy tile's +# streets sit at z≈41-55 m); pointing it at the physical environment makes it measure heights +# above the LIDAR ground model instead. +# -------------------------------------------------------------------------------------------- +[Config TerrainTwoRayInterference] +extends = TerrainTwoRay +description = "TwoRayInterference with antenna heights above the LIDAR terrain" +*.radioMedium.pathLoss.typename = "TwoRayInterference" +*.radioMedium.pathLoss.physicalEnvironmentModule = "physicalEnvironment" diff --git a/src/inet/physicallayer/wireless/common/contract/packetlevel/ITracingObstacleLoss.h b/src/inet/physicallayer/wireless/common/contract/packetlevel/ITracingObstacleLoss.h index 4a07b8ba248..a2d360ea545 100644 --- a/src/inet/physicallayer/wireless/common/contract/packetlevel/ITracingObstacleLoss.h +++ b/src/inet/physicallayer/wireless/common/contract/packetlevel/ITracingObstacleLoss.h @@ -20,6 +20,12 @@ class INET_API ITracingObstacleLoss : public IObstacleLoss public: class INET_API ObstaclePenetratedEvent : public cObject { public: + /** + * The penetrated physical object, or nullptr for obstructions that + * are not physical objects (e.g. terrain). When an object is present, + * the intersection points and normals below are in the object's local + * coordinate frame; when it is nullptr, they are world coordinates. + */ const physicalenvironment::IPhysicalObject *object; const Coord intersection1; const Coord intersection2; diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc index 8817f0ff412..d454231bcd8 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc @@ -8,6 +8,7 @@ #include "inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h" #include +#include #include "inet/environment/ground/PointCloudGround.h" @@ -21,11 +22,38 @@ void TerrainObstacleLoss::initialize(int stage) { if (stage == INITSTAGE_LOCAL) { physicalEnvironment.reference(this, "physicalEnvironmentModule", true); + const char *modeString = par("mode"); + if (!strcmp(modeString, "los")) + mode = LOS; + else if (!strcmp(modeString, "fresnel")) + mode = FRESNEL; + else + throw cRuntimeError("Unknown mode '%s'", modeString); sampleStep = par("sampleStep"); logLinkEvents = par("logLinkEvents"); } } +double TerrainObstacleLoss::computeKnifeEdgeLoss(double v) +{ + // ITU-R P.526 approximate single knife-edge diffraction loss, valid for + // v > -0.78; below that the obstruction is far enough outside the first + // Fresnel zone to have no effect (the curve is ~0 dB at the cutoff). + if (v <= -0.78) + return 0; + return 6.9 + 20 * std::log10(std::sqrt((v - 0.1) * (v - 0.1) + 1) + v - 0.1); +} + +void TerrainObstacleLoss::emitObstaclePenetrated(const Coord& rayPoint, const Coord& groundPoint, double lossFactor) const +{ + // No physical object is associated with terrain: coordinates are world-frame + // (see ObstaclePenetratedEvent). The intersection segment connects the direct + // ray to the terrain surface at the critical point, making the obstruction + // (or near-grazing pinch) visible in the tracing obstacle loss visualizers. + ObstaclePenetratedEvent event(nullptr, rayPoint, groundPoint, getHeightfield()->getNormal(groundPoint.x, groundPoint.y), Coord(0, 0, 1), lossFactor); + const_cast(this)->emit(obstaclePenetratedSignal, &event); +} + const cModule *TerrainObstacleLoss::findNodeAt(const Coord& position) const { if (nodes.empty()) { @@ -93,7 +121,7 @@ std::ostream& TerrainObstacleLoss::printToStream(std::ostream& stream, int level { stream << "TerrainObstacleLoss"; if (level <= PRINT_LEVEL_TRACE) - stream << EV_FIELD(sampleStep); + stream << EV_FIELD(mode, mode == LOS ? "los" : "fresnel") << EV_FIELD(sampleStep); return stream; } @@ -103,25 +131,61 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans double dx = receptionPosition.x - transmissionPosition.x; double dy = receptionPosition.y - transmissionPosition.y; double dz = receptionPosition.z - transmissionPosition.z; - double distance = std::sqrt(dx * dx + dy * dy); + double horizontalDistance = std::sqrt(dx * dx + dy * dy); + double totalDistance = transmissionPosition.distance(receptionPosition); double step = sampleStep > 0 ? sampleStep : hf->getCellSize(); - int numSamples = std::min(4096, std::max(2, (int)(distance / step) + 1)); - // test interior samples only, so an antenna standing on the surface does not occlude itself + int numSamples = std::min(4096, std::max(2, (int)(horizontalDistance / step) + 1)); + double waveLength = SPEED_OF_LIGHT / frequency.get(); + // Walk interior samples only, so an antenna standing on the surface does not + // occlude itself. LOS mode stops at the first obstruction; FRESNEL mode scans + // the whole profile for the worst diffraction parameter v (the deepest + // intrusion into the first Fresnel zone, in zone-normalized units), which the + // single knife-edge curve J(v) then maps to a graded loss. bool blocked = false; + double worstV = -INFINITY; + double worstT = NaN; + double worstGroundZ = NaN; for (int i = 1; i < numSamples - 1; i++) { double t = (double)i / (numSamples - 1); double groundZ = hf->getElevation(transmissionPosition.x + dx * t, transmissionPosition.y + dy * t); if (std::isnan(groundZ)) continue; // outside the data extent: nothing to block there double rayZ = transmissionPosition.z + dz * t; - if (groundZ > rayZ) { - blocked = true; // terrain or building blocks the direct ray - break; + if (mode == LOS) { + if (groundZ > rayZ) { + blocked = true; // terrain or building blocks the direct ray + worstT = t; + worstGroundZ = groundZ; + break; + } + } + else { + double d1 = totalDistance * t; + double d2 = totalDistance - d1; + double v = (groundZ - rayZ) * std::sqrt(2 * totalDistance / (waveLength * d1 * d2)); + if (v > worstV) { + worstV = v; + worstT = t; + worstGroundZ = groundZ; + } } } + double lossFactor; + if (mode == LOS) + lossFactor = blocked ? 0 : 1; + else { + double lossDb = std::isfinite(worstV) ? computeKnifeEdgeLoss(worstV) : 0; + lossFactor = std::pow(10.0, -lossDb / 10.0); + blocked = worstV > 0; // the direct ray itself is geometrically obstructed + } + if (lossFactor < 1) { + Coord rayPoint(transmissionPosition.x + dx * worstT, transmissionPosition.y + dy * worstT, transmissionPosition.z + dz * worstT); + Coord groundPoint(rayPoint.x, rayPoint.y, worstGroundZ); + emitObstaclePenetrated(rayPoint, groundPoint, lossFactor); + } if (logLinkEvents) logLineOfSightChange(transmissionPosition, receptionPosition, blocked); - return blocked ? 0 : 1; + return lossFactor; } } // namespace physicallayer diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h index 33315a2706b..6c108cc1de4 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h @@ -11,30 +11,35 @@ #include #include -#include "inet/common/Module.h" #include "inet/common/ModuleRefByPar.h" #include "inet/common/geometry/common/Heightfield.h" #include "inet/environment/contract/IPhysicalEnvironment.h" #include "inet/mobility/contract/IMobility.h" -#include "inet/physicallayer/wireless/common/contract/packetlevel/IObstacleLoss.h" +#include "inet/physicallayer/wireless/common/base/packetlevel/TracingObstacleLossBase.h" namespace inet { namespace physicallayer { /** - * Minimal terrain occlusion: blocks a transmission entirely (loss factor 0) - * when the terrain surface — as rasterized by the PointCloudGround ground - * module — rises above the straight line between the transmitter and the - * receiver; otherwise the transmission is unaffected (loss factor 1). - * Endpoints are excluded from the test so antennas standing on the surface - * do not occlude themselves. See TerrainObstacleLoss.ned. + * Terrain occlusion from the PointCloudGround heightfield. In "los" mode a + * transmission is blocked entirely (loss factor 0) when the terrain rises + * above the straight line between transmitter and receiver. In "fresnel" mode + * the loss is graded: the worst first-Fresnel-zone clearance along the path is + * mapped through the ITU-R P.526 single knife-edge curve J(v), so links + * degrade smoothly as terrain approaches and then crosses the direct ray. + * Emits ITracingObstacleLoss::obstaclePenetratedSignal at the critical + * terrain point (with no associated physical object, world coordinates). + * See TerrainObstacleLoss.ned. */ -class INET_API TerrainObstacleLoss : public Module, public IObstacleLoss +class INET_API TerrainObstacleLoss : public TracingObstacleLossBase { protected: + enum Mode { LOS, FRESNEL }; + ModuleRefByPar physicalEnvironment; mutable const Heightfield *heightfield = nullptr; // resolved lazily: the ground module may initialize after this one + Mode mode = LOS; double sampleStep = 0; bool logLinkEvents = false; @@ -47,8 +52,16 @@ class INET_API TerrainObstacleLoss : public Module, public IObstacleLoss const Heightfield *getHeightfield() const; const cModule *findNodeAt(const Coord& position) const; void logLineOfSightChange(const Coord& transmissionPosition, const Coord& receptionPosition, bool blocked) const; + void emitObstaclePenetrated(const Coord& rayPoint, const Coord& groundPoint, double lossFactor) const; public: + /** + * Single knife-edge diffraction loss J(v) in dB from ITU-R P.526 for the + * dimensionless diffraction parameter v; 0 dB for v <= -0.78 (the curve + * is continuous at the cutoff). + */ + static double computeKnifeEdgeLoss(double v); + virtual std::ostream& printToStream(std::ostream& stream, int level, int evFlags = 0) const override; virtual double computeObstacleLoss(Hz frequency, const Coord& transmissionPosition, const Coord& receptionPosition) const override; }; diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned index c7893d85c5a..7ebc477f10a 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned @@ -7,28 +7,38 @@ package inet.physicallayer.wireless.common.obstacleloss; -import inet.common.Module; -import inet.physicallayer.wireless.common.contract.packetlevel.IObstacleLoss; +import inet.physicallayer.wireless.common.base.packetlevel.TracingObstacleLossBase; // -// Minimal terrain occlusion model: a transmission is blocked entirely (loss -// factor 0) when the terrain surface rises above the straight line between -// the transmitter and the receiver, and is unaffected (loss factor 1) -// otherwise. The surface comes from the physical environment's ground -// module, which must be a ~PointCloudGround (e.g. a LIDAR scan, so buildings -// present in the scan occlude links passing behind them). The line is sampled -// at the ground heightfield's cell size by default; endpoints are excluded so -// antennas standing on the surface do not occlude themselves. +// Terrain occlusion model driven by the physical environment's ground module, +// which must be a ~PointCloudGround (e.g. a LIDAR scan, so buildings present +// in the scan occlude links passing behind them). The straight line between +// the transmitter and the receiver is sampled at the ground heightfield's +// cell size by default; endpoints are excluded so antennas standing on the +// surface do not occlude themselves. // -// This is the binary "los" mode of the planned terrain occlusion support; -// graded Fresnel-zone attenuation and knife-edge diffraction are future work. +// Two modes: +// - "los" (default): binary — the transmission is blocked entirely (loss +// factor 0) when the terrain surface rises above the direct line, and is +// unaffected (loss factor 1) otherwise. +// - "fresnel": graded — the worst first-Fresnel-zone clearance along the path +// is mapped through the ITU-R P.526 single knife-edge diffraction curve +// J(v): no loss while the zone is clear by more than ~0.55 of its radius, +// about 6 dB with the terrain right at the direct ray, deepening smoothly +// (and frequency-dependently) as the ray sinks below the terrain. Links +// thus degrade progressively instead of dropping at an invisible boundary. // -module TerrainObstacleLoss extends Module like IObstacleLoss +// Emits obstaclePenetrated events at the critical terrain point (with no +// associated physical object, world coordinates), so tracing obstacle loss +// visualizers can display where links are obstructed. +// +// Multi-edge (Deygout) diffraction is future work. +// +module TerrainObstacleLoss extends TracingObstacleLossBase { parameters: - string physicalEnvironmentModule = default("physicalEnvironment"); // the physical environment whose ground provides the terrain + string mode @enum("los", "fresnel") = default("los"); // binary line-of-sight blocking, or graded Fresnel/knife-edge attenuation double sampleStep @unit(m) = default(0m); // spacing of the line-of-sight samples; 0 = the ground heightfield's cell size bool logLinkEvents = default(true); // log line-of-sight transitions per node pair (established / lost / re-established / cannot be established) @class(TerrainObstacleLoss); - @display("i=block/control"); } diff --git a/src/inet/visualizer/canvas/physicallayer/TracingObstacleLossCanvasVisualizer.cc b/src/inet/visualizer/canvas/physicallayer/TracingObstacleLossCanvasVisualizer.cc index 5ac4508328a..28c95580620 100644 --- a/src/inet/visualizer/canvas/physicallayer/TracingObstacleLossCanvasVisualizer.cc +++ b/src/inet/visualizer/canvas/physicallayer/TracingObstacleLossCanvasVisualizer.cc @@ -59,8 +59,10 @@ const TracingObstacleLossVisualizerBase::ObstacleLossVisualization *TracingObsta auto normal1 = obstaclePenetratedEvent->normal1; auto normal2 = obstaclePenetratedEvent->normal2; auto loss = obstaclePenetratedEvent->loss; - const RotationMatrix rotation(object->getOrientation().toEulerAngles()); - const Coord& position = object->getPosition(); + // A nullptr object (e.g. terrain from TerrainObstacleLoss) means the intersection + // points and normals are already world-frame: identity rotation, no offset. + const RotationMatrix rotation = object != nullptr ? RotationMatrix(object->getOrientation().toEulerAngles()) : RotationMatrix(); + const Coord position = object != nullptr ? object->getPosition() : Coord::ZERO; const Coord rotatedIntersection1 = rotation.rotateVector(intersection1); const Coord rotatedIntersection2 = rotation.rotateVector(intersection2); double intersectionDistance = intersection2.distance(intersection1); diff --git a/src/inet/visualizer/osg/physicallayer/TracingObstacleLossOsgVisualizer.cc b/src/inet/visualizer/osg/physicallayer/TracingObstacleLossOsgVisualizer.cc index 8ea0d5e7559..848e889a4f6 100644 --- a/src/inet/visualizer/osg/physicallayer/TracingObstacleLossOsgVisualizer.cc +++ b/src/inet/visualizer/osg/physicallayer/TracingObstacleLossOsgVisualizer.cc @@ -49,8 +49,10 @@ const TracingObstacleLossVisualizerBase::ObstacleLossVisualization *TracingObsta auto normal1 = obstaclePenetratedEvent->normal1; auto normal2 = obstaclePenetratedEvent->normal2; // TODO display auto loss = obstaclePenetratedEvent->loss; - const RotationMatrix rotation(object->getOrientation().toEulerAngles()); - const Coord& position = object->getPosition(); + // A nullptr object (e.g. terrain from TerrainObstacleLoss) means the intersection + // points and normals are already world-frame: identity rotation, no offset. + const RotationMatrix rotation = object != nullptr ? RotationMatrix(object->getOrientation().toEulerAngles()) : RotationMatrix(); + const Coord position = object != nullptr ? object->getPosition() : Coord::ZERO; const Coord rotatedIntersection1 = rotation.rotateVector(intersection1); const Coord rotatedIntersection2 = rotation.rotateVector(intersection2); double intersectionDistance = intersection2.distance(intersection1); diff --git a/src/inet/visualizer/vsg/physicallayer/TracingObstacleLossVsgVisualizer.cc b/src/inet/visualizer/vsg/physicallayer/TracingObstacleLossVsgVisualizer.cc index 6661695886f..93b58bf7602 100644 --- a/src/inet/visualizer/vsg/physicallayer/TracingObstacleLossVsgVisualizer.cc +++ b/src/inet/visualizer/vsg/physicallayer/TracingObstacleLossVsgVisualizer.cc @@ -50,8 +50,10 @@ TracingObstacleLossVsgVisualizer::createObstacleLossVisualization(const ITracing auto normal2 = obstaclePenetratedEvent->normal2; // TODO: display obstaclePenetratedEvent->loss as a label or color - const RotationMatrix rotation(object->getOrientation().toEulerAngles()); - const Coord& position = object->getPosition(); + // A nullptr object (e.g. terrain from TerrainObstacleLoss) means the intersection + // points and normals are already world-frame: identity rotation, no offset. + const RotationMatrix rotation = object != nullptr ? RotationMatrix(object->getOrientation().toEulerAngles()) : RotationMatrix(); + const Coord position = object != nullptr ? object->getPosition() : Coord::ZERO; const Coord rotatedIntersection1 = rotation.rotateVector(intersection1) + position; const Coord rotatedIntersection2 = rotation.rotateVector(intersection2) + position; diff --git a/tests/unit/TerrainObstacleLoss_1.test b/tests/unit/TerrainObstacleLoss_1.test new file mode 100644 index 00000000000..518dac50613 --- /dev/null +++ b/tests/unit/TerrainObstacleLoss_1.test @@ -0,0 +1,57 @@ +%description: +Test TerrainObstacleLoss::computeKnifeEdgeLoss: the ITU-R P.526 single +knife-edge diffraction curve J(v). Checks the no-loss region below the +v = -0.78 cutoff, continuity at the cutoff, the ~6 dB grazing value at v = 0, +reference points in the shadow region, and monotonicity. + +%includes: +#include +#include "inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h" + +using namespace inet; +using namespace inet::physicallayer; + +%global: + +#define CHECK(label, cond) EV << label << ": " << ((cond) ? "OK" : "FAIL") << "\n"; + +%activity: + +// --- no loss while the obstruction stays outside ~0.55 of the first Fresnel zone --- +CHECK("clear v=-2", TerrainObstacleLoss::computeKnifeEdgeLoss(-2) == 0); +CHECK("clear v=-0.9", TerrainObstacleLoss::computeKnifeEdgeLoss(-0.9) == 0); + +// --- the curve is continuous at the v=-0.78 cutoff (~0 dB just above it) --- +CHECK("continuity at cutoff", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(-0.78 + 1e-6)) < 0.02); + +// --- grazing incidence: terrain exactly at the direct ray costs ~6 dB --- +CHECK("grazing v=0", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(0) - 6.03) < 0.05); + +// --- reference points in the shadow region --- +CHECK("shadow v=1", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(1) - 13.93) < 0.1); +CHECK("shadow v=2", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(2) - 19.04) < 0.1); + +// --- strictly monotone increasing through the transition and shadow regions --- +{ + bool monotone = true; + double previous = TerrainObstacleLoss::computeKnifeEdgeLoss(-0.7); + for (double v = -0.6; v <= 4.0 + 1e-9; v += 0.1) { + double loss = TerrainObstacleLoss::computeKnifeEdgeLoss(v); + if (loss <= previous) + monotone = false; + previous = loss; + } + CHECK("monotone", monotone); +} + +EV << "done\n"; + +%contains: stdout +clear v=-2: OK +clear v=-0.9: OK +continuity at cutoff: OK +grazing v=0: OK +shadow v=1: OK +shadow v=2: OK +monotone: OK +done From 06b5d85f6cb13c0babe38de35744f64bfaf3c804 Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 14:46:42 +0300 Subject: [PATCH 12/17] physicallayer: TerrainObstacleLoss: selectable obstruction marker geometry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New markerStyle parameter for the emitted obstaclePenetrated events: - "depth" (default): vertical segment from the direct ray up to the terrain surface at the worst point — its length shows how deeply the link is buried below the skyline. - "ray": the chord of the direct ray below the terrain surface (the contiguous obstructed sample run around the worst point), matching the look of physical-object obstacle loss models; near-grazing links in fresnel mode draw a one-sample tick along the ray at the pinch point. The geometry is chosen by the emitter, so all tracing visualizers (canvas/OSG/VSG) support both styles unchanged. Physics untouched: same-seed 60 s TerrainFresnel/TerrainOcclusion app counts identical to before. Demo configs show one of each style. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 6 ++ .../obstacleloss/TerrainObstacleLoss.cc | 63 +++++++++++++++---- .../common/obstacleloss/TerrainObstacleLoss.h | 4 +- .../obstacleloss/TerrainObstacleLoss.ned | 11 +++- 4 files changed, 70 insertions(+), 14 deletions(-) diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index 49ebbc03f67..f43cf39e758 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -230,6 +230,8 @@ description = "Two-ray ground reflection over the real LIDAR surface (PointCloud extends = TerrainTwoRay description = "LIDAR terrain blocks non-line-of-sight links (TerrainObstacleLoss over PointCloudGround)" *.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" +# "ray" markers: red chords along the signal path where it passes through buildings +*.radioMedium.obstacleLoss.markerStyle = "ray" *.visualizer.obstacleLossVisualizer.displayIntersections = true *.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" *.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s @@ -245,6 +247,10 @@ extends = TerrainTwoRay description = "graded Fresnel/knife-edge terrain attenuation (TerrainObstacleLoss mode=fresnel)" *.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" *.radioMedium.obstacleLoss.mode = "fresnel" +# "depth" markers (the default): vertical gauges from the direct ray up to the blocking +# surface — their length shows how deeply each link is buried. Set to "ray" for chords +# along the signal path instead (see TerrainOcclusion). +*.radioMedium.obstacleLoss.markerStyle = "depth" *.visualizer.obstacleLossVisualizer.displayIntersections = true *.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" *.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc index d454231bcd8..71f7a727c4b 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc @@ -29,6 +29,13 @@ void TerrainObstacleLoss::initialize(int stage) mode = FRESNEL; else throw cRuntimeError("Unknown mode '%s'", modeString); + const char *markerStyleString = par("markerStyle"); + if (!strcmp(markerStyleString, "depth")) + markerStyle = DEPTH; + else if (!strcmp(markerStyleString, "ray")) + markerStyle = RAY; + else + throw cRuntimeError("Unknown markerStyle '%s'", markerStyleString); sampleStep = par("sampleStep"); logLinkEvents = par("logLinkEvents"); } @@ -44,13 +51,14 @@ double TerrainObstacleLoss::computeKnifeEdgeLoss(double v) return 6.9 + 20 * std::log10(std::sqrt((v - 0.1) * (v - 0.1) + 1) + v - 0.1); } -void TerrainObstacleLoss::emitObstaclePenetrated(const Coord& rayPoint, const Coord& groundPoint, double lossFactor) const +void TerrainObstacleLoss::emitObstaclePenetrated(const Coord& intersection1, const Coord& intersection2, double lossFactor) const { // No physical object is associated with terrain: coordinates are world-frame - // (see ObstaclePenetratedEvent). The intersection segment connects the direct - // ray to the terrain surface at the critical point, making the obstruction - // (or near-grazing pinch) visible in the tracing obstacle loss visualizers. - ObstaclePenetratedEvent event(nullptr, rayPoint, groundPoint, getHeightfield()->getNormal(groundPoint.x, groundPoint.y), Coord(0, 0, 1), lossFactor); + // (see ObstaclePenetratedEvent). The intersection segment geometry depends on + // markerStyle — see the NED documentation. + const Heightfield *hf = getHeightfield(); + ObstaclePenetratedEvent event(nullptr, intersection1, intersection2, + hf->getNormal(intersection1.x, intersection1.y), hf->getNormal(intersection2.x, intersection2.y), lossFactor); const_cast(this)->emit(obstaclePenetratedSignal, &event); } @@ -143,7 +151,7 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans // single knife-edge curve J(v) then maps to a graded loss. bool blocked = false; double worstV = -INFINITY; - double worstT = NaN; + int worstIndex = -1; double worstGroundZ = NaN; for (int i = 1; i < numSamples - 1; i++) { double t = (double)i / (numSamples - 1); @@ -154,7 +162,7 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans if (mode == LOS) { if (groundZ > rayZ) { blocked = true; // terrain or building blocks the direct ray - worstT = t; + worstIndex = i; worstGroundZ = groundZ; break; } @@ -165,7 +173,7 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans double v = (groundZ - rayZ) * std::sqrt(2 * totalDistance / (waveLength * d1 * d2)); if (v > worstV) { worstV = v; - worstT = t; + worstIndex = i; worstGroundZ = groundZ; } } @@ -178,10 +186,41 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans lossFactor = std::pow(10.0, -lossDb / 10.0); blocked = worstV > 0; // the direct ray itself is geometrically obstructed } - if (lossFactor < 1) { - Coord rayPoint(transmissionPosition.x + dx * worstT, transmissionPosition.y + dy * worstT, transmissionPosition.z + dz * worstT); - Coord groundPoint(rayPoint.x, rayPoint.y, worstGroundZ); - emitObstaclePenetrated(rayPoint, groundPoint, lossFactor); + if (lossFactor < 1 && worstIndex > 0) { + auto rayPointAt = [&](double tt) { + return Coord(transmissionPosition.x + dx * tt, transmissionPosition.y + dy * tt, transmissionPosition.z + dz * tt); + }; + Coord intersection1, intersection2; + if (markerStyle == RAY) { + if (blocked) { + // the chord of the direct ray below the terrain surface: the contiguous + // obstructed run of samples around the worst point + auto isObstructedAt = [&](int i) { + Coord p = rayPointAt((double)i / (numSamples - 1)); + double groundZ = hf->getElevation(p.x, p.y); + return !std::isnan(groundZ) && groundZ > p.z; + }; + int first = worstIndex, last = worstIndex; + while (first > 1 && isObstructedAt(first - 1)) + first--; + while (last < numSamples - 2 && isObstructedAt(last + 1)) + last++; + intersection1 = rayPointAt((double)first / (numSamples - 1)); + intersection2 = rayPointAt((double)last / (numSamples - 1)); + } + else { + // near-grazing attenuation without geometric obstruction (fresnel mode): + // a one-sample-long tick along the ray at the pinch point + intersection1 = rayPointAt((worstIndex - 0.5) / (numSamples - 1)); + intersection2 = rayPointAt((worstIndex + 0.5) / (numSamples - 1)); + } + } + else { // DEPTH: vertical gauge from the direct ray up to the terrain surface + Coord rayPoint = rayPointAt((double)worstIndex / (numSamples - 1)); + intersection1 = rayPoint; + intersection2 = Coord(rayPoint.x, rayPoint.y, worstGroundZ); + } + emitObstaclePenetrated(intersection1, intersection2, lossFactor); } if (logLinkEvents) logLineOfSightChange(transmissionPosition, receptionPosition, blocked); diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h index 6c108cc1de4..06ee21109fc 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h @@ -36,10 +36,12 @@ class INET_API TerrainObstacleLoss : public TracingObstacleLossBase { protected: enum Mode { LOS, FRESNEL }; + enum MarkerStyle { DEPTH, RAY }; ModuleRefByPar physicalEnvironment; mutable const Heightfield *heightfield = nullptr; // resolved lazily: the ground module may initialize after this one Mode mode = LOS; + MarkerStyle markerStyle = DEPTH; double sampleStep = 0; bool logLinkEvents = false; @@ -52,7 +54,7 @@ class INET_API TerrainObstacleLoss : public TracingObstacleLossBase const Heightfield *getHeightfield() const; const cModule *findNodeAt(const Coord& position) const; void logLineOfSightChange(const Coord& transmissionPosition, const Coord& receptionPosition, bool blocked) const; - void emitObstaclePenetrated(const Coord& rayPoint, const Coord& groundPoint, double lossFactor) const; + void emitObstaclePenetrated(const Coord& intersection1, const Coord& intersection2, double lossFactor) const; public: /** diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned index 7ebc477f10a..c47d67321eb 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned @@ -30,7 +30,15 @@ import inet.physicallayer.wireless.common.base.packetlevel.TracingObstacleLossBa // // Emits obstaclePenetrated events at the critical terrain point (with no // associated physical object, world coordinates), so tracing obstacle loss -// visualizers can display where links are obstructed. +// visualizers can display where links are obstructed. The marker geometry is +// selectable via markerStyle: +// - "depth" (default): a vertical segment from the direct ray up to the +// terrain surface at the worst point — its length shows how deeply the +// link is buried below the skyline. +// - "ray": the chord of the direct ray below the terrain surface (entry to +// exit around the worst point), matching the look of physical-object +// obstacle loss models; near-grazing links (fresnel mode, zone intruded +// but ray clear) draw a short tick along the ray at the pinch point. // // Multi-edge (Deygout) diffraction is future work. // @@ -38,6 +46,7 @@ module TerrainObstacleLoss extends TracingObstacleLossBase { parameters: string mode @enum("los", "fresnel") = default("los"); // binary line-of-sight blocking, or graded Fresnel/knife-edge attenuation + string markerStyle @enum("depth", "ray") = default("depth"); // geometry of the emitted obstruction markers (see module description) double sampleStep @unit(m) = default(0m); // spacing of the line-of-sight samples; 0 = the ground heightfield's cell size bool logLinkEvents = default(true); // log line-of-sight transitions per node pair (established / lost / re-established / cannot be established) @class(TerrainObstacleLoss); From 8cf3be745b5fb8c7a86b81052b799092d4f9db2d Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 15:42:26 +0300 Subject: [PATCH 13/17] physicallayer: TerrainObstacleLoss: multi-edge Deygout diffraction mode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mode="diffraction": extends the single knife-edge "fresnel" mode to several ridges via the Deygout construction — charge the dominant edge along the profile, then recurse into the sub-paths on either side of it, summing J(v) over up to maxDiffractionEdges (default 3) edges. Behind a row of buildings this predicts a deeper, more realistic shadow than the single dominant edge, while still fading rather than cutting off like "los". The dominant edge reuses the geometry the fresnel scan already computes, so fresnel is exactly Deygout truncated to one edge. computeDeygoutLoss is a static free-standing routine, unit-tested in tests/unit/TerrainObstacleLoss_1.test (single-edge equivalence to J(v), clear path, second-ridge superadditivity, edge-count cap). Demo config TerrainDiffraction. Verified headless: los/fresnel/two-ray regressions unchanged (597/279); diffraction differs from fresnel in 196 radio-level scalars (extra loss on multi-ridge paths) while app-level echo counts coincide on this single-tower tile. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 16 +++++ .../obstacleloss/TerrainObstacleLoss.cc | 69 ++++++++++++++++--- .../common/obstacleloss/TerrainObstacleLoss.h | 12 +++- .../obstacleloss/TerrainObstacleLoss.ned | 8 ++- tests/unit/TerrainObstacleLoss_1.test | 47 ++++++++++++- 5 files changed, 141 insertions(+), 11 deletions(-) diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index f43cf39e758..faef55d216a 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -255,6 +255,22 @@ description = "graded Fresnel/knife-edge terrain attenuation (TerrainObstacleLos *.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" *.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s +# -------------------------------------------------------------------------------------------- +# Multi-edge diffraction: the Deygout construction sums the knife-edge loss over the dominant +# ridge on the path plus the secondary ridges on either side of it (up to maxDiffractionEdges). +# Behind a cluster of buildings this predicts a deeper, more realistic shadow than the single +# edge of "fresnel", while still fading rather than cutting off like "los". +# -------------------------------------------------------------------------------------------- +[Config TerrainDiffraction] +extends = TerrainTwoRay +description = "multi-edge Deygout terrain diffraction (TerrainObstacleLoss mode=diffraction)" +*.radioMedium.obstacleLoss.typename = "TerrainObstacleLoss" +*.radioMedium.obstacleLoss.mode = "diffraction" +*.radioMedium.obstacleLoss.markerStyle = "depth" +*.visualizer.obstacleLossVisualizer.displayIntersections = true +*.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" +*.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s + # -------------------------------------------------------------------------------------------- # Phase-correct two-ray path loss over the real terrain: TwoRayInterference normally assumes # flat ground at z=0 (its antenna "heights" would be ~340 m here, since the Sandy tile's diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc index 71f7a727c4b..4d3eff12151 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc @@ -27,8 +27,11 @@ void TerrainObstacleLoss::initialize(int stage) mode = LOS; else if (!strcmp(modeString, "fresnel")) mode = FRESNEL; + else if (!strcmp(modeString, "diffraction")) + mode = DIFFRACTION; else throw cRuntimeError("Unknown mode '%s'", modeString); + maxDiffractionEdges = par("maxDiffractionEdges"); const char *markerStyleString = par("markerStyle"); if (!strcmp(markerStyleString, "depth")) markerStyle = DEPTH; @@ -51,6 +54,40 @@ double TerrainObstacleLoss::computeKnifeEdgeLoss(double v) return 6.9 + 20 * std::log10(std::sqrt((v - 0.1) * (v - 0.1) + 1) + v - 0.1); } +double TerrainObstacleLoss::computeDeygoutLoss(const std::vector& z, double unitDistance, double waveLength, int i0, int i1, int depth) +{ + // Deygout multi-edge diffraction: find the dominant knife edge (the interior + // sample with the largest diffraction parameter v, measured above the chord + // joining the two sub-path endpoints), charge its J(v), then recurse into the + // paths on either side of it — up to `depth` edges total. This captures the + // graded shadowing behind a row of buildings, where each successive ridge adds + // its own diffraction loss, better than a single dominant edge alone. + if (depth <= 0 || i1 - i0 < 2) + return 0; + double za = z[i0], zb = z[i1]; + int span = i1 - i0; + double bestV = -INFINITY; + int bestIndex = -1; + for (int i = i0 + 1; i < i1; i++) { + if (std::isinf(z[i])) + continue; // no terrain data at this sample: nothing to diffract over + double d1 = (i - i0) * unitDistance; + double d2 = (i1 - i) * unitDistance; + double chord = za + (zb - za) * (double)(i - i0) / span; + double h = z[i] - chord; + double v = h * std::sqrt(2.0 * (d1 + d2) / (waveLength * d1 * d2)); + if (v > bestV) { + bestV = v; + bestIndex = i; + } + } + if (bestIndex < 0 || bestV <= -0.78) + return 0; // no edge intrudes far enough into the Fresnel zone to matter + return computeKnifeEdgeLoss(bestV) + + computeDeygoutLoss(z, unitDistance, waveLength, i0, bestIndex, depth - 1) + + computeDeygoutLoss(z, unitDistance, waveLength, bestIndex, i1, depth - 1); +} + void TerrainObstacleLoss::emitObstaclePenetrated(const Coord& intersection1, const Coord& intersection2, double lossFactor) const { // No physical object is associated with terrain: coordinates are world-frame @@ -128,8 +165,10 @@ const Heightfield *TerrainObstacleLoss::getHeightfield() const std::ostream& TerrainObstacleLoss::printToStream(std::ostream& stream, int level, int evFlags) const { stream << "TerrainObstacleLoss"; - if (level <= PRINT_LEVEL_TRACE) - stream << EV_FIELD(mode, mode == LOS ? "los" : "fresnel") << EV_FIELD(sampleStep); + if (level <= PRINT_LEVEL_TRACE) { + const char *modeName = mode == LOS ? "los" : mode == FRESNEL ? "fresnel" : "diffraction"; + stream << EV_FIELD(mode, modeName) << EV_FIELD(sampleStep); + } return stream; } @@ -145,19 +184,29 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans int numSamples = std::min(4096, std::max(2, (int)(horizontalDistance / step) + 1)); double waveLength = SPEED_OF_LIGHT / frequency.get(); // Walk interior samples only, so an antenna standing on the surface does not - // occlude itself. LOS mode stops at the first obstruction; FRESNEL mode scans - // the whole profile for the worst diffraction parameter v (the deepest - // intrusion into the first Fresnel zone, in zone-normalized units), which the - // single knife-edge curve J(v) then maps to a graded loss. + // occlude itself. LOS mode stops at the first obstruction; the graded modes + // scan the whole profile for the worst diffraction parameter v (the deepest + // intrusion into the first Fresnel zone, in zone-normalized units) — which is + // also the dominant Deygout edge, hence the marker point and the FRESNEL loss. + // DIFFRACTION additionally records the full profile so the Deygout recursion + // can sum the loss over several ridges. bool blocked = false; double worstV = -INFINITY; int worstIndex = -1; double worstGroundZ = NaN; + std::vector profile; + if (mode == DIFFRACTION) { + profile.assign(numSamples, -INFINITY); + profile.front() = transmissionPosition.z; + profile.back() = receptionPosition.z; + } for (int i = 1; i < numSamples - 1; i++) { double t = (double)i / (numSamples - 1); double groundZ = hf->getElevation(transmissionPosition.x + dx * t, transmissionPosition.y + dy * t); if (std::isnan(groundZ)) - continue; // outside the data extent: nothing to block there + continue; // outside the data extent: nothing to block there (profile stays -inf) + if (mode == DIFFRACTION) + profile[i] = groundZ; double rayZ = transmissionPosition.z + dz * t; if (mode == LOS) { if (groundZ > rayZ) { @@ -182,7 +231,11 @@ double TerrainObstacleLoss::computeObstacleLoss(Hz frequency, const Coord& trans if (mode == LOS) lossFactor = blocked ? 0 : 1; else { - double lossDb = std::isfinite(worstV) ? computeKnifeEdgeLoss(worstV) : 0; + double lossDb; + if (mode == DIFFRACTION) + lossDb = computeDeygoutLoss(profile, horizontalDistance / (numSamples - 1), waveLength, 0, numSamples - 1, maxDiffractionEdges); + else + lossDb = std::isfinite(worstV) ? computeKnifeEdgeLoss(worstV) : 0; lossFactor = std::pow(10.0, -lossDb / 10.0); blocked = worstV > 0; // the direct ray itself is geometrically obstructed } diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h index 06ee21109fc..e4004725ad9 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h @@ -35,7 +35,7 @@ namespace physicallayer { class INET_API TerrainObstacleLoss : public TracingObstacleLossBase { protected: - enum Mode { LOS, FRESNEL }; + enum Mode { LOS, FRESNEL, DIFFRACTION }; enum MarkerStyle { DEPTH, RAY }; ModuleRefByPar physicalEnvironment; @@ -43,6 +43,7 @@ class INET_API TerrainObstacleLoss : public TracingObstacleLossBase Mode mode = LOS; MarkerStyle markerStyle = DEPTH; double sampleStep = 0; + int maxDiffractionEdges = 3; bool logLinkEvents = false; // line-of-sight bookkeeping for logging: last known blocked/clear state per node pair @@ -64,6 +65,15 @@ class INET_API TerrainObstacleLoss : public TracingObstacleLossBase */ static double computeKnifeEdgeLoss(double v); + /** + * Total multi-edge diffraction loss in dB over the terrain profile z (heights + * sampled at uniform horizontal spacing unitDistance, with z[i0] and z[i1] the + * path endpoints) using the Deygout construction: charge the dominant knife + * edge in [i0, i1], then recurse into the two sub-paths on either side of it, + * up to depth edges. z[i] == -infinity marks a sample with no terrain data. + */ + static double computeDeygoutLoss(const std::vector& z, double unitDistance, double waveLength, int i0, int i1, int depth); + virtual std::ostream& printToStream(std::ostream& stream, int level, int evFlags = 0) const override; virtual double computeObstacleLoss(Hz frequency, const Coord& transmissionPosition, const Coord& receptionPosition) const override; }; diff --git a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned index c47d67321eb..5de1329f72b 100644 --- a/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned @@ -27,6 +27,11 @@ import inet.physicallayer.wireless.common.base.packetlevel.TracingObstacleLossBa // about 6 dB with the terrain right at the direct ray, deepening smoothly // (and frequency-dependently) as the ray sinks below the terrain. Links // thus degrade progressively instead of dropping at an invisible boundary. +// - "diffraction": multi-edge — the Deygout construction charges the dominant +// knife edge along the profile, then recurses into the sub-paths on either +// side of it (up to maxDiffractionEdges edges), summing the losses. This +// models the added attenuation behind a row of ridges/buildings, where the +// single-edge "fresnel" mode would underestimate the shadow. // // Emits obstaclePenetrated events at the critical terrain point (with no // associated physical object, world coordinates), so tracing obstacle loss @@ -45,9 +50,10 @@ import inet.physicallayer.wireless.common.base.packetlevel.TracingObstacleLossBa module TerrainObstacleLoss extends TracingObstacleLossBase { parameters: - string mode @enum("los", "fresnel") = default("los"); // binary line-of-sight blocking, or graded Fresnel/knife-edge attenuation + string mode @enum("los", "fresnel", "diffraction") = default("los"); // binary line-of-sight blocking, graded single knife-edge (Fresnel), or multi-edge Deygout diffraction string markerStyle @enum("depth", "ray") = default("depth"); // geometry of the emitted obstruction markers (see module description) double sampleStep @unit(m) = default(0m); // spacing of the line-of-sight samples; 0 = the ground heightfield's cell size + int maxDiffractionEdges = default(3); // in "diffraction" mode, the maximum number of Deygout knife edges summed along a path bool logLinkEvents = default(true); // log line-of-sight transitions per node pair (established / lost / re-established / cannot be established) @class(TerrainObstacleLoss); } diff --git a/tests/unit/TerrainObstacleLoss_1.test b/tests/unit/TerrainObstacleLoss_1.test index 518dac50613..ffff3af2adc 100644 --- a/tests/unit/TerrainObstacleLoss_1.test +++ b/tests/unit/TerrainObstacleLoss_1.test @@ -2,10 +2,13 @@ Test TerrainObstacleLoss::computeKnifeEdgeLoss: the ITU-R P.526 single knife-edge diffraction curve J(v). Checks the no-loss region below the v = -0.78 cutoff, continuity at the cutoff, the ~6 dB grazing value at v = 0, -reference points in the shadow region, and monotonicity. +reference points in the shadow region, and monotonicity. Also tests the +Deygout multi-edge construction (computeDeygoutLoss): single-edge equivalence, +clear path, superadditivity of a second ridge, and the edge-count cap. %includes: #include +#include #include "inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h" using namespace inet; @@ -44,6 +47,43 @@ CHECK("shadow v=2", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(2) - 19. CHECK("monotone", monotone); } +// --- Deygout: a single knife edge equals J(v) of that edge --- +{ + double lambda = 0.125, d = 100.0; + std::vector z = {0.0, 30.0, 0.0}; // one peak, endpoints at ground + double v = 30.0 * std::sqrt(2.0 * (d + d) / (lambda * d * d)); // dominant edge at index 1 + double single = TerrainObstacleLoss::computeKnifeEdgeLoss(v); + double deygout = TerrainObstacleLoss::computeDeygoutLoss(z, d, lambda, 0, 2, 3); + CHECK("deygout single edge", std::fabs(deygout - single) < 1e-9); +} + +// --- Deygout: a path clear of every sample has no loss --- +{ + std::vector z = {0.0, -10.0, -5.0, -8.0, 0.0}; + CHECK("deygout clear", TerrainObstacleLoss::computeDeygoutLoss(z, 50.0, 0.125, 0, 4, 3) == 0); +} + +// --- Deygout: a second ridge adds loss beyond the dominant edge alone --- +{ + double lambda = 0.125, d = 50.0; + std::vector one = {0.0, 40.0, 0.0, 0.0, 0.0}; // single dominant ridge + std::vector two = {0.0, 40.0, 0.0, 25.0, 0.0}; // plus a secondary ridge + double lossOne = TerrainObstacleLoss::computeDeygoutLoss(one, d, lambda, 0, 4, 3); + double lossTwo = TerrainObstacleLoss::computeDeygoutLoss(two, d, lambda, 0, 4, 3); + CHECK("deygout second ridge adds loss", lossTwo > lossOne + 0.5); +} + +// --- Deygout: the edge cap limits how many edges are charged --- +{ + double lambda = 0.125, d = 50.0; + std::vector z = {0.0, 40.0, 0.0, 25.0, 0.0}; + double capped = TerrainObstacleLoss::computeDeygoutLoss(z, d, lambda, 0, 4, 1); // dominant edge only + double full = TerrainObstacleLoss::computeDeygoutLoss(z, d, lambda, 0, 4, 3); + CHECK("deygout edge cap limits loss", capped < full); + double v1 = 40.0 * std::sqrt(2.0 * (4 * d) / (lambda * (1 * d) * (3 * d))); // dominant edge at index 1 + CHECK("deygout cap equals dominant edge", std::fabs(capped - TerrainObstacleLoss::computeKnifeEdgeLoss(v1)) < 1e-9); +} + EV << "done\n"; %contains: stdout @@ -54,4 +94,9 @@ grazing v=0: OK shadow v=1: OK shadow v=2: OK monotone: OK +deygout single edge: OK +deygout clear: OK +deygout second ridge adds loss: OK +deygout edge cap limits loss: OK +deygout cap equals dominant edge: OK done From 4f9a33b82730fc051d5af8ffb6bf99a9c95081de Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 16:21:33 +0300 Subject: [PATCH 14/17] visualizer: vsg: render a ground module's heightfield as a shaded mesh MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New createTerrainMeshFromHeightfield() + a lit terrain-mesh pipeline (position + colour + normal, TRIANGLE_LIST, two-sided): builds an elevation-colored, diffuse-shaded triangle surface from a Heightfield in its own simulation coordinates, so the displayed ground is exactly the DSM the radio models sample — no separate display transform that could drift from the physics. Cells with no data drop the triangles that touch them. Wired via a new SceneVsgVisualizerBase.groundModel parameter naming a PointCloudGround module (built by INITSTAGE_LAST, when the scene floor is created); it takes precedence over the sceneModel point cloud. Demo config TerrainMesh renders the vsglidar tile's physics heightfield. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 11 ++ .../vsg/base/SceneVsgVisualizerBase.cc | 17 ++- .../vsg/base/SceneVsgVisualizerBase.ned | 1 + src/inet/visualizer/vsg/util/VsgUtils.cc | 137 ++++++++++++++++++ src/inet/visualizer/vsg/util/VsgUtils.h | 7 + 5 files changed, 172 insertions(+), 1 deletion(-) diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index faef55d216a..7a0a6b1a008 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -255,6 +255,17 @@ description = "graded Fresnel/knife-edge terrain attenuation (TerrainObstacleLos *.visualizer.obstacleLossVisualizer.fadeOutMode = "simulationTime" *.visualizer.obstacleLossVisualizer.fadeOutTime = 1.5s +# -------------------------------------------------------------------------------------------- +# Ground-mesh render: instead of the raw LIDAR point cloud, draw the physics ground module's +# own heightfield as a shaded, elevation-colored surface. What you see is exactly the surface +# the radio models sample (same rasterized DSM, same de-spiking) — the display and the physics +# can no longer disagree. Uses the same PointCloudGround the two-ray/occlusion configs rely on. +# -------------------------------------------------------------------------------------------- +[Config TerrainMesh] +extends = TerrainTwoRay +description = "render the physics heightfield as a shaded ground surface (groundModel)" +*.visualizer.sceneVisualizer.groundModel = "physicalEnvironment.ground" + # -------------------------------------------------------------------------------------------- # Multi-edge diffraction: the Deygout construction sums the knife-edge loss over the dominant # ridge on the path plus the secondary ridges on either side of it (up to maxDiffractionEdges). diff --git a/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.cc b/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.cc index bd9cef72e40..e751ffef203 100644 --- a/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.cc +++ b/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.cc @@ -10,6 +10,7 @@ #include #include "inet/common/ModuleAccess.h" +#include "inet/environment/ground/PointCloudGround.h" #include "inet/visualizer/vsg/scene/NetworkNodeVsgVisualizer.h" #include "inet/visualizer/vsg/util/VsgScene.h" #include "inet/visualizer/vsg/util/VsgUtils.h" @@ -72,8 +73,22 @@ void SceneVsgVisualizerBase::initializeSceneFloor() if (sceneBounds.getMin() == sceneBounds.getMax()) return; auto scene = inet::vsg::TopLevelScene::getSimulationScene(visualizationTargetModule); + const char *groundModel = par("groundModel"); const char *sceneModel = par("sceneModel"); - if (sceneModel != nullptr && *sceneModel != '\0') { + if (groundModel != nullptr && *groundModel != '\0') { + // Render the physics ground's own heightfield as a shaded mesh, so the displayed surface is + // exactly the one the radio models sample. The path (relative to the network) must name a + // PointCloudGround module; its heightfield is already built at this init stage (INITSTAGE_LAST). + auto module = getSimulation()->getSystemModule()->getModuleByPath(groundModel); + auto ground = dynamic_cast(module); + if (ground == nullptr) + throw cRuntimeError("groundModel '%s' does not resolve to a PointCloudGround module", groundModel); + const Heightfield& heightfield = ground->getHeightfield(); + if (!heightfield.isValid()) + throw cRuntimeError("groundModel '%s' has no valid heightfield", groundModel); + scene->addChild(inet::vsg::createTerrainMeshFromHeightfield(heightfield)); + } + else if (sceneModel != nullptr && *sceneModel != '\0') { // Load an external terrain model (a PLY point cloud, e.g. a LIDAR scan) as the ground, in place // of the flat quad. The path is resolved relative to the working directory (the example dir). auto terrain = inet::vsg::createTerrainFromPLY(sceneModel, sceneBounds.getMin(), sceneBounds.getMax()); diff --git a/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.ned b/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.ned index aa3eb59368c..3c16842cd1f 100644 --- a/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.ned +++ b/src/inet/visualizer/vsg/base/SceneVsgVisualizerBase.ned @@ -32,6 +32,7 @@ simple SceneVsgVisualizerBase extends SceneVisualizerBase double sceneImageSize @unit(m) = default(1m); // Scene texture size bool sceneShading = default(true); // Exponential shading for scene color, enabled by default string sceneModel = default(""); // Path to an external terrain model (a PLY point cloud, e.g. a LIDAR scan) to use as the ground instead of the flat quad; recentred and aspect-fit into the scene bounds. Resolved relative to the working directory. Empty = flat floor. + string groundModel = default(""); // Network-relative path of a PointCloudGround module whose heightfield is rendered as a shaded surface mesh, so the displayed ground is exactly the one the physics samples. Takes precedence over sceneModel. Empty = use sceneModel/flat floor. @class(SceneVsgVisualizerBase); } diff --git a/src/inet/visualizer/vsg/util/VsgUtils.cc b/src/inet/visualizer/vsg/util/VsgUtils.cc index 45a40cc3755..ad36e5add0a 100644 --- a/src/inet/visualizer/vsg/util/VsgUtils.cc +++ b/src/inet/visualizer/vsg/util/VsgUtils.cc @@ -15,6 +15,7 @@ #include #include +#include "inet/common/geometry/common/Heightfield.h" #include "inet/common/geometry/common/PlyPointCloudReader.h" namespace inet { @@ -767,6 +768,142 @@ ref_ptr createTerrainFromPLY(const std::string& path, const Coord& sceneMi return stateGroup; } +// --------------------------------------------------------------------------------------------- +// Ground mesh from a Heightfield (renders the SAME surface the physics uses) +// --------------------------------------------------------------------------------------------- + +static const char *terrainMeshVertexShader = R"(#version 450 +layout(push_constant) uniform PushConstants { mat4 projection; mat4 modelView; } pc; +layout(location = 0) in vec3 vsg_Vertex; +layout(location = 1) in vec4 vsg_Color; +layout(location = 2) in vec3 vsg_Normal; +layout(location = 0) out vec4 vColor; +layout(location = 1) out vec3 vWorldNormal; +void main() { + vColor = vsg_Color; + vWorldNormal = vsg_Normal; // heightfield normals are already world/up-referenced (no model rotation) + gl_Position = (pc.projection * pc.modelView) * vec4(vsg_Vertex, 1.0); +} +)"; + +static const char *terrainMeshFragmentShader = R"(#version 450 +layout(location = 0) in vec4 vColor; +layout(location = 1) in vec3 vWorldNormal; +layout(location = 0) out vec4 fragColor; +void main() { + vec3 n = normalize(vWorldNormal); + vec3 lightDir = normalize(vec3(0.35, 0.35, 0.87)); // upper-front, mirrors the scene key light + float diffuse = max(dot(n, lightDir), 0.0); + float shade = 0.40 + 0.60 * diffuse; // ambient floor + diffuse, so buildings read as 3D + fragColor = vec4(vColor.rgb * shade, vColor.a); +} +)"; + +// Cached opaque terrain-mesh pipeline (position + colour + normal, TRIANGLE_LIST, two-sided, depth test + write). +static ref_ptr getTerrainMeshPipeline() +{ + static ref_ptr& pipeline = *(new ref_ptr()); + if (pipeline) return pipeline; + + auto vertexShader = ShaderStage::create(VK_SHADER_STAGE_VERTEX_BIT, "main", terrainMeshVertexShader); + auto fragmentShader = ShaderStage::create(VK_SHADER_STAGE_FRAGMENT_BIT, "main", terrainMeshFragmentShader); + ShaderStages shaderStages{vertexShader, fragmentShader}; + auto shaderCompiler = ShaderCompiler::create(); + if (!shaderCompiler->supported() || !shaderCompiler->compile(shaderStages)) + throw cRuntimeError("groundModel terrain-mesh rendering needs a VSG built with the GLSL compiler (glslang)"); + + auto pipelineLayout = PipelineLayout::create(DescriptorSetLayouts{}, + PushConstantRanges{{VK_SHADER_STAGE_VERTEX_BIT, 0, 128}}); + VertexInputState::Bindings vertexBindings{ + {0, sizeof(::vsg::vec3), VK_VERTEX_INPUT_RATE_VERTEX}, + {1, sizeof(::vsg::vec4), VK_VERTEX_INPUT_RATE_VERTEX}, + {2, sizeof(::vsg::vec3), VK_VERTEX_INPUT_RATE_VERTEX}}; + VertexInputState::Attributes vertexAttributes{ + {0, 0, VK_FORMAT_R32G32B32_SFLOAT, 0}, + {1, 1, VK_FORMAT_R32G32B32A32_SFLOAT, 0}, + {2, 2, VK_FORMAT_R32G32B32_SFLOAT, 0}}; + auto inputAssembly = InputAssemblyState::create(); + inputAssembly->topology = VK_PRIMITIVE_TOPOLOGY_TRIANGLE_LIST; + auto rasterization = RasterizationState::create(); + rasterization->cullMode = VK_CULL_MODE_NONE; // viewable from above or below + auto depthStencil = DepthStencilState::create(); // defaults: depth test + write on (opaque) + auto colorBlend = ColorBlendState::create(); + VkPipelineColorBlendAttachmentState att{}; + att.blendEnable = VK_FALSE; // opaque terrain + att.colorWriteMask = VK_COLOR_COMPONENT_R_BIT | VK_COLOR_COMPONENT_G_BIT | VK_COLOR_COMPONENT_B_BIT | VK_COLOR_COMPONENT_A_BIT; + colorBlend->attachments = ColorBlendState::ColorBlendAttachments{att}; + + GraphicsPipelineStates pipelineStates{ + VertexInputState::create(vertexBindings, vertexAttributes), + inputAssembly, rasterization, MultisampleState::create(), depthStencil, colorBlend}; + pipeline = GraphicsPipeline::create(pipelineLayout, shaderStages, pipelineStates); + return pipeline; +} + +ref_ptr createTerrainMeshFromHeightfield(const Heightfield& heightfield) +{ + // The heightfield is already in simulation coordinates (the ground module baked + // its transform in), so its cells map 1:1 onto the scene the nodes live in — the + // rendered surface is exactly the surface the physics samples. Sample one vertex + // per cell center; a cell with no data (NaN) drops the triangles that touch it. + int nx = heightfield.getNumCellsX(), ny = heightfield.getNumCellsY(); + if (nx < 2 || ny < 2) + return ::vsg::Group::create(); + double cell = heightfield.getCellSize(); + double x0 = heightfield.getMinX() + cell * 0.5, y0 = heightfield.getMinY() + cell * 0.5; + + std::vector zs(nx * ny); + double minZ = INFINITY, maxZ = -INFINITY; + for (int iy = 0; iy < ny; iy++) + for (int ix = 0; ix < nx; ix++) { + double z = heightfield.getElevation(x0 + ix * cell, y0 + iy * cell); + zs[iy * nx + ix] = (float)z; + if (std::isfinite(z)) { minZ = std::min(minZ, z); maxZ = std::max(maxZ, z); } + } + double dz = (std::isfinite(minZ) && maxZ > minZ) ? maxZ - minZ : 1.0; + + auto vertices = vec3Array::create(nx * ny); + auto colors = vec4Array::create(nx * ny); + auto normals = vec3Array::create(nx * ny); + for (int iy = 0; iy < ny; iy++) + for (int ix = 0; ix < nx; ix++) { + int i = iy * nx + ix; + double x = x0 + ix * cell, y = y0 + iy * cell; + float z = zs[i]; + vertices->set(i, ::vsg::vec3((float)x, (float)y, std::isfinite(z) ? z : (float)minZ)); + colors->set(i, elevationColor(std::isfinite(z) ? (z - minZ) / dz : 0.0)); + Coord n = heightfield.getNormal(x, y); + normals->set(i, std::isfinite(n.z) ? ::vsg::vec3((float)n.x, (float)n.y, (float)n.z) : ::vsg::vec3(0.0f, 0.0f, 1.0f)); + } + + // two triangles per cell, skipping any quad with a missing (NaN) corner + std::vector idx; + idx.reserve((size_t)(nx - 1) * (ny - 1) * 6); + auto valid = [&](int ix, int iy) { return std::isfinite(zs[iy * nx + ix]); }; + for (int iy = 0; iy < ny - 1; iy++) + for (int ix = 0; ix < nx - 1; ix++) { + if (!valid(ix, iy) || !valid(ix + 1, iy) || !valid(ix, iy + 1) || !valid(ix + 1, iy + 1)) + continue; + uint32_t a = iy * nx + ix, b = iy * nx + ix + 1, c = (iy + 1) * nx + ix, d = (iy + 1) * nx + ix + 1; + idx.insert(idx.end(), {a, c, b, b, c, d}); + } + if (idx.empty()) + return ::vsg::Group::create(); + auto indices = uintArray::create(idx.size()); + for (size_t i = 0; i < idx.size(); i++) + indices->set(i, idx[i]); + + auto stateGroup = StateGroup::create(); + stateGroup->add(BindGraphicsPipeline::create(getTerrainMeshPipeline())); + auto commands = Commands::create(); + DataList arrays{vertices, colors, normals}; + commands->addChild(BindVertexBuffers::create(0, arrays)); + commands->addChild(BindIndexBuffer::create(indices)); + commands->addChild(DrawIndexed::create(indices->size(), 1, 0, 0, 0)); + stateGroup->addChild(commands); + return stateGroup; +} + // --------------------------------------------------------------------------------------------- // Vertex-array builders // --------------------------------------------------------------------------------------------- diff --git a/src/inet/visualizer/vsg/util/VsgUtils.h b/src/inet/visualizer/vsg/util/VsgUtils.h index 8e705166d27..577f7c88faf 100644 --- a/src/inet/visualizer/vsg/util/VsgUtils.h +++ b/src/inet/visualizer/vsg/util/VsgUtils.h @@ -39,6 +39,8 @@ namespace inet { +class Heightfield; + namespace vsg { using namespace ::vsg; @@ -106,6 +108,11 @@ ref_ptr createSphereWaveShader(double innerRadius, double outerRadius, con // and return it as a coloured POINT_LIST node, recentred and aspect-fit into [sceneMin, sceneMax] and // coloured by elevation when the file has no colours. For a LIDAR/terrain "ground". Empty group on failure. ref_ptr createTerrainFromPLY(const std::string& path, const Coord& sceneMin, const Coord& sceneMax); + +// Renders a Heightfield (a ground module's digital surface model) as a shaded, elevation-colored +// triangle mesh, in the heightfield's own (simulation) coordinates — so the displayed surface is +// exactly the one the physics samples. Cells with no data drop the triangles that touch them. +ref_ptr createTerrainMeshFromHeightfield(const Heightfield& heightfield); ref_ptr createQuad(const Coord& min, const Coord& max, const cFigure::Color& color, double opacity = 1.0); ref_ptr createPolygon(const std::vector& points, const cFigure::Color& color, double opacity = 1.0, const Coord& translation = Coord::ZERO); ref_ptr createArrowhead(const Coord& start, const Coord& end, const cFigure::Color& color, double width = 10.0, double height = 20.0, double opacity = 1.0); From 32e3b12885dddedb7113027d31aae70f8cda8cce Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 17:07:52 +0300 Subject: [PATCH 15/17] common: Heightfield: de-spike against the second-highest neighbor MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The de-spike previously clamped only cells taller than EVERY neighbor by the threshold, which left two-cell spikes (each shielded by the other) and any spike shallower than the threshold above the local grade. Compare each cell to its SECOND-highest neighbor instead, over a few erosion passes: this catches isolated and two-cell spikes and shallower ones, while genuine structures are preserved — any cell on a real surface or edge has at least two neighbors at its height, so its second-highest neighbor is high and it is left alone. Verified against the vsglidar tile: the real buildings (e.g. 154 m / 146 m / the 289 m tower) are byte-identical across thresholds 0/20/50, while stray returns are pulled to the surrounding grade. vsglidar demo drops despikeThreshold to 20 m to also clear the shallower plain spikes; the swarm physics regressions are unchanged (the noise was never on a link path). Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 1 + .../common/geometry/common/Heightfield.cc | 58 ++++++++++++------- 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index 7a0a6b1a008..2dd0e77affd 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -149,6 +149,7 @@ description = "Two-ray ground reflection over the real LIDAR surface (PointCloud *.physicalEnvironment.ground.typename = "PointCloudGround" *.physicalEnvironment.ground.file = "terrain.ply" *.physicalEnvironment.ground.transformMode = "fit" +*.physicalEnvironment.ground.despikeThreshold = 20m # pull stray LIDAR returns down to their surroundings (birds/noise sit ~40-90m over the ~45m grade); real multi-cell buildings keep their height *.physicalEnvironment.ground.fitMinX = 0m *.physicalEnvironment.ground.fitMinY = 0m *.physicalEnvironment.ground.fitMinZ = 0m diff --git a/src/inet/common/geometry/common/Heightfield.cc b/src/inet/common/geometry/common/Heightfield.cc index 80c1da5b5fd..ef9c13f0c66 100644 --- a/src/inet/common/geometry/common/Heightfield.cc +++ b/src/inet/common/geometry/common/Heightfield.cc @@ -7,6 +7,7 @@ #include "inet/common/geometry/common/Heightfield.h" +#include #include #include @@ -55,33 +56,46 @@ void Heightfield::buildFromPoints(const std::vector& xs, const std::vect cell = z; } - // clamp isolated single-cell spikes (stray LIDAR returns) that exceed ALL - // non-NaN neighbors by more than the threshold; real structures survive - // because their edge cells have same-height neighbors + // clamp spikes (stray LIDAR returns: birds, atmospheric noise) that stand more + // than the threshold above their SECOND-highest non-NaN neighbor, pulling each + // down to that level. Testing the second-highest (rather than the highest, as + // before) also catches two-cell spikes — where each cell's single tall neighbor + // is the other half of the spike — while genuine structures survive: any cell on + // a real edge or surface has at least two neighbors sharing its height, so its + // second-highest neighbor is high and it is left alone. A few erosion passes peel + // wider spike clusters from the outside in. if (despikeThreshold > 0) { - std::vector> clamps; - for (int iy = 0; iy < numCellsY; iy++) { - for (int ix = 0; ix < numCellsX; ix++) { - float v = getCell(ix, iy); - if (std::isnan(v)) - continue; - float maxNeighbor = -std::numeric_limits::infinity(); - int n = 0; - for (int dy = -1; dy <= 1; dy++) { - for (int dx = -1; dx <= 1; dx++) { - if (dx == 0 && dy == 0) continue; - int jx = ix + dx, jy = iy + dy; - if (jx < 0 || jx >= numCellsX || jy < 0 || jy >= numCellsY) continue; - float w = getCell(jx, jy); - if (!std::isnan(w)) { maxNeighbor = std::max(maxNeighbor, w); n++; } + for (int pass = 0; pass < 3; pass++) { + std::vector> clamps; + for (int iy = 0; iy < numCellsY; iy++) { + for (int ix = 0; ix < numCellsX; ix++) { + float v = getCell(ix, iy); + if (std::isnan(v)) + continue; + float highest = -std::numeric_limits::infinity(); + float secondHighest = -std::numeric_limits::infinity(); + int n = 0; + for (int dy = -1; dy <= 1; dy++) { + for (int dx = -1; dx <= 1; dx++) { + if (dx == 0 && dy == 0) continue; + int jx = ix + dx, jy = iy + dy; + if (jx < 0 || jx >= numCellsX || jy < 0 || jy >= numCellsY) continue; + float w = getCell(jx, jy); + if (std::isnan(w)) continue; + n++; + if (w > highest) { secondHighest = highest; highest = w; } + else if (w > secondHighest) secondHighest = w; + } } + if (n >= 3 && v > secondHighest + despikeThreshold) + clamps.emplace_back((size_t)(iy * numCellsX + ix), secondHighest); } - if (n >= 3 && v > maxNeighbor + despikeThreshold) - clamps.emplace_back((size_t)(iy * numCellsX + ix), maxNeighbor); } + if (clamps.empty()) + break; + for (auto& c : clamps) + cells[c.first] = c.second; } - for (auto& c : clamps) - cells[c.first] = c.second; } // bounded hole filling: NaN cells take the average of their non-NaN 8-neighbors From 62588ec550a9d926721e4e7259cb67bc957af10b Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 17:15:28 +0300 Subject: [PATCH 16/17] examples: vsglidar: add TerrainMeshOcclusion (mesh render + live occlusion) Combines the shaded ground-mesh render with binary line-of-sight blocking and the moving swarm, so occlusion is watchable: as a drone passes behind a building its data-link arrow to the gcs breaks and fades, then returns when it comes back into view. TerrainMesh alone extends TerrainTwoRay (no obstacle loss), so there the arrows draw straight through the buildings. Co-Authored-By: Claude Fable 5 --- examples/visualizer/vsglidar/omnetpp.ini | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/examples/visualizer/vsglidar/omnetpp.ini b/examples/visualizer/vsglidar/omnetpp.ini index 2dd0e77affd..3e467963c99 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -267,6 +267,18 @@ extends = TerrainTwoRay description = "render the physics heightfield as a shaded ground surface (groundModel)" *.visualizer.sceneVisualizer.groundModel = "physicalEnvironment.ground" +# -------------------------------------------------------------------------------------------- +# The occlusion story, made watchable: the shaded terrain mesh (so the buildings are solid +# forms) PLUS binary line-of-sight blocking, with the swarm flying among the buildings. As a +# drone passes behind a building its direct ray to the gcs is cut, its pings stop arriving, and +# the data-link arrow pair fades within ~1.5s; when it comes back into view the arrows return. +# (TerrainMesh alone has no obstacle loss, so there the arrows draw straight through buildings.) +# -------------------------------------------------------------------------------------------- +[Config TerrainMeshOcclusion] +extends = TerrainOcclusion +description = "swarm over the shaded mesh: data-link arrows break as drones move behind buildings" +*.visualizer.sceneVisualizer.groundModel = "physicalEnvironment.ground" + # -------------------------------------------------------------------------------------------- # Multi-edge diffraction: the Deygout construction sums the knife-edge loss over the dominant # ridge on the path plus the secondary ridges on either side of it (up to maxDiffractionEdges). From ca589d400d2c41db74c11ee8423a50d339cb2003 Mon Sep 17 00:00:00 2001 From: tabgab Date: Mon, 6 Jul 2026 17:37:55 +0300 Subject: [PATCH 17/17] tests: strengthen Heightfield/TerrainObstacleLoss unit tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Heightfield: FlatGround-equivalence check over flat terrain (constant elevation + (0,0,1) normal at many points, grid-aligned or not) — the flat-PLY == FlatGround cross-check at the layer PointCloudGround delegates to; plus a two-cell spike case for the second-highest-neighbor de-spike. - TerrainObstacleLoss: a hand-derived two-edge Deygout worked example (profile [0,30,0,22,0], lambda 0.06 m) asserting the total equals J(v1)+J(v') to 1e-6 (~76 dB), with the derivation documented, and the one-edge cap. Co-Authored-By: Claude Fable 5 --- tests/unit/Heightfield_1.test | 43 +++++++++++++++++++++++++-- tests/unit/TerrainObstacleLoss_1.test | 30 ++++++++++++++++++- 2 files changed, 70 insertions(+), 3 deletions(-) diff --git a/tests/unit/Heightfield_1.test b/tests/unit/Heightfield_1.test index 85b53e98021..4d3b570e999 100644 --- a/tests/unit/Heightfield_1.test +++ b/tests/unit/Heightfield_1.test @@ -1,7 +1,8 @@ %description: -Test Heightfield: bilinear elevation on a flat plane, normals on flat and +Test Heightfield: bilinear elevation on a flat plane, FlatGround equivalence +over flat terrain (the flat-PLY == FlatGround cross-check), normals on flat and ramp surfaces, elevation profile across a step (building wall), out-of-extent -NaN behavior, and the maxCells guard. +NaN behavior, the maxCells guard, and one- and two-cell spike removal. %includes: #include @@ -42,6 +43,27 @@ static double zStep(double x, double) { return x < 10 ? 0 : 30; } CHECK("out of extent NaN", std::isnan(hf.getElevation(-100, -100))); } +// --- flat terrain is FlatGround-equivalent. PointCloudGround.computeGroundProjection returns +// (x, y, heightfield.getElevation(x,y)) and computeGroundNormal returns heightfield.getNormal(), +// while FlatGround returns (x, y, elevation) and (0,0,1). So over a flat tile the heightfield +// must reproduce FlatGround exactly at every point, grid-aligned or not — that is the +// flat-PLY == FlatGround cross-check, at the layer PointCloudGround actually delegates to. --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 40, 0, 40, 1, zFlat); // flat plane at z=5, i.e. FlatGround(elevation=5) + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 2); + bool elevationEquiv = true, normalEquiv = true; + double probes[][2] = {{5, 5}, {12.3, 7.8}, {20, 20}, {33.33, 4.5}, {2.5, 38}, {38.9, 38.9}}; + for (auto& p : probes) { + if (std::fabs(hf.getElevation(p[0], p[1]) - 5) > 1e-6) elevationEquiv = false; // FlatGround z = 5 + Coord n = hf.getNormal(p[0], p[1]); + if (n.z < 1 - 1e-6 || std::fabs(n.x) > 1e-6 || std::fabs(n.y) > 1e-6) normalEquiv = false; // FlatGround (0,0,1) + } + CHECK("flatground elevation equivalence", elevationEquiv); + CHECK("flatground normal equivalence", normalEquiv); +} + // --- ramp z=x: elevation increases along x, normal tilts against the slope --- { std::vector xs, ys, zs; @@ -91,6 +113,19 @@ static double zStep(double x, double) { return x < 10 ? 0 : 30; } CHECK("building survives despike", std::fabs(hf.getElevation(17, 10) - 30) < 1e-6); } +// --- despike catches a TWO-cell spike (each cell shields the other from the "taller than every +// neighbor" test, but not from the second-highest-neighbor test), still sparing the building --- +{ + std::vector xs, ys, zs; + makeGrid(xs, ys, zs, 0, 20, 0, 20, 0.5, zStep); // building at x>=10, z=30 + xs.push_back(4); ys.push_back(4); zs.push_back(400); // two adjacent stray returns (one grid cell apart) + xs.push_back(4.5); ys.push_back(4); zs.push_back(400); + Heightfield hf; + hf.buildFromPoints(xs, ys, zs, 0.5, 67108864, 8, 50); + CHECK("two-cell spike clamped", hf.getElevation(4, 4) < 10 && hf.getElevation(4.5, 4) < 10); + CHECK("building survives two-cell despike", std::fabs(hf.getElevation(17, 10) - 30) < 1e-6); +} + // --- degenerate input rejected --- { std::vector xs(3, 1), ys(3, 1), zs(3, 1); // all points at the same x/y @@ -109,6 +144,8 @@ EV << ".\n"; flat elevation: OK flat normal: OK out of extent NaN: OK +flatground elevation equivalence: OK +flatground normal equivalence: OK ramp monotonic: OK ramp normal tilt: OK step low side: OK @@ -117,5 +154,7 @@ profile spans the wall: OK maxCells guard: OK spike clamped: OK building survives despike: OK +two-cell spike clamped: OK +building survives two-cell despike: OK degenerate rejected: OK . diff --git a/tests/unit/TerrainObstacleLoss_1.test b/tests/unit/TerrainObstacleLoss_1.test index ffff3af2adc..5de4ddda96c 100644 --- a/tests/unit/TerrainObstacleLoss_1.test +++ b/tests/unit/TerrainObstacleLoss_1.test @@ -4,7 +4,8 @@ knife-edge diffraction curve J(v). Checks the no-loss region below the v = -0.78 cutoff, continuity at the cutoff, the ~6 dB grazing value at v = 0, reference points in the shadow region, and monotonicity. Also tests the Deygout multi-edge construction (computeDeygoutLoss): single-edge equivalence, -clear path, superadditivity of a second ridge, and the edge-count cap. +clear path, superadditivity of a second ridge, the edge-count cap, and a +hand-derived two-edge worked example (total == J(v1) + J(v') to 1e-6). %includes: #include @@ -84,6 +85,30 @@ CHECK("shadow v=2", std::fabs(TerrainObstacleLoss::computeKnifeEdgeLoss(2) - 19. CHECK("deygout cap equals dominant edge", std::fabs(capped - TerrainObstacleLoss::computeKnifeEdgeLoss(v1)) < 1e-9); } +// --- Deygout worked example (two edges, hand-derived). Profile z = [0, 30, 0, 22, 0] at 50 m +// spacing, lambda = 0.06 m. Deygout charges the dominant edge, then recurses on the sub-path +// beyond it and charges that dominant edge, and so on: +// * Dominant edge over the whole path [0..4] is index 1 (h = 30 m above the flat endpoint +// chord; d1 = 50 m, d2 = 150 m) -> v1 = 30 * sqrt(2*(d1+d2)/(lambda*d1*d2)) ~= 28.28. +// * Recurse on [1..4]. The endpoint chord runs from z=30 (index 1) to z=0 (index 4); at +// index 3 it is 30 * (1 - 2/3) = 10 m, so that edge stands h' = 22 - 10 = 12 m above it, +// with d1' = 100 m, d2' = 50 m -> v' = 12 * sqrt(2*(d1'+d2')/(lambda*d1'*d2')) = 12.0. +// (Index 2 sits below its chord, so it is not an edge; [1..3] and [3..4] add nothing.) +// Total loss must be exactly J(v1) + J(v'); depth-capped to one edge it must be J(v1) alone. --- +{ + double lambda = 0.06, d = 50.0; + std::vector prof = {0.0, 30.0, 0.0, 22.0, 0.0}; + double v1 = 30.0 * std::sqrt(2.0 * (1 * d + 3 * d) / (lambda * (1 * d) * (3 * d))); // ~28.28 + double chordAt3 = 30.0 + (0.0 - 30.0) * (double)(3 - 1) / (4 - 1); // 10 m + double vPrime = (22.0 - chordAt3) * std::sqrt(2.0 * (2 * d + 1 * d) / (lambda * (2 * d) * (1 * d))); // 12.0 + double expected = TerrainObstacleLoss::computeKnifeEdgeLoss(v1) + TerrainObstacleLoss::computeKnifeEdgeLoss(vPrime); + double total = TerrainObstacleLoss::computeDeygoutLoss(prof, d, lambda, 0, 4, 3); + CHECK("deygout worked example total", std::fabs(total - expected) < 1e-6); + CHECK("deygout worked example magnitude ~76 dB", total > 74.0 && total < 79.0); + double capped = TerrainObstacleLoss::computeDeygoutLoss(prof, d, lambda, 0, 4, 1); + CHECK("deygout worked example one edge", std::fabs(capped - TerrainObstacleLoss::computeKnifeEdgeLoss(v1)) < 1e-6); +} + EV << "done\n"; %contains: stdout @@ -99,4 +124,7 @@ deygout clear: OK deygout second ridge adds loss: OK deygout edge cap limits loss: OK deygout cap equals dominant edge: OK +deygout worked example total: OK +deygout worked example magnitude ~76 dB: OK +deygout worked example one edge: OK done