PostGIS  2.1.10dev-r@@SVN_REVISION@@
cu_measures.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: cu_measures.c 13242 2015-02-19 00:22:22Z pramsey $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include "CUnit/Basic.h"
17 
18 #include "liblwgeom_internal.h"
19 #include "cu_tester.h"
20 #include "measures.h"
21 #include "lwtree.h"
22 
23 static LWGEOM* lwgeom_from_text(const char *str)
24 {
26  if( LW_FAILURE == lwgeom_parse_wkt(&r, (char*)str, LW_PARSER_CHECK_NONE) )
27  return NULL;
28  return r.geom;
29 }
30 
31 #define DIST2DTEST(str1, str2, res) do_test_mindistance2d_tolerance(str1, str2, res, __LINE__)
32 
33 static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line)
34 {
35  LWGEOM *lw1;
36  LWGEOM *lw2;
37  double distance;
38  char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n";
39  char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n %s\n\n";
40 
43 
44  if ( ! lw1 )
45  {
46  printf(msg2, in1);
47  exit(1);
48  }
49  if ( ! lw2 )
50  {
51  printf(msg2, in2);
52  exit(1);
53  }
54 
55  distance = lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0);
56  lwgeom_free(lw1);
57  lwgeom_free(lw2);
58 
59  if ( fabs(distance - expected_res) > 0.00001 )
60  {
61  printf(msg1, distance, expected_res, line);
62  CU_FAIL();
63  }
64  else
65  {
66  CU_PASS();
67  }
68 
69 }
70 
72 {
73  /*
74  ** Simple case.
75  */
76  DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
77 
78  /*
79  ** Point vs Geometry Collection.
80  */
81  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
82 
83  /*
84  ** Point vs Geometry Collection Collection.
85  */
86  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
87 
88  /*
89  ** Point vs Geometry Collection Collection Collection.
90  */
91  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
92 
93  /*
94  ** Point vs Geometry Collection Collection Collection Multipoint.
95  */
96  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
97 
98  /*
99  ** Geometry Collection vs Geometry Collection
100  */
101  DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
102 
103  /*
104  ** Geometry Collection Collection vs Geometry Collection Collection
105  */
106  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
107 
108  /*
109  ** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint
110  */
111  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0);
112 
113  /*
114  ** Linestring vs its start point
115  */
116  DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
117 
118  /*
119  ** Linestring vs its end point
120  */
121  DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
122 
123  /*
124  ** Linestring vs its start point (tricky number, see #1459)
125  */
126  DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
127 
128  /*
129  ** Linestring vs its end point (tricky number, see #1459)
130  */
131  DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
132 
133  /*
134  * Circular string and point
135  */
136  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1);
137  DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1);
138 
139  /*
140  * Circular string and Circular string
141  */
142  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1);
143 
144  /*
145  * CurvePolygon and Point
146  */
147  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)))";
148  DIST2DTEST(cs1, "POINT(3 14)", 1);
149  DIST2DTEST(cs1, "POINT(3 8)", 0);
150  DIST2DTEST(cs1, "POINT(6 5)", 1);
151  DIST2DTEST(cs1, "POINT(6 4)", 0);
152 
153  /*
154  * CurvePolygon and Linestring
155  */
156  DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484);
157  DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0);
158  DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0);
159  DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786);
160  DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913);
161 
162  /*
163  * CurvePolygon and Polygon
164  */
165  DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415);
166  DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0);
167  DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0);
168 
169  /*
170  * CurvePolygon and CurvePolygon
171  */
172  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666);
173  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0);
174 
175  /*
176  * MultiSurface and CurvePolygon
177  */
178  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)))";
179  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
180  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
181  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
182  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
183  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
184  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
185  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
186  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
187 
188  /*
189  * MultiCurve and Linestring
190  */
191  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);
192 
193 }
194 
196 {
197  LWPOLY *poly;
198  POINT2D p;
199  RECT_NODE* tree;
200  int result;
201  int boundary = 0;
202 
203  /* square */
204  poly = (LWPOLY*)lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
205  tree = rect_tree_new(poly->rings[0]);
206 
207  /* inside square */
208  boundary = 0;
209  p.x = 0.5;
210  p.y = 0.5;
211  result = rect_tree_contains_point(tree, &p, &boundary);
212  CU_ASSERT_NOT_EQUAL(result, 0);
213 
214  /* outside square */
215  boundary = 0;
216  p.x = 1.5;
217  p.y = 0.5;
218  result = rect_tree_contains_point(tree, &p, &boundary);
219  CU_ASSERT_EQUAL(result, 0);
220 
221  rect_tree_free(tree);
222  lwpoly_free(poly);
223 
224  /* ziggy zaggy horizontal saw tooth polygon */
225  poly = (LWPOLY*)lwgeom_from_wkt("POLYGON((0 0, 1 3, 2 0, 3 3, 4 0, 4 5, 0 5, 0 0))", LW_PARSER_CHECK_NONE);
226  tree = rect_tree_new(poly->rings[0]);
227 
228  /* not in, left side */
229  boundary = 0;
230  p.x = -0.5;
231  p.y = 0.5;
232  result = rect_tree_contains_point(tree, &p, &boundary);
233  CU_ASSERT_EQUAL(result, 0);
234 
235  /* not in, right side */
236  boundary = 0;
237  p.x = 3.0;
238  p.y = 1.0;
239  result = rect_tree_contains_point(tree, &p, &boundary);
240  CU_ASSERT_EQUAL(result, 0);
241 
242  /* inside */
243  boundary = 0;
244  p.x = 2.0;
245  p.y = 1.0;
246  result = rect_tree_contains_point(tree, &p, &boundary);
247  CU_ASSERT_NOT_EQUAL(result, 0);
248 
249  /* on left border */
250  boundary = 0;
251  p.x = 0.0;
252  p.y = 1.0;
253  result = rect_tree_contains_point(tree, &p, &boundary);
254  CU_ASSERT_EQUAL(boundary, 1);
255 
256  /* on right border */
257  boundary = 0;
258  p.x = 4.0;
259  p.y = 0.0;
260  result = rect_tree_contains_point(tree, &p, &boundary);
261  CU_ASSERT_EQUAL(boundary, 1);
262 
263  /* on tooth concave */
264  boundary = 0;
265  p.x = 3.0;
266  p.y = 3.0;
267  result = rect_tree_contains_point(tree, &p, &boundary);
268  CU_ASSERT_EQUAL(boundary, 1);
269 
270  /* on tooth convex */
271  boundary = 0;
272  p.x = 2.0;
273  p.y = 0.0;
274  result = rect_tree_contains_point(tree, &p, &boundary);
275  CU_ASSERT_EQUAL(boundary, 1);
276 
277  rect_tree_free(tree);
278  lwpoly_free(poly);
279 
280  /* ziggy zaggy vertical saw tooth polygon */
281  poly = (LWPOLY*)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);
282  tree = rect_tree_new(poly->rings[0]);
283 
284  /* not in, left side */
285  boundary = 0;
286  p.x = -0.5;
287  p.y = 3.5;
288  result = rect_tree_contains_point(tree, &p, &boundary);
289  CU_ASSERT_EQUAL(result, 0);
290 
291  /* not in, right side */
292  boundary = 0;
293  p.x = 6.0;
294  p.y = 2.2;
295  result = rect_tree_contains_point(tree, &p, &boundary);
296  CU_ASSERT_EQUAL(result, 0);
297 
298  /* inside */
299  boundary = 0;
300  p.x = 3.0;
301  p.y = 2.0;
302  result = rect_tree_contains_point(tree, &p, &boundary);
303  CU_ASSERT_NOT_EQUAL(result, 0);
304 
305  /* on bottom border */
306  boundary = 0;
307  p.x = 1.0;
308  p.y = 0.0;
309  result = rect_tree_contains_point(tree, &p, &boundary);
310  CU_ASSERT_EQUAL(boundary, 1);
311 
312  /* on top border */
313  boundary = 0;
314  p.x = 3.0;
315  p.y = 6.0;
316  result = rect_tree_contains_point(tree, &p, &boundary);
317  CU_ASSERT_EQUAL(boundary, 1);
318 
319  /* on tooth concave */
320  boundary = 0;
321  p.x = 3.0;
322  p.y = 1.0;
323  result = rect_tree_contains_point(tree, &p, &boundary);
324  CU_ASSERT_EQUAL(boundary, 1);
325 
326  /* on tooth convex */
327  boundary = 0;
328  p.x = 0.0;
329  p.y = 2.0;
330  result = rect_tree_contains_point(tree, &p, &boundary);
331  CU_ASSERT_EQUAL(boundary, 1);
332 
333  /* on tooth convex */
334  boundary = 0;
335  p.x = 0.0;
336  p.y = 6.0;
337  result = rect_tree_contains_point(tree, &p, &boundary);
338  CU_ASSERT_EQUAL(boundary, 1);
339 
340  rect_tree_free(tree);
341  lwpoly_free(poly);
342 
343 }
344 
346 {
347  LWPOLY *poly1, *poly2;
348  RECT_NODE *tree1, *tree2;
349  int result;
350 
351  /* total overlap, A == B */
352  poly1 = (LWPOLY*)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);
353  poly2 = (LWPOLY*)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);
354  tree1 = rect_tree_new(poly1->rings[0]);
355  tree2 = rect_tree_new(poly2->rings[0]);
356  result = rect_tree_intersects_tree(tree1, tree2);
357  CU_ASSERT_EQUAL(result, LW_TRUE);
358  lwpoly_free(poly1);
359  lwpoly_free(poly2);
360  rect_tree_free(tree1);
361  rect_tree_free(tree2);
362 
363  /* hiding between the tines of the comb */
364  poly1 = (LWPOLY*)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);
365  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 0.4 0.7, 0.3 0.7))", LW_PARSER_CHECK_NONE);
366  tree1 = rect_tree_new(poly1->rings[0]);
367  tree2 = rect_tree_new(poly2->rings[0]);
368  result = rect_tree_intersects_tree(tree1, tree2);
369  CU_ASSERT_EQUAL(result, LW_FALSE);
370  lwpoly_free(poly1);
371  lwpoly_free(poly2);
372  rect_tree_free(tree1);
373  rect_tree_free(tree2);
374 
375  /* between the tines, but with a corner overlapping */
376  poly1 = (LWPOLY*)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);
377  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 1.3 0.3, 0.3 0.7))", LW_PARSER_CHECK_NONE);
378  tree1 = rect_tree_new(poly1->rings[0]);
379  tree2 = rect_tree_new(poly2->rings[0]);
380  result = rect_tree_intersects_tree(tree1, tree2);
381  CU_ASSERT_EQUAL(result, LW_TRUE);
382  lwpoly_free(poly1);
383  lwpoly_free(poly2);
384  rect_tree_free(tree1);
385  rect_tree_free(tree2);
386 
387  /* Just touching the top left corner of the comb */
388  poly1 = (LWPOLY*)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);
389  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((-1 5, 0 5, 0 7, -1 7, -1 5))", LW_PARSER_CHECK_NONE);
390  tree1 = rect_tree_new(poly1->rings[0]);
391  tree2 = rect_tree_new(poly2->rings[0]);
392  result = rect_tree_intersects_tree(tree1, tree2);
393  CU_ASSERT_EQUAL(result, LW_TRUE);
394  lwpoly_free(poly1);
395  lwpoly_free(poly2);
396  rect_tree_free(tree1);
397  rect_tree_free(tree2);
398 
399 
400 }
401 
402 static void
404 {
405  LWGEOM *linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
406  LWGEOM *lineout = lwgeom_segmentize2d(linein, 5);
407  char *strout = lwgeom_to_ewkt(lineout);
408  CU_ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0)");
409  lwfree(strout);
410  lwgeom_free(linein);
411  lwgeom_free(lineout);
412 }
413 
414 static void
416 {
417  LWGEOM *geom = NULL;
418  LWGEOM *out = NULL;
419  double measure = 105.0;
420  char *str;
421 
422  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
423  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
424  out = lwgeom_locate_along(geom, measure, 0.0);
425  str = lwgeom_to_wkt(out, WKT_ISO, 8, NULL);
426  lwgeom_free(geom);
427  lwgeom_free(out);
428  CU_ASSERT_STRING_EQUAL("MULTIPOINT M (55.226131 55.226131 105)", str);
429  lwfree(str);
430 
431  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
432  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);
433  out = lwgeom_locate_along(geom, measure, 0.0);
434  str = lwgeom_to_wkt(out, WKT_ISO, 8, NULL);
435  lwgeom_free(geom);
436  lwgeom_free(out);
437  CU_ASSERT_STRING_EQUAL("MULTIPOINT M (55.226131 55.226131 105)", str);
438  lwfree(str);
439 }
440 
441 static void
443 {
444  /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
445  DISTPTS dl;
446  POINT2D P, A1, A2, A3;
447  int rv;
448 
449 
450  /* Point within unit semicircle, 0.5 units from arc */
451  A1.x = -1; A1.y = 0;
452  A2.x = 0 ; A2.y = 1;
453  A3.x = 1 ; A3.y = 0;
454  P.x = 0 ; P.y = 0.5;
455 
457  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
458  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
459 
460  /* Point outside unit semicircle, 0.5 units from arc */
461  P.x = 0 ; P.y = 1.5;
463  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
464  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
465 
466  /* Point outside unit semicircle, sqrt(2) units from arc end point*/
467  P.x = 0 ; P.y = -1;
469  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
470  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
471 
472  /* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
473  P.x = 1 ; P.y = 1;
475  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
476  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
477 
478  /* Point on unit semicircle midpoint */
479  P.x = 0 ; P.y = 1;
481  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
482  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
483 
484  /* Point on unit semicircle endpoint */
485  P.x = 1 ; P.y = 0;
487  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
488  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
489 
490  /* Point inside closed circle */
491  P.x = 0 ; P.y = 0.5;
492  A2.x = 1; A2.y = 0;
493  A3 = A1;
495  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
496  //printf("distance %g\n", dl.distance);
497  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
498 }
499 
500 static void
502 {
503  /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */
504 
505  DISTPTS dl;
506  POINT2D A1, A2, B1, B2, B3;
507  int rv;
508 
509  /* Unit semicircle */
510  B1.x = -1; B1.y = 0;
511  B2.x = 0 ; B2.y = 1;
512  B3.x = 1 ; B3.y = 0;
513 
514  /* Edge above the unit semicircle */
516  A1.x = -2; A1.y = 2;
517  A2.x = 2 ; A2.y = 2;
518  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
519  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
520 
521  /* Edge to the right of the unit semicircle */
523  A1.x = 2; A1.y = -2;
524  A2.x = 2; A2.y = 2;
525  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
526  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
527 
528  /* Edge to the left of the unit semicircle */
530  A1.x = -2; A1.y = -2;
531  A2.x = -2; A2.y = 2;
532  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
533  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
534 
535  /* Edge within the unit semicircle */
537  A1.x = 0; A1.y = 0;
538  A2.x = 0; A2.y = 0.5;
539  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
540  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
541 
542  /* Edge grazing the unit semicircle */
544  A1.x = -2; A1.y = 1;
545  A2.x = 2; A2.y = 1;
546  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
547  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);
548 
549  /* Line grazing the unit semicircle, but edge not */
551  A1.x = 1; A1.y = 1;
552  A2.x = 2; A2.y = 1;
553  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
554  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
555 
556  /* Edge intersecting the unit semicircle */
558  A1.x = 0; A1.y = 0;
559  A2.x = 2; A2.y = 2;
560  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
561  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
562 
563  /* Line intersecting the unit semicircle, but edge not */
565  A1.x = -1; A1.y = 1;
566  A2.x = -2; A2.y = 2;
567  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
568  //printf("distance %g\n", dl.distance);
569  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
570 }
571 
572 static void
574 {
575  /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
576  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
577  DISTPTS *dl) */
578  DISTPTS dl;
579  POINT2D A1, A2, A3, B1, B2, B3;
580  int rv;
581 
582  /* Unit semicircle at 0,0 */
583  B1.x = -1; B1.y = 0;
584  B2.x = 0 ; B2.y = 1;
585  B3.x = 1 ; B3.y = 0;
586 
587  /* Arc above the unit semicircle */
589  A1.x = -1; A1.y = 3;
590  A2.x = 0 ; A2.y = 2;
591  A3.x = 1 ; A3.y = 3;
592  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
593  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
594 
595  /* Arc grazes the unit semicircle */
597  A1.x = -1; A1.y = 2;
598  A2.x = 0 ; A2.y = 1;
599  A3.x = 1 ; A3.y = 2;
600  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
601  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
602 
603  /* Circles intersect, but arcs do not */
605  A1.x = -1; A1.y = 1;
606  A2.x = 0; A2.y = 2;
607  A3.x = 1; A3.y = 1;
608  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
609  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);
610 
611  /* Circles and arcs intersect */
613  A1.x = -1; A1.y = 1;
614  A2.x = 0; A2.y = 0;
615  A3.x = 1; A3.y = 1;
616  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
617  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
618 
619  /* inscribed and closest on arcs */
621  A1.x = -0.5; A1.y = 0.0;
622  A2.x = 0.0; A2.y = 0.5;
623  A3.x = 0.5; A3.y = 0.0;
624  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
625  //printf("distance %g\n", dl.distance);
626  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
627 
628  /* inscribed and closest not on arcs */
630  A1.x = -0.5; A1.y = 0.0;
631  A2.x = 0.0; A2.y = -0.5;
632  A3.x = 0.5; A3.y = 0.0;
633  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
634  //printf("distance %g\n", dl.distance);
635  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
636 }
637 
638 static void
640 {
641 /* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
642 
643  POINT2D A1, A2, A3;
644  double d;
645 
646  /* Unit semicircle at 0,0 */
647  A1.x = -1; A1.y = 0;
648  A2.x = 0 ; A2.y = 1;
649  A3.x = 1 ; A3.y = 0;
650 
651  /* Arc above the unit semicircle */
652  d = lw_arc_length(&A1, &A2, &A3);
653  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
654  d = lw_arc_length(&A3, &A2, &A1);
655  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
656 
657  /* Unit semicircle at 0,0 */
658  A1.x = 0; A1.y = 1;
659  A2.x = 1; A2.y = 0;
660  A3.x = 0; A3.y = -1;
661 
662  /* Arc to right of the unit semicircle */
663  d = lw_arc_length(&A1, &A2, &A3);
664  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
665  d = lw_arc_length(&A3, &A2, &A1);
666  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
667 
668  /* Unit 3/4 circle at 0,0 */
669  A1.x = -1; A1.y = 0;
670  A2.x = 1; A2.y = 0;
671  A3.x = 0; A3.y = -1;
672 
673  /* Arc to right of the unit semicircle */
674  d = lw_arc_length(&A1, &A2, &A3);
675  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);
676  d = lw_arc_length(&A3, &A2, &A1);
677  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);
678 }
679 
680 static void
682 {
683  /* lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) */
684  DISTPTS dl;
685  int rv;
686  LWLINE *lwline;
687  POINT2D P;
688 
689  /* Unit semi-circle above X axis */
690  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
691 
692  /* Point at origin */
693  P.x = P.y = 0;
695  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
696  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
697 
698  /* Point above arc on Y axis */
699  P.y = 2;
701  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
702  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
703 
704  /* Point 45 degrees off arc, 2 radii from center */
705  P.y = P.x = 2 * cos(M_PI/4);
707  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
708  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
709 
710  /* Four unit semi-circles surrounding the 2x2 box around origin */
711  lwline_free(lwline);
712  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)"));
713 
714  /* Point at origin */
715  P.x = P.y = 0;
717  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
718  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
719 
720  /* Point on box edge */
721  P.x = -1; P.y = 0;
723  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
724  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
725 
726  /* Point within a semicircle lobe */
727  P.x = -1.5; P.y = 0;
729  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
730  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
731 
732  /* Point outside a semicircle lobe */
733  P.x = -2.5; P.y = 0;
735  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
736  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
737 
738  /* Point outside a semicircle lobe */
739  P.y = -2.5; P.x = 0;
741  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
742  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
743 
744  /* Point outside a semicircle lobe */
745  P.y = 2; P.x = 1;
747  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
748  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1.0, 0.000001);
749 
750  /* Clean up */
751  lwline_free(lwline);
752 }
753 
754 static void
756 {
757  /* int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl) */
758  DISTPTS dl;
759  int rv;
760  LWLINE *lwline1;
761  LWLINE *lwline2;
762 
763  /* Unit semi-circle above X axis */
764  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
765 
766  /* Line above top of semi-circle */
767  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
769  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
770  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
771 
772  /* Reversed arguments, should fail */
775  rv = lw_dist2d_ptarray_ptarrayarc(lwline1->points, lwline2->points, &dl);
776  //printf("%s\n", cu_error_msg);
777  CU_ASSERT_STRING_EQUAL("lw_dist2d_ptarray_ptarrayarc called with non-arc input", cu_error_msg);
778 
779  lwline_free(lwline2);
780 
781  /* Line along side of semi-circle */
782  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
784  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
785  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
786 
787  /* Four unit semi-circles surrounding the 2x2 box around origin */
788  lwline_free(lwline1);
789  lwline_free(lwline2);
790  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)"));
791  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2.5 -3, -2.5 -2, -2.5 2, -2.5 3)"));
792  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
793  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
794 
795  lwline_free(lwline2);
796  lwline_free(lwline1);
797 }
798 
799 /*
800 ** Used by test harness to register the tests in this file.
801 */
802 void measures_suite_setup(void);
804 {
805  CU_pSuite suite = CU_add_suite("Measures", NULL, NULL);
817 }
static void test_lw_dist2d_seg_arc(void)
Definition: cu_measures.c:501
double lwgeom_mindistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:173
char * r
Definition: cu_in_wkt.c:25
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:655
Datum boundary(PG_FUNCTION_ARGS)
void lwfree(void *mem)
Definition: lwutil.c:190
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
#define DIST_MIN
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
void lwline_free(LWLINE *line)
Definition: lwline.c:63
static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line)
Definition: cu_measures.c:33
static LWGEOM * lwgeom_from_text(const char *str)
Definition: cu_measures.c:23
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:844
char ** result
Definition: liblwgeom.h:218
static void test_lw_dist2d_arc_arc(void)
Definition: cu_measures.c:573
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:1438
RECT_NODE * rect_tree_new(const POINTARRAY *pa)
Build a tree of nodes from a point array, one node per edge, and each with an associated measure rang...
Definition: lwtree.c:151
#define LW_FAILURE
Definition: liblwgeom.h:54
static void test_lw_dist2d_pt_arc(void)
Definition: cu_measures.c:442
static void test_lwgeom_segmentize2d(void)
Definition: cu_measures.c:403
static void test_lw_dist2d_pt_ptarrayarc(void)
Definition: cu_measures.c:681
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1706
double x
Definition: liblwgeom.h:284
int rect_tree_contains_point(const RECT_NODE *node, const POINT2D *pt, int *on_boundary)
Definition: lwtree.c:34
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
#define WKT_ISO
Definition: liblwgeom.h:1776
#define LW_FALSE
Definition: liblwgeom.h:52
int lwgeom_parse_wkt(LWGEOM_PARSER_RESULT *parser_result, char *wktstr, int parse_flags)
Parse a WKT geometry string into an LWGEOM structure.
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
void cu_error_msg_reset()
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM...
Definition: liblwgeom.h:1713
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:1143
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1395
#define PG_ADD_TEST(suite, testfunc)
POINTARRAY ** rings
Definition: liblwgeom.h:413
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:1039
double y
Definition: liblwgeom.h:284
Datum distance(PG_FUNCTION_ARGS)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:89
static void test_rect_tree_contains_point(void)
Definition: cu_measures.c:195
#define DIST2DTEST(str1, str2, res)
Definition: cu_measures.c:31
void measures_suite_setup(void)
Definition: cu_measures.c:803
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.
Note that p1 and p2 are pointers into an independent POINTARRAY, do not free them.
Definition: lwtree.h:4
double distance
Definition: measures.h:23
static void test_lwgeom_locate_along(void)
Definition: cu_measures.c:415
LWGEOM * lwgeom_segmentize2d(LWGEOM *line, double dist)
Definition: lwgeom.c:624
int rect_tree_intersects_tree(const RECT_NODE *n1, const RECT_NODE *n2)
Definition: lwtree.c:55
static void test_mindistance2d_tolerance(void)
Definition: cu_measures.c:71
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:1248
static void test_lw_dist2d_ptarray_ptarrayarc(void)
Definition: cu_measures.c:755
void rect_tree_free(RECT_NODE *node)
Recurse from top of node tree and free all children.
Definition: lwtree.c:18
Structure used in distance-calculations.
Definition: measures.h:21
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:27
char cu_error_msg[MAX_CUNIT_ERROR_LENGTH+1]
static void test_rect_tree_intersects_tree(void)
Definition: cu_measures.c:345
static void test_lw_arc_length(void)
Definition: cu_measures.c:639
POINTARRAY * points
Definition: liblwgeom.h:378