PostGIS  2.3.8dev-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 #include "liblwgeom_internal.h"
27 #include "lwgeom_log.h"
28 #include "measures3d.h"
29 
30 static int
31 segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
32 {
33  double m1 = p1->m;
34  double m2 = p2->m;
35  double mprop;
36 
37  /* M is out of range, no new point generated. */
38  if ( (m < FP_MIN(m1,m2)) || (m > FP_MAX(m1,m2)) )
39  {
40  return LW_FALSE;
41  }
42 
43  if( m1 == m2 )
44  {
45  /* Degenerate case: same M on both points.
46  If they are the same point we just return one of them. */
47  if ( p4d_same(p1,p2) )
48  {
49  *pn = *p1;
50  return LW_TRUE;
51  }
52  /* If the points are different we split the difference */
53  mprop = 0.5;
54  }
55  else
56  {
57  mprop = (m - m1) / (m2 - m1);
58  }
59 
60  /* M is in range, new point to be generated. */
61  pn->x = p1->x + (p2->x - p1->x) * mprop;
62  pn->y = p1->y + (p2->y - p1->y) * mprop;
63  pn->z = p1->z + (p2->z - p1->z) * mprop;
64  pn->m = m;
65 
66  /* Offset to the left or right, if necessary. */
67  if ( offset != 0.0 )
68  {
69  double theta = atan2(p2->y - p1->y, p2->x - p1->x);
70  pn->x -= sin(theta) * offset;
71  pn->y += cos(theta) * offset;
72  }
73 
74  return LW_TRUE;
75 }
76 
77 
78 static POINTARRAY*
79 ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
80 {
81  int i;
82  POINT4D p1, p2, pn;
83  POINTARRAY *dpa = NULL;
84 
85  /* Can't do anything with degenerate point arrays */
86  if ( ! pa || pa->npoints < 2 ) return NULL;
87 
88  /* Walk through each segment in the point array */
89  for ( i = 1; i < pa->npoints; i++ )
90  {
91  getPoint4d_p(pa, i-1, &p1);
92  getPoint4d_p(pa, i, &p2);
93 
94  /* No derived point? Move to next segment. */
95  if ( segment_locate_along(&p1, &p2, m, offset, &pn) == LW_FALSE )
96  continue;
97 
98  /* No pointarray, make a fresh one */
99  if ( dpa == NULL )
101 
102  /* Add our new point to the array */
103  ptarray_append_point(dpa, &pn, 0);
104  }
105 
106  return dpa;
107 }
108 
109 static LWMPOINT*
110 lwline_locate_along(const LWLINE *lwline, double m, double offset)
111 {
112  POINTARRAY *opa = NULL;
113  LWMPOINT *mp = NULL;
114  LWGEOM *lwg = lwline_as_lwgeom(lwline);
115  int hasz, hasm, srid;
116 
117  /* Return degenerates upwards */
118  if ( ! lwline ) return NULL;
119 
120  /* Create empty return shell */
121  srid = lwgeom_get_srid(lwg);
122  hasz = lwgeom_has_z(lwg);
123  hasm = lwgeom_has_m(lwg);
124 
125  if ( hasm )
126  {
127  /* Find points along */
128  opa = ptarray_locate_along(lwline->points, m, offset);
129  }
130  else
131  {
132  LWLINE *lwline_measured = lwline_measured_from_lwline(lwline, 0.0, 1.0);
133  opa = ptarray_locate_along(lwline_measured->points, m, offset);
134  lwline_free(lwline_measured);
135  }
136 
137  /* Return NULL as EMPTY */
138  if ( ! opa )
139  return lwmpoint_construct_empty(srid, hasz, hasm);
140 
141  /* Convert pointarray into a multipoint */
142  mp = lwmpoint_construct(srid, opa);
143  ptarray_free(opa);
144  return mp;
145 }
146 
147 static LWMPOINT*
148 lwmline_locate_along(const LWMLINE *lwmline, double m, double offset)
149 {
150  LWMPOINT *lwmpoint = NULL;
151  LWGEOM *lwg = lwmline_as_lwgeom(lwmline);
152  int i, j;
153 
154  /* Return degenerates upwards */
155  if ( (!lwmline) || (lwmline->ngeoms < 1) ) return NULL;
156 
157  /* Construct return */
159 
160  /* Locate along each sub-line */
161  for ( i = 0; i < lwmline->ngeoms; i++ )
162  {
163  LWMPOINT *along = lwline_locate_along(lwmline->geoms[i], m, offset);
164  if ( along )
165  {
166  if ( ! lwgeom_is_empty((LWGEOM*)along) )
167  {
168  for ( j = 0; j < along->ngeoms; j++ )
169  {
170  lwmpoint_add_lwpoint(lwmpoint, along->geoms[j]);
171  }
172  }
173  /* Free the containing geometry, but leave the sub-geometries around */
174  along->ngeoms = 0;
175  lwmpoint_free(along);
176  }
177  }
178  return lwmpoint;
179 }
180 
181 static LWMPOINT*
182 lwpoint_locate_along(const LWPOINT *lwpoint, double m, double offset)
183 {
184  double point_m = lwpoint_get_m(lwpoint);
185  LWGEOM *lwg = lwpoint_as_lwgeom(lwpoint);
187  if ( FP_EQUALS(m, point_m) )
188  {
189  lwmpoint_add_lwpoint(r, lwpoint_clone(lwpoint));
190  }
191  return r;
192 }
193 
194 static LWMPOINT*
195 lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset)
196 {
197  LWGEOM *lwg = lwmpoint_as_lwgeom(lwin);
198  LWMPOINT *lwout = NULL;
199  int i;
200 
201  /* Construct return */
203 
204  for ( i = 0; i < lwin->ngeoms; i++ )
205  {
206  double point_m = lwpoint_get_m(lwin->geoms[i]);
207  if ( FP_EQUALS(m, point_m) )
208  {
209  lwmpoint_add_lwpoint(lwout, lwpoint_clone(lwin->geoms[i]));
210  }
211  }
212 
213  return lwout;
214 }
215 
216 LWGEOM*
217 lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
218 {
219  if ( ! lwin ) return NULL;
220 
221  if ( ! lwgeom_has_m(lwin) )
222  lwerror("Input geometry does not have a measure dimension");
223 
224  switch (lwin->type)
225  {
226  case POINTTYPE:
227  return (LWGEOM*)lwpoint_locate_along((LWPOINT*)lwin, m, offset);
228  case MULTIPOINTTYPE:
229  return (LWGEOM*)lwmpoint_locate_along((LWMPOINT*)lwin, m, offset);
230  case LINETYPE:
231  return (LWGEOM*)lwline_locate_along((LWLINE*)lwin, m, offset);
232  case MULTILINETYPE:
233  return (LWGEOM*)lwmline_locate_along((LWMLINE*)lwin, m, offset);
234  /* Only line types supported right now */
235  /* TO DO: CurveString, CompoundCurve, MultiCurve */
236  /* TO DO: Point, MultiPoint */
237  default:
238  lwerror("Only linear geometries are supported, %s provided.",lwtype_name(lwin->type));
239  return NULL;
240  }
241  return NULL;
242 }
243 
251 double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
252 {
253  if ( ! p )
254  {
255  lwerror("Null input geometry.");
256  return 0.0;
257  }
258 
259  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
260  {
261  lwerror("Cannot extract %c ordinate.", ordinate);
262  return 0.0;
263  }
264 
265  if ( ordinate == 'X' )
266  return p->x;
267  if ( ordinate == 'Y' )
268  return p->y;
269  if ( ordinate == 'Z' )
270  return p->z;
271  if ( ordinate == 'M' )
272  return p->m;
273 
274  /* X */
275  return p->x;
276 
277 }
278 
283 void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
284 {
285  if ( ! p )
286  {
287  lwerror("Null input geometry.");
288  return;
289  }
290 
291  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
292  {
293  lwerror("Cannot set %c ordinate.", ordinate);
294  return;
295  }
296 
297  LWDEBUGF(4, " setting ordinate %c to %g", ordinate, value);
298 
299  switch ( ordinate )
300  {
301  case 'X':
302  p->x = value;
303  return;
304  case 'Y':
305  p->y = value;
306  return;
307  case 'Z':
308  p->z = value;
309  return;
310  case 'M':
311  p->m = value;
312  return;
313  }
314 }
315 
321 int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
322 {
323  static char* dims = "XYZM";
324  double p1_value = lwpoint_get_ordinate(p1, ordinate);
325  double p2_value = lwpoint_get_ordinate(p2, ordinate);
326  double proportion;
327  int i = 0;
328 
329  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
330  {
331  lwerror("Cannot set %c ordinate.", ordinate);
332  return 0;
333  }
334 
335  if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
336  FP_MAX(p1_value, p2_value) < interpolation_value )
337  {
338  lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
339  return 0;
340  }
341 
342  proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
343 
344  for ( i = 0; i < 4; i++ )
345  {
346  double newordinate = 0.0;
347  if ( dims[i] == 'Z' && ! hasz ) continue;
348  if ( dims[i] == 'M' && ! hasm ) continue;
349  p1_value = lwpoint_get_ordinate(p1, dims[i]);
350  p2_value = lwpoint_get_ordinate(p2, dims[i]);
351  newordinate = p1_value + proportion * (p2_value - p1_value);
352  lwpoint_set_ordinate(p, dims[i], newordinate);
353  LWDEBUGF(4, " clip ordinate(%c) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", dims[i], p1_value, p2_value, proportion, newordinate );
354  }
355 
356  return 1;
357 }
358 
359 
364 lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
365 {
366  LWCOLLECTION *lwgeom_out = NULL;
367  char hasz, hasm;
368  POINT4D p4d;
369  double ordinate_value;
370 
371  /* Nothing to do with NULL */
372  if ( ! point )
373  lwerror("Null input geometry.");
374 
375  /* Ensure 'from' is less than 'to'. */
376  if ( to < from )
377  {
378  double t = from;
379  from = to;
380  to = t;
381  }
382 
383  /* Read Z/M info */
384  hasz = lwgeom_has_z(lwpoint_as_lwgeom(point));
385  hasm = lwgeom_has_m(lwpoint_as_lwgeom(point));
386 
387  /* Prepare return object */
388  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, point->srid, hasz, hasm);
389 
390  /* Test if ordinate is in range */
391  lwpoint_getPoint4d_p(point, &p4d);
392  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
393  if ( from <= ordinate_value && to >= ordinate_value )
394  {
395  LWPOINT *lwp = lwpoint_clone(point);
396  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
397  }
398 
399  /* Set the bbox, if necessary */
400  if ( lwgeom_out->bbox )
401  {
402  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
403  lwgeom_add_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  int 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_drop_bbox((LWGEOM*)lwgeom_out);
460  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
461  }
462 
463  return lwgeom_out;
464 }
465 
470 lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
471 {
472  LWCOLLECTION *lwgeom_out = NULL;
473 
474  if ( ! mline )
475  {
476  lwerror("Null input geometry.");
477  return NULL;
478  }
479 
480  if ( mline->ngeoms == 1)
481  {
482  lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
483  }
484  else
485  {
486  LWCOLLECTION *col;
487  char hasz = lwgeom_has_z(lwmline_as_lwgeom(mline));
488  char hasm = lwgeom_has_m(lwmline_as_lwgeom(mline));
489  int i, j;
490  char homogeneous = 1;
491  size_t geoms_size = 0;
492  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
493  FLAGS_SET_Z(lwgeom_out->flags, hasz);
494  FLAGS_SET_M(lwgeom_out->flags, hasm);
495  for ( i = 0; i < mline->ngeoms; i ++ )
496  {
497  col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
498  if ( col )
499  {
500  /* Something was left after the clip. */
501  if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
502  {
503  geoms_size += 16;
504  if ( lwgeom_out->geoms )
505  {
506  lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
507  }
508  else
509  {
510  lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
511  }
512  }
513  for ( j = 0; j < col->ngeoms; j++ )
514  {
515  lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
516  lwgeom_out->ngeoms++;
517  }
518  if ( col->type != mline->type )
519  {
520  homogeneous = 0;
521  }
522  /* Shallow free the struct, leaving the geoms behind. */
523  if ( col->bbox ) lwfree(col->bbox);
524  lwfree(col->geoms);
525  lwfree(col);
526  }
527  }
528  if ( lwgeom_out->bbox )
529  {
530  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
531  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
532  }
533 
534  if ( ! homogeneous )
535  {
536  lwgeom_out->type = COLLECTIONTYPE;
537  }
538  }
539 
540  return lwgeom_out;
541 
542 }
543 
544 
550 lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
551 {
552 
553  POINTARRAY *pa_in = NULL;
554  LWCOLLECTION *lwgeom_out = NULL;
555  POINTARRAY *dp = NULL;
556  int i;
557  int added_last_point = 0;
558  POINT4D *p = NULL, *q = NULL, *r = NULL;
559  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
560  char hasz = lwgeom_has_z(lwline_as_lwgeom(line));
561  char hasm = lwgeom_has_m(lwline_as_lwgeom(line));
562  char dims = FLAGS_NDIMS(line->flags);
563 
564  /* Null input, nothing we can do. */
565  if ( ! line )
566  {
567  lwerror("Null input geometry.");
568  return NULL;
569  }
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_drop_bbox((LWGEOM*)lwgeom_out);
765  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
766  }
767 
768  return lwgeom_out;
769 
770 }
771 
773 lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
774 {
775  LWCOLLECTION *out_col;
776  LWCOLLECTION *out_offset;
777  int i;
778 
779  if ( ! lwin )
780  lwerror("lwgeom_clip_to_ordinate_range: null input geometry!");
781 
782  switch ( lwin->type )
783  {
784  case LINETYPE:
785  out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to);
786  break;
787  case MULTILINETYPE:
788  out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to);
789  break;
790  case MULTIPOINTTYPE:
791  out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to);
792  break;
793  case POINTTYPE:
794  out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to);
795  break;
796  default:
797  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
798  return NULL;;
799  }
800 
801  /* Stop if result is NULL */
802  if ( out_col == NULL )
803  lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL");
804 
805  /* Return if we aren't going to offset the result */
806  if ( FP_EQUALS(offset, 0.0) || lwgeom_is_empty(lwcollection_as_lwgeom(out_col)) )
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  int 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
1050 ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from)
1051 {
1052  int 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  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  /* We use lwgeom_calculate_gbox() instead of lwgeom_get_gbox() */
1104  /* because we cannot afford the float rounding inaccuracy when */
1105  /* we compare the ranges for overlap below */
1106  lwgeom_calculate_gbox(g1, &gbox1);
1107  lwgeom_calculate_gbox(g2, &gbox2);
1108 
1109  /*
1110  * Find overlapping M range
1111  * WARNING: may be larger than the real one
1112  */
1113 
1114  tmin = FP_MAX(gbox1.mmin, gbox2.mmin);
1115  tmax = FP_MIN(gbox1.mmax, gbox2.mmax);
1116 
1117  if ( tmax < tmin )
1118  {
1119  LWDEBUG(1, "Inputs never exist at the same time");
1120  return -2;
1121  }
1122 
1123  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1124 
1125  /*
1126  * Collect M values in common time range from inputs
1127  */
1128 
1129  mvals = lwalloc( sizeof(double) *
1130  ( l1->points->npoints + l2->points->npoints ) );
1131 
1132  /* TODO: also clip the lines ? */
1133  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1134  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1135 
1136  /* Sort values in ascending order */
1137  qsort(mvals, nmvals, sizeof(double), compare_double);
1138 
1139  /* Remove duplicated values */
1140  nmvals = uniq(mvals, nmvals);
1141 
1142  if ( nmvals < 2 )
1143  {
1144  {
1145  /* there's a single time, must be that one... */
1146  double t0 = mvals[0];
1147  POINT4D p0, p1;
1148  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1149  if ( mindist )
1150  {
1151  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1152  {
1153  lwfree(mvals);
1154  lwerror("Could not find point with M=%g on first geom", t0);
1155  return -1;
1156  }
1157  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1158  {
1159  lwfree(mvals);
1160  lwerror("Could not find point with M=%g on second geom", t0);
1161  return -1;
1162  }
1163  *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1);
1164  }
1165  lwfree(mvals);
1166  return t0;
1167  }
1168  }
1169 
1170  /*
1171  * For each consecutive pair of measures, compute time of closest point
1172  * approach and actual distance between points at that time
1173  */
1174  mintime = tmin;
1175  for (i=1; i<nmvals; ++i)
1176  {
1177  double t0 = mvals[i-1];
1178  double t1 = mvals[i];
1179  double t;
1180  POINT4D p0, p1, q0, q1;
1181  int seg;
1182  double dist2;
1183 
1184  // lwnotice("T %g-%g", t0, t1);
1185 
1186  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1187  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1188  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1189 
1190  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1191  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1192  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1193 
1194  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1195  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1196  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1197 
1198  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1199  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1200  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1201 
1202  t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1203 
1204  /*
1205  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1206  p0.x, p0.y, p0.z,
1207  q0.x, q0.y, q0.z, t);
1208  */
1209 
1210  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1211  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1212  ( q0.z - p0.z ) * ( q0.z - p0.z );
1213  if ( dist2 < mindist2 )
1214  {
1215  mindist2 = dist2;
1216  mintime = t;
1217  // lwnotice("MINTIME: %g", mintime);
1218  }
1219  }
1220 
1221  /*
1222  * Release memory
1223  */
1224 
1225  lwfree(mvals);
1226 
1227  if ( mindist )
1228  {
1229  *mindist = sqrt(mindist2);
1230  }
1231  /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
1232 
1233  return mintime;
1234 }
1235 
1236 int
1237 lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
1238 {
1239  LWLINE *l1, *l2;
1240  int i;
1241  GBOX gbox1, gbox2;
1242  double tmin, tmax;
1243  double *mvals;
1244  int nmvals = 0;
1245  double maxdist2 = maxdist * maxdist;
1246  int within = LW_FALSE;
1247 
1248  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1249  {
1250  lwerror("Both input geometries must have a measure dimension");
1251  return LW_FALSE;
1252  }
1253 
1254  l1 = lwgeom_as_lwline(g1);
1255  l2 = lwgeom_as_lwline(g2);
1256 
1257  if ( ! l1 || ! l2 )
1258  {
1259  lwerror("Both input geometries must be linestrings");
1260  return LW_FALSE;
1261  }
1262 
1263  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1264  {
1265  /* TODO: return distance between these two points */
1266  lwerror("Both input lines must have at least 2 points");
1267  return LW_FALSE;
1268  }
1269 
1270  /* We use lwgeom_calculate_gbox() instead of lwgeom_get_gbox() */
1271  /* because we cannot afford the float rounding inaccuracy when */
1272  /* we compare the ranges for overlap below */
1273  lwgeom_calculate_gbox(g1, &gbox1);
1274  lwgeom_calculate_gbox(g2, &gbox2);
1275 
1276  /*
1277  * Find overlapping M range
1278  * WARNING: may be larger than the real one
1279  */
1280 
1281  tmin = FP_MAX(gbox1.mmin, gbox2.mmin);
1282  tmax = FP_MIN(gbox1.mmax, gbox2.mmax);
1283 
1284  if ( tmax < tmin )
1285  {
1286  LWDEBUG(1, "Inputs never exist at the same time");
1287  return LW_FALSE;
1288  }
1289 
1290  /*
1291  * Collect M values in common time range from inputs
1292  */
1293 
1294  mvals = lwalloc( sizeof(double) *
1295  ( l1->points->npoints + l2->points->npoints ) );
1296 
1297  /* TODO: also clip the lines ? */
1298  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1299  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1300 
1301  /* Sort values in ascending order */
1302  qsort(mvals, nmvals, sizeof(double), compare_double);
1303 
1304  /* Remove duplicated values */
1305  nmvals = uniq(mvals, nmvals);
1306 
1307  if ( nmvals < 2 )
1308  {
1309  /* there's a single time, must be that one... */
1310  double t0 = mvals[0];
1311  POINT4D p0, p1;
1312  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1313  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1314  {
1315  lwnotice("Could not find point with M=%g on first geom", t0);
1316  return LW_FALSE;
1317  }
1318  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1319  {
1320  lwnotice("Could not find point with M=%g on second geom", t0);
1321  return LW_FALSE;
1322  }
1323  if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist )
1324  within = LW_TRUE;
1325  lwfree(mvals);
1326  return within;
1327  }
1328 
1329  /*
1330  * For each consecutive pair of measures, compute time of closest point
1331  * approach and actual distance between points at that time
1332  */
1333  for (i=1; i<nmvals; ++i)
1334  {
1335  double t0 = mvals[i-1];
1336  double t1 = mvals[i];
1337 #if POSTGIS_DEBUG_LEVEL >= 1
1338  double t;
1339 #endif
1340  POINT4D p0, p1, q0, q1;
1341  int seg;
1342  double dist2;
1343 
1344  // lwnotice("T %g-%g", t0, t1);
1345 
1346  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1347  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1348  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1349 
1350  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1351  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1352  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1353 
1354  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1355  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1356  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1357 
1358  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1359  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1360  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1361 
1362 #if POSTGIS_DEBUG_LEVEL >= 1
1363  t =
1364 #endif
1365  segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1366 
1367  /*
1368  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1369  p0.x, p0.y, p0.z,
1370  q0.x, q0.y, q0.z, t);
1371  */
1372 
1373  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1374  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1375  ( q0.z - p0.z ) * ( q0.z - p0.z );
1376  if ( dist2 <= maxdist2 )
1377  {
1378  LWDEBUGF(1, "Within distance %g at time %g, breaking", sqrt(dist2), t);
1379  within = LW_TRUE;
1380  break;
1381  }
1382  }
1383 
1384  /*
1385  * Release memory
1386  */
1387 
1388  lwfree(mvals);
1389 
1390  return within;
1391 }
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
uint8_t type
Definition: liblwgeom.h:406
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...
#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
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:665
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:1310
#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)
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
type
Definition: ovdump.py:41
#define FP_EQUALS(A, B)
double mmax
Definition: liblwgeom.h:298
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
int value
Definition: genraster.py:61
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