PostGIS  2.3.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 550 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().

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