lwpoly_to_points()

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

extern

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

{
 
        typedef struct CellCorner {
                double x;
                double y;
        } CellCorner;
 
        CellCorner* cells;
        uint32_t num_cells = 0;
 
        double area, bbox_area, bbox_width, bbox_height, area_ratio;
        GBOX bbox;
        const LWGEOM* lwgeom = (LWGEOM*)lwpoly;
        uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
        double sample_cell_size;
        uint32_t i, j, n;
        uint32_t iterations = 0;
        uint32_t npoints_generated = 0;
        uint32_t npoints_tested = 0;
        GEOSGeometry* g;
        const GEOSPreparedGeometry* gprep;
        LWMPOINT* mpt;
        int32_t srid = lwgeom_get_srid(lwgeom);
        int done = 0;
 
        if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
        {
                lwerror("%s: only polygons supported", __func__);
                return NULL;
        }
 
        if (npoints == 0 || lwgeom_is_empty(lwgeom)) return NULL;
 
        if (!lwpoly->bbox)
                lwgeom_calculate_gbox(lwgeom, &bbox);
        else
                bbox = *(lwpoly->bbox);
 
        area = lwpoly_area(lwpoly);
        bbox_width = bbox.xmax - bbox.xmin;
        bbox_height = bbox.ymax - bbox.ymin;
        bbox_area = bbox_width * bbox_height;
 
        if (area == 0.0 || bbox_area == 0.0)
        {
                lwerror("%s: zero area input polygon, TBD", __func__);
                return NULL;
        }
 
        /* Gross up our test set a bit to increase odds of getting coverage
         * in one pass. For narrow inputs, it is possible we will get
         * no hits at all. A fall-back approach might use a random stab-line
         * to find pairs of edges that bound interior area, and generate a point
         * between those edges, in the case that this gridding system
         * does not generate the desired number of points */
        area_ratio = FP_MIN(bbox_area / area, 10000.0);
        sample_npoints = npoints * area_ratio;
 
        /* We're going to generate points using a sample grid as described
         * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html
         * to try and get a more uniform "random" set of points. So we have to figure
         *  out how to stick a grid into our box */
        sample_sqrt = ceil(sqrt(sample_npoints));
 
        /* Calculate the grids we're going to randomize within */
        if (bbox_width > bbox_height)
        {
                sample_width = sample_sqrt;
                sample_height = ceil((double)sample_npoints / (double)sample_width);
                sample_cell_size = bbox_width / sample_width;
        }
        else
        {
                sample_height = sample_sqrt;
                sample_width = ceil((double)sample_npoints / (double)sample_height);
                sample_cell_size = bbox_height / sample_height;
        }
 
        /* Prepare the polygon for fast true/false testing */
        initGEOS(lwnotice, lwgeom_geos_error);
        g = (GEOSGeometry*)LWGEOM2GEOS(lwgeom, 0);
        if (!g)
        {
                lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
                return NULL;
        }
        gprep = GEOSPrepare(g);
 
        /* START_PROFILING(); */
 
        /* Now we fill in an array of cells, only using cells that actually
         * intersect the geometry of interest, and finally weshuffle that array,
         * so we can visit the cells in random order to avoid visual ugliness
         * caused by visiting them sequentially */
        cells = lwalloc(sizeof(CellCorner) * sample_height * sample_width);
 
        for (i = 0; i < sample_width; i++)
        {
                for (j = 0; j < sample_height; j++)
                {
                        int intersects = 0;
                        double xmin = bbox.xmin + i * sample_cell_size;
                        double ymin = bbox.ymin + j * sample_cell_size;
                        GEOSGeometry *gcell = lwpoly_to_points_make_cell(xmin, ymin, sample_cell_size);
                        intersects = GEOSPreparedIntersects(gprep, gcell);
                        GEOSGeom_destroy(gcell);
                        if (intersects == 2)
                        {
                                GEOSPreparedGeom_destroy(gprep);
                                GEOSGeom_destroy(g);
                                lwerror("%s: GEOS exception on GEOSPreparedIntersects: %s", __func__, lwgeom_geos_errmsg);
                                return NULL;
                        }
                        if (intersects == 1)
                        {
                                cells[num_cells].x = xmin;
                                cells[num_cells].y = ymin;
                                num_cells++;
                        }
                }
        }
 
        /* STOP_PROFILING(CellGeneration); */
 
        /* Initiate random number generator.
         * Repeatable numbers are generated with seed values >= 1.
         * When seed is zero and has not previously been set, it is based on
         * Unix time (seconds) and process ID. */
        lwrandom_set_seed(seed);
 
        /* Fisher-Yates shuffle */
        n = num_cells;
        if (n > 1)
        {
                for (i = n - 1; i > 0; i--)
                {
                        CellCorner tmp;
                        j = (uint32_t)(lwrandom_uniform() * (i + 1));
                        memcpy(     &tmp, cells + j, sizeof(CellCorner));
                        memcpy(cells + j, cells + i, sizeof(CellCorner));
                        memcpy(cells + i,      &tmp, sizeof(CellCorner));
                }
        }
 
        /* Get an empty array of points ready to return */
        mpt = lwmpoint_construct_empty(srid, 0, 0);
 
        /* Start testing points */
        while (npoints_generated < npoints)
        {
                iterations++;
                for (i = 0; i < num_cells; i++)
                {
                        int contains = 0;
                        double y = cells[i].y;
                        double x = cells[i].x;
                        x += lwrandom_uniform() * sample_cell_size;
                        y += lwrandom_uniform() * sample_cell_size;
                        if (x >= bbox.xmax || y >= bbox.ymax) continue;
 
#if POSTGIS_GEOS_VERSION >= 31200
                        contains = GEOSPreparedIntersectsXY(gprep, x, y);
#else
                        {
                                GEOSGeometry *gpt;
                                GEOSCoordSequence *gseq = GEOSCoordSeq_create(1, 2);;
                                GEOSCoordSeq_setXY(gseq, 0, x, y);
                                gpt = GEOSGeom_createPoint(gseq);
                                contains = GEOSPreparedIntersects(gprep, gpt);
                                GEOSGeom_destroy(gpt);
                        }
#endif
                        if (contains == 2)
                        {
                                GEOSPreparedGeom_destroy(gprep);
                                GEOSGeom_destroy(g);
                                lwerror("%s: GEOS exception on GEOSPreparedIntersects: %s", __func__, lwgeom_geos_errmsg);
                                return NULL;
                        }
                        if (contains == 1)
                        {
                                mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
                                npoints_generated++;
                                if (npoints_generated == npoints)
                                {
                                        done = 1;
                                        break;
                                }
                        }
 
                        /* Short-circuit check for ctrl-c occasionally */
                        npoints_tested++;
                        if (npoints_tested % 10000 == 0)
                                LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
 
                        if (done) break;
                }
                if (done || iterations > 100) break;
        }
 
        GEOSPreparedGeom_destroy(gprep);
        GEOSGeom_destroy(g);
        lwfree(cells);
 
        return mpt;
}

References LWPOLY::bbox, contains(), FP_MIN, 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(), lwpoly_to_points_make_cell(), lwrandom_set_seed(), lwrandom_uniform(), POLYGONTYPE, GBOX::xmax, GBOX::xmin, 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: