PostGIS  2.5.0dev-r@@SVN_REVISION@@
lwlinearreferencing.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2015 Sandro Santilli <strk@kbt.io>
22  * Copyright (C) 2011 Paul Ramsey
23  *
24  **********************************************************************/
25 
26 
27 #include "liblwgeom_internal.h"
28 #include "lwgeom_log.h"
29 #include "measures3d.h"
30 
31 static int
32 segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
33 {
34  double m1 = p1->m;
35  double m2 = p2->m;
36  double mprop;
37 
38  /* M is out of range, no new point generated. */
39  if ( (m < FP_MIN(m1,m2)) || (m > FP_MAX(m1,m2)) )
40  {
41  return LW_FALSE;
42  }
43 
44  if( m1 == m2 )
45  {
46  /* Degenerate case: same M on both points.
47  If they are the same point we just return one of them. */
48  if ( p4d_same(p1,p2) )
49  {
50  *pn = *p1;
51  return LW_TRUE;
52  }
53  /* If the points are different we split the difference */
54  mprop = 0.5;
55  }
56  else
57  {
58  mprop = (m - m1) / (m2 - m1);
59  }
60 
61  /* M is in range, new point to be generated. */
62  pn->x = p1->x + (p2->x - p1->x) * mprop;
63  pn->y = p1->y + (p2->y - p1->y) * mprop;
64  pn->z = p1->z + (p2->z - p1->z) * mprop;
65  pn->m = m;
66 
67  /* Offset to the left or right, if necessary. */
68  if ( offset != 0.0 )
69  {
70  double theta = atan2(p2->y - p1->y, p2->x - p1->x);
71  pn->x -= sin(theta) * offset;
72  pn->y += cos(theta) * offset;
73  }
74 
75  return LW_TRUE;
76 }
77 
78 
79 static POINTARRAY*
80 ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
81 {
82  uint32_t i;
83  POINT4D p1, p2, pn;
84  POINTARRAY *dpa = NULL;
85 
86  /* Can't do anything with degenerate point arrays */
87  if ( ! pa || pa->npoints < 2 ) return NULL;
88 
89  /* Walk through each segment in the point array */
90  for ( i = 1; i < pa->npoints; i++ )
91  {
92  getPoint4d_p(pa, i-1, &p1);
93  getPoint4d_p(pa, i, &p2);
94 
95  /* No derived point? Move to next segment. */
96  if ( segment_locate_along(&p1, &p2, m, offset, &pn) == LW_FALSE )
97  continue;
98 
99  /* No pointarray, make a fresh one */
100  if ( dpa == NULL )
102 
103  /* Add our new point to the array */
104  ptarray_append_point(dpa, &pn, 0);
105  }
106 
107  return dpa;
108 }
109 
110 static LWMPOINT*
111 lwline_locate_along(const LWLINE *lwline, double m, double offset)
112 {
113  POINTARRAY *opa = NULL;
114  LWMPOINT *mp = NULL;
115  LWGEOM *lwg = lwline_as_lwgeom(lwline);
116  int hasz, hasm, srid;
117 
118  /* Return degenerates upwards */
119  if ( ! lwline ) return NULL;
120 
121  /* Create empty return shell */
122  srid = lwgeom_get_srid(lwg);
123  hasz = lwgeom_has_z(lwg);
124  hasm = lwgeom_has_m(lwg);
125 
126  if ( hasm )
127  {
128  /* Find points along */
129  opa = ptarray_locate_along(lwline->points, m, offset);
130  }
131  else
132  {
133  LWLINE *lwline_measured = lwline_measured_from_lwline(lwline, 0.0, 1.0);
134  opa = ptarray_locate_along(lwline_measured->points, m, offset);
135  lwline_free(lwline_measured);
136  }
137 
138  /* Return NULL as EMPTY */
139  if ( ! opa )
140  return lwmpoint_construct_empty(srid, hasz, hasm);
141 
142  /* Convert pointarray into a multipoint */
143  mp = lwmpoint_construct(srid, opa);
144  ptarray_free(opa);
145  return mp;
146 }
147 
148 static LWMPOINT*
149 lwmline_locate_along(const LWMLINE *lwmline, double m, double offset)
150 {
151  LWMPOINT *lwmpoint = NULL;
152  LWGEOM *lwg = lwmline_as_lwgeom(lwmline);
153  uint32_t i, j;
154 
155  /* Return degenerates upwards */
156  if ( (!lwmline) || (lwmline->ngeoms < 1) ) return NULL;
157 
158  /* Construct return */
160 
161  /* Locate along each sub-line */
162  for ( i = 0; i < lwmline->ngeoms; i++ )
163  {
164  LWMPOINT *along = lwline_locate_along(lwmline->geoms[i], m, offset);
165  if ( along )
166  {
167  if ( ! lwgeom_is_empty((LWGEOM*)along) )
168  {
169  for ( j = 0; j < along->ngeoms; j++ )
170  {
171  lwmpoint_add_lwpoint(lwmpoint, along->geoms[j]);
172  }
173  }
174  /* Free the containing geometry, but leave the sub-geometries around */
175  along->ngeoms = 0;
176  lwmpoint_free(along);
177  }
178  }
179  return lwmpoint;
180 }
181 
182 static LWMPOINT*
183 lwpoint_locate_along(const LWPOINT *lwpoint, double m, double offset)
184 {
185  double point_m = lwpoint_get_m(lwpoint);
186  LWGEOM *lwg = lwpoint_as_lwgeom(lwpoint);
188  if ( FP_EQUALS(m, point_m) )
189  {
190  lwmpoint_add_lwpoint(r, lwpoint_clone(lwpoint));
191  }
192  return r;
193 }
194 
195 static LWMPOINT*
196 lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset)
197 {
198  LWGEOM *lwg = lwmpoint_as_lwgeom(lwin);
199  LWMPOINT *lwout = NULL;
200  uint32_t i;
201 
202  /* Construct return */
204 
205  for ( i = 0; i < lwin->ngeoms; i++ )
206  {
207  double point_m = lwpoint_get_m(lwin->geoms[i]);
208  if ( FP_EQUALS(m, point_m) )
209  {
210  lwmpoint_add_lwpoint(lwout, lwpoint_clone(lwin->geoms[i]));
211  }
212  }
213 
214  return lwout;
215 }
216 
217 LWGEOM*
218 lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
219 {
220  if ( ! lwin ) return NULL;
221 
222  if ( ! lwgeom_has_m(lwin) )
223  lwerror("Input geometry does not have a measure dimension");
224 
225  switch (lwin->type)
226  {
227  case POINTTYPE:
228  return (LWGEOM*)lwpoint_locate_along((LWPOINT*)lwin, m, offset);
229  case MULTIPOINTTYPE:
230  return (LWGEOM*)lwmpoint_locate_along((LWMPOINT*)lwin, m, offset);
231  case LINETYPE:
232  return (LWGEOM*)lwline_locate_along((LWLINE*)lwin, m, offset);
233  case MULTILINETYPE:
234  return (LWGEOM*)lwmline_locate_along((LWMLINE*)lwin, m, offset);
235  /* Only line types supported right now */
236  /* TO DO: CurveString, CompoundCurve, MultiCurve */
237  /* TO DO: Point, MultiPoint */
238  default:
239  lwerror("Only linear geometries are supported, %s provided.",lwtype_name(lwin->type));
240  return NULL;
241  }
242  return NULL;
243 }
244 
252 double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
253 {
254  if ( ! p )
255  {
256  lwerror("Null input geometry.");
257  return 0.0;
258  }
259 
260  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
261  {
262  lwerror("Cannot extract %c ordinate.", ordinate);
263  return 0.0;
264  }
265 
266  if ( ordinate == 'X' )
267  return p->x;
268  if ( ordinate == 'Y' )
269  return p->y;
270  if ( ordinate == 'Z' )
271  return p->z;
272  if ( ordinate == 'M' )
273  return p->m;
274 
275  /* X */
276  return p->x;
277 
278 }
279 
284 void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
285 {
286  if ( ! p )
287  {
288  lwerror("Null input geometry.");
289  return;
290  }
291 
292  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
293  {
294  lwerror("Cannot set %c ordinate.", ordinate);
295  return;
296  }
297 
298  LWDEBUGF(4, " setting ordinate %c to %g", ordinate, value);
299 
300  switch ( ordinate )
301  {
302  case 'X':
303  p->x = value;
304  return;
305  case 'Y':
306  p->y = value;
307  return;
308  case 'Z':
309  p->z = value;
310  return;
311  case 'M':
312  p->m = value;
313  return;
314  }
315 }
316 
322 int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
323 {
324  static char* dims = "XYZM";
325  double p1_value = lwpoint_get_ordinate(p1, ordinate);
326  double p2_value = lwpoint_get_ordinate(p2, ordinate);
327  double proportion;
328  int i = 0;
329 
330  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
331  {
332  lwerror("Cannot set %c ordinate.", ordinate);
333  return 0;
334  }
335 
336  if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
337  FP_MAX(p1_value, p2_value) < interpolation_value )
338  {
339  lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
340  return 0;
341  }
342 
343  proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
344 
345  for ( i = 0; i < 4; i++ )
346  {
347  double newordinate = 0.0;
348  if ( dims[i] == 'Z' && ! hasz ) continue;
349  if ( dims[i] == 'M' && ! hasm ) continue;
350  p1_value = lwpoint_get_ordinate(p1, dims[i]);
351  p2_value = lwpoint_get_ordinate(p2, dims[i]);
352  newordinate = p1_value + proportion * (p2_value - p1_value);
353  lwpoint_set_ordinate(p, dims[i], newordinate);
354  LWDEBUGF(4, " clip ordinate(%c) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", dims[i], p1_value, p2_value, proportion, newordinate );
355  }
356 
357  return 1;
358 }
359 
360 
365 lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
366 {
367  LWCOLLECTION *lwgeom_out = NULL;
368  char hasz, hasm;
369  POINT4D p4d;
370  double ordinate_value;
371 
372  /* Nothing to do with NULL */
373  if ( ! point )
374  lwerror("Null input geometry.");
375 
376  /* Ensure 'from' is less than 'to'. */
377  if ( to < from )
378  {
379  double t = from;
380  from = to;
381  to = t;
382  }
383 
384  /* Read Z/M info */
385  hasz = lwgeom_has_z(lwpoint_as_lwgeom(point));
386  hasm = lwgeom_has_m(lwpoint_as_lwgeom(point));
387 
388  /* Prepare return object */
389  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, point->srid, hasz, hasm);
390 
391  /* Test if ordinate is in range */
392  lwpoint_getPoint4d_p(point, &p4d);
393  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
394  if ( from <= ordinate_value && to >= ordinate_value )
395  {
396  LWPOINT *lwp = lwpoint_clone(point);
397  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
398  }
399 
400  /* Set the bbox, if necessary */
401  if ( lwgeom_out->bbox )
402  {
403  lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
404  }
405 
406  return lwgeom_out;
407 }
408 
409 
410 
415 lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
416 {
417  LWCOLLECTION *lwgeom_out = NULL;
418  char hasz, hasm;
419  uint32_t i;
420 
421  /* Nothing to do with NULL */
422  if ( ! mpoint )
423  lwerror("Null input geometry.");
424 
425  /* Ensure 'from' is less than 'to'. */
426  if ( to < from )
427  {
428  double t = from;
429  from = to;
430  to = t;
431  }
432 
433  /* Read Z/M info */
434  hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint));
435  hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint));
436 
437  /* Prepare return object */
438  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, mpoint->srid, hasz, hasm);
439 
440  /* For each point, is its ordinate value between from and to? */
441  for ( i = 0; i < mpoint->ngeoms; i ++ )
442  {
443  POINT4D p4d;
444  double ordinate_value;
445 
446  lwpoint_getPoint4d_p(mpoint->geoms[i], &p4d);
447  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
448 
449  if ( from <= ordinate_value && to >= ordinate_value )
450  {
451  LWPOINT *lwp = lwpoint_clone(mpoint->geoms[i]);
452  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
453  }
454  }
455 
456  /* Set the bbox, if necessary */
457  if ( lwgeom_out->bbox )
458  {
459  lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
460  }
461 
462  return lwgeom_out;
463 }
464 
469 lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
470 {
471  LWCOLLECTION *lwgeom_out = NULL;
472 
473  if ( ! mline )
474  {
475  lwerror("Null input geometry.");
476  return NULL;
477  }
478 
479  if ( mline->ngeoms == 1)
480  {
481  lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
482  }
483  else
484  {
485  LWCOLLECTION *col;
486  char hasz = lwgeom_has_z(lwmline_as_lwgeom(mline));
487  char hasm = lwgeom_has_m(lwmline_as_lwgeom(mline));
488  uint32_t i, j;
489  char homogeneous = 1;
490  size_t geoms_size = 0;
491  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
492  FLAGS_SET_Z(lwgeom_out->flags, hasz);
493  FLAGS_SET_M(lwgeom_out->flags, hasm);
494  for ( i = 0; i < mline->ngeoms; i ++ )
495  {
496  col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
497  if ( col )
498  {
499  /* Something was left after the clip. */
500  if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
501  {
502  geoms_size += 16;
503  if ( lwgeom_out->geoms )
504  {
505  lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
506  }
507  else
508  {
509  lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
510  }
511  }
512  for ( j = 0; j < col->ngeoms; j++ )
513  {
514  lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
515  lwgeom_out->ngeoms++;
516  }
517  if ( col->type != mline->type )
518  {
519  homogeneous = 0;
520  }
521  /* Shallow free the struct, leaving the geoms behind. */
522  if ( col->bbox ) lwfree(col->bbox);
523  lwfree(col->geoms);
524  lwfree(col);
525  }
526  }
527  if ( lwgeom_out->bbox )
528  {
529  lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
530  }
531 
532  if ( ! homogeneous )
533  {
534  lwgeom_out->type = COLLECTIONTYPE;
535  }
536  }
537 
538  return lwgeom_out;
539 
540 }
541 
542 
548 lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
549 {
550 
551  POINTARRAY *pa_in = NULL;
552  LWCOLLECTION *lwgeom_out = NULL;
553  POINTARRAY *dp = NULL;
554  uint32_t i;
555  int added_last_point = 0;
556  POINT4D *p = NULL, *q = NULL, *r = NULL;
557  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
558  char hasz, hasm;
559  char dims;
560 
561  /* Null input, nothing we can do. */
562  if ( ! line )
563  {
564  lwerror("Null input geometry.");
565  return NULL;
566  }
567  hasz = lwgeom_has_z(lwline_as_lwgeom(line));
568  hasm = lwgeom_has_m(lwline_as_lwgeom(line));
569  dims = FLAGS_NDIMS(line->flags);
570 
571  /* Ensure 'from' is less than 'to'. */
572  if ( to < from )
573  {
574  double t = from;
575  from = to;
576  to = t;
577  }
578 
579  LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
580  LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
581 
582  /* Asking for an ordinate we don't have. Error. */
583  if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
584  {
585  lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
586  return NULL;
587  }
588 
589  /* Prepare our working point objects. */
590  p = lwalloc(sizeof(POINT4D));
591  q = lwalloc(sizeof(POINT4D));
592  r = lwalloc(sizeof(POINT4D));
593 
594  /* Construct a collection to hold our outputs. */
595  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
596 
597  /* Get our input point array */
598  pa_in = line->points;
599 
600  for ( i = 0; i < pa_in->npoints; i++ )
601  {
602  LWDEBUGF(4, "Point #%d", i);
603  LWDEBUGF(4, "added_last_point %d", added_last_point);
604  if ( i > 0 )
605  {
606  *q = *p;
607  ordinate_value_q = ordinate_value_p;
608  }
609  getPoint4d_p(pa_in, i, p);
610  ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
611  LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
612  LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
613 
614  /* Is this point inside the ordinate range? Yes. */
615  if ( ordinate_value_p >= from && ordinate_value_p <= to )
616  {
617  LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
618 
619  if ( ! added_last_point )
620  {
621  LWDEBUG(4," new ptarray required");
622  /* We didn't add the previous point, so this is a new segment.
623  * Make a new point array. */
624  dp = ptarray_construct_empty(hasz, hasm, 32);
625 
626  /* We're transiting into the range so add an interpolated
627  * point at the range boundary.
628  * If we're on a boundary and crossing from the far side,
629  * we also need an interpolated point. */
630  if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
631  ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
632  ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
633  ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
634  {
635  double interpolation_value;
636  (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
637  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
639  LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
640  }
641  }
642  /* Add the current vertex to the point array. */
644  if ( ordinate_value_p == from || ordinate_value_p == to )
645  {
646  added_last_point = 2; /* Added on boundary. */
647  }
648  else
649  {
650  added_last_point = 1; /* Added inside range. */
651  }
652  }
653  /* Is this point inside the ordinate range? No. */
654  else
655  {
656  LWDEBUGF(4, " added_last_point (%d)", added_last_point);
657  if ( added_last_point == 1 )
658  {
659  /* We're transiting out of the range, so add an interpolated point
660  * to the point array at the range boundary. */
661  double interpolation_value;
662  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
663  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
665  LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
666  }
667  else if ( added_last_point == 2 )
668  {
669  /* We're out and the last point was on the boundary.
670  * If the last point was the near boundary, nothing to do.
671  * If it was the far boundary, we need an interpolated point. */
672  if ( from != to && (
673  (ordinate_value_q == from && ordinate_value_p > from) ||
674  (ordinate_value_q == to && ordinate_value_p < to) ) )
675  {
676  double interpolation_value;
677  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
678  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
680  LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
681  }
682  }
683  else if ( i && ordinate_value_q < from && ordinate_value_p > to )
684  {
685  /* We just hopped over the whole range, from bottom to top,
686  * so we need to add *two* interpolated points! */
687  dp = ptarray_construct(hasz, hasm, 2);
688  /* Interpolate lower point. */
689  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
690  ptarray_set_point4d(dp, 0, r);
691  /* Interpolate upper point. */
692  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
693  ptarray_set_point4d(dp, 1, r);
694  }
695  else if ( i && ordinate_value_q > to && ordinate_value_p < from )
696  {
697  /* We just hopped over the whole range, from top to bottom,
698  * so we need to add *two* interpolated points! */
699  dp = ptarray_construct(hasz, hasm, 2);
700  /* Interpolate upper point. */
701  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
702  ptarray_set_point4d(dp, 0, r);
703  /* Interpolate lower point. */
704  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
705  ptarray_set_point4d(dp, 1, r);
706  }
707  /* We have an extant point-array, save it out to a multi-line. */
708  if ( dp )
709  {
710  LWDEBUG(4, "saving pointarray to multi-line (1)");
711 
712  /* Only one point, so we have to make an lwpoint to hold this
713  * and set the overall output type to a generic collection. */
714  if ( dp->npoints == 1 )
715  {
716  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
717  lwgeom_out->type = COLLECTIONTYPE;
718  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
719 
720  }
721  else
722  {
723  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
724  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
725  }
726 
727  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
728  dp = NULL;
729  }
730  added_last_point = 0;
731 
732  }
733  }
734 
735  /* Still some points left to be saved out. */
736  if ( dp && dp->npoints > 0 )
737  {
738  LWDEBUG(4, "saving pointarray to multi-line (2)");
739  LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
740  LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
741 
742  if ( dp->npoints == 1 )
743  {
744  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
745  lwgeom_out->type = COLLECTIONTYPE;
746  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
747  }
748  else
749  {
750  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
751  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
752  }
753 
754  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
755  dp = NULL;
756  }
757 
758  lwfree(p);
759  lwfree(q);
760  lwfree(r);
761 
762  if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
763  {
764  lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
765  }
766 
767  return lwgeom_out;
768 
769 }
770 
772 lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
773 {
774  LWCOLLECTION *out_col;
775  LWCOLLECTION *out_offset;
776  uint32_t i;
777 
778  if ( ! lwin )
779  lwerror("lwgeom_clip_to_ordinate_range: null input geometry!");
780 
781  switch ( lwin->type )
782  {
783  case LINETYPE:
784  out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to);
785  break;
786  case MULTILINETYPE:
787  out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to);
788  break;
789  case MULTIPOINTTYPE:
790  out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to);
791  break;
792  case POINTTYPE:
793  out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to);
794  break;
795  default:
796  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
797  return NULL;;
798  }
799 
800  /* Stop if result is NULL */
801  if ( out_col == NULL )
802  lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL");
803 
804  /* Return if we aren't going to offset the result */
805  if (FP_IS_ZERO(offset) ||
807  return out_col;
808 
809  /* Construct a collection to hold our outputs. */
810  /* Things get ugly: GEOS offset drops Z's and M's so we have to drop ours */
811  out_offset = lwcollection_construct_empty(MULTILINETYPE, lwin->srid, 0, 0);
812 
813  /* Try and offset the linear portions of the return value */
814  for ( i = 0; i < out_col->ngeoms; i++ )
815  {
816  int type = out_col->geoms[i]->type;
817  if ( type == POINTTYPE )
818  {
819  lwnotice("lwgeom_clip_to_ordinate_range cannot offset a clipped point");
820  continue;
821  }
822  else if ( type == LINETYPE )
823  {
824  /* lwgeom_offsetcurve(line, offset, quadsegs, joinstyle (round), mitrelimit) */
825  LWGEOM *lwoff = lwgeom_offsetcurve(lwgeom_as_lwline(out_col->geoms[i]), offset, 8, 1, 5.0);
826  if ( ! lwoff )
827  {
828  lwerror("lwgeom_offsetcurve returned null");
829  }
830  lwcollection_add_lwgeom(out_offset, lwoff);
831  }
832  else
833  {
834  lwerror("lwgeom_clip_to_ordinate_range found an unexpected type (%s) in the offset routine",lwtype_name(type));
835  }
836  }
837 
838  return out_offset;
839 }
840 
842 lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset)
843 {
844  if ( ! lwgeom_has_m(lwin) )
845  lwerror("Input geometry does not have a measure dimension");
846 
847  return lwgeom_clip_to_ordinate_range(lwin, 'M', from, to, offset);
848 }
849 
850 double
851 lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
852 {
853  POINT4D p, p_proj;
854  double ret = 0.0;
855 
856  if ( ! lwin )
857  lwerror("lwgeom_interpolate_point: null input geometry!");
858 
859  if ( ! lwgeom_has_m(lwin) )
860  lwerror("Input geometry does not have a measure dimension");
861 
862  if ( lwgeom_is_empty(lwin) || lwpoint_is_empty(lwpt) )
863  lwerror("Input geometry is empty");
864 
865  switch ( lwin->type )
866  {
867  case LINETYPE:
868  {
869  LWLINE *lwline = lwgeom_as_lwline(lwin);
870  lwpoint_getPoint4d_p(lwpt, &p);
871  ret = ptarray_locate_point(lwline->points, &p, NULL, &p_proj);
872  ret = p_proj.m;
873  break;
874  }
875  default:
876  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
877  }
878  return ret;
879 }
880 
881 /*
882  * Time of closest point of approach
883  *
884  * Given two vectors (p1-p2 and q1-q2) and
885  * a time range (t1-t2) return the time in which
886  * a point p is closest to a point q on their
887  * respective vectors, and the actual points
888  *
889  * Here we use algorithm from softsurfer.com
890  * that can be found here
891  * http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
892  *
893  * @param p0 start of first segment, will be set to actual
894  * closest point of approach on segment.
895  * @param p1 end of first segment
896  * @param q0 start of second segment, will be set to actual
897  * closest point of approach on segment.
898  * @param q1 end of second segment
899  * @param t0 start of travel time
900  * @param t1 end of travel time
901  *
902  * @return time of closest point of approach
903  *
904  */
905 static double
907  POINT4D* q0, const POINT4D* q1,
908  double t0, double t1)
909 {
910  POINT3DZ pv; /* velocity of p, aka u */
911  POINT3DZ qv; /* velocity of q, aka v */
912  POINT3DZ dv; /* velocity difference */
913  POINT3DZ w0; /* vector between first points */
914 
915  /*
916  lwnotice("FROM %g,%g,%g,%g -- %g,%g,%g,%g",
917  p0->x, p0->y, p0->z, p0->m,
918  p1->x, p1->y, p1->z, p1->m);
919  lwnotice(" TO %g,%g,%g,%g -- %g,%g,%g,%g",
920  q0->x, q0->y, q0->z, q0->m,
921  q1->x, q1->y, q1->z, q1->m);
922  */
923 
924  /* PV aka U */
925  pv.x = ( p1->x - p0->x );
926  pv.y = ( p1->y - p0->y );
927  pv.z = ( p1->z - p0->z );
928  /*lwnotice("PV: %g, %g, %g", pv.x, pv.y, pv.z);*/
929 
930  /* QV aka V */
931  qv.x = ( q1->x - q0->x );
932  qv.y = ( q1->y - q0->y );
933  qv.z = ( q1->z - q0->z );
934  /*lwnotice("QV: %g, %g, %g", qv.x, qv.y, qv.z);*/
935 
936  dv.x = pv.x - qv.x;
937  dv.y = pv.y - qv.y;
938  dv.z = pv.z - qv.z;
939  /*lwnotice("DV: %g, %g, %g", dv.x, dv.y, dv.z);*/
940 
941  double dv2 = DOT(dv,dv);
942  /*lwnotice("DOT: %g", dv2);*/
943 
944  if ( dv2 == 0.0 )
945  {
946  /* Distance is the same at any time, we pick the earliest */
947  return t0;
948  }
949 
950  /* Distance at any given time, with t0 */
951  w0.x = ( p0->x - q0->x );
952  w0.y = ( p0->y - q0->y );
953  w0.z = ( p0->z - q0->z );
954 
955  /*lwnotice("W0: %g, %g, %g", w0.x, w0.y, w0.z);*/
956 
957  /* Check that at distance dt w0 is distance */
958 
959  /* This is the fraction of measure difference */
960  double t = -DOT(w0,dv) / dv2;
961  /*lwnotice("CLOSEST TIME (fraction): %g", t);*/
962 
963  if ( t > 1.0 )
964  {
965  /* Getting closer as we move to the end */
966  /*lwnotice("Converging");*/
967  t = 1;
968  }
969  else if ( t < 0.0 )
970  {
971  /*lwnotice("Diverging");*/
972  t = 0;
973  }
974 
975  /* Interpolate the actual points now */
976 
977  p0->x += pv.x * t;
978  p0->y += pv.y * t;
979  p0->z += pv.z * t;
980 
981  q0->x += qv.x * t;
982  q0->y += qv.y * t;
983  q0->z += qv.z * t;
984 
985  t = t0 + (t1 - t0) * t;
986  /*lwnotice("CLOSEST TIME (real): %g", t);*/
987 
988  return t;
989 }
990 
991 static int
992 ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals)
993 {
994  POINT4D pbuf;
995  uint32_t i, n=0;
996  for (i=0; i<pa->npoints; ++i)
997  {
998  getPoint4d_p(pa, i, &pbuf); /* could be optimized */
999  if ( pbuf.m >= tmin && pbuf.m <= tmax )
1000  mvals[n++] = pbuf.m;
1001  }
1002  return n;
1003 }
1004 
1005 static int
1006 compare_double(const void *pa, const void *pb)
1007 {
1008  double a = *((double *)pa);
1009  double b = *((double *)pb);
1010  if ( a < b )
1011  return -1;
1012  else if ( a > b )
1013  return 1;
1014  else
1015  return 0;
1016 }
1017 
1018 /* Return number of elements in unique array */
1019 static int
1020 uniq(double *vals, int nvals)
1021 {
1022  int i, last=0;
1023  for (i=1; i<nvals; ++i)
1024  {
1025  // lwnotice("(I%d):%g", i, vals[i]);
1026  if ( vals[i] != vals[last] )
1027  {
1028  vals[++last] = vals[i];
1029  // lwnotice("(O%d):%g", last, vals[last]);
1030  }
1031  }
1032  return last+1;
1033 }
1034 
1035 /*
1036  * Find point at a given measure
1037  *
1038  * The function assumes measures are linear so that always a single point
1039  * is returned for a single measure.
1040  *
1041  * @param pa the point array to perform search on
1042  * @param m the measure to search for
1043  * @param p the point to write result into
1044  * @param from the segment number to start from
1045  *
1046  * @return the segment number the point was found into
1047  * or -1 if given measure was out of the known range.
1048  */
1049 static int
1051 {
1052  uint32_t i = from;
1053  POINT4D p1, p2;
1054 
1055  /* Walk through each segment in the point array */
1056  getPoint4d_p(pa, i, &p1);
1057  for ( i = from+1; i < pa->npoints; i++ )
1058  {
1059  getPoint4d_p(pa, i, &p2);
1060 
1061  if ( segment_locate_along(&p1, &p2, m, 0, p) == LW_TRUE )
1062  return i-1; /* found */
1063 
1064  p1 = p2;
1065  }
1066 
1067  return -1; /* not found */
1068 }
1069 
1070 double
1071 lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
1072 {
1073  LWLINE *l1, *l2;
1074  int i;
1075  const GBOX *gbox1, *gbox2;
1076  double tmin, tmax;
1077  double *mvals;
1078  int nmvals = 0;
1079  double mintime;
1080  double mindist2 = FLT_MAX; /* minimum distance, squared */
1081 
1082  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1083  {
1084  lwerror("Both input geometries must have a measure dimension");
1085  return -1;
1086  }
1087 
1088  l1 = lwgeom_as_lwline(g1);
1089  l2 = lwgeom_as_lwline(g2);
1090 
1091  if ( ! l1 || ! l2 )
1092  {
1093  lwerror("Both input geometries must be linestrings");
1094  return -1;
1095  }
1096 
1097  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1098  {
1099  lwerror("Both input lines must have at least 2 points");
1100  return -1;
1101  }
1102 
1103  /* WARNING: these ranges may be wider than real ones */
1104  gbox1 = lwgeom_get_bbox(g1);
1105  gbox2 = lwgeom_get_bbox(g2);
1106 
1107  assert(gbox1); /* or the npoints check above would have failed */
1108  assert(gbox2); /* or the npoints check above would have failed */
1109 
1110  /*
1111  * Find overlapping M range
1112  * WARNING: may be larger than the real one
1113  */
1114 
1115  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1116  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1117 
1118  if ( tmax < tmin )
1119  {
1120  LWDEBUG(1, "Inputs never exist at the same time");
1121  return -2;
1122  }
1123 
1124  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1125 
1126  /*
1127  * Collect M values in common time range from inputs
1128  */
1129 
1130  mvals = lwalloc( sizeof(double) *
1131  ( l1->points->npoints + l2->points->npoints ) );
1132 
1133  /* TODO: also clip the lines ? */
1134  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1135  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1136 
1137  /* Sort values in ascending order */
1138  qsort(mvals, nmvals, sizeof(double), compare_double);
1139 
1140  /* Remove duplicated values */
1141  nmvals = uniq(mvals, nmvals);
1142 
1143  if ( nmvals < 2 )
1144  {
1145  {
1146  /* there's a single time, must be that one... */
1147  double t0 = mvals[0];
1148  POINT4D p0, p1;
1149  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1150  if ( mindist )
1151  {
1152  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1153  {
1154  lwfree(mvals);
1155  lwerror("Could not find point with M=%g on first geom", t0);
1156  return -1;
1157  }
1158  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1159  {
1160  lwfree(mvals);
1161  lwerror("Could not find point with M=%g on second geom", t0);
1162  return -1;
1163  }
1164  *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1);
1165  }
1166  lwfree(mvals);
1167  return t0;
1168  }
1169  }
1170 
1171  /*
1172  * For each consecutive pair of measures, compute time of closest point
1173  * approach and actual distance between points at that time
1174  */
1175  mintime = tmin;
1176  for (i=1; i<nmvals; ++i)
1177  {
1178  double t0 = mvals[i-1];
1179  double t1 = mvals[i];
1180  double t;
1181  POINT4D p0, p1, q0, q1;
1182  int seg;
1183  double dist2;
1184 
1185  // lwnotice("T %g-%g", t0, t1);
1186 
1187  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1188  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1189  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1190 
1191  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1192  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1193  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1194 
1195  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1196  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1197  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1198 
1199  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1200  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1201  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1202 
1203  t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1204 
1205  /*
1206  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1207  p0.x, p0.y, p0.z,
1208  q0.x, q0.y, q0.z, t);
1209  */
1210 
1211  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1212  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1213  ( q0.z - p0.z ) * ( q0.z - p0.z );
1214  if ( dist2 < mindist2 )
1215  {
1216  mindist2 = dist2;
1217  mintime = t;
1218  // lwnotice("MINTIME: %g", mintime);
1219  }
1220  }
1221 
1222  /*
1223  * Release memory
1224  */
1225 
1226  lwfree(mvals);
1227 
1228  if ( mindist )
1229  {
1230  *mindist = sqrt(mindist2);
1231  }
1232  /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
1233 
1234  return mintime;
1235 }
1236 
1237 int
1238 lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
1239 {
1240  LWLINE *l1, *l2;
1241  int i;
1242  const GBOX *gbox1, *gbox2;
1243  double tmin, tmax;
1244  double *mvals;
1245  int nmvals = 0;
1246  double maxdist2 = maxdist * maxdist;
1247  int within = LW_FALSE;
1248 
1249  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1250  {
1251  lwerror("Both input geometries must have a measure dimension");
1252  return LW_FALSE;
1253  }
1254 
1255  l1 = lwgeom_as_lwline(g1);
1256  l2 = lwgeom_as_lwline(g2);
1257 
1258  if ( ! l1 || ! l2 )
1259  {
1260  lwerror("Both input geometries must be linestrings");
1261  return LW_FALSE;
1262  }
1263 
1264  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1265  {
1266  /* TODO: return distance between these two points */
1267  lwerror("Both input lines must have at least 2 points");
1268  return LW_FALSE;
1269  }
1270 
1271  /* WARNING: these ranges may be wider than real ones */
1272  gbox1 = lwgeom_get_bbox(g1);
1273  gbox2 = lwgeom_get_bbox(g2);
1274 
1275  assert(gbox1); /* or the npoints check above would have failed */
1276  assert(gbox2); /* or the npoints check above would have failed */
1277 
1278  /*
1279  * Find overlapping M range
1280  * WARNING: may be larger than the real one
1281  */
1282 
1283  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1284  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1285 
1286  if ( tmax < tmin )
1287  {
1288  LWDEBUG(1, "Inputs never exist at the same time");
1289  return LW_FALSE;
1290  }
1291 
1292  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1293 
1294  /*
1295  * Collect M values in common time range from inputs
1296  */
1297 
1298  mvals = lwalloc( sizeof(double) *
1299  ( l1->points->npoints + l2->points->npoints ) );
1300 
1301  /* TODO: also clip the lines ? */
1302  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1303  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1304 
1305  /* Sort values in ascending order */
1306  qsort(mvals, nmvals, sizeof(double), compare_double);
1307 
1308  /* Remove duplicated values */
1309  nmvals = uniq(mvals, nmvals);
1310 
1311  if ( nmvals < 2 )
1312  {
1313  /* there's a single time, must be that one... */
1314  double t0 = mvals[0];
1315  POINT4D p0, p1;
1316  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1317  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1318  {
1319  lwnotice("Could not find point with M=%g on first geom", t0);
1320  return LW_FALSE;
1321  }
1322  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1323  {
1324  lwnotice("Could not find point with M=%g on second geom", t0);
1325  return LW_FALSE;
1326  }
1327  if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist )
1328  within = LW_TRUE;
1329  lwfree(mvals);
1330  return within;
1331  }
1332 
1333  /*
1334  * For each consecutive pair of measures, compute time of closest point
1335  * approach and actual distance between points at that time
1336  */
1337  for (i=1; i<nmvals; ++i)
1338  {
1339  double t0 = mvals[i-1];
1340  double t1 = mvals[i];
1341 #if POSTGIS_DEBUG_LEVEL >= 1
1342  double t;
1343 #endif
1344  POINT4D p0, p1, q0, q1;
1345  int seg;
1346  double dist2;
1347 
1348  // lwnotice("T %g-%g", t0, t1);
1349 
1350  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1351  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1352  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1353 
1354  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1355  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1356  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1357 
1358  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1359  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1360  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1361 
1362  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1363  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1364  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1365 
1366 #if POSTGIS_DEBUG_LEVEL >= 1
1367  t =
1368 #endif
1369  segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1370 
1371  /*
1372  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1373  p0.x, p0.y, p0.z,
1374  q0.x, q0.y, q0.z, t);
1375  */
1376 
1377  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1378  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1379  ( q0.z - p0.z ) * ( q0.z - p0.z );
1380  if ( dist2 <= maxdist2 )
1381  {
1382  LWDEBUGF(1, "Within distance %g at time %g, breaking", sqrt(dist2), t);
1383  within = LW_TRUE;
1384  break;
1385  }
1386  }
1387 
1388  /*
1389  * Release memory
1390  */
1391 
1392  lwfree(mvals);
1393 
1394  return within;
1395 }
double x
Definition: liblwgeom.h:351
#define LINETYPE
Definition: liblwgeom.h:85
double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
Find the measure value at the location on the line closest to the point.
LWLINE * lwline_measured_from_lwline(const LWLINE *lwline, double m_start, double m_end)
Add a measure dimension to a line, interpolating linearly from the start to the end value...
Definition: lwline.c:388
void lwgeom_refresh_bbox(LWGEOM *lwgeom)
Drop current bbox and calculate a fresh one.
Definition: lwgeom.c:691
uint8_t type
Definition: liblwgeom.h:476
double z
Definition: liblwgeom.h:333
double y
Definition: liblwgeom.h:333
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
Given a POINT4D and an ordinate number, return the value of the ordinate.
LWCOLLECTION * lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
Clip an input POINT between two values, on any ordinate input.
double m
Definition: liblwgeom.h:351
uint32_t ngeoms
Definition: liblwgeom.h:467
LWCOLLECTION * lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
Clip an input MULTIPOINT between two values, on any ordinate input.
double x
Definition: liblwgeom.h:333
char * r
Definition: cu_in_wkt.c:24
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
static LWMPOINT * lwline_locate_along(const LWLINE *lwline, double m, double offset)
uint8_t type
Definition: liblwgeom.h:502
static int segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
static int uniq(double *vals, int nvals)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:549
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:328
static int compare_double(const void *pa, const void *pb)
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
LWMPOINT * lwmpoint_construct(int srid, const POINTARRAY *pa)
Definition: lwmpoint.c:52
LWGEOM * lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
void lwline_free(LWLINE *line)
Definition: lwline.c:76
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:425
GBOX * bbox
Definition: liblwgeom.h:504
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:239
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:909
LWGEOM * lwmline_as_lwgeom(const LWMLINE *obj)
Definition: lwgeom.c:289
int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
Given two points, a dimensionality, an ordinate, and an interpolation value generate a new point that...
int32_t srid
Definition: liblwgeom.h:420
#define FP_MIN(A, B)
uint32_t ngeoms
Definition: liblwgeom.h:506
LWCOLLECTION * lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
Take in a LINESTRING and return a MULTILINESTRING of those portions of the LINESTRING between the fro...
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:923
static POINTARRAY * ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
int32_t srid
Definition: liblwgeom.h:398
#define FP_IS_ZERO(A)
LWCOLLECTION * lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset)
Determine the segments along a measured line that fall within the m-range given.
#define DOT(u, v)
Definition: measures3d.h:31
unsigned int uint32_t
Definition: uthash.h:78
uint8_t flags
Definition: liblwgeom.h:503
double ptarray_locate_point(const POINTARRAY *pa, const POINT4D *pt, double *dist, POINT4D *p_located)
Definition: ptarray.c:1298
#define FLAGS_SET_Z(flags, value)
Definition: liblwgeom.h:145
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:72
static int ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals)
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:329
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_FALSE, then a duplicate point will not be added.
Definition: ptarray.c:156
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
LWGEOM ** geoms
Definition: liblwgeom.h:508
void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
Given a point, ordinate number and value, set that ordinate on the point.
static LWMPOINT * lwmline_locate_along(const LWMLINE *lwmline, double m, double offset)
static int ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, uint32_t from)
const GBOX * lwgeom_get_bbox(const LWGEOM *lwgeom)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:727
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition: lwpoint.c:57
LWPOINT ** geoms
Definition: liblwgeom.h:469
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:39
LWCOLLECTION * lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
Clip an input MULTILINESTRING between two values, on any ordinate input.
double z
Definition: liblwgeom.h:351
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:169
int32_t srid
Definition: liblwgeom.h:409
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
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.
uint32_t ngeoms
Definition: liblwgeom.h:480
LWLINE ** geoms
Definition: liblwgeom.h:482
double lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
Find the time of closest point of approach.
double lwpoint_get_m(const LWPOINT *point)
Definition: lwpoint.c:107
int lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
Is the closest point of approach within a distance ?
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
static LWMPOINT * lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset)
double mmin
Definition: liblwgeom.h:297
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:792
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:237
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int lwpoint_is_empty(const LWPOINT *point)
Definition: lwpoint.c:291
LWCOLLECTION * lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
Given a geometry clip based on the from/to range of one of its ordinates (x, y, z, m).
uint8_t type
Definition: liblwgeom.h:395
#define FP_EQUALS(A, B)
double mmax
Definition: liblwgeom.h:298
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:129
void * lwalloc(size_t size)
Definition: lwutil.c:229
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1386
double y
Definition: liblwgeom.h:351
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
uint8_t flags
Definition: liblwgeom.h:418
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
int32_t srid
Definition: liblwgeom.h:466
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:930
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
static double segments_tcpa(POINT4D *p0, const POINT4D *p1, POINT4D *q0, const POINT4D *q1, double t0, double t1)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
int p4d_same(const POINT4D *p1, const POINT4D *p2)
Definition: lwalgorithm.c:31
static LWMPOINT * lwpoint_locate_along(const LWPOINT *lwpoint, double m, double offset)
#define FP_MAX(A, B)
int32_t srid
Definition: liblwgeom.h:479
#define FLAGS_SET_M(flags, value)
Definition: liblwgeom.h:146
POINTARRAY * points
Definition: liblwgeom.h:421
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:299
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:294
uint32_t npoints
Definition: liblwgeom.h:370