PostGIS  3.1.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  ** "Fast path" case
213  */
214  DIST2DTEST("LINESTRING(10 0,11 1,12 2,13 3,14 4,15 5,16 6)",
215  "LINESTRING(1 1.5,2 3,3 4.5,4 6,5 7.5,6 9)",
216  8.3205, default_accepted_error);
217 
218  /*
219  ** Ticket 4326
220  */
221  DIST2DTEST(
222  "CURVEPOLYGON(CIRCULARSTRING(7874821 8715927,8907663 8715927,8844683 7750316,7937800 7750316,7874821 8715927))",
223  "POINT(5433865 8243495)", 2271704.2698450615, default_accepted_error);
224 }
225 
226 static void
227 test_mindistance3d_tolerance(void)
228 {
229  double default_accepted_error = 0.00001;
230  double zero_accepted_error = 0.0;
231  /* 2D [Z=0] should work just the same */
232  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5, default_accepted_error);
233  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5, default_accepted_error);
234  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0, default_accepted_error);
235  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))", 5.0, default_accepted_error);
236  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0))))", 5.0, default_accepted_error);
237  DIST3DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0, default_accepted_error);
238  DIST3DTEST("GEOMETRYCOLLECTION(POINT(0 0 0))", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0, default_accepted_error);
239  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0 0)))",
240  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))",
241  5.0, default_accepted_error);
242  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0 0)))",
243  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4 0)))",
244  5.0, default_accepted_error);
245  DIST3DTEST("LINESTRING(-2 0 0, -0.2 0 0)", "POINT(-2 0 0)", 0, zero_accepted_error);
246  DIST3DTEST("LINESTRING(-0.2 0 0, -2 0 0)", "POINT(-2 0 0)", 0, zero_accepted_error);
247  DIST3DTEST("LINESTRING(-1e-8 0 0, -0.2 0 0)", "POINT(-1e-8 0 0)", 0, zero_accepted_error);
248  DIST3DTEST("LINESTRING(-0.2 0 0, -1e-8 0 0)", "POINT(-1e-8 0 0)", 0, zero_accepted_error);
249 
250  /* Tests around intersections */
251  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);
252  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);
253  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);
254 
255  /* Same but triangles */
256  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);
257  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);
258  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);
259 
260  /* Triangle to triangle*/
261  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);
262 
263  /* Line in polygon */
264  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);
265 
266  /* Line has a point in the same plane as the polygon but isn't the closest*/
267  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);
268 
269  /* This is an invalid polygon since it defines just a line */
270  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);
271  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);
272  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);
273  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);
274 
275  /* A box in a box: two solids, one inside another */
276  DIST3DTEST(
277  "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)))",
278  "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)))",
279  0, zero_accepted_error);
280 
281  /* A box in a box with a hat: two solids, one inside another, Z ray up is hitting hat */
282  DIST3DTEST(
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  "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)))",
285  0, zero_accepted_error);
286 
287  /* Same but as TIN */
288  DIST3DTEST(
289  "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)))",
290  "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)))",
291  0, zero_accepted_error);
292 
293  /* Same but both are TIN */
294  DIST3DTEST(
295  "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)))",
296  "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)))",
297  0, zero_accepted_error);
298 
299  /* Point inside TIN */
300  DIST3DTEST(
301  "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)))",
302  "POINT(0 0 0)",
303  0, zero_accepted_error);
304 
305  /* A point hits vertical Z edge */
306  DIST3DTEST(
307  "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)))",
308  "POINT(0 0 0)",
309  0, zero_accepted_error);
310 
311  /* A point in the middle of a hole of extruded polygon */
312  DIST3DTEST(
313  "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)))",
314  "POINT(0 0 1)",
315  1, zero_accepted_error);
316 
317  /* A point at the face of a hole of extruded polygon */
318  DIST3DTEST(
319  "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)))",
320  "POINT(1 1 2)",
321  0, zero_accepted_error);
322 
323  /* A point at the face of a hole of extruded polygon */
324  DIST3DTEST(
325  "LINESTRING Z (-27974.1264 -110211.5032 148.9768,-27975.4229 -110210.9441 149.0093)",
326  "LINESTRING Z (-27995.4183 -110201.8041 149.3354,-27975.4229 -110210.9441 149.0093)",
327  0, zero_accepted_error);
328 }
329 
330 static int tree_pt(RECT_NODE *tree, double x, double y)
331 {
332  POINT2D pt;
333  pt.x = x; pt.y = y;
334  return rect_tree_contains_point(tree, &pt);
335 }
336 
337 static void test_rect_tree_contains_point(void)
338 {
339  LWGEOM *poly;
340  RECT_NODE* tree;
341 
342  /**********************************************************************
343  * curvepolygon
344  */
345  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);
346  tree = rect_tree_from_lwgeom(poly);
347  // char *wkt = rect_tree_to_wkt(tree);
348  // printf("%s\n", wkt);
349  // lwfree(wkt);
350  // return;
351 
352  /* in hole, within arc stroke */
353  CU_ASSERT_EQUAL(tree_pt(tree, 5, 7.5), 0);
354  /* inside */
355  CU_ASSERT_EQUAL(tree_pt(tree, 8, 9), 1);
356  /* outside */
357  CU_ASSERT_EQUAL(tree_pt(tree, -1, 5), 0);
358  /* outside */
359  CU_ASSERT_EQUAL(tree_pt(tree, -1, 7.5), 0);
360  /* outside, within arc stroke */
361  CU_ASSERT_EQUAL(tree_pt(tree, 0.2, 7.5), 0);
362  /* inside, within arc stroke */
363  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
364  /* inside, crossing arc stroke */
365  CU_ASSERT_EQUAL(tree_pt(tree, 2, 7.5), 1);
366  /* touching hole corner */
367  CU_ASSERT_EQUAL(tree_pt(tree, 7, 7), 1);
368 
369  rect_tree_free(tree);
370  lwgeom_free(poly);
371 
372 
373  /**********************************************************************
374  * polygon with hole and concavities
375  */
376  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);
377  tree = rect_tree_from_lwgeom(poly);
378 
379  /* inside, many grazings */
380  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
381  /* inside */
382  CU_ASSERT_EQUAL(tree_pt(tree, 3, 5.5), 1);
383  /* outside */
384  CU_ASSERT_EQUAL(tree_pt(tree, -3, 5.5), 0);
385  /* touching interior ring */
386  CU_ASSERT_EQUAL(tree_pt(tree, 4, 4), 1);
387  CU_ASSERT_EQUAL(tree_pt(tree, 6, 6), 1);
388  /* touching interior ring */
389  CU_ASSERT_EQUAL(tree_pt(tree, 4.5, 4), 1);
390  /* touching exterior ring */
391  CU_ASSERT_EQUAL(tree_pt(tree, 8, 0), 1);
392  CU_ASSERT_EQUAL(tree_pt(tree, 9, 0), 1);
393  CU_ASSERT_EQUAL(tree_pt(tree, 10, 1), 1);
394  CU_ASSERT_EQUAL(tree_pt(tree, 9.5, 1), 1);
395  CU_ASSERT_EQUAL(tree_pt(tree, 0, 10), 1);
396  /* touching grazing spike */
397  CU_ASSERT_EQUAL(tree_pt(tree, 1, 6), 1);
398  /* outide, many grazings */
399  CU_ASSERT_EQUAL(tree_pt(tree, -1, 6), 0);
400  /* within hole */
401  CU_ASSERT_EQUAL(tree_pt(tree, 5, 5), 0);
402  /* within */
403  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 4), 1);
404  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 6), 1);
405  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 9), 1);
406 
407  rect_tree_free(tree);
408  lwgeom_free(poly);
409 
410  /**********************************************************************
411  * square
412  */
413  poly = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
414  tree = rect_tree_from_lwgeom(poly);
415 
416  /* inside square */
417  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
418  /* outside square */
419  CU_ASSERT_EQUAL(tree_pt(tree, 1.5, 0.5), 0);
420  /* outside square grazing some edges */
421  CU_ASSERT_EQUAL(tree_pt(tree, -1, 1), 0);
422  /* inside square on corner */
423  CU_ASSERT_EQUAL(tree_pt(tree, 1, 1), 1);
424  /* inside square on top edge */
425  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 1), 1);
426  /* inside square on side edge */
427  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0.5), 1);
428 
429  rect_tree_free(tree);
430  lwgeom_free(poly);
431 
432  /* ziggy zaggy horizontal saw tooth polygon */
433  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);
434  tree = rect_tree_from_lwgeom(poly);
435 
436  /* not in, left side */
437  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 0.5), 0);
438  /* not in, right side */
439  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 0);
440  /* inside */
441  CU_ASSERT_EQUAL(tree_pt(tree, 2, 1), 1);
442  /* on left border */
443  CU_ASSERT_EQUAL(tree_pt(tree, 0, 1), 1);
444  /* on left border, grazing */
445  CU_ASSERT_EQUAL(tree_pt(tree, 0, 3), 1);
446  /* on right border */
447  CU_ASSERT_EQUAL(tree_pt(tree, 4, 0), 1);
448  /* on tooth concave */
449  CU_ASSERT_EQUAL(tree_pt(tree, 3, 3), 1);
450  /* on tooth convex */
451  CU_ASSERT_EQUAL(tree_pt(tree, 2, 0), 1);
452 
453  rect_tree_free(tree);
454  lwgeom_free(poly);
455 
456  /**********************************************************************
457  * ziggy zaggy vertical saw tooth polygon
458  */
459  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);
460  tree = rect_tree_from_lwgeom(poly);
461 
462  /* not in, left side */
463  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 3.5), 0);
464  /* not in, right side */
465  CU_ASSERT_EQUAL(tree_pt(tree, 6.0, 2.2), 0);
466  /* inside */
467  CU_ASSERT_EQUAL(tree_pt(tree, 3, 2), 1);
468  /* on bottom border */
469  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0), 1);
470  /* on top border */
471  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
472  /* on tooth concave */
473  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 1);
474  /* on tooth convex */
475  CU_ASSERT_EQUAL(tree_pt(tree, 0, 2), 1);
476  /* on tooth convex */
477  CU_ASSERT_EQUAL(tree_pt(tree, 0, 6), 1);
478 
479  rect_tree_free(tree);
480  lwgeom_free(poly);
481 
482 }
483 
484 static int tree_inter(const char *wkt1, const char *wkt2)
485 {
486  LWGEOM *g1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
487  LWGEOM *g2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
488  RECT_NODE *t1 = rect_tree_from_lwgeom(g1);
489  RECT_NODE *t2 = rect_tree_from_lwgeom(g2);
490  int result = rect_tree_intersects_tree(t1, t2);
491  rect_tree_free(t1);
492  rect_tree_free(t2);
493  lwgeom_free(g1);
494  lwgeom_free(g2);
495  return result;
496 }
497 
498 static void test_rect_tree_intersects_tree(void)
499 {
500  /* total overlap, A == B */
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 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))"),
504  LW_TRUE
505  );
506 
507  /* hiding between the tines of the comb */
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, 0.4 0.7, 0.3 0.7))"),
511  LW_FALSE
512  );
513 
514  /* between the tines, but with a corner overlapping */
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((0.3 0.7, 0.3 0.8, 0.4 0.8, 1.3 0.3, 0.3 0.7))"),
518  LW_TRUE
519  );
520 
521  /* Just touching the top left corner of the comb */
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  "POLYGON((-1 5, 0 5, 0 7, -1 7, -1 5))"),
525  LW_TRUE
526  );
527 
528  /* Contained, 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((1 2, 3 2)),POINT(1 2))"),
532  LW_TRUE
533  );
534 
535  /* 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(5 3))"),
539  LW_TRUE
540  );
541 
542  /* Not Touching, 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((6 3, 8 4)),POINT(1 3.5))"),
546  LW_FALSE
547  );
548 
549  /* Crossing, complex */
550  CU_ASSERT_EQUAL(tree_inter(
551  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
552  "GEOMETRYCOLLECTION(MULTILINESTRING((1.5 4.1, 1.6 2)),POINT(1 3.5))"),
553  LW_TRUE
554  );
555 }
556 
557 
558 static double
559 test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
560 {
561  LWGEOM *lw1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
562  LWGEOM *lw2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
563  RECT_NODE *n1 = rect_tree_from_lwgeom(lw1);
564  RECT_NODE *n2 = rect_tree_from_lwgeom(lw2);
565 
566  // rect_tree_printf(n1, 0);
567  // rect_tree_printf(n2, 0);
568  //
569  // printf("%s\n", rect_tree_to_wkt(n1));
570  // printf("%s\n", rect_tree_to_wkt(n2));
571 
572  double dist = rect_tree_distance_tree(n1, n2, 0.0);
573  // printf("%g\n", dist);
574  rect_tree_free(n1);
575  rect_tree_free(n2);
576  lwgeom_free(lw1);
577  lwgeom_free(lw2);
578  return dist;
579 }
580 
581 #define TDT(w1, w2, d) CU_ASSERT_DOUBLE_EQUAL(test_rect_tree_distance_tree_case(w1, w2), d, 0.00001);
582 
583 static void
584 test_rect_tree_distance_tree(void)
585 {
586  const char *wkt;
587 
588  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)))";
589  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);
590 
591  wkt = "POINT(0 0)";
592  TDT(wkt, "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
593  TDT(wkt, "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
594  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
595  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
596  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
597 
598  TDT("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
599  TDT("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
600  TDT("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
601  TDT("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
602 
603  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)))";
604  TDT(wkt, "POINT(3 14)", 1);
605  TDT(wkt, "POINT(3 8)", 0);
606  TDT(wkt, "POINT(6 5)", 1);
607  TDT(wkt, "POINT(6 4)", 0);
608 
609  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)))";
610  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
611  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
612  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
613  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
614  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
615  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
616  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
617  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
618 
619  wkt = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,5 0,0 0)))";
620  TDT(wkt, "POINT(3 0)", 0.0);
621  TDT(wkt, "POINT(5 0)", 0.0);
622  TDT(wkt, "POINT(7 0)", 2.0);
623  TDT(wkt, "POINT(2.5 3.5)", 1.0);
624 
625  wkt = "POINT(0 0)";
626  TDT(wkt, "POINT(0 1)", 1.0);
627  TDT(wkt, "POINT(1 0)", 1.0);
628 
629  wkt = "LINESTRING(0 0,1 0)";
630  TDT(wkt, "LINESTRING(1 0,1 1)", 0.0);
631  TDT(wkt, "LINESTRING(0 1,1 1)", 1.0);
632 
633  wkt = "POLYGON((0 0,0 1,1 1,1 0,0 0))";
634  TDT(wkt, "POINT(2 2)", sqrt(2));
635  TDT(wkt, "POINT(0.5 0.5)", 0);
636  TDT(wkt, "POINT(1 1)", 0);
637 
638  wkt = "POLYGON((0 0,0 10,10 10,10 0,0 0), (4 4,4 6,6 6,6 4,4 4))";
639  TDT(wkt, "POINT(5 5)", 1);
640  TDT(wkt, "POLYGON((5 5,5 5.5,5.5 5.5,5.5 5, 5 5))", 0.5);
641 }
642 
643 
644 static void
645 test_lwgeom_segmentize2d(void)
646 {
647  LWGEOM *linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
648  LWGEOM *lineout = lwgeom_segmentize2d(linein, 5);
649  char *strout = lwgeom_to_ewkt(lineout);
650  ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0)");
651  lwfree(strout);
652  lwgeom_free(linein);
653  lwgeom_free(lineout);
654 
655  /* test that segmentize is proportional - request every 6, get every 5 */
656  linein = lwgeom_from_wkt("LINESTRING(0 0, 20 0)", LW_PARSER_CHECK_NONE);
657  lineout = lwgeom_segmentize2d(linein, 6);
658  strout = lwgeom_to_ewkt(lineout);
659  ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0,15 0,20 0)");
660  lwfree(strout);
661  lwgeom_free(linein);
662  lwgeom_free(lineout);
663 
664  /* test too many segments */
665  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
666  lineout = lwgeom_segmentize2d(linein, 1e-100);
667  CU_ASSERT_EQUAL(lineout, NULL);
668  lwgeom_free(linein);
669 
670  /* test interruption */
671 
672  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
673  lwgeom_request_interrupt();
674  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
675  CU_ASSERT_EQUAL(lineout, NULL);
676  lwgeom_free(linein);
677 
678  linein = lwgeom_from_wkt("MULTILINESTRING((0 0,10 0),(20 0, 30 0))", LW_PARSER_CHECK_NONE);
679  lwgeom_request_interrupt();
680  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
681  CU_ASSERT_EQUAL(lineout, NULL);
682  lwgeom_free(linein);
683 
684  linein = lwgeom_from_wkt(
685  "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)))"
686  , LW_PARSER_CHECK_NONE);
687  lwgeom_request_interrupt();
688  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
689  CU_ASSERT_EQUAL(lineout, NULL);
690  lwgeom_free(linein);
691 
692  linein = lwgeom_from_wkt(
693  "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))"
694  , LW_PARSER_CHECK_NONE);
695  CU_ASSERT_FATAL(linein != NULL);
696  lwgeom_request_interrupt();
697  lineout = lwgeom_segmentize2d(linein, INT32_MAX);
698  CU_ASSERT_EQUAL(lineout, NULL);
699  lwgeom_free(linein);
700 
701  linein = lwgeom_from_wkt("LINESTRING(20 0, 30 0)", LW_PARSER_CHECK_NONE);
702  CU_ASSERT_FATAL(linein != NULL);
703  /* NOT INTERRUPTED */
704  lineout = lwgeom_segmentize2d(linein, 5);
705  CU_ASSERT_NOT_EQUAL_FATAL(lineout, NULL);
706  strout = lwgeom_to_ewkt(lineout);
707  ASSERT_STRING_EQUAL(strout, "LINESTRING(20 0,25 0,30 0)");
708  lwfree(strout);
709  lwgeom_free(linein);
710  lwgeom_free(lineout);
711 }
712 
713 static void
714 test_lwgeom_locate_along(void)
715 {
716  LWGEOM *geom = NULL;
717  LWGEOM *out = NULL;
718  double measure = 105.0;
719  char *str;
720 
721  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
722  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
723  out = lwgeom_locate_along(geom, measure, 0.0);
724  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
725  lwgeom_free(geom);
726  lwgeom_free(out);
727  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
728  lwfree(str);
729 
730  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
731  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);
732  out = lwgeom_locate_along(geom, measure, 0.0);
733  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
734  lwgeom_free(geom);
735  lwgeom_free(out);
736  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
737  lwfree(str);
738 }
739 
740 static void
741 test_lw_dist2d_pt_arc(void)
742 {
743  /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
744  DISTPTS dl;
745  POINT2D P, A1, A2, A3;
746  int rv;
747 
748 
749  /* Point within unit semicircle, 0.5 units from arc */
750  A1.x = -1; A1.y = 0;
751  A2.x = 0 ; A2.y = 1;
752  A3.x = 1 ; A3.y = 0;
753  P.x = 0 ; P.y = 0.5;
754 
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, 0.5 units from arc */
761  P.x = 0 ; P.y = 1.5;
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, 0.5, 0.000001);
766 
767  /* Point outside unit semicircle, sqrt(2) units from arc end point*/
768  P.x = 0 ; 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), 0.000001);
773 
774  /* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
775  P.x = 1 ; 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, sqrt(2.0)-1, 0.000001);
780 
781  /* Point on unit semicircle midpoint */
782  P.x = 0 ; P.y = 1;
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 unit semicircle endpoint */
789  P.x = 1 ; 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, 0, 0.000001);
794 
795  /* Point on semicircle center */
796  P.x = 0 ; P.y = 0;
797  lw_dist2d_distpts_init(&dl, DIST_MIN);
798  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
799  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
800  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
801 
802  /* Point inside closed circle */
803  P.x = 0 ; P.y = 0.5;
804  A2.x = 1; A2.y = 0;
805  A3 = A1;
806  lw_dist2d_distpts_init(&dl, DIST_MIN);
807  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
808  //printf("distance %g\n", dl.distance);
809  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
810  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
811 }
812 
813 static void
814 test_lw_dist2d_seg_arc(void)
815 {
816  /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */
817 
818  DISTPTS dl;
819  POINT2D A1, A2, B1, B2, B3;
820  int rv;
821 
822  /* Unit semicircle */
823  B1.x = -1; B1.y = 0;
824  B2.x = 0 ; B2.y = 1;
825  B3.x = 1 ; B3.y = 0;
826 
827  /* Edge above the unit semicircle */
828  lw_dist2d_distpts_init(&dl, DIST_MIN);
829  A1.x = -2; A1.y = 2;
830  A2.x = 2 ; A2.y = 2;
831  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
832  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
833  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
834 
835  /* Edge to the right of the unit semicircle */
836  lw_dist2d_distpts_init(&dl, DIST_MIN);
837  A1.x = 2; A1.y = -2;
838  A2.x = 2; A2.y = 2;
839  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
840  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
841  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
842 
843  /* Edge to the left of the unit semicircle */
844  lw_dist2d_distpts_init(&dl, DIST_MIN);
845  A1.x = -2; A1.y = -2;
846  A2.x = -2; A2.y = 2;
847  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
848  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
849  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
850 
851  /* Edge within the unit semicircle */
852  lw_dist2d_distpts_init(&dl, DIST_MIN);
853  A1.x = 0; A1.y = 0;
854  A2.x = 0; A2.y = 0.5;
855  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
856  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
857  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
858 
859  /* Edge grazing the unit semicircle */
860  lw_dist2d_distpts_init(&dl, DIST_MIN);
861  A1.x = -2; A1.y = 1;
862  A2.x = 2; A2.y = 1;
863  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
864  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
865  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);
866 
867  /* Line grazing the unit semicircle, but edge not */
868  lw_dist2d_distpts_init(&dl, DIST_MIN);
869  A1.x = 1; A1.y = 1;
870  A2.x = 2; A2.y = 1;
871  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
872  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
873  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
874 
875  /* Edge intersecting the unit semicircle */
876  lw_dist2d_distpts_init(&dl, DIST_MIN);
877  A1.x = 0; A1.y = 0;
878  A2.x = 2; A2.y = 2;
879  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
880  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
881  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
882 
883  /* Line intersecting the unit semicircle, but edge not */
884  lw_dist2d_distpts_init(&dl, DIST_MIN);
885  A1.x = -1; A1.y = 1;
886  A2.x = -2; A2.y = 2;
887  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
888  //printf("distance %g\n", dl.distance);
889  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
890  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
891 }
892 
893 static void
894 test_lw_dist2d_arc_arc(void)
895 {
896  /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
897  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
898  DISTPTS *dl) */
899  DISTPTS dl;
900  POINT2D A1, A2, A3, B1, B2, B3;
901  int rv;
902 
903  /* Ticket #4326 */
904  lw_dist2d_distpts_init(&dl, DIST_MIN);
905  A1.x = -1.0; A1.y = 4.0;
906  A2.x = 0.0; A2.y = 5.0;
907  A3.x = 1.0; A3.y = 4.0;
908  B1.x = 1.0; B1.y = 6.0;
909  B2.x = 6.0; B2.y = 1.0;
910  B3.x = 9.0; B3.y = 7.0;
911  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
912  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
913  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0475666, 0.000001);
914 
915  /* Unit semicircle at 0,0 */
916  B1.x = -1; B1.y = 0;
917  B2.x = 0 ; B2.y = 1;
918  B3.x = 1 ; B3.y = 0;
919 
920  /* Arc above the unit semicircle */
921  lw_dist2d_distpts_init(&dl, DIST_MIN);
922  A1.x = -1; A1.y = 3;
923  A2.x = 0 ; A2.y = 2;
924  A3.x = 1 ; A3.y = 3;
925  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
926  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
927  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
928 
929  /* Arc grazes the unit semicircle */
930  lw_dist2d_distpts_init(&dl, DIST_MIN);
931  A1.x = -1; A1.y = 2;
932  A2.x = 0 ; A2.y = 1;
933  A3.x = 1 ; A3.y = 2;
934  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
935  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
936  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
937 
938  /* Circles intersect, but arcs do not */
939  lw_dist2d_distpts_init(&dl, DIST_MIN);
940  A1.x = -1; A1.y = 1;
941  A2.x = 0; A2.y = 2;
942  A3.x = 1; A3.y = 1;
943  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
944  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
945  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);
946 
947  /* Circles and arcs intersect */
948  lw_dist2d_distpts_init(&dl, DIST_MIN);
949  A1.x = -1; A1.y = 1;
950  A2.x = 0; A2.y = 0;
951  A3.x = 1; A3.y = 1;
952  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
953  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
954  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
955 
956  /* Concentric: and fully parallel */
957  lw_dist2d_distpts_init(&dl, DIST_MIN);
958  A1.x = -2.0; A1.y = 0.0;
959  A2.x = 0.0; A2.y = 2.0;
960  A3.x = 2.0; A3.y = 0.0;
961  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
962  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
963  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
964 
965  /* Concentric: with A fully included in B's range */
966  lw_dist2d_distpts_init(&dl, DIST_MIN);
967  A1.x = -0.5 / sqrt(2.0); A1.y = 0.5 / sqrt(2.0);
968  A2.x = 0.0; A2.y = 0.5;
969  A3.x = 0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
970  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
971  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
972  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
973 
974  /* Concentric: with A partially included in B's range */
975  lw_dist2d_distpts_init(&dl, DIST_MIN);
976  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
977  A2.x = -0.5; A2.y = 0.0;
978  A3.x = -0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
979  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
980  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
981  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
982 
983  /* Concentric: with A and B without parallel segments */
984  lw_dist2d_distpts_init(&dl, DIST_MIN);
985  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
986  A2.x = 0.0; A2.y = -0.5;
987  A3.x = 0.5 / sqrt(2.0); A3.y = -0.5 / sqrt(2.0);
988  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
989  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
990  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.736813, 0.000001);
991 
992  /* Concentric: Arcs are the same */
993  lw_dist2d_distpts_init(&dl, DIST_MIN);
994  A1.x = -1.0; A1.y = 0.0;
995  A2.x = 0.0; A2.y = 1.0;
996  A3.x = 1.0; A3.y = 0.0;
997  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
998  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
999  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0, 0.000001);
1000 
1001  /* Check different orientations */
1002  B1.x = -10.0; B1.y = 0.0;
1003  B2.x = 0.0 ; B2.y = 10.0;
1004  B3.x = 10.0 ; B3.y = 0.0;
1005 
1006  lw_dist2d_distpts_init(&dl, DIST_MIN);
1007  A1.x = -22.0; A1.y = 0.0;
1008  A2.x = -17.0; A2.y = -5.0;
1009  A3.x = -12.0; A3.y = 0.0;
1010  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1011  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1012  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 2.0, 0.000001);
1013 
1014  lw_dist2d_distpts_init(&dl, DIST_MIN);
1015  A1.x = -19.0; A1.y = 0.0;
1016  A2.x = -14.0; A2.y = -5.0;
1017  A3.x = - 9.0; A3.y = 0.0;
1018  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1019  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1020  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1021 
1022  lw_dist2d_distpts_init(&dl, DIST_MIN);
1023  A1.x = -9.0; A1.y = 0.0;
1024  A2.x = -4.0; A2.y = -5.0;
1025  A3.x = 1.0; A3.y = 0.0;
1026  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1027  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1028  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1029 
1030  lw_dist2d_distpts_init(&dl, DIST_MIN);
1031  A1.x = -1.0; A1.y = 0.0;
1032  A2.x = 4.0; A2.y = -5.0;
1033  A3.x = 9.0; A3.y = 0.0;
1034  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1035  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1036  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1037 
1038  lw_dist2d_distpts_init(&dl, DIST_MIN);
1039  A1.x = 1.0; A1.y = 0.0;
1040  A2.x = 6.0; A2.y = -5.0;
1041  A3.x = 11.0; A3.y = 0.0;
1042  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1043  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1044  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1045 
1046  lw_dist2d_distpts_init(&dl, DIST_MIN);
1047  A1.x = 11.0; A1.y = 0.0;
1048  A2.x = 16.0; A2.y = -5.0;
1049  A3.x = 21.0; A3.y = 0.0;
1050  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1051  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1052  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1053 
1054 
1055  lw_dist2d_distpts_init(&dl, DIST_MIN);
1056  A1.x = -15.0; A1.y = -6.0;
1057  A2.x = -10.0; A2.y = -1.0;
1058  A3.x = - 5.0; A3.y = -6.0;
1059  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1060  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1061  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
1062 
1063  lw_dist2d_distpts_init(&dl, DIST_MIN);
1064  A1.x = -5.0; A1.y = 0.0;
1065  A2.x = 0.0; A2.y = 5.0;
1066  A3.x = 5.0; A3.y = 0.0;
1067  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1068  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1069  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
1070 
1071  lw_dist2d_distpts_init(&dl, DIST_MIN);
1072  A1.x = -5.0; A1.y = 0.0;
1073  A2.x = 0.0; A2.y = -5.0;
1074  A3.x = 5.0; A3.y = 0.0;
1075  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
1076  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1077  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
1078 
1079 }
1080 
1081 static void
1082 test_lw_arc_length(void)
1083 {
1084 /* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
1085 
1086  POINT2D A1, A2, A3;
1087  double d;
1088 
1089  /* Unit semicircle at 0,0 */
1090  A1.x = -1; A1.y = 0;
1091  A2.x = 0 ; A2.y = 1;
1092  A3.x = 1 ; A3.y = 0;
1093 
1094  /* Arc above the unit semicircle */
1095  d = lw_arc_length(&A1, &A2, &A3);
1096  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1097  d = lw_arc_length(&A3, &A2, &A1);
1098  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1099 
1100  /* Unit semicircle at 0,0 */
1101  A1.x = 0; A1.y = 1;
1102  A2.x = 1; A2.y = 0;
1103  A3.x = 0; A3.y = -1;
1104 
1105  /* Arc to right of the unit semicircle */
1106  d = lw_arc_length(&A1, &A2, &A3);
1107  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1108  d = lw_arc_length(&A3, &A2, &A1);
1109  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1110 
1111  /* Unit 3/4 circle at 0,0 */
1112  A1.x = -1; A1.y = 0;
1113  A2.x = 1; A2.y = 0;
1114  A3.x = 0; A3.y = -1;
1115 
1116  /* Arc to right of the unit semicircle */
1117  d = lw_arc_length(&A1, &A2, &A3);
1118  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1119  d = lw_arc_length(&A3, &A2, &A1);
1120  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1121 }
1122 
1123 static void
1124 test_lw_dist2d_pt_ptarrayarc(void)
1125 {
1126  /* lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) */
1127  DISTPTS dl;
1128  int rv;
1129  LWLINE *lwline;
1130  POINT2D P;
1131 
1132  /* Unit semi-circle above X axis */
1133  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1134 
1135  /* Point at origin */
1136  P.x = P.y = 0;
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 above arc on Y axis */
1143  P.y = 2;
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  /* Point 45 degrees off arc, 2 radii from center */
1150  P.y = P.x = 2 * cos(M_PI_4);
1151  lw_dist2d_distpts_init(&dl, DIST_MIN);
1152  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1153  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1154  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1155 
1156  /* Four unit semi-circles surrounding the 2x2 box around origin */
1157  lwline_free(lwline);
1158  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)"));
1159 
1160  /* Point at origin */
1161  P.x = 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, sqrt(2.0), 0.000001);
1166 
1167  /* Point on box edge */
1168  P.x = -1; 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, 1, 0.000001);
1173 
1174  /* Point within a semicircle lobe */
1175  P.x = -1.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.x = -2.5; P.y = 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.5; P.x = 0;
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, 0.5, 0.000001);
1194 
1195  /* Point outside a semicircle lobe */
1196  P.y = 2; P.x = 1;
1197  lw_dist2d_distpts_init(&dl, DIST_MIN);
1198  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1199  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1200  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1.0, 0.000001);
1201 
1202  /* Clean up */
1203  lwline_free(lwline);
1204 }
1205 
1206 static void
1207 test_lw_dist2d_ptarray_ptarrayarc(void)
1208 {
1209  /* int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl) */
1210  DISTPTS dl;
1211  int rv;
1212  LWLINE *lwline1;
1213  LWLINE *lwline2;
1214 
1215  /* Unit semi-circle above X axis */
1216  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1217 
1218  /* Line above top of semi-circle */
1219  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
1220  lw_dist2d_distpts_init(&dl, DIST_MIN);
1221  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1222  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1223  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1224 
1225  /* Reversed arguments, should fail */
1226  lw_dist2d_distpts_init(&dl, DIST_MIN);
1227  cu_error_msg_reset();
1228  rv = lw_dist2d_ptarray_ptarrayarc(lwline1->points, lwline2->points, &dl);
1229  //printf("%s\n", cu_error_msg);
1230  CU_ASSERT_EQUAL( rv, LW_FAILURE );
1231  ASSERT_STRING_EQUAL(
1232  cu_error_msg,
1233  "lw_dist2d_ptarray_ptarrayarc called with non-arc input");
1234 
1235  lwline_free(lwline2);
1236 
1237  /* Line along side of semi-circle */
1238  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
1239  lw_dist2d_distpts_init(&dl, DIST_MIN);
1240  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1241  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1242  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1243 
1244  /* Four unit semi-circles surrounding the 2x2 box around origin */
1245  lwline_free(lwline1);
1246  lwline_free(lwline2);
1247  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)"));
1248  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2.5 -3, -2.5 -2, -2.5 2, -2.5 3)"));
1249  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1250  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1251  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1252 
1253  lwline_free(lwline2);
1254  lwline_free(lwline1);
1255 }
1256 
1257 static void
1258 test_lwgeom_tcpa(void)
1259 {
1260  LWGEOM *g1, *g2;
1261  double m, dist;
1262 
1263  /* Invalid input, lack of dimensions */
1264 
1265  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1266  g2 = lwgeom_from_wkt("LINESTRING (0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1267  m = lwgeom_tcpa(g1, g2, NULL);
1268  lwgeom_free(g1);
1269  lwgeom_free(g2);
1270  ASSERT_DOUBLE_EQUAL(m, -1.0);
1271  ASSERT_STRING_EQUAL(
1272  cu_error_msg,
1273  "Both input geometries must have a measure dimension");
1274 
1275  /* Invalid input, not linestrings */
1276 
1277  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1278  g2 = lwgeom_from_wkt("POINT M (0 0 2)", LW_PARSER_CHECK_NONE);
1279  m = lwgeom_tcpa(g1, g2, NULL);
1280  lwgeom_free(g1);
1281  lwgeom_free(g2);
1282  ASSERT_DOUBLE_EQUAL(m, -1.0);
1283  ASSERT_STRING_EQUAL(cu_error_msg,
1284  "Both input geometries must be linestrings");
1285 
1286  /* Invalid input, too short linestring */
1287 
1288  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1289  g2 = lwgeom_from_wkt("LINESTRING M(2 0 1)", LW_PARSER_CHECK_NONE);
1290  dist = -77;
1291  m = lwgeom_tcpa(g1, g2, &dist);
1292  lwgeom_free(g1);
1293  lwgeom_free(g2);
1294  ASSERT_DOUBLE_EQUAL(dist, -77.0); /* not touched */
1295  ASSERT_DOUBLE_EQUAL(m, -1.0);
1296  ASSERT_STRING_EQUAL(
1297  cu_error_msg,
1298  "Both input lines must have at least 2 points" /* should be accepted
1299  ? */
1300 
1301  );
1302 
1303  /* Invalid input, empty linestring */
1304 
1305  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1306  g2 = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1307  m = lwgeom_tcpa(g1, g2, NULL);
1308  lwgeom_free(g1);
1309  lwgeom_free(g2);
1310  ASSERT_DOUBLE_EQUAL(m, -1.0);
1311  ASSERT_STRING_EQUAL(cu_error_msg,
1312  "Both input lines must have at least 2 points"
1313 
1314  );
1315 
1316  /* Timeranges do not overlap */
1317 
1318  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1319  g2 = lwgeom_from_wkt("LINESTRING M(0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1320  m = lwgeom_tcpa(g1, g2, NULL);
1321  lwgeom_free(g1);
1322  lwgeom_free(g2);
1323  ASSERT_DOUBLE_EQUAL(m, -2.0); /* means timeranges do not overlap */
1324 
1325  /* One of the tracks is still, the other passes to that point */
1326 
1327  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1328  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1329  dist = -1;
1330  m = lwgeom_tcpa(g1, g2, &dist);
1331  ASSERT_DOUBLE_EQUAL(m, 1.0);
1332  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1333  CU_ASSERT( lwgeom_cpa_within(g1, g2, 0.0) == LW_TRUE );
1334  lwgeom_free(g1);
1335  lwgeom_free(g2);
1336 
1337  /* One of the tracks is still, the other passes at 10 meters from point */
1338 
1339  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1340  g2 = lwgeom_from_wkt("LINESTRING M(-10 10 1, 10 10 5)", LW_PARSER_CHECK_NONE);
1341  dist = -1;
1342  m = lwgeom_tcpa(g1, g2, &dist);
1343  ASSERT_DOUBLE_EQUAL(m, 3.0);
1344  ASSERT_DOUBLE_EQUAL(dist, 10.0);
1345  CU_ASSERT( lwgeom_cpa_within(g1, g2, 11.0) == LW_TRUE );
1346  CU_ASSERT( lwgeom_cpa_within(g1, g2, 10.0) == LW_TRUE );
1347  CU_ASSERT( lwgeom_cpa_within(g1, g2, 9.0) == LW_FALSE );
1348  lwgeom_free(g1);
1349  lwgeom_free(g2);
1350 
1351  /* Equal tracks, 2d */
1352 
1353  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1354  g2 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1355  dist = -1;
1356  m = lwgeom_tcpa(g1, g2, &dist);
1357  lwgeom_free(g1);
1358  lwgeom_free(g2);
1359  ASSERT_DOUBLE_EQUAL(m, 10.0);
1360  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1361 
1362  /* Reversed tracks, 2d */
1363 
1364  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1365  g2 = lwgeom_from_wkt("LINESTRING M(10 0 10, 0 0 20)", LW_PARSER_CHECK_NONE);
1366  dist = -1;
1367  m = lwgeom_tcpa(g1, g2, &dist);
1368  lwgeom_free(g1);
1369  lwgeom_free(g2);
1370  ASSERT_DOUBLE_EQUAL(m, 15.0);
1371  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1372 
1373  /* Parallel tracks, same speed, 2d */
1374 
1375  g1 = lwgeom_from_wkt("LINESTRING M(2 0 10, 12 0 20)", LW_PARSER_CHECK_NONE);
1376  g2 = lwgeom_from_wkt("LINESTRING M(13 0 10, 23 0 20)", LW_PARSER_CHECK_NONE);
1377  dist = -1;
1378  m = lwgeom_tcpa(g1, g2, &dist);
1379  lwgeom_free(g1);
1380  lwgeom_free(g2);
1381  ASSERT_DOUBLE_EQUAL(m, 10.0);
1382  ASSERT_DOUBLE_EQUAL(dist, 11.0);
1383 
1384  /* Parallel tracks, different speed (g2 gets closer as time passes), 2d */
1385 
1386  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1387  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 9 0 20)", LW_PARSER_CHECK_NONE);
1388  dist = -1;
1389  m = lwgeom_tcpa(g1, g2, &dist);
1390  lwgeom_free(g1);
1391  lwgeom_free(g2);
1392  ASSERT_DOUBLE_EQUAL(m, 20.0);
1393  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1394 
1395  /* Parallel tracks, different speed (g2 left behind as time passes), 2d */
1396 
1397  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1398  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 6 0 20)", LW_PARSER_CHECK_NONE);
1399  dist = -1;
1400  m = lwgeom_tcpa(g1, g2, &dist);
1401  lwgeom_free(g1);
1402  lwgeom_free(g2);
1403  ASSERT_DOUBLE_EQUAL(m, 10.0);
1404  ASSERT_DOUBLE_EQUAL(dist, 2.0);
1405 
1406  /* Tracks, colliding, 2d */
1407 
1408  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1409  g2 = lwgeom_from_wkt("LINESTRING M(5 -8 0, 5 8 10)", LW_PARSER_CHECK_NONE);
1410  dist = -1;
1411  m = lwgeom_tcpa(g1, g2, &dist);
1412  lwgeom_free(g1);
1413  lwgeom_free(g2);
1414  ASSERT_DOUBLE_EQUAL(m, 5.0);
1415  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1416 
1417  /* Tracks crossing, NOT colliding, 2d */
1418 
1419  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1420  g2 = lwgeom_from_wkt("LINESTRING M(8 -5 0, 8 5 10)", LW_PARSER_CHECK_NONE);
1421  dist = -1;
1422  m = lwgeom_tcpa(g1, g2, &dist);
1423  lwgeom_free(g1);
1424  lwgeom_free(g2);
1425  ASSERT_DOUBLE_EQUAL(m, 6.5);
1426  ASSERT_DOUBLE_EQUAL(rint(dist*100), 212.0);
1427 
1428  /* Same origin, different direction, 2d */
1429 
1430  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1431  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, -100 0 10)", LW_PARSER_CHECK_NONE);
1432  dist = -1;
1433  m = lwgeom_tcpa(g1, g2, &dist);
1434  lwgeom_free(g1);
1435  lwgeom_free(g2);
1436  ASSERT_DOUBLE_EQUAL(m, 1.0);
1437  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1438 
1439  /* Same ending, different direction, 2d */
1440 
1441  g1 = lwgeom_from_wkt("LINESTRING M(10 0 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1442  g2 = lwgeom_from_wkt("LINESTRING M(0 -100 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1443  dist = -1;
1444  m = lwgeom_tcpa(g1, g2, &dist);
1445  lwgeom_free(g1);
1446  lwgeom_free(g2);
1447  ASSERT_DOUBLE_EQUAL(m, 10.0);
1448  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1449 
1450  /* Converging tracks, 3d */
1451 
1452  g1 = lwgeom_from_wkt("LINESTRING ZM(0 0 0 10, 10 0 0 20)", LW_PARSER_CHECK_NONE);
1453  g2 = lwgeom_from_wkt("LINESTRING ZM(0 0 8 10, 10 0 5 20)", LW_PARSER_CHECK_NONE);
1454  dist = -1;
1455  m = lwgeom_tcpa(g1, g2, &dist);
1456  lwgeom_free(g1);
1457  lwgeom_free(g2);
1458  ASSERT_DOUBLE_EQUAL(m, 20.0);
1459  ASSERT_DOUBLE_EQUAL(dist, 5.0);
1460 
1461  /* G1 stops at t=1 until t=4 to let G2 pass by, then continues */
1462  /* G2 passes at 1 meter from G1 t=3 */
1463 
1464  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 0 1 1, 0 1 4, 0 10 13)", LW_PARSER_CHECK_NONE);
1465  g2 = lwgeom_from_wkt("LINESTRING M(-10 2 0, 0 2 3, 12 2 13)", LW_PARSER_CHECK_NONE);
1466  dist = -1;
1467  m = lwgeom_tcpa(g1, g2, &dist);
1468  lwgeom_free(g1);
1469  lwgeom_free(g2);
1470  ASSERT_DOUBLE_EQUAL(m, 3.0);
1471  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1472 
1473  /* Test for https://trac.osgeo.org/postgis/ticket/3136 */
1474 
1475  g1 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,2 0 1432291536) ", LW_PARSER_CHECK_NONE);
1476  g2 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,1 0 1432291466.25,2 0 1432291500)", LW_PARSER_CHECK_NONE);
1477  dist = -1;
1478  m = lwgeom_tcpa(g1, g2, &dist);
1479  lwgeom_free(g1);
1480  lwgeom_free(g2);
1481  ASSERT_DOUBLE_EQUAL(m, 1432291464.0);
1482  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1483 
1484  /* Tracks share a single point in time */
1485 
1486  g1 = lwgeom_from_wkt("LINESTRINGM(0 0 0, 1 0 2)", LW_PARSER_CHECK_NONE);
1487  g2 = lwgeom_from_wkt("LINESTRINGM(0 0 2, 1 0 3)", LW_PARSER_CHECK_NONE);
1488  dist = -1;
1489  m = lwgeom_tcpa(g1, g2, &dist);
1490  lwgeom_free(g1);
1491  lwgeom_free(g2);
1492  ASSERT_DOUBLE_EQUAL(m, 2.0);
1493  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1494 
1495 }
1496 
1497 static void
1498 test_lwgeom_is_trajectory(void)
1499 {
1500  LWGEOM *g;
1501  int ret;
1502 
1503  g = lwgeom_from_wkt("POINT M(0 0 1)", LW_PARSER_CHECK_NONE);
1504  ret = lwgeom_is_trajectory(g);
1505  lwgeom_free(g);
1506  ASSERT_INT_EQUAL(ret, LW_FALSE); /* not a linestring */
1507 
1508  g = lwgeom_from_wkt("LINESTRING Z(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1509  ret = lwgeom_is_trajectory(g);
1510  lwgeom_free(g);
1511  ASSERT_INT_EQUAL(ret, LW_FALSE); /* no measure */
1512 
1513  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1514  ret = lwgeom_is_trajectory(g);
1515  lwgeom_free(g);
1516  ASSERT_INT_EQUAL(ret, LW_FALSE); /* same measure in two points */
1517 
1518  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 0)", LW_PARSER_CHECK_NONE);
1519  ret = lwgeom_is_trajectory(g);
1520  lwgeom_free(g);
1521  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1522 
1523  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 1 0 3, 2 2 2)", LW_PARSER_CHECK_NONE);
1524  ret = lwgeom_is_trajectory(g);
1525  lwgeom_free(g);
1526  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1527 
1528  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 2)", LW_PARSER_CHECK_NONE);
1529  ret = lwgeom_is_trajectory(g);
1530  lwgeom_free(g);
1531  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (still) */
1532 
1533  g = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1534  ret = lwgeom_is_trajectory(g);
1535  lwgeom_free(g);
1536  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (empty) */
1537 
1538  g = lwgeom_from_wkt("LINESTRING M (0 0 1)", LW_PARSER_CHECK_NONE);
1539  ret = lwgeom_is_trajectory(g);
1540  lwgeom_free(g);
1541  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (corner case) */
1542 }
1543 
1544 /*
1545 ** Used by test harness to register the tests in this file.
1546 */
1547 void measures_suite_setup(void);
1548 void measures_suite_setup(void)
1549 {
1550  CU_pSuite suite = CU_add_suite("measures", NULL, NULL);
1551  PG_ADD_TEST(suite, test_mindistance2d_tolerance);
1552  PG_ADD_TEST(suite, test_mindistance3d_tolerance);
1553  PG_ADD_TEST(suite, test_rect_tree_contains_point);
1554  PG_ADD_TEST(suite, test_rect_tree_intersects_tree);
1555  PG_ADD_TEST(suite, test_lwgeom_segmentize2d);
1556  PG_ADD_TEST(suite, test_lwgeom_locate_along);
1557  PG_ADD_TEST(suite, test_lw_dist2d_pt_arc);
1558  PG_ADD_TEST(suite, test_lw_dist2d_seg_arc);
1559  PG_ADD_TEST(suite, test_lw_dist2d_arc_arc);
1560  PG_ADD_TEST(suite, test_lw_arc_length);
1561  PG_ADD_TEST(suite, test_lw_dist2d_pt_ptarrayarc);
1562  PG_ADD_TEST(suite, test_lw_dist2d_ptarray_ptarrayarc);
1563  PG_ADD_TEST(suite, test_lwgeom_tcpa);
1564  PG_ADD_TEST(suite, test_lwgeom_is_trajectory);
1565  PG_ADD_TEST(suite, test_rect_tree_distance_tree);
1566 }
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:1207
test_lw_arc_length
static void test_lw_arc_length(void)
Definition: cu_measures.c:1082
test_rect_tree_intersects_tree
static void test_rect_tree_intersects_tree(void)
Definition: cu_measures.c:498
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:1498
measures_suite_setup
void measures_suite_setup(void)
Definition: cu_measures.c:1548
tree_inter
static int tree_inter(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:484
test_lwgeom_tcpa
static void test_lwgeom_tcpa(void)
Definition: cu_measures.c:1258
test_lw_dist2d_seg_arc
static void test_lw_dist2d_seg_arc(void)
Definition: cu_measures.c:814
test_rect_tree_distance_tree
static void test_rect_tree_distance_tree(void)
Definition: cu_measures.c:584
TDT
#define TDT(w1, w2, d)
Definition: cu_measures.c:581
test_lw_dist2d_pt_ptarrayarc
static void test_lw_dist2d_pt_ptarrayarc(void)
Definition: cu_measures.c:1124
test_lw_dist2d_arc_arc
static void test_lw_dist2d_arc_arc(void)
Definition: cu_measures.c:894
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:337
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:714
test_lwgeom_segmentize2d
static void test_lwgeom_segmentize2d(void)
Definition: cu_measures.c:645
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:741
test_rect_tree_distance_tree_case
static double test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:559
test_mindistance3d_tolerance
static void test_mindistance3d_tolerance(void)
Definition: cu_measures.c:227
tree_pt
static int tree_pt(RECT_NODE *tree, double x, double y)
Definition: cu_measures.c:330
DIST2DTEST
#define DIST2DTEST(str1, str2, res, accepted_error)
Definition: cu_measures.c:33
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
cu_error_msg_reset
void cu_error_msg_reset()
Definition: liblwgeom/cunit/cu_tester.c:306
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:32
ASSERT_INT_EQUAL
#define ASSERT_INT_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:38
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:21
ASSERT_STRING_EQUAL
#define ASSERT_STRING_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:62
lwgeom_as_lwline
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:162
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:664
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:2085
lwgeom_segmentize2d
LWGEOM * lwgeom_segmentize2d(const LWGEOM *line, double dist)
Definition: lwgeom.c:754
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:2473
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:548
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:2155
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:704
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:44
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:475
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:497
POINT2D::y
double y
Definition: liblwgeom.h:404
POINT2D::x
double x
Definition: liblwgeom.h:404
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:2093