PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ ptarray_segmentize_sphere()

static POINTARRAY* ptarray_segmentize_sphere ( const POINTARRAY pa_in,
double  max_seg_length 

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

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

Definition at line 1597 of file lwgeodetic.c.

References geog2cart(), geographic_point_init(), getPoint4d_p(), LW_TRUE, lwerror(), POINTARRAY::npoints, p4d_same(), ptarray_append_point(), ptarray_construct_empty(), ptarray_has_m(), ptarray_has_z(), ptarray_segmentize_sphere_edge_recursive(), sphere_distance(), POINT4D::x, and POINT4D::y.

Referenced by lwgeom_segmentize_sphere().

1598 {
1599  POINTARRAY *pa_out;
1600  int hasz = ptarray_has_z(pa_in);
1601  int hasm = ptarray_has_m(pa_in);
1602  POINT4D p1, p2;
1603  POINT3D q1, q2;
1604  GEOGRAPHIC_POINT g1, g2;
1605  uint32_t i;
1607  /* Just crap out on crazy input */
1608  if ( ! pa_in )
1609  lwerror("%s: null input pointarray", __func__);
1610  if ( max_seg_length <= 0.0 )
1611  lwerror("%s: maximum segment length must be positive", __func__);
1613  /* Empty starting array */
1614  pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints);
1616  /* Simple loop per edge */
1617  for (i = 1; i < pa_in->npoints; i++)
1618  {
1619  getPoint4d_p(pa_in, i-1, &p1);
1620  getPoint4d_p(pa_in, i, &p2);
1621  geographic_point_init(p1.x, p1.y, &g1);
1622  geographic_point_init(p2.x, p2.y, &g2);
1624  /* Skip duplicate points (except in case of 2-point lines!) */
1625  if ((pa_in->npoints > 2) && p4d_same(&p1, &p2))
1626  continue;
1628  /* How long is this edge? */
1629  double d = sphere_distance(&g1, &g2);
1631  if (d > max_seg_length)
1632  {
1633  geog2cart(&g1, &q1);
1634  geog2cart(&g2, &q2);
1635  /* 3-d end points, XYZM end point, current edge size, min edge size */
1636  ptarray_segmentize_sphere_edge_recursive(&q1, &q2, &p1, &p2, d, max_seg_length, pa_out);
1637  }
1638  /* If we don't segmentize, we need to add first point manually */
1639  else
1640  {
1641  ptarray_append_point(pa_out, &p1, LW_TRUE);
1642  }
1643  }
1644  /* Always add the last point */
1645  ptarray_append_point(pa_out, &p2, LW_TRUE);
1646  return pa_out;
1647 }
double x
Definition: liblwgeom.h:354
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:917
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
unsigned int uint32_t
Definition: uthash.h:78
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_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:373
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:171
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
double y
Definition: liblwgeom.h:354
static int ptarray_segmentize_sphere_edge_recursive(const POINT3D *p1, const POINT3D *p2, const POINT4D *v1, const POINT4D *v2, double d, double max_seg_length, POINTARRAY *pa)
Definition: lwgeodetic.c:1546
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
int p4d_same(const POINT4D *p1, const POINT4D *p2)
Definition: lwalgorithm.c:31
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: