PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ pta_unstroke()

LWGEOM * pta_unstroke ( const POINTARRAY points,
int32_t  srid 
)

Definition at line 975 of file lwstroke.c.

976 {
977  int i = 0, j, k;
978  POINT4D a1, a2, a3, b;
979  POINT4D first, center;
980  char *edges_in_arcs;
981  int found_arc = LW_FALSE;
982  int current_arc = 1;
983  int num_edges;
984  int edge_type; /* non-zero if edge is part of an arc */
985  int start, end;
986  LWCOLLECTION *outcol;
987  /* Minimum number of edges, per quadrant, required to define an arc */
988  const unsigned int min_quad_edges = 2;
989 
990  /* Die on null input */
991  if ( ! points )
992  lwerror("pta_unstroke called with null pointarray");
993 
994  /* Null on empty input? */
995  if ( points->npoints == 0 )
996  return NULL;
997 
998  /* We can't desegmentize anything shorter than four points */
999  if ( points->npoints < 4 )
1000  {
1001  /* Return a linestring here*/
1002  lwerror("pta_unstroke needs implementation for npoints < 4");
1003  }
1004 
1005  /* Allocate our result array of vertices that are part of arcs */
1006  num_edges = points->npoints - 1;
1007  edges_in_arcs = lwalloc(num_edges + 1);
1008  memset(edges_in_arcs, 0, num_edges + 1);
1009 
1010  /* We make a candidate arc of the first two edges, */
1011  /* And then see if the next edge follows it */
1012  while( i < num_edges-2 )
1013  {
1014  unsigned int arc_edges;
1015  double num_quadrants;
1016  double angle;
1017 
1018  found_arc = LW_FALSE;
1019  /* Make candidate arc */
1020  getPoint4d_p(points, i , &a1);
1021  getPoint4d_p(points, i+1, &a2);
1022  getPoint4d_p(points, i+2, &a3);
1023  memcpy(&first, &a1, sizeof(POINT4D));
1024 
1025  for( j = i+3; j < num_edges+1; j++ )
1026  {
1027  LWDEBUGF(4, "i=%d, j=%d", i, j);
1028  getPoint4d_p(points, j, &b);
1029  /* Does this point fall on our candidate arc? */
1030  if ( pt_continues_arc(&a1, &a2, &a3, &b) )
1031  {
1032  /* Yes. Mark this edge and the two preceding it as arc components */
1033  LWDEBUGF(4, "pt_continues_arc #%d", current_arc);
1034  found_arc = LW_TRUE;
1035  for ( k = j-1; k > j-4; k-- )
1036  edges_in_arcs[k] = current_arc;
1037  }
1038  else
1039  {
1040  /* No. So we're done with this candidate arc */
1041  LWDEBUG(4, "pt_continues_arc = false");
1042  current_arc++;
1043  break;
1044  }
1045 
1046  memcpy(&a1, &a2, sizeof(POINT4D));
1047  memcpy(&a2, &a3, sizeof(POINT4D));
1048  memcpy(&a3, &b, sizeof(POINT4D));
1049  }
1050  /* Jump past all the edges that were added to the arc */
1051  if ( found_arc )
1052  {
1053  /* Check if an arc was composed by enough edges to be
1054  * really considered an arc
1055  * See http://trac.osgeo.org/postgis/ticket/2420
1056  */
1057  arc_edges = j - 1 - i;
1058  LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
1059  if ( first.x == b.x && first.y == b.y ) {
1060  LWDEBUG(4, "arc is a circle");
1061  num_quadrants = 4;
1062  }
1063  else {
1064  lw_arc_center((POINT2D*)&first, (POINT2D*)&b, (POINT2D*)&a1, (POINT2D*)&center);
1065  angle = lw_arc_angle((POINT2D*)&first, (POINT2D*)&center, (POINT2D*)&b);
1066  int p2_side = lw_segment_side((POINT2D*)&first, (POINT2D*)&a1, (POINT2D*)&b);
1067  if ( p2_side >= 0 ) angle = -angle;
1068 
1069  if ( angle < 0 ) angle = 2 * M_PI + angle;
1070  num_quadrants = ( 4 * angle ) / ( 2 * M_PI );
1071  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);
1072  }
1073  /* a1 is first point, b is last point */
1074  if ( arc_edges < min_quad_edges * num_quadrants ) {
1075  LWDEBUGF(4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants);
1076  for ( k = j-1; k >= i; k-- )
1077  edges_in_arcs[k] = 0;
1078  }
1079 
1080  i = j-1;
1081  }
1082  else
1083  {
1084  /* Mark this edge as a linear edge */
1085  edges_in_arcs[i] = 0;
1086  i = i+1;
1087  }
1088  }
1089 
1090 #if POSTGIS_DEBUG_LEVEL > 3
1091  {
1092  char *edgestr = lwalloc(num_edges+1);
1093  for ( i = 0; i < num_edges; i++ )
1094  {
1095  if ( edges_in_arcs[i] )
1096  edgestr[i] = 48 + edges_in_arcs[i];
1097  else
1098  edgestr[i] = '.';
1099  }
1100  edgestr[num_edges] = 0;
1101  LWDEBUGF(3, "edge pattern %s", edgestr);
1102  lwfree(edgestr);
1103  }
1104 #endif
1105 
1106  start = 0;
1107  edge_type = edges_in_arcs[0];
1108  outcol = lwcollection_construct_empty(COMPOUNDTYPE, srid, ptarray_has_z(points), ptarray_has_m(points));
1109  for( i = 1; i < num_edges; i++ )
1110  {
1111  if( edge_type != edges_in_arcs[i] )
1112  {
1113  end = i - 1;
1114  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1115  start = i;
1116  edge_type = edges_in_arcs[i];
1117  }
1118  }
1119  lwfree(edges_in_arcs); /* not needed anymore */
1120 
1121  /* Roll out last item */
1122  end = num_edges - 1;
1123  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1124 
1125  /* Strip down to singleton if only one entry */
1126  if ( outcol->ngeoms == 1 )
1127  {
1128  LWGEOM *outgeom = outcol->geoms[0];
1129  outcol->ngeoms = 0; lwcollection_free(outcol);
1130  return outgeom;
1131  }
1132  return lwcollection_as_lwgeom(outcol);
1133 }
#define LW_FALSE
Definition: liblwgeom.h:108
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:291
#define COMPOUNDTYPE
Definition: liblwgeom.h:124
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:107
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:229
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:65
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:869
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:891
static LWGEOM * geom_from_pa(const POINTARRAY *pa, int32_t srid, int is_arc, int start, int end)
Definition: lwstroke.c:965
uint32_t ngeoms
Definition: liblwgeom.h:566
LWGEOM ** geoms
Definition: liblwgeom.h:561
double y
Definition: liblwgeom.h:376
double x
Definition: liblwgeom.h:376
double x
Definition: liblwgeom.h:400
double y
Definition: liblwgeom.h:400
uint32_t npoints
Definition: liblwgeom.h:413

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: