1630 double area, bbox_area, bbox_width, bbox_height;
1633 uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
1634 double sample_cell_size;
1636 uint32_t iterations = 0;
1637 uint32_t npoints_generated = 0;
1638 uint32_t npoints_tested = 0;
1640 const GEOSPreparedGeometry* gprep;
1642 GEOSCoordSequence* gseq;
1647 const size_t size = 2 *
sizeof(int);
1648 char tmp[2 *
sizeof(int)];
1649 const size_t stride = 2 *
sizeof(int);
1653 lwerror(
"%s: only polygons supported", __func__);
1662 bbox = *(lwpoly->bbox);
1665 bbox_width = bbox.
xmax - bbox.
xmin;
1666 bbox_height = bbox.
ymax - bbox.
ymin;
1667 bbox_area = bbox_width * bbox_height;
1669 if (area == 0.0 || bbox_area == 0.0)
1671 lwerror(
"%s: zero area input polygon, TBD", __func__);
1676 sample_npoints = npoints * bbox_area / area;
1681 sample_sqrt = lround(sqrt(sample_npoints));
1682 if (sample_sqrt == 0)
1686 if (bbox_width > bbox_height)
1688 sample_width = sample_sqrt;
1689 sample_height = ceil((
double)sample_npoints / (
double)sample_width);
1690 sample_cell_size = bbox_width / sample_width;
1694 sample_height = sample_sqrt;
1695 sample_width = ceil((
double)sample_npoints / (
double)sample_height);
1696 sample_cell_size = bbox_height / sample_height;
1707 gprep = GEOSPrepare(g);
1721 cells =
lwalloc(2 *
sizeof(
int) * sample_height * sample_width);
1722 for (i = 0; i < sample_width; i++)
1724 for (j = 0; j < sample_height; j++)
1726 cells[2 * (i * sample_height + j)] = i;
1727 cells[2 * (i * sample_height + j) + 1] = j;
1732 n = sample_height * sample_width;
1735 for (i = n - 1; i > 0; i--)
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);
1746 while (npoints_generated < npoints)
1749 for (i = 0; i < sample_width * sample_height; i++)
1752 double y = bbox.
ymin + cells[2 * i] * sample_cell_size;
1753 double x = bbox.
xmin + cells[2 * i + 1] * sample_cell_size;
1756 if (
x >= bbox.
xmax ||
y >= bbox.
ymax)
continue;
1758 gseq = GEOSCoordSeq_create(1, 2);
1759 #if POSTGIS_GEOS_VERSION < 38
1760 GEOSCoordSeq_setX(gseq, 0,
x);
1761 GEOSCoordSeq_setY(gseq, 0,
y);
1763 GEOSCoordSeq_setXY(gseq, 0,
x,
y);
1765 gpt = GEOSGeom_createPoint(gseq);
1767 contains = GEOSPreparedIntersects(gprep, gpt);
1769 GEOSGeom_destroy(gpt);
1773 GEOSPreparedGeom_destroy(gprep);
1774 GEOSGeom_destroy(g);
1780 npoints_generated++;
1782 if (npoints_generated == npoints)
1791 if (npoints_tested % 10000 == 0)
1792 LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g);
return NULL);
1796 if (done || iterations > 100)
break;
1799 GEOSPreparedGeom_destroy(gprep);
1800 GEOSGeom_destroy(g);
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.
LWPOINT * lwpoint_make2d(int32_t srid, double x, double y)
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
LWMPOINT * lwmpoint_construct_empty(int32_t srid, char hasz, char hasm)
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...
void * lwalloc(size_t size)
#define LW_ON_INTERRUPT(x)
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
static uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
void lwrandom_set_seed(int32_t seed)
double lwrandom_uniform(void)
Datum contains(PG_FUNCTION_ARGS)