diff --git a/vortex-geo/Cargo.toml b/vortex-geo/Cargo.toml index 82d5eb520d9..0b4262fcf90 100644 --- a/vortex-geo/Cargo.toml +++ b/vortex-geo/Cargo.toml @@ -74,5 +74,9 @@ harness = false name = "collect" harness = false +[[bench]] +name = "convex_hull" +harness = false + [lints] workspace = true diff --git a/vortex-geo/benches/convex_hull.rs b/vortex-geo/benches/convex_hull.rs new file mode 100644 index 00000000000..0c3d5fdc684 --- /dev/null +++ b/vortex-geo/benches/convex_hull.rs @@ -0,0 +1,96 @@ +// SPDX-License-Identifier: Apache-2.0 +// SPDX-FileCopyrightText: Copyright the Vortex contributors + +//! Microbenchmarks for native `ST_ConvexHull` over `MultiPoint` rows. +//! +//! The cases separate ordinary small hulls, larger point sets, and strict null propagation. They +//! execute the result to its canonical polygon representation. +//! +//! Run with `cargo bench -p vortex-geo --bench convex_hull`. + +#![expect(clippy::unwrap_used)] + +use std::f64::consts::TAU; +use std::sync::LazyLock; + +use divan::Bencher; +use divan::counter::ItemsCount; +use mimalloc::MiMalloc; +use vortex_array::ArrayRef; +use vortex_array::Canonical; +use vortex_array::ExecutionCtx; +use vortex_array::IntoArray; +use vortex_array::VortexSessionExecute; +use vortex_array::arrays::MaskedArray; +use vortex_array::validity::Validity; +use vortex_geo::scalar_fn::convex_hull::GeoConvexHull; +use vortex_geo::test_harness::geo_session; +use vortex_geo::test_harness::multipoint_column; +use vortex_session::VortexSession; + +// Scalar function execution allocates its output inside the timed region, so use the vendored +// allocator instead of measuring glibc differences between CodSpeed runner images. +#[global_allocator] +static GLOBAL: MiMalloc = MiMalloc; + +static SESSION: LazyLock = LazyLock::new(geo_session); + +const ROWS: usize = 512; + +fn main() { + divan::main(); +} + +fn multipoints(points_per_row: usize) -> ArrayRef { + multipoint_column( + (0..ROWS) + .map(|row| { + (0..points_per_row) + .map(|point| { + let angle = TAU * point as f64 / points_per_row as f64; + let radius = 10.0 + ((row + point) % 7) as f64; + (radius * angle.cos(), radius * angle.sin()) + }) + .collect() + }) + .collect(), + ) + .unwrap() +} + +fn hulls(input: &ArrayRef, ctx: &mut ExecutionCtx) -> ArrayRef { + GeoConvexHull::try_new_array(input.clone()) + .unwrap() + .into_array() + .execute::(ctx) + .unwrap() + .into_array() +} + +fn bench_hulls(bencher: Bencher, input: ArrayRef) { + let mut ctx = SESSION.create_execution_ctx(); + bencher + .counter(ItemsCount::new(ROWS)) + .bench_local(|| hulls(&input, &mut ctx)); +} + +#[divan::bench] +fn eight_points(bencher: Bencher) { + bench_hulls(bencher, multipoints(8)); +} + +#[divan::bench] +fn sixty_four_points(bencher: Bencher) { + bench_hulls(bencher, multipoints(64)); +} + +#[divan::bench] +fn nullable_eight_points(bencher: Bencher) { + let input = MaskedArray::try_new( + multipoints(8), + Validity::from_iter((0..ROWS).map(|row| !row.is_multiple_of(8))), + ) + .unwrap() + .into_array(); + bench_hulls(bencher, input); +} diff --git a/vortex-geo/src/extension/polygon.rs b/vortex-geo/src/extension/polygon.rs index 9a74c3ce7b3..9b825dfb1bf 100644 --- a/vortex-geo/src/extension/polygon.rs +++ b/vortex-geo/src/extension/polygon.rs @@ -13,9 +13,11 @@ use arrow_schema::Field; use arrow_schema::extension::ExtensionType; use geo_traits::to_geo::ToGeoGeometry; use geo_types::Geometry; +use geoarrow::array::GeoArrowArray; use geoarrow::array::GeoArrowArrayAccessor; use geoarrow::array::IntoArrow; use geoarrow::array::PolygonArray; +use geoarrow::array::PolygonBuilder; use geoarrow::datatypes::CoordType; use geoarrow::datatypes::PolygonType; use prost::Message; @@ -24,6 +26,7 @@ use vortex_array::ExecutionCtx; use vortex_array::IntoArray; use vortex_array::arrays::ExtensionArray; use vortex_array::arrays::extension::ExtensionArrayExt; +use vortex_array::builtins::ArrayBuiltins; use vortex_array::dtype::DType; use vortex_array::dtype::Nullability; use vortex_array::dtype::extension::ExtDType; @@ -114,6 +117,25 @@ fn polygon_type(geo_metadata: &GeoMetadata, dimension: Dimension) -> PolygonType PolygonType::new(dimension.into(), geoarrow_metadata(geo_metadata)) } +/// Build a native 2-D [`Polygon`] array from row-oriented `geo_types` polygons. +pub(crate) fn build_polygon_array( + polygons: &[Option>], + metadata: GeoMetadata, + nullability: Nullability, +) -> VortexResult { + let polygons = + PolygonBuilder::from_nullable_polygons(polygons, polygon_type(&metadata, Dimension::Xy)) + .finish(); + let storage_dtype = polygon_storage_dtype(Dimension::Xy, nullability); + let storage = ArrayRef::from_arrow( + polygons.to_array_ref().as_ref(), + nullability == Nullability::Nullable, + )? + .cast(storage_dtype.clone())?; + let ext_dtype = ExtDType::::try_new(metadata, storage_dtype)?; + Ok(ExtensionArray::try_new(ext_dtype.erased(), storage)?.into_array()) +} + /// Decode `Polygon` storage (`List>`) to `geo_types` polygons, for the geo scalar /// functions. CRS does not affect planar geometry ops, so default metadata is used. pub(crate) fn polygon_geometries( diff --git a/vortex-geo/src/lib.rs b/vortex-geo/src/lib.rs index 0a9afecb85e..901b8223947 100644 --- a/vortex-geo/src/lib.rs +++ b/vortex-geo/src/lib.rs @@ -24,6 +24,7 @@ use crate::prune::GeoIntersectsPrune; use crate::scalar_fn::area::GeoArea; use crate::scalar_fn::collect::GeoCollect; use crate::scalar_fn::contains::GeoContains; +use crate::scalar_fn::convex_hull::GeoConvexHull; use crate::scalar_fn::distance::GeoDistance; use crate::scalar_fn::envelope::GeoEnvelope; use crate::scalar_fn::intersects::GeoIntersects; @@ -70,6 +71,7 @@ pub fn initialize(session: &VortexSession) { // Register the geometry scalar functions. session.scalar_fns().register(GeoArea); session.scalar_fns().register(GeoCollect); + session.scalar_fns().register(GeoConvexHull); session.scalar_fns().register(GeoEnvelope); session.scalar_fns().register(GeoContains); session.scalar_fns().register(GeoDistance); diff --git a/vortex-geo/src/scalar_fn/convex_hull.rs b/vortex-geo/src/scalar_fn/convex_hull.rs new file mode 100644 index 00000000000..b87e147f434 --- /dev/null +++ b/vortex-geo/src/scalar_fn/convex_hull.rs @@ -0,0 +1,353 @@ +// SPDX-License-Identifier: Apache-2.0 +// SPDX-FileCopyrightText: Copyright the Vortex contributors + +//! `ST_ConvexHull`: the planar convex hull of each native `MultiPoint`. + +use geo::ConvexHull; +use vortex_array::ArrayRef; +use vortex_array::ExecutionCtx; +use vortex_array::IntoArray; +use vortex_array::arrays::ConstantArray; +use vortex_array::arrays::ScalarFnArray; +use vortex_array::dtype::DType; +use vortex_array::dtype::extension::ExtDType; +use vortex_array::dtype::extension::ExtDTypeRef; +use vortex_array::expr::Expression; +use vortex_array::expr::union_child_validities; +use vortex_array::scalar_fn::Arity; +use vortex_array::scalar_fn::ChildName; +use vortex_array::scalar_fn::EmptyOptions; +use vortex_array::scalar_fn::ExecutionArgs; +use vortex_array::scalar_fn::ScalarFnId; +use vortex_array::scalar_fn::ScalarFnVTable; +use vortex_array::scalar_fn::TypedScalarFnInstance; +use vortex_error::VortexResult; +use vortex_error::vortex_bail; +use vortex_error::vortex_ensure; +use vortex_mask::AllOr; +use vortex_mask::Mask; +use vortex_session::VortexSession; +use vortex_session::registry::CachedId; + +use crate::extension::MultiPoint; +use crate::extension::Polygon; +use crate::extension::build_polygon_array; +use crate::extension::coordinate::Dimension; +use crate::extension::geometries; +use crate::extension::polygon_storage_dtype; +use crate::extension::single_geometry; +use crate::scalar_fn::execute::Execution; +use crate::scalar_fn::execute::Operand; +use crate::scalar_fn::execute::dispatch_unary; + +/// Resolve the strict native `MultiPoint -> Polygon` overload. +fn convex_hull_dtype(dtypes: &[DType]) -> VortexResult { + vortex_ensure!( + dtypes.len() == 1, + "geo: convex_hull requires exactly one MultiPoint operand, got {}", + dtypes.len() + ); + let Some(input) = dtypes[0].as_extension_opt() else { + vortex_bail!( + "geo: convex_hull operand {} is not a native MultiPoint", + dtypes[0] + ); + }; + vortex_ensure!( + input.is::(), + "geo: convex_hull operand {} is not a native MultiPoint", + dtypes[0] + ); + + Ok(ExtDType::::try_new( + input.metadata::().clone(), + polygon_storage_dtype(Dimension::Xy, dtypes[0].nullability()), + )? + .erased()) +} + +/// Compute hulls for the valid rows and scatter them into a full-length native polygon array. +fn convex_hull_array( + array: ArrayRef, + valid: &Mask, + output_dtype: &ExtDTypeRef, + ctx: &mut ExecutionCtx, +) -> VortexResult { + let decoded = geometries(&array.filter(valid.clone())?, ctx)?; + let hulls = decoded.iter().map(ConvexHull::convex_hull); + let polygons = match valid.indices() { + AllOr::All => hulls.map(Some).collect(), + AllOr::None => vec![None; array.len()], + AllOr::Some(rows) => { + let mut polygons = vec![None; array.len()]; + for (&row, hull) in rows.iter().zip(hulls) { + polygons[row] = Some(hull); + } + polygons + } + }; + build_polygon_array( + &polygons, + output_dtype.metadata::().clone(), + output_dtype.nullability(), + ) +} + +/// Execute convex hull after shared unary shape and null dispatch. +fn execute_convex_hull( + execution: Execution<1>, + output_dtype: &ExtDTypeRef, + ctx: &mut ExecutionCtx, +) -> VortexResult { + match execution.operands { + [Operand::Constant(scalar)] => { + let hull = single_geometry(&scalar, ctx)?.convex_hull(); + let output = build_polygon_array( + &[Some(hull)], + output_dtype.metadata::().clone(), + output_dtype.nullability(), + )?; + Ok(ConstantArray::new(output.execute_scalar(0, ctx)?, execution.len).into_array()) + } + [Operand::Column(array)] => convex_hull_array(array, &execution.valid, output_dtype, ctx), + } +} + +/// Compute the two-dimensional convex hull of each native `MultiPoint` as a native `Polygon`. +/// Empty, single-point, and collinear inputs remain typed polygons with degenerate exterior rings. +#[derive(Debug, Clone, Default, PartialEq, Eq, Hash)] +pub struct GeoConvexHull; + +impl GeoConvexHull { + /// A lazy `ScalarFnArray` computing a polygon hull for each native `MultiPoint` row. + pub fn try_new_array(array: ArrayRef) -> VortexResult { + ScalarFnArray::try_new( + TypedScalarFnInstance::new(GeoConvexHull, EmptyOptions).erased(), + vec![array], + ) + } +} + +impl ScalarFnVTable for GeoConvexHull { + type Options = EmptyOptions; + + fn id(&self) -> ScalarFnId { + static ID: CachedId = CachedId::new("vortex.geo.convex_hull"); + *ID + } + + fn serialize(&self, _: &Self::Options) -> VortexResult>> { + Ok(Some(vec![])) + } + + fn deserialize(&self, _: &[u8], _: &VortexSession) -> VortexResult { + Ok(EmptyOptions) + } + + fn arity(&self, _: &Self::Options) -> Arity { + Arity::Exact(1) + } + + fn child_name(&self, _: &Self::Options, child_idx: usize) -> ChildName { + match child_idx { + 0 => ChildName::from("multipoint"), + _ => unreachable!("convex_hull has exactly one child"), + } + } + + fn return_dtype(&self, _: &Self::Options, dtypes: &[DType]) -> VortexResult { + Ok(DType::Extension(convex_hull_dtype(dtypes)?)) + } + + fn execute( + &self, + _: &Self::Options, + args: &dyn ExecutionArgs, + ctx: &mut ExecutionCtx, + ) -> VortexResult { + let input = args.get(0)?; + let output_dtype = convex_hull_dtype(std::slice::from_ref(input.dtype()))?; + dispatch_unary( + &input, + DType::Extension(output_dtype.clone()), + |execution, ctx| execute_convex_hull(execution, &output_dtype, ctx), + ctx, + ) + } + + fn validity( + &self, + _: &Self::Options, + expression: &Expression, + ) -> VortexResult> { + union_child_validities(expression) + } + + fn is_strict(&self, _: &Self::Options) -> bool { + true + } + + fn is_fallible(&self, _: &Self::Options) -> bool { + false + } +} + +#[cfg(test)] +mod tests { + use rstest::rstest; + use vortex_array::ArrayRef; + use vortex_array::Columnar; + use vortex_array::IntoArray; + use vortex_array::VortexSessionExecute; + use vortex_array::arrays::ConstantArray; + use vortex_array::arrays::ListArray; + use vortex_array::arrays::MaskedArray; + use vortex_array::arrays::PrimitiveArray; + use vortex_array::assert_arrays_eq; + use vortex_array::scalar_fn::EmptyOptions; + use vortex_array::scalar_fn::ScalarFnVTable; + use vortex_array::validity::Validity; + use vortex_error::VortexResult; + use vortex_error::vortex_err; + + use super::GeoConvexHull; + use crate::scalar_fn::area::GeoArea; + use crate::scalar_fn::collect::GeoCollect; + use crate::test_harness::multipoint_column; + use crate::test_harness::point_column; + use crate::test_harness::polygon_column; + + #[test] + fn computes_polygon_hulls() -> VortexResult<()> { + let input = multipoint_column(vec![vec![ + (0.0, 0.0), + (2.0, 0.0), + (2.0, 2.0), + (0.0, 2.0), + (1.0, 1.0), + ]])?; + let expected = polygon_column(vec![vec![vec![ + (2.0, 0.0), + (2.0, 2.0), + (0.0, 2.0), + (0.0, 0.0), + (2.0, 0.0), + ]]])?; + let result = GeoConvexHull::try_new_array(input)?.into_array(); + let mut ctx = vortex_array::array_session().create_execution_ctx(); + + assert_arrays_eq!(result, expected, &mut ctx); + Ok(()) + } + + #[rstest] + #[case::empty(vec![], vec![])] + #[case::one_point( + vec![(1.0, 2.0)], + vec![vec![(1.0, 2.0), (1.0, 2.0)]] + )] + #[case::collinear( + vec![(0.0, 0.0), (1.0, 1.0), (2.0, 2.0)], + vec![vec![(0.0, 0.0), (2.0, 2.0), (0.0, 0.0)]] + )] + fn degenerate_hulls_remain_polygons( + #[case] points: Vec<(f64, f64)>, + #[case] expected_rings: Vec>, + ) -> VortexResult<()> { + let input = multipoint_column(vec![points])?; + let expected = polygon_column(vec![expected_rings])?; + let result = GeoConvexHull::try_new_array(input)?.into_array(); + let mut ctx = vortex_array::array_session().create_execution_ctx(); + + assert_arrays_eq!(result, expected, &mut ctx); + Ok(()) + } + + #[test] + fn propagates_nulls() -> VortexResult<()> { + let input = MaskedArray::try_new( + multipoint_column(vec![ + vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)], + vec![(2.0, 2.0)], + ])?, + Validity::from_iter([true, false]), + )? + .into_array(); + let expected = MaskedArray::try_new( + polygon_column(vec![ + vec![vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (0.0, 0.0)]], + vec![vec![(2.0, 2.0), (2.0, 2.0)]], + ])?, + Validity::from_iter([true, false]), + )? + .into_array(); + let result = GeoConvexHull::try_new_array(input)?.into_array(); + let mut ctx = vortex_array::array_session().create_execution_ctx(); + + assert_arrays_eq!(result, expected, &mut ctx); + Ok(()) + } + + #[test] + fn constant_remains_constant() -> VortexResult<()> { + let mut ctx = vortex_array::array_session().create_execution_ctx(); + let scalar = multipoint_column(vec![vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]])? + .execute_scalar(0, &mut ctx)?; + let input = ConstantArray::new(scalar, 3).into_array(); + + let result = GeoConvexHull::try_new_array(input)?.into_array(); + let Columnar::Constant(constant) = result.clone().execute::(&mut ctx)? else { + return Err(vortex_err!( + "convex_hull of a constant should remain constant" + )); + }; + assert_eq!(constant.len(), 3); + let hull = vec![vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (0.0, 0.0)]]; + let expected = polygon_column(vec![hull.clone(), hull.clone(), hull])?; + assert_arrays_eq!(result, expected, &mut ctx); + Ok(()) + } + + #[test] + fn collect_hull_area_pipeline() -> VortexResult<()> { + let points = point_column( + vec![0.0, 2.0, 2.0, 0.0, 1.0, 0.0, 1.0, 2.0], + vec![0.0, 0.0, 2.0, 2.0, 1.0, 0.0, 1.0, 2.0], + )?; + let point_lists = ListArray::try_new( + points, + PrimitiveArray::from_iter([0_u32, 5, 8]).into_array(), + Validity::NonNullable, + )? + .into_array(); + + let collected = GeoCollect::try_new_array(point_lists)?.into_array(); + let hulls = GeoConvexHull::try_new_array(collected)?.into_array(); + let areas = GeoArea::try_new_array(hulls)?.into_array(); + let expected = PrimitiveArray::from_iter([4.0_f64, 0.0]).into_array(); + let mut ctx = vortex_array::array_session().create_execution_ctx(); + + assert_arrays_eq!(areas, expected, &mut ctx); + Ok(()) + } + + #[rstest] + #[case::none(0)] + #[case::two(2)] + fn rejects_wrong_arity(#[case] arity: usize) -> VortexResult<()> { + let dtype = multipoint_column(vec![vec![]])?.dtype().clone(); + assert!( + GeoConvexHull + .return_dtype(&EmptyOptions, &vec![dtype; arity]) + .is_err() + ); + Ok(()) + } + + #[test] + fn rejects_non_multipoint_input() -> VortexResult<()> { + let input: ArrayRef = point_column(vec![0.0], vec![0.0])?; + assert!(GeoConvexHull::try_new_array(input).is_err()); + Ok(()) + } +} diff --git a/vortex-geo/src/scalar_fn/mod.rs b/vortex-geo/src/scalar_fn/mod.rs index ef5447fad91..99fe5d28528 100644 --- a/vortex-geo/src/scalar_fn/mod.rs +++ b/vortex-geo/src/scalar_fn/mod.rs @@ -6,6 +6,7 @@ pub mod area; pub mod collect; pub mod contains; +pub mod convex_hull; pub mod distance; pub mod envelope; mod execute;