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 = "" 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..3e467963c99 100644 --- a/examples/visualizer/vsglidar/omnetpp.ini +++ b/examples/visualizer/vsglidar/omnetpp.ini @@ -124,3 +124,185 @@ 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.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 +*.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: 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 = "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. +# 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" +# "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 + +# -------------------------------------------------------------------------------------------- +# 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" +# "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 + +# -------------------------------------------------------------------------------------------- +# 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" + +# -------------------------------------------------------------------------------------------- +# 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). +# 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 +# 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/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/common/geometry/common/Heightfield.cc b/src/inet/common/geometry/common/Heightfield.cc new file mode 100644 index 00000000000..ef9c13f0c66 --- /dev/null +++ b/src/inet/common/geometry/common/Heightfield.cc @@ -0,0 +1,218 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/common/geometry/common/Heightfield.h" + +#include +#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 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) { + 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 (clamps.empty()) + break; + 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/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); +} 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 new file mode 100644 index 00000000000..4d3eff12151 --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.cc @@ -0,0 +1,285 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +#include "inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h" + +#include +#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); + const char *modeString = par("mode"); + if (!strcmp(modeString, "los")) + 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; + else if (!strcmp(markerStyleString, "ray")) + markerStyle = RAY; + else + throw cRuntimeError("Unknown markerStyle '%s'", markerStyleString); + 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); +} + +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 + // (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); +} + +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) { + const char *modeName = mode == LOS ? "los" : mode == FRESNEL ? "fresnel" : "diffraction"; + stream << EV_FIELD(mode, modeName) << 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 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)(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; 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 (profile stays -inf) + if (mode == DIFFRACTION) + profile[i] = groundZ; + double rayZ = transmissionPosition.z + dz * t; + if (mode == LOS) { + if (groundZ > rayZ) { + blocked = true; // terrain or building blocks the direct ray + worstIndex = i; + 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; + worstIndex = i; + worstGroundZ = groundZ; + } + } + } + double lossFactor; + if (mode == LOS) + lossFactor = blocked ? 0 : 1; + else { + 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 + } + 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); + return lossFactor; +} + +} // 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..e4004725ad9 --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.h @@ -0,0 +1,85 @@ +// +// 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/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/base/packetlevel/TracingObstacleLossBase.h" + +namespace inet { + +namespace physicallayer { + +/** + * 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 TracingObstacleLossBase +{ + protected: + enum Mode { LOS, FRESNEL, DIFFRACTION }; + 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; + int maxDiffractionEdges = 3; + 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; + void emitObstaclePenetrated(const Coord& intersection1, const Coord& intersection2, 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); + + /** + * 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; +}; + +} // 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..5de1329f72b --- /dev/null +++ b/src/inet/physicallayer/wireless/common/obstacleloss/TerrainObstacleLoss.ned @@ -0,0 +1,59 @@ +// +// Copyright (C) 2026 OpenSim Ltd. +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// + + +package inet.physicallayer.wireless.common.obstacleloss; + +import inet.physicallayer.wireless.common.base.packetlevel.TracingObstacleLossBase; + +// +// 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. +// +// 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. +// - "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 +// 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. +// +module TerrainObstacleLoss extends TracingObstacleLossBase +{ + parameters: + 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/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) } 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/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/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" 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/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/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. 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/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/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(); diff --git a/src/inet/visualizer/vsg/util/VsgUtils.cc b/src/inet/visualizer/vsg/util/VsgUtils.cc index 12723570835..ad36e5add0a 100644 --- a/src/inet/visualizer/vsg/util/VsgUtils.cc +++ b/src/inet/visualizer/vsg/util/VsgUtils.cc @@ -15,6 +15,9 @@ #include #include +#include "inet/common/geometry/common/Heightfield.h" +#include "inet/common/geometry/common/PlyPointCloudReader.h" + namespace inet { namespace vsg { @@ -708,24 +711,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 +726,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 +753,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) --- @@ -864,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); 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? diff --git a/tests/unit/Heightfield_1.test b/tests/unit/Heightfield_1.test new file mode 100644 index 00000000000..4d3b570e999 --- /dev/null +++ b/tests/unit/Heightfield_1.test @@ -0,0 +1,160 @@ +%description: +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, the maxCells guard, and one- and two-cell spike removal. + +%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))); +} + +// --- 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; + 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); +} + +// --- 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 + 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 +flatground elevation equivalence: OK +flatground normal equivalence: 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 +two-cell spike clamped: OK +building survives two-cell 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 +. diff --git a/tests/unit/TerrainObstacleLoss_1.test b/tests/unit/TerrainObstacleLoss_1.test new file mode 100644 index 00000000000..5de4ddda96c --- /dev/null +++ b/tests/unit/TerrainObstacleLoss_1.test @@ -0,0 +1,130 @@ +%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. Also tests the +Deygout multi-edge construction (computeDeygoutLoss): single-edge equivalence, +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 +#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); +} + +// --- 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); +} + +// --- 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 +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 +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 +deygout worked example total: OK +deygout worked example magnitude ~76 dB: OK +deygout worked example one edge: OK +done