PostGIS  2.2.8dev-r@@SVN_REVISION@@

◆ lwline_clip_to_ordinate_range()

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 537 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_add_bbox(), lwgeom_drop_bbox(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_to_ewkt(), lwline_as_lwgeom(), lwline_construct(), lwpoint_as_lwgeom(), lwpoint_construct(), lwpoint_get_ordinate(), MULTILINETYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, point_interpolate(), ptarray_append_point(), ptarray_construct(), ptarray_construct_empty(), ptarray_set_point4d(), r, LWGEOM::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().

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 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:536
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:214
int npoints
Definition: liblwgeom.h:355
uint8_t type
Definition: liblwgeom.h:487
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
GBOX * bbox
Definition: liblwgeom.h:489
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
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
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:836
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:586
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
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:599
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:98
void * lwalloc(size_t size)
Definition: lwutil.c:199
#define MULTILINETYPE
Definition: liblwgeom.h:74
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
uint8_t flags
Definition: liblwgeom.h:403
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:136
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
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
POINTARRAY * points
Definition: liblwgeom.h:406
Here is the call graph for this function:
Here is the caller graph for this function: