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