PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ lwpoly_to_points()

LWMPOINT* lwpoly_to_points ( const LWPOLY poly,
int  npoints 
)

Definition at line 1668 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, shuffle(), pixval::x, GBOX::xmax, GBOX::xmin, pixval::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_to_points(), and lwmpoly_to_points().

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