 Home Current PostGIS Blog Articles

Search

Browse

Email me  ## COGO: Finding centre and radius of a curve defined by three points (PostGIS)

Wednesday February 29 2012 at 06:25

KeywordsCOGO find centre radius center circle three 3 points
Summary

This article describes a pgPlSql function that, given three points, will return the centre and radius of the circle so defined. This works only for projected data.

Recently I had need to convert a PL/SQL Oracle Spatial function I created years ago called FindCircle to SQL Server 2008 for use in another project. That function was original work already released to the public domain as part of my free COGO package for Oracle. Here is that function
for SQL Server.

Note that I have a schema call cogo in which I create functions like this. You can use anything you like.

1. /** ----------------------------------------------------------------------------------------
2.   * @function   : FindCircle
3.   * @precis     : Function that determines if three points form a circle. If so a table containing
4.   *               centre and radius is returned. If not, a null table is returned.
5.   * @version    : 1.0
6.   * @param      : p_pt1        : First point in curve
7.   * @param      : p_pt2        : Second point in curve
8.   * @param      : p_pt3        : Third point in curve
9.   * @return     : geometry     : In which X,Y ordinates are the centre X, Y and the Z being the radius of found circle
10.   *                              or NULL if three points do not form a circle.
11.   * @history    : Simon Greener - Feb 2012 - Original coding.
12.   * @copyright  : Simon Greener @ 2012
14. **/
15. CREATE OR REPLACE FUNCTION Find_Circle(p_pt1 geometry, p_pt2 geometry, p_pt3 geometry)
16.   RETURNS geometry AS
17. \$BODY\$
18. DECLARE
19.    v_Centre geometry;
21.    v_CX     NUMERIC;
22.    v_CY     NUMERIC;
23.    v_dA     NUMERIC;
24.    v_dB     NUMERIC;
25.    v_dC     NUMERIC;
26.    v_dD     NUMERIC;
27.    v_dE     NUMERIC;
28.    v_dF     NUMERIC;
29.    v_dG     NUMERIC;
30. BEGIN
31.    IF ( p_pt1 IS NULL OR p_pt2 IS NULL OR p_pt3 IS NULL ) THEN
32.       RAISE EXCEPTION 'All supplied points must be not null.';
33.       RETURN NULL;
34.    END IF;
35.    IF ( ST_GeometryType(p_pt1) <> 'ST_Point' OR
36.         ST_GeometryType(p_pt1) <> 'ST_Point' OR
37.         ST_GeometryType(p_pt1) <> 'ST_Point' ) THEN
38.       RAISE EXCEPTION 'All supplied geometries must be points.';
39.       RETURN NULL;
40.    END IF;
41.    v_dA := ST_X(p_pt2) - ST_X(p_pt1);
42.    v_dB := ST_Y(p_pt2) - ST_Y(p_pt1);
43.    v_dC := ST_X(p_pt3) - ST_X(p_pt1);
44.    v_dD := ST_Y(p_pt3) - ST_Y(p_pt1);
45.    v_dE := v_dA * (ST_X(p_pt1) + ST_X(p_pt2)) + v_dB * (ST_Y(p_pt1) + ST_Y(p_pt2));
46.    v_dF := v_dC * (ST_X(p_pt1) + ST_X(p_pt3)) + v_dD * (ST_Y(p_pt1) + ST_Y(p_pt3));
47.    v_dG := 2.0  * (v_dA * (ST_Y(p_pt3) - ST_Y(p_pt2)) - v_dB * (ST_X(p_pt3) - ST_X(p_pt2)));
48.    -- If v_dG is zero then the three points are collinear and no finite-radius
49.    -- circle through them exists.
50.    IF ( v_dG = 0 ) THEN
51.       RETURN NULL;
52.    ELSE
53.       v_CX := (v_dD * v_dE - v_dB * v_dF) / v_dG;
54.       v_CY := (v_dA * v_dF - v_dC * v_dE) / v_dG;
55.       v_Radius := SQRT(POWER(ST_X(p_pt1) - v_CX,2) + POWER(ST_Y(p_pt1) - v_CY,2) );
56.    END IF;
58. END;
59. \$BODY\$
60.   LANGUAGE plpgsql VOLATILE STRICT;
p.The function can be used as follows:
1. SELECT ST_X(c.circle) AS CX, ST_Y(c.circle) AS CY, ST_Z(c.circle) AS radius
2.   FROM (SELECT Find_Circle(ST_MakePoint(0.0,0.0),
3.                            ST_MakePoint(10.0,0.0),
4.                            ST_MakePoint(10.0,10.0)) AS circle
5.         ) AS c;

Resulting in ….

5 5 7.07106781186548

Or can be used as follows…

1. SELECT ST_AsEWKT(Find_Circle(ST_SetSrid(ST_MakePoint(525000, 5202000),32615),
2.                              ST_SetSrid(ST_MakePoint(525500, 5202000),32615),
3.                              ST_SetSrid(ST_MakePoint(525500, 5202500),32615))) AS circle;

Resulting in ….

circle
SRID=32615;POINT (525250 5202250 353.553390593274)

I hope this is of help to someone.

Simon          ### Comment

 Name Remember E-mail Message Textile Help 