diff --git a/src/pg/sql/09_voronoi.sql b/src/pg/sql/09_voronoi.sql index 51d22d4..20be4c9 100644 --- a/src/pg/sql/09_voronoi.sql +++ b/src/pg/sql/09_voronoi.sql @@ -1,3 +1,69 @@ +/** ---------------------------------------------------------------------------------------- + * @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 + * @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. + * @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(p_pt1 geometry, p_pt2 geometry, 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_pt1) <> 'ST_Point' OR + ST_GeometryType(p_pt1) <> '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; + +-- ============================================================================================= +-- +-- CDB_Voronoi +-- +-- ============================================================================================= CREATE OR REPLACE FUNCTION CDB_voronoi( IN geomin geometry[], IN buffer numeric DEFAULT 0.5, @@ -143,63 +209,4 @@ BEGIN 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 - * @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. - * @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(p_pt1 geometry, p_pt2 geometry, 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_pt1) <> 'ST_Point' OR - ST_GeometryType(p_pt1) <> '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; +