216 lines
6.5 KiB
PL/PgSQL
216 lines
6.5 KiB
PL/PgSQL
-- =============================================================================================
|
|
--
|
|
-- CDB_Voronoi
|
|
--
|
|
-- =============================================================================================
|
|
CREATE OR REPLACE FUNCTION CDB_voronoi(
|
|
IN geomin geometry[],
|
|
IN buffer numeric DEFAULT 0.5,
|
|
IN tolerance numeric DEFAULT 1e-9
|
|
)
|
|
RETURNS geometry AS $$
|
|
DECLARE
|
|
geomout geometry;
|
|
BEGIN
|
|
WITH
|
|
convexhull_1 as (
|
|
SELECT
|
|
ST_ConvexHull(ST_Collect(geomin)) as g,
|
|
buffer * |/ st_area(ST_ConvexHull(ST_Collect(geomin))) as r
|
|
),
|
|
clipper as(
|
|
SELECT
|
|
st_buffer(ST_MinimumBoundingCircle(a.g), buffer*a.r) as g
|
|
FROM convexhull_1 a
|
|
),
|
|
env0 as (
|
|
SELECT
|
|
(st_dumppoints(st_expand(a.g, buffer*a.r))).geom as e
|
|
FROM convexhull_1 a
|
|
),
|
|
env as (
|
|
SELECT
|
|
array_agg(env0.e) as e
|
|
FROM env0
|
|
),
|
|
sample AS (
|
|
SELECT
|
|
ST_Collect(geomin || env.e) as geom
|
|
FROM env
|
|
),
|
|
convexhull as (
|
|
SELECT
|
|
ST_ConvexHull(ST_Collect(geomin)) as cg
|
|
),
|
|
tin as (
|
|
SELECT
|
|
ST_Dump(ST_DelaunayTriangles(geom, tolerance, 0)) as gd
|
|
FROM
|
|
sample
|
|
),
|
|
tin_polygons as (
|
|
SELECT
|
|
(gd).Path as id,
|
|
(gd).Geom as pg,
|
|
ST_Centroid(ST_MinimumBoundingCircle((gd).Geom, 180)) as ct
|
|
FROM tin
|
|
),
|
|
tin_lines as (
|
|
SELECT
|
|
id,
|
|
ST_ExteriorRing(pg) as lg
|
|
FROM tin_polygons
|
|
),
|
|
tin_nodes as (
|
|
SELECT
|
|
id,
|
|
ST_PointN(lg,1) p1,
|
|
ST_PointN(lg,2) p2,
|
|
ST_PointN(lg,3) p3
|
|
FROM tin_lines
|
|
),
|
|
tin_edges AS (
|
|
SELECT
|
|
p.id,
|
|
UNNEST(ARRAY[
|
|
ST_MakeLine(n.p1,n.p2) ,
|
|
ST_MakeLine(n.p2,n.p3) ,
|
|
ST_MakeLine(n.p3,n.p1)]) as Edge,
|
|
ST_Force2D(cdb_crankshaft._Find_Circle(n.p1,n.p2,n.p3)) as ct,
|
|
CASE WHEN st_distance(p.ct, ST_ExteriorRing(p.pg)) < tolerance THEN
|
|
TRUE
|
|
ELSE FALSE END AS ctx,
|
|
p.pg,
|
|
ST_within(p.ct, convexhull.cg) as ctin
|
|
FROM
|
|
tin_polygons p,
|
|
tin_nodes n,
|
|
convexhull
|
|
WHERE p.id = n.id
|
|
),
|
|
voro_nodes as (
|
|
SELECT
|
|
CASE WHEN x.ctx = TRUE THEN
|
|
ST_Centroid(x.edge)
|
|
ELSE
|
|
x.ct
|
|
END as xct,
|
|
CASE WHEN y.id is null THEN
|
|
CASE WHEN x.ctin = TRUE THEN
|
|
ST_SetSRID(ST_MakePoint(
|
|
ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * (1+buffer)),
|
|
ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * (1+buffer))
|
|
), ST_SRID(x.ct))
|
|
END
|
|
ELSE
|
|
y.ct
|
|
END as yct
|
|
FROM
|
|
tin_edges x
|
|
LEFT OUTER JOIN
|
|
tin_edges y
|
|
ON x.id <> y.id AND ST_Equals(x.edge, y.edge)
|
|
),
|
|
voro_edges as(
|
|
SELECT
|
|
ST_LineMerge(ST_Collect(ST_MakeLine(xct, yct))) as v
|
|
FROM
|
|
voro_nodes
|
|
),
|
|
voro_cells as(
|
|
SELECT
|
|
ST_Polygonize(
|
|
ST_Node(
|
|
ST_LineMerge(
|
|
ST_Union(v, ST_ExteriorRing(
|
|
ST_Convexhull(v)
|
|
)
|
|
)
|
|
)
|
|
)
|
|
) as g
|
|
FROM
|
|
voro_edges
|
|
),
|
|
voro_set as(
|
|
SELECT
|
|
(st_dump(v.g)).geom as g
|
|
FROM voro_cells v
|
|
)
|
|
SELECT
|
|
st_collect(ST_intersection(c.g, v.g))
|
|
-- ST_intersection(c.g, v.g)
|
|
INTO geomout
|
|
FROM
|
|
voro_set v,
|
|
clipper c;
|
|
RETURN geomout;
|
|
END;
|
|
$$ language plpgsql IMMUTABLE;
|
|
|
|
/** ----------------------------------------------------------------------------------------
|
|
* @function : FindCircle
|
|
* @precis : Function that determines if three points form a circle. If so a table containing
|
|
* centre and radius is returned. If not, a null table is returned.
|
|
* @version : 1.0.1
|
|
* @param : p_pt1 : First point in curve
|
|
* @param : p_pt2 : Second point in curve
|
|
* @param : p_pt3 : Third point in curve
|
|
* @return : geometry : In which X,Y ordinates are the centre X, Y and the Z being the radius of found circle
|
|
* or NULL if three points do not form a circle.
|
|
* @history : Simon Greener - Feb 2012 - Original coding.
|
|
* Rafa de la Torre - Aug 2016 - Small fix for type checking
|
|
* @copyright : Simon Greener @ 2012
|
|
* Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
|
|
**/
|
|
CREATE OR REPLACE FUNCTION _Find_Circle(
|
|
IN p_pt1 geometry,
|
|
IN p_pt2 geometry,
|
|
IN p_pt3 geometry)
|
|
RETURNS geometry AS
|
|
$BODY$
|
|
DECLARE
|
|
v_Centre geometry;
|
|
v_radius NUMERIC;
|
|
v_CX NUMERIC;
|
|
v_CY NUMERIC;
|
|
v_dA NUMERIC;
|
|
v_dB NUMERIC;
|
|
v_dC NUMERIC;
|
|
v_dD NUMERIC;
|
|
v_dE NUMERIC;
|
|
v_dF NUMERIC;
|
|
v_dG NUMERIC;
|
|
BEGIN
|
|
IF ( p_pt1 IS NULL OR p_pt2 IS NULL OR p_pt3 IS NULL ) THEN
|
|
RAISE EXCEPTION 'All supplied points must be not null.';
|
|
RETURN NULL;
|
|
END IF;
|
|
IF ( ST_GeometryType(p_pt1) <> 'ST_Point' OR
|
|
ST_GeometryType(p_pt2) <> 'ST_Point' OR
|
|
ST_GeometryType(p_pt3) <> 'ST_Point' ) THEN
|
|
RAISE EXCEPTION 'All supplied geometries must be points.';
|
|
RETURN NULL;
|
|
END IF;
|
|
v_dA := ST_X(p_pt2) - ST_X(p_pt1);
|
|
v_dB := ST_Y(p_pt2) - ST_Y(p_pt1);
|
|
v_dC := ST_X(p_pt3) - ST_X(p_pt1);
|
|
v_dD := ST_Y(p_pt3) - ST_Y(p_pt1);
|
|
v_dE := v_dA * (ST_X(p_pt1) + ST_X(p_pt2)) + v_dB * (ST_Y(p_pt1) + ST_Y(p_pt2));
|
|
v_dF := v_dC * (ST_X(p_pt1) + ST_X(p_pt3)) + v_dD * (ST_Y(p_pt1) + ST_Y(p_pt3));
|
|
v_dG := 2.0 * (v_dA * (ST_Y(p_pt3) - ST_Y(p_pt2)) - v_dB * (ST_X(p_pt3) - ST_X(p_pt2)));
|
|
-- If v_dG is zero then the three points are collinear and no finite-radius
|
|
-- circle through them exists.
|
|
IF ( v_dG = 0 ) THEN
|
|
RETURN NULL;
|
|
ELSE
|
|
v_CX := (v_dD * v_dE - v_dB * v_dF) / v_dG;
|
|
v_CY := (v_dA * v_dF - v_dC * v_dE) / v_dG;
|
|
v_Radius := SQRT(POWER(ST_X(p_pt1) - v_CX,2) + POWER(ST_Y(p_pt1) - v_CY,2) );
|
|
END IF;
|
|
RETURN ST_SetSRID(ST_MakePoint(v_CX, v_CY, v_radius),ST_Srid(p_pt1));
|
|
END;
|
|
$BODY$
|
|
LANGUAGE plpgsql VOLATILE STRICT;
|
|
|