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