PostGIS  2.5.0dev-r@@SVN_REVISION@@
LWCOLLECTION* lwline_clip_to_ordinate_range ( const LWLINE line,
char  ordinate,
double  from,
double  to 
)

Clip a line based on the from/to range of one of its ordinates.

Use for m- and z- clipping

Clip a line based on the from/to range of one of its ordinates.

Definition at line 548 of file lwlinearreferencing.c.

References LWCOLLECTION::bbox, COLLECTIONTYPE, LWLINE::flags, FLAGS_NDIMS, getPoint4d_p(), LW_FALSE, lwalloc(), lwcollection_add_lwgeom(), lwcollection_construct_empty(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_refresh_bbox(), lwgeom_to_ewkt(), lwline_as_lwgeom(), lwline_construct(), lwpoint_as_lwgeom(), lwpoint_construct(), lwpoint_get_ordinate(), MULTILINETYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, point_interpolate(), LWLINE::points, ptarray_append_point(), ptarray_construct(), ptarray_construct_empty(), ptarray_set_point4d(), r, LWLINE::srid, and LWCOLLECTION::type.

Referenced by lwgeom_clip_to_ordinate_range(), lwmline_clip_to_ordinate_range(), test_lwline_clip(), test_lwline_clip_big(), and test_lwmline_clip().

549 {
550 
551  POINTARRAY *pa_in = NULL;
552  LWCOLLECTION *lwgeom_out = NULL;
553  POINTARRAY *dp = NULL;
554  uint32_t i;
555  int added_last_point = 0;
556  POINT4D *p = NULL, *q = NULL, *r = NULL;
557  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
558  char hasz = lwgeom_has_z(lwline_as_lwgeom(line));
559  char hasm = lwgeom_has_m(lwline_as_lwgeom(line));
560  char dims = FLAGS_NDIMS(line->flags);
561 
562  /* Null input, nothing we can do. */
563  if ( ! line )
564  {
565  lwerror("Null input geometry.");
566  return NULL;
567  }
568 
569  /* Ensure 'from' is less than 'to'. */
570  if ( to < from )
571  {
572  double t = from;
573  from = to;
574  to = t;
575  }
576 
577  LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
578  LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
579 
580  /* Asking for an ordinate we don't have. Error. */
581  if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
582  {
583  lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
584  return NULL;
585  }
586 
587  /* Prepare our working point objects. */
588  p = lwalloc(sizeof(POINT4D));
589  q = lwalloc(sizeof(POINT4D));
590  r = lwalloc(sizeof(POINT4D));
591 
592  /* Construct a collection to hold our outputs. */
593  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
594 
595  /* Get our input point array */
596  pa_in = line->points;
597 
598  for ( i = 0; i < pa_in->npoints; i++ )
599  {
600  LWDEBUGF(4, "Point #%d", i);
601  LWDEBUGF(4, "added_last_point %d", added_last_point);
602  if ( i > 0 )
603  {
604  *q = *p;
605  ordinate_value_q = ordinate_value_p;
606  }
607  getPoint4d_p(pa_in, i, p);
608  ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
609  LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
610  LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
611 
612  /* Is this point inside the ordinate range? Yes. */
613  if ( ordinate_value_p >= from && ordinate_value_p <= to )
614  {
615  LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
616 
617  if ( ! added_last_point )
618  {
619  LWDEBUG(4," new ptarray required");
620  /* We didn't add the previous point, so this is a new segment.
621  * Make a new point array. */
622  dp = ptarray_construct_empty(hasz, hasm, 32);
623 
624  /* We're transiting into the range so add an interpolated
625  * point at the range boundary.
626  * If we're on a boundary and crossing from the far side,
627  * we also need an interpolated point. */
628  if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
629  ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
630  ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
631  ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
632  {
633  double interpolation_value;
634  (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
635  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
637  LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
638  }
639  }
640  /* Add the current vertex to the point array. */
642  if ( ordinate_value_p == from || ordinate_value_p == to )
643  {
644  added_last_point = 2; /* Added on boundary. */
645  }
646  else
647  {
648  added_last_point = 1; /* Added inside range. */
649  }
650  }
651  /* Is this point inside the ordinate range? No. */
652  else
653  {
654  LWDEBUGF(4, " added_last_point (%d)", added_last_point);
655  if ( added_last_point == 1 )
656  {
657  /* We're transiting out of the range, so add an interpolated point
658  * to the point array at the range boundary. */
659  double interpolation_value;
660  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
661  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
663  LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
664  }
665  else if ( added_last_point == 2 )
666  {
667  /* We're out and the last point was on the boundary.
668  * If the last point was the near boundary, nothing to do.
669  * If it was the far boundary, we need an interpolated point. */
670  if ( from != to && (
671  (ordinate_value_q == from && ordinate_value_p > from) ||
672  (ordinate_value_q == to && ordinate_value_p < to) ) )
673  {
674  double interpolation_value;
675  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
676  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
678  LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
679  }
680  }
681  else if ( i && ordinate_value_q < from && ordinate_value_p > to )
682  {
683  /* We just hopped over the whole range, from bottom to top,
684  * so we need to add *two* interpolated points! */
685  dp = ptarray_construct(hasz, hasm, 2);
686  /* Interpolate lower point. */
687  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
688  ptarray_set_point4d(dp, 0, r);
689  /* Interpolate upper point. */
690  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
691  ptarray_set_point4d(dp, 1, r);
692  }
693  else if ( i && ordinate_value_q > to && ordinate_value_p < from )
694  {
695  /* We just hopped over the whole range, from top to bottom,
696  * so we need to add *two* interpolated points! */
697  dp = ptarray_construct(hasz, hasm, 2);
698  /* Interpolate upper point. */
699  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
700  ptarray_set_point4d(dp, 0, r);
701  /* Interpolate lower point. */
702  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
703  ptarray_set_point4d(dp, 1, r);
704  }
705  /* We have an extant point-array, save it out to a multi-line. */
706  if ( dp )
707  {
708  LWDEBUG(4, "saving pointarray to multi-line (1)");
709 
710  /* Only one point, so we have to make an lwpoint to hold this
711  * and set the overall output type to a generic collection. */
712  if ( dp->npoints == 1 )
713  {
714  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
715  lwgeom_out->type = COLLECTIONTYPE;
716  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
717 
718  }
719  else
720  {
721  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
722  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
723  }
724 
725  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
726  dp = NULL;
727  }
728  added_last_point = 0;
729 
730  }
731  }
732 
733  /* Still some points left to be saved out. */
734  if ( dp && dp->npoints > 0 )
735  {
736  LWDEBUG(4, "saving pointarray to multi-line (2)");
737  LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
738  LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
739 
740  if ( dp->npoints == 1 )
741  {
742  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
743  lwgeom_out->type = COLLECTIONTYPE;
744  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
745  }
746  else
747  {
748  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
749  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
750  }
751 
752  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
753  dp = NULL;
754  }
755 
756  lwfree(p);
757  lwfree(q);
758  lwfree(r);
759 
760  if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
761  {
762  lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
763  }
764 
765  return lwgeom_out;
766 
767 }
void lwgeom_refresh_bbox(LWGEOM *lwgeom)
Drop current bbox and calculate a fresh one.
Definition: lwgeom.c:691
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.
char * r
Definition: cu_in_wkt.c:24
void lwfree(void *mem)
Definition: lwutil.c:244
uint8_t type
Definition: liblwgeom.h:502
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:549
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:428
GBOX * bbox
Definition: liblwgeom.h:504
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
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:420
uint32_t ngeoms
Definition: liblwgeom.h:506
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
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:329
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
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:129
void * lwalloc(size_t size)
Definition: lwutil.c:229
#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:418
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:930
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
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: