diff --git a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/CutProfile.java b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/CutProfile.java index 3dcca28dd..ed5be1e3a 100644 --- a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/CutProfile.java +++ b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/CutProfile.java @@ -277,7 +277,8 @@ public List computePts2DGround() { } // keep track of the obstacle under our current position. If -1 there is only ground below int overObstacleIndex = getCutPoints().get(0).getBuildingId(); - for (CutPoint cut : getCutPoints()) { + for (int i=0; i < pts.size(); i++) { + CutPoint cut = pts.get(i); if (cut.getType() != GROUND_EFFECT) { Coordinate coordinate; if (BUILDING.equals(cut.getType()) || WALL.equals(cut.getType())) { diff --git a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/ProfileBuilder.java b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/ProfileBuilder.java index f1ee13860..118788373 100644 --- a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/ProfileBuilder.java +++ b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/ProfileBuilder.java @@ -23,6 +23,7 @@ import org.locationtech.jts.index.ItemVisitor; import org.locationtech.jts.index.strtree.STRtree; import org.locationtech.jts.math.Vector2D; +import org.locationtech.jts.math.Vector3D; import org.locationtech.jts.operation.distance.DistanceOp; import org.locationtech.jts.triangulate.quadedge.Vertex; import org.noise_planet.noisemodelling.pathfinder.delaunay.LayerDelaunay; @@ -1311,35 +1312,21 @@ public void addTopoCutPts(Coordinate p1, Coordinate p2, CutProfile profile, bool } else { LOGGER.warn(String.format(Locale.ROOT, "Propagation out of the DEM area from %s to %s", p1.toString(), p2.toString())); + return; } profile.hasTopographyIntersection = !freeField; - // Remove unnecessary points - ArrayList retainedCoordinates = new ArrayList<>(coordinates.size()); - for(int i =0; i < coordinates.size(); i++) { - // Always add first and last points - Coordinate previous; - Coordinate current = coordinates.get(i); - Coordinate next; - if(retainedCoordinates.isEmpty()) { - previous = new Coordinate(p1.x, p1.y, getZGround(p1)); - } else { - previous = retainedCoordinates.get(retainedCoordinates.size() - 1); - } - if(i == coordinates.size() - 1) { - next = new Coordinate(p2.x, p2.y, getZGround(p2)); - } else { - next = coordinates.get(i + 1); - } + // avoid resizing of array by reserving memory + profile.reservePoints(coordinates.size()); + for(int idPoint = 1; idPoint < coordinates.size() - 1; idPoint++) { + final Coordinate previous = coordinates.get(idPoint - 1); + final Coordinate current = coordinates.get(idPoint); + final Coordinate next = coordinates.get(idPoint+1); // Do not add topographic points which are simply the linear interpolation between two points + // triangulation add a lot of interpolated lines from line segment DEM if(CGAlgorithms3D.distancePointSegment(current, previous, next) >= DELTA) { - retainedCoordinates.add(coordinates.get(i)); + profile.addTopoCutPt(current, idPoint); } } - // Feed profile - profile.reservePoints(retainedCoordinates.size()); - for(int i =0; i < retainedCoordinates.size(); i++) { - profile.addTopoCutPt(retainedCoordinates.get(i), i); - } } /** @@ -1390,7 +1377,7 @@ boolean findClosestTriangleIntersection(LineSegment segment, final Coordinate in } /** - * Fetch all intersections with TIN + * Fetch all intersections with TIN. For simplification only plane change are pushed. * @param p1 first point * @param p2 second point * @param stopAtObstacleOverSourceReceiver Stop fetching intersections if the segment p1-p2 is intersecting with TIN @@ -1421,6 +1408,8 @@ public boolean fetchTopographicProfile(List outputPoints,Coordinate // Add p1 coordinate Coordinate[] vertices = getTriangleVertices(curTriP1); outputPoints.add(new Coordinate(p1.x, p1.y, Vertex.interpolateZ(p1, vertices[0], vertices[1], vertices[2]))); + Vector3D previousTriangleNormal = null; + boolean freeField = true; while (navigationTri != -1) { navigationHistory.add(navigationTri); Coordinate intersectionPt = new Coordinate(); @@ -1433,11 +1422,18 @@ public boolean fetchTopographicProfile(List outputPoints,Coordinate // Found next triangle (if propaTri >= 0) // extract X,Y,Z values of intersection with triangle segment if(!Double.isNaN(intersectionPt.z)) { - outputPoints.add(intersectionPt); - if(stopAtObstacleOverSourceReceiver) { - Coordinate closestPointOnPropagationLine = propaLine.closestPoint(intersectionPt); - double interpolatedZ = Vertex.interpolateZ(closestPointOnPropagationLine, propaLine.p0, propaLine.p1); - if(interpolatedZ < intersectionPt.z) { + Coordinate[] trianglePoints = getTriangle(propaTri); + final Vector3D triangleNormal = JTSUtility.getTriangleNormal(trianglePoints[0], trianglePoints[1], trianglePoints[2]); + // We do not push coplanar intersection points + if(previousTriangleNormal == null || Math.abs(computeNormalsAngle(triangleNormal, previousTriangleNormal)) > epsilon) { + outputPoints.add(intersectionPt); + previousTriangleNormal = triangleNormal; + } + Coordinate closestPointOnPropagationLine = propaLine.closestPoint(intersectionPt); + double interpolatedZ = Vertex.interpolateZ(closestPointOnPropagationLine, propaLine.p0, propaLine.p1); + if(interpolatedZ < intersectionPt.z) { + freeField = false; + if(stopAtObstacleOverSourceReceiver) { return false; } } @@ -1445,7 +1441,16 @@ public boolean fetchTopographicProfile(List outputPoints,Coordinate } navigationTri = propaTri; } - return true; + return freeField; + } + + /** + * @param normal1 Normalized vector 1 + * @param normal2 Normalized vector 2 + * @return The angle between the two normals + */ + private double computeNormalsAngle(Vector3D normal1, Vector3D normal2) { + return Math.acos(normal1.dot(normal2)); } /** diff --git a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/utils/geometry/JTSUtility.java b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/utils/geometry/JTSUtility.java index 2e2f74e36..9c50ace8a 100644 --- a/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/utils/geometry/JTSUtility.java +++ b/noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/utils/geometry/JTSUtility.java @@ -17,9 +17,11 @@ import org.locationtech.jts.geom.LineSegment; import org.locationtech.jts.geom.LineString; import org.locationtech.jts.math.Vector2D; +import org.locationtech.jts.math.Vector3D; import java.util.ArrayList; import java.util.List; +import java.util.Vector; import java.util.concurrent.atomic.AtomicReference; /** @@ -186,6 +188,14 @@ public static boolean dotInTri(Coordinate p, Coordinate a, Coordinate b, return dotInTri(p, a, b, c, null); } + public static Vector3D getTriangleNormal(Coordinate p1, Coordinate p2, Coordinate p3) { + Vector3D a = new Vector3D(p1, p2); + Vector3D b = new Vector3D(p1, p3); + return new Vector3D(a.getY() * b.getZ() - a.getZ() * b.getY(), + a.getZ() * b.getX() - a.getX() * b.getZ(), + a.getX()*b.getY() - a.getY() * b.getX()).normalize(); + } + /** * Fast dot in triangle test * http://www.blackpawn.com/texts/pointinpoly/default.html @@ -235,6 +245,17 @@ public static boolean dotInTri(Coordinate p, Coordinate a, Coordinate b, * @return X Z projected points */ public static List getNewCoordinateSystem(List listPoints) { + return getNewCoordinateSystem(listPoints, 0); + } + /** + * ChangeCoordinateSystem, use original coordinate in 3D to change into a new markland in 2D + * with new x' computed by algorithm and y' is original height of point. + * @param listPoints X Y Z points, all should be on the same plane as first and last points. + * @param tolerance Simplify the point list by not adding points where the distance from the line segments + * formed from the previous and the next point is inferior to this tolerance (remove intermediate collinear points) + * @return X Z projected points + */ + public static List getNewCoordinateSystem(List listPoints, double tolerance) { if(listPoints.isEmpty()) { return new ArrayList<>(); } @@ -245,6 +266,20 @@ public static List getNewCoordinateSystem(List listPoint // Get 2D distance newCoord.add(new Coordinate(newCoord.get(idPoint - 1).x + pt.distance(listPoints.get(idPoint - 1)), pt.z)); } + if(tolerance > 0) { + // remove collinear points using tolerance + for (int idPoint = 1; idPoint < newCoord.size() - 1;) { + final Coordinate previous = newCoord.get(idPoint - 1); + final Coordinate current = newCoord.get(idPoint); + final Coordinate next = newCoord.get(idPoint+1); + final LineSegment lineSegment = new LineSegment(previous, next); + if(lineSegment.distance(current) < tolerance) { + newCoord.remove(idPoint); + } else { + idPoint++; + } + } + } return newCoord; } diff --git a/noisemodelling-pathfinder/src/test/java/org/noise_planet/noisemodelling/pathfinder/PathFinderTest.java b/noisemodelling-pathfinder/src/test/java/org/noise_planet/noisemodelling/pathfinder/PathFinderTest.java index 016ebf3dd..db5159be0 100644 --- a/noisemodelling-pathfinder/src/test/java/org/noise_planet/noisemodelling/pathfinder/PathFinderTest.java +++ b/noisemodelling-pathfinder/src/test/java/org/noise_planet/noisemodelling/pathfinder/PathFinderTest.java @@ -16,6 +16,10 @@ import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.LineSegment; +import org.locationtech.jts.math.Plane3D; +import org.locationtech.jts.math.Vector3D; +import org.locationtech.jts.operation.distance3d.PlanarPolygon3D; import org.noise_planet.noisemodelling.pathfinder.cnossos.CnossosPath; //import org.noise_planet.noisemodelling.pathfinder.path.CnossosPathParameters; import org.noise_planet.noisemodelling.pathfinder.path.Scene; @@ -235,6 +239,46 @@ public static void addTopographicTC5Model(ProfileBuilder profileBuilder) { .addTopographicLine(205, -5, 10, 225, -20, 0); } + + public static void addTopographicTC23Model(ProfileBuilder profileBuilder) { + // Create parallel lines for the slope edge because unit test table values are rounded and + // the rounding make the lines non-parallel + + // base line left + Vector3D v1 = new Vector3D(new Coordinate(46.27, 36.28, 0), + new Coordinate(59.6, -9.87, 0)); + + // original top segment (left side) + Vector3D v2 = new Vector3D(new Coordinate(54.68, 37.59, 5), + new Coordinate(67.35, -6.83, 5)); + + + // base line right + Vector3D v3 = new Vector3D(new Coordinate(63.71, 41.16, 0), + new Coordinate(76.84, -5.28, 0)); + + Vector3D parallelPoint1 = new Vector3D(new Coordinate(54.68, 37.59, 5)).add(v1.normalize().divide(1/v2.length())); + Vector3D parallelPoint2 = new Vector3D(new Coordinate(55.93, 37.93, 5)).add(v3.normalize().divide(1/v2.length())); + + + profileBuilder.addTopographicLine(30, -14, 0, 122, -14, 0)// 1 + .addTopographicLine(122, -14, 0, 122, 45, 0)// 2 + .addTopographicLine(122, 45, 0, 30, 45, 0)// 3 + .addTopographicLine(30, 45, 0, 30, -14, 0)// 4 + .addTopographicLine(59.6, -9.87, 0, 76.84, -5.28, 0)// 5 + .addTopographicLine(76.84, -5.28, 0, 63.71, 41.16, 0)// 6 + .addTopographicLine(63.71, 41.16, 0, 46.27, 36.28, 0)// 7 + .addTopographicLine(46.27, 36.28, 0, 59.6, -9.87, 0)// 8 + .addTopographicLine(46.27, 36.28, 0, 54.68, 37.59, 5)// 9 + .addTopographicLine(54.68, 37.59, 5, parallelPoint2.getX(), parallelPoint2.getY(), 5)// 10 + .addTopographicLine(55.93, 37.93, 5, 63.71, 41.16, 0)// 11 + .addTopographicLine(59.6, -9.87, 0, parallelPoint1.getX(), parallelPoint1.getY(), 5)// 12 + .addTopographicLine(parallelPoint1.getX(), parallelPoint1.getY(), 5, parallelPoint2.getX(), parallelPoint2.getY(), 5)// 13 + .addTopographicLine(parallelPoint2.getX(), parallelPoint2.getY(), 5, 76.84, -5.28, 0)// 14 + .addTopographicLine(54.68, 37.59, 5, parallelPoint1.getX(), parallelPoint1.getY(), 5)// 15 + .addTopographicLine(55.93, 37.93, 5, parallelPoint2.getX(), parallelPoint2.getY(), 5); // 16 + } + /** * Test TC05 -- Ground with spatially varying heights and acoustic properties */ @@ -1329,19 +1373,20 @@ public void TC18() { //Run computation computeRays.run(propDataOut); - CutProfile cutProfile = computeRays.getData().profileBuilder.getProfile(rayData.sourceGeometries.get(0).getCoordinate(), rayData.receivers.get(0), computeRays.getData().gS, false); - List result = cutProfile.computePts2DGround(); - - // Expected Values - /* Table 193 */ + /* Table 193 Z Profile SR */ List expectedZ_profile = new ArrayList<>(); expectedZ_profile.add(new Coordinate(0.0, 0.0)); expectedZ_profile.add(new Coordinate(112.41, 0.0)); expectedZ_profile.add(new Coordinate(178.84, 10)); expectedZ_profile.add(new Coordinate(194.16, 10)); + CutProfile cutProfile = propDataOut.getPropagationPaths().get(0).getCutProfile(); + List result = cutProfile.computePts2DGround(); + assertZProfil(expectedZ_profile, result); + + /* Table 194 */ double [][] segmentsMeanPlanes0 = new double[][]{ // a b zs zr dp Gp Gp' @@ -1356,7 +1401,8 @@ public void TC18() { }; - assertZProfil(expectedZ_profile, result); + + // S-R (not the rayleigh segments SO OR) assertPlanes(segmentsMeanPlanes0, propDataOut.getPropagationPaths().get(0).getSRSegment()); // Check reflexion mean planes @@ -1751,23 +1797,6 @@ public void TC23() { new Coordinate(118, 10, 8), new Coordinate(83, 10, 8)},buildingsAbs) // Ground Surface - - .addTopographicLine(30, -14, 0, 122, -14, 0)// 1 - .addTopographicLine(122, -14, 0, 122, 45, 0)// 2 - .addTopographicLine(122, 45, 0, 30, 45, 0)// 3 - .addTopographicLine(30, 45, 0, 30, -14, 0)// 4 - .addTopographicLine(59.6, -9.87, 0, 76.84, -5.28, 0)// 5 - .addTopographicLine(76.84, -5.28, 0, 63.71, 41.16, 0)// 6 - .addTopographicLine(63.71, 41.16, 0, 46.27, 36.28, 0)// 7 - .addTopographicLine(46.27, 36.28, 0, 59.6, -9.87, 0)// 8 - .addTopographicLine(46.27, 36.28, 0, 54.68, 37.59, 5)// 9 - .addTopographicLine(54.68, 37.59, 5, 55.93, 37.93, 5)// 10 - .addTopographicLine(55.93, 37.93, 5, 63.71, 41.16, 0)// 11 - .addTopographicLine(59.6, -9.87, 0, 67.35, -6.83, 5)// 12 - .addTopographicLine(67.35, -6.83, 5, 68.68, -6.49, 5)// 13 - .addTopographicLine(68.68, -6.49, 5, 76.84, -5.28, 0)// 14 - .addTopographicLine(54.68, 37.59, 5, 67.35, -6.83, 5)// 15 - .addTopographicLine(55.93, 37.93, 5, 68.68, -6.49, 5)// 16 .addGroundEffect(factory.createPolygon(new Coordinate[]{ new Coordinate(59.6, -9.87, 0), // 5 new Coordinate(76.84, -5.28, 0), // 5-6 @@ -1782,8 +1811,8 @@ public void TC23() { new Coordinate(30, 45, 0), // 7-8 new Coordinate(30, -14, 0) }), 0.); + addTopographicTC23Model(builder); builder.finishFeeding(); - //.finishFeeding(); //Propagation data building Scene rayData = new ProfileBuilderDecorator(builder) @@ -1848,6 +1877,8 @@ public void TC24() { //Create obstruction test object ProfileBuilder builder = new ProfileBuilder(); + + builder.addBuilding(new Coordinate[]{ new Coordinate(75, 34, 9), new Coordinate(110, 34, 9), @@ -1859,23 +1890,6 @@ public void TC24() { new Coordinate(118, 10, 6), new Coordinate(83, 10, 6)},buildingsAbs) // Ground Surface - - .addTopographicLine(30, -14, 0, 122, -14, 0)// 1 - .addTopographicLine(122, -14, 0, 122, 45, 0)// 2 - .addTopographicLine(122, 45, 0, 30, 45, 0)// 3 - .addTopographicLine(30, 45, 0, 30, -14, 0)// 4 - .addTopographicLine(59.6, -9.87, 0, 76.84, -5.28, 0)// 5 - .addTopographicLine(76.84, -5.28, 0, 63.71, 41.16, 0)// 6 - .addTopographicLine(63.71, 41.16, 0, 46.27, 36.28, 0)// 7 - .addTopographicLine(46.27, 36.28, 0, 59.6, -9.87, 0)// 8 - .addTopographicLine(46.27, 36.28, 0, 54.68, 37.59, 5)// 9 - .addTopographicLine(54.68, 37.59, 5, 55.93, 37.93, 5)// 10 - .addTopographicLine(55.93, 37.93, 5, 63.71, 41.16, 0)// 11 - .addTopographicLine(59.6, -9.87, 0, 67.35, -6.83, 5)// 12 - .addTopographicLine(67.35, -6.83, 5, 68.68, -6.49, 5)// 13 - .addTopographicLine(68.68, -6.49, 5, 76.84, -5.28, 0)// 14 - .addTopographicLine(54.68, 37.59, 5, 67.35, -6.83, 5)// 15 - .addTopographicLine(55.93, 37.93, 5, 68.68, -6.49, 5)// 16 .addGroundEffect(factory.createPolygon(new Coordinate[]{ new Coordinate(59.6, -9.87, 0), // 5 new Coordinate(76.84, -5.28, 0), // 5-6 @@ -1891,8 +1905,8 @@ public void TC24() { new Coordinate(30, -14, 0) }), 0.); builder.setzBuildings(true); + addTopographicTC23Model(builder); builder.finishFeeding(); - //.finishFeeding(); //Propagation data building Scene rayData = new ProfileBuilderDecorator(builder) @@ -1911,27 +1925,26 @@ public void TC24() { //Run computation computeRays.run(propDataOut); - computeRays.run(propDataOut); - - CutProfile cutProfile = computeRays.getData().profileBuilder.getProfile(rayData.sourceGeometries.get(0).getCoordinate(), rayData.receivers.get(0), computeRays.getData().gS, false); - List result = cutProfile.computePts2DGround(); - - // Expected Values /* Table 279 */ List expectedZ_profile = new ArrayList<>(); expectedZ_profile.add(new Coordinate(0.0, 0.0)); expectedZ_profile.add(new Coordinate(14.46, 0.0)); - expectedZ_profile.add(new Coordinate(19.03, 2.64)); expectedZ_profile.add(new Coordinate(23.03, 5.0)); expectedZ_profile.add(new Coordinate(24.39, 5.0)); - expectedZ_profile.add(new Coordinate(28.40, 2.65)); expectedZ_profile.add(new Coordinate(32.85, 0.0)); expectedZ_profile.add(new Coordinate(45.10, 0.0)); + expectedZ_profile.add(new Coordinate(45.10, 6.0)); + expectedZ_profile.add(new Coordinate(60.58, 6.0)); expectedZ_profile.add(new Coordinate(60.58, 0.0)); expectedZ_profile.add(new Coordinate(68.15, 0.0)); + + + List result = propDataOut.getPropagationPaths().get(0).getCutProfile().computePts2DGround(); + assertZProfil(expectedZ_profile,result); + /* Table 280 */ double [][] segmentsMeanPlanes0 = new double[][]{ // a b zs zr dp Gp Gp' @@ -1939,7 +1952,7 @@ public void TC24() { {0.0, 0.0, 6.0, 4.0, 7.57, 0.00, NaN} }; assertPlanes(segmentsMeanPlanes0,propDataOut.getPropagationPaths().get(0).getSegmentList()); - assertZProfil(expectedZ_profile,result); + }