fix tc23 and tc24 dem as topographic lines over beam was not parallel…

… so artifacts were generated. The cause of non parallel is rounding on line table in the cnossos document.

noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/profilebuilder/CutProfile.java

@@ -277,7 +277,8 @@ public List<Coordinate> 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())) {

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<Coordinate> 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<Coordinate> 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,19 +1422,35 @@ public boolean fetchTopographicProfile(List<Coordinate> 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;
                         }
                     }
                 }
             }
             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));
     }
 
     /**

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<Coordinate> getNewCoordinateSystem(List<Coordinate> 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<Coordinate> getNewCoordinateSystem(List<Coordinate> listPoints, double tolerance) {
         if(listPoints.isEmpty()) {
             return new ArrayList<>();
         }
@@ -245,6 +266,20 @@ public static List<Coordinate> getNewCoordinateSystem(List<Coordinate> 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;
     }