PostGIS  3.1.6dev-r@@SVN_REVISION@@

◆ lwpoly_to_points()

LWMPOINT* lwpoly_to_points ( const LWPOLY poly,
uint32_t  npoints,
int32_t  seed 
)

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

1629 {
1630  double area, bbox_area, bbox_width, bbox_height;
1631  GBOX bbox;
1632  const LWGEOM* lwgeom = (LWGEOM*)lwpoly;
1633  uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
1634  double sample_cell_size;
1635  uint32_t i, j, n;
1636  uint32_t iterations = 0;
1637  uint32_t npoints_generated = 0;
1638  uint32_t npoints_tested = 0;
1639  GEOSGeometry* g;
1640  const GEOSPreparedGeometry* gprep;
1641  GEOSGeometry* gpt;
1642  GEOSCoordSequence* gseq;
1643  LWMPOINT* mpt;
1644  int32_t srid = lwgeom_get_srid(lwgeom);
1645  int done = 0;
1646  int* cells;
1647  const size_t size = 2 * sizeof(int);
1648  char tmp[2 * sizeof(int)];
1649  const size_t stride = 2 * sizeof(int);
1650 
1651  if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
1652  {
1653  lwerror("%s: only polygons supported", __func__);
1654  return NULL;
1655  }
1656 
1657  if (npoints == 0 || lwgeom_is_empty(lwgeom)) return NULL;
1658 
1659  if (!lwpoly->bbox)
1660  lwgeom_calculate_gbox(lwgeom, &bbox);
1661  else
1662  bbox = *(lwpoly->bbox);
1663 
1664  area = lwpoly_area(lwpoly);
1665  bbox_width = bbox.xmax - bbox.xmin;
1666  bbox_height = bbox.ymax - bbox.ymin;
1667  bbox_area = bbox_width * bbox_height;
1668 
1669  if (area == 0.0 || bbox_area == 0.0)
1670  {
1671  lwerror("%s: zero area input polygon, TBD", __func__);
1672  return NULL;
1673  }
1674 
1675  /* Gross up our test set a bit to increase odds of getting coverage in one pass */
1676  sample_npoints = npoints * bbox_area / area;
1677 
1678  /* We're going to generate points using a sample grid as described
1679  * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html to try and get a more uniform
1680  * "random" set of points. So we have to figure out how to stick a grid into our box */
1681  sample_sqrt = lround(sqrt(sample_npoints));
1682  if (sample_sqrt == 0)
1683  sample_sqrt = 1;
1684 
1685  /* Calculate the grids we're going to randomize within */
1686  if (bbox_width > bbox_height)
1687  {
1688  sample_width = sample_sqrt;
1689  sample_height = ceil((double)sample_npoints / (double)sample_width);
1690  sample_cell_size = bbox_width / sample_width;
1691  }
1692  else
1693  {
1694  sample_height = sample_sqrt;
1695  sample_width = ceil((double)sample_npoints / (double)sample_height);
1696  sample_cell_size = bbox_height / sample_height;
1697  }
1698 
1699  /* Prepare the polygon for fast true/false testing */
1700  initGEOS(lwnotice, lwgeom_geos_error);
1701  g = (GEOSGeometry*)LWGEOM2GEOS(lwgeom, 0);
1702  if (!g)
1703  {
1704  lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
1705  return NULL;
1706  }
1707  gprep = GEOSPrepare(g);
1708 
1709  /* Get an empty multi-point ready to return */
1710  mpt = lwmpoint_construct_empty(srid, 0, 0);
1711 
1712  /* Initiate random number generator.
1713  * Repeatable numbers are generated with seed values >= 1.
1714  * When seed is zero and has not previously been set, it is based on
1715  * Unix time (seconds) and process ID. */
1716  lwrandom_set_seed(seed);
1717 
1718  /* Now we fill in an array of cells, and then shuffle that array, */
1719  /* so we can visit the cells in random order to avoid visual ugliness */
1720  /* caused by visiting them sequentially */
1721  cells = lwalloc(2 * sizeof(int) * sample_height * sample_width);
1722  for (i = 0; i < sample_width; i++)
1723  {
1724  for (j = 0; j < sample_height; j++)
1725  {
1726  cells[2 * (i * sample_height + j)] = i;
1727  cells[2 * (i * sample_height + j) + 1] = j;
1728  }
1729  }
1730 
1731  /* Fisher-Yates shuffle */
1732  n = sample_height * sample_width;
1733  if (n > 1)
1734  {
1735  for (i = n - 1; i > 0; i--)
1736  {
1737  size_t j = (size_t)(lwrandom_uniform() * (i + 1));
1738 
1739  memcpy(tmp, (char *)cells + j * stride, size);
1740  memcpy((char *)cells + j * stride, (char *)cells + i * stride, size);
1741  memcpy((char *)cells + i * stride, tmp, size);
1742  }
1743  }
1744 
1745  /* Start testing points */
1746  while (npoints_generated < npoints)
1747  {
1748  iterations++;
1749  for (i = 0; i < sample_width * sample_height; i++)
1750  {
1751  int contains = 0;
1752  double y = bbox.ymin + cells[2 * i] * sample_cell_size;
1753  double x = bbox.xmin + cells[2 * i + 1] * sample_cell_size;
1754  x += lwrandom_uniform() * sample_cell_size;
1755  y += lwrandom_uniform() * sample_cell_size;
1756  if (x >= bbox.xmax || y >= bbox.ymax) continue;
1757 
1758  gseq = GEOSCoordSeq_create(1, 2);
1759 #if POSTGIS_GEOS_VERSION < 38
1760  GEOSCoordSeq_setX(gseq, 0, x);
1761  GEOSCoordSeq_setY(gseq, 0, y);
1762 #else
1763  GEOSCoordSeq_setXY(gseq, 0, x, y);
1764 #endif
1765  gpt = GEOSGeom_createPoint(gseq);
1766 
1767  contains = GEOSPreparedIntersects(gprep, gpt);
1768 
1769  GEOSGeom_destroy(gpt);
1770 
1771  if (contains == 2)
1772  {
1773  GEOSPreparedGeom_destroy(gprep);
1774  GEOSGeom_destroy(g);
1775  lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
1776  return NULL;
1777  }
1778  if (contains)
1779  {
1780  npoints_generated++;
1781  mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
1782  if (npoints_generated == npoints)
1783  {
1784  done = 1;
1785  break;
1786  }
1787  }
1788 
1789  /* Short-circuit check for ctrl-c occasionally */
1790  npoints_tested++;
1791  if (npoints_tested % 10000 == 0)
1792  LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
1793 
1794  if (done) break;
1795  }
1796  if (done || iterations > 100) break;
1797  }
1798 
1799  GEOSPreparedGeom_destroy(gprep);
1800  GEOSGeom_destroy(g);
1801  lwfree(cells);
1802 
1803  return mpt;
1804 }
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
void lwgeom_geos_error(const char *fmt,...)
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:910
LWPOINT * lwpoint_make2d(int32_t srid, double x, double y)
Definition: lwpoint.c:163
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
void lwfree(void *mem)
Definition: lwutil.c:242
#define POLYGONTYPE
Definition: liblwgeom.h:118
LWMPOINT * lwmpoint_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwmpoint.c:39
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:738
void * lwalloc(size_t size)
Definition: lwutil.c:227
#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:434
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
static uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwinline.h:145
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
void lwrandom_set_seed(int32_t seed)
Definition: lwrandom.c:48
double lwrandom_uniform(void)
Definition: lwrandom.c:94
Datum contains(PG_FUNCTION_ARGS)
double ymax
Definition: liblwgeom.h:371
double xmax
Definition: liblwgeom.h:369
double ymin
Definition: liblwgeom.h:370
double xmin
Definition: liblwgeom.h:368

References 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(), lwrandom_set_seed(), lwrandom_uniform(), 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: