From d60efbc12c99a82f9b5e090d25115f6952c4fcce Mon Sep 17 00:00:00 2001
From: Regina Obe <lr@pcorp.us>
Date: Wed, 4 Oct 2017 16:19:29 +0000
Subject: [PATCH] =?utf8?q?ST=5FAngle=20function=20from=20R=C3=A9mi=20Cura?=
 =?utf8?q?=20Closes=20https://github.com/postgis/postgis/pull/97=20Closes?=
 =?utf8?q?=20#3876?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

git-svn-id: http://svn.osgeo.org/postgis/trunk@15884 b70326c6-7e19-0410-871a-916f4a2858ee
---
 NEWS                             |   2 +-
 doc/reference_measure.xml        |  66 ++++++++++++++++
 postgis/lwgeom_functions_basic.c | 126 +++++++++++++++++++++++++++++++
 postgis/postgis.sql.in           |  19 +++++
 regress/lwgeom_regress.sql       |  34 +++++++++
 regress/lwgeom_regress_expected  |  11 +++
 6 files changed, 257 insertions(+), 1 deletion(-)

diff --git a/NEWS b/NEWS
index 6b7a8335f..56ec4bfe5 100644
--- a/NEWS
+++ b/NEWS
@@ -1,7 +1,7 @@
 PostGIS 2.5.0
 2018/xx/xx
 * New Features *
-
+  - #3876, ST_Angle function (Rémi Cura)
 
 PostGIS 2.4.0
 2017/09/30
diff --git a/doc/reference_measure.xml b/doc/reference_measure.xml
index bf855a158..c376f02ac 100644
--- a/doc/reference_measure.xml
+++ b/doc/reference_measure.xml
@@ -819,6 +819,72 @@ SELECT degrees(ST_Azimuth(ST_Point(25, 45), ST_Point(75, 100))) AS degA_B,
 			<para><xref linkend="ST_Point" />, <xref linkend="ST_Translate" />, <xref linkend="ST_Project" />, <ulink url="http://www.postgresql.org/docs/current/interactive/functions-math.html">PostgreSQL Math Functions</ulink></para>
 		</refsection>
 
+  </refentry>
+  
+  	<refentry id="ST_Angle">
+		<refnamediv>
+		  <refname>ST_Angle</refname>
+
+		  <refpurpose>Returns the angle between 3 points, or between 2 vectors (4 points or 2 lines).</refpurpose>
+		</refnamediv>
+		<refsynopsisdiv>
+		  <funcsynopsis>
+			<funcprototype>
+			  <funcdef>float <function>ST_Angle</function></funcdef>
+			  <paramdef><type>geometry </type><parameter>point1</parameter></paramdef>
+			  <paramdef><type>geometry </type><parameter>point2</parameter></paramdef>
+			  <paramdef><type>geometry </type><parameter>point3</parameter></paramdef>
+			  <paramdef choice="opt"><type>geometry </type><parameter>point4</parameter></paramdef>
+			</funcprototype> 
+			<funcprototype>
+			  <funcdef>float <function>ST_Angle</function></funcdef>
+			  <paramdef><type>geometry </type><parameter>line1</parameter></paramdef>
+			  <paramdef><type>geometry </type><parameter>line2</parameter></paramdef>
+			</funcprototype>
+		  </funcsynopsis>
+		</refsynopsisdiv>
+		<refsection>
+			<title>Description</title>
+
+			<para> For 3 points, computes the angle measured clockwise of P1P2P3.
+			If input are 2 lines, get first and last point of the lines as 4 points.
+			For 4 points,compute the angle measured clockwise of P1P2,P3P4.
+			Results are always positive, between 0 and 2*Pi radians.
+			
+			Uses azimuth of pairs or points.
+			</para>
+			<para>ST_Angle(P1,P2,P3) = ST_Angle(P2,P1,P2,P3)</para> 
+			<para>Result is in radian and can be converted to degrees using a built-in PostgreSQL function degrees(), as shown in the example.</para> 
+			<para>Availability: 2.5.0</para> 
+		</refsection>
+
+		<refsection>
+		<title>Examples</title>
+		<para>Geometry Azimuth in degrees </para>
+<programlisting>
+	WITH rand AS (
+		SELECT s, random() * 2 * PI() AS rad1
+			, random() * 2 * PI() AS rad2
+		FROM  generate_series(1,2,2) AS s
+	)
+	 , points AS ( 
+		SELECT s, rad1,rad2, ST_MakePoint(cos1+s,sin1+s) as p1, ST_MakePoint(s,s) AS p2, ST_MakePoint(cos2+s,sin2+s) as p3
+		FROM rand 
+			,cos(rad1) cos1, sin(rad1) sin1 
+			,cos(rad2) cos2, sin(rad2) sin2
+	)
+	SELECT s, ST_AsText(ST_SnapToGrid(ST_MakeLine(ARRAY[p1,p2,p3]),0.001)) AS line
+		, degrees(ST_Angle(p1,p2,p3)) as computed_angle
+		, round(degrees(2*PI()-rad2 -2*PI()+rad1+2*PI()))::int%360 AS reference 
+		, round(degrees(2*PI()-rad2 -2*PI()+rad1+2*PI()))::int%360 AS reference 
+	FROM points ;
+
+1 | line | computed_angle | reference
+------------------+------------------
+1 | LINESTRING(1.511 1.86,1 1,0.896 0.005) | 155.27033848688 | 155
+
+</programlisting> 
+		</refsection> 
   </refentry>
 
   <refentry id="ST_Centroid">
diff --git a/postgis/lwgeom_functions_basic.c b/postgis/lwgeom_functions_basic.c
index 0f9d6ac97..a81794f7e 100644
--- a/postgis/lwgeom_functions_basic.c
+++ b/postgis/lwgeom_functions_basic.c
@@ -102,6 +102,7 @@ Datum LWGEOM_setpoint_linestring(PG_FUNCTION_ARGS);
 Datum LWGEOM_asEWKT(PG_FUNCTION_ARGS);
 Datum LWGEOM_hasBBOX(PG_FUNCTION_ARGS);
 Datum LWGEOM_azimuth(PG_FUNCTION_ARGS);
+Datum LWGEOM_angle(PG_FUNCTION_ARGS); 
 Datum LWGEOM_affine(PG_FUNCTION_ARGS);
 Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS);
 Datum optimistic_overlap(PG_FUNCTION_ARGS);
@@ -2420,6 +2421,131 @@ Datum LWGEOM_azimuth(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT8(result);
 }
 
+/**
+ * Compute the angle defined by 3 points or the angle between 2 vectors
+ * defined by 4 points
+ * given Point geometries.
+ * @return NULL on exception (same point).
+ * 		Return radians otherwise (always positive).
+ */
+PG_FUNCTION_INFO_V1(LWGEOM_angle);
+Datum LWGEOM_angle(PG_FUNCTION_ARGS)
+{
+	GSERIALIZED * seri_geoms[4];
+	LWGEOM *geom_unser;
+	LWPOINT *lwpoint;
+	POINT2D points[4]; 
+	double az1,az2 ;
+	double result;
+	int srids[4];
+	int i = 0 ; 
+	int j = 0;
+	int err_code = 0; 
+	int n_args = PG_NARGS(); 
+
+	/* no deserialize, checking for common error first*/
+	for(i=0; i<n_args; i++)
+	{
+		seri_geoms[i] = PG_GETARG_GSERIALIZED_P(i); 
+		if (gserialized_is_empty(seri_geoms[i]) )
+		{/* empty geom */
+			if (i==3)
+			{
+				n_args = 3 ;
+			}
+			else
+			{
+				err_code = 1 ;
+				break ;
+			}
+		} else
+		{
+			if(gserialized_get_type(seri_geoms[i]) != POINTTYPE)
+			{/* geom type */
+				err_code = 2 ;
+				break;
+			}
+			else
+			{
+				srids[i] = gserialized_get_srid(seri_geoms[i]) ; 
+				if(srids[0] != srids[i])
+				{/* error on srid*/
+					err_code = 3 ;
+					break;
+				}
+			}
+		}
+	}
+	if (err_code >0)
+		switch (err_code){
+			default: /*always executed*/
+			for (j=0;j<=i;j++)
+				PG_FREE_IF_COPY(seri_geoms[j], j);
+			
+			case 1:
+			lwpgerror("Empty geometry"); 
+			PG_RETURN_NULL() ;
+			break;
+			
+			case 2:
+			lwpgerror("Argument must be POINT geometries");
+			PG_RETURN_NULL();
+			break;
+			
+			case 3:
+			lwpgerror("Operation on mixed SRID geometries");
+			PG_RETURN_NULL();
+			break;
+		}
+	/* extract points */
+	for(i=0; i<n_args; i++)
+	{ 
+		geom_unser = lwgeom_from_gserialized(seri_geoms[i]) ; 
+		lwpoint = lwgeom_as_lwpoint(geom_unser); 
+		if (!lwpoint)
+		{ 
+			for (j=0;j<n_args;j++)
+				PG_FREE_IF_COPY(seri_geoms[j], j);
+			lwpgerror("Error unserializing geometry"); 
+			PG_RETURN_NULL() ;
+		} 
+		
+		if ( ! getPoint2d_p(lwpoint->point, 0, &points[i]) )
+		{ 
+			/* // can't free serialized geom, it might be needed by lw
+			for (j=0;j<n_args;j++)
+				PG_FREE_IF_COPY(seri_geoms[j], j); */
+			lwpgerror("Error extracting point");
+			PG_RETURN_NULL();
+		}
+		/* lwfree(geom_unser);don't do, lw may rely on this memory
+		lwpoint_free(lwpoint); dont do , this memory is needed ! */ 
+	}
+	/* // can't free serialized geom, it might be needed by lw
+	for (j=0;j<n_args;j++)
+		PG_FREE_IF_COPY(seri_geoms[j], j); */
+	
+	/* compute azimuth for the 2 pairs of points
+	 * note that angle is not defined identically for 3 points or 4 points*/ 
+	if (n_args == 3)
+	{/* we rely on azimuth to complain if points are identical */
+		if ( ! azimuth_pt_pt(&points[0], &points[1], &az1) )
+			PG_RETURN_NULL(); 
+		if ( ! azimuth_pt_pt(&points[2], &points[1], &az2) )
+			PG_RETURN_NULL(); 
+	} else
+	{
+		if ( ! azimuth_pt_pt(&points[0], &points[1], &az1) )
+			PG_RETURN_NULL();
+		if ( ! azimuth_pt_pt(&points[2], &points[3], &az2) )
+			PG_RETURN_NULL(); 
+	}
+	result = az2-az1 ;
+	result += (result<0) * 2 * M_PI ; /* we dont want negative angle*/
+	PG_RETURN_FLOAT8(result);
+}
+
+
 /*
  * optimistic_overlap(Polygon P1, Multipolygon MP2, double dist)
  * returns true if P1 overlaps MP2
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index d9df91f94..8775fc6e2 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -1290,7 +1290,15 @@ CREATE OR REPLACE FUNCTION ST_azimuth(geom1 geometry, geom2 geometry)
 	RETURNS float8
 	AS 'MODULE_PATHNAME', 'LWGEOM_azimuth'
 	LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
+	
+-- Availability: 2.5.0
+CREATE OR REPLACE FUNCTION ST_Angle(pt1 geometry, pt2 geometry, pt3 geometry, pt4 geometry default 'POINT EMPTY'::geometry)
+	RETURNS float8
+	AS 'MODULE_PATHNAME', 'LWGEOM_angle'
+	LANGUAGE 'c' IMMUTABLE STRICT;
 
+	 
+	 
 ------------------------------------------------------------------------
 -- MISC
 ------------------------------------------------------------------------
@@ -5974,6 +5982,17 @@ CREATE OR REPLACE FUNCTION ST_AsX3D(geom geometry, maxdecimaldigits integer DEFA
 	AS $$SELECT @extschema@._ST_AsX3D(3,$1,$2,$3,'');$$
 	LANGUAGE 'sql' IMMUTABLE  _PARALLEL;
 
+
+-----------------------------------------------------------------------
+-- ST_Angle
+-----------------------------------------------------------------------
+-- Availability: 2.3.0
+-- has to be here because need ST_StartPoint
+CREATE OR REPLACE FUNCTION ST_Angle(line1 geometry, line2 geometry)
+	RETURNS float8 AS 'SELECT ST_Angle(St_StartPoint($1), ST_EndPoint($1), St_StartPoint($2), ST_EndPoint($2))'
+	LANGUAGE 'sql' IMMUTABLE STRICT;
+
+
 -- make views and spatial_ref_sys public viewable --
 GRANT SELECT ON TABLE geography_columns TO public;
 GRANT SELECT ON TABLE geometry_columns TO public;
diff --git a/regress/lwgeom_regress.sql b/regress/lwgeom_regress.sql
index f66a4766b..f7f149d3b 100644
--- a/regress/lwgeom_regress.sql
+++ b/regress/lwgeom_regress.sql
@@ -135,3 +135,37 @@ SELECT 'BoundingDiagonal5', ST_AsEwkt(ST_BoundingDiagonal(
 SELECT 'BoundingDiagonal6', ST_AsEwkt(ST_BoundingDiagonal(
     'SRID=3857;POLYGON M EMPTY'::geometry
 ));
+
+
+--- ST_Azimuth
+SELECT 'ST_Azimuth_regular' , round(ST_Azimuth(geom1,geom2)::numeric,4)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, CAST('POINT(1 0)' AS geometry) AS geom2 ; 
+SELECT 'ST_Azimuth_same_point' , ST_Azimuth(geom1,geom1)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1 ; 
+SELECT 'ST_Azimuth_mixed_srid' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('POINT(1 0)',4326) AS geom2; 
+SELECT 'ST_Azimuth_not_point' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('LINESTRING(1 0 ,2 0)',4326) AS geom2; 
+SELECT 'ST_Azimuth_null_geom' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('EMPTY') AS geom2; 
+
+--- ST_Angle( points)
+SELECT 'ST_Angle_4_pts', St_Angle(p1,p2,p3,p4)
+	FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+	, ST_GeomFromtext('POINT(1 0)') AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+SELECT 'ST_Angle_4_pts', St_Angle(p1,p2,p3,p4)
+	FROM ST_GeomFromtext('POINT(2 0)') AS p1, ST_GeomFromtext('POINT(1 0)') AS p2
+	, ST_GeomFromtext('POINT(1 -1)') AS p3, ST_GeomFromtext('POINT(0 0)') AS p4 ;
+SELECT 'ST_Angle_3_pts', St_Angle(p1,p2,p3)
+	FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+	, ST_GeomFromtext('POINT(1 0)') AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ; 
+SELECT 'ST_Angle_mixed_srid', St_Angle(p1,p2,p3,p4)
+	FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+	, ST_GeomFromtext('POINT(1 0)',4326) AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+SELECT 'ST_Angle_empty' , St_Angle(p1,p2,p3,p4)
+	FROM ST_GeomFromtext('POINT EMPTY') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+	, ST_GeomFromtext('POINT(1 0)',4326) AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+--- ST_Angle( lines)
+SELECT 'ST_Angle_2_lines', St_Angle(l1,l2)
+	FROM ST_GeomFromtext('LINESTRING(0 1,0 0)') AS l1
+	, ST_GeomFromtext('LINESTRING(1 0, 2 0)') AS l2 ;
diff --git a/regress/lwgeom_regress_expected b/regress/lwgeom_regress_expected
index 5e20f9c69..4c9408eee 100644
--- a/regress/lwgeom_regress_expected
+++ b/regress/lwgeom_regress_expected
@@ -21,3 +21,14 @@ BoundingDiagonal3|SRID=4326;LINESTRING(999999986991104 999999986991104,1.0000000
 BoundingDiagonal4|SRID=3857;LINESTRING(-1 -2 -8 2,1 2 3 9)
 BoundingDiagonal5|SRID=3857;LINESTRINGM(4 4 0,5 4 1)
 BoundingDiagonal6|SRID=3857;LINESTRINGM EMPTY
+ST_Azimuth_regular|2.3562
+ST_Azimuth_same_point|
+ERROR:  Operation on mixed SRID geometries
+ERROR:  Argument must be POINT geometries
+ERROR:  parse error - invalid geometry
+ST_Angle_4_pts|4.71238898038469
+ST_Angle_4_pts|0.785398163397448
+ST_Angle_3_pts|1.5707963267949
+ERROR:  Operation on mixed SRID geometries
+ERROR:  Empty geometry
+ST_Angle_2_lines|4.71238898038469
-- 
2.50.1