Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

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

Wednesday February 29 2012 at 17: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
  13.   *               Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
  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;
  20.    v_radius NUMERIC;
  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;
  57.    RETURN ST_SetSRID(ST_MakePoint(v_CX, v_CY, v_radius),ST_Srid(p_pt1));
  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 ….

cx cy radius
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

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment