PostGIS  2.5.0dev-r@@SVN_REVISION@@
POINTARRAY* lwline_interpolate_points ( const LWLINE line,
double  length_fraction,
char  repeat 
)

Interpolate one or more points along a line.

Definition at line 542 of file lwline.c.

References distance2d_pt_pt(), getPoint2d_cp(), getPoint4d(), getPoint4d_p(), interpolate_point4d(), lwgeom_has_m(), lwgeom_has_z(), lwline_as_lwgeom(), lwline_is_empty(), POINTARRAY::npoints, LWLINE::points, ptarray_construct(), ptarray_construct_empty(), ptarray_length_2d(), and ptarray_set_point4d().

Referenced by LWGEOM_line_interpolate_point(), and test_lwline_interpolate_points().

542  {
543  POINT4D pt;
544  uint32_t i;
545  uint32_t points_to_interpolate;
546  uint32_t points_found = 0;
547  double length;
548  double length_fraction_increment = length_fraction;
549  double length_fraction_consumed = 0;
550  char has_z = (char) lwgeom_has_z(lwline_as_lwgeom(line));
551  char has_m = (char) lwgeom_has_m(lwline_as_lwgeom(line));
552  const POINTARRAY* ipa = line->points;
553  POINTARRAY* opa;
554 
555  /* Empty.InterpolatePoint == Point Empty */
556  if ( lwline_is_empty(line) )
557  {
558  return ptarray_construct_empty(has_z, has_m, 0);
559  }
560 
561  /* If distance is one of the two extremes, return the point on that
562  * end rather than doing any computations
563  */
564  if ( length_fraction == 0.0 || length_fraction == 1.0 )
565  {
566  if ( length_fraction == 0.0 )
567  getPoint4d_p(ipa, 0, &pt);
568  else
569  getPoint4d_p(ipa, ipa->npoints-1, &pt);
570 
571  opa = ptarray_construct(has_z, has_m, 1);
572  ptarray_set_point4d(opa, 0, &pt);
573 
574  return opa;
575  }
576 
577  /* Interpolate points along the line */
578  length = ptarray_length_2d(ipa);
579  points_to_interpolate = repeat ? (uint32_t) floor(1 / length_fraction) : 1;
580  opa = ptarray_construct(has_z, has_m, points_to_interpolate);
581 
582  const POINT2D* p1 = getPoint2d_cp(ipa, 0);
583  for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
584  {
585  const POINT2D* p2 = getPoint2d_cp(ipa, i+1);
586  double segment_length_frac = distance2d_pt_pt(p1, p2) / length;
587 
588  /* If our target distance is before the total length we've seen
589  * so far. create a new point some distance down the current
590  * segment.
591  */
592  while ( length_fraction < length_fraction_consumed + segment_length_frac && points_found < points_to_interpolate )
593  {
594  POINT4D p1_4d = getPoint4d(ipa, i);
595  POINT4D p2_4d = getPoint4d(ipa, i+1);
596 
597  double segment_fraction = (length_fraction - length_fraction_consumed) / segment_length_frac;
598  interpolate_point4d(&p1_4d, &p2_4d, &pt, segment_fraction);
599  ptarray_set_point4d(opa, points_found++, &pt);
600  length_fraction += length_fraction_increment;
601  }
602 
603  length_fraction_consumed += segment_length_frac;
604 
605  p1 = p2;
606  }
607 
608  /* Return the last point on the line. This shouldn't happen, but
609  * could if there's some floating point rounding errors. */
610  if (points_found < points_to_interpolate) {
611  getPoint4d_p(ipa, ipa->npoints - 1, &pt);
612  ptarray_set_point4d(opa, points_found, &pt);
613  }
614 
615  return opa;
616 }
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
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it's 3d)
Definition: ptarray.c:1682
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:425
int lwline_is_empty(const LWLINE *line)
Definition: lwline.c:511
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2312
void interpolate_point4d(const POINT4D *A, const POINT4D *B, POINT4D *I, double F)
Find interpolation point I between point A and point B so that the len(AI) == len(AB)*F and I falls o...
Definition: lwgeom_api.c:704
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:923
unsigned int uint32_t
Definition: uthash.h:78
POINT4D getPoint4d(const POINTARRAY *pa, uint32_t n)
Definition: lwgeom_api.c:96
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:329
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:930
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
POINTARRAY * points
Definition: liblwgeom.h:421
uint32_t npoints
Definition: liblwgeom.h:370

Here is the call graph for this function:

Here is the caller graph for this function: