PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ pta_unstroke()

LWGEOM * pta_unstroke ( const POINTARRAY points,
int  srid 
)

Definition at line 884 of file lwstroke.c.

885 {
886  int i = 0, j, k;
887  POINT4D a1, a2, a3, b;
888  POINT4D first, center;
889  char *edges_in_arcs;
890  int found_arc = LW_FALSE;
891  int current_arc = 1;
892  int num_edges;
893  int edge_type; /* non-zero if edge is part of an arc */
894  int start, end;
895  LWCOLLECTION *outcol;
896  /* Minimum number of edges, per quadrant, required to define an arc */
897  const unsigned int min_quad_edges = 2;
898 
899  /* Die on null input */
900  if ( ! points )
901  lwerror("pta_unstroke called with null pointarray");
902 
903  /* Null on empty input? */
904  if ( points->npoints == 0 )
905  return NULL;
906 
907  /* We can't desegmentize anything shorter than four points */
908  if ( points->npoints < 4 )
909  {
910  /* Return a linestring here*/
911  lwerror("pta_unstroke needs implementation for npoints < 4");
912  }
913 
914  /* Allocate our result array of vertices that are part of arcs */
915  num_edges = points->npoints - 1;
916  edges_in_arcs = lwalloc(num_edges + 1);
917  memset(edges_in_arcs, 0, num_edges + 1);
918 
919  /* We make a candidate arc of the first two edges, */
920  /* And then see if the next edge follows it */
921  while( i < num_edges-2 )
922  {
923  unsigned int arc_edges;
924  double num_quadrants;
925  double angle;
926 
927  found_arc = LW_FALSE;
928  /* Make candidate arc */
929  getPoint4d_p(points, i , &a1);
930  getPoint4d_p(points, i+1, &a2);
931  getPoint4d_p(points, i+2, &a3);
932  memcpy(&first, &a1, sizeof(POINT4D));
933 
934  for( j = i+3; j < num_edges+1; j++ )
935  {
936  LWDEBUGF(4, "i=%d, j=%d", i, j);
937  getPoint4d_p(points, j, &b);
938  /* Does this point fall on our candidate arc? */
939  if ( pt_continues_arc(&a1, &a2, &a3, &b) )
940  {
941  /* Yes. Mark this edge and the two preceding it as arc components */
942  LWDEBUGF(4, "pt_continues_arc #%d", current_arc);
943  found_arc = LW_TRUE;
944  for ( k = j-1; k > j-4; k-- )
945  edges_in_arcs[k] = current_arc;
946  }
947  else
948  {
949  /* No. So we're done with this candidate arc */
950  LWDEBUG(4, "pt_continues_arc = false");
951  current_arc++;
952  break;
953  }
954 
955  memcpy(&a1, &a2, sizeof(POINT4D));
956  memcpy(&a2, &a3, sizeof(POINT4D));
957  memcpy(&a3, &b, sizeof(POINT4D));
958  }
959  /* Jump past all the edges that were added to the arc */
960  if ( found_arc )
961  {
962  /* Check if an arc was composed by enough edges to be
963  * really considered an arc
964  * See http://trac.osgeo.org/postgis/ticket/2420
965  */
966  arc_edges = j - 1 - i;
967  LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
968  if ( first.x == b.x && first.y == b.y ) {
969  LWDEBUG(4, "arc is a circle");
970  num_quadrants = 4;
971  }
972  else {
973  lw_arc_center((POINT2D*)&first, (POINT2D*)&b, (POINT2D*)&a1, (POINT2D*)&center);
974  angle = lw_arc_angle((POINT2D*)&first, (POINT2D*)&center, (POINT2D*)&b);
975  int p2_side = lw_segment_side((POINT2D*)&first, (POINT2D*)&a1, (POINT2D*)&b);
976  if ( p2_side >= 0 ) angle = -angle;
977 
978  if ( angle < 0 ) angle = 2 * M_PI + angle;
979  num_quadrants = ( 4 * angle ) / ( 2 * M_PI );
980  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);
981  }
982  /* a1 is first point, b is last point */
983  if ( arc_edges < min_quad_edges * num_quadrants ) {
984  LWDEBUGF(4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants);
985  for ( k = j-1; k >= i; k-- )
986  edges_in_arcs[k] = 0;
987  }
988 
989  i = j-1;
990  }
991  else
992  {
993  /* Mark this edge as a linear edge */
994  edges_in_arcs[i] = 0;
995  i = i+1;
996  }
997  }
998 
999 #if POSTGIS_DEBUG_LEVEL > 3
1000  {
1001  char *edgestr = lwalloc(num_edges+1);
1002  for ( i = 0; i < num_edges; i++ )
1003  {
1004  if ( edges_in_arcs[i] )
1005  edgestr[i] = 48 + edges_in_arcs[i];
1006  else
1007  edgestr[i] = '.';
1008  }
1009  edgestr[num_edges] = 0;
1010  LWDEBUGF(3, "edge pattern %s", edgestr);
1011  lwfree(edgestr);
1012  }
1013 #endif
1014 
1015  start = 0;
1016  edge_type = edges_in_arcs[0];
1017  outcol = lwcollection_construct_empty(COMPOUNDTYPE, srid, ptarray_has_z(points), ptarray_has_m(points));
1018  for( i = 1; i < num_edges; i++ )
1019  {
1020  if( edge_type != edges_in_arcs[i] )
1021  {
1022  end = i - 1;
1023  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1024  start = i;
1025  edge_type = edges_in_arcs[i];
1026  }
1027  }
1028  lwfree(edges_in_arcs); /* not needed anymore */
1029 
1030  /* Roll out last item */
1031  end = num_edges - 1;
1032  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1033 
1034  /* Strip down to singleton if only one entry */
1035  if ( outcol->ngeoms == 1 )
1036  {
1037  LWGEOM *outgeom = outcol->geoms[0];
1038  outcol->ngeoms = 0; lwcollection_free(outcol);
1039  return outgeom;
1040  }
1041  return lwcollection_as_lwgeom(outcol);
1042 }
#define LW_FALSE
Definition: liblwgeom.h:77
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:300
#define COMPOUNDTYPE
Definition: liblwgeom.h:93
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
void lwfree(void *mem)
Definition: lwutil.c:244
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:356
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:123
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
void * lwalloc(size_t size)
Definition: lwutil.c:229
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
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:228
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
#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 LWGEOM * geom_from_pa(const POINTARRAY *pa, int srid, int is_arc, int start, int end)
Definition: lwstroke.c:874
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:778
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:800
uint32_t ngeoms
Definition: liblwgeom.h:510
LWGEOM ** geoms
Definition: liblwgeom.h:512
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:355
double y
Definition: liblwgeom.h:355
uint32_t npoints
Definition: liblwgeom.h:374

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: