PostGIS  2.2.7dev-r@@SVN_REVISION@@
static POINTARRAY* ptarray_segmentize_sphere ( const POINTARRAY pa_in,
double  max_seg_length 
)
static

Create a new point array with no segment longer than the input segment length (expressed in radians!)

Parameters
pa_in- input point array pointer
max_seg_length- maximum output segment length in radians

Definition at line 1531 of file lwgeodetic.c.

References cart2geog(), geog2cart(), geographic_point_init(), getPoint4d_p(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_FALSE, LW_TRUE, lwerror(), POINT4D::m, normalize(), POINTARRAY::npoints, p4d_same(), ptarray_append_point(), ptarray_construct_empty(), ptarray_has_m(), ptarray_has_z(), rad2deg, sphere_distance(), POINT3D::x, POINT4D::x, POINT3D::y, POINT4D::y, POINT3D::z, and POINT4D::z.

Referenced by lwgeom_segmentize_sphere().

1532 {
1533  POINTARRAY *pa_out;
1534  int hasz = ptarray_has_z(pa_in);
1535  int hasm = ptarray_has_m(pa_in);
1536  int pa_in_offset = 0; /* input point offset */
1537  POINT4D p1, p2, p;
1538  POINT3D q1, q2, q, qn;
1539  GEOGRAPHIC_POINT g1, g2, g;
1540  double d;
1541 
1542  /* Just crap out on crazy input */
1543  if ( ! pa_in )
1544  lwerror("ptarray_segmentize_sphere: null input pointarray");
1545  if ( max_seg_length <= 0.0 )
1546  lwerror("ptarray_segmentize_sphere: maximum segment length must be positive");
1547 
1548  /* Empty starting array */
1549  pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints);
1550 
1551  /* Add first point */
1552  getPoint4d_p(pa_in, pa_in_offset, &p1);
1553  ptarray_append_point(pa_out, &p1, LW_FALSE);
1554  geographic_point_init(p1.x, p1.y, &g1);
1555  pa_in_offset++;
1556 
1557  while ( pa_in_offset < pa_in->npoints )
1558  {
1559  getPoint4d_p(pa_in, pa_in_offset, &p2);
1560  geographic_point_init(p2.x, p2.y, &g2);
1561 
1562  /* Skip duplicate points (except in case of 2-point lines!) */
1563  if ( (pa_in->npoints > 2) && p4d_same(&p1, &p2) )
1564  {
1565  /* Move one offset forward */
1566  p1 = p2;
1567  g1 = g2;
1568  pa_in_offset++;
1569  continue;
1570  }
1571 
1572  /* How long is this edge? */
1573  d = sphere_distance(&g1, &g2);
1574 
1575  /* We need to segmentize this edge */
1576  if ( d > max_seg_length )
1577  {
1578  int nsegs = 1 + d / max_seg_length;
1579  int i;
1580  double dx, dy, dz, dzz = 0, dmm = 0;
1581 
1582  geog2cart(&g1, &q1);
1583  geog2cart(&g2, &q2);
1584 
1585  dx = (q2.x - q1.x) / nsegs;
1586  dy = (q2.y - q1.y) / nsegs;
1587  dz = (q2.z - q1.z) / nsegs;
1588 
1589  /* The independent Z/M values on the ptarray */
1590  if ( hasz ) dzz = (p2.z - p1.z) / nsegs;
1591  if ( hasm ) dmm = (p2.m - p1.m) / nsegs;
1592 
1593  q = q1;
1594  p = p1;
1595 
1596  for ( i = 0; i < nsegs - 1; i++ )
1597  {
1598  /* Move one increment forwards */
1599  q.x += dx; q.y += dy; q.z += dz;
1600  qn = q;
1601  normalize(&qn);
1602 
1603  /* Back to spherical coordinates */
1604  cart2geog(&qn, &g);
1605  /* Back to lon/lat */
1606  p.x = rad2deg(g.lon);
1607  p.y = rad2deg(g.lat);
1608  if ( hasz )
1609  p.z += dzz;
1610  if ( hasm )
1611  p.m += dmm;
1612  ptarray_append_point(pa_out, &p, LW_FALSE);
1613  }
1614 
1615  ptarray_append_point(pa_out, &p2, LW_FALSE);
1616  }
1617  /* This edge is already short enough */
1618  else
1619  {
1620  ptarray_append_point(pa_out, &p2, (pa_in->npoints==2)?LW_TRUE:LW_FALSE);
1621  }
1622 
1623  /* Move one offset forward */
1624  p1 = p2;
1625  g1 = g2;
1626  pa_in_offset++;
1627  }
1628 
1629  return pa_out;
1630 }
double x
Definition: liblwgeom.h:336
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:898
double m
Definition: liblwgeom.h:336
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:565
double y
Definition: liblwgeom.h:324
int npoints
Definition: liblwgeom.h:355
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesion coordinates on unit sphere to spherical coordinates.
Definition: lwgeodetic.c:364
double x
Definition: liblwgeom.h:324
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:32
double z
Definition: liblwgeom.h:324
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
#define rad2deg(r)
Definition: lwgeodetic.h:60
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
double z
Definition: liblwgeom.h:336
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:354
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:156
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
double y
Definition: liblwgeom.h:336
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
int p4d_same(const POINT4D *p1, const POINT4D *p2)
Definition: lwalgorithm.c:17

Here is the call graph for this function:

Here is the caller graph for this function: