PostGIS  3.0.6dev-r@@SVN_REVISION@@
cu_measures.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright (C) 2015 Sandro Santilli <strk@kbt.io>
7  * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
8  *
9  * This is free software; you can redistribute and/or modify it under
10  * the terms of the GNU General Public Licence. See the COPYING file.
11  *
12  **********************************************************************/
13 
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include "CUnit/Basic.h"
18 
19 #include "liblwgeom_internal.h"
20 #include "cu_tester.h"
21 #include "measures.h"
22 #include "measures3d.h"
23 #include "lwtree.h"
24 
25 static LWGEOM* lwgeom_from_text(const char *str)
26 {
27  LWGEOM_PARSER_RESULT r;
28  if( LW_FAILURE == lwgeom_parse_wkt(&r, (char*)str, LW_PARSER_CHECK_NONE) )
29  return NULL;
30  return r.geom;
31 }
32 
33 #define DIST2DTEST(str1, str2, res, accepted_error) \
34  do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance2d_tolerance, accepted_error);\
35  do_test_mindistance_tolerance(str2, str1, res, __LINE__, lwgeom_mindistance2d_tolerance, accepted_error)
36 #define DIST3DTEST(str1, str2, res, accepted_error) \
37  do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance3d_tolerance, accepted_error);\
38  do_test_mindistance_tolerance(str2, str1, res, __LINE__, lwgeom_mindistance3d_tolerance, accepted_error)
39 
40 static void
41 do_test_mindistance_tolerance(char *in1,
42  char *in2,
43  double expected_res,
44  int line,
45  double (*distancef)(const LWGEOM *, const LWGEOM *, double),
46  double accepted_error)
47 {
48  LWGEOM *lw1;
49  LWGEOM *lw2;
50  double distance;
51  char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n";
52  char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n %s\n\n";
53 
54  lw1 = lwgeom_from_wkt(in1, LW_PARSER_CHECK_NONE);
55  lw2 = lwgeom_from_wkt(in2, LW_PARSER_CHECK_NONE);
56 
57  if ( ! lw1 )
58  {
59  printf(msg2, in1);
60  exit(1);
61  }
62  if ( ! lw2 )
63  {
64  printf(msg2, in2);
65  exit(1);
66  }
67 
68  FLAGS_SET_SOLID(lw1->flags, 1);
69  FLAGS_SET_SOLID(lw2->flags, 1);
70 
71  distance = distancef(lw1, lw2, 0.0);
72  lwgeom_free(lw1);
73  lwgeom_free(lw2);
74 
75  if ( fabs(distance - expected_res) > accepted_error )
76  {
77  printf(msg1, distance, expected_res, line);
78  CU_FAIL();
79  }
80  else
81  {
82  CU_PASS();
83  }
84 
85 }
86 
87 static void test_mindistance2d_tolerance(void)
88 {
89  double default_accepted_error = 0.00001;
90  double zero_accepted_error = 0.0;
91  /*
92  ** Simple case.
93  */
94  DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5, default_accepted_error);
95 
96  /*
97  ** Point vs Geometry Collection.
98  */
99  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0, default_accepted_error);
100 
101  /*
102  ** Point vs Geometry Collection Collection.
103  */
104  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0, default_accepted_error);
105 
106  /*
107  ** Point vs Geometry Collection Collection Collection.
108  */
109  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0, default_accepted_error);
110 
111  /*
112  ** Point vs Geometry Collection Collection Collection Multipoint.
113  */
114  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0, default_accepted_error);
115 
116  /*
117  ** Geometry Collection vs Geometry Collection
118  */
119  DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0, default_accepted_error);
120 
121  /*
122  ** Geometry Collection Collection vs Geometry Collection Collection
123  */
124  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0, default_accepted_error);
125 
126  /*
127  ** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint
128  */
129  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0, default_accepted_error);
130 
131  /*
132  ** Linestring vs its start point
133  */
134  DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0, zero_accepted_error);
135 
136  /*
137  ** Linestring vs its end point
138  */
139  DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0, zero_accepted_error);
140 
141  /*
142  ** Linestring vs its start point (tricky number, see #1459)
143  */
144  DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0, zero_accepted_error);
145 
146  /*
147  ** Linestring vs its end point (tricky number, see #1459)
148  */
149  DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0, zero_accepted_error);
150 
151  /*
152  * Circular string and point
153  */
154  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1, default_accepted_error);
155  DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1, default_accepted_error);
156 
157  /*
158  * Circular string and Circular string
159  */
160  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1, default_accepted_error);
161 
162  /*
163  * CurvePolygon and Point
164  */
165  static char *cs1 = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))";
166  DIST2DTEST(cs1, "POINT(3 14)", 1, default_accepted_error);
167  DIST2DTEST(cs1, "POINT(3 8)", 0, default_accepted_error);
168  DIST2DTEST(cs1, "POINT(6 5)", 1, default_accepted_error);
169  DIST2DTEST(cs1, "POINT(6 4)", 0, default_accepted_error);
170 
171  /*
172  * CurvePolygon and Linestring
173  */
174  DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484, default_accepted_error);
175  DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0, default_accepted_error);
176  DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0, default_accepted_error);
177  DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786, default_accepted_error);
178  DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913, default_accepted_error);
179 
180  /*
181  * CurvePolygon and Polygon
182  */
183  DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415, default_accepted_error);
184  DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0, default_accepted_error);
185  DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0, default_accepted_error);
186 
187  /*
188  * CurvePolygon and CurvePolygon
189  */
190  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666, default_accepted_error);
191  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0, default_accepted_error);
192 
193  /*
194  * MultiSurface and CurvePolygon
195  */
196  static char *cs2 = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))";
197  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1, default_accepted_error);
198  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0, default_accepted_error);
199  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0, default_accepted_error);
200  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0, default_accepted_error);
201  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5, default_accepted_error);
202  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0, default_accepted_error);
203  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0, default_accepted_error);
204  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555, default_accepted_error);
205 
206  /*
207  * MultiCurve and Linestring
208  */
209  DIST2DTEST("LINESTRING(0.5 1,0.5 3)", "MULTICURVE(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3),(0 0, 0 5))", 0.5, default_accepted_error);
210 
211  /*
212  ** Ticket 4326
213  */
214  DIST2DTEST(
215  "CURVEPOLYGON(CIRCULARSTRING(7874821 8715927,8907663 8715927,8844683 7750316,7937800 7750316,7874821 8715927))",
216  "POINT(5433865 8243495)", 2271704.2698450615, default_accepted_error);
217 }
218 
219 static void
220 test_mindistance3d_tolerance(void)
221 {
222  double default_accepted_error = 0.00001;
223  double zero_accepted_error = 0.0;
224  /* 2D [Z=0] should work just the same */
225  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5, default_accepted_error);
226  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5, default_accepted_error);
227  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0, default_accepted_error);
228  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))", 5.0, default_accepted_error);
229  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0))))", 5.0, default_accepted_error);
230  DIST3DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0, default_accepted_error);
231  DIST3DTEST("GEOMETRYCOLLECTION(POINT(0 0 0))", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0, default_accepted_error);
232  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0 0)))",
233  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))",
234  5.0, default_accepted_error);
235  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0 0)))",
236  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4 0)))",
237  5.0, default_accepted_error);
238  DIST3DTEST("LINESTRING(-2 0 0, -0.2 0 0)", "POINT(-2 0 0)", 0, zero_accepted_error);
239  DIST3DTEST("LINESTRING(-0.2 0 0, -2 0 0)", "POINT(-2 0 0)", 0, zero_accepted_error);
240  DIST3DTEST("LINESTRING(-1e-8 0 0, -0.2 0 0)", "POINT(-1e-8 0 0)", 0, zero_accepted_error);
241  DIST3DTEST("LINESTRING(-0.2 0 0, -1e-8 0 0)", "POINT(-1e-8 0 0)", 0, zero_accepted_error);
242 
243  /* Tests around intersections */
244  DIST3DTEST("LINESTRING(1 0 0 , 2 0 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0, default_accepted_error);
245  DIST3DTEST("LINESTRING(1 1 1 , 2 1 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 0.0, zero_accepted_error);
246  DIST3DTEST("LINESTRING(1 1 1 , 2 1 1)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0, default_accepted_error);
247 
248  /* Same but triangles */
249  DIST3DTEST("LINESTRING(1 0 0 , 2 0 0)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 1.0, default_accepted_error);
250  DIST3DTEST("LINESTRING(1 1 1 , 2 1 0)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 0.0, zero_accepted_error);
251  DIST3DTEST("LINESTRING(1 1 1 , 2 1 1)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 1.0, default_accepted_error);
252 
253  /* Triangle to triangle*/
254  DIST3DTEST("TRIANGLE((-1 1 0, -2 1 0, -1 2 0, -1 1 0))", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 2.0, default_accepted_error);
255 
256  /* Line in polygon */
257  DIST3DTEST("LINESTRING(1 1 1 , 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 1, 0 0 0))", 0.0, zero_accepted_error);
258 
259  /* Line has a point in the same plane as the polygon but isn't the closest*/
260  DIST3DTEST("LINESTRING(-10000 -10000 0, 0 0 1)", "POLYGON((0 0 0, 1 0 0, 1 1 0, 0 1 0, 0 0 0))", 1, default_accepted_error);
261 
262  /* This is an invalid polygon since it defines just a line */
263  DIST3DTEST("LINESTRING(1 1 1, 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0, zero_accepted_error);
264  DIST3DTEST("TRIANGLE((1 1 1, 2 2 2, 3 3 3, 1 1 1))", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0, zero_accepted_error);
265  DIST3DTEST("POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", "TRIANGLE((1 1 1, 2 2 2, 3 3 3, 1 1 1))", 0, zero_accepted_error);
266  DIST3DTEST("TRIANGLE((0 0 0, 2 2 2, 3 3 3, 0 0 0))", "LINESTRING(1 1 1, 2 2 2)", 0, zero_accepted_error);
267 
268  /* A box in a box: two solids, one inside another */
269  DIST3DTEST(
270  "POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
271  "POLYHEDRALSURFACE Z (((-1 -1 -1,-1 2 -1,2 2 -1,2 -1 -1,-1 -1 -1)),((-1 -1 2,2 -1 2,2 2 2,-1 2 2,-1 -1 2)),((-1 -1 -1,-1 -1 2,-1 2 2,-1 2 -1,-1 -1 -1)),((2 -1 -1,2 2 -1,2 2 2,2 -1 2,2 -1 -1)),((-1 -1 -1,2 -1 -1,2 -1 2,-1 -1 2,-1 -1 -1)),((-1 2 -1,-1 2 2,2 2 2,2 2 -1,-1 2 -1)))",
272  0, zero_accepted_error);
273 
274  /* A box in a box with a hat: two solids, one inside another, Z ray up is hitting hat */
275  DIST3DTEST(
276  "POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
277  "POLYHEDRALSURFACE Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
278  0, zero_accepted_error);
279 
280  /* Same but as TIN */
281  DIST3DTEST(
282  "TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
283  "POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
284  0, zero_accepted_error);
285 
286  /* Same but both are TIN */
287  DIST3DTEST(
288  "TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
289  "TIN Z (((0 0 0,0 1 0,1 1 0,0 0 0)),((1 0 0,0 0 0,1 1 0,1 0 0)),((0 1 1,1 0 1,1 1 1,0 1 1)),((0 1 1,0 0 1,1 0 1,0 1 1)),((0 0 0,0 0 1,0 1 1,0 0 0)),((0 1 0,0 0 0,0 1 1,0 1 0)),((1 0 1,1 1 0,1 1 1,1 0 1)),((1 0 1,1 0 0,1 1 0,1 0 1)),((0 0 1,1 0 0,1 0 1,0 0 1)),((0 0 1,0 0 0,1 0 0,0 0 1)),((0 1 0,0 1 1,1 1 1,0 1 0)),((1 1 0,0 1 0,1 1 1,1 1 0)))",
290  0, zero_accepted_error);
291 
292  /* Point inside TIN */
293  DIST3DTEST(
294  "TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
295  "POINT(0 0 0)",
296  0, zero_accepted_error);
297 
298  /* A point hits vertical Z edge */
299  DIST3DTEST(
300  "POLYHEDRALSURFACE Z (((0 -1 1,-1 -1 1,-1 -1 -1,0 -1 -1,1 -1 -1,0 -1 2,0 -1 1)),((0 1 1,0 1 2,1 1 -1,0 1 -1,-1 1 -1,-1 1 1,0 1 1)),((0 -1 1,0 1 1,-1 1 1,-1 -1 1,0 -1 1)),((-1 -1 1,-1 1 1,-1 1 -1,-1 -1 -1,-1 -1 1)),((-1 -1 -1,-1 1 -1,0 1 -1,0 -1 -1,-1 -1 -1)),((0 -1 -1,0 1 -1,1 1 -1,1 -1 -1,0 -1 -1)),((1 -1 -1,1 1 -1,0 1 2,0 -1 2,1 -1 -1)),((0 -1 2,0 1 2,0 1 1,0 -1 1,0 -1 2)))",
301  "POINT(0 0 0)",
302  0, zero_accepted_error);
303 
304  /* A point in the middle of a hole of extruded polygon */
305  DIST3DTEST(
306  "POLYHEDRALSURFACE Z (((-3 -3 0,-3 3 0,3 3 0,3 -3 0,-3 -3 0),(-1 -1 0,1 -1 0,1 1 0,-1 1 0,-1 -1 0)),((-3 -3 3,3 -3 3,3 3 3,-3 3 3,-3 -3 3),(-1 -1 3,-1 1 3,1 1 3,1 -1 3,-1 -1 3)),((-3 -3 0,-3 -3 3,-3 3 3,-3 3 0,-3 -3 0)),((-3 3 0,-3 3 3,3 3 3,3 3 0,-3 3 0)),((3 3 0,3 3 3,3 -3 3,3 -3 0,3 3 0)),((3 -3 0,3 -3 3,-3 -3 3,-3 -3 0,3 -3 0)),((-1 -1 0,-1 -1 3,1 -1 3,1 -1 0,-1 -1 0)),((1 -1 0,1 -1 3,1 1 3,1 1 0,1 -1 0)),((1 1 0,1 1 3,-1 1 3,-1 1 0,1 1 0)),((-1 1 0,-1 1 3,-1 -1 3,-1 -1 0,-1 1 0)))",
307  "POINT(0 0 1)",
308  1, zero_accepted_error);
309 
310  /* A point at the face of a hole of extruded polygon */
311  DIST3DTEST(
312  "POLYHEDRALSURFACE Z (((-3 -3 0,-3 3 0,3 3 0,3 -3 0,-3 -3 0),(-1 -1 0,1 -1 0,1 1 0,-1 1 0,-1 -1 0)),((-3 -3 3,3 -3 3,3 3 3,-3 3 3,-3 -3 3),(-1 -1 3,-1 1 3,1 1 3,1 -1 3,-1 -1 3)),((-3 -3 0,-3 -3 3,-3 3 3,-3 3 0,-3 -3 0)),((-3 3 0,-3 3 3,3 3 3,3 3 0,-3 3 0)),((3 3 0,3 3 3,3 -3 3,3 -3 0,3 3 0)),((3 -3 0,3 -3 3,-3 -3 3,-3 -3 0,3 -3 0)),((-1 -1 0,-1 -1 3,1 -1 3,1 -1 0,-1 -1 0)),((1 -1 0,1 -1 3,1 1 3,1 1 0,1 -1 0)),((1 1 0,1 1 3,-1 1 3,-1 1 0,1 1 0)),((-1 1 0,-1 1 3,-1 -1 3,-1 -1 0,-1 1 0)))",
313  "POINT(1 1 2)",
314  0, zero_accepted_error);
315 
316  /* A point at the face of a hole of extruded polygon */
317  DIST3DTEST(
318  "LINESTRING Z (-27974.1264 -110211.5032 148.9768,-27975.4229 -110210.9441 149.0093)",
319  "LINESTRING Z (-27995.4183 -110201.8041 149.3354,-27975.4229 -110210.9441 149.0093)",
320  0, zero_accepted_error);
321 }
322 
323 static int tree_pt(RECT_NODE *tree, double x, double y)
324 {
325  POINT2D pt;
326  pt.x = x; pt.y = y;
327  return rect_tree_contains_point(tree, &pt);
328 }
329 
330 static void test_rect_tree_contains_point(void)
331 {
332  LWGEOM *poly;
333  RECT_NODE* tree;
334 
335  /**********************************************************************
336  * curvepolygon
337  */
338  poly = lwgeom_from_wkt("CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,1 5,0 10),(0 10,10 10,10 0, 0 0)),COMPOUNDCURVE(CIRCULARSTRING(3 7,5 8,7 7),(7 7,7 3,3 3, 3 7)))", LW_PARSER_CHECK_NONE);
339  tree = rect_tree_from_lwgeom(poly);
340  // char *wkt = rect_tree_to_wkt(tree);
341  // printf("%s\n", wkt);
342  // lwfree(wkt);
343  // return;
344 
345  /* in hole, within arc stroke */
346  CU_ASSERT_EQUAL(tree_pt(tree, 5, 7.5), 0);
347  /* inside */
348  CU_ASSERT_EQUAL(tree_pt(tree, 8, 9), 1);
349  /* outside */
350  CU_ASSERT_EQUAL(tree_pt(tree, -1, 5), 0);
351  /* outside */
352  CU_ASSERT_EQUAL(tree_pt(tree, -1, 7.5), 0);
353  /* outside, within arc stroke */
354  CU_ASSERT_EQUAL(tree_pt(tree, 0.2, 7.5), 0);
355  /* inside, within arc stroke */
356  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
357  /* inside, crossing arc stroke */
358  CU_ASSERT_EQUAL(tree_pt(tree, 2, 7.5), 1);
359  /* touching hole corner */
360  CU_ASSERT_EQUAL(tree_pt(tree, 7, 7), 1);
361 
362  rect_tree_free(tree);
363  lwgeom_free(poly);
364 
365 
366  /**********************************************************************
367  * polygon with hole and concavities
368  */
369  poly = lwgeom_from_wkt("POLYGON((0 0,0 10,10 10,10 0,9 0,9 9,8 6,8 0,2 0,2 9,1 6,1 0,0 0),(4 4,4 6,6 6,6 4,4 4))", LW_PARSER_CHECK_NONE);
370  tree = rect_tree_from_lwgeom(poly);
371 
372  /* inside, many grazings */
373  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
374  /* inside */
375  CU_ASSERT_EQUAL(tree_pt(tree, 3, 5.5), 1);
376  /* outside */
377  CU_ASSERT_EQUAL(tree_pt(tree, -3, 5.5), 0);
378  /* touching interior ring */
379  CU_ASSERT_EQUAL(tree_pt(tree, 4, 4), 1);
380  CU_ASSERT_EQUAL(tree_pt(tree, 6, 6), 1);
381  /* touching interior ring */
382  CU_ASSERT_EQUAL(tree_pt(tree, 4.5, 4), 1);
383  /* touching exterior ring */
384  CU_ASSERT_EQUAL(tree_pt(tree, 8, 0), 1);
385  CU_ASSERT_EQUAL(tree_pt(tree, 9, 0), 1);
386  CU_ASSERT_EQUAL(tree_pt(tree, 10, 1), 1);
387  CU_ASSERT_EQUAL(tree_pt(tree, 9.5, 1), 1);
388  CU_ASSERT_EQUAL(tree_pt(tree, 0, 10), 1);
389  /* touching grazing spike */
390  CU_ASSERT_EQUAL(tree_pt(tree, 1, 6), 1);
391  /* outide, many grazings */
392  CU_ASSERT_EQUAL(tree_pt(tree, -1, 6), 0);
393  /* within hole */
394  CU_ASSERT_EQUAL(tree_pt(tree, 5, 5), 0);
395  /* within */
396  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 4), 1);
397  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 6), 1);
398  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 9), 1);
399 
400  rect_tree_free(tree);
401  lwgeom_free(poly);
402 
403  /**********************************************************************
404  * square
405  */
406  poly = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
407  tree = rect_tree_from_lwgeom(poly);
408 
409  /* inside square */
410  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
411  /* outside square */
412  CU_ASSERT_EQUAL(tree_pt(tree, 1.5, 0.5), 0);
413  /* outside square grazing some edges */
414  CU_ASSERT_EQUAL(tree_pt(tree, -1, 1), 0);
415  /* inside square on corner */
416  CU_ASSERT_EQUAL(tree_pt(tree, 1, 1), 1);
417  /* inside square on top edge */
418  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 1), 1);
419  /* inside square on side edge */
420  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0.5), 1);
421 
422  rect_tree_free(tree);
423  lwgeom_free(poly);
424 
425  /* ziggy zaggy horizontal saw tooth polygon */
426  poly = lwgeom_from_wkt("POLYGON((0 0, 1 3, 2 0, 3 3, 4 0, 4 5, 0 5, 0 0))", LW_PARSER_CHECK_NONE);
427  tree = rect_tree_from_lwgeom(poly);
428 
429  /* not in, left side */
430  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 0.5), 0);
431  /* not in, right side */
432  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 0);
433  /* inside */
434  CU_ASSERT_EQUAL(tree_pt(tree, 2, 1), 1);
435  /* on left border */
436  CU_ASSERT_EQUAL(tree_pt(tree, 0, 1), 1);
437  /* on left border, grazing */
438  CU_ASSERT_EQUAL(tree_pt(tree, 0, 3), 1);
439  /* on right border */
440  CU_ASSERT_EQUAL(tree_pt(tree, 4, 0), 1);
441  /* on tooth concave */
442  CU_ASSERT_EQUAL(tree_pt(tree, 3, 3), 1);
443  /* on tooth convex */
444  CU_ASSERT_EQUAL(tree_pt(tree, 2, 0), 1);
445 
446  rect_tree_free(tree);
447  lwgeom_free(poly);
448 
449  /**********************************************************************
450  * ziggy zaggy vertical saw tooth polygon
451  */
452  poly = lwgeom_from_wkt("POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))", LW_PARSER_CHECK_NONE);
453  tree = rect_tree_from_lwgeom(poly);
454 
455  /* not in, left side */
456  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 3.5), 0);
457  /* not in, right side */
458  CU_ASSERT_EQUAL(tree_pt(tree, 6.0, 2.2), 0);
459  /* inside */
460  CU_ASSERT_EQUAL(tree_pt(tree, 3, 2), 1);
461  /* on bottom border */
462  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0), 1);
463  /* on top border */
464  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
465  /* on tooth concave */
466  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 1);
467  /* on tooth convex */
468  CU_ASSERT_EQUAL(tree_pt(tree, 0, 2), 1);
469  /* on tooth convex */
470  CU_ASSERT_EQUAL(tree_pt(tree, 0, 6), 1);
471 
472  rect_tree_free(tree);
473  lwgeom_free(poly);
474 
475 }
476 
477 static int tree_inter(const char *wkt1, const char *wkt2)
478 {
479  LWGEOM *g1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
480  LWGEOM *g2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
481  RECT_NODE *t1 = rect_tree_from_lwgeom(g1);
482  RECT_NODE *t2 = rect_tree_from_lwgeom(g2);
483  int result = rect_tree_intersects_tree(t1, t2);
484  rect_tree_free(t1);
485  rect_tree_free(t2);
486  lwgeom_free(g1);
487  lwgeom_free(g2);
488  return result;
489 }
490 
491 static void test_rect_tree_intersects_tree(void)
492 {
493  /* total overlap, A == B */
494  CU_ASSERT_EQUAL(tree_inter(
495  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
496  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))"),
497  LW_TRUE
498  );
499 
500  /* hiding between the tines of the comb */
501  CU_ASSERT_EQUAL(tree_inter(
502  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
503  "POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 0.4 0.7, 0.3 0.7))"),
504  LW_FALSE
505  );
506 
507  /* between the tines, but with a corner overlapping */
508  CU_ASSERT_EQUAL(tree_inter(
509  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
510  "POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 1.3 0.3, 0.3 0.7))"),
511  LW_TRUE
512  );
513 
514  /* Just touching the top left corner of the comb */
515  CU_ASSERT_EQUAL(tree_inter(
516  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
517  "POLYGON((-1 5, 0 5, 0 7, -1 7, -1 5))"),
518  LW_TRUE
519  );
520 
521  /* Contained, complex */
522  CU_ASSERT_EQUAL(tree_inter(
523  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
524  "GEOMETRYCOLLECTION(MULTILINESTRING((1 2, 3 2)),POINT(1 2))"),
525  LW_TRUE
526  );
527 
528  /* Touching, complex */
529  CU_ASSERT_EQUAL(tree_inter(
530  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
531  "GEOMETRYCOLLECTION(MULTILINESTRING((6 3, 8 4)),POINT(5 3))"),
532  LW_TRUE
533  );
534 
535  /* Not Touching, complex */
536  CU_ASSERT_EQUAL(tree_inter(
537  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
538  "GEOMETRYCOLLECTION(MULTILINESTRING((6 3, 8 4)),POINT(1 3.5))"),
539  LW_FALSE
540  );
541 
542  /* Crossing, complex */
543  CU_ASSERT_EQUAL(tree_inter(
544  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
545  "GEOMETRYCOLLECTION(MULTILINESTRING((1.5 4.1, 1.6 2)),POINT(1 3.5))"),
546  LW_TRUE
547  );
548 }
549 
550 
551 static double
552 test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
553 {
554  LWGEOM *lw1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
555  LWGEOM *lw2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
556  RECT_NODE *n1 = rect_tree_from_lwgeom(lw1);
557  RECT_NODE *n2 = rect_tree_from_lwgeom(lw2);
558 
559  // rect_tree_printf(n1, 0);
560  // rect_tree_printf(n2, 0);
561  //
562  // printf("%s\n", rect_tree_to_wkt(n1));
563  // printf("%s\n", rect_tree_to_wkt(n2));
564 
565  double dist = rect_tree_distance_tree(n1, n2, 0.0);
566  // printf("%g\n", dist);
567  rect_tree_free(n1);
568  rect_tree_free(n2);
569  lwgeom_free(lw1);
570  lwgeom_free(lw2);
571  return dist;
572 }
573 
574 #define TDT(w1, w2, d) CU_ASSERT_DOUBLE_EQUAL(test_rect_tree_distance_tree_case(w1, w2), d, 0.00001);
575 
576 static void
577 test_rect_tree_distance_tree(void)
578 {
579  const char *wkt;
580 
581  wkt = "MULTIPOLYGON(((-123.35702791281 48.4232302445918,-123.35689654493 48.4237265810249,-123.354053908057 48.4234039978588,-123.35417179975 48.4229151379279,-123.354369811539 48.4220987102936,-123.355779071731 48.4222571534228,-123.357238860904 48.4224209369449,-123.35702791281 48.4232302445918)))";
582  TDT(wkt, "MULTIPOLYGON(((-123.353452578038 48.4259519079838,-123.35072012771 48.4256699150083,-123.347337809991 48.4254740864963,-123.347469111645 48.4245757659326,-123.349409235923 48.4246224093429,-123.349966167324 48.4246562342604,-123.353650661317 48.4250703224683,-123.353452578038 48.4259519079838)))", 0.0017144228293396);
583 
584  wkt = "POINT(0 0)";
585  TDT(wkt, "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
586  TDT(wkt, "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
587  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
588  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
589  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
590 
591  TDT("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
592  TDT("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
593  TDT("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
594  TDT("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
595 
596  wkt = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))";
597  TDT(wkt, "POINT(3 14)", 1);
598  TDT(wkt, "POINT(3 8)", 0);
599  TDT(wkt, "POINT(6 5)", 1);
600  TDT(wkt, "POINT(6 4)", 0);
601 
602  wkt = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))";
603  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
604  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
605  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
606  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
607  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
608  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
609  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
610  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
611 
612  wkt = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,5 0,0 0)))";
613  TDT(wkt, "POINT(3 0)", 0.0);
614  TDT(wkt, "POINT(5 0)", 0.0);
615  TDT(wkt, "POINT(7 0)", 2.0);
616  TDT(wkt, "POINT(2.5 3.5)", 1.0);
617 
618  wkt = "POINT(0 0)";
619  TDT(wkt, "POINT(0 1)", 1.0);
620  TDT(wkt, "POINT(1 0)", 1.0);
621 
622  wkt = "LINESTRING(0 0,1 0)";
623  TDT(wkt, "LINESTRING(1 0,1 1)", 0.0);
624  TDT(wkt, "LINESTRING(0 1,1 1)", 1.0);
625 
626  wkt = "POLYGON((0 0,0 1,1 1,1 0,0 0))";
627  TDT(wkt, "POINT(2 2)", sqrt(2));
628  TDT(wkt, "POINT(0.5 0.5)", 0);
629  TDT(wkt, "POINT(1 1)", 0);
630 
631  wkt = "POLYGON((0 0,0 10,10 10,10 0,0 0), (4 4,4 6,6 6,6 4,4 4))";
632  TDT(wkt, "POINT(5 5)", 1);
633  TDT(wkt, "POLYGON((5 5,5 5.5,5.5 5.5,5.5 5, 5 5))", 0.5);
634 }
635 
636 
637 static void
638 test_lwgeom_segmentize2d(void)
639 {
640  LWGEOM *linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
641  LWGEOM *lineout = lwgeom_segmentize2d(linein, 5);
642  char *strout = lwgeom_to_ewkt(lineout);
643  ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0)");
644  lwfree(strout);
645  lwgeom_free(linein);
646  lwgeom_free(lineout);
647 
648  /* test that segmentize is proportional - request every 6, get every 5 */
649  linein = lwgeom_from_wkt("LINESTRING(0 0, 20 0)", LW_PARSER_CHECK_NONE);
650  lineout = lwgeom_segmentize2d(linein, 6);
651  strout = lwgeom_to_ewkt(lineout);
652  ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0,15 0,20 0)");
653  lwfree(strout);
654  lwgeom_free(linein);
655  lwgeom_free(lineout);
656 
657  /* test too many segments */
658  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
659  lineout = lwgeom_segmentize2d(linein, 1e-100);
660  CU_ASSERT_EQUAL(lineout, NULL);
661  lwgeom_free(linein);
662 
663  /* test interruption */
664 
665  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
666  lwgeom_request_interrupt();
667  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
668  CU_ASSERT_EQUAL(lineout, NULL);
669  lwgeom_free(linein);
670 
671  linein = lwgeom_from_wkt("MULTILINESTRING((0 0,10 0),(20 0, 30 0))", LW_PARSER_CHECK_NONE);
672  lwgeom_request_interrupt();
673  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
674  CU_ASSERT_EQUAL(lineout, NULL);
675  lwgeom_free(linein);
676 
677  linein = lwgeom_from_wkt(
678  "MULTIPOLYGON(((0 0,20 0,20 20,0 20,0 0),(2 2,2 4,4 4,4 2,2 2),(6 6,6 8,8 8,8 6,6 6)),((40 0,40 20,60 20,60 0,40 0),(42 2,42 4,44 4,44 2,42 2)))"
679  , LW_PARSER_CHECK_NONE);
680  lwgeom_request_interrupt();
681  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
682  CU_ASSERT_EQUAL(lineout, NULL);
683  lwgeom_free(linein);
684 
685  linein = lwgeom_from_wkt(
686  "GEOMETRYCOLLECTION(MULTIPOLYGON(((0 0,20 0,20 20,0 20,0 0),(2 2,2 4,4 4,4 2,2 2),(6 6,6 8,8 8,8 6,6 6)),((40 0,40 20,60 20,60 0,40 0),(42 2,42 4,44 4,44 2,42 2))),MULTILINESTRING((0 0,10 0),(20 0, 30 0)),MULTIPOINT(0 0, 3 4))"
687  , LW_PARSER_CHECK_NONE);
688  CU_ASSERT_FATAL(linein != NULL);
689  lwgeom_request_interrupt();
690  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
691  CU_ASSERT_EQUAL(lineout, NULL);
692  lwgeom_free(linein);
693 
694  linein = lwgeom_from_wkt("LINESTRING(20 0, 30 0)", LW_PARSER_CHECK_NONE);
695  CU_ASSERT_FATAL(linein != NULL);
696  /* NOT INTERRUPTED */
697  lineout = lwgeom_segmentize2d(linein, 5);
698  CU_ASSERT_NOT_EQUAL_FATAL(lineout, NULL);
699  strout = lwgeom_to_ewkt(lineout);
700  ASSERT_STRING_EQUAL(strout, "LINESTRING(20 0,25 0,30 0)");
701  lwfree(strout);
702  lwgeom_free(linein);
703  lwgeom_free(lineout);
704 }
705 
706 static void
707 test_lwgeom_locate_along(void)
708 {
709  LWGEOM *geom = NULL;
710  LWGEOM *out = NULL;
711  double measure = 105.0;
712  char *str;
713 
714  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
715  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
716  out = lwgeom_locate_along(geom, measure, 0.0);
717  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
718  lwgeom_free(geom);
719  lwgeom_free(out);
720  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
721  lwfree(str);
722 
723  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
724  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 3 4 2, 9 4 3), (1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
725  out = lwgeom_locate_along(geom, measure, 0.0);
726  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
727  lwgeom_free(geom);
728  lwgeom_free(out);
729  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
730  lwfree(str);
731 }
732 
733 static void
734 test_lw_dist2d_pt_arc(void)
735 {
736  /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
737  DISTPTS dl;
738  POINT2D P, A1, A2, A3;
739  int rv;
740 
741 
742  /* Point within unit semicircle, 0.5 units from arc */
743  A1.x = -1; A1.y = 0;
744  A2.x = 0 ; A2.y = 1;
745  A3.x = 1 ; A3.y = 0;
746  P.x = 0 ; P.y = 0.5;
747 
748  lw_dist2d_distpts_init(&dl, DIST_MIN);
749  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
750  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
751  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
752 
753  /* Point outside unit semicircle, 0.5 units from arc */
754  P.x = 0 ; P.y = 1.5;
755  lw_dist2d_distpts_init(&dl, DIST_MIN);
756  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
757  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
758  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
759 
760  /* Point outside unit semicircle, sqrt(2) units from arc end point*/
761  P.x = 0 ; P.y = -1;
762  lw_dist2d_distpts_init(&dl, DIST_MIN);
763  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
764  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
765  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
766 
767  /* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
768  P.x = 1 ; P.y = 1;
769  lw_dist2d_distpts_init(&dl, DIST_MIN);
770  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
771  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
772  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
773 
774  /* Point on unit semicircle midpoint */
775  P.x = 0 ; P.y = 1;
776  lw_dist2d_distpts_init(&dl, DIST_MIN);
777  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
778  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
779  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
780 
781  /* Point on unit semicircle endpoint */
782  P.x = 1 ; P.y = 0;
783  lw_dist2d_distpts_init(&dl, DIST_MIN);
784  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
785  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
786  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
787 
788  /* Point on semicircle center */
789  P.x = 0 ; P.y = 0;
790  lw_dist2d_distpts_init(&dl, DIST_MIN);
791  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
792  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
793  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
794 
795  /* Point inside closed circle */
796  P.x = 0 ; P.y = 0.5;
797  A2.x = 1; A2.y = 0;
798  A3 = A1;
799  lw_dist2d_distpts_init(&dl, DIST_MIN);
800  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
801  //printf("distance %g\n", dl.distance);
802  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
803  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
804 }
805 
806 static void
807 test_lw_dist2d_seg_arc(void)
808 {
809  /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */
810 
811  DISTPTS dl;
812  POINT2D A1, A2, B1, B2, B3;
813  int rv;
814 
815  /* Unit semicircle */
816  B1.x = -1; B1.y = 0;
817  B2.x = 0 ; B2.y = 1;
818  B3.x = 1 ; B3.y = 0;
819 
820  /* Edge above the unit semicircle */
821  lw_dist2d_distpts_init(&dl, DIST_MIN);
822  A1.x = -2; A1.y = 2;
823  A2.x = 2 ; A2.y = 2;
824  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
825  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
826  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
827 
828  /* Edge to the right of the unit semicircle */
829  lw_dist2d_distpts_init(&dl, DIST_MIN);
830  A1.x = 2; A1.y = -2;
831  A2.x = 2; A2.y = 2;
832  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
833  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
834  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
835 
836  /* Edge to the left of the unit semicircle */
837  lw_dist2d_distpts_init(&dl, DIST_MIN);
838  A1.x = -2; A1.y = -2;
839  A2.x = -2; A2.y = 2;
840  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
841  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
842  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
843 
844  /* Edge within the unit semicircle */
845  lw_dist2d_distpts_init(&dl, DIST_MIN);
846  A1.x = 0; A1.y = 0;
847  A2.x = 0; A2.y = 0.5;
848  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
849  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
850  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
851 
852  /* Edge grazing the unit semicircle */
853  lw_dist2d_distpts_init(&dl, DIST_MIN);
854  A1.x = -2; A1.y = 1;
855  A2.x = 2; A2.y = 1;
856  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
857  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
858  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);
859 
860  /* Line grazing the unit semicircle, but edge not */
861  lw_dist2d_distpts_init(&dl, DIST_MIN);
862  A1.x = 1; A1.y = 1;
863  A2.x = 2; A2.y = 1;
864  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
865  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
866  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
867 
868  /* Edge intersecting the unit semicircle */
869  lw_dist2d_distpts_init(&dl, DIST_MIN);
870  A1.x = 0; A1.y = 0;
871  A2.x = 2; A2.y = 2;
872  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
873  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
874  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
875 
876  /* Line intersecting the unit semicircle, but edge not */
877  lw_dist2d_distpts_init(&dl, DIST_MIN);
878  A1.x = -1; A1.y = 1;
879  A2.x = -2; A2.y = 2;
880  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
881  //printf("distance %g\n", dl.distance);
882  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
883  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
884 }
885 
886 static void
887 test_lw_dist2d_arc_arc(void)
888 {
889  /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
890  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
891  DISTPTS *dl) */
892  DISTPTS dl;
893  POINT2D A1, A2, A3, B1, B2, B3;
894  int rv;
895 
896  /* Ticket #4326 */
897  lw_dist2d_distpts_init(&dl, DIST_MIN);
898  A1.x = -1.0; A1.y = 4.0;
899  A2.x = 0.0; A2.y = 5.0;
900  A3.x = 1.0; A3.y = 4.0;
901  B1.x = 1.0; B1.y = 6.0;
902  B2.x = 6.0; B2.y = 1.0;
903  B3.x = 9.0; B3.y = 7.0;
904  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
905  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
906  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0475666, 0.000001);
907 
908  /* Unit semicircle at 0,0 */
909  B1.x = -1; B1.y = 0;
910  B2.x = 0 ; B2.y = 1;
911  B3.x = 1 ; B3.y = 0;
912 
913  /* Arc above the unit semicircle */
914  lw_dist2d_distpts_init(&dl, DIST_MIN);
915  A1.x = -1; A1.y = 3;
916  A2.x = 0 ; A2.y = 2;
917  A3.x = 1 ; A3.y = 3;
918  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
919  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
920  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
921 
922  /* Arc grazes the unit semicircle */
923  lw_dist2d_distpts_init(&dl, DIST_MIN);
924  A1.x = -1; A1.y = 2;
925  A2.x = 0 ; A2.y = 1;
926  A3.x = 1 ; A3.y = 2;
927  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
928  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
929  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
930 
931  /* Circles intersect, but arcs do not */
932  lw_dist2d_distpts_init(&dl, DIST_MIN);
933  A1.x = -1; A1.y = 1;
934  A2.x = 0; A2.y = 2;
935  A3.x = 1; A3.y = 1;
936  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
937  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
938  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);
939 
940  /* Circles and arcs intersect */
941  lw_dist2d_distpts_init(&dl, DIST_MIN);
942  A1.x = -1; A1.y = 1;
943  A2.x = 0; A2.y = 0;
944  A3.x = 1; A3.y = 1;
945  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
946  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
947  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
948 
949  /* Concentric: and fully parallel */
950  lw_dist2d_distpts_init(&dl, DIST_MIN);
951  A1.x = -2.0; A1.y = 0.0;
952  A2.x = 0.0; A2.y = 2.0;
953  A3.x = 2.0; A3.y = 0.0;
954  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
955  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
956  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
957 
958  /* Concentric: with A fully included in B's range */
959  lw_dist2d_distpts_init(&dl, DIST_MIN);
960  A1.x = -0.5 / sqrt(2.0); A1.y = 0.5 / sqrt(2.0);
961  A2.x = 0.0; A2.y = 0.5;
962  A3.x = 0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
963  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
964  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
965  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
966 
967  /* Concentric: with A partially included in B's range */
968  lw_dist2d_distpts_init(&dl, DIST_MIN);
969  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
970  A2.x = -0.5; A2.y = 0.0;
971  A3.x = -0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
972  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
973  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
974  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
975 
976  /* Concentric: with A and B without parallel segments */
977  lw_dist2d_distpts_init(&dl, DIST_MIN);
978  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
979  A2.x = 0.0; A2.y = -0.5;
980  A3.x = 0.5 / sqrt(2.0); A3.y = -0.5 / sqrt(2.0);
981  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
982  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
983  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.736813, 0.000001);
984 
985  /* Concentric: Arcs are the same */
986  lw_dist2d_distpts_init(&dl, DIST_MIN);
987  A1.x = -1.0; A1.y = 0.0;
988  A2.x = 0.0; A2.y = 1.0;
989  A3.x = 1.0; A3.y = 0.0;
990  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
991  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
992  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0, 0.000001);
993 
994  /* Check different orientations */
995  B1.x = -10.0; B1.y = 0.0;
996  B2.x = 0.0 ; B2.y = 10.0;
997  B3.x = 10.0 ; B3.y = 0.0;
998 
999  lw_dist2d_distpts_init(&dl, DIST_MIN);
1000  A1.x = -22.0; A1.y = 0.0;
1001  A2.x = -17.0; A2.y = -5.0;
1002  A3.x = -12.0; A3.y = 0.0;
1003  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1004  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1005  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 2.0, 0.000001);
1006 
1007  lw_dist2d_distpts_init(&dl, DIST_MIN);
1008  A1.x = -19.0; A1.y = 0.0;
1009  A2.x = -14.0; A2.y = -5.0;
1010  A3.x = - 9.0; A3.y = 0.0;
1011  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1012  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1013  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1014 
1015  lw_dist2d_distpts_init(&dl, DIST_MIN);
1016  A1.x = -9.0; A1.y = 0.0;
1017  A2.x = -4.0; A2.y = -5.0;
1018  A3.x = 1.0; A3.y = 0.0;
1019  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1020  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1021  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1022 
1023  lw_dist2d_distpts_init(&dl, DIST_MIN);
1024  A1.x = -1.0; A1.y = 0.0;
1025  A2.x = 4.0; A2.y = -5.0;
1026  A3.x = 9.0; A3.y = 0.0;
1027  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1028  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1029  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1030 
1031  lw_dist2d_distpts_init(&dl, DIST_MIN);
1032  A1.x = 1.0; A1.y = 0.0;
1033  A2.x = 6.0; A2.y = -5.0;
1034  A3.x = 11.0; A3.y = 0.0;
1035  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1036  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1037  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1038 
1039  lw_dist2d_distpts_init(&dl, DIST_MIN);
1040  A1.x = 11.0; A1.y = 0.0;
1041  A2.x = 16.0; A2.y = -5.0;
1042  A3.x = 21.0; A3.y = 0.0;
1043  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1044  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1045  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1046 
1047 
1048  lw_dist2d_distpts_init(&dl, DIST_MIN);
1049  A1.x = -15.0; A1.y = -6.0;
1050  A2.x = -10.0; A2.y = -1.0;
1051  A3.x = - 5.0; A3.y = -6.0;
1052  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1053  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1054  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1055 
1056  lw_dist2d_distpts_init(&dl, DIST_MIN);
1057  A1.x = -5.0; A1.y = 0.0;
1058  A2.x = 0.0; A2.y = 5.0;
1059  A3.x = 5.0; A3.y = 0.0;
1060  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1061  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1062  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
1063 
1064  lw_dist2d_distpts_init(&dl, DIST_MIN);
1065  A1.x = -5.0; A1.y = 0.0;
1066  A2.x = 0.0; A2.y = -5.0;
1067  A3.x = 5.0; A3.y = 0.0;
1068  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1069  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1070  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
1071 
1072 }
1073 
1074 static void
1075 test_lw_arc_length(void)
1076 {
1077 /* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
1078 
1079  POINT2D A1, A2, A3;
1080  double d;
1081 
1082  /* Unit semicircle at 0,0 */
1083  A1.x = -1; A1.y = 0;
1084  A2.x = 0 ; A2.y = 1;
1085  A3.x = 1 ; A3.y = 0;
1086 
1087  /* Arc above the unit semicircle */
1088  d = lw_arc_length(&A1, &A2, &A3);
1089  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1090  d = lw_arc_length(&A3, &A2, &A1);
1091  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1092 
1093  /* Unit semicircle at 0,0 */
1094  A1.x = 0; A1.y = 1;
1095  A2.x = 1; A2.y = 0;
1096  A3.x = 0; A3.y = -1;
1097 
1098  /* Arc to right of the unit semicircle */
1099  d = lw_arc_length(&A1, &A2, &A3);
1100  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1101  d = lw_arc_length(&A3, &A2, &A1);
1102  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1103 
1104  /* Unit 3/4 circle at 0,0 */
1105  A1.x = -1; A1.y = 0;
1106  A2.x = 1; A2.y = 0;
1107  A3.x = 0; A3.y = -1;
1108 
1109  /* Arc to right of the unit semicircle */
1110  d = lw_arc_length(&A1, &A2, &A3);
1111  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1112  d = lw_arc_length(&A3, &A2, &A1);
1113  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1114 }
1115 
1116 static void
1117 test_lw_dist2d_pt_ptarrayarc(void)
1118 {
1119  /* lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) */
1120  DISTPTS dl;
1121  int rv;
1122  LWLINE *lwline;
1123  POINT2D P;
1124 
1125  /* Unit semi-circle above X axis */
1126  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1127 
1128  /* Point at origin */
1129  P.x = P.y = 0;
1130  lw_dist2d_distpts_init(&dl, DIST_MIN);
1131  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1132  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1133  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1134 
1135  /* Point above arc on Y axis */
1136  P.y = 2;
1137  lw_dist2d_distpts_init(&dl, DIST_MIN);
1138  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1139  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1140  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1141 
1142  /* Point 45 degrees off arc, 2 radii from center */
1143  P.y = P.x = 2 * cos(M_PI_4);
1144  lw_dist2d_distpts_init(&dl, DIST_MIN);
1145  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1146  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1147  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1148 
1149  /* Four unit semi-circles surrounding the 2x2 box around origin */
1150  lwline_free(lwline);
1151  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
1152 
1153  /* Point at origin */
1154  P.x = P.y = 0;
1155  lw_dist2d_distpts_init(&dl, DIST_MIN);
1156  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1157  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1158  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
1159 
1160  /* Point on box edge */
1161  P.x = -1; P.y = 0;
1162  lw_dist2d_distpts_init(&dl, DIST_MIN);
1163  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1164  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1165  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1166 
1167  /* Point within a semicircle lobe */
1168  P.x = -1.5; P.y = 0;
1169  lw_dist2d_distpts_init(&dl, DIST_MIN);
1170  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1171  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1172  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1173 
1174  /* Point outside a semicircle lobe */
1175  P.x = -2.5; P.y = 0;
1176  lw_dist2d_distpts_init(&dl, DIST_MIN);
1177  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1178  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1179  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1180 
1181  /* Point outside a semicircle lobe */
1182  P.y = -2.5; P.x = 0;
1183  lw_dist2d_distpts_init(&dl, DIST_MIN);
1184  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1185  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1186  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1187 
1188  /* Point outside a semicircle lobe */
1189  P.y = 2; P.x = 1;
1190  lw_dist2d_distpts_init(&dl, DIST_MIN);
1191  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1192  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1193  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1.0, 0.000001);
1194 
1195  /* Clean up */
1196  lwline_free(lwline);
1197 }
1198 
1199 static void
1200 test_lw_dist2d_ptarray_ptarrayarc(void)
1201 {
1202  /* int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl) */
1203  DISTPTS dl;
1204  int rv;
1205  LWLINE *lwline1;
1206  LWLINE *lwline2;
1207 
1208  /* Unit semi-circle above X axis */
1209  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1210 
1211  /* Line above top of semi-circle */
1212  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
1213  lw_dist2d_distpts_init(&dl, DIST_MIN);
1214  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1215  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1216  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1217 
1218  /* Reversed arguments, should fail */
1219  lw_dist2d_distpts_init(&dl, DIST_MIN);
1220  cu_error_msg_reset();
1221  rv = lw_dist2d_ptarray_ptarrayarc(lwline1->points, lwline2->points, &dl);
1222  //printf("%s\n", cu_error_msg);
1223  CU_ASSERT_EQUAL( rv, LW_FAILURE );
1224  ASSERT_STRING_EQUAL(
1225  cu_error_msg,
1226  "lw_dist2d_ptarray_ptarrayarc called with non-arc input");
1227 
1228  lwline_free(lwline2);
1229 
1230  /* Line along side of semi-circle */
1231  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
1232  lw_dist2d_distpts_init(&dl, DIST_MIN);
1233  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1234  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1235  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1236 
1237  /* Four unit semi-circles surrounding the 2x2 box around origin */
1238  lwline_free(lwline1);
1239  lwline_free(lwline2);
1240  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
1241  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2.5 -3, -2.5 -2, -2.5 2, -2.5 3)"));
1242  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1243  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1244  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1245 
1246  lwline_free(lwline2);
1247  lwline_free(lwline1);
1248 }
1249 
1250 static void
1251 test_lwgeom_tcpa(void)
1252 {
1253  LWGEOM *g1, *g2;
1254  double m, dist;
1255 
1256  /* Invalid input, lack of dimensions */
1257 
1258  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1259  g2 = lwgeom_from_wkt("LINESTRING (0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1260  m = lwgeom_tcpa(g1, g2, NULL);
1261  lwgeom_free(g1);
1262  lwgeom_free(g2);
1263  ASSERT_DOUBLE_EQUAL(m, -1.0);
1264  ASSERT_STRING_EQUAL(
1265  cu_error_msg,
1266  "Both input geometries must have a measure dimension");
1267 
1268  /* Invalid input, not linestrings */
1269 
1270  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1271  g2 = lwgeom_from_wkt("POINT M (0 0 2)", LW_PARSER_CHECK_NONE);
1272  m = lwgeom_tcpa(g1, g2, NULL);
1273  lwgeom_free(g1);
1274  lwgeom_free(g2);
1275  ASSERT_DOUBLE_EQUAL(m, -1.0);
1276  ASSERT_STRING_EQUAL(cu_error_msg,
1277  "Both input geometries must be linestrings");
1278 
1279  /* Invalid input, too short linestring */
1280 
1281  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1282  g2 = lwgeom_from_wkt("LINESTRING M(2 0 1)", LW_PARSER_CHECK_NONE);
1283  dist = -77;
1284  m = lwgeom_tcpa(g1, g2, &dist);
1285  lwgeom_free(g1);
1286  lwgeom_free(g2);
1287  ASSERT_DOUBLE_EQUAL(dist, -77.0); /* not touched */
1288  ASSERT_DOUBLE_EQUAL(m, -1.0);
1289  ASSERT_STRING_EQUAL(
1290  cu_error_msg,
1291  "Both input lines must have at least 2 points" /* should be accepted
1292  ? */
1293 
1294  );
1295 
1296  /* Invalid input, empty linestring */
1297 
1298  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1299  g2 = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1300  m = lwgeom_tcpa(g1, g2, NULL);
1301  lwgeom_free(g1);
1302  lwgeom_free(g2);
1303  ASSERT_DOUBLE_EQUAL(m, -1.0);
1304  ASSERT_STRING_EQUAL(cu_error_msg,
1305  "Both input lines must have at least 2 points"
1306 
1307  );
1308 
1309  /* Timeranges do not overlap */
1310 
1311  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1312  g2 = lwgeom_from_wkt("LINESTRING M(0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1313  m = lwgeom_tcpa(g1, g2, NULL);
1314  lwgeom_free(g1);
1315  lwgeom_free(g2);
1316  ASSERT_DOUBLE_EQUAL(m, -2.0); /* means timeranges do not overlap */
1317 
1318  /* One of the tracks is still, the other passes to that point */
1319 
1320  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1321  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1322  dist = -1;
1323  m = lwgeom_tcpa(g1, g2, &dist);
1324  ASSERT_DOUBLE_EQUAL(m, 1.0);
1325  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1326  CU_ASSERT( lwgeom_cpa_within(g1, g2, 0.0) == LW_TRUE );
1327  lwgeom_free(g1);
1328  lwgeom_free(g2);
1329 
1330  /* One of the tracks is still, the other passes at 10 meters from point */
1331 
1332  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1333  g2 = lwgeom_from_wkt("LINESTRING M(-10 10 1, 10 10 5)", LW_PARSER_CHECK_NONE);
1334  dist = -1;
1335  m = lwgeom_tcpa(g1, g2, &dist);
1336  ASSERT_DOUBLE_EQUAL(m, 3.0);
1337  ASSERT_DOUBLE_EQUAL(dist, 10.0);
1338  CU_ASSERT( lwgeom_cpa_within(g1, g2, 11.0) == LW_TRUE );
1339  CU_ASSERT( lwgeom_cpa_within(g1, g2, 10.0) == LW_TRUE );
1340  CU_ASSERT( lwgeom_cpa_within(g1, g2, 9.0) == LW_FALSE );
1341  lwgeom_free(g1);
1342  lwgeom_free(g2);
1343 
1344  /* Equal tracks, 2d */
1345 
1346  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1347  g2 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1348  dist = -1;
1349  m = lwgeom_tcpa(g1, g2, &dist);
1350  lwgeom_free(g1);
1351  lwgeom_free(g2);
1352  ASSERT_DOUBLE_EQUAL(m, 10.0);
1353  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1354 
1355  /* Reversed tracks, 2d */
1356 
1357  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1358  g2 = lwgeom_from_wkt("LINESTRING M(10 0 10, 0 0 20)", LW_PARSER_CHECK_NONE);
1359  dist = -1;
1360  m = lwgeom_tcpa(g1, g2, &dist);
1361  lwgeom_free(g1);
1362  lwgeom_free(g2);
1363  ASSERT_DOUBLE_EQUAL(m, 15.0);
1364  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1365 
1366  /* Parallel tracks, same speed, 2d */
1367 
1368  g1 = lwgeom_from_wkt("LINESTRING M(2 0 10, 12 0 20)", LW_PARSER_CHECK_NONE);
1369  g2 = lwgeom_from_wkt("LINESTRING M(13 0 10, 23 0 20)", LW_PARSER_CHECK_NONE);
1370  dist = -1;
1371  m = lwgeom_tcpa(g1, g2, &dist);
1372  lwgeom_free(g1);
1373  lwgeom_free(g2);
1374  ASSERT_DOUBLE_EQUAL(m, 10.0);
1375  ASSERT_DOUBLE_EQUAL(dist, 11.0);
1376 
1377  /* Parallel tracks, different speed (g2 gets closer as time passes), 2d */
1378 
1379  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1380  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 9 0 20)", LW_PARSER_CHECK_NONE);
1381  dist = -1;
1382  m = lwgeom_tcpa(g1, g2, &dist);
1383  lwgeom_free(g1);
1384  lwgeom_free(g2);
1385  ASSERT_DOUBLE_EQUAL(m, 20.0);
1386  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1387 
1388  /* Parallel tracks, different speed (g2 left behind as time passes), 2d */
1389 
1390  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1391  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 6 0 20)", LW_PARSER_CHECK_NONE);
1392  dist = -1;
1393  m = lwgeom_tcpa(g1, g2, &dist);
1394  lwgeom_free(g1);
1395  lwgeom_free(g2);
1396  ASSERT_DOUBLE_EQUAL(m, 10.0);
1397  ASSERT_DOUBLE_EQUAL(dist, 2.0);
1398 
1399  /* Tracks, colliding, 2d */
1400 
1401  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1402  g2 = lwgeom_from_wkt("LINESTRING M(5 -8 0, 5 8 10)", LW_PARSER_CHECK_NONE);
1403  dist = -1;
1404  m = lwgeom_tcpa(g1, g2, &dist);
1405  lwgeom_free(g1);
1406  lwgeom_free(g2);
1407  ASSERT_DOUBLE_EQUAL(m, 5.0);
1408  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1409 
1410  /* Tracks crossing, NOT colliding, 2d */
1411 
1412  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1413  g2 = lwgeom_from_wkt("LINESTRING M(8 -5 0, 8 5 10)", LW_PARSER_CHECK_NONE);
1414  dist = -1;
1415  m = lwgeom_tcpa(g1, g2, &dist);
1416  lwgeom_free(g1);
1417  lwgeom_free(g2);
1418  ASSERT_DOUBLE_EQUAL(m, 6.5);
1419  ASSERT_DOUBLE_EQUAL(rint(dist*100), 212.0);
1420 
1421  /* Same origin, different direction, 2d */
1422 
1423  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1424  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, -100 0 10)", LW_PARSER_CHECK_NONE);
1425  dist = -1;
1426  m = lwgeom_tcpa(g1, g2, &dist);
1427  lwgeom_free(g1);
1428  lwgeom_free(g2);
1429  ASSERT_DOUBLE_EQUAL(m, 1.0);
1430  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1431 
1432  /* Same ending, different direction, 2d */
1433 
1434  g1 = lwgeom_from_wkt("LINESTRING M(10 0 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1435  g2 = lwgeom_from_wkt("LINESTRING M(0 -100 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1436  dist = -1;
1437  m = lwgeom_tcpa(g1, g2, &dist);
1438  lwgeom_free(g1);
1439  lwgeom_free(g2);
1440  ASSERT_DOUBLE_EQUAL(m, 10.0);
1441  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1442 
1443  /* Converging tracks, 3d */
1444 
1445  g1 = lwgeom_from_wkt("LINESTRING ZM(0 0 0 10, 10 0 0 20)", LW_PARSER_CHECK_NONE);
1446  g2 = lwgeom_from_wkt("LINESTRING ZM(0 0 8 10, 10 0 5 20)", LW_PARSER_CHECK_NONE);
1447  dist = -1;
1448  m = lwgeom_tcpa(g1, g2, &dist);
1449  lwgeom_free(g1);
1450  lwgeom_free(g2);
1451  ASSERT_DOUBLE_EQUAL(m, 20.0);
1452  ASSERT_DOUBLE_EQUAL(dist, 5.0);
1453 
1454  /* G1 stops at t=1 until t=4 to let G2 pass by, then continues */
1455  /* G2 passes at 1 meter from G1 t=3 */
1456 
1457  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 0 1 1, 0 1 4, 0 10 13)", LW_PARSER_CHECK_NONE);
1458  g2 = lwgeom_from_wkt("LINESTRING M(-10 2 0, 0 2 3, 12 2 13)", LW_PARSER_CHECK_NONE);
1459  dist = -1;
1460  m = lwgeom_tcpa(g1, g2, &dist);
1461  lwgeom_free(g1);
1462  lwgeom_free(g2);
1463  ASSERT_DOUBLE_EQUAL(m, 3.0);
1464  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1465 
1466  /* Test for https://trac.osgeo.org/postgis/ticket/3136 */
1467 
1468  g1 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,2 0 1432291536) ", LW_PARSER_CHECK_NONE);
1469  g2 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,1 0 1432291466.25,2 0 1432291500)", LW_PARSER_CHECK_NONE);
1470  dist = -1;
1471  m = lwgeom_tcpa(g1, g2, &dist);
1472  lwgeom_free(g1);
1473  lwgeom_free(g2);
1474  ASSERT_DOUBLE_EQUAL(m, 1432291464.0);
1475  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1476 
1477  /* Tracks share a single point in time */
1478 
1479  g1 = lwgeom_from_wkt("LINESTRINGM(0 0 0, 1 0 2)", LW_PARSER_CHECK_NONE);
1480  g2 = lwgeom_from_wkt("LINESTRINGM(0 0 2, 1 0 3)", LW_PARSER_CHECK_NONE);
1481  dist = -1;
1482  m = lwgeom_tcpa(g1, g2, &dist);
1483  lwgeom_free(g1);
1484  lwgeom_free(g2);
1485  ASSERT_DOUBLE_EQUAL(m, 2.0);
1486  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1487 
1488 }
1489 
1490 static void
1491 test_lwgeom_is_trajectory(void)
1492 {
1493  LWGEOM *g;
1494  int ret;
1495 
1496  g = lwgeom_from_wkt("POINT M(0 0 1)", LW_PARSER_CHECK_NONE);
1497  ret = lwgeom_is_trajectory(g);
1498  lwgeom_free(g);
1499  ASSERT_INT_EQUAL(ret, LW_FALSE); /* not a linestring */
1500 
1501  g = lwgeom_from_wkt("LINESTRING Z(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1502  ret = lwgeom_is_trajectory(g);
1503  lwgeom_free(g);
1504  ASSERT_INT_EQUAL(ret, LW_FALSE); /* no measure */
1505 
1506  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1507  ret = lwgeom_is_trajectory(g);
1508  lwgeom_free(g);
1509  ASSERT_INT_EQUAL(ret, LW_FALSE); /* same measure in two points */
1510 
1511  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 0)", LW_PARSER_CHECK_NONE);
1512  ret = lwgeom_is_trajectory(g);
1513  lwgeom_free(g);
1514  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1515 
1516  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 1 0 3, 2 2 2)", LW_PARSER_CHECK_NONE);
1517  ret = lwgeom_is_trajectory(g);
1518  lwgeom_free(g);
1519  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1520 
1521  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 2)", LW_PARSER_CHECK_NONE);
1522  ret = lwgeom_is_trajectory(g);
1523  lwgeom_free(g);
1524  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (still) */
1525 
1526  g = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1527  ret = lwgeom_is_trajectory(g);
1528  lwgeom_free(g);
1529  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (empty) */
1530 
1531  g = lwgeom_from_wkt("LINESTRING M (0 0 1)", LW_PARSER_CHECK_NONE);
1532  ret = lwgeom_is_trajectory(g);
1533  lwgeom_free(g);
1534  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (corner case) */
1535 }
1536 
1537 /*
1538 ** Used by test harness to register the tests in this file.
1539 */
1540 void measures_suite_setup(void);
1541 void measures_suite_setup(void)
1542 {
1543  CU_pSuite suite = CU_add_suite("measures", NULL, NULL);
1544  PG_ADD_TEST(suite, test_mindistance2d_tolerance);
1545  PG_ADD_TEST(suite, test_mindistance3d_tolerance);
1546  PG_ADD_TEST(suite, test_rect_tree_contains_point);
1547  PG_ADD_TEST(suite, test_rect_tree_intersects_tree);
1548  PG_ADD_TEST(suite, test_lwgeom_segmentize2d);
1549  PG_ADD_TEST(suite, test_lwgeom_locate_along);
1550  PG_ADD_TEST(suite, test_lw_dist2d_pt_arc);
1551  PG_ADD_TEST(suite, test_lw_dist2d_seg_arc);
1552  PG_ADD_TEST(suite, test_lw_dist2d_arc_arc);
1553  PG_ADD_TEST(suite, test_lw_arc_length);
1554  PG_ADD_TEST(suite, test_lw_dist2d_pt_ptarrayarc);
1555  PG_ADD_TEST(suite, test_lw_dist2d_ptarray_ptarrayarc);
1556  PG_ADD_TEST(suite, test_lwgeom_tcpa);
1557  PG_ADD_TEST(suite, test_lwgeom_is_trajectory);
1558  PG_ADD_TEST(suite, test_rect_tree_distance_tree);
1559 }
r
char * r
Definition: cu_in_wkt.c:24
test_lw_dist2d_ptarray_ptarrayarc
static void test_lw_dist2d_ptarray_ptarrayarc(void)
Definition: cu_measures.c:1200
test_lw_arc_length
static void test_lw_arc_length(void)
Definition: cu_measures.c:1075
test_rect_tree_intersects_tree
static void test_rect_tree_intersects_tree(void)
Definition: cu_measures.c:491
test_mindistance2d_tolerance
static void test_mindistance2d_tolerance(void)
Definition: cu_measures.c:87
test_lwgeom_is_trajectory
static void test_lwgeom_is_trajectory(void)
Definition: cu_measures.c:1491
measures_suite_setup
void measures_suite_setup(void)
Definition: cu_measures.c:1541
tree_inter
static int tree_inter(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:477
test_lwgeom_tcpa
static void test_lwgeom_tcpa(void)
Definition: cu_measures.c:1251
test_lw_dist2d_seg_arc
static void test_lw_dist2d_seg_arc(void)
Definition: cu_measures.c:807
test_rect_tree_distance_tree
static void test_rect_tree_distance_tree(void)
Definition: cu_measures.c:577
TDT
#define TDT(w1, w2, d)
Definition: cu_measures.c:574
test_lw_dist2d_pt_ptarrayarc
static void test_lw_dist2d_pt_ptarrayarc(void)
Definition: cu_measures.c:1117
test_lw_dist2d_arc_arc
static void test_lw_dist2d_arc_arc(void)
Definition: cu_measures.c:887
DIST3DTEST
#define DIST3DTEST(str1, str2, res, accepted_error)
Definition: cu_measures.c:36
test_rect_tree_contains_point
static void test_rect_tree_contains_point(void)
Definition: cu_measures.c:330
do_test_mindistance_tolerance
static void do_test_mindistance_tolerance(char *in1, char *in2, double expected_res, int line, double(*distancef)(const LWGEOM *, const LWGEOM *, double), double accepted_error)
Definition: cu_measures.c:41
test_lwgeom_locate_along
static void test_lwgeom_locate_along(void)
Definition: cu_measures.c:707
test_lwgeom_segmentize2d
static void test_lwgeom_segmentize2d(void)
Definition: cu_measures.c:638
lwgeom_from_text
static LWGEOM * lwgeom_from_text(const char *str)
Definition: cu_measures.c:25
test_lw_dist2d_pt_arc
static void test_lw_dist2d_pt_arc(void)
Definition: cu_measures.c:734
test_rect_tree_distance_tree_case
static double test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:552
test_mindistance3d_tolerance
static void test_mindistance3d_tolerance(void)
Definition: cu_measures.c:220
tree_pt
static int tree_pt(RECT_NODE *tree, double x, double y)
Definition: cu_measures.c:323
DIST2DTEST
#define DIST2DTEST(str1, str2, res, accepted_error)
Definition: cu_measures.c:33
cu_error_msg_reset
void cu_error_msg_reset()
Definition: liblwgeom/cunit/cu_tester.c:304
cu_error_msg
char cu_error_msg[MAX_CUNIT_ERROR_LENGTH+1]
Definition: liblwgeom/cunit/cu_tester.c:19
ASSERT_DOUBLE_EQUAL
#define ASSERT_DOUBLE_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:31
ASSERT_INT_EQUAL
#define ASSERT_INT_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:37
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:20
ASSERT_STRING_EQUAL
#define ASSERT_STRING_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:43
lwgeom_as_lwline
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:108
lwgeom_request_interrupt
void lwgeom_request_interrupt(void)
Request interruption of any running code.
Definition: lwgeom_api.c:670
lwgeom_locate_along
LWGEOM * lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
Determine the location(s) along a measured line where m occurs and return as a multipoint.
Definition: lwlinearreferencing.c:219
LW_FAILURE
#define LW_FAILURE
Definition: liblwgeom.h:110
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1138
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2060
lwgeom_segmentize2d
LWGEOM * lwgeom_segmentize2d(const LWGEOM *line, double dist)
Definition: lwgeom.c:753
LW_SUCCESS
#define LW_SUCCESS
Definition: liblwgeom.h:111
lwgeom_is_trajectory
int lwgeom_is_trajectory(const LWGEOM *geom)
Return LW_TRUE or LW_FALSE depending on whether or not a geometry is a linestring with measure value ...
Definition: lwgeom.c:2462
lwgeom_parse_wkt
int lwgeom_parse_wkt(LWGEOM_PARSER_RESULT *parser_result, char *wktstr, int parse_flags)
Parse a WKT geometry string into an LWGEOM structure.
Definition: lwin_wkt_parse.c:122
lwgeom_to_ewkt
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:547
lwgeom_tcpa
double lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
Find the time of closest point of approach.
Definition: lwlinearreferencing.c:1144
lwfree
void lwfree(void *mem)
Definition: lwutil.c:242
WKT_ISO
#define WKT_ISO
Definition: liblwgeom.h:2130
lwgeom_to_wkt
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:676
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:905
LW_TRUE
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:107
FLAGS_SET_SOLID
#define FLAGS_SET_SOLID(flags, value)
Definition: liblwgeom.h:191
lwline_free
void lwline_free(LWLINE *line)
Definition: lwline.c:67
lwgeom_cpa_within
int lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
Is the closest point of approach within a distance ?
Definition: lwlinearreferencing.c:1311
lw_arc_length
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns the length of a circular arc segment.
Definition: lwalgorithm.c:119
liblwgeom_internal.h
str
#define str(s)
Definition: lwgeom_functions_basic.c:43
INT32_MAX
#define INT32_MAX
Definition: lwin_wkt_lex.c:334
rect_tree_from_lwgeom
RECT_NODE * rect_tree_from_lwgeom(const LWGEOM *lwgeom)
Create a tree index on top an LWGEOM.
Definition: lwtree.c:861
rect_tree_contains_point
int rect_tree_contains_point(RECT_NODE *node, const POINT2D *pt)
Definition: lwtree.c:398
rect_tree_free
void rect_tree_free(RECT_NODE *node)
Recurse from top of node tree and free all children.
Definition: lwtree.c:69
distance
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
rect_tree_distance_tree
double rect_tree_distance_tree(RECT_NODE *n1, RECT_NODE *n2, double threshold)
Return the distance between two RECT_NODE trees.
Definition: lwtree.c:1354
rect_tree_intersects_tree
int rect_tree_intersects_tree(RECT_NODE *n1, RECT_NODE *n2)
Test if two RECT_NODE trees intersect one another.
Definition: lwtree.c:990
lw_dist2d_pt_arc
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1512
lw_dist2d_seg_arc
int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Calculate the shortest distance between an arc and an edge.
Definition: measures.c:1362
lw_dist2d_distpts_init
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:64
lw_dist2d_arc_arc
int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Definition: measures.c:1575
lw_dist2d_ptarray_ptarrayarc
int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
Test each segment of pa against each arc of pb for distance.
Definition: measures.c:1257
lw_dist2d_pt_ptarrayarc
int lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl)
Search all the arcs of pointarray to see which one is closest to p1 Returns minimum distance between ...
Definition: measures.c:1159
DIST_MIN
#define DIST_MIN
Definition: measures.h:44
cu_tester.h
DISTPTS::distance
double distance
Definition: measures.h:51
DISTPTS
Structure used in distance-calculations.
Definition: measures.h:50
LWGEOM::flags
lwflags_t flags
Definition: liblwgeom.h:447
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:469
POINT2D::y
double y
Definition: liblwgeom.h:376
POINT2D::x
double x
Definition: liblwgeom.h:376
rect_node
Definition: lwtree.h:66
struct_lwgeom_parser_result
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM.
Definition: liblwgeom.h:2068