Skip to content

Commit

Permalink
Repair major diffraction issues (horizontal through berms or topograp…
Browse files Browse the repository at this point in the history 
…hy and vertical diffraction was sometimes not computed)
  • Loading branch information
@pierromond
pierromond committed Jul 20, 2023
1 parent 274b4c4 commit e2526fb
Show file tree
Hide file tree
Showing 3 changed files with 159 additions and 44 deletions.
37 changes: 37 additions & 0 deletions ...bc/src/test/java/org/noise_planet/noisemodelling/jdbc/EvaluateAttenuationCnossosTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5380,6 +5380,43 @@ public void TC28() {
computeRays.setThreadCount(1);
computeRays.run(propDataOut);

//Expected values
//Path0 : vertical plane
double[] expectedAlphaAtm = new double[]{0.12, 0.41, 1.04, 1.93, 3.66, 9.66, 32.77, 116.88};
double[] expectedAAtm = new double[]{0.12, 0.41, 1.04, 1.93, 3.66, 9.68, 32.81, 117.03};
double[] expectedADiv = new double[]{71.01, 71.01, 71.01, 71.01, 71.01, 71.01, 71.01, 71.01};
double[] expectedABoundaryH = new double[]{14.85, 17.56, 20.52, 22.31, 22.31, 22.31, 22.32, 22.32};
double[] expectedABoundaryF = new double[]{7.48, 10.11, 12.94, 15.86, 18.82, 19.78, 19.78, 19.78};
double[] expectedLH = new double[]{64.02, 61.02, 57.43, 54.75, 53.01, 47.00, 23.86, -60.35};
double[] expectedLF = new double[]{71.39, 68.46, 65.00, 61.20, 56.51, 49.54, 26.40, -57.82};
double[] expectedL = new double[]{69.11, 66.17, 62.69, 59.08, 55.10, 48.45, 25.31, -58.90};
double[] expectedLA = new double[]{42.91, 50.07, 54.09, 55.88, 55.10, 49.65, 26.31, -60.00};

PropagationPath proPath = propDataOut.getPropagationPaths().get(0);

double[] actualAlphaAtm = propDataOut.genericMeteoData.getAlpha_atmo();
double[] actualAAtm = proPath.absorptionData.aAtm;
double[] actualADiv = proPath.absorptionData.aDiv;
double[] actualABoundaryH = proPath.absorptionData.aBoundaryH;
double[] actualABoundaryF = proPath.absorptionData.aBoundaryF;
double[] actualLH = addArray(proPath.absorptionData.aGlobalH, SOUND_POWER_LEVELS);
double[] actualLF = addArray(proPath.absorptionData.aGlobalF, SOUND_POWER_LEVELS);
double[] actualL = addArray(proPath.absorptionData.aGlobal, SOUND_POWER_LEVELS);
double[] actualLA = addArray(actualL, A_WEIGHTING);

/* assertDoubleArrayEquals("AlphaAtm - vertical plane", expectedAlphaAtm, actualAlphaAtm, ERROR_EPSILON_LOWEST);
assertDoubleArrayEquals("AAtm - vertical plane", expectedAAtm, actualAAtm, ERROR_EPSILON_LOWEST);
assertDoubleArrayEquals("ADiv - vertical plane", expectedADiv, actualADiv, ERROR_EPSILON_LOWEST);
assertDoubleArrayEquals("ABoundaryH - vertical plane", expectedABoundaryH, actualABoundaryH, ERROR_EPSILON_VERY_HIGH);
assertDoubleArrayEquals("ABoundaryF - vertical plane", expectedABoundaryF, actualABoundaryF, ERROR_EPSILON_VERY_HIGH);
assertDoubleArrayEquals("LH - vertical plane", expectedLH, actualLH, ERROR_EPSILON_VERY_HIGH);
assertDoubleArrayEquals("LF - vertical plane", expectedLF, actualLF, ERROR_EPSILON_VERY_HIGH);
assertDoubleArrayEquals("L - vertical plane", expectedL, actualL, ERROR_EPSILON_VERY_HIGH);
assertDoubleArrayEquals("LA - vertical plane", expectedLA, actualLA, ERROR_EPSILON_VERY_HIGH);*/




double[] L = addArray(propDataOut.getVerticesSoundLevel().get(0).value, new double[]{150-26.2,150-16.1,150-8.6,150-3.2,150,150+1.2,150+1.0,150-1.1});
assertArrayEquals( new double[]{43.56,50.59,54.49,56.14,55.31,49.77,23.37,-59.98},L, ERROR_EPSILON_VERY_HIGH);

Expand Down
104 changes: 61 additions & 43 deletions ...thfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/ComputeCnossosRays.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@
import org.locationtech.jts.index.ItemVisitor;
import org.locationtech.jts.math.Vector2D;
import org.locationtech.jts.math.Vector3D;
import org.locationtech.jts.simplify.DouglasPeuckerSimplifier;
import org.locationtech.jts.triangulate.quadedge.Vertex;
import org.noise_planet.noisemodelling.pathfinder.utils.ProfilerThread;
import org.noise_planet.noisemodelling.pathfinder.utils.ReceiverStatsMetric;
Expand All @@ -55,7 +56,6 @@
import java.util.concurrent.atomic.AtomicInteger;
import java.util.stream.Collectors;

import static java.lang.Double.isInfinite;
import static java.lang.Double.isNaN;
import static java.lang.Math.*;
import static org.noise_planet.noisemodelling.pathfinder.ComputeCnossosRays.ComputationSide.LEFT;
Expand Down Expand Up @@ -617,6 +617,10 @@ public PropagationPath computeVEdgeDiffraction(Coordinate rcvCoord, Coordinate s
for(int i=0; i<coordinates.size()-1; i++) {
ProfileBuilder.CutProfile profile = data.profileBuilder.getProfile(coordinates.get(i), coordinates.get(i+1), data.gS);
profile.setSrcOrientation(orientation);
if (!profile.isFreeField()){
path = null;
return path;
}
double dist = dist2D(coordinates.get(i), coordinates.get(i+1));
g+=profile.getGPath()*dist;
d+=dist;
Expand Down Expand Up @@ -737,22 +741,25 @@ public PropagationPath computeHEdgeDiffraction(ProfileBuilder.CutProfile cutProf
if(pts2D.size() != cutPts.size()) {
throw new IllegalArgumentException("The two arrays size should be the same");
}
//Remove aligned cut points
//Remove aligned cut points thanks to jts DouglasPeuckerSimplifier algo
List<ProfileBuilder.CutPoint> newCutPts = new ArrayList<>(cutPts.size());
List<Coordinate> newPts2D = new ArrayList<>(pts2D.size());
newCutPts.add(cutPts.get(0));
newPts2D.add(pts2D.get(0));
for(int i=0; i<pts2D.size()-2; i++) {
Coordinate c0 = pts2D.get(i);
Coordinate c1 = pts2D.get(i+1);
Coordinate c2 = pts2D.get(i+2);
if(new LineSegment(c0, c2).distance(c1) >= 0.1) {
newPts2D.add(c1);
newCutPts.add(cutPts.get(i+1));
Geometry lineString = new GeometryFactory().createLineString(pts2D.toArray(new Coordinate[0]));
List<Coordinate> newPts2D = List.of(DouglasPeuckerSimplifier.simplify(lineString, 0.5*cutProfile.getDistanceToSR()).getCoordinates());
List<Integer> matchingIndices = new ArrayList<>();
for (int i = 0; i < pts2D.size(); i++) {
Coordinate coordinate1 = pts2D.get(i);
for (int j = 0; j < newPts2D.size(); j++) {
Coordinate coordinate2 = newPts2D.get(j);
if (coordinate1.equals(coordinate2)) {
matchingIndices.add(i);
break; // Move to the next coordinate in list1 once a match is found
}
}
}
newPts2D.add(pts2D.get(pts2D.size()-1));
newCutPts.add(cutPts.get(cutPts.size() - 1));
for (int index : matchingIndices) {
newCutPts.add(cutPts.get(index));
}

pts2D = newPts2D;
cutPts = newCutPts;
double[] meanPlane = JTSUtility.getMeanPlaneCoefficients(pts2D.toArray(new Coordinate[0]));
Expand All @@ -765,39 +772,50 @@ public PropagationPath computeHEdgeDiffraction(ProfileBuilder.CutProfile cutProf
propagationPath.setCutPoints(cutPts);
LineSegment srcRcvLine = new LineSegment(firstPts2D, lastPts2D);
List<Coordinate> pts = new ArrayList<>();
pts.add(firstPts2D);
for(int i=1; i<pts2D.size(); i++) {
Coordinate pt = pts2D.get(i);
double frac = srcRcvLine.segmentFraction(pt);
double y = -Double.MAX_VALUE;
for(int j=i+1; j<pts2D.size(); j++) {
y = max(y, srcRcvLine.p0.y + frac*(pts2D.get(j).y-srcRcvLine.p0.y));

// Extract the first and last points to define the line segment
Coordinate firstPt = pts2D.get(0);
Coordinate lastPt = pts2D.get(pts2D.size() - 1);

// Compute the slope and y-intercept of the line segment
double slope = (lastPt.y - firstPt.y) / (lastPt.x - firstPt.x);
double yIntercept = firstPt.y - slope * firstPt.x;

// Filter out points that are below the line segment
List<Coordinate> filteredCoordinates = new ArrayList<>();
for (Coordinate coord : pts2D) {
double lineY = slope * coord.x + yIntercept;
if (coord.y >= lineY-0.001) { // espilon to avoir float issues
filteredCoordinates.add(coord);
}
if(y <= pt.y){
pts.add(pt);
srcRcvLine = new LineSegment(pt, lastPts2D);

//Filter point to only keep hull.
List<Coordinate> toRemove = new ArrayList<>();
//check if last-1 point is under or not the surrounding points
for(int j = pts.size()-2; j > 0; j--) {
Coordinate jPt = pts.get(j);
if(jPt.y==Double.MAX_VALUE || isInfinite(jPt.y)) {
toRemove.add(jPt);
}
//line between last point and previous-1 point
else {
LineSegment lineRm = new LineSegment(pts.get(j - 1), pt);
double fracRm = lineRm.segmentFraction(jPt);
double zRm = lineRm.p0.z + fracRm * (lineRm.p1.z - lineRm.p0.z);
if (zRm >= jPt.z) {
toRemove.add(jPt);
}
}
}

// Compute the convex hull using JTS
GeometryFactory geomFactory = new GeometryFactory();
Coordinate[] coordsArray = filteredCoordinates.toArray(new Coordinate[0]);
ConvexHull convexHull = new ConvexHull(coordsArray, geomFactory);
Coordinate[] convexHullCoords = convexHull.getConvexHull().getCoordinates();
Arrays.sort(convexHullCoords, Comparator.comparingDouble(coord -> coord.x));

CoordinateSequence coordSequence = geomFactory.getCoordinateSequenceFactory().create(convexHullCoords);
Geometry geom = geomFactory.createLineString(coordSequence);
Geometry uniqueGeom = geom.union(); // Removes duplicate coordinates
convexHullCoords = uniqueGeom.getCoordinates();

// Convert the result back to your format (List<Point2D> pts)
List<Coordinate> convexHullPoints = new ArrayList<>();
if (convexHullCoords.length ==3){
convexHullPoints = Arrays.asList(convexHullCoords);
}else {
for (int j = 0; j < convexHullCoords.length; j++) {
// Check if the y-coordinate is valid (not equal to Double.MAX_VALUE and not infinite)
if (convexHullCoords[j].y == Double.MAX_VALUE || Double.isInfinite(convexHullCoords[j].y)) {
continue; // Skip this point as it's not part of the hull
}
pts.removeAll(toRemove);
convexHullPoints.add(convexHullCoords[j]);
}
}
pts = convexHullPoints;

double e = 0;
Coordinate src = null;
Expand Down
62 changes: 61 additions & 1 deletion ...g-pathfinder/src/main/java/org/noise_planet/noisemodelling/pathfinder/ProfileBuilder.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1647,6 +1647,7 @@ public static class CutProfile {
/** True if contains a ground effect cutting point. */
private Boolean hasGroundEffectInter = false;
private Boolean isFreeField;
private double distanceToSR = 0;
private Orientation srcOrientation;

/**
Expand Down Expand Up @@ -1747,6 +1748,11 @@ public CutPoint getSource() {
return source;
}

/**
* get Distance of the not free field point to the Source-Receiver Segement
* @return
*/
public double getDistanceToSR(){return distanceToSR;}
/**
* Retrieve the profile receiver.
* @return The profile receiver.
Expand Down Expand Up @@ -1869,15 +1875,69 @@ else if(cut.getType() == TOPOGRAPHY) {
for(CutPoint pt : ptsWithouGroundEffect) {
double frac = (pt.coordinate.x-s.x)/(r.x-s.x);
double z = source.getCoordinate().z + frac * (receiver.getCoordinate().z-source.getCoordinate().z);
if(z < pt.getCoordinate().z && !pt.isCorner()) {
double[] distanceSRpt = distance3D(source.getCoordinate(), receiver.getCoordinate(), pt.getCoordinate());
if(distanceSRpt[0]>0 && distanceSRpt[1]>0 && !pt.isCorner()) {
isFreeField = false;
distanceToSR = distanceSRpt[0];
break;
}
}
}
return isFreeField;
}

/**
* 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 + ", " +
Expand Down

0 comments on commit e2526fb

Please sign in to comment.