voronoi

syntax

voronoi(point_data_item, extent_unit)

definition

voronoi(point_data_item, extent_unit) results in a Polygon Data item with, for each point in the point_data_item Argument, the Voronoi cell (also called Thiessen polygon) of that point: the set of locations that are at least as close to it as to any other point. The result has the same Domain unit and the same Values unit as the point_data_item.

The cells of the points on the convex hull are unbounded, so the extent_unit argument supplies the rectangle they are clipped to: its Range is used as the clipping rectangle.

description

If every point lies inside extent_unit, the resulting cells partition that rectangle exactly: they do not overlap and they leave no gaps, so the areas add up to the area of the rectangle. That makes voronoi the usual way to turn a set of measurement locations into an area partition, for instance to assign each part of a region to its nearest facility.

The Voronoi diagram is the dual of the Delaunay triangulation, which triangualize returns: the cells of two points share a boundary exactly when those points are connected by a Delaunay edge.

applies to

Attribute point_data_item with a point Value type and Single Composition
Unit extent_unit with the same point Value type and a defined Range

conditions

The extent_unit must have a proper range: a range whose lower bound is strictly below its upper bound in both dimensions. Without it the unbounded cells cannot be closed and the operator reports an error.
The Values unit of point_data_item and the extent_unit must be compatible.
A point outside the extent_unit range gets an empty cell.
[[null Undefined]] points get an empty cell.
Of a group of coinciding points, one keeps the cell and the others get an empty cell. Their true cells would coincide, so returning them all would produce overlapping polygons that silently double-count area.
Collinear point sets are allowed and give the expected slabs.

since version

20.12

performance

O(m log m), dominated by the Delaunay triangulation that the cells are derived from - see triangualize for that part. Each cell is then built by clipping the extent rectangle with the perpendicular bisector between the point and each of its Delaunay neighbours, which is O(deg²) per point; the average degree of a Delaunay vertex is below 6, so in practice this stage is linear.

Half-plane clipping is used rather than joining circumcentres because it needs nothing but the neighbour relation: hull points need no unbounded-ray special case, and a collinear point set - whose triangulation has no triangles at all, and therefore no circumcentres - still yields the correct cells.

Memory is O(m): about 200 bytes per point for the triangulation, plus the neighbour lists, plus the resulting polygons.

example

unit<ipoint> extent := range(ipoint, point_yx(0i, 0i), point_yx(1000i, 1000i));

unit<uint32> station: nrofrows = 4
{
    attribute<extent>  geometry     : [yx(250,250), yx(250,750), yx(750,250), yx(750,750)];
    attribute<extent>  cell   (poly) := voronoi(geometry, extent);
    attribute<float64> cell_area     := area(dpolygon(cell), float64);
}

geometry	cell_area
yx(250,250)	250000
yx(250,750)	250000
yx(750,250)	250000
yx(750,750)	250000

domain station, nr of rows = 4

The four points sit at the centres of the quadrants of the extent, so each cell is one 500 x 500 quadrant and the four areas add up to the 1000 x 1000 extent.