diff --git a/README.md b/README.md index 2c64bfde..02c45d17 100644 --- a/README.md +++ b/README.md @@ -140,6 +140,12 @@ If --username and/or --dbname are not specified the current username is used as --keep_projection (Default: false) Keep projection of input data + --sortby (Default: AREA) Sort features by AREA/VOLUME + AREA is the default and faster (uses ST_Area). + VOLUME is slower (uses 3D bounding-box volume) and is useful when + geometries have relatively large Z components (e.g. infrastructure / + vertical walls). + --subdivision (Default: QUADTREE) Subdivision schema QUADTREE/OCTREE --help Display this help screen. @@ -220,10 +226,11 @@ For styling see [styling 3D Tiles](styling.md) Input geometries must be of type LineString/MultilineString/Polygon/MultiPolygon/PolyhedralSurface (with z values). When the geometry is not triangulated, pg2b3dm will perform triangulation. Geometries with interior rings are supported. -For large datasets create a spatial index on the geometry column: +For large datasets create the following indexes: ``` -psql> CREATE INDEX ON the_table USING gist(st_centroid(st_envelope(geom_triangle))); +psql> CREATE INDEX ON the_table USING gist (st_centroid(st_envelope(geom))); +psql> CREATE INDEX ON the_table using btree(md5(st_asbinary(geom)::text)); ``` When there the spatial index is not present the following warning is shown. @@ -264,7 +271,11 @@ When the input geometries are distributed in a flat area (like buildings in a ci OCTREE is used when the input geometries are distributed in a cube-like area. -Most features are supported when using OCTREE subdivision, except LOD support; +Most features are supported when using OCTREE subdivision, except + +- LOD support; + +- Update boundingboxes when using explicit tiling; ## Query parameter diff --git a/md5_queries.md b/md5_queries.md new file mode 100644 index 00000000..826f860d --- /dev/null +++ b/md5_queries.md @@ -0,0 +1,134 @@ +## Queries for MD5 + + +### Initial + +1] Get bounding box whole table (1.9 s) + +```sql +SELECT st_xmin(geom1),st_ymin(geom1), st_xmax(geom1), st_ymax(geom1), st_zmin(geom1), st_zmax(geom1) FROM (select st_transform(ST_3DExtent(geom), 4979) as geom1 from bertt.nantes_reconstructed_buildings ) as t +``` +Result: + +``` +-1.8471041030488762 47.14626298148698 -1.1473131952502678 47.62268076404559 34.427586472817715 475.03764899302183 +``` + +## Tile 0_0_0.glb + +2] Count geometries in bounding box (0.2s) + +```sql +SELECT COUNT(geom) FROM bertt.nantes_reconstructed_buildings WHERE ST_Centroid(ST_Envelope(geom)) && st_transform(ST_MakeEnvelope(-1.847105103048876, 47.14626198148698, -1.1473121952502678, 47.62268176404559, 4326), 5698) +``` + +Result: 385856 + +3] Get geometries for tile 0_0_0.glb - 1000 largest geometries in whole table (2 s) + +```sql +SELECT ST_AsBinary(st_transform(geom, 4978)), id , MD5(ST_AsBinary(geom)::text) as geom_hash FROM bertt.nantes_reconstructed_buildings where ST_Centroid(ST_Envelope(geom)) && st_transform(ST_MakeEnvelope(-1.847105103048876, 47.14626198148698, -1.1473121952502678, 47.62268176404559, 4326), 5698) ORDER BY ST_Area(ST_Envelope(geom)) DESC LIMIT 1000 +``` + +md5 hashes (for example '9759cdee666f512a0c13df8245b667f9') are remembered to be excluded in higher level (z) tile + +potential improvement: make exception for first tile on z=0 - do not filter on envelope (all features are included) + +## Tile 1_0_0.glb (level 1, x=0, y=0) + +4] Filter the hashes from previous list to only geometries within tile 1_0_0 + +```sql + SELECT MD5(ST_AsBinary(geom)::text) as geom_hash + FROM bertt.nantes_reconstructed_buildings + WHERE MD5(ST_AsBinary(geom)::text) = ANY($1) + AND ST_Within( + ST_Centroid(ST_Envelope(geom)), + ST_Transform(ST_MakeEnvelope($2, $3, $4, $5, 4326), 5698) + ) +``` + +Note: Using parameterized query with array parameter instead of string concatenation. + +5] Count geometries in bounding box on level 1 excluding the geometries from tile 0_0_0.glb, including only the geometries within the tile + +```sql +SELECT COUNT(geom) FROM bertt.nantes_reconstructed_buildings WHERE ST_Centroid(ST_Envelope(geom)) && st_transform(ST_MakeEnvelope(-1.847105103048876, 47.14626198148698, -1.497208649149572, 47.384471872766284, 4326), 5698) AND MD5(ST_AsBinary(geom)::text) != ALL($1) +``` + +Note: Using parameterized query with array parameter instead of string concatenation. + +Result: 235787 + +6] Get geometries for tile 1_0_0.glb - 1000 largest geometries in tile 1_0_0 + +```sql +SELECT ST_AsBinary(st_transform(geom, 4978)), id , MD5(ST_AsBinary(geom)::text) as geom_hash FROM bertt.nantes_reconstructed_buildings where ST_Centroid(ST_Envelope(geom)) && st_transform(ST_MakeEnvelope(-1.847105103048876, 47.14626198148698, -1.497208649149572, 47.384471872766284, 4326), 5698) AND MD5(ST_AsBinary(geom)::text) != ALL($1) ORDER BY ST_Area(ST_Envelope(geom)) DESC LIMIT 1000 +``` + +Note: Using parameterized query with array parameter instead of string concatenation. + +## Issue + +List of hashes can get long (maximum z*1000 items). Previously this was handled with string concatenation which could lead to performance issues and potential SQL injection vulnerabilities. + +**Solution**: Now using parameterized queries with PostgreSQL's `= ANY()` and `!= ALL()` operators for better performance and security. + +## Spatial indexing + + Recommended Indexes + + 1. Spatial Index with MD5 Hash (Composite) + + CREATE INDEX idx_geom_centroid_hash ON the_table + USING btree(MD5(ST_AsBinary(geom_triangle)::text)); + + 2. Spatial Index (GIST) - Still Required + + CREATE INDEX idx_geom_centroid_spatial ON the_table + USING gist(ST_Centroid(ST_Envelope(geom_triangle))); + + Rationale + + The queries now use three main patterns: + + 1] Spatial filtering with MD5 hash exclusion (GetGeometrySubset): WHERE ST_Centroid(ST_Envelope(geom_triangle)) && + AND MD5(ST_AsBinary(geom_triangle)::text) != ALL($1) + + 2] MD5 hash filtering with spatial validation (FilterHashesByEnvelope): WHERE MD5(ST_AsBinary(geom_triangle)::text) = ANY($1) + AND ST_Within(ST_Centroid(ST_Envelope(geom_triangle)), ) + + 3] Hash-only filtering (GetGeometriesBoundingBox): WHERE MD5(ST_AsBinary(geom_triangle)::text) = ANY($1) + + Performance Notes: + + 1] The GIST spatial index handles the ST_Centroid(ST_Envelope(geom_triangle)) predicates + + 2] The MD5 hash BTREE index handles the MD5(ST_AsBinary(geom_triangle)::text) = ANY/!= ALL predicates + + 3] PostgreSQL will use both indexes (bitmap index scan) for queries with both predicates + + 4] Using parameterized queries with ANY/ALL operators provides better performance than string-concatenated IN/NOT IN clauses + + Optional: Materialized Hash Column + +## Solution + +The hash filtering now uses PostgreSQL's `= ANY(@param)` operator with array parameters instead of string concatenation: + +1. **Hash Inclusion (IN clause)**: Changed from `MD5(...) IN ('hash1', 'hash2', ...)` to `MD5(...) = ANY(@hashes)` with parameterized array +2. **Hash Exclusion (NOT IN clause)**: Changed from `MD5(...) NOT IN ('hash1', 'hash2', ...)` to `MD5(...) != ALL(@excludeHashes)` with parameterized array + +Benefits: +- Eliminates SQL injection risk (even though MD5 hashes are predictable) +- Better performance with large hash lists +- Cleaner, more maintainable code +- Proper use of parameterized queries + +## Todo + +- ~~idea: make a temporary blacklist table with the to be exluded hashes?~~ (Solved using parameterized arrays) + +- idea: force use of id column (longs)? + +- ~~Other solutions?~~ (Implemented using `ANY` and `ALL` operators) diff --git a/src/README.md b/src/README.md new file mode 100644 index 00000000..e69de29b diff --git a/src/b3dm.tileset.tests/CesiumTilerTests.cs b/src/b3dm.tileset.tests/CesiumTilerTests.cs index 72a80782..c9cfd923 100644 --- a/src/b3dm.tileset.tests/CesiumTilerTests.cs +++ b/src/b3dm.tileset.tests/CesiumTilerTests.cs @@ -1,7 +1,6 @@ using System; using System.Collections.Generic; using System.IO; -using B3dm.Tileset.settings; using NUnit.Framework; using pg2b3dm; using subtree; diff --git a/src/b3dm.tileset.tests/GeometryRepositoryOrderByTests.cs b/src/b3dm.tileset.tests/GeometryRepositoryOrderByTests.cs new file mode 100644 index 00000000..1f2b973e --- /dev/null +++ b/src/b3dm.tileset.tests/GeometryRepositoryOrderByTests.cs @@ -0,0 +1,29 @@ +using NUnit.Framework; + +namespace B3dm.Tileset.Tests; + +public class GeometryRepositoryOrderByTests +{ + [Test] + public void GetOrderBy_Area_UsesEnvelopeArea() + { + var sql = GeometryRepository.GetOrderBy("geom", SortBy.AREA); + Assert.That(sql, Does.Contain("ST_Area(ST_Envelope(geom))")); + Assert.That(sql, Does.Contain("DESC")); + } + + [Test] + public void GetOrderBy_Volume_UsesZExtents() + { + var sql = GeometryRepository.GetOrderBy("geom", SortBy.VOLUME); + Assert.That(sql, Does.Contain("ST_ZMax(geom)")); + Assert.That(sql, Does.Contain("ST_ZMin(geom)")); + } + + [Test] + public void TilingSettings_DefaultSortBy_IsArea() + { + var settings = new TilingSettings(); + Assert.That(settings.SortBy, Is.EqualTo(SortBy.AREA)); + } +} diff --git a/src/b3dm.tileset.tests/GeometryRepositoryWhereTests.cs b/src/b3dm.tileset.tests/GeometryRepositoryWhereTests.cs new file mode 100644 index 00000000..0213ace7 --- /dev/null +++ b/src/b3dm.tileset.tests/GeometryRepositoryWhereTests.cs @@ -0,0 +1,61 @@ +using NUnit.Framework; +using Wkx; + +namespace B3dm.Tileset.Tests; + +public class GeometryRepositoryWhereTests +{ + [Test] + public void GetWhere_2D_UseSourceEpsgDirectly() + { + var from = new Point(100000.0, 400000.0); + var to = new Point(200000.0, 500000.0); + + var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698); + + Assert.That(sql, Does.Contain("5698")); + Assert.That(sql, Does.Not.Contain("4326")); + Assert.That(sql, Does.Not.Contain("ST_Transform")); + Assert.That(sql, Does.Contain("ST_MakeEnvelope")); + } + + [Test] + public void GetWhere_2D_Epsg4326_UseSourceEpsgDirectly() + { + var from = new Point(-75.8, 38.4); + var to = new Point(-75.0, 39.8); + + var sql = GeometryRepository.GetWhere("geom", from, to, "", 4326); + + Assert.That(sql, Does.Contain("4326")); + Assert.That(sql, Does.Not.Contain("ST_Transform")); + Assert.That(sql, Does.Contain("ST_MakeEnvelope")); + } + + [Test] + public void GetWhere_3D_UseSourceEpsgDirectly() + { + var from = new Point(100000.0, 400000.0, 200.0); + var to = new Point(200000.0, 500000.0, 300.0); + + var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698); + + Assert.That(sql, Does.Contain("5698")); + Assert.That(sql, Does.Not.Contain("4979")); + Assert.That(sql, Does.Not.Contain("4326")); + Assert.That(sql, Does.Not.Contain("ST_Transform")); + Assert.That(sql, Does.Contain("ST_3DMakeBox")); + } + + [Test] + public void GetWhere_3D_ContainsCentroidAndSrid() + { + var from = new Point(100000.0, 400000.0, 200.0); + var to = new Point(200000.0, 500000.0, 300.0); + + var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698); + + Assert.That(sql, Does.Contain("st_setsrid")); + Assert.That(sql, Does.Contain("ST_3DIntersects")); + } +} diff --git a/src/b3dm.tileset/BoundingBoxRepository.cs b/src/b3dm.tileset/BoundingBoxRepository.cs index e64b7e59..7a719bbc 100644 --- a/src/b3dm.tileset/BoundingBoxRepository.cs +++ b/src/b3dm.tileset/BoundingBoxRepository.cs @@ -1,4 +1,5 @@ -using System.Data; +using System; +using System.Data; using Wkx; namespace B3dm.Tileset; @@ -16,6 +17,31 @@ public static (BoundingBox bbox, double zmin, double zmax) GetBoundingBoxForTabl return bbox3d; } + public static (BoundingBox bbox, double zmin, double zmax) GetBoundingBoxAs4979(IDbConnection conn,(BoundingBox bbox, double zmin, double zmax) bboxTable, int sourceEpsg) + { + if (sourceEpsg == 4979 || sourceEpsg == 4326) { + return bboxTable; + } + var bbox = bboxTable.bbox; + var sqlBounds = FormattableString.Invariant($@"SELECT st_xmin(geom1),st_ymin(geom1), st_xmax(geom1), st_ymax(geom1), st_zmin(geom1), st_zmax(geom1) +FROM ( + SELECT ST_3DExtent( + ST_Transform( + ST_SetSRID( + ST_3DMakeBox( + ST_MakePoint({bbox.XMin}, {bbox.YMin}, {bboxTable.zmin}), + ST_MakePoint({bbox.XMax}, {bbox.YMax}, {bboxTable.zmax}) + ), + {sourceEpsg} + ), + 4979 + ) + ) AS geom1 +) AS t"); + return GetBounds(conn, sqlBounds); + } + + private static (BoundingBox, double, double) GetBounds(IDbConnection conn, string sql) { conn.Open(); diff --git a/src/b3dm.tileset/CesiumTiler.cs b/src/b3dm.tileset/CesiumTiler.cs index 53f1c232..d5a5d964 100644 --- a/src/b3dm.tileset/CesiumTiler.cs +++ b/src/b3dm.tileset/CesiumTiler.cs @@ -3,7 +3,6 @@ using System.IO; using System.Linq; using B3dm.Tileset.Extensions; -using B3dm.Tileset.settings; using Newtonsoft.Json; using subtree; using Wkx; diff --git a/src/b3dm.tileset/FeatureCountRepository.cs b/src/b3dm.tileset/FeatureCountRepository.cs index 389a92d7..c0f2e203 100644 --- a/src/b3dm.tileset/FeatureCountRepository.cs +++ b/src/b3dm.tileset/FeatureCountRepository.cs @@ -1,23 +1,36 @@ -using Npgsql; +using System.Collections.Generic; +using System.Linq; +using Npgsql; using Wkx; namespace B3dm.Tileset; public static class FeatureCountRepository { - public static int CountFeaturesInBox(NpgsqlConnection conn, string geometry_table, string geometry_column, Point from, Point to, string query, int source_epsg, bool keepProjection = false) + public static int CountFeaturesInBox(NpgsqlConnection conn, string geometry_table, string geometry_column, Point from, Point to, string query, int source_epsg, HashSet excludeHashes = null) { var select = $"COUNT({geometry_column})"; - var where = GeometryRepository.GetWhere(geometry_column, from, to, query, source_epsg, keepProjection); + var where = GeometryRepository.GetWhere(geometry_column, from, to, query, source_epsg); + + // Add hash exclusion filter using parameterized query + if (excludeHashes != null && excludeHashes.Count > 0) { + where += $" AND MD5(ST_AsBinary({geometry_column})::text) != ALL(@excludeHashes)"; + } var sql = $"SELECT {select} FROM {geometry_table} WHERE {where}"; conn.Open(); - var cmd = new NpgsqlCommand(sql, conn); - var reader = cmd.ExecuteReader(); - reader.Read(); - var count = reader.GetInt32(0); - reader.Close(); - conn.Close(); - return count; + try { + using var cmd = new NpgsqlCommand(sql, conn); + if (excludeHashes != null && excludeHashes.Count > 0) { + cmd.Parameters.AddWithValue("excludeHashes", excludeHashes.ToArray()); + } + using var reader = cmd.ExecuteReader(); + reader.Read(); + var count = reader.GetInt32(0); + return count; + } + finally { + conn.Close(); + } } } diff --git a/src/b3dm.tileset/GeometryRepository.cs b/src/b3dm.tileset/GeometryRepository.cs index f9297b46..31029240 100644 --- a/src/b3dm.tileset/GeometryRepository.cs +++ b/src/b3dm.tileset/GeometryRepository.cs @@ -13,47 +13,69 @@ namespace B3dm.Tileset; public static class GeometryRepository { + /// /// Returns double array with 6 bounding box coordinates, xmin, ymin, xmax, ymax, zmin, zmax /// - public static double[] GetGeometriesBoundingBox(NpgsqlConnection conn, string geometry_table, string geometry_column, int epsg, Tile t, string query = "", bool keepProjection = false) + public static double[] GetGeometriesBoundingBox(NpgsqlConnection conn, string geometry_table, string geometry_column, int epsg, Tile t, HashSet tileHashes, string query = "", bool keepProjection = false) { var sqlSelect = keepProjection? $"select st_Asbinary(st_3dextent({geometry_column})) ": $"select st_Asbinary(st_3dextent(st_transform({geometry_column}, 4979))) "; - var sqlWhere = GetWhere(geometry_column, new Point(t.BoundingBox[0], t.BoundingBox[1]), new Point(t.BoundingBox[2], t.BoundingBox[3]), query, epsg, keepProjection); + + var sqlWhere = $" MD5(ST_AsBinary({geometry_column})::text) = ANY(@hashes)"; var sql = $"{sqlSelect} from {geometry_table} where {sqlWhere}"; conn.Open(); - var cmd = new NpgsqlCommand(sql, conn); - var reader = cmd.ExecuteReader(); - reader.Read(); - var stream = reader.GetStream(0); - var geometry = Geometry.Deserialize(stream); - var result = BBox3D.GetBoundingBoxPoints(geometry); - - reader.Close(); - conn.Close(); + try { + using var cmd = new NpgsqlCommand(sql, conn); + cmd.Parameters.AddWithValue("hashes", tileHashes.ToArray()); + using var reader = cmd.ExecuteReader(); + reader.Read(); + var stream = reader.GetStream(0); + var geometry = Geometry.Deserialize(stream); + var result = BBox3D.GetBoundingBoxPoints(geometry); - return result; + return result; + } + finally { + conn.Close(); + } } - public static List GetGeometrySubset(NpgsqlConnection conn, string geometry_table, string geometry_column, double[] bbox, int source_epsg, int target_srs, string shaderColumn = "", string attributesColumns = "", string query = "", string radiusColumn = "", bool keepProjection = false) + public static List GetGeometrySubset(NpgsqlConnection conn, string geometry_table, string geometry_column, double[] bbox, int source_epsg, int target_srs, string shaderColumn = "", string attributesColumns = "", string query = "", string radiusColumn = "", HashSet excludeHashes = null, int? maxFeatures = null, SortBy sortBy = SortBy.AREA) { var sqlselect = GetSqlSelect(geometry_column, shaderColumn, attributesColumns, radiusColumn, target_srs); var sqlFrom = "FROM " + geometry_table; - - // todo: fix unit test when there is no z var points = GetPoints(bbox); - var sqlWhere = GetWhere(geometry_column, points.fromPoint, points.toPoint, query, source_epsg, keepProjection); - var sql = sqlselect + sqlFrom + " where " + sqlWhere; + var sqlWhere = GetWhere(geometry_column, points.fromPoint, points.toPoint, query, source_epsg); + + // Add hash exclusion filter using parameterized query + if (excludeHashes != null && excludeHashes.Count > 0) { + sqlWhere += $" AND MD5(ST_AsBinary({geometry_column})::text) != ALL(@excludeHashes)"; + } + + var sqlOrderBy = GetOrderBy(geometry_column, sortBy); + var sqlLimit = maxFeatures.HasValue ? $" LIMIT {maxFeatures.Value}" : ""; + var sql = sqlselect + sqlFrom + " where " + sqlWhere + sqlOrderBy + sqlLimit; - var geometries = GetGeometries(conn, shaderColumn, attributesColumns, sql, radiusColumn); - return geometries; + conn.Open(); + try { + using var cmd = new NpgsqlCommand(sql, conn); + if (excludeHashes != null && excludeHashes.Count > 0) { + cmd.Parameters.AddWithValue("excludeHashes", excludeHashes.ToArray()); + } + + var geometries = GetGeometries(cmd, shaderColumn, attributesColumns, radiusColumn, geometry_column); + return geometries; + } + finally { + conn.Close(); + } } - public static string GetWhere(string geometry_column, Point from, Point to, string query, int source_epsg, bool keepProjection) + public static string GetWhere(string geometry_column, Point from, Point to, string query, int source_epsg) { var fromX = from.X.Value.ToString(CultureInfo.InvariantCulture); var fromY = from.Y.Value.ToString(CultureInfo.InvariantCulture); @@ -64,22 +86,14 @@ public static string GetWhere(string geometry_column, Point from, Point to, stri var where = ""; if (!hasZ) { - where = keepProjection ? - $"ST_Centroid(ST_Envelope({geometry_column})) && ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, {source_epsg}) {query}" : - $"ST_Centroid(ST_Envelope({geometry_column})) && st_transform(ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, 4326), {source_epsg}) {query}"; + where = $"ST_Centroid(ST_Envelope({geometry_column})) && ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, {source_epsg}) {query}"; } else { - var fromBox = keepProjection ? - $"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})" : - $"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), 4979)"; - var toBox = keepProjection ? - $"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})" : - $"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), 4979)"; + var fromBox = $"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})"; + var toBox = $"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})"; var geom = $"st_setsrid(st_makepoint((st_xmin({geometry_column}) + st_xmax({geometry_column}))/2,(st_ymin({geometry_column}) + st_ymax({geometry_column}))/2, (st_zmin({geometry_column}) + st_zmax({geometry_column}))/2), {source_epsg})"; - where = keepProjection ? - $"ST_3DIntersects({geom}, ST_3DMakeBox({fromBox}, {toBox})) {query}" : - $"ST_3DIntersects({geom}, st_transform(ST_3DMakeBox({fromBox}, {toBox}), {source_epsg})) {query}"; + where = $"ST_3DIntersects({geom}, ST_3DMakeBox({fromBox}, {toBox})) {query}"; } return where; @@ -98,6 +112,8 @@ public static string GetSqlSelect(string geometry_column, string shaderColumn, s if (radiusColumn != String.Empty) { sqlselect = $"{sqlselect}, {radiusColumn} "; } + // Add MD5 hash of geometry + sqlselect = $"{sqlselect}, MD5(ST_AsBinary({geometry_column})::text) as geom_hash "; return sqlselect; } @@ -107,15 +123,22 @@ public static string GetGeometryColumn(string geometry_column, int target_srs) return $"st_transform({geometry_column}, {target_srs})"; } - public static List GetGeometries(NpgsqlConnection conn, string shaderColumn, string attributesColumns, string sql, string radiusColumn) + public static string GetOrderBy(string geometry_column, SortBy sortBy) + { + if (sortBy == SortBy.VOLUME) { + return $" ORDER BY (ST_XMax({geometry_column}) - ST_XMin({geometry_column})) *(ST_YMax({geometry_column}) - ST_YMin({geometry_column})) *(ST_ZMax({geometry_column}) - ST_ZMin({geometry_column})) DESC"; + } + return $" ORDER BY ST_Area(ST_Envelope({geometry_column})) DESC"; + } + + public static List GetGeometries(NpgsqlCommand cmd, string shaderColumn, string attributesColumns, string radiusColumn, string geometry_column = "") { var geometries = new List(); - conn.Open(); - var cmd = new NpgsqlCommand(sql, conn); var reader = cmd.ExecuteReader(); var attributesColumnIds = new Dictionary(); var shadersColumnId = int.MinValue; var radiusColumnId = int.MinValue; + var hashColumnId = int.MinValue; if (attributesColumns != String.Empty) { var attributesColumnsList = attributesColumns.Split(',').ToList(); @@ -133,6 +156,11 @@ public static List GetGeometries(NpgsqlConnection conn, string s radiusColumnId = FindField(reader, radiusColumn).Value; } } + // Find hash column + var hashFld = FindField(reader, "geom_hash"); + if (hashFld.HasValue) { + hashColumnId = hashFld.Value; + } var batchId = 0; while (reader.Read()) { @@ -154,13 +182,15 @@ public static List GetGeometries(NpgsqlConnection conn, string s var radius = reader.GetFieldValue(radiusColumnId); geometryRecord.Radius = Convert.ToSingle(radius); } + if (hashColumnId != int.MinValue) { + geometryRecord.Hash = reader.GetString(hashColumnId); + } geometries.Add(geometryRecord); batchId++; } reader.Close(); - conn.Close(); return geometries; } diff --git a/src/b3dm.tileset/OctreeTiler.cs b/src/b3dm.tileset/OctreeTiler.cs index ac3dd9a8..c42ed92e 100644 --- a/src/b3dm.tileset/OctreeTiler.cs +++ b/src/b3dm.tileset/OctreeTiler.cs @@ -1,9 +1,7 @@ -using System; -using System.Collections.Generic; +using System.Collections.Generic; using System.IO; -using B3dm.Tileset.settings; +using System.Linq; using Npgsql; -using pg2b3dm; using subtree; using Wkx; @@ -28,14 +26,23 @@ public OctreeTiler(string connectionString, InputTable inputTable, TilingSetting public List GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile, List tiles) { - return GenerateTiles3D(bbox, level, tile, tiles, null); + return GenerateTiles3D(bbox, level, tile, tiles, null, null); } public List GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile, List tiles, Dictionary tileBounds) { + return GenerateTiles3D(bbox, level, tile, tiles, tileBounds, null); + } + + public List GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile, List tiles, Dictionary tileBounds, HashSet processedGeometries) + { + if (processedGeometries == null) { + processedGeometries = new HashSet(); + } + var where = inputTable.GetQueryClause(); - var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin, bbox.ZMin), new Point(bbox.XMax, bbox.YMax, bbox.ZMax), where, inputTable.EPSGCode, tilingSettings.KeepProjection); + var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin, bbox.ZMin), new Point(bbox.XMax, bbox.YMax, bbox.ZMax), where, inputTable.EPSGCode, processedGeometries); if (numberOfFeatures == 0) { var t2 = new Tile3D(level, tile.X, tile.Y, tile.Z); t2.Available = false; @@ -46,6 +53,9 @@ public List GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile, } } else if (numberOfFeatures > tilingSettings.MaxFeaturesPerTile) { + // First, create a tile with the largest geometries up to MaxFeaturesPerTile for this level + var localProcessedGeometries = CreateTileForLargestGeometries3D(bbox, level, tile, tiles, tileBounds, where, processedGeometries); + level++; for (var x = 0; x < 2; x++) { for (var y = 0; y < 2; y++) { @@ -64,46 +74,77 @@ public List GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile, var zend = z_start + dz; var bbox3d = new BoundingBox3D(xstart, ystart, z_start, xend, yend, zend); + var bboxOctant = new BoundingBox(xstart, ystart, xend, yend); + var new_tile = new Tile3D(level, tile.X * 2 + x, tile.Y * 2 + y, tile.Z * 2 + z); - GenerateTiles3D(bbox3d, level, new_tile, tiles, tileBounds); + GenerateTiles3D(bbox3d, level, new_tile, tiles, tileBounds, localProcessedGeometries); } } } } else { - var boundingBox = new BoundingBox(bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax); + CreateTile3D(bbox, level, tile, tiles, tileBounds, where, processedGeometries); + } - int target_srs = 4978; + return tiles; - if (tilingSettings.KeepProjection) { - target_srs = inputTable.EPSGCode; - } + } - var bbox1 = new double[] { bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax, bbox.ZMin, bbox.ZMax }; - var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, tilingSettings.KeepProjection); + private HashSet CreateTileForLargestGeometries3D(BoundingBox3D bbox, int level, Tile3D tile, List tiles, Dictionary tileBounds, string where, HashSet processedGeometries) + { + // clone processedGeometries to avoid modifying the original set in recursive calls + var localProcessedGeometries = new HashSet(processedGeometries); - if (geometries.Count > 0) { + // Get the largest geometries (up to MaxFeaturesPerTile) for this tile at this level + int target_srs = tilingSettings.KeepProjection ? inputTable.EPSGCode : 4978; - if (!tilingSettings.SkipCreateTiles) { - var bytes = TileWriter.ToTile(geometries, tilesetSettings.Translation, copyright: tilesetSettings.Copyright, addOutlines: stylingSettings.AddOutlines, defaultColor: stylingSettings.DefaultColor, defaultMetallicRoughness: stylingSettings.DefaultMetallicRoughness, doubleSided: stylingSettings.DoubleSided, defaultAlphaMode: stylingSettings.DefaultAlphaMode, alphaCutoff: stylingSettings.AlphaCutoff, createGltf: tilingSettings.CreateGltf); - var file = $"{tilesetSettings.OutputSettings.ContentFolder}{Path.AltDirectorySeparatorChar}{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}.glb"; - Console.Write($"\rCreating tile: {file} "); - File.WriteAllBytes($"{file}", bytes); - } - tile.Available = true; - } - else { - tile.Available = false; - } + var bbox1 = new double[] { bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax, bbox.ZMin, bbox.ZMax }; + var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy); - tiles.Add(tile); - if (tileBounds != null) { - var key = $"{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}"; - tileBounds[key] = bbox; + if (geometriesToProcess.Count > 0) { + foreach (var geom in geometriesToProcess.Where(geom => !string.IsNullOrEmpty(geom.Hash))) { + localProcessedGeometries.Add(geom.Hash); } + + var file = $"{tilesetSettings.OutputSettings.ContentFolder}{Path.AltDirectorySeparatorChar}{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}.glb"; + TileCreationHelper.WriteTileIfNeeded(geometriesToProcess, tilesetSettings.Translation, stylingSettings, tilesetSettings.Copyright, tilingSettings.CreateGltf, tilingSettings.SkipCreateTiles, file, file); + tile.Available = true; } - return tiles; + tiles.Add(tile); + if (tileBounds != null) { + var key = $"{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}"; + tileBounds[key] = bbox; + } + + return localProcessedGeometries; + } + + private void CreateTile3D(BoundingBox3D bbox, int level, Tile3D tile, List tiles, Dictionary tileBounds, string where, HashSet processedGeometries) + { + int target_srs = tilingSettings.KeepProjection ? inputTable.EPSGCode : 4978; + var bbox1 = new double[] { bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax, bbox.ZMin, bbox.ZMax }; + var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy); + + if (geometries.Count > 0) { + // Collect hashes of processed geometries + foreach (var geom in geometries.Where(geom => !string.IsNullOrEmpty(geom.Hash))) { + processedGeometries.Add(geom.Hash); + } + + var file = $"{tilesetSettings.OutputSettings.ContentFolder}{Path.AltDirectorySeparatorChar}{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}.glb"; + TileCreationHelper.WriteTileIfNeeded(geometries, tilesetSettings.Translation, stylingSettings, tilesetSettings.Copyright, tilingSettings.CreateGltf, tilingSettings.SkipCreateTiles, file, file); + tile.Available = true; + } + else { + tile.Available = false; + } + + tiles.Add(tile); + if (tileBounds != null) { + var key = $"{tile.Level}_{tile.Z}_{tile.X}_{tile.Y}"; + tileBounds[key] = bbox; + } } } diff --git a/src/b3dm.tileset/OutputDirectoryCreator.cs b/src/b3dm.tileset/OutputDirectoryCreator.cs index 8f1c8f69..998b9a10 100644 --- a/src/b3dm.tileset/OutputDirectoryCreator.cs +++ b/src/b3dm.tileset/OutputDirectoryCreator.cs @@ -1,5 +1,5 @@ using System.IO; -using B3dm.Tileset.settings; +using B3dm.Tileset; namespace pg2b3dm; diff --git a/src/b3dm.tileset/QuadtreeTiler.cs b/src/b3dm.tileset/QuadtreeTiler.cs index 637f0c95..d1ae083b 100644 --- a/src/b3dm.tileset/QuadtreeTiler.cs +++ b/src/b3dm.tileset/QuadtreeTiler.cs @@ -4,7 +4,6 @@ using System.Linq; using B3dm.Tileset; using B3dm.Tileset.Extensions; -using B3dm.Tileset.settings; using Npgsql; using subtree; using Wkx; @@ -23,8 +22,10 @@ public class QuadtreeTiler private readonly bool skipCreateTiles; private readonly StylingSettings stylingSettings; private InputTable inputTable; + private readonly bool useImplicitTiling; + private readonly SortBy sortBy; - public QuadtreeTiler(string connectionString, InputTable inputTable, StylingSettings stylingSettings, int maxFeaturesPerTile, double[] translation, string outputFolder, List lods, string copyright = "", bool skipCreateTiles = false) + public QuadtreeTiler(string connectionString, InputTable inputTable, StylingSettings stylingSettings, int maxFeaturesPerTile, double[] translation, string outputFolder, List lods, string copyright = "", bool skipCreateTiles = false, bool useImplicitTiling = false, SortBy sortBy = SortBy.AREA) { this.conn = new NpgsqlConnection(connectionString); this.inputTable = inputTable; @@ -36,10 +37,16 @@ public QuadtreeTiler(string connectionString, InputTable inputTable, StylingSett this.copyright = copyright; this.skipCreateTiles = skipCreateTiles; this.stylingSettings = stylingSettings; + this.useImplicitTiling = useImplicitTiling; + this.sortBy = sortBy; } - public List GenerateTiles(BoundingBox bbox, Tile tile, List tiles, int lod = 0, bool createGltf = false, bool keepProjection = false) + public List GenerateTiles(BoundingBox bbox, Tile tile, List tiles, int lod = 0, bool createGltf = false, bool keepProjection = false, HashSet processedGeometries = null) { + if (processedGeometries == null) { + processedGeometries = new HashSet(); + } + var where = inputTable.GetQueryClause(); var lodquery = LodQuery.GetLodQuery(inputTable.LodColumn, lod); @@ -48,15 +55,15 @@ public List GenerateTiles(BoundingBox bbox, Tile tile, List tiles, i where += $" and {lodquery}"; } - var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin), new Point(bbox.XMax, bbox.YMax), where, source_epsg, keepProjection); + var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin), new Point(bbox.XMax, bbox.YMax), where, source_epsg, processedGeometries); if (numberOfFeatures == 0) { tile.Available = false; tiles.Add(tile); } else if (numberOfFeatures > maxFeaturesPerTile) { - tile.Available = false; - tiles.Add(tile); + // First, get the largest geometries up to maxFeaturesPerTile for this level + var localProcessedGeometries = CreateTileForLargestGeometries(bbox, tile, tiles, where, lod, createGltf, keepProjection, processedGeometries); var z = tile.Z + 1; @@ -74,11 +81,39 @@ public List GenerateTiles(BoundingBox bbox, Tile tile, List tiles, i var bboxQuad = new BoundingBox(xstart, ystart, xend, yend); var new_tile = new Tile(z, tile.X * 2 + x, tile.Y * 2 + y); new_tile.BoundingBox = bboxQuad.ToArray(); - GenerateTiles(bboxQuad, new_tile, tiles, lod, createGltf, keepProjection); + + GenerateTiles(bboxQuad, new_tile, tiles, lod, createGltf, keepProjection, localProcessedGeometries); } } } else { + CreateTile(bbox, tile, tiles, where, lod, createGltf, keepProjection, processedGeometries); + } + + return tiles; + } + + private HashSet CreateTileForLargestGeometries(BoundingBox bbox, Tile tile, List tiles, string where, int lod, bool createGltf, bool keepProjection, HashSet processedGeometries) + { + // clone processedIds to avoid modifying the original set in recursive calls + var localProcessedGeometries = new HashSet(processedGeometries); + var tileHashes = new HashSet(); + + // Get the largest geometries (up to maxFeaturesPerTile) for this tile at this level + tile.Available = false; + tile.BoundingBox = bbox.ToArray(); + + int target_srs = keepProjection ? source_epsg : 4978; + + var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, maxFeaturesPerTile, sortBy: sortBy); + + if (geometriesToProcess.Count > 0) { + + // Collect hashes of processed geometries + foreach (var geom in geometriesToProcess.Where(geom => !string.IsNullOrEmpty(geom.Hash))) { + localProcessedGeometries.Add(geom.Hash); + tileHashes.Add(geom.Hash); + } var file = $"{tile.Z}_{tile.X}_{tile.Y}"; if (inputTable.LodColumn != String.Empty) { @@ -90,55 +125,89 @@ public List GenerateTiles(BoundingBox bbox, Tile tile, List tiles, i Console.Write($"\rCreating tile: {file} "); tile.ContentUri = file; - int target_srs = 4978; + tile.Lod = lod; - if(keepProjection) { - target_srs = source_epsg; - } + var outputPath = $"{outputFolder}{Path.AltDirectorySeparatorChar}{file}"; + TileCreationHelper.WriteTileIfNeeded(geometriesToProcess, translation, stylingSettings, copyright, createGltf, skipCreateTiles, outputPath, file); - byte[] bytes = null; + ProcessLodLevels(bbox, tile, lod, createGltf, keepProjection, localProcessedGeometries); + if (!useImplicitTiling) { + UpdateTileBoundingBox(tile, tileHashes, where, keepProjection); + } - var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, keepProjection); - // var scale = new double[] { 1, 1, 1 }; - if (geometries.Count > 0) { + tile.Available = true; + } + + tiles.Add(tile); + return localProcessedGeometries; + } - tile.Lod = lod; + private void CreateTile(BoundingBox bbox, Tile tile, List tiles, string where, int lod, bool createGltf, bool keepProjection, HashSet processedGeometries) + { + tile.BoundingBox = bbox.ToArray(); + var tileHashes = new HashSet(); - if (!skipCreateTiles) { - bytes = TileWriter.ToTile(geometries, translation, copyright: copyright, addOutlines: stylingSettings.AddOutlines, defaultColor: stylingSettings.DefaultColor, defaultMetallicRoughness: stylingSettings.DefaultMetallicRoughness, doubleSided: stylingSettings.DoubleSided, defaultAlphaMode: stylingSettings.DefaultAlphaMode, alphaCutoff: stylingSettings.AlphaCutoff, createGltf: createGltf); - File.WriteAllBytes($"{outputFolder}{Path.AltDirectorySeparatorChar}{file}", bytes); - } - if (inputTable.LodColumn != String.Empty) { - if (lod < lods.Max()) { - // take the next lod - var currentIndex = lods.FindIndex(p => p == lod); - var nextIndex = currentIndex + 1; - var nextLod = lods[nextIndex]; - // make a copy of the tile - var t2 = new Tile(tile.Z, tile.X, tile.Y); - t2.BoundingBox = tile.BoundingBox; - var lodNextTiles = GenerateTiles(bbox, t2, new List(), nextLod, createGltf, keepProjection); - tile.Children = lodNextTiles; - }; - } + var file = $"{tile.Z}_{tile.X}_{tile.Y}"; + if (inputTable.LodColumn != String.Empty) { + file += $"_{lod}"; + } - // next code is used to fix geometries that have centroid in the tile, but some parts outside... - var bbox_geometries = GeometryRepository.GetGeometriesBoundingBox(conn, inputTable.TableName, inputTable.GeometryColumn, source_epsg, tile, where, keepProjection); - var bbox_tile = new double[] { bbox_geometries[0], bbox_geometries[1], bbox_geometries[2], bbox_geometries[3] }; - tile.BoundingBox = bbox_tile; - tile.ZMin = bbox_geometries[4]; - tile.ZMax = bbox_geometries[5]; + var ext = createGltf ? ".glb" : ".b3dm"; + file += ext; + Console.Write($"\rCreating tile: {file} "); + tile.ContentUri = file; - tile.Available = true; + int target_srs = keepProjection ? source_epsg : 4978; - if (skipCreateTiles) { tile.Available = true; } + var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, sortBy: sortBy); + + if (geometries.Count > 0) { + // Collect hashes of processed geometries + foreach (var geom in geometries.Where(g => !string.IsNullOrEmpty(g.Hash))) { + tileHashes.Add(geom.Hash); + processedGeometries.Add(geom.Hash); } - else { - tile.Available = false; + + tile.Lod = lod; + + var outputPath = $"{outputFolder}{Path.AltDirectorySeparatorChar}{file}"; + TileCreationHelper.WriteTileIfNeeded(geometries, translation, stylingSettings, copyright, createGltf, skipCreateTiles, outputPath, file); + + ProcessLodLevels(bbox, tile, lod, createGltf, keepProjection, processedGeometries); + if (!useImplicitTiling && tileHashes.Count > 0) { + UpdateTileBoundingBox(tile, tileHashes, where, keepProjection); } - tiles.Add(tile); + + tile.Available = true; + } + else { + tile.Available = false; } + tiles.Add(tile); + } - return tiles; + private void ProcessLodLevels(BoundingBox bbox, Tile tile, int lod, bool createGltf, bool keepProjection, HashSet processedGeometries) + { + if (inputTable.LodColumn != String.Empty && lod < lods.Max()) { + // take the next lod + var currentIndex = lods.FindIndex(p => p == lod); + var nextIndex = currentIndex + 1; + var nextLod = lods[nextIndex]; + // make a copy of the tile + var t2 = new Tile(tile.Z, tile.X, tile.Y); + t2.BoundingBox = tile.BoundingBox; + var lodNextTiles = GenerateTiles(bbox, t2, new List(), nextLod, createGltf, keepProjection, processedGeometries); + tile.Children = lodNextTiles; + } + } + + private void UpdateTileBoundingBox(Tile tile, HashSet tileHashes, string where, bool keepProjection) + { + // next code is used to fix geometries that have centroid in the tile, but some parts outside... + var bbox_geometries = GeometryRepository.GetGeometriesBoundingBox(conn, inputTable.TableName, inputTable.GeometryColumn, source_epsg, tile, tileHashes, where, keepProjection); + var bbox_tile = new double[] { bbox_geometries[0], bbox_geometries[1], bbox_geometries[2], bbox_geometries[3] }; + tile.BoundingBox = bbox_tile; + tile.ZMin = bbox_geometries[4]; + tile.ZMax = bbox_geometries[5]; } } diff --git a/src/b3dm.tileset/SpatialIndexChecker.cs b/src/b3dm.tileset/SpatialIndexChecker.cs index c4d95ae3..eec6d978 100644 --- a/src/b3dm.tileset/SpatialIndexChecker.cs +++ b/src/b3dm.tileset/SpatialIndexChecker.cs @@ -4,9 +4,22 @@ namespace B3dm.Tileset; public static class SpatialIndexChecker { public static bool HasSpatialIndex(NpgsqlConnection conn, string geometry_table, string geometry_column) + { + var indexDef = $"%st_centroid(st_envelope({geometry_column}))%"; + + return HasIndex(conn, geometry_table, indexDef); + } + + public static bool HasMd5Index(NpgsqlConnection conn, string geometry_table, string geometry_column) + { + var indexDef = $"%md5((st_asbinary({geometry_column}))::text)%"; + return HasIndex(conn, geometry_table, indexDef); + } + + private static bool HasIndex(NpgsqlConnection conn, string geometry_table, string indexDef) { var schema = "public"; - if(geometry_table.Contains('.')) { + if (geometry_table.Contains('.')) { var items = geometry_table.Split('.', 2); schema = items[0]; geometry_table = items[1]; @@ -20,7 +33,7 @@ public static bool HasSpatialIndex(NpgsqlConnection conn, string geometry_table, cmd.Parameters.AddWithValue("schema", schema); cmd.Parameters.AddWithValue("geometry_table", geometry_table); - cmd.Parameters.AddWithValue("index", "%st_centroid(st_envelope(" + geometry_column + "%))"); + cmd.Parameters.AddWithValue("index", indexDef); var reader = cmd.ExecuteReader(); reader.Read(); diff --git a/src/b3dm.tileset/TileCreationHelper.cs b/src/b3dm.tileset/TileCreationHelper.cs new file mode 100644 index 00000000..de2c59d4 --- /dev/null +++ b/src/b3dm.tileset/TileCreationHelper.cs @@ -0,0 +1,21 @@ +using System; +using System.Collections.Generic; +using System.IO; +using pg2b3dm; +using Wkb2Gltf; + +namespace B3dm.Tileset; + +public static class TileCreationHelper +{ + public static void WriteTileIfNeeded(List geometries, double[] translation, StylingSettings stylingSettings, string copyright, bool createGltf, bool skipCreateTiles, string outputPath, string displayName) + { + if (skipCreateTiles) { + return; + } + + var bytes = TileWriter.ToTile(geometries, translation, copyright: copyright, addOutlines: stylingSettings.AddOutlines, defaultColor: stylingSettings.DefaultColor, defaultMetallicRoughness: stylingSettings.DefaultMetallicRoughness, doubleSided: stylingSettings.DoubleSided, defaultAlphaMode: stylingSettings.DefaultAlphaMode, alphaCutoff: stylingSettings.AlphaCutoff, createGltf: createGltf); + Console.Write($"\rCreating tile: {displayName} "); + File.WriteAllBytes(outputPath, bytes); + } +} diff --git a/src/b3dm.tileset/settings/InputTable.cs b/src/b3dm.tileset/settings/InputTable.cs index 98aa45e1..f69d953b 100644 --- a/src/b3dm.tileset/settings/InputTable.cs +++ b/src/b3dm.tileset/settings/InputTable.cs @@ -1,4 +1,4 @@ -namespace B3dm.Tileset.settings; +namespace B3dm.Tileset; public class InputTable { diff --git a/src/b3dm.tileset/settings/OutputSettings.cs b/src/b3dm.tileset/settings/OutputSettings.cs index 3b5bb00a..c4bb224b 100644 --- a/src/b3dm.tileset/settings/OutputSettings.cs +++ b/src/b3dm.tileset/settings/OutputSettings.cs @@ -1,4 +1,4 @@ -namespace B3dm.Tileset.settings; +namespace B3dm.Tileset; public class OutputSettings { diff --git a/src/b3dm.tileset/settings/SortBy.cs b/src/b3dm.tileset/settings/SortBy.cs new file mode 100644 index 00000000..c39f187d --- /dev/null +++ b/src/b3dm.tileset/settings/SortBy.cs @@ -0,0 +1,11 @@ +using Newtonsoft.Json; +using Newtonsoft.Json.Converters; + +namespace B3dm.Tileset; + +[JsonConverter(typeof(StringEnumConverter))] +public enum SortBy +{ + AREA, + VOLUME +} diff --git a/src/b3dm.tileset/settings/StylingSettings.cs b/src/b3dm.tileset/settings/StylingSettings.cs index 25747b0d..1d24108f 100644 --- a/src/b3dm.tileset/settings/StylingSettings.cs +++ b/src/b3dm.tileset/settings/StylingSettings.cs @@ -1,6 +1,6 @@ using SharpGLTF.Materials; -namespace B3dm.Tileset.settings; +namespace B3dm.Tileset; public class StylingSettings { diff --git a/src/b3dm.tileset/settings/TilesetSettings.cs b/src/b3dm.tileset/settings/TilesetSettings.cs index 71fb39ed..f7f4d0b8 100644 --- a/src/b3dm.tileset/settings/TilesetSettings.cs +++ b/src/b3dm.tileset/settings/TilesetSettings.cs @@ -1,6 +1,6 @@ using System; -namespace B3dm.Tileset.settings; +namespace B3dm.Tileset; public class TilesetSettings { diff --git a/src/b3dm.tileset/settings/TilingSettings.cs b/src/b3dm.tileset/settings/TilingSettings.cs index e4892981..b6e56671 100644 --- a/src/b3dm.tileset/settings/TilingSettings.cs +++ b/src/b3dm.tileset/settings/TilingSettings.cs @@ -1,6 +1,6 @@ using System.Collections.Generic; -namespace B3dm.Tileset.settings; +namespace B3dm.Tileset; public class TilingSettings { @@ -14,6 +14,8 @@ public class TilingSettings public int MaxFeaturesPerTile { get; set; } = 1000; + public SortBy SortBy { get; set; } = SortBy.AREA; + public bool UseImplicitTiling { get; set; } = true; public List Lods { get; set; } diff --git a/src/pg2b3dm.database.tests/UnitTest1.cs b/src/pg2b3dm.database.tests/UnitTest1.cs index 39314603..a7a3b434 100644 --- a/src/pg2b3dm.database.tests/UnitTest1.cs +++ b/src/pg2b3dm.database.tests/UnitTest1.cs @@ -1,6 +1,5 @@ using B3dm.Tileset; using B3dm.Tileset.Extensions; -using B3dm.Tileset.settings; using DotNet.Testcontainers.Builders; using Npgsql; using subtree; @@ -17,7 +16,7 @@ public class UnitTest1 public async Task Setup() { _containerPostgres = new PostgreSqlBuilder() - .WithImage("postgis/postgis:16-3.4-alpine") + .WithImage("postgis/postgis:18-3.6-alpine") .WithWaitStrategy(Wait.ForUnixContainer().UntilInternalTcpPortIsAvailable(5432)) .Build(); await _containerPostgres.StartAsync(); @@ -44,9 +43,12 @@ public void TestArvieuxBuildingsOctree() OutputDirectoryCreator.GetFolders("output"); var connectionString = _containerPostgres.GetConnectionString(); var conn = new NpgsqlConnection(connectionString); - var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom"); + // Get bbox in source projection (EPSG:5698) for spatial queries + var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom", true); + // Get WGS84 bbox separately for ECEF translation + var bbox_wgs84 = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom", false); - var center_wgs84 = bbox_table.bbox.GetCenter(); + var center_wgs84 = bbox_wgs84.bbox.GetCenter(); var translation = SpatialConverter.GeodeticToEcef((double)center_wgs84.X!, (double)center_wgs84.Y!, 0); var trans = new double[] { translation.X, translation.Y, }; @@ -78,7 +80,7 @@ public void TestArvieuxBuildingsOctree() var implicitTiler = new OctreeTiler(connectionString, inputTable, tilingSettings, stylingSettings, tilesetSettings); var tiles = implicitTiler.GenerateTiles3D(boundingBox3D, 0, new Tile3D(0, 0, 0, 0 ), new List()); - Assert.That(tiles.Count, Is.EqualTo(36)); + Assert.That(tiles.Count, Is.EqualTo(25)); } [Test] @@ -119,7 +121,7 @@ public void TestArvieuxBuildingsOctreeKeepProjection() var implicitTiler = new OctreeTiler(connectionString, inputTable, tilingSettings, stylingSettings, tilesetSettings); var tiles = implicitTiler.GenerateTiles3D(boundingBox3D, 0, new Tile3D(0, 0, 0, 0), new List()); - Assert.That(tiles.Count, Is.EqualTo(36)); + Assert.That(tiles.Count, Is.EqualTo(25)); } @@ -168,12 +170,13 @@ public void ImplicitTilingTest() trans, "output/content", new List() { 0 }, - skipCreateTiles: true); + skipCreateTiles: true, + useImplicitTiling: true); var tiles = implicitTiler.GenerateTiles( bbox_wgs84.bbox, new Tile(0, 0, 0), new List()); - Assert.That(tiles.Count, Is.EqualTo(29)); + Assert.That(tiles.Count, Is.EqualTo(17)); } [Test] @@ -201,12 +204,13 @@ public void LodTest() trans, "output/content", new List() { 0, 1 }, - skipCreateTiles: true); + skipCreateTiles: true, + useImplicitTiling: true); var tiles = implicitTiler.GenerateTiles( bbox_wgs84.bbox, new Tile(0, 0, 0), new List()); - Assert.That(tiles.Count, Is.EqualTo(145)); + Assert.That(tiles.Count, Is.EqualTo(89)); } @@ -231,7 +235,8 @@ public void GeometryTest() trans, "output/content", new List() { 0 }, - skipCreateTiles: false); + skipCreateTiles: false, + useImplicitTiling: true); var tile = new Tile(0, 0, 0) { BoundingBox = bbox_wgs84.bbox.ToArray() }; diff --git a/src/pg2b3dm.sln b/src/pg2b3dm.sln index 63549f8b..d1a1488a 100644 --- a/src/pg2b3dm.sln +++ b/src/pg2b3dm.sln @@ -24,7 +24,8 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "doc", "doc", "{0403ED2E-B5A ..\..\..\..\..\Users\bertt\Desktop\delaware.png = ..\..\..\..\..\Users\bertt\Desktop\delaware.png ..\..\..\..\..\Users\bert\Desktop\demo_pg2b3dm.gif = ..\..\..\..\..\Users\bert\Desktop\demo_pg2b3dm.gif ..\getting_started.md = ..\getting_started.md - C:\Users\bert\Desktop\test.csv\materials_for_features.csv = C:\Users\bert\Desktop\test.csv\materials_for_features.csv + ..\dataprocessing\materials_for_features.csv = ..\dataprocessing\materials_for_features.csv + ..\md5_queries.md = ..\md5_queries.md ..\README.md = ..\README.md release_notes_0.10.md = release_notes_0.10.md release_notes_0.11.md = release_notes_0.11.md diff --git a/src/pg2b3dm/Options.cs b/src/pg2b3dm/Options.cs index 47c279e6..64cdf313 100644 --- a/src/pg2b3dm/Options.cs +++ b/src/pg2b3dm/Options.cs @@ -83,6 +83,9 @@ public class Options [Option("keep_projection", Required = false, Default = false, HelpText = "Keep projection")] public bool? KeepProjection { get; set; } + [Option("sortby", Required = false, Default = SortBy.AREA, HelpText = "Sort features by AREA (default) or VOLUME")] + public SortBy SortBy { get; set; } + [Option("subdivision", Required=false, Default=SubdivisionScheme.QUADTREE, HelpText = "Subdivision schema QUADTREE/OCTREE" )] public SubdivisionScheme subdivisionScheme { get; set; } } diff --git a/src/pg2b3dm/Program.cs b/src/pg2b3dm/Program.cs index 81244c7e..e2cc340b 100644 --- a/src/pg2b3dm/Program.cs +++ b/src/pg2b3dm/Program.cs @@ -9,7 +9,6 @@ using subtree; using B3dm.Tileset.Extensions; using SharpGLTF.Schema2; -using B3dm.Tileset.settings; namespace pg2b3dm; @@ -104,18 +103,33 @@ static void Main(string[] args) Console.WriteLine($"Spatial index detected on {inputTable.TableName}.{inputTable.GeometryColumn}"); } + // Check md5 index + var hasMd5Index = SpatialIndexChecker.HasMd5Index(conn, inputTable.TableName, inputTable.GeometryColumn); + if (!hasMd5Index) { + Console.WriteLine(); + Console.WriteLine("-----------------------------------------------------------------------------"); + Console.WriteLine($"WARNING: No md5 index detected on {inputTable.TableName}.{inputTable.GeometryColumn}"); + Console.WriteLine("Fix: add a md5 index, for example: "); + Console.WriteLine($"'CREATE INDEX ON {inputTable.TableName} using btree(md5(st_asbinary({inputTable.GeometryColumn})::text))'"); + Console.WriteLine("-----------------------------------------------------------------------------"); + Console.WriteLine(); + } + else { + Console.WriteLine($"Md5 index detected on {inputTable.TableName}.{inputTable.GeometryColumn}"); + } + var skipCreateTiles = (bool)o.SkipCreateTiles; Console.WriteLine("Skip create tiles: " + skipCreateTiles); Console.WriteLine($"Query bounding box of {inputTable.TableName}.{inputTable.GeometryColumn}..."); var where = (inputTable.Query != string.Empty ? $" where {inputTable.Query}" : String.Empty); - var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, inputTable.TableName, inputTable.GeometryColumn, keepProjection, where); + var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, inputTable.TableName, inputTable.GeometryColumn, true, where); var bbox = bbox_table.bbox; var zmin = bbox_table.zmin; var zmax = bbox_table.zmax; - var proj = keepProjection ? $"EPSG:{source_epsg}" : $"EPSG:4326 (WGS84)"; + var proj = $"EPSG:{source_epsg}"; Console.WriteLine($"Bounding box for {inputTable.TableName}.{inputTable.GeometryColumn} ({proj}): " + $"{Math.Round(bbox.XMin, 8)}, {Math.Round(bbox.YMin, 8)}, " + $"{Math.Round(bbox.XMax, 8)}, {Math.Round(bbox.YMax, 8)}"); @@ -134,14 +148,23 @@ static void Main(string[] args) Console.WriteLine($"Attribute columns: {att}"); var center = bbox.GetCenter(); - Console.WriteLine($"Center ({proj}): {center.X}, {center.Y}"); + var center_z = (zmin + zmax) / 2; + Console.WriteLine($"Center ({proj}): {center.X}, {center.Y}, {center_z}"); Tiles3DExtensions.RegisterExtensions(); - var translation = keepProjection ? - [(double)center.X, (double)center.Y, 0] : - Translation.ToEcef(center); - Console.WriteLine($"Translation: {String.Join(',', translation)}"); + double[] translation = [(double)center.X, (double)center.Y, center_z]; + + var bbox_wgs84 = bbox; // fallback, only set when !keepProjection + if(!keepProjection) { + // Convert bbox to EPSG:4979 for ECEF translation and bounding volume region + var bbox_4979 = BoundingBoxRepository.GetBoundingBoxAs4979(conn, bbox_table, source_epsg); + bbox_wgs84 = bbox_4979.bbox; + var center_wgs84 = bbox_wgs84.GetCenter(); + var p = new Wkx.Point((double)center_wgs84.X, (double)center_wgs84.Y, center_z); + translation = Translation.ToEcef(p); + } + Console.WriteLine($"Translation (EPSG:4978): {String.Join(',', translation)}"); var lodcolumn = o.LodColumn; var addOutlines = (bool)o.AddOutlines; @@ -160,6 +183,7 @@ static void Main(string[] args) Console.WriteLine($"Refinement: {refinement}"); Console.WriteLine($"Keep projection: {keepProjection}"); Console.WriteLine($"Subdivision scheme: {subdivisionScheme}"); + Console.WriteLine($"Sort by: {o.SortBy}"); if (keepProjection && !useImplicitTiling) { Console.WriteLine("Warning: keepProjection is only supported with implicit tiling."); @@ -182,7 +206,7 @@ static void Main(string[] args) var rootBoundingVolumeRegion = keepProjection ? bbox.ToRegion(zmin, zmax) : - bbox.ToRadians().ToRegion(zmin, zmax); + bbox_wgs84.ToRadians().ToRegion(zmin, zmax); Console.WriteLine($"Maximum features per tile: " + maxFeaturesPerTile); @@ -222,6 +246,7 @@ static void Main(string[] args) tilingSettings.UseImplicitTiling = useImplicitTiling; tilingSettings.SkipCreateTiles = skipCreateTiles; tilingSettings.MaxFeaturesPerTile = maxFeaturesPerTile; + tilingSettings.SortBy = o.SortBy; tilingSettings.Lods = lods; if (subdivisionScheme == SubdivisionScheme.QUADTREE) { @@ -278,7 +303,7 @@ private static void QuadtreeTile(string connectionString, InputTable inputTable, tile.BoundingBox = bbox.ToArray(); var outputSettings = tilesetSettings.OutputSettings; - var quadtreeTiler = new QuadtreeTiler(connectionString, inputTable, stylingSettings, tilingSettings.MaxFeaturesPerTile, tilesetSettings.Translation, outputSettings.ContentFolder, tilingSettings.Lods, tilesetSettings.Copyright, tilingSettings.SkipCreateTiles); + var quadtreeTiler = new QuadtreeTiler(connectionString, inputTable, stylingSettings, tilingSettings.MaxFeaturesPerTile, tilesetSettings.Translation, outputSettings.ContentFolder, tilingSettings.Lods, tilesetSettings.Copyright, tilingSettings.SkipCreateTiles, tilingSettings.UseImplicitTiling, sortBy: tilingSettings.SortBy); var tiles = quadtreeTiler.GenerateTiles(bbox, tile, new List(), inputTable.LodColumn != string.Empty ? tilingSettings.Lods.First() : 0, tilingSettings.CreateGltf, tilingSettings.KeepProjection); Console.WriteLine(); Console.WriteLine("Tiles created: " + tiles.Count(tile => tile.Available)); diff --git a/src/wkb2gltf.core/GeometryRecord.cs b/src/wkb2gltf.core/GeometryRecord.cs index 13c1ba9a..6e3c1617 100644 --- a/src/wkb2gltf.core/GeometryRecord.cs +++ b/src/wkb2gltf.core/GeometryRecord.cs @@ -20,6 +20,8 @@ public GeometryRecord(int batchId) public float? Radius { get; set; } + public string Hash { get; set; } + public List GetTriangles(double[] translation = null, double[] scale = null) { var triangles = GeometryProcessor.GetTriangles(Geometry, BatchId, translation, scale, Shader, Radius);