PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ lwpoly_to_points()

LWMPOINT* lwpoly_to_points ( const LWPOLY lwpoly,
uint32_t  npoints 
)

Definition at line 1486 of file liblwgeom/lwgeom_geos.c.

1487 {
1488  double area, bbox_area, bbox_width, bbox_height;
1489  GBOX bbox;
1490  const LWGEOM* lwgeom = (LWGEOM*)lwpoly;
1491  uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
1492  double sample_cell_size;
1493  uint32_t i, j, n;
1494  uint32_t iterations = 0;
1495  uint32_t npoints_generated = 0;
1496  uint32_t npoints_tested = 0;
1497  GEOSGeometry* g;
1498  const GEOSPreparedGeometry* gprep;
1499  GEOSGeometry* gpt;
1500  GEOSCoordSequence* gseq;
1501  LWMPOINT* mpt;
1502  int srid = lwgeom_get_srid(lwgeom);
1503  int done = 0;
1504  int* cells;
1505  const size_t size = 2 * sizeof(int);
1506  char tmp[2 * sizeof(int)];
1507  const size_t stride = 2 * sizeof(int);
1508 
1509  if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
1510  {
1511  lwerror("%s: only polygons supported", __func__);
1512  return NULL;
1513  }
1514 
1515  if (npoints == 0 || lwgeom_is_empty(lwgeom)) return NULL;
1516 
1517  if (!lwpoly->bbox)
1518  lwgeom_calculate_gbox(lwgeom, &bbox);
1519  else
1520  bbox = *(lwpoly->bbox);
1521 
1522  area = lwpoly_area(lwpoly);
1523  bbox_width = bbox.xmax - bbox.xmin;
1524  bbox_height = bbox.ymax - bbox.ymin;
1525  bbox_area = bbox_width * bbox_height;
1526 
1527  if (area == 0.0 || bbox_area == 0.0)
1528  {
1529  lwerror("%s: zero area input polygon, TBD", __func__);
1530  return NULL;
1531  }
1532 
1533  /* Gross up our test set a bit to increase odds of getting coverage in one pass */
1534  sample_npoints = npoints * bbox_area / area;
1535 
1536  /* We're going to generate points using a sample grid as described
1537  * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html to try and get a more uniform
1538  * "random" set of points. So we have to figure out how to stick a grid into our box */
1539  sample_sqrt = lround(sqrt(sample_npoints));
1540  if (sample_sqrt == 0) sample_sqrt = 1;
1541 
1542  /* Calculate the grids we're going to randomize within */
1543  if (bbox_width > bbox_height)
1544  {
1545  sample_width = sample_sqrt;
1546  sample_height = ceil((double)sample_npoints / (double)sample_width);
1547  sample_cell_size = bbox_width / sample_width;
1548  }
1549  else
1550  {
1551  sample_height = sample_sqrt;
1552  sample_width = ceil((double)sample_npoints / (double)sample_height);
1553  sample_cell_size = bbox_height / sample_height;
1554  }
1555 
1556  /* Prepare the polygon for fast true/false testing */
1557  initGEOS(lwnotice, lwgeom_geos_error);
1558  g = (GEOSGeometry*)LWGEOM2GEOS(lwgeom, 0);
1559  if (!g)
1560  {
1561  lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
1562  return NULL;
1563  }
1564  gprep = GEOSPrepare(g);
1565 
1566  /* Get an empty multi-point ready to return */
1567  mpt = lwmpoint_construct_empty(srid, 0, 0);
1568 
1569  /* Init random number generator */
1570  srand(time(NULL));
1571 
1572  /* Now we fill in an array of cells, and then shuffle that array, */
1573  /* so we can visit the cells in random order to avoid visual ugliness */
1574  /* caused by visiting them sequentially */
1575  cells = lwalloc(2 * sizeof(int) * sample_height * sample_width);
1576  for (i = 0; i < sample_width; i++)
1577  {
1578  for (j = 0; j < sample_height; j++)
1579  {
1580  cells[2 * (i * sample_height + j)] = i;
1581  cells[2 * (i * sample_height + j) + 1] = j;
1582  }
1583  }
1584 
1585  /* shuffle */
1586  n = sample_height * sample_width;
1587  if (n > 1)
1588  {
1589  for (i = 0; i < n - 1; ++i)
1590  {
1591  size_t rnd = (size_t)rand();
1592  size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
1593 
1594  memcpy(tmp, (char *)cells + j * stride, size);
1595  memcpy((char *)cells + j * stride, (char *)cells + i * stride, size);
1596  memcpy((char *)cells + i * stride, tmp, size);
1597  }
1598  }
1599 
1600  /* Start testing points */
1601  while (npoints_generated < npoints)
1602  {
1603  iterations++;
1604  for (i = 0; i < sample_width * sample_height; i++)
1605  {
1606  int contains = 0;
1607  double y = bbox.ymin + cells[2 * i] * sample_cell_size;
1608  double x = bbox.xmin + cells[2 * i + 1] * sample_cell_size;
1609  x += rand() * sample_cell_size / RAND_MAX;
1610  y += rand() * sample_cell_size / RAND_MAX;
1611  if (x >= bbox.xmax || y >= bbox.ymax) continue;
1612 
1613  gseq = GEOSCoordSeq_create(1, 2);
1614 #if POSTGIS_GEOS_VERSION < 38
1615  GEOSCoordSeq_setX(gseq, 0, x);
1616  GEOSCoordSeq_setY(gseq, 0, y);
1617 #else
1618  GEOSCoordSeq_setXY(gseq, 0, x, y);
1619 #endif
1620  gpt = GEOSGeom_createPoint(gseq);
1621 
1622  contains = GEOSPreparedIntersects(gprep, gpt);
1623 
1624  GEOSGeom_destroy(gpt);
1625 
1626  if (contains == 2)
1627  {
1628  GEOSPreparedGeom_destroy(gprep);
1629  GEOSGeom_destroy(g);
1630  lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
1631  return NULL;
1632  }
1633  if (contains)
1634  {
1635  npoints_generated++;
1636  mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
1637  if (npoints_generated == npoints)
1638  {
1639  done = 1;
1640  break;
1641  }
1642  }
1643 
1644  /* Short-circuit check for ctrl-c occasionally */
1645  npoints_tested++;
1646  if (npoints_tested % 10000 == 0)
1647  LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
1648 
1649  if (done) break;
1650  }
1651  if (done || iterations > 100) break;
1652  }
1653 
1654  GEOSPreparedGeom_destroy(gprep);
1655  GEOSGeom_destroy(g);
1656  lwfree(cells);
1657 
1658  return mpt;
1659 }
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
void lwgeom_geos_error(const char *fmt,...)
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:916
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:923
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
void lwfree(void *mem)
Definition: lwutil.c:244
#define POLYGONTYPE
Definition: liblwgeom.h:87
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwgeom.c:1393
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:746
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:39
void * lwalloc(size_t size)
Definition: lwutil.c:229
#define LW_ON_INTERRUPT(x)
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
Definition: lwpoly.c:441
Datum area(PG_FUNCTION_ARGS)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
Datum contains(PG_FUNCTION_ARGS)
double ymax
Definition: liblwgeom.h:298
double xmax
Definition: liblwgeom.h:296
double ymin
Definition: liblwgeom.h:297
double xmin
Definition: liblwgeom.h:295
GBOX * bbox
Definition: liblwgeom.h:456
unsigned int uint32_t
Definition: uthash.h:78

References area(), LWPOLY::bbox, contains(), LW_ON_INTERRUPT, lwalloc(), lwerror(), lwfree(), LWGEOM2GEOS(), lwgeom_calculate_gbox(), lwgeom_geos_errmsg, lwgeom_geos_error(), lwgeom_get_srid(), lwgeom_get_type(), lwgeom_is_empty(), lwmpoint_add_lwpoint(), lwmpoint_construct_empty(), lwnotice(), lwpoint_make2d(), lwpoly_area(), POLYGONTYPE, pixval::x, GBOX::xmax, GBOX::xmin, pixval::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_to_points(), and lwmpoly_to_points().

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