Update to DuckDB V1.2.2 + Add Spatial Join Operator

.github/workflows/StableDistributionPipeline.yml

@@ -5,13 +5,13 @@ name: Stable Extension Distribution Pipeline
 on:
   pull_request:
     branches:
-      - v1.2.1
+      -
     paths-ignore:
       - '**/README.md'
       - 'doc/**'
   push:
     branches:
-      - v1.2.1
+      - v1.2.2
     paths-ignore:
       - '**/README.md'
       - 'doc/**'
@@ -26,7 +26,7 @@ jobs:
     name: Build extension binaries
     uses: duckdb/extension-ci-tools/.github/workflows/_extension_distribution.yml@v1.2.1
     with:
-      duckdb_version: v1.2.1
+      duckdb_version: v1.2.2
       extension_name: spatial
       ci_tools_version: v1.2.1
       vcpkg_commit: 5e5d0e1cd7785623065e77eff011afdeec1a3574
@@ -39,7 +39,7 @@ jobs:
     uses: duckdb/extension-ci-tools/.github/workflows/_extension_deploy.yml@v1.2.1
     secrets: inherit
     with:
-      duckdb_version: v1.2.1
+      duckdb_version: v1.2.2
       ci_tools_version: v1.2.1
       extension_name: spatial
       deploy_latest: ${{ startsWith(github.ref, 'refs/tags/v') || github.ref == 'refs/heads/v1.2.1' }}

Makefile

@@ -1,7 +1,7 @@
 PROJ_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST))))
 
 # Configuration of extension
-EXT_NAME=excel
+EXT_NAME=spatial
 EXT_CONFIG=${PROJ_DIR}extension_config.cmake
 
 # Include the Makefile from extension-ci-tools

src/sgl/sgl.hpp

@@ -1037,7 +1037,7 @@ inline void affine_transform(sgl::allocator *alloc, sgl::geometry *geom, const s
 			}
 
 			// Now, apply the transformation
-			vertex_xyzm old_vertex = {0, 0, 0, 0};
+			vertex_xyzm old_vertex = {0, 0, 1, 1};
 			for (uint32_t i = 0; i < vertex_count; i++) {
 				memcpy(&old_vertex, old_vertex_data + i * vertex_width, vertex_width);
                 auto new_vertex = matrix->apply_xyz(old_vertex);
@@ -1781,6 +1781,83 @@ inline double haversine_distance(const double lat1_p, const double lon1_p, const
 	return R * c;
 }
 
+//----------------------------------------------------------------------------------------------------------------------
+// Hilbert Curve Encoding
+// From (Public Domain): https://github.com/rawrunprotected/hilbert_curves
+//----------------------------------------------------------------------------------------------------------------------
+inline uint32_t hilbert_interleave(uint32_t x) {
+	x = (x | (x << 8)) & 0x00FF00FF;
+	x = (x | (x << 4)) & 0x0F0F0F0F;
+	x = (x | (x << 2)) & 0x33333333;
+	x = (x | (x << 1)) & 0x55555555;
+	return x;
+}
+
+inline uint32_t hilbert_encode(uint32_t n, uint32_t x, uint32_t y) {
+	x = x << (16 - n);
+	y = y << (16 - n);
+
+	// Initial prefix scan round, prime with x and y
+	uint32_t a = x ^ y;
+	uint32_t b = 0xFFFF ^ a;
+	uint32_t c = 0xFFFF ^ (x | y);
+	uint32_t d = x & (y ^ 0xFFFF);
+	uint32_t A = a | (b >> 1);
+	uint32_t B = (a >> 1) ^ a;
+	uint32_t C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
+	uint32_t D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
+
+	a = A;
+	b = B;
+	c = C;
+	d = D;
+	A = ((a & (a >> 2)) ^ (b & (b >> 2)));
+	B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
+	C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
+	D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
+
+	a = A;
+	b = B;
+	c = C;
+	d = D;
+	A = ((a & (a >> 4)) ^ (b & (b >> 4)));
+	B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
+	C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
+	D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
+
+	// Final round and projection
+	a = A;
+	b = B;
+	c = C;
+	d = D;
+	C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
+	D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
+
+	// Undo transformation prefix scan
+	a = C ^ (C >> 1);
+	b = D ^ (D >> 1);
+
+	// Recover index bits
+	uint32_t i0 = x ^ y;
+	uint32_t i1 = b | (0xFFFF ^ (i0 | a));
+
+	return ((hilbert_interleave(i1) << 1) | hilbert_interleave(i0)) >> (32 - 2 * n);
+}
+
+inline uint32_t hilbert_f32_to_u32(float f) {
+	if (std::isnan(f)) {
+		return 0xFFFFFFFF;
+	}
+	uint32_t res;
+	memcpy(&res, &f, sizeof(res));
+	if ((res & 0x80000000) != 0) {
+		res ^= 0xFFFFFFFF;
+	} else {
+		res |= 0x80000000;
+	}
+	return res;
+}
+
 } // namespace util
 
 } // namespace sgl

src/spatial/CMakeLists.txt

@@ -2,11 +2,11 @@ add_subdirectory(util)
 add_subdirectory(modules)
 add_subdirectory(geometry)
 add_subdirectory(index)
+add_subdirectory(operators)
 
 set(EXTENSION_SOURCES
     ${EXTENSION_SOURCES}
     ${CMAKE_CURRENT_SOURCE_DIR}/spatial_extension.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/spatial_types.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/spatial_optimizers.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/spatial_geoarrow.cpp
 PARENT_SCOPE)

src/spatial/geometry/geometry_type.hpp

@@ -1,11 +1,12 @@
 #pragma once
 
+#include "duckdb/common/string_util.hpp"
+#include "duckdb/common/types/string_type.hpp"
 #include "spatial/geometry/bbox.hpp"
 #include "spatial/geometry/geometry_properties.hpp"
 #include "spatial/util/cursor.hpp"
 
-#include "duckdb/common/types/string_type.hpp"
-#include "duckdb/common/string_util.hpp"
+#include <spatial/util/math.hpp>
 
 namespace duckdb {
 
@@ -98,7 +99,7 @@ class geometry_t {
 		return props;
 	}
 
-	bool TryGetCachedBounds(Box2D<double> &bbox) const {
+	bool TryGetCachedBounds(Box2D<float> &bbox) const {
 		Cursor cursor(data);
 
 		// Read the header
@@ -135,12 +136,12 @@ class geometry_t {
 				return false;
 			}
 
-			auto x = cursor.Read<double>();
-			auto y = cursor.Read<double>();
-			bbox.min.x = x;
-			bbox.min.y = y;
-			bbox.max.x = x;
-			bbox.max.y = y;
+			const auto x = cursor.Read<double>();
+			const auto y = cursor.Read<double>();
+			bbox.min.x = MathUtil::DoubleToFloatDown(x);
+			bbox.min.y = MathUtil::DoubleToFloatDown(y);
+			bbox.max.x = MathUtil::DoubleToFloatUp(x);
+			bbox.max.y = MathUtil::DoubleToFloatUp(y);
 			return true;
 		}
 		return false;

src/spatial/index/rtree/rtree_index.cpp

@@ -166,18 +166,12 @@ ErrorData RTreeIndex::Insert(IndexLock &lock, DataChunk &input, Vector &rowid_ve
 
 		const auto rowid = rowid_data[i];
 
-		Box2D<double> box_2d;
-		if (!geom_data[i].TryGetCachedBounds(box_2d)) {
+		Box2D<float> bbox;
+		if (!geom_data[i].TryGetCachedBounds(bbox)) {
 			valid_buffer[i] = false;
 			continue;
 		}
 
-		Box2D<float> bbox;
-		bbox.min.x = MathUtil::DoubleToFloatDown(box_2d.min.x);
-		bbox.min.y = MathUtil::DoubleToFloatDown(box_2d.min.y);
-		bbox.max.x = MathUtil::DoubleToFloatUp(box_2d.max.x);
-		bbox.max.y = MathUtil::DoubleToFloatUp(box_2d.max.y);
-
 		entry_buffer[i] = {RTree::MakeRowId(rowid), bbox};
 		valid_buffer[i] = true;
 	}
@@ -227,17 +221,11 @@ void RTreeIndex::Delete(IndexLock &lock, DataChunk &input, Vector &rowid_vec) {
 		auto &geom = UnifiedVectorFormat::GetData<geometry_t>(geom_format)[geom_idx];
 		auto &rowid = UnifiedVectorFormat::GetData<row_t>(rowid_format)[rowid_idx];
 
-		Box2D<double> raw_bounds;
-		if (!geom.TryGetCachedBounds(raw_bounds)) {
+		Box2D<float> approx_bounds;
+		if (!geom.TryGetCachedBounds(approx_bounds)) {
 			continue;
 		}
 
-		Box2D<float> approx_bounds;
-		approx_bounds.min.x = MathUtil::DoubleToFloatDown(raw_bounds.min.x);
-		approx_bounds.min.y = MathUtil::DoubleToFloatDown(raw_bounds.min.y);
-		approx_bounds.max.x = MathUtil::DoubleToFloatUp(raw_bounds.max.x);
-		approx_bounds.max.y = MathUtil::DoubleToFloatUp(raw_bounds.max.y);
-
 		RTreeEntry new_entry = {RTree::MakeRowId(rowid), approx_bounds};
 		tree->Delete(new_entry);
 	}

src/spatial/index/rtree/rtree_index_plan_scan.cpp

@@ -81,20 +81,12 @@ class RTreeIndexScanOptimizer : public OptimizerExtension {
 		return true;
 	}
 
-	static bool TryGetBoundingBox(const Value &value, Box2D<float> &bbox_f) {
+	static bool TryGetBoundingBox(const Value &value, Box2D<float> &bbox) {
 		const auto str = value.GetValueUnsafe<string_t>();
 		const geometry_t blob(str);
-
-		Box2D<double> bbox;
 		if (!blob.TryGetCachedBounds(bbox)) {
 			return false;
 		}
-
-		bbox_f.min.x = MathUtil::DoubleToFloatDown(bbox.min.x);
-		bbox_f.min.y = MathUtil::DoubleToFloatDown(bbox.min.y);
-		bbox_f.max.x = MathUtil::DoubleToFloatUp(bbox.max.x);
-		bbox_f.max.y = MathUtil::DoubleToFloatUp(bbox.max.y);
-
 		return true;
 	}
 

src/spatial/modules/geos/geos_module.cpp

@@ -284,7 +284,6 @@ struct ST_Buffer {
 
 			    const auto buffer = geom.get_buffer_style(radius, segments, cap_style, join_style, mitre_limit);
 			    return lstate.Serialize(result, buffer);
-			    ;
 		    });
 	}
 

src/spatial/modules/main/spatial_functions_aggregate.cpp

@@ -43,24 +43,30 @@ struct ExtentAggFunction {
 	}
 
 	template <class INPUT_TYPE, class STATE, class OP>
-	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &) {
-		Box2D<double> bbox;
-		if (!state.is_set) {
-			if (input.TryGetCachedBounds(bbox)) {
+	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &aggregate) {
+
+		// TODO: Vectorize this so we dont clear the arena after each row
+		sgl::geometry geom;
+		Serde::Deserialize(geom, aggregate.input.allocator, input.GetDataUnsafe(), input.GetSize());
+
+		auto bbox = sgl::box_xy::smallest();
+		if (sgl::ops::try_get_extent_xy(&geom, &bbox)) {
+
+			if (!state.is_set) {
 				state.is_set = true;
 				state.xmin = bbox.min.x;
 				state.xmax = bbox.max.x;
 				state.ymin = bbox.min.y;
 				state.ymax = bbox.max.y;
-			}
-		} else {
-			if (input.TryGetCachedBounds(bbox)) {
+			} else {
 				state.xmin = std::min(state.xmin, bbox.min.x);
 				state.xmax = std::max(state.xmax, bbox.max.x);
 				state.ymin = std::min(state.ymin, bbox.min.y);
 				state.ymax = std::max(state.ymax, bbox.max.y);
 			}
 		}
+
+		aggregate.input.allocator.Reset();
 	}
 
 	template <class INPUT_TYPE, class STATE, class OP>
@@ -101,8 +107,7 @@ struct ExtentAggFunction {
 			Serde::Serialize(bbox, blob.GetDataWriteable(), size);
 			blob.Finalize();
 
-			// TODO: dont use geometry_t here
-			target = geometry_t(blob);
+			target = blob;
 		}
 	}
 
@@ -137,7 +142,7 @@ static constexpr const char *DOC_ALIAS_DESCRIPTION = R"(
 void RegisterSpatialAggregateFunctions(DatabaseInstance &db) {
 
 	// TODO: Dont use geometry_t here
-	const auto agg = AggregateFunction::UnaryAggregate<ExtentAggState, geometry_t, geometry_t, ExtentAggFunction>(
+	const auto agg = AggregateFunction::UnaryAggregate<ExtentAggState, string_t, string_t, ExtentAggFunction>(
 	    GeoTypes::GEOMETRY(), GeoTypes::GEOMETRY());
 
 	FunctionBuilder::RegisterAggregate(db, "ST_Extent_Agg", [&](AggregateFunctionBuilder &func) {