PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ lwpoly_to_points()

LWMPOINT* lwpoly_to_points ( const LWPOLY lwpoly,
int  npoints 
)

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

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().

1681 {
1682  double area, bbox_area, bbox_width, bbox_height;
1683  GBOX bbox;
1684  const LWGEOM *lwgeom = (LWGEOM*)lwpoly;
1685  int sample_npoints, sample_sqrt, sample_width, sample_height;
1686  double sample_cell_size;
1687  int i, j, n;
1688  int iterations = 0;
1689  int npoints_generated = 0;
1690  int npoints_tested = 0;
1691  GEOSGeometry *g;
1692  const GEOSPreparedGeometry *gprep;
1693  GEOSGeometry *gpt;
1694  GEOSCoordSequence *gseq;
1695  LWMPOINT *mpt;
1696  int srid = lwgeom_get_srid(lwgeom);
1697  int done = 0;
1698  int *cells;
1699  const size_t size = 2*sizeof(int);
1700  char tmp[2*sizeof(int)];
1701  const size_t stride = 2*sizeof(int);
1702 
1703 
1704  if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
1705  {
1706  lwerror("%s: only polygons supported", __func__);
1707  return NULL;
1708  }
1709 
1710  if (npoints == 0 || lwgeom_is_empty(lwgeom))
1711  {
1712  return NULL;
1713  // return lwmpoint_construct_empty(lwgeom_get_srid(poly), lwgeom_has_z(poly), lwgeom_has_m(poly));
1714  }
1715 
1716  if (!lwpoly->bbox)
1717  {
1718  lwgeom_calculate_gbox(lwgeom, &bbox);
1719  }
1720  else
1721  {
1722  bbox = *(lwpoly->bbox);
1723  }
1724  area = lwpoly_area(lwpoly);
1725  bbox_width = bbox.xmax - bbox.xmin;
1726  bbox_height = bbox.ymax - bbox.ymin;
1727  bbox_area = bbox_width * bbox_height;
1728 
1729  if (area == 0.0 || bbox_area == 0.0)
1730  {
1731  lwerror("%s: zero area input polygon, TBD", __func__);
1732  return NULL;
1733  }
1734 
1735  /* Gross up our test set a bit to increase odds of getting */
1736  /* coverage in one pass */
1737  sample_npoints = npoints * bbox_area / area;
1738 
1739  /* We're going to generate points using a sample grid */
1740  /* as described http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html */
1741  /* to try and get a more uniform "random" set of points */
1742  /* So we have to figure out how to stick a grid into our box */
1743  sample_sqrt = lround(sqrt(sample_npoints));
1744  if (sample_sqrt == 0)
1745  sample_sqrt = 1;
1746 
1747  /* Calculate the grids we're going to randomize within */
1748  if (bbox_width > bbox_height)
1749  {
1750  sample_width = sample_sqrt;
1751  sample_height = ceil((double)sample_npoints / (double)sample_width);
1752  sample_cell_size = bbox_width / sample_width;
1753  }
1754  else
1755  {
1756  sample_height = sample_sqrt;
1757  sample_width = ceil((double)sample_npoints / (double)sample_height);
1758  sample_cell_size = bbox_height / sample_height;
1759  }
1760 
1761  /* Prepare the polygon for fast true/false testing */
1762  initGEOS(lwnotice, lwgeom_geos_error);
1763  g = (GEOSGeometry *)LWGEOM2GEOS(lwgeom, 0);
1764  if (!g)
1765  {
1766  lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
1767  return NULL;
1768  }
1769  gprep = GEOSPrepare(g);
1770 
1771  /* Get an empty multi-point ready to return */
1772  mpt = lwmpoint_construct_empty(srid, 0, 0);
1773 
1774  /* Init random number generator */
1775  srand(time(NULL));
1776 
1777  /* Now we fill in an array of cells, and then shuffle that array, */
1778  /* so we can visit the cells in random order to avoid visual ugliness */
1779  /* caused by visiting them sequentially */
1780  cells = lwalloc(2*sizeof(int)*sample_height*sample_width);
1781  for (i = 0; i < sample_width; i++)
1782  {
1783  for (j = 0; j < sample_height; j++)
1784  {
1785  cells[2*(i*sample_height+j)] = i;
1786  cells[2*(i*sample_height+j)+1] = j;
1787  }
1788  }
1789 
1790  /* shuffle */
1791  {
1792  n = sample_height*sample_width;
1793  if (n > 1) {
1794  for (i = 0; i < n - 1; ++i) {
1795  size_t rnd = (size_t) rand();
1796  size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
1797 
1798  memcpy(tmp, (char *)cells + j * stride, size);
1799  memcpy((char *)cells + j * stride, (char *)cells + i * stride, size);
1800  memcpy((char *)cells + i * stride, tmp, size);
1801  }
1802  }
1803  }
1804 
1805 
1806  /* Start testing points */
1807  while (npoints_generated < npoints)
1808  {
1809  iterations++;
1810  for (i = 0; i < sample_width*sample_height; i++)
1811  {
1812  int contains = 0;
1813  double y = bbox.ymin + cells[2*i] * sample_cell_size;
1814  double x = bbox.xmin + cells[2*i+1] * sample_cell_size;
1815  x += rand() * sample_cell_size / RAND_MAX;
1816  y += rand() * sample_cell_size / RAND_MAX;
1817  if (x >= bbox.xmax || y >= bbox.ymax)
1818  continue;
1819 
1820  gseq = GEOSCoordSeq_create(1, 2);
1821  GEOSCoordSeq_setX(gseq, 0, x);
1822  GEOSCoordSeq_setY(gseq, 0, y);
1823  gpt = GEOSGeom_createPoint(gseq);
1824 
1825  contains = GEOSPreparedIntersects(gprep, gpt);
1826 
1827  GEOSGeom_destroy(gpt);
1828 
1829  if (contains == 2)
1830  {
1831  GEOSPreparedGeom_destroy(gprep);
1832  GEOSGeom_destroy(g);
1833  lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
1834  return NULL;
1835  }
1836  if (contains)
1837  {
1838  npoints_generated++;
1839  mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
1840  if (npoints_generated == npoints)
1841  {
1842  done = 1;
1843  break;
1844  }
1845  }
1846 
1847  /* Short-circuit check for ctrl-c occasionally */
1848  npoints_tested++;
1849  if (npoints_tested % 10000 == 0)
1850  {
1851  LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
1852  }
1853 
1854  if (done) break;
1855  }
1856  if (done || iterations > 100) break;
1857  }
1858 
1859  GEOSPreparedGeom_destroy(gprep);
1860  GEOSGeom_destroy(g);
1861  lwfree(cells);
1862 
1863  return mpt;
1864 }
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:878
#define POLYGONTYPE
Definition: liblwgeom.h:87
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
Datum area(PG_FUNCTION_ARGS)
double xmax
Definition: liblwgeom.h:293
Datum contains(PG_FUNCTION_ARGS)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:871
#define LW_ON_INTERRUPT(x)
GBOX * bbox
Definition: liblwgeom.h:453
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:701
double ymin
Definition: liblwgeom.h:294
double xmin
Definition: liblwgeom.h:292
void lwgeom_geos_error(const char *fmt,...)
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
Definition: lwpoly.c:524
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:39
double ymax
Definition: liblwgeom.h:295
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
void * lwalloc(size_t size)
Definition: lwutil.c:229
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:1346
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
Here is the call graph for this function:
Here is the caller graph for this function: