diff --git a/c/sedona-geos/benches/geos-functions.rs b/c/sedona-geos/benches/geos-functions.rs index 3dca26dfe1..d535031405 100644 --- a/c/sedona-geos/benches/geos-functions.rs +++ b/c/sedona-geos/benches/geos-functions.rs @@ -28,9 +28,6 @@ fn criterion_benchmark(c: &mut Criterion) { benchmark::scalar(c, &f, "geos", "st_area", Polygon(10)); benchmark::scalar(c, &f, "geos", "st_area", Polygon(500)); - benchmark::scalar(c, &f, "geos", "st_boundary", Polygon(10)); - benchmark::scalar(c, &f, "geos", "st_boundary", Polygon(500)); - benchmark::scalar( c, &f, diff --git a/c/sedona-geos/src/lib.rs b/c/sedona-geos/src/lib.rs index b77236e74f..036187ccbe 100644 --- a/c/sedona-geos/src/lib.rs +++ b/c/sedona-geos/src/lib.rs @@ -23,7 +23,6 @@ mod geos_to_wkb; mod overlay; pub mod register; mod st_area; -mod st_boundary; mod st_buffer; mod st_buildarea; mod st_centroid; diff --git a/c/sedona-geos/src/register.rs b/c/sedona-geos/src/register.rs index f28e2a131f..c016f90d38 100644 --- a/c/sedona-geos/src/register.rs +++ b/c/sedona-geos/src/register.rs @@ -40,7 +40,6 @@ macro_rules! define_aggregate_kernels { pub fn scalar_kernels() -> Vec<(&'static str, Vec)> { define_scalar_kernels!( "st_area" => crate::st_area::st_area_impl, - "st_boundary" => crate::st_boundary::st_boundary_impl, "st_buildarea" => crate::st_buildarea::st_build_area_impl, "st_buffer" => crate::st_buffer::st_buffer_impl, "st_buffer" => crate::st_buffer::st_buffer_style_impl, diff --git a/c/sedona-geos/src/st_boundary.rs b/c/sedona-geos/src/st_boundary.rs deleted file mode 100644 index 4c41c4574d..0000000000 --- a/c/sedona-geos/src/st_boundary.rs +++ /dev/null @@ -1,393 +0,0 @@ -// Licensed to the Apache Software Foundation (ASF) under one -// or more contributor license agreements. See the NOTICE file -// distributed with this work for additional information -// regarding copyright ownership. The ASF licenses this file -// to you under the Apache License, Version 2.0 (the -// "License"); you may not use this file except in compliance -// with the License. You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, -// software distributed under the License is distributed on an -// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY -// KIND, either express or implied. See the License for the -// specific language governing permissions and limitations -// under the License. - -use std::sync::Arc; - -use arrow_array::builder::BinaryBuilder; -use datafusion_common::{error::Result, DataFusionError}; -use datafusion_expr::ColumnarValue; -use geos::{Geom, Geometry, GeometryTypes}; -use sedona_expr::{ - item_crs::ItemCrsKernel, - scalar_udf::{ScalarKernelRef, SedonaScalarKernel}, -}; -use sedona_geometry::wkb_factory::WKB_MIN_PROBABLE_BYTES; -use sedona_schema::{ - datatypes::{SedonaType, WKB_GEOGRAPHY, WKB_GEOMETRY}, - matchers::ArgMatcher, -}; - -use crate::executor::GeosExecutor; -use crate::geos_to_wkb::write_geos_geometry; - -/// ST_Boundary() implementation using the geos crate -pub fn st_boundary_impl() -> Vec { - ItemCrsKernel::wrap_impl(vec![ - Arc::new(STBoundary { - matcher: ArgMatcher::new(vec![ArgMatcher::is_geometry()], WKB_GEOMETRY), - }), - Arc::new(STBoundary { - matcher: ArgMatcher::new(vec![ArgMatcher::is_geography()], WKB_GEOGRAPHY), - }), - ]) -} - -#[derive(Debug)] -struct STBoundary { - matcher: ArgMatcher, -} - -impl SedonaScalarKernel for STBoundary { - fn return_type(&self, args: &[SedonaType]) -> datafusion_common::Result> { - self.matcher.match_args(args) - } - - fn invoke_batch( - &self, - arg_types: &[SedonaType], - args: &[ColumnarValue], - ) -> datafusion_common::Result { - let executor = GeosExecutor::new(arg_types, args); - let mut builder = BinaryBuilder::with_capacity( - executor.num_iterations(), - WKB_MIN_PROBABLE_BYTES * executor.num_iterations(), - ); - - executor.execute_wkb_void(|maybe_wkb| { - match maybe_wkb { - Some(wkb) => { - invoke_scalar(wkb, &mut builder)?; - builder.append_value([]); - } - _ => builder.append_null(), - } - Ok(()) - })?; - - executor.finish(Arc::new(builder.finish())) - } -} - -fn invoke_scalar(geos_geom: &geos::Geometry, writer: &mut BinaryBuilder) -> Result<()> { - let result_geom = geos_boundary(geos_geom)?; - - write_geos_geometry(&result_geom, writer)?; - Ok(()) -} - -/// For simple geometries, it calls `geometry.boundary()`. -/// For a `GeometryCollection`, it recursively computes the boundary for each -/// component and then re-aggregates the resulting geometry types (Point, LineString, etc.) -/// into their respective Multi-geometry types (MultiPoint, MultiLineString, etc.) -/// before forming the final `GeometryCollection`. This aggregation step is crucial -/// for adhering to OGC specifications for `GeometryCollection` boundaries. -fn geos_boundary(geometry: &impl Geom) -> Result { - if geometry - .geometry_type() - .map_err(|e| DataFusionError::Execution(format!("Failed to get geometry type: {e}")))? - == GeometryTypes::GeometryCollection - { - let num_geometries = geometry.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of geometries: {e}")) - })?; - - let mut empty_collections: Vec = Vec::new(); - let mut points: Vec = Vec::new(); - let mut lines: Vec = Vec::new(); - let mut polygons: Vec = Vec::new(); - - for i in 0..num_geometries { - let child_geom = geometry.get_geometry_n(i).map_err(|e| { - DataFusionError::Execution(format!("Failed to get {}th child geometry: {e}", i + 1)) - })?; - - // Recursively calculate the boundary of the child geometry - let child_boundary = geos_boundary(&child_geom)?; - - // Collect components based on whether they are an empty GeometryCollection - if is_empty_geometry_collection(&child_boundary)? { - empty_collections.push(child_boundary); - } else { - // Collect and group non-empty boundary components (Points, LineStrings, etc.) - collect_boundary_components( - &child_boundary, - &mut points, - &mut lines, - &mut polygons, - )?; - } - } - - let mut result_components: Vec = Vec::new(); - - // Aggregate the result - result_components.extend(empty_collections); - - if points.len() == 1 { - result_components.push(points.into_iter().next().unwrap()); - } else if !points.is_empty() { - let multi_point = Geometry::create_multipoint(points).map_err(|e| { - DataFusionError::Execution(format!("Failed to create multipoint: {e}")) - })?; - result_components.push(multi_point); - } - - if lines.len() == 1 { - result_components.push(lines.into_iter().next().unwrap()); - } else if !lines.is_empty() { - let multi_line = Geometry::create_multiline_string(lines).map_err(|e| { - DataFusionError::Execution(format!("Failed to create multilinestring: {e}")) - })?; - result_components.push(multi_line); - } - - if polygons.len() == 1 { - result_components.push(polygons.into_iter().next().unwrap()); - } else if !polygons.is_empty() { - let multi_polygon = Geometry::create_multipolygon(polygons).map_err(|e| { - DataFusionError::Execution(format!("Failed to create multipolygon: {e}")) - })?; - result_components.push(multi_polygon); - } - - if result_components.len() == 1 { - Ok(Geom::clone(result_components.first().unwrap()) - .map_err(|e| DataFusionError::Execution(e.to_string()))?) - } else { - Geometry::create_geometry_collection(result_components).map_err(|e| { - DataFusionError::Execution(format!("Failed to create geometry collection: {e}")) - }) - } - } else { - // For simple geometries, use the standard geos boundary function - geometry - .boundary() - .map_err(|e| DataFusionError::Execution(format!("Failed to calculate boundary: {e}"))) - } -} - -/// Checks if a geometry is an empty `GeometryCollection`. -fn is_empty_geometry_collection(geom: &Geometry) -> Result { - if geom - .geometry_type() - .map_err(|e| DataFusionError::Execution(format!("Failed to get geometry type: {e}")))? - == GeometryTypes::GeometryCollection - { - let num = geom.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of geometries: {e}")) - })?; - Ok(num == 0) - } else { - Ok(false) - } -} - -/// Recursively collects and groups individual boundary components (Point, LineString, etc.) -/// from a given boundary geometry (which could be a `GeometryCollection` itself) -/// into mutable vectors based on their type. -fn collect_boundary_components( - boundary: &Geometry, - points: &mut Vec, - lines: &mut Vec, - polygons: &mut Vec, -) -> Result<()> { - match boundary - .geometry_type() - .map_err(|e| DataFusionError::Execution(format!("Failed to get geometry type: {e}")))? - { - // Recurse into sub-geometries if it's a collection - GeometryTypes::GeometryCollection => { - let num_geoms = boundary.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of geometries: {e}")) - })?; - - (0..num_geoms).try_for_each(|i| { - let component = boundary - .get_geometry_n(i) - .and_then(|g| Geom::clone(&g)) - .map_err(|e| { - DataFusionError::Execution(format!( - "Failed to process {}th geometry: {e}", - i + 1 - )) - })?; - collect_boundary_components(&component, points, lines, polygons) - })?; - } - // Collect simple, single-part components - GeometryTypes::Point => { - points.push( - Geom::clone(boundary).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - GeometryTypes::LineString => { - lines.push( - Geom::clone(boundary).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - GeometryTypes::Polygon => { - polygons.push( - Geom::clone(boundary).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - // Decompose Multi-geometries and collect their parts - GeometryTypes::MultiPoint => { - let num_points = boundary.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of points: {e}")) - })?; - for i in 0..num_points { - let point = boundary.get_geometry_n(i).map_err(|e| { - DataFusionError::Execution(format!("Failed to get {}th point: {e}", i + 1)) - })?; - points.push( - Geom::clone(&point).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - } - GeometryTypes::MultiLineString => { - let num_lines = boundary.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of linestrings: {e}")) - })?; - for i in 0..num_lines { - let line = boundary.get_geometry_n(i).map_err(|e| { - DataFusionError::Execution(format!("Failed to get {}th linestring: {e}", i + 1)) - })?; - lines.push( - Geom::clone(&line).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - } - GeometryTypes::MultiPolygon => { - let num_polygons = boundary.get_num_geometries().map_err(|e| { - DataFusionError::Execution(format!("Failed to get number of polygons: {e}")) - })?; - for i in 0..num_polygons { - let polygon = boundary.get_geometry_n(i).map_err(|e| { - DataFusionError::Execution(format!("Failed to get {}th polygon: {e}", i + 1)) - })?; - polygons.push( - Geom::clone(&polygon).map_err(|e| DataFusionError::Execution(e.to_string()))?, - ); - } - } - // Ignore other types (e.g., empty geometries) - _ => {} - } - Ok(()) -} - -#[cfg(test)] -mod tests { - use rstest::rstest; - use sedona_expr::scalar_udf::SedonaScalarUDF; - use sedona_schema::datatypes::{WKB_GEOGRAPHY_ITEM_CRS, WKB_GEOMETRY_ITEM_CRS}; - use sedona_testing::testers::ScalarUdfTester; - - use super::*; - - #[rstest] - fn udf(#[values(WKB_GEOMETRY, WKB_GEOGRAPHY)] sedona_type: SedonaType) { - let udf = SedonaScalarUDF::from_impl("st_boundary", st_boundary_impl()); - let tester = ScalarUdfTester::new(udf.into(), vec![sedona_type.clone()]); - tester.assert_return_type(sedona_type.clone()); - - let result = tester - .invoke_scalar( - "GEOMETRYCOLLECTION(LINESTRING(1 1,2 2),GEOMETRYCOLLECTION(POLYGON((3 3,4 4,5 5,3 3)),GEOMETRYCOLLECTION(LINESTRING(6 6,7 7),POLYGON((8 8,9 9,10 10,8 8)))))", - ) - .unwrap(); - tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION(MULTIPOINT((1 1),(2 2),(6 6),(7 7)),MULTILINESTRING((3 3,4 4,5 5,3 3),(8 8,9 9,10 10,8 8)))"); - - let result = tester - .invoke_scalar("LINESTRING(100 150,50 60, 70 80, 160 170)") - .unwrap(); - tester.assert_scalar_result_equals(result, "MULTIPOINT((100 150),(160 170))"); - - let result = tester - .invoke_scalar( - "POLYGON (( 10 130, 50 190, 110 190, 140 150, 150 80, 100 10, 20 40, 10 130 ), ( 70 40, 100 50, 120 80, 80 110, 50 90, 70 40 ))" - ) - .unwrap(); - tester.assert_scalar_result_equals( - result, - "MULTILINESTRING((10 130,50 190,110 190,140 150,150 80,100 10,20 40,10 130), (70 40,100 50,120 80,80 110,50 90,70 40))" - ); - - let result = tester - .invoke_scalar("MULTILINESTRING ((10 10, 20 20), (30 30, 40 40, 30 30))") - .unwrap(); - tester.assert_scalar_result_equals(result, "MULTIPOINT (10 10, 20 20)"); - - let result = tester.invoke_scalar("GEOMETRYCOLLECTION(MULTIPOINT(-2 3, -2 2), LINESTRING(5 5, 10 10), POLYGON((-7 4.2, -7.1 5, -7.1 4.3, -7 4.2)))").unwrap(); - tester.assert_scalar_result_equals( - result, - "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION EMPTY, MULTIPOINT(5 5, 10 10), LINESTRING(-7 4.2, -7.1 5, -7.1 4.3, -7 4.2))" - ); - - let result = tester.invoke_scalar("POINT (10 20)").unwrap(); - tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); - - let result = tester - .invoke_scalar("MULTIPOINT (5 5, 10 10, 15 15)") - .unwrap(); - tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); - - let result = tester - .invoke_scalar("LINESTRING (0 0, 1 1, 0 1, 0 0)") - .unwrap(); - tester.assert_scalar_result_equals(result, "MULTIPOINT EMPTY"); - - let result = tester - .invoke_scalar("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))") - .unwrap(); - tester.assert_scalar_result_equals(result, "LINESTRING(0 0,0 10,10 10,10 0,0 0)"); - - let result = tester - .invoke_scalar( - "MULTIPOLYGON (((0 0, 0 1, 1 1, 1 0, 0 0)), ((10 10, 10 11, 11 11, 11 10, 10 10)))", - ) - .unwrap(); - tester.assert_scalar_result_equals( - result, - "MULTILINESTRING((0 0,0 1,1 1,1 0,0 0),(10 10,10 11,11 11,11 10,10 10))", - ); - - let result = tester - .invoke_scalar("GEOMETRYCOLLECTION(POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0)), GEOMETRYCOLLECTION(LINESTRING(10 10, 10 20)))") - .unwrap(); - tester.assert_scalar_result_equals( - result, - "GEOMETRYCOLLECTION(MULTIPOINT((10 10),(10 20)), LINESTRING(0 0,0 1,1 1,1 0,0 0))", - ); - - let result = tester.invoke_scalar("GEOMETRYCOLLECTION EMPTY").unwrap(); - tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); - } - - #[rstest] - fn udf_invoke_item_crs( - #[values(WKB_GEOMETRY_ITEM_CRS.clone(), WKB_GEOGRAPHY_ITEM_CRS.clone())] - sedona_type: SedonaType, - ) { - let udf = SedonaScalarUDF::from_impl("st_boundary", st_boundary_impl()); - let tester = ScalarUdfTester::new(udf.into(), vec![sedona_type.clone()]); - tester.assert_return_type(sedona_type); - - let result = tester.invoke_scalar("POINT (1 3)").unwrap(); - tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); - } -} diff --git a/rust/sedona-functions/benches/native-functions.rs b/rust/sedona-functions/benches/native-functions.rs index ff90660ace..5487d6af27 100644 --- a/rust/sedona-functions/benches/native-functions.rs +++ b/rust/sedona-functions/benches/native-functions.rs @@ -27,6 +27,9 @@ fn criterion_benchmark(c: &mut Criterion) { benchmark::scalar(c, &f, "native", "st_astext", Point); benchmark::scalar(c, &f, "native", "st_astext", LineString(10)); + benchmark::scalar(c, &f, "native", "st_boundary", Polygon(10)); + benchmark::scalar(c, &f, "native", "st_boundary", Polygon(500)); + benchmark::scalar(c, &f, "native", "st_dimension", Point); benchmark::scalar(c, &f, "native", "st_dimension", LineString(10)); diff --git a/rust/sedona-functions/src/lib.rs b/rust/sedona-functions/src/lib.rs index c848aa01d7..9a34a2ab27 100644 --- a/rust/sedona-functions/src/lib.rs +++ b/rust/sedona-functions/src/lib.rs @@ -28,6 +28,7 @@ mod st_asbinary; mod st_asewkb; mod st_astext; mod st_azimuth; +mod st_boundary; mod st_collect_agg; mod st_dimension; mod st_dump; diff --git a/rust/sedona-functions/src/register.rs b/rust/sedona-functions/src/register.rs index c01d2e1196..72ad3ecfca 100644 --- a/rust/sedona-functions/src/register.rs +++ b/rust/sedona-functions/src/register.rs @@ -46,6 +46,7 @@ pub fn default_function_set() -> FunctionSet { crate::st_asewkb::st_asewkb_udf, crate::st_astext::st_astext_udf, crate::st_azimuth::st_azimuth_udf, + crate::st_boundary::st_boundary_udf, crate::st_dimension::st_dimension_udf, crate::st_dump::st_dump_udf, crate::st_envelope::st_envelope_udf, diff --git a/rust/sedona-functions/src/st_boundary.rs b/rust/sedona-functions/src/st_boundary.rs new file mode 100644 index 0000000000..e2a4a55b5c --- /dev/null +++ b/rust/sedona-functions/src/st_boundary.rs @@ -0,0 +1,459 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use std::{collections::HashMap, io::Write, sync::Arc}; + +use arrow_array::builder::BinaryBuilder; +use datafusion_common::{error::Result, exec_datafusion_err}; +use datafusion_expr::{ColumnarValue, Volatility}; +use geo_traits::{ + CoordTrait, GeometryCollectionTrait, GeometryTrait, LineStringTrait, MultiLineStringTrait, + MultiPolygonTrait, PolygonTrait, +}; +use sedona_expr::{ + item_crs::ItemCrsKernel, + scalar_udf::{SedonaScalarKernel, SedonaScalarUDF}, +}; +use sedona_geometry::{ + error::SedonaGeometryError, + wkb_factory::{ + write_wkb_geometrycollection_header, write_wkb_linestring_header, + write_wkb_multilinestring_header, write_wkb_multipoint_header, write_wkb_point_header, + WKB_MIN_PROBABLE_BYTES, + }, +}; +use sedona_schema::{ + datatypes::{SedonaType, WKB_GEOGRAPHY, WKB_GEOMETRY}, + matchers::ArgMatcher, +}; +use wkb::{ + reader::{Coord, LinearRing, Wkb}, + Endianness, +}; + +use crate::executor::WkbExecutor; + +/// ST_Boundary() scalar UDF implementation using geo-traits over WKB +pub fn st_boundary_udf() -> SedonaScalarUDF { + SedonaScalarUDF::new( + "st_boundary", + ItemCrsKernel::wrap_impl(vec![ + Arc::new(STBoundary { + matcher: ArgMatcher::new(vec![ArgMatcher::is_geometry()], WKB_GEOMETRY), + }), + Arc::new(STBoundary { + matcher: ArgMatcher::new(vec![ArgMatcher::is_geography()], WKB_GEOGRAPHY), + }), + ]), + Volatility::Immutable, + ) +} + +#[derive(Debug)] +struct STBoundary { + matcher: ArgMatcher, +} + +impl SedonaScalarKernel for STBoundary { + fn return_type(&self, args: &[SedonaType]) -> datafusion_common::Result> { + self.matcher.match_args(args) + } + + fn invoke_batch( + &self, + arg_types: &[SedonaType], + args: &[ColumnarValue], + ) -> datafusion_common::Result { + let executor = WkbExecutor::new(arg_types, args); + let mut builder = BinaryBuilder::with_capacity( + executor.num_iterations(), + WKB_MIN_PROBABLE_BYTES * executor.num_iterations(), + ); + + executor.execute_wkb_void(|maybe_wkb| { + match maybe_wkb { + Some(wkb) => { + invoke_scalar(wkb, &mut builder) + .map_err(|e| exec_datafusion_err!("ST_Boundary error: {e}"))?; + builder.append_value([]); + } + _ => builder.append_null(), + } + Ok(()) + })?; + + executor.finish(Arc::new(builder.finish())) + } +} + +fn invoke_scalar(geom: &Wkb, writer: &mut impl Write) -> Result<(), SedonaGeometryError> { + boundary(geom)?.write(writer) +} + +#[derive(Debug)] +enum Boundary<'a> { + EmptyCollection(geo_traits::Dimensions), + Point(Coord<'a>), + MultiPoint(geo_traits::Dimensions, Vec>), + Line(LinearRing<'a>), + MultiLine(geo_traits::Dimensions, Vec>), + Collection(geo_traits::Dimensions, Vec>), +} + +impl Boundary<'_> { + fn write(&self, writer: &mut impl Write) -> Result<(), SedonaGeometryError> { + match self { + Self::EmptyCollection(dims) => write_wkb_geometrycollection_header(writer, *dims, 0), + Self::Point(coord) => { + write_wkb_point_header(writer, coord.dim())?; + write_coord(writer, coord) + } + Self::MultiPoint(dims, coords) => { + write_wkb_multipoint_header(writer, *dims, coords.len())?; + for coord in coords { + write_wkb_point_header(writer, coord.dim())?; + write_coord(writer, coord)?; + } + Ok(()) + } + Self::Line(ring) => write_ring(writer, ring), + Self::MultiLine(dims, rings) => { + write_wkb_multilinestring_header(writer, *dims, rings.len())?; + for ring in rings { + write_ring(writer, ring)?; + } + Ok(()) + } + Self::Collection(dims, components) => { + write_wkb_geometrycollection_header(writer, *dims, components.len())?; + for component in components { + component.write(writer)?; + } + Ok(()) + } + } + } +} + +fn boundary<'a>(geom: &'a Wkb<'a>) -> Result, SedonaGeometryError> { + let dims = geom.dim(); + match geom.as_type() { + geo_traits::GeometryType::Point(_) | geo_traits::GeometryType::MultiPoint(_) => { + Ok(Boundary::EmptyCollection(dims)) + } + geo_traits::GeometryType::LineString(line) => { + Ok(Boundary::MultiPoint(dims, line_boundary(line))) + } + geo_traits::GeometryType::MultiLineString(lines) => { + Ok(Boundary::MultiPoint(dims, multiline_boundary(lines))) + } + geo_traits::GeometryType::Polygon(polygon) => Ok(polygon_boundary(polygon, dims)), + geo_traits::GeometryType::MultiPolygon(polygons) => { + let mut rings = Vec::new(); + for polygon in polygons.polygons() { + collect_polygon_rings(polygon, &mut rings); + } + Ok(Boundary::MultiLine(dims, rings)) + } + geo_traits::GeometryType::GeometryCollection(collection) => { + geometry_collection_boundary(collection, dims) + } + _ => Err(SedonaGeometryError::Invalid( + "Unsupported geometry type for boundary operation".to_string(), + )), + } +} + +fn line_boundary<'a>(line: &wkb::reader::LineString<'a>) -> Vec> { + match line.num_coords() { + 0 => vec![], + count => { + let first = line.coord(0).unwrap(); + let last = line.coord(count - 1).unwrap(); + if coord_key(&first) == coord_key(&last) { + vec![] + } else { + vec![first, last] + } + } + } +} + +fn multiline_boundary<'a>(lines: &wkb::reader::MultiLineString<'a>) -> Vec> { + let mut indices = HashMap::new(); + let mut endpoints: Vec<(Coord<'a>, bool)> = Vec::new(); + + for line in lines.line_strings() { + if line.num_coords() == 0 { + continue; + } + + let first = line.coord(0).unwrap(); + let last = line.coord(line.num_coords() - 1).unwrap(); + toggle_endpoint(first, &mut indices, &mut endpoints); + toggle_endpoint(last, &mut indices, &mut endpoints); + } + + endpoints + .into_iter() + .filter_map(|(coord, is_boundary)| is_boundary.then_some(coord)) + .collect() +} + +fn toggle_endpoint<'a>( + coord: Coord<'a>, + indices: &mut HashMap<(u64, u64), usize>, + endpoints: &mut Vec<(Coord<'a>, bool)>, +) { + let key = coord_key(&coord); + if let Some(index) = indices.get(&key) { + endpoints[*index].1 = !endpoints[*index].1; + } else { + indices.insert(key, endpoints.len()); + endpoints.push((coord, true)); + } +} + +fn coord_key(coord: &impl CoordTrait) -> (u64, u64) { + // Normalize signed zero because -0 and +0 are the same topological position. + let x = if coord.x() == 0.0 { 0.0 } else { coord.x() }; + let y = if coord.y() == 0.0 { 0.0 } else { coord.y() }; + (x.to_bits(), y.to_bits()) +} + +fn polygon_boundary<'a>( + polygon: &wkb::reader::Polygon<'a>, + dims: geo_traits::Dimensions, +) -> Boundary<'a> { + let mut rings = Vec::with_capacity(polygon.num_interiors() + 1); + collect_polygon_rings(polygon, &mut rings); + if rings.len() == 1 { + Boundary::Line(rings[0]) + } else { + Boundary::MultiLine(dims, rings) + } +} + +fn collect_polygon_rings<'a>(polygon: &wkb::reader::Polygon<'a>, rings: &mut Vec>) { + if let Some(exterior) = polygon.exterior() { + rings.push(*exterior); + } + rings.extend(polygon.interiors().copied()); +} + +fn geometry_collection_boundary<'a>( + collection: &'a wkb::reader::GeometryCollection<'a>, + dims: geo_traits::Dimensions, +) -> Result, SedonaGeometryError> { + let mut empty_collections = Vec::new(); + let mut points = Vec::new(); + let mut lines = Vec::new(); + + for geom in collection.geometries() { + let child = boundary(geom)?; + if matches!(child, Boundary::EmptyCollection(_)) { + empty_collections.push(child); + } else { + collect_components(child, &mut points, &mut lines); + } + } + + let mut components = empty_collections; + match points.len() { + 0 => {} + 1 => components.push(Boundary::Point(points[0])), + _ => components.push(Boundary::MultiPoint(dims, points)), + } + match lines.len() { + 0 => {} + 1 => components.push(Boundary::Line(lines[0])), + _ => components.push(Boundary::MultiLine(dims, lines)), + } + + Ok(match components.len() { + 0 => Boundary::EmptyCollection(dims), + 1 => components.pop().unwrap(), + _ => Boundary::Collection(dims, components), + }) +} + +fn collect_components<'a>( + boundary: Boundary<'a>, + points: &mut Vec>, + lines: &mut Vec>, +) { + match boundary { + Boundary::Point(point) => points.push(point), + Boundary::MultiPoint(_, mut child_points) => points.append(&mut child_points), + Boundary::Line(line) => lines.push(line), + Boundary::MultiLine(_, mut child_lines) => lines.append(&mut child_lines), + Boundary::Collection(_, components) => { + for component in components { + // Empty collections nested inside a non-empty collection are ignored, + // matching the existing GEOS-backed aggregation behaviour. + if !matches!(component, Boundary::EmptyCollection(_)) { + collect_components(component, points, lines); + } + } + } + Boundary::EmptyCollection(_) => {} + } +} + +fn write_ring(writer: &mut impl Write, ring: &LinearRing<'_>) -> Result<(), SedonaGeometryError> { + write_wkb_linestring_header(writer, ring.dim(), ring.num_coords())?; + write_coords(writer, ring.coords_slice(), ring.byte_order()) +} + +fn write_coord(writer: &mut impl Write, coord: &Coord<'_>) -> Result<(), SedonaGeometryError> { + write_coords(writer, coord.coord_slice(), coord.byte_order()) +} + +fn write_coords( + writer: &mut impl Write, + coords: &[u8], + byte_order: Endianness, +) -> Result<(), SedonaGeometryError> { + if matches!(byte_order, Endianness::LittleEndian) { + writer.write_all(coords)?; + } else { + for ordinate in coords.as_chunks::<{ size_of::() }>().0 { + let mut little_endian = [0; size_of::()]; + little_endian.copy_from_slice(ordinate); + little_endian.reverse(); + writer.write_all(&little_endian)?; + } + } + Ok(()) +} + +#[cfg(test)] +mod tests { + use rstest::rstest; + use sedona_schema::datatypes::{WKB_GEOGRAPHY_ITEM_CRS, WKB_GEOMETRY_ITEM_CRS}; + use sedona_testing::testers::ScalarUdfTester; + + use super::*; + + #[rstest] + fn udf(#[values(WKB_GEOMETRY, WKB_GEOGRAPHY)] sedona_type: SedonaType) { + let tester = ScalarUdfTester::new(st_boundary_udf().into(), vec![sedona_type.clone()]); + tester.assert_return_type(sedona_type.clone()); + + let result = tester + .invoke_scalar( + "GEOMETRYCOLLECTION(LINESTRING(1 1,2 2),GEOMETRYCOLLECTION(POLYGON((3 3,4 4,5 5,3 3)),GEOMETRYCOLLECTION(LINESTRING(6 6,7 7),POLYGON((8 8,9 9,10 10,8 8)))))", + ) + .unwrap(); + tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION(MULTIPOINT((1 1),(2 2),(6 6),(7 7)),MULTILINESTRING((3 3,4 4,5 5,3 3),(8 8,9 9,10 10,8 8)))"); + + let result = tester + .invoke_scalar("LINESTRING(100 150,50 60, 70 80, 160 170)") + .unwrap(); + tester.assert_scalar_result_equals(result, "MULTIPOINT((100 150),(160 170))"); + + let result = tester + .invoke_scalar( + "POLYGON (( 10 130, 50 190, 110 190, 140 150, 150 80, 100 10, 20 40, 10 130 ), ( 70 40, 100 50, 120 80, 80 110, 50 90, 70 40 ))" + ) + .unwrap(); + tester.assert_scalar_result_equals( + result, + "MULTILINESTRING((10 130,50 190,110 190,140 150,150 80,100 10,20 40,10 130), (70 40,100 50,120 80,80 110,50 90,70 40))" + ); + + let result = tester + .invoke_scalar("MULTILINESTRING ((10 10, 20 20), (30 30, 40 40, 30 30))") + .unwrap(); + tester.assert_scalar_result_equals(result, "MULTIPOINT (10 10, 20 20)"); + + // Endpoints shared by an even number of components are not in the + // boundary (the OGC mod-2 boundary node rule). + let result = tester + .invoke_scalar("MULTILINESTRING ((0 0, 1 1), (1 1, 2 2))") + .unwrap(); + tester.assert_scalar_result_equals(result, "MULTIPOINT (0 0, 2 2)"); + + let result = tester.invoke_scalar("GEOMETRYCOLLECTION(MULTIPOINT(-2 3, -2 2), LINESTRING(5 5, 10 10), POLYGON((-7 4.2, -7.1 5, -7.1 4.3, -7 4.2)))").unwrap(); + tester.assert_scalar_result_equals( + result, + "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION EMPTY, MULTIPOINT(5 5, 10 10), LINESTRING(-7 4.2, -7.1 5, -7.1 4.3, -7 4.2))" + ); + + let result = tester.invoke_scalar("POINT (10 20)").unwrap(); + tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); + + let result = tester + .invoke_scalar("MULTIPOINT (5 5, 10 10, 15 15)") + .unwrap(); + tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); + + let result = tester + .invoke_scalar("LINESTRING (0 0, 1 1, 0 1, 0 0)") + .unwrap(); + tester.assert_scalar_result_equals(result, "MULTIPOINT EMPTY"); + + let result = tester + .invoke_scalar("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))") + .unwrap(); + tester.assert_scalar_result_equals(result, "LINESTRING(0 0,0 10,10 10,10 0,0 0)"); + + let result = tester + .invoke_scalar( + "MULTIPOLYGON (((0 0, 0 1, 1 1, 1 0, 0 0)), ((10 10, 10 11, 11 11, 11 10, 10 10)))", + ) + .unwrap(); + tester.assert_scalar_result_equals( + result, + "MULTILINESTRING((0 0,0 1,1 1,1 0,0 0),(10 10,10 11,11 11,11 10,10 10))", + ); + + let result = tester + .invoke_scalar("GEOMETRYCOLLECTION(POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0)), GEOMETRYCOLLECTION(LINESTRING(10 10, 10 20)))") + .unwrap(); + tester.assert_scalar_result_equals( + result, + "GEOMETRYCOLLECTION(MULTIPOINT((10 10),(10 20)), LINESTRING(0 0,0 1,1 1,1 0,0 0))", + ); + + let result = tester.invoke_scalar("GEOMETRYCOLLECTION EMPTY").unwrap(); + tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); + + let result = tester.invoke_scalar("POLYGON EMPTY").unwrap(); + tester.assert_scalar_result_equals(result, "MULTILINESTRING EMPTY"); + + let result = tester.invoke_scalar("MULTIPOLYGON EMPTY").unwrap(); + tester.assert_scalar_result_equals(result, "MULTILINESTRING EMPTY"); + + let result = tester + .invoke_scalar("LINESTRING Z (0 0 1, 1 1 2, 0 0 3)") + .unwrap(); + tester.assert_scalar_result_equals(result, "MULTIPOINT Z EMPTY"); + } + + #[rstest] + fn udf_invoke_item_crs( + #[values(WKB_GEOMETRY_ITEM_CRS.clone(), WKB_GEOGRAPHY_ITEM_CRS.clone())] + sedona_type: SedonaType, + ) { + let tester = ScalarUdfTester::new(st_boundary_udf().into(), vec![sedona_type.clone()]); + tester.assert_return_type(sedona_type); + + let result = tester.invoke_scalar("POINT (1 3)").unwrap(); + tester.assert_scalar_result_equals(result, "GEOMETRYCOLLECTION EMPTY"); + } +}