Size-Based Geometry Prioritization for Tile Generation

README.md

@@ -220,10 +220,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({inputTable.GeometryColumn})::text));
 ```
 
 When there the spatial index is not present the following warning is shown.

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)) && <envelope>
+      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)), <envelope>)
+
+  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)

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, bool keepProjection = false, HashSet<string> excludeHashes = null)
     {
         var select = $"COUNT({geometry_column})";
         var where = GeometryRepository.GetWhere(geometry_column, from, to, query, source_epsg, keepProjection);
+
+        // 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();
+        }
     }
 }

src/b3dm.tileset/GeometryRepository.cs

@@ -13,32 +13,79 @@ namespace B3dm.Tileset;
 
 public static class GeometryRepository
 {
+    public static HashSet<string> FilterHashesByEnvelope(NpgsqlConnection conn, string tableName, string geometryColumn, BoundingBox bbox, int source_epsg, HashSet<string> geometryHashes, bool keepProjection)
+    {
+        if (geometryHashes.Count == 0) {
+            return new HashSet<string>();
+        }
+
+        var filteredHashes = new HashSet<string>();
+
+        conn.Open();
+        try {
+            var envelope = keepProjection ?
+                $"ST_MakeEnvelope(@xmin, @ymin, @xmax, @ymax, {source_epsg})" :
+                $"ST_Transform(ST_MakeEnvelope(@xmin, @ymin, @xmax, @ymax, 4326), {source_epsg})";
+
+            var query = $@"
+            SELECT MD5(ST_AsBinary({geometryColumn})::text) as geom_hash
+            FROM {tableName}
+            WHERE MD5(ST_AsBinary({geometryColumn})::text) = ANY(@hashes)
+            AND ST_Within(
+                ST_Centroid(ST_Envelope({geometryColumn})),
+                {envelope}
+            )";
+
+            using var cmd = new NpgsqlCommand(query, conn);
+            cmd.Parameters.AddWithValue("hashes", geometryHashes.ToArray());
+            cmd.Parameters.AddWithValue("xmin", bbox.XMin);
+            cmd.Parameters.AddWithValue("ymin", bbox.YMin);
+            cmd.Parameters.AddWithValue("xmax", bbox.XMax);
+            cmd.Parameters.AddWithValue("ymax", bbox.YMax);
+
+            using var reader = cmd.ExecuteReader();
+            while (reader.Read()) {
+                var hash = reader.GetString(0);
+                filteredHashes.Add(hash);
+            }
+
+            return filteredHashes;
+        }
+        finally {
+            conn.Close();
+        }
+    }
+
     /// <summary>
     /// Returns double array with 6 bounding box coordinates, xmin, ymin, xmax, ymax, zmin, zmax
     /// </summary>
-    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<string> 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<WkbSerializer>(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<WkbSerializer>(stream);
+            var result = BBox3D.GetBoundingBoxPoints(geometry);
 
-        return result;
+            return result;
+        }
+        finally {
+            conn.Close();
+        }
     }
 
-    public static List<GeometryRecord> 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<GeometryRecord> 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, HashSet<string> excludeHashes = null, int? maxFeatures = null)
     {
         var sqlselect = GetSqlSelect(geometry_column, shaderColumn, attributesColumns, radiusColumn, target_srs);
         var sqlFrom = "FROM " + geometry_table;
@@ -47,10 +94,29 @@ public static List<GeometryRecord> GetGeometrySubset(NpgsqlConnection conn, stri
         var points = GetPoints(bbox);
 
         var sqlWhere = GetWhere(geometry_column, points.fromPoint, points.toPoint, query, source_epsg, keepProjection);
-        var sql = sqlselect + sqlFrom + " where " + sqlWhere;
+
+        // 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);
+        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)
@@ -98,6 +164,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 +175,19 @@ public static string GetGeometryColumn(string geometry_column, int target_srs)
         return $"st_transform({geometry_column}, {target_srs})";
     }
 
-    public static List<GeometryRecord> GetGeometries(NpgsqlConnection conn, string shaderColumn, string attributesColumns, string sql, string radiusColumn)
+    public static string GetOrderBy(string geometry_column)
+    {
+        return $" ORDER BY ST_Area(ST_Envelope({geometry_column})) DESC";
+    }
+
+    public static List<GeometryRecord> GetGeometries(NpgsqlCommand cmd, string shaderColumn, string attributesColumns, string radiusColumn, string geometry_column = "")
     {
         var geometries = new List<GeometryRecord>();
-        conn.Open();
-        var cmd = new NpgsqlCommand(sql, conn);
         var reader = cmd.ExecuteReader();
         var attributesColumnIds = new Dictionary<string, int>();
         var shadersColumnId = int.MinValue;
         var radiusColumnId = int.MinValue;
+        var hashColumnId = int.MinValue;
 
         if (attributesColumns != String.Empty) {
             var attributesColumnsList = attributesColumns.Split(',').ToList();
@@ -133,6 +205,11 @@ public static List<GeometryRecord> 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 +231,15 @@ public static List<GeometryRecord> GetGeometries(NpgsqlConnection conn, string s
                 var radius = reader.GetFieldValue<object>(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;
     }