fix freefield intersection test

noisemodelling-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/AttenuationCnossosTest.java

@@ -392,13 +392,13 @@ public void northSouthTest() {
  new Coordinate(-5, 20),
  new Coordinate(5, 20),
  new Coordinate(5, 10)
- }, 0.0)
+ }, 1.0)
  .addBuilding(new Coordinate[]{
  new Coordinate(-5, -10 ),
  new Coordinate(-5, -20 ),
  new Coordinate(5, -20),
  new Coordinate(5, -10)
- }, 0.0)
+ }, 1.0)
  .finishFeeding();
 
  //Propagation data building

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

@@ -15,7 +15,6 @@
 import java.util.ArrayList;
 import java.util.Collections;
 import java.util.List;
-import java.util.stream.Collectors;
 
 //import static org.noise_planet.noisemodelling.pathfinder.utils.geometry.JTSUtility.dist2D;
 import static org.noise_planet.noisemodelling.pathfinder.profilebuilder.ProfileBuilder.IntersectionType.*;
@@ -29,12 +28,10 @@ public class CutProfile {
  /** Receiver cut point. */
  CutPoint receiver;
  //TODO cache has intersection properties
- /** True if contains a building cutting point. */
- Boolean hasBuildingInter = false;
- /** True if contains a topography cutting point. */
- Boolean hasTopographyInter = false;
- /** True if contains a ground effect cutting point. */
- Boolean hasGroundEffectInter = false;
+ /** True if Source-Receiver linestring is below building intersection */
+ Boolean hasBuildingIntersection = false;
+ /** True if Source-Receiver linestring is below topography cutting point. */
+ Boolean hasTopographyIntersection = false;
  Boolean isFreeField;
  double distanceToSR = 0;
  Orientation srcOrientation;
@@ -69,7 +66,6 @@ public CutPoint addBuildingCutPt(Coordinate coord, int buildingId, int wallId, b
  cut.buildingId = buildingId;
  cut.wallId = wallId;
  pts.add(cut);
- hasBuildingInter = true;
  return cut;
  }
 
@@ -82,7 +78,6 @@ public CutPoint addWallCutPt(Coordinate coord, int id, boolean corner) {
  CutPoint wallPoint = new CutPoint(coord, ProfileBuilder.IntersectionType.WALL, id, corner);
  wallPoint.wallId = id;
  pts.add(wallPoint);
- hasBuildingInter = true;
  return wallPoint;
  }
 
@@ -95,7 +90,6 @@ public void addWallCutPt(Coordinate coord, int id, boolean corner, List<Double>
  pts.add(new CutPoint(coord, ProfileBuilder.IntersectionType.WALL, id, corner));
  pts.get(pts.size()-1).wallId = id;
  pts.get(pts.size()-1).setWallAlpha(alphas);
- hasBuildingInter = true;
  }
 
  /**
@@ -105,7 +99,6 @@ public void addWallCutPt(Coordinate coord, int id, boolean corner, List<Double>
  */
  public void addTopoCutPt(Coordinate coord, int id) {
  pts.add(new CutPoint(coord, TOPOGRAPHY, id));
- hasTopographyInter = true;
  }
 
  /**
@@ -125,7 +118,6 @@ public CutPoint addGroundCutPt(Coordinate coordinate, int id, double groundCoeff
  CutPoint pt = new CutPoint(coordinate, ProfileBuilder.IntersectionType.GROUND_EFFECT, id);
  pt.setGroundCoef(groundCoefficient);
  pts.add(pt);
- hasGroundEffectInter = true;
  return pt;
  }
 
@@ -192,18 +184,13 @@ public Orientation getSrcOrientation(){
  }
 
  public boolean intersectBuilding(){
- return hasBuildingInter;
+ return hasBuildingIntersection;
  }
 
  public boolean intersectTopography(){
- return hasTopographyInter;
+ return hasTopographyIntersection;
  }
 
- public boolean intersectGroundEffect(){
- return hasGroundEffectInter;
- }
-
-
  /**
  * compute the path between two points
  * @param p0
@@ -235,66 +222,21 @@ public double getGPath() {
  * @return
  */
  public boolean isFreeField() {
- return !hasBuildingInter && !hasTopographyInter;
+ return !hasBuildingIntersection && !hasTopographyIntersection;
  }
 
- /**
- * Get distance between a segment (p1,p2) and a point (point) with point perpendicular to (p1,p2)
- * @param p1
- * @param p2
- * @param point
- * @return distance in meters
- */
- private static double[] distance3D(Coordinate p1, Coordinate p2, Coordinate point) {
- double[] DistanceInfo = new double[2];
- double x1 = p1.getX();
- double y1 = p1.getY();
- double z1 = p1.getZ();
-
- double x2 = p2.getX();
- double y2 = p2.getY();
- double z2 = p2.getZ();
-
- double x0 = point.getX();
- double y0 = point.getY();
- double z0 = point.getZ();
-
- // Vector representing the LineSegment
- double dx = x2 - x1;
- double dy = y2 - y1;
- double dz = z2 - z1;
-
- // Vector from the start point of the LineSegment to the Point
- double px = x0 - x1;
- double py = y0 - y1;
- double pz = z0 - z1;
-
- // Compute the dot product of the vectors
- double dotProduct = dx * px + dy * py + dz * pz;
-
- // Calculate the projection of the Point onto the LineSegment
- double t = dotProduct / (dx * dx + dy * dy + dz * dz);
-
- // Calculate the closest point on the LineSegment to the Point
- double closestX = x1 + t * dx;
- double closestY = y1 + t * dy;
- double closestZ = z1 + t * dz;
-
- // Calculate the distance between the closest point and the Point
- double distance = Math.sqrt((x0 - closestX) * (x0 - closestX)
- + (y0 - closestY) * (y0 - closestY)
- + (z0 - closestZ) * (z0 - closestZ));
- double sign = z0 - closestZ;
- DistanceInfo[0]=distance;
- DistanceInfo[1]=sign;
- return DistanceInfo;
- }
 
  @Override
  public String toString() {
- return "CutProfile{" + "pts=" + pts + ", source=" + source + ", receiver=" + receiver + ", " +
- "hasBuildingInter=" + hasBuildingInter + ", hasTopographyInter=" + hasTopographyInter + ", " +
- "hasGroundEffectInter=" + hasGroundEffectInter + ", isFreeField=" + isFreeField + ", " +
- "srcOrientation=" + srcOrientation + '}';
+ return "CutProfile{" +
+ "pts=" + pts +
+ ", source=" + source +
+ ", receiver=" + receiver +
+ ", hasBuildingIntersection=" + hasBuildingIntersection +
+ ", hasTopographyIntersection=" + hasTopographyIntersection +
+ ", isFreeField=" + isFreeField +
+ ", distanceToSR=" + distanceToSR +
+ ", srcOrientation=" + srcOrientation +
+ '}';
  }
 }

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

@@ -996,7 +996,7 @@ public CutProfile getProfile(Coordinate sourceCoordinate, Coordinate receiverCoo
  //Fetch topography evolution between sourceCoordinate and receiverCoordinate
  if(topoTree != null) {
  addTopoCutPts(sourceCoordinate, receiverCoordinate, profile, stopAtObstacleOverSourceReceiver);
- if(stopAtObstacleOverSourceReceiver && profile.hasTopographyInter) {
+ if(stopAtObstacleOverSourceReceiver && profile.hasTopographyIntersection) {
  return profile;
  }
  } else {
@@ -1008,7 +1008,7 @@ public CutProfile getProfile(Coordinate sourceCoordinate, Coordinate receiverCoo
  if(rtree != null) {
  LineSegment fullLine = new LineSegment(sourceCoordinate, receiverCoordinate);
  addGroundBuildingCutPts(fullLine, profile, stopAtObstacleOverSourceReceiver);
- if(stopAtObstacleOverSourceReceiver && profile.hasBuildingInter) {
+ if(stopAtObstacleOverSourceReceiver && profile.hasBuildingIntersection) {
  return profile;
  }
  }
@@ -1046,6 +1046,11 @@ public CutProfile getProfile(Coordinate sourceCoordinate, Coordinate receiverCoo
  new Coordinate(previousZGround.coordinate.x, previousZGround.coordinate.y,
  previousZGround.getzGround()),
  new Coordinate(nextPoint.coordinate.x, nextPoint.coordinate.y, nextPoint.getzGround()));
+ if(Double.isNaN(cutPoint.coordinate.z)) {
+ // Bottom of walls are set to NaN z because it can be computed here at low cost
+ // (without fetch dem r-tree)
+ cutPoint.coordinate.setZ(cutPoint.zGround);
+ }
  } else {
  // we have an update on Z ground
  previousZGround = cutPoint;
@@ -1083,6 +1088,7 @@ public int getIntersectingGroundAbsorption(Geometry query) {
  */
  private void addGroundBuildingCutPts(LineSegment fullLine, CutProfile profile, boolean stopAtObstacleOverSourceReceiver) {
  Vector2D directionAfter = Vector2D.create(fullLine.p0, fullLine.p1).normalize().multiply(MILLIMETER);
+ Vector2D directionBefore = directionAfter.negate();
  // Collect all objects where envelope intersects all sub-segments of fullLine
  Set<Integer> processed = new HashSet<>();
 
@@ -1118,15 +1124,21 @@ private void addGroundBuildingCutPts(LineSegment fullLine, CutProfile profile, b
  Vector2D exteriorVector = facetVector.rotate(LEFT_SIDE).normalize().multiply(MILLIMETER);
  Coordinate exteriorPoint = exteriorVector.add(Vector2D.create(intersection)).toCoordinate();
  CutPoint exteriorPointCutPoint = profile.addBuildingCutPt(exteriorPoint, facetLine.originId, i, false);
- if (topoTree == null) {
- exteriorPointCutPoint.coordinate.setZ(0.0);
- } else {
- exteriorPointCutPoint.coordinate.setZ(getZGround(exteriorPointCutPoint));
- pt.zGround = exteriorPointCutPoint.coordinate.z;
- exteriorPointCutPoint.zGround = exteriorPointCutPoint.coordinate.z;
+ exteriorPointCutPoint.coordinate.setZ(NaN);
+ double zRayReceiverSource = Vertex.interpolateZ(intersection,fullLine.p0, fullLine.p1);
+ if(zRayReceiverSource <= intersection.z) {
+ profile.hasBuildingIntersection = true;
  }
  } else if (facetLine.type == IntersectionType.WALL) {
+ profile.addWallCutPt(Vector2D.create(intersection).add(directionBefore).toCoordinate(),
+ facetLine.originId, false, facetLine.alphas);
  profile.addWallCutPt(intersection, facetLine.originId, false, facetLine.alphas);
+ profile.addWallCutPt(Vector2D.create(intersection).add(directionAfter).toCoordinate(),
+ facetLine.originId, false, facetLine.alphas);
+ double zRayReceiverSource = Vertex.interpolateZ(intersection,fullLine.p0, fullLine.p1);
+ if(zRayReceiverSource <= intersection.z) {
+ profile.hasBuildingIntersection = true;
+ }
  } else if (facetLine.type == GROUND_EFFECT) {
  // we hit the border of a ground effect
  // we need to add a new point with the new value of the ground effect
@@ -1283,7 +1295,7 @@ public int getTriangleIdByCoordinate(Coordinate pt) {
  public void addTopoCutPts(Coordinate p1, Coordinate p2, CutProfile profile, boolean stopAtObstacleOverSourceReceiver) {
  List<Coordinate> coordinates = new ArrayList<>();
  boolean freeField = fetchTopographicProfile(coordinates, p1, p2, stopAtObstacleOverSourceReceiver);
- profile.hasTopographyInter = !freeField;
+ profile.hasTopographyIntersection = !freeField;
  // Remove unnecessary points
  ArrayList<Coordinate> retainedCoordinates = new ArrayList<>(coordinates.size());
  for(int i =0; i < coordinates.size(); i++) {

noisemodelling-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/utils/geometry/JTSUtility.java

@@ -34,6 +34,59 @@ public class JTSUtility {
  **/
  private JTSUtility() {}
 
+
+ /**
+ * Get distance between a segment (p1,p2) and a point (point) with point perpendicular to (p1,p2)
+ * @param p1 origin segment
+ * @param p2 destination segment
+ * @param point reference point
+ * @return DistanceInfo[0]=distance; DistanceInfo[1]=sign;
+ */
+ private static double[] distance3D(Coordinate p1, Coordinate p2, Coordinate point) {
+ double[] DistanceInfo = new double[2];
+ double x1 = p1.getX();
+ double y1 = p1.getY();
+ double z1 = p1.getZ();
+
+ double x2 = p2.getX();
+ double y2 = p2.getY();
+ double z2 = p2.getZ();
+
+ double x0 = point.getX();
+ double y0 = point.getY();
+ double z0 = point.getZ();
+
+ // Vector representing the LineSegment
+ double dx = x2 - x1;
+ double dy = y2 - y1;
+ double dz = z2 - z1;
+
+ // Vector from the start point of the LineSegment to the Point
+ double px = x0 - x1;
+ double py = y0 - y1;
+ double pz = z0 - z1;
+
+ // Compute the dot product of the vectors
+ double dotProduct = dx * px + dy * py + dz * pz;
+
+ // Calculate the projection of the Point onto the LineSegment
+ double t = dotProduct / (dx * dx + dy * dy + dz * dz);
+
+ // Calculate the closest point on the LineSegment to the Point
+ double closestX = x1 + t * dx;
+ double closestY = y1 + t * dy;
+ double closestZ = z1 + t * dz;
+
+ // Calculate the distance between the closest point and the Point
+ double distance = Math.sqrt((x0 - closestX) * (x0 - closestX)
+ + (y0 - closestY) * (y0 - closestY)
+ + (z0 - closestZ) * (z0 - closestZ));
+ double sign = z0 - closestZ;
+ DistanceInfo[0]=distance;
+ DistanceInfo[1]=sign;
+ return DistanceInfo;
+ }
+
  /**
  * Compute a and b linear function of the line p1 p2
  * @param p1 p1