PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ pta_unstroke()

LWGEOM * pta_unstroke ( const POINTARRAY points,
int32_t  srid 
)

Definition at line 987 of file lwstroke.c.

988 {
989  int i = 0, j, k;
990  POINT4D a1, a2, a3, b;
991  POINT4D first, center;
992  char *edges_in_arcs;
993  int found_arc = LW_FALSE;
994  int current_arc = 1;
995  int num_edges;
996  int edge_type; /* non-zero if edge is part of an arc */
997  int start, end;
998  LWCOLLECTION *outcol;
999  /* Minimum number of edges, per quadrant, required to define an arc */
1000  const unsigned int min_quad_edges = 2;
1001 
1002  /* Die on null input */
1003  if ( ! points )
1004  lwerror("pta_unstroke called with null pointarray");
1005 
1006  /* Null on empty input? */
1007  if ( points->npoints == 0 )
1008  return NULL;
1009 
1010  /* We can't desegmentize anything shorter than four points */
1011  if ( points->npoints < 4 )
1012  {
1013  /* Return a linestring here*/
1014  lwerror("pta_unstroke needs implementation for npoints < 4");
1015  }
1016 
1017  /* Allocate our result array of vertices that are part of arcs */
1018  num_edges = points->npoints - 1;
1019  edges_in_arcs = lwalloc(num_edges + 1);
1020  memset(edges_in_arcs, 0, num_edges + 1);
1021 
1022  /* We make a candidate arc of the first two edges, */
1023  /* And then see if the next edge follows it */
1024  while( i < num_edges-2 )
1025  {
1026  unsigned int arc_edges;
1027  double num_quadrants;
1028  double angle;
1029 
1030  found_arc = LW_FALSE;
1031  /* Make candidate arc */
1032  getPoint4d_p(points, i , &a1);
1033  getPoint4d_p(points, i+1, &a2);
1034  getPoint4d_p(points, i+2, &a3);
1035  memcpy(&first, &a1, sizeof(POINT4D));
1036 
1037  for( j = i+3; j < num_edges+1; j++ )
1038  {
1039  LWDEBUGF(4, "i=%d, j=%d", i, j);
1040  getPoint4d_p(points, j, &b);
1041  /* Does this point fall on our candidate arc? */
1042  if ( pt_continues_arc(&a1, &a2, &a3, &b) )
1043  {
1044  /* Yes. Mark this edge and the two preceding it as arc components */
1045  LWDEBUGF(4, "pt_continues_arc #%d", current_arc);
1046  found_arc = LW_TRUE;
1047  for ( k = j-1; k > j-4; k-- )
1048  edges_in_arcs[k] = current_arc;
1049  }
1050  else
1051  {
1052  /* No. So we're done with this candidate arc */
1053  LWDEBUG(4, "pt_continues_arc = false");
1054  current_arc++;
1055  break;
1056  }
1057 
1058  memcpy(&a1, &a2, sizeof(POINT4D));
1059  memcpy(&a2, &a3, sizeof(POINT4D));
1060  memcpy(&a3, &b, sizeof(POINT4D));
1061  }
1062  /* Jump past all the edges that were added to the arc */
1063  if ( found_arc )
1064  {
1065  /* Check if an arc was composed by enough edges to be
1066  * really considered an arc
1067  * See http://trac.osgeo.org/postgis/ticket/2420
1068  */
1069  arc_edges = j - 1 - i;
1070  LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
1071  if ( first.x == b.x && first.y == b.y ) {
1072  LWDEBUG(4, "arc is a circle");
1073  num_quadrants = 4;
1074  }
1075  else {
1076  lw_arc_center((POINT2D*)&first, (POINT2D*)&b, (POINT2D*)&a1, (POINT2D*)&center);
1077  angle = lw_arc_angle((POINT2D*)&first, (POINT2D*)&center, (POINT2D*)&b);
1078  int p2_side = lw_segment_side((POINT2D*)&first, (POINT2D*)&a1, (POINT2D*)&b);
1079  if ( p2_side >= 0 ) angle = -angle;
1080 
1081  if ( angle < 0 ) angle = 2 * M_PI + angle;
1082  num_quadrants = ( 4 * angle ) / ( 2 * M_PI );
1083  LWDEBUGF(4, "arc angle (%g %g, %g %g, %g %g) is %g (side is %d), quadrants:%g", first.x, first.y, center.x, center.y, b.x, b.y, angle, p2_side, num_quadrants);
1084  }
1085  /* a1 is first point, b is last point */
1086  if ( arc_edges < min_quad_edges * num_quadrants ) {
1087  LWDEBUGF(4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants);
1088  for ( k = j-1; k >= i; k-- )
1089  edges_in_arcs[k] = 0;
1090  }
1091 
1092  i = j-1;
1093  }
1094  else
1095  {
1096  /* Mark this edge as a linear edge */
1097  edges_in_arcs[i] = 0;
1098  i = i+1;
1099  }
1100  }
1101 
1102 #if POSTGIS_DEBUG_LEVEL > 3
1103  {
1104  char *edgestr = lwalloc(num_edges+1);
1105  for ( i = 0; i < num_edges; i++ )
1106  {
1107  if ( edges_in_arcs[i] )
1108  edgestr[i] = 48 + edges_in_arcs[i];
1109  else
1110  edgestr[i] = '.';
1111  }
1112  edgestr[num_edges] = 0;
1113  LWDEBUGF(3, "edge pattern %s", edgestr);
1114  lwfree(edgestr);
1115  }
1116 #endif
1117 
1118  start = 0;
1119  edge_type = edges_in_arcs[0];
1120  outcol = lwcollection_construct_empty(COMPOUNDTYPE, srid, ptarray_has_z(points), ptarray_has_m(points));
1121  for( i = 1; i < num_edges; i++ )
1122  {
1123  if( edge_type != edges_in_arcs[i] )
1124  {
1125  end = i - 1;
1126  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1127  start = i;
1128  edge_type = edges_in_arcs[i];
1129  }
1130  }
1131  lwfree(edges_in_arcs); /* not needed anymore */
1132 
1133  /* Roll out last item */
1134  end = num_edges - 1;
1135  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1136 
1137  /* Strip down to singleton if only one entry */
1138  if ( outcol->ngeoms == 1 )
1139  {
1140  LWGEOM *outgeom = outcol->geoms[0];
1141  outcol->ngeoms = 0; lwcollection_free(outcol);
1142  return outgeom;
1143  }
1144  return lwcollection_as_lwgeom(outcol);
1145 }
#define LW_FALSE
Definition: liblwgeom.h:94
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:309
#define COMPOUNDTYPE
Definition: liblwgeom.h:110
void lwfree(void *mem)
Definition: lwutil.c:242
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
Definition: lwcollection.c:92
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:357
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:125
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:188
void * lwalloc(size_t size)
Definition: lwutil.c:227
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:93
double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
Determines the center of the circle defined by the three given points.
Definition: lwalgorithm.c:226
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:37
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:62
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:44
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static double lw_arc_angle(const POINT2D *a, const POINT2D *b, const POINT2D *c)
Return ABC angle in radians TODO: move to lwalgorithm.
Definition: lwstroke.c:881
static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
Returns LW_TRUE if b is on the arc formed by a1/a2/a3, but not within that portion already described ...
Definition: lwstroke.c:903
static LWGEOM * geom_from_pa(const POINTARRAY *pa, int32_t srid, int is_arc, int start, int end)
Definition: lwstroke.c:977
uint32_t ngeoms
Definition: liblwgeom.h:580
LWGEOM ** geoms
Definition: liblwgeom.h:575
double y
Definition: liblwgeom.h:390
double x
Definition: liblwgeom.h:390
double x
Definition: liblwgeom.h:414
double y
Definition: liblwgeom.h:414
uint32_t npoints
Definition: liblwgeom.h:427

References COMPOUNDTYPE, geom_from_pa(), LWCOLLECTION::geoms, getPoint4d_p(), lw_arc_angle(), lw_arc_center(), LW_FALSE, lw_segment_side(), LW_TRUE, lwalloc(), lwcollection_add_lwgeom(), lwcollection_as_lwgeom(), lwcollection_construct_empty(), lwcollection_free(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), LWCOLLECTION::ngeoms, POINTARRAY::npoints, pt_continues_arc(), ptarray_has_m(), ptarray_has_z(), POINT2D::x, POINT4D::x, POINT2D::y, and POINT4D::y.

Referenced by lwline_unstroke(), and lwpolygon_unstroke().

Here is the call graph for this function:
Here is the caller graph for this function: