PostGIS  2.3.7dev-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  int 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  int 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  int 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_drop_bbox((LWGEOM*)lwgeom_out);
404  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
405  }
406 
407  return lwgeom_out;
408 }
409 
410 
411 
416 lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
417 {
418  LWCOLLECTION *lwgeom_out = NULL;
419  char hasz, hasm;
420  int i;
421 
422  /* Nothing to do with NULL */
423  if ( ! mpoint )
424  lwerror("Null input geometry.");
425 
426  /* Ensure 'from' is less than 'to'. */
427  if ( to < from )
428  {
429  double t = from;
430  from = to;
431  to = t;
432  }
433 
434  /* Read Z/M info */
435  hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint));
436  hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint));
437 
438  /* Prepare return object */
439  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, mpoint->srid, hasz, hasm);
440 
441  /* For each point, is its ordinate value between from and to? */
442  for ( i = 0; i < mpoint->ngeoms; i ++ )
443  {
444  POINT4D p4d;
445  double ordinate_value;
446 
447  lwpoint_getPoint4d_p(mpoint->geoms[i], &p4d);
448  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
449 
450  if ( from <= ordinate_value && to >= ordinate_value )
451  {
452  LWPOINT *lwp = lwpoint_clone(mpoint->geoms[i]);
453  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
454  }
455  }
456 
457  /* Set the bbox, if necessary */
458  if ( lwgeom_out->bbox )
459  {
460  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
461  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
462  }
463 
464  return lwgeom_out;
465 }
466 
471 lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
472 {
473  LWCOLLECTION *lwgeom_out = NULL;
474 
475  if ( ! mline )
476  {
477  lwerror("Null input geometry.");
478  return NULL;
479  }
480 
481  if ( mline->ngeoms == 1)
482  {
483  lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
484  }
485  else
486  {
487  LWCOLLECTION *col;
488  char hasz = lwgeom_has_z(lwmline_as_lwgeom(mline));
489  char hasm = lwgeom_has_m(lwmline_as_lwgeom(mline));
490  int i, j;
491  char homogeneous = 1;
492  size_t geoms_size = 0;
493  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
494  FLAGS_SET_Z(lwgeom_out->flags, hasz);
495  FLAGS_SET_M(lwgeom_out->flags, hasm);
496  for ( i = 0; i < mline->ngeoms; i ++ )
497  {
498  col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
499  if ( col )
500  {
501  /* Something was left after the clip. */
502  if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
503  {
504  geoms_size += 16;
505  if ( lwgeom_out->geoms )
506  {
507  lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
508  }
509  else
510  {
511  lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
512  }
513  }
514  for ( j = 0; j < col->ngeoms; j++ )
515  {
516  lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
517  lwgeom_out->ngeoms++;
518  }
519  if ( col->type != mline->type )
520  {
521  homogeneous = 0;
522  }
523  /* Shallow free the struct, leaving the geoms behind. */
524  if ( col->bbox ) lwfree(col->bbox);
525  lwfree(col->geoms);
526  lwfree(col);
527  }
528  }
529  if ( lwgeom_out->bbox )
530  {
531  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
532  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
533  }
534 
535  if ( ! homogeneous )
536  {
537  lwgeom_out->type = COLLECTIONTYPE;
538  }
539  }
540 
541  return lwgeom_out;
542 
543 }
544 
545 
551 lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
552 {
553 
554  POINTARRAY *pa_in = NULL;
555  LWCOLLECTION *lwgeom_out = NULL;
556  POINTARRAY *dp = NULL;
557  int i;
558  int added_last_point = 0;
559  POINT4D *p = NULL, *q = NULL, *r = NULL;
560  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
561  char hasz = lwgeom_has_z(lwline_as_lwgeom(line));
562  char hasm = lwgeom_has_m(lwline_as_lwgeom(line));
563  char dims = FLAGS_NDIMS(line->flags);
564 
565  /* Null input, nothing we can do. */
566  if ( ! line )
567  {
568  lwerror("Null input geometry.");
569  return NULL;
570  }
571 
572  /* Ensure 'from' is less than 'to'. */
573  if ( to < from )
574  {
575  double t = from;
576  from = to;
577  to = t;
578  }
579 
580  LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
581  LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
582 
583  /* Asking for an ordinate we don't have. Error. */
584  if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
585  {
586  lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
587  return NULL;
588  }
589 
590  /* Prepare our working point objects. */
591  p = lwalloc(sizeof(POINT4D));
592  q = lwalloc(sizeof(POINT4D));
593  r = lwalloc(sizeof(POINT4D));
594 
595  /* Construct a collection to hold our outputs. */
596  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
597 
598  /* Get our input point array */
599  pa_in = line->points;
600 
601  for ( i = 0; i < pa_in->npoints; i++ )
602  {
603  LWDEBUGF(4, "Point #%d", i);
604  LWDEBUGF(4, "added_last_point %d", added_last_point);
605  if ( i > 0 )
606  {
607  *q = *p;
608  ordinate_value_q = ordinate_value_p;
609  }
610  getPoint4d_p(pa_in, i, p);
611  ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
612  LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
613  LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
614 
615  /* Is this point inside the ordinate range? Yes. */
616  if ( ordinate_value_p >= from && ordinate_value_p <= to )
617  {
618  LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
619 
620  if ( ! added_last_point )
621  {
622  LWDEBUG(4," new ptarray required");
623  /* We didn't add the previous point, so this is a new segment.
624  * Make a new point array. */
625  dp = ptarray_construct_empty(hasz, hasm, 32);
626 
627  /* We're transiting into the range so add an interpolated
628  * point at the range boundary.
629  * If we're on a boundary and crossing from the far side,
630  * we also need an interpolated point. */
631  if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
632  ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
633  ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
634  ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
635  {
636  double interpolation_value;
637  (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
638  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
640  LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
641  }
642  }
643  /* Add the current vertex to the point array. */
645  if ( ordinate_value_p == from || ordinate_value_p == to )
646  {
647  added_last_point = 2; /* Added on boundary. */
648  }
649  else
650  {
651  added_last_point = 1; /* Added inside range. */
652  }
653  }
654  /* Is this point inside the ordinate range? No. */
655  else
656  {
657  LWDEBUGF(4, " added_last_point (%d)", added_last_point);
658  if ( added_last_point == 1 )
659  {
660  /* We're transiting out of the range, so add an interpolated point
661  * to the point array at the range boundary. */
662  double interpolation_value;
663  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
664  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
666  LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
667  }
668  else if ( added_last_point == 2 )
669  {
670  /* We're out and the last point was on the boundary.
671  * If the last point was the near boundary, nothing to do.
672  * If it was the far boundary, we need an interpolated point. */
673  if ( from != to && (
674  (ordinate_value_q == from && ordinate_value_p > from) ||
675  (ordinate_value_q == to && ordinate_value_p < to) ) )
676  {
677  double interpolation_value;
678  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
679  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
681  LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
682  }
683  }
684  else if ( i && ordinate_value_q < from && ordinate_value_p > to )
685  {
686  /* We just hopped over the whole range, from bottom to top,
687  * so we need to add *two* interpolated points! */
688  dp = ptarray_construct(hasz, hasm, 2);
689  /* Interpolate lower point. */
690  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
691  ptarray_set_point4d(dp, 0, r);
692  /* Interpolate upper point. */
693  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
694  ptarray_set_point4d(dp, 1, r);
695  }
696  else if ( i && ordinate_value_q > to && ordinate_value_p < from )
697  {
698  /* We just hopped over the whole range, from top to bottom,
699  * so we need to add *two* interpolated points! */
700  dp = ptarray_construct(hasz, hasm, 2);
701  /* Interpolate upper point. */
702  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
703  ptarray_set_point4d(dp, 0, r);
704  /* Interpolate lower point. */
705  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
706  ptarray_set_point4d(dp, 1, r);
707  }
708  /* We have an extant point-array, save it out to a multi-line. */
709  if ( dp )
710  {
711  LWDEBUG(4, "saving pointarray to multi-line (1)");
712 
713  /* Only one point, so we have to make an lwpoint to hold this
714  * and set the overall output type to a generic collection. */
715  if ( dp->npoints == 1 )
716  {
717  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
718  lwgeom_out->type = COLLECTIONTYPE;
719  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
720 
721  }
722  else
723  {
724  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
725  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
726  }
727 
728  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
729  dp = NULL;
730  }
731  added_last_point = 0;
732 
733  }
734  }
735 
736  /* Still some points left to be saved out. */
737  if ( dp && dp->npoints > 0 )
738  {
739  LWDEBUG(4, "saving pointarray to multi-line (2)");
740  LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
741  LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
742 
743  if ( dp->npoints == 1 )
744  {
745  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
746  lwgeom_out->type = COLLECTIONTYPE;
747  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
748  }
749  else
750  {
751  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
752  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
753  }
754 
755  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
756  dp = NULL;
757  }
758 
759  lwfree(p);
760  lwfree(q);
761  lwfree(r);
762 
763  if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
764  {
765  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
766  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
767  }
768 
769  return lwgeom_out;
770 
771 }
772 
774 lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
775 {
776  LWCOLLECTION *out_col;
777  LWCOLLECTION *out_offset;
778  int i;
779 
780  if ( ! lwin )
781  lwerror("lwgeom_clip_to_ordinate_range: null input geometry!");
782 
783  switch ( lwin->type )
784  {
785  case LINETYPE:
786  out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to);
787  break;
788  case MULTILINETYPE:
789  out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to);
790  break;
791  case MULTIPOINTTYPE:
792  out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to);
793  break;
794  case POINTTYPE:
795  out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to);
796  break;
797  default:
798  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
799  return NULL;;
800  }
801 
802  /* Stop if result is NULL */
803  if ( out_col == NULL )
804  lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL");
805 
806  /* Return if we aren't going to offset the result */
807  if ( FP_EQUALS(offset, 0.0) || lwgeom_is_empty(lwcollection_as_lwgeom(out_col)) )
808  return out_col;
809 
810  /* Construct a collection to hold our outputs. */
811  /* Things get ugly: GEOS offset drops Z's and M's so we have to drop ours */
812  out_offset = lwcollection_construct_empty(MULTILINETYPE, lwin->srid, 0, 0);
813 
814  /* Try and offset the linear portions of the return value */
815  for ( i = 0; i < out_col->ngeoms; i++ )
816  {
817  int type = out_col->geoms[i]->type;
818  if ( type == POINTTYPE )
819  {
820  lwnotice("lwgeom_clip_to_ordinate_range cannot offset a clipped point");
821  continue;
822  }
823  else if ( type == LINETYPE )
824  {
825  /* lwgeom_offsetcurve(line, offset, quadsegs, joinstyle (round), mitrelimit) */
826  LWGEOM *lwoff = lwgeom_offsetcurve(lwgeom_as_lwline(out_col->geoms[i]), offset, 8, 1, 5.0);
827  if ( ! lwoff )
828  {
829  lwerror("lwgeom_offsetcurve returned null");
830  }
831  lwcollection_add_lwgeom(out_offset, lwoff);
832  }
833  else
834  {
835  lwerror("lwgeom_clip_to_ordinate_range found an unexpected type (%s) in the offset routine",lwtype_name(type));
836  }
837  }
838 
839  return out_offset;
840 }
841 
843 lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset)
844 {
845  if ( ! lwgeom_has_m(lwin) )
846  lwerror("Input geometry does not have a measure dimension");
847 
848  return lwgeom_clip_to_ordinate_range(lwin, 'M', from, to, offset);
849 }
850 
851 double
852 lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
853 {
854  POINT4D p, p_proj;
855  double ret = 0.0;
856 
857  if ( ! lwin )
858  lwerror("lwgeom_interpolate_point: null input geometry!");
859 
860  if ( ! lwgeom_has_m(lwin) )
861  lwerror("Input geometry does not have a measure dimension");
862 
863  if ( lwgeom_is_empty(lwin) || lwpoint_is_empty(lwpt) )
864  lwerror("Input geometry is empty");
865 
866  switch ( lwin->type )
867  {
868  case LINETYPE:
869  {
870  LWLINE *lwline = lwgeom_as_lwline(lwin);
871  lwpoint_getPoint4d_p(lwpt, &p);
872  ret = ptarray_locate_point(lwline->points, &p, NULL, &p_proj);
873  ret = p_proj.m;
874  break;
875  }
876  default:
877  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
878  }
879  return ret;
880 }
881 
882 /*
883  * Time of closest point of approach
884  *
885  * Given two vectors (p1-p2 and q1-q2) and
886  * a time range (t1-t2) return the time in which
887  * a point p is closest to a point q on their
888  * respective vectors, and the actual points
889  *
890  * Here we use algorithm from softsurfer.com
891  * that can be found here
892  * http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
893  *
894  * @param p0 start of first segment, will be set to actual
895  * closest point of approach on segment.
896  * @param p1 end of first segment
897  * @param q0 start of second segment, will be set to actual
898  * closest point of approach on segment.
899  * @param q1 end of second segment
900  * @param t0 start of travel time
901  * @param t1 end of travel time
902  *
903  * @return time of closest point of approach
904  *
905  */
906 static double
908  POINT4D* q0, const POINT4D* q1,
909  double t0, double t1)
910 {
911  POINT3DZ pv; /* velocity of p, aka u */
912  POINT3DZ qv; /* velocity of q, aka v */
913  POINT3DZ dv; /* velocity difference */
914  POINT3DZ w0; /* vector between first points */
915 
916  /*
917  lwnotice("FROM %g,%g,%g,%g -- %g,%g,%g,%g",
918  p0->x, p0->y, p0->z, p0->m,
919  p1->x, p1->y, p1->z, p1->m);
920  lwnotice(" TO %g,%g,%g,%g -- %g,%g,%g,%g",
921  q0->x, q0->y, q0->z, q0->m,
922  q1->x, q1->y, q1->z, q1->m);
923  */
924 
925  /* PV aka U */
926  pv.x = ( p1->x - p0->x );
927  pv.y = ( p1->y - p0->y );
928  pv.z = ( p1->z - p0->z );
929  /*lwnotice("PV: %g, %g, %g", pv.x, pv.y, pv.z);*/
930 
931  /* QV aka V */
932  qv.x = ( q1->x - q0->x );
933  qv.y = ( q1->y - q0->y );
934  qv.z = ( q1->z - q0->z );
935  /*lwnotice("QV: %g, %g, %g", qv.x, qv.y, qv.z);*/
936 
937  dv.x = pv.x - qv.x;
938  dv.y = pv.y - qv.y;
939  dv.z = pv.z - qv.z;
940  /*lwnotice("DV: %g, %g, %g", dv.x, dv.y, dv.z);*/
941 
942  double dv2 = DOT(dv,dv);
943  /*lwnotice("DOT: %g", dv2);*/
944 
945  if ( dv2 == 0.0 )
946  {
947  /* Distance is the same at any time, we pick the earliest */
948  return t0;
949  }
950 
951  /* Distance at any given time, with t0 */
952  w0.x = ( p0->x - q0->x );
953  w0.y = ( p0->y - q0->y );
954  w0.z = ( p0->z - q0->z );
955 
956  /*lwnotice("W0: %g, %g, %g", w0.x, w0.y, w0.z);*/
957 
958  /* Check that at distance dt w0 is distance */
959 
960  /* This is the fraction of measure difference */
961  double t = -DOT(w0,dv) / dv2;
962  /*lwnotice("CLOSEST TIME (fraction): %g", t);*/
963 
964  if ( t > 1.0 )
965  {
966  /* Getting closer as we move to the end */
967  /*lwnotice("Converging");*/
968  t = 1;
969  }
970  else if ( t < 0.0 )
971  {
972  /*lwnotice("Diverging");*/
973  t = 0;
974  }
975 
976  /* Interpolate the actual points now */
977 
978  p0->x += pv.x * t;
979  p0->y += pv.y * t;
980  p0->z += pv.z * t;
981 
982  q0->x += qv.x * t;
983  q0->y += qv.y * t;
984  q0->z += qv.z * t;
985 
986  t = t0 + (t1 - t0) * t;
987  /*lwnotice("CLOSEST TIME (real): %g", t);*/
988 
989  return t;
990 }
991 
992 static int
993 ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals)
994 {
995  POINT4D pbuf;
996  int i, n=0;
997  for (i=0; i<pa->npoints; ++i)
998  {
999  getPoint4d_p(pa, i, &pbuf); /* could be optimized */
1000  if ( pbuf.m >= tmin && pbuf.m <= tmax )
1001  mvals[n++] = pbuf.m;
1002  }
1003  return n;
1004 }
1005 
1006 static int
1007 compare_double(const void *pa, const void *pb)
1008 {
1009  double a = *((double *)pa);
1010  double b = *((double *)pb);
1011  if ( a < b )
1012  return -1;
1013  else if ( a > b )
1014  return 1;
1015  else
1016  return 0;
1017 }
1018 
1019 /* Return number of elements in unique array */
1020 static int
1021 uniq(double *vals, int nvals)
1022 {
1023  int i, last=0;
1024  for (i=1; i<nvals; ++i)
1025  {
1026  // lwnotice("(I%d):%g", i, vals[i]);
1027  if ( vals[i] != vals[last] )
1028  {
1029  vals[++last] = vals[i];
1030  // lwnotice("(O%d):%g", last, vals[last]);
1031  }
1032  }
1033  return last+1;
1034 }
1035 
1036 /*
1037  * Find point at a given measure
1038  *
1039  * The function assumes measures are linear so that always a single point
1040  * is returned for a single measure.
1041  *
1042  * @param pa the point array to perform search on
1043  * @param m the measure to search for
1044  * @param p the point to write result into
1045  * @param from the segment number to start from
1046  *
1047  * @return the segment number the point was found into
1048  * or -1 if given measure was out of the known range.
1049  */
1050 static int
1051 ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from)
1052 {
1053  int i = from;
1054  POINT4D p1, p2;
1055 
1056  /* Walk through each segment in the point array */
1057  getPoint4d_p(pa, i, &p1);
1058  for ( i = from+1; i < pa->npoints; i++ )
1059  {
1060  getPoint4d_p(pa, i, &p2);
1061 
1062  if ( segment_locate_along(&p1, &p2, m, 0, p) == LW_TRUE )
1063  return i-1; /* found */
1064 
1065  p1 = p2;
1066  }
1067 
1068  return -1; /* not found */
1069 }
1070 
1071 double
1072 lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
1073 {
1074  LWLINE *l1, *l2;
1075  int i;
1076  const GBOX *gbox1, *gbox2;
1077  double tmin, tmax;
1078  double *mvals;
1079  int nmvals = 0;
1080  double mintime;
1081  double mindist2 = FLT_MAX; /* minimum distance, squared */
1082 
1083  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1084  {
1085  lwerror("Both input geometries must have a measure dimension");
1086  return -1;
1087  }
1088 
1089  l1 = lwgeom_as_lwline(g1);
1090  l2 = lwgeom_as_lwline(g2);
1091 
1092  if ( ! l1 || ! l2 )
1093  {
1094  lwerror("Both input geometries must be linestrings");
1095  return -1;
1096  }
1097 
1098  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1099  {
1100  lwerror("Both input lines must have at least 2 points");
1101  return -1;
1102  }
1103 
1104  /* WARNING: these ranges may be wider than real ones */
1105  gbox1 = lwgeom_get_bbox(g1);
1106  gbox2 = lwgeom_get_bbox(g2);
1107 
1108  assert(gbox1); /* or the npoints check above would have failed */
1109  assert(gbox2); /* or the npoints check above would have failed */
1110 
1111  /*
1112  * Find overlapping M range
1113  * WARNING: may be larger than the real one
1114  */
1115 
1116  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1117  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1118 
1119  if ( tmax < tmin )
1120  {
1121  LWDEBUG(1, "Inputs never exist at the same time");
1122  return -2;
1123  }
1124 
1125  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1126 
1127  /*
1128  * Collect M values in common time range from inputs
1129  */
1130 
1131  mvals = lwalloc( sizeof(double) *
1132  ( l1->points->npoints + l2->points->npoints ) );
1133 
1134  /* TODO: also clip the lines ? */
1135  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1136  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1137 
1138  /* Sort values in ascending order */
1139  qsort(mvals, nmvals, sizeof(double), compare_double);
1140 
1141  /* Remove duplicated values */
1142  nmvals = uniq(mvals, nmvals);
1143 
1144  if ( nmvals < 2 )
1145  {
1146  {
1147  /* there's a single time, must be that one... */
1148  double t0 = mvals[0];
1149  POINT4D p0, p1;
1150  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1151  if ( mindist )
1152  {
1153  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1154  {
1155  lwfree(mvals);
1156  lwerror("Could not find point with M=%g on first geom", t0);
1157  return -1;
1158  }
1159  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1160  {
1161  lwfree(mvals);
1162  lwerror("Could not find point with M=%g on second geom", t0);
1163  return -1;
1164  }
1165  *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1);
1166  }
1167  lwfree(mvals);
1168  return t0;
1169  }
1170  }
1171 
1172  /*
1173  * For each consecutive pair of measures, compute time of closest point
1174  * approach and actual distance between points at that time
1175  */
1176  mintime = tmin;
1177  for (i=1; i<nmvals; ++i)
1178  {
1179  double t0 = mvals[i-1];
1180  double t1 = mvals[i];
1181  double t;
1182  POINT4D p0, p1, q0, q1;
1183  int seg;
1184  double dist2;
1185 
1186  // lwnotice("T %g-%g", t0, t1);
1187 
1188  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1189  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1190  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1191 
1192  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1193  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1194  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1195 
1196  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1197  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1198  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1199 
1200  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1201  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1202  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1203 
1204  t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1205 
1206  /*
1207  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1208  p0.x, p0.y, p0.z,
1209  q0.x, q0.y, q0.z, t);
1210  */
1211 
1212  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1213  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1214  ( q0.z - p0.z ) * ( q0.z - p0.z );
1215  if ( dist2 < mindist2 )
1216  {
1217  mindist2 = dist2;
1218  mintime = t;
1219  // lwnotice("MINTIME: %g", mintime);
1220  }
1221  }
1222 
1223  /*
1224  * Release memory
1225  */
1226 
1227  lwfree(mvals);
1228 
1229  if ( mindist )
1230  {
1231  *mindist = sqrt(mindist2);
1232  }
1233  /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
1234 
1235  return mintime;
1236 }
1237 
1238 int
1239 lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
1240 {
1241  LWLINE *l1, *l2;
1242  int i;
1243  const GBOX *gbox1, *gbox2;
1244  double tmin, tmax;
1245  double *mvals;
1246  int nmvals = 0;
1247  double maxdist2 = maxdist * maxdist;
1248  int within = LW_FALSE;
1249 
1250  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1251  {
1252  lwerror("Both input geometries must have a measure dimension");
1253  return LW_FALSE;
1254  }
1255 
1256  l1 = lwgeom_as_lwline(g1);
1257  l2 = lwgeom_as_lwline(g2);
1258 
1259  if ( ! l1 || ! l2 )
1260  {
1261  lwerror("Both input geometries must be linestrings");
1262  return LW_FALSE;
1263  }
1264 
1265  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1266  {
1267  /* TODO: return distance between these two points */
1268  lwerror("Both input lines must have at least 2 points");
1269  return LW_FALSE;
1270  }
1271 
1272  /* WARNING: these ranges may be wider than real ones */
1273  gbox1 = lwgeom_get_bbox(g1);
1274  gbox2 = lwgeom_get_bbox(g2);
1275 
1276  assert(gbox1); /* or the npoints check above would have failed */
1277  assert(gbox2); /* or the npoints check above would have failed */
1278 
1279  /*
1280  * Find overlapping M range
1281  * WARNING: may be larger than the real one
1282  */
1283 
1284  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1285  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1286 
1287  if ( tmax < tmin )
1288  {
1289  LWDEBUG(1, "Inputs never exist at the same time");
1290  return LW_FALSE;
1291  }
1292 
1293  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1294 
1295  /*
1296  * Collect M values in common time range from inputs
1297  */
1298 
1299  mvals = lwalloc( sizeof(double) *
1300  ( l1->points->npoints + l2->points->npoints ) );
1301 
1302  /* TODO: also clip the lines ? */
1303  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1304  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1305 
1306  /* Sort values in ascending order */
1307  qsort(mvals, nmvals, sizeof(double), compare_double);
1308 
1309  /* Remove duplicated values */
1310  nmvals = uniq(mvals, nmvals);
1311 
1312  if ( nmvals < 2 )
1313  {
1314  /* there's a single time, must be that one... */
1315  double t0 = mvals[0];
1316  POINT4D p0, p1;
1317  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1318  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1319  {
1320  lwnotice("Could not find point with M=%g on first geom", t0);
1321  return LW_FALSE;
1322  }
1323  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1324  {
1325  lwnotice("Could not find point with M=%g on second geom", t0);
1326  return LW_FALSE;
1327  }
1328  if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist )
1329  within = LW_TRUE;
1330  lwfree(mvals);
1331  return within;
1332  }
1333 
1334  /*
1335  * For each consecutive pair of measures, compute time of closest point
1336  * approach and actual distance between points at that time
1337  */
1338  for (i=1; i<nmvals; ++i)
1339  {
1340  double t0 = mvals[i-1];
1341  double t1 = mvals[i];
1342 #if POSTGIS_DEBUG_LEVEL >= 1
1343  double t;
1344 #endif
1345  POINT4D p0, p1, q0, q1;
1346  int seg;
1347  double dist2;
1348 
1349  // lwnotice("T %g-%g", t0, t1);
1350 
1351  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1352  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1353  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1354 
1355  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1356  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1357  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1358 
1359  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1360  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1361  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1362 
1363  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1364  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1365  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1366 
1367 #if POSTGIS_DEBUG_LEVEL >= 1
1368  t =
1369 #endif
1370  segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1371 
1372  /*
1373  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1374  p0.x, p0.y, p0.z,
1375  q0.x, q0.y, q0.z, t);
1376  */
1377 
1378  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1379  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1380  ( q0.z - p0.z ) * ( q0.z - p0.z );
1381  if ( dist2 <= maxdist2 )
1382  {
1383  LWDEBUGF(1, "Within distance %g at time %g, breaking", sqrt(dist2), t);
1384  within = LW_TRUE;
1385  break;
1386  }
1387  }
1388 
1389  /*
1390  * Release memory
1391  */
1392 
1393  lwfree(mvals);
1394 
1395  return within;
1396 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
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:396
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
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:89
void lwfree(void *mem)
Definition: lwutil.c:242
static LWMPOINT * lwline_locate_along(const LWLINE *lwline, double m, double offset)
int npoints
Definition: liblwgeom.h:370
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:482
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:330
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
GBOX * bbox
Definition: liblwgeom.h:504
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:835
LWGEOM * lwmline_as_lwgeom(const LWMLINE *obj)
Definition: lwgeom.c:222
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)
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:849
static POINTARRAY * ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:599
int32_t srid
Definition: liblwgeom.h:398
int ngeoms
Definition: liblwgeom.h:480
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
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:1303
#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:216
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
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
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:221
#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)
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:653
static int ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from)
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:102
int32_t srid
Definition: liblwgeom.h:409
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.
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:95
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:795
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:235
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int lwpoint_is_empty(const LWPOINT *point)
Definition: lwpoint.c:273
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:612
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:267
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:111
void * lwalloc(size_t size)
Definition: lwutil.c:227
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:1310
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
int ngeoms
Definition: liblwgeom.h:467
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:856
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:102
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:244
#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:232
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:227