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