Skip to content

Commit

Permalink
Check contouring keeping Z value
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolas-f
nicolas-f committed Oct 5, 2023
1 parent ec29b63 commit 2159475
Show file tree
Hide file tree
Showing 3 changed files with 162 additions and 63 deletions.
59 changes: 47 additions & 12 deletions noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/BezierContouring.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
*/
package org.noise_planet.noisemodelling.jdbc;

import org.h2gis.utilities.GeometryTableUtilities;
import org.h2gis.utilities.JDBCUtilities;
import org.h2gis.utilities.TableLocation;
import org.h2gis.utilities.jts_utils.Contouring;
Expand Down Expand Up @@ -53,6 +54,8 @@ public class BezierContouring {
List<Double> isoLevels;
List<String> isoLabels;
boolean smooth = true;

boolean mergeTriangles = true;
double smoothCoefficient = 1.0;
double deltaPoints = 0.5; // minimal distance between bezier points
double epsilon = 0.05;
Expand Down Expand Up @@ -96,6 +99,20 @@ public double getSmoothCoefficient() {
return smoothCoefficient;
}

/**
* @return True if triangles will be merged using isolevel attribute
*/
public boolean isMergeTriangles() {
return mergeTriangles;
}

/**
* @param mergeTriangles True if triangles will be merged using isolevel attribute, Z ordinate will be lost
*/
public void setMergeTriangles(boolean mergeTriangles) {
this.mergeTriangles = mergeTriangles;
}

/**
* @param smoothCoefficient Coefficient of polygons smoothing [0-1]
*/
Expand Down Expand Up @@ -181,11 +198,11 @@ static List<Coordinate> curve4(Coordinate anchor1, Coordinate control1, Coordina

static Coordinate[] generateBezierCurves(Coordinate[] coordinates, Quadtree segmentTree, double pointsDelta) {
ArrayList<Coordinate> pts = new ArrayList<>();
pts.add(coordinates[0]);
pts.add(new Coordinate(coordinates[0].x, coordinates[0].y));
for(int i = 0; i < coordinates.length - 1; i++) {
final int i2 = i + 1;
Coordinate p1 = coordinates[i];
Coordinate p2 = coordinates[i2];
Coordinate p1 = new Coordinate(coordinates[i].x, coordinates[i].y);
Coordinate p2 = new Coordinate(coordinates[i2].x, coordinates[i2].y);

Segment segment = new Segment(p1, p2);
List<Segment> segments = (List<Segment>)segmentTree.query(segment.getEnvelope());
Expand Down Expand Up @@ -409,7 +426,7 @@ void processCell(Connection connection, int cellId, Map<Short, ArrayList<Geometr
+ "(cell_id, the_geom, ISOLVL, ISOLABEL) VALUES (?, ?, ?, ?);")) {
for (Map.Entry<Short, ArrayList<Geometry>> entry : polys.entrySet()) {
ArrayList<Polygon> polygons = new ArrayList<>();
if(!smooth) {
if(!smooth && mergeTriangles) {
// Merge triangles
try {
CascadedPolygonUnion union = new CascadedPolygonUnion(entry.getValue());
Expand Down Expand Up @@ -454,12 +471,21 @@ public void createTable(Connection connection) throws SQLException {
Map<Short, ArrayList<Geometry>> polyMap = new HashMap<>();
int lastCellId = -1;
try(Statement st = connection.createStatement()) {
String geometryType = "GEOMETRY(POLYGONZ,"+srid+")";
if(smooth || mergeTriangles) {
geometryType = "GEOMETRY(POLYGON,"+srid+")";
}
st.execute("DROP TABLE IF EXISTS " + TableLocation.parse(outputTable));
st.execute("CREATE TABLE " + TableLocation.parse(outputTable) + "(PK SERIAL, CELL_ID INTEGER, THE_GEOM GEOMETRY, ISOLVL INTEGER, ISOLABEL VARCHAR);");
String query = "SELECT CELL_ID, ST_X(p1.the_geom) xa,ST_Y(p1.the_geom) ya,ST_X(p2.the_geom) xb,ST_Y(p2.the_geom) yb,ST_X(p3.the_geom) xc,ST_Y(p3.the_geom) yc, p1."+pointTableField+" lvla, p2."+pointTableField+" lvlb, p3."+pointTableField+" lvlc FROM "+triangleTable+" t, "+pointTable+" p1,"+pointTable+" p2,"+pointTable+" p3 WHERE t.PK_1 = p1."+pkField+" and t.PK_2 = p2."+pkField+" AND t.PK_3 = p3."+pkField+" order by cell_id;";
st.execute("CREATE TABLE " + TableLocation.parse(outputTable) +
"(PK SERIAL, CELL_ID INTEGER, THE_GEOM "+geometryType+", ISOLVL INTEGER, ISOLABEL VARCHAR);");
String query = "SELECT CELL_ID, ST_X(p1.the_geom) xa,ST_Y(p1.the_geom) ya, ST_Z(p1.the_geom) za," +
"ST_X(p2.the_geom) xb,ST_Y(p2.the_geom) yb, ST_Z(p2.the_geom) zb," +
"ST_X(p3.the_geom) xc,ST_Y(p3.the_geom) yc, ST_Z(p3.the_geom) zc," +
" p1."+pointTableField+" lvla, p2."+pointTableField+" lvlb, p3."+pointTableField+" lvlc FROM "+triangleTable+" t, "+pointTable+" p1,"+pointTable+" p2,"+pointTable+" p3 WHERE t.PK_1 = p1."+pkField+" and t.PK_2 = p2."+pkField+" AND t.PK_3 = p3."+pkField+" order by cell_id;";
try(ResultSet rs = st.executeQuery(query)) {
// Cache columns index
int xa = 0, xb = 0, xc = 0, ya = 0, yb = 0, yc = 0, lvla = 0, lvlb = 0, lvlc = 0, cell_id = 0;
int xa = 0, xb = 0, xc = 0, ya = 0, yb = 0, yc = 0, za = 0, zb = 1, zc = 1, lvla = 0, lvlb = 0,
lvlc = 0, cell_id = 0;
ResultSetMetaData resultSetMetaData = rs.getMetaData();
for (int columnId = 1; columnId <= resultSetMetaData.getColumnCount(); columnId++) {
switch (resultSetMetaData.getColumnLabel(columnId).toUpperCase()) {
Expand All @@ -481,6 +507,15 @@ public void createTable(Connection connection) throws SQLException {
case "YC":
yc = columnId;
break;
case "ZA":
za = columnId;
break;
case "ZB":
zb = columnId;
break;
case "ZC":
zc = columnId;
break;
case "LVLA":
lvla = columnId;
break;
Expand All @@ -495,8 +530,8 @@ public void createTable(Connection connection) throws SQLException {
break;
}
}
if (xa == 0 || xb == 0 || xc == 0 || ya == 0 || yb == 0 || yc == 0 || lvla == 0 || lvlb == 0 ||
lvlc == 0 || cell_id == 0) {
if (xa == 0 || xb == 0 || xc == 0 || ya == 0 || yb == 0 || yc == 0 || za == 0 || zb == 0 || zc == 0
|| lvla == 0 || lvlb == 0 || lvlc == 0 || cell_id == 0) {
throw new SQLException("Missing field in input tables");
}
while(rs.next()) {
Expand All @@ -508,9 +543,9 @@ public void createTable(Connection connection) throws SQLException {
}
lastCellId = cellId;
// Split current triangle
Coordinate a = new Coordinate(rs.getDouble(xa), rs.getDouble(ya));
Coordinate b = new Coordinate(rs.getDouble(xb), rs.getDouble(yb));
Coordinate c = new Coordinate(rs.getDouble(xc), rs.getDouble(yc));
Coordinate a = new Coordinate(rs.getDouble(xa), rs.getDouble(ya), rs.getDouble(za));
Coordinate b = new Coordinate(rs.getDouble(xb), rs.getDouble(yb), rs.getDouble(zb));
Coordinate c = new Coordinate(rs.getDouble(xc), rs.getDouble(yc), rs.getDouble(zc));
// Fetch data
TriMarkers triMarkers = new TriMarkers(a, b, c, dbaToW(rs.getDouble(lvla)),
dbaToW(rs.getDouble(lvlb)),
Expand Down
103 changes: 56 additions & 47 deletions noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/TriangleNoiseMap.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -299,6 +299,60 @@ public void setEpsilon(double epsilon) {
this.epsilon = epsilon;
}

public static void generateResultTable(Connection connection, String receiverTableName, String trianglesTableName,
AtomicInteger receiverPK, List<Coordinate> vertices,
GeometryFactory geometryFactory, List<Triangle> triangles, int cellI,
int cellJ, int gridDim) throws SQLException {

if(!JDBCUtilities.tableExists(connection, receiverTableName)) {
Statement st = connection.createStatement();
st.execute("CREATE TABLE "+TableLocation.parse(receiverTableName)+"(pk serial NOT NULL, the_geom geometry not null, PRIMARY KEY (PK))");
}
if(!JDBCUtilities.tableExists(connection, trianglesTableName)) {
Statement st = connection.createStatement();
st.execute("CREATE TABLE "+TableLocation.parse(trianglesTableName)+"(pk serial NOT NULL, the_geom geometry , PK_1 integer not null, PK_2 integer not null, PK_3 integer not null, cell_id integer not null, PRIMARY KEY (PK))");
}
int receiverPkOffset = receiverPK.get();
// Add vertices to receivers
PreparedStatement ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(receiverTableName)+" VALUES (?, ?);");
int batchSize = 0;
for(Coordinate v : vertices) {
ps.setInt(1, receiverPK.getAndAdd(1));
ps.setObject(2, geometryFactory.createPoint(v));
ps.addBatch();
batchSize++;
if (batchSize >= BATCH_MAX_SIZE) {
ps.executeBatch();
ps.clearBatch();
batchSize = 0;
}
}
if (batchSize > 0) {
ps.executeBatch();
}
// Add triangles
ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(trianglesTableName)+"(the_geom, PK_1, PK_2, PK_3, CELL_ID) VALUES (?, ?, ?, ?, ?);");
batchSize = 0;
for(Triangle t : triangles) {
ps.setObject(1, geometryFactory.createPolygon(new Coordinate[]{vertices.get(t.getA()),
vertices.get(t.getB()), vertices.get(t.getC()), vertices.get(t.getA())}));
ps.setInt(2, t.getA() + receiverPkOffset);
ps.setInt(3, t.getC() + receiverPkOffset);
ps.setInt(4, t.getB() + receiverPkOffset);
ps.setInt(5, cellI * gridDim + cellJ);
ps.addBatch();
batchSize++;
if (batchSize >= BATCH_MAX_SIZE) {
ps.executeBatch();
ps.clearBatch();
batchSize = 0;
}
}
if (batchSize > 0) {
ps.executeBatch();
}
}

public void generateReceivers(Connection connection, int cellI, int cellJ, String receiverTableName, String trianglesTableName, AtomicInteger receiverPK) throws SQLException, LayerDelaunayError, IOException {

int ij = cellI * gridDim + cellJ + 1;
Expand Down Expand Up @@ -367,53 +421,8 @@ public void generateReceivers(Connection connection, int cellI, int cellJ, Strin
}
nbreceivers += vertices.size();

if(!JDBCUtilities.tableExists(connection, receiverTableName)) {
Statement st = connection.createStatement();
st.execute("CREATE TABLE "+TableLocation.parse(receiverTableName)+"(pk serial NOT NULL, the_geom geometry not null, PRIMARY KEY (PK))");
}
if(!JDBCUtilities.tableExists(connection, trianglesTableName)) {
Statement st = connection.createStatement();
st.execute("CREATE TABLE "+TableLocation.parse(trianglesTableName)+"(pk serial NOT NULL, the_geom geometry , PK_1 integer not null, PK_2 integer not null, PK_3 integer not null, cell_id integer not null, PRIMARY KEY (PK))");
}
int receiverPkOffset = receiverPK.get();
// Add vertices to receivers
PreparedStatement ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(receiverTableName)+" VALUES (?, ?);");
int batchSize = 0;
for(Coordinate v : vertices) {
ps.setInt(1, receiverPK.getAndAdd(1));
ps.setObject(2, geometryFactory.createPoint(v));
ps.addBatch();
batchSize++;
if (batchSize >= BATCH_MAX_SIZE) {
ps.executeBatch();
ps.clearBatch();
batchSize = 0;
}
}
if (batchSize > 0) {
ps.executeBatch();
}
// Add triangles
ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(trianglesTableName)+"(the_geom, PK_1, PK_2, PK_3, CELL_ID) VALUES (?, ?, ?, ?, ?);");
batchSize = 0;
for(Triangle t : triangles) {
ps.setObject(1, geometryFactory.createPolygon(new Coordinate[]{vertices.get(t.getA()),
vertices.get(t.getB()), vertices.get(t.getC()), vertices.get(t.getA())}));
ps.setInt(2, t.getA() + receiverPkOffset);
ps.setInt(3, t.getC() + receiverPkOffset);
ps.setInt(4, t.getB() + receiverPkOffset);
ps.setInt(5, cellI * gridDim + cellJ);
ps.addBatch();
batchSize++;
if (batchSize >= BATCH_MAX_SIZE) {
ps.executeBatch();
ps.clearBatch();
batchSize = 0;
}
}
if (batchSize > 0) {
ps.executeBatch();
}
generateResultTable(connection, receiverTableName, trianglesTableName, receiverPK, vertices, geometryFactory,
triangles, cellI, cellJ, gridDim);
}

public double getRoadWidth() {
Expand Down
63 changes: 59 additions & 4 deletions ...ing-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/BezierContouringJDBCTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,28 @@
import org.h2gis.functions.io.geojson.GeoJsonRead;
import org.h2gis.functions.io.shp.SHPRead;
import org.h2gis.functions.io.shp.SHPWrite;
import org.h2gis.functions.spatial.mesh.DelaunayData;
import org.h2gis.utilities.JDBCUtilities;
import org.h2gis.utilities.SpatialResultSet;
import org.h2gis.utilities.TableLocation;
import org.h2gis.utilities.TableUtilities;
import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.noise_planet.noisemodelling.pathfinder.LayerDelaunayError;
import org.noise_planet.noisemodelling.pathfinder.LayerTinfour;

import java.io.IOException;
import java.sql.Connection;
import java.sql.SQLException;
import java.sql.Statement;
import java.nio.file.Paths;
import java.sql.*;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.concurrent.atomic.AtomicInteger;

import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;

public class BezierContouringJDBCTest {
Expand Down Expand Up @@ -70,6 +79,52 @@ public void testBezierContouring() throws SQLException, IOException {
assertTrue(fieldValues.contains("8"));
assertTrue(fieldValues.contains("9"));

SHPWrite.exportTable(connection, "target/contouring.shp", "CONTOURING_NOISE_MAP","UTF-8",true);
}

@Test
public void testContouring3D() throws SQLException, IOException, LayerDelaunayError {
// Will create elevation iso from DEM table
GeoJsonRead.importTable(connection, Paths.get(Paths.get(System.getProperty("user.dir")).getParent().toString(),
"wps_scripts/src/test/resources/org/noise_planet/noisemodelling/wps/dem.geojson").toString());
LayerTinfour delaunayTool = new LayerTinfour();
try (PreparedStatement st = connection.prepareStatement(
"SELECT the_geom FROM DEM")) {
try (SpatialResultSet rs = st.executeQuery().unwrap(SpatialResultSet.class)) {
while (rs.next()) {
Geometry pt = rs.getGeometry();
if(pt != null) {
delaunayTool.addVertex(pt.getCoordinate());
}
}
}
}
delaunayTool.processDelaunay();
TriangleNoiseMap.generateResultTable(connection, "RECEIVERS", "TRIANGLES",
new AtomicInteger(), delaunayTool.getVertices(), new GeometryFactory(), delaunayTool.getTriangles(),
0, 0, 1);
try(Statement st = connection.createStatement()) {
st.execute("ALTER TABLE RECEIVERS ADD COLUMN HEIGHT FLOAT");
st.execute("UPDATE RECEIVERS SET HEIGHT = ST_Z(THE_GEOM)");
}
long start = System.currentTimeMillis();
BezierContouring bezierContouring = new BezierContouring(Arrays.asList(0.,5.,10.,15.,20.,25.,30.,35.), 2154);
bezierContouring.setPointTable("RECEIVERS");
bezierContouring.setPointTableField("HEIGHT");
bezierContouring.setSmooth(false);
bezierContouring.setMergeTriangles(false);
bezierContouring.createTable(connection);
System.out.println("Contouring done in " + (System.currentTimeMillis() - start) + " ms");

assertTrue(JDBCUtilities.tableExists(connection, "CONTOURING_NOISE_MAP"));

// Check Z values in CONTOURING_NOISE_MAP
try(Statement st = connection.createStatement()) {
try(ResultSet rs = st.executeQuery("SELECT MAX(ST_ZMAX(THE_GEOM)) MAXZ, MIN(ST_ZMIN(THE_GEOM)) MINZ FROM CONTOURING_NOISE_MAP")) {
assertTrue(rs.next());
assertEquals(33.2, rs.getDouble("MAXZ"), 0.01);
assertEquals(-1.79, rs.getDouble("MINZ"), 0.01);
}
}

}
}

0 comments on commit 2159475

Please sign in to comment.