PostGIS  3.0.0dev-r@@SVN_REVISION@@

◆ lwline_interpolate_points()

POINTARRAY* lwline_interpolate_points ( const LWLINE line,
double  length_fraction,
char  repeat 
)

Interpolate one or more points along a line.

Definition at line 534 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().

534  {
535  POINT4D pt;
536  uint32_t i;
537  uint32_t points_to_interpolate;
538  uint32_t points_found = 0;
539  double length;
540  double length_fraction_increment = length_fraction;
541  double length_fraction_consumed = 0;
542  char has_z = (char) lwgeom_has_z(lwline_as_lwgeom(line));
543  char has_m = (char) lwgeom_has_m(lwline_as_lwgeom(line));
544  const POINTARRAY* ipa = line->points;
545  POINTARRAY* opa;
546 
547  /* Empty.InterpolatePoint == Point Empty */
548  if ( lwline_is_empty(line) )
549  {
550  return ptarray_construct_empty(has_z, has_m, 0);
551  }
552 
553  /* If distance is one of the two extremes, return the point on that
554  * end rather than doing any computations
555  */
556  if ( length_fraction == 0.0 || length_fraction == 1.0 )
557  {
558  if ( length_fraction == 0.0 )
559  getPoint4d_p(ipa, 0, &pt);
560  else
561  getPoint4d_p(ipa, ipa->npoints-1, &pt);
562 
563  opa = ptarray_construct(has_z, has_m, 1);
564  ptarray_set_point4d(opa, 0, &pt);
565 
566  return opa;
567  }
568 
569  /* Interpolate points along the line */
570  length = ptarray_length_2d(ipa);
571  points_to_interpolate = repeat ? (uint32_t) floor(1 / length_fraction) : 1;
572  opa = ptarray_construct(has_z, has_m, points_to_interpolate);
573 
574  const POINT2D* p1 = getPoint2d_cp(ipa, 0);
575  for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
576  {
577  const POINT2D* p2 = getPoint2d_cp(ipa, i+1);
578  double segment_length_frac = distance2d_pt_pt(p1, p2) / length;
579 
580  /* If our target distance is before the total length we've seen
581  * so far. create a new point some distance down the current
582  * segment.
583  */
584  while ( length_fraction < length_fraction_consumed + segment_length_frac && points_found < points_to_interpolate )
585  {
586  POINT4D p1_4d = getPoint4d(ipa, i);
587  POINT4D p2_4d = getPoint4d(ipa, i+1);
588 
589  double segment_fraction = (length_fraction - length_fraction_consumed) / segment_length_frac;
590  interpolate_point4d(&p1_4d, &p2_4d, &pt, segment_fraction);
591  ptarray_set_point4d(opa, points_found++, &pt);
592  length_fraction += length_fraction_increment;
593  }
594 
595  length_fraction_consumed += segment_length_frac;
596 
597  p1 = p2;
598  }
599 
600  /* Return the last point on the line. This shouldn't happen, but
601  * could if there's some floating point rounding errors. */
602  if (points_found < points_to_interpolate) {
603  getPoint4d_p(ipa, ipa->npoints - 1, &pt);
604  ptarray_set_point4d(opa, points_found, &pt);
605  }
606 
607  return opa;
608 }
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:50
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:58
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it&#39;s 3d)
Definition: ptarray.c:1676
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:404
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2313
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:683
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:914
static const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwinline.h:91
unsigned int uint32_t
Definition: uthash.h:78
POINT4D getPoint4d(const POINTARRAY *pa, uint32_t n)
Definition: lwgeom_api.c:95
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:321
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:112
int lwline_is_empty(const LWLINE *line)
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:921
POINTARRAY * points
Definition: liblwgeom.h:424
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: