PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ lwgeom_subdivide_recursive()

static int lwgeom_subdivide_recursive ( const LWGEOM geom,
uint8_t  dimension,
uint32_t  maxvertices,
uint32_t  depth,
LWCOLLECTION col 
)
static

Definition at line 2265 of file lwgeom.c.

2266 {
2267  const uint32_t maxdepth = 50;
2268  GBOX clip, subbox1, subbox2;
2269  uint32_t nvertices = 0;
2270  uint32_t i, n = 0;
2271  uint32_t split_ordinate;
2272  double width;
2273  double height;
2274  double pivot = DBL_MAX;
2275  double center = DBL_MAX;
2276  LWPOLY *lwpoly = NULL;
2277  LWGEOM *clipped;
2278  const GBOX *box_in;
2279  if (!geom) return 0;
2280  box_in = lwgeom_get_bbox(geom);
2281  if (!box_in) return 0;
2282  gbox_duplicate(box_in, &clip);
2283  width = clip.xmax - clip.xmin;
2284  height = clip.ymax - clip.ymin;
2285 
2286  if ( geom->type == POLYHEDRALSURFACETYPE || geom->type == TINTYPE )
2287  lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(geom->type));
2288 
2289  if ( width == 0.0 && height == 0.0 )
2290  {
2291  if ( geom->type == POINTTYPE && dimension == 0)
2292  {
2294  return 1;
2295  }
2296  else
2297  return 0;
2298  }
2299 
2300  if (width == 0.0)
2301  {
2302  clip.xmax += FP_TOLERANCE;
2303  clip.xmin -= FP_TOLERANCE;
2304  width = 2 * FP_TOLERANCE;
2305  }
2306  if (height == 0.0)
2307  {
2308  clip.ymax += FP_TOLERANCE;
2309  clip.ymin -= FP_TOLERANCE;
2310  height = 2 * FP_TOLERANCE;
2311  }
2312 
2313  /* Always just recurse into collections */
2314  if ( lwgeom_is_collection(geom) && geom->type != MULTIPOINTTYPE )
2315  {
2316  LWCOLLECTION *incol = (LWCOLLECTION*)geom;
2317  int n = 0;
2318  /* Don't increment depth yet, since we aren't actually
2319  * subdividing geometries yet */
2320  for ( i = 0; i < incol->ngeoms; i++ )
2321  n += lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col);
2322  return n;
2323  }
2324 
2325  if (lwgeom_dimension(geom) < dimension)
2326  {
2327  /* We've hit a lower dimension object produced by clipping at
2328  * a shallower recursion level. Ignore it. */
2329  return 0;
2330  }
2331 
2332  /* But don't go too far. 2^50 ~= 10^15, that's enough subdivision */
2333  /* Just add what's left */
2334  if ( depth > maxdepth )
2335  {
2337  return 1;
2338  }
2339 
2340  nvertices = lwgeom_count_vertices(geom);
2341 
2342  /* Skip empties entirely */
2343  if (nvertices == 0) return 0;
2344 
2345  /* If it is under the vertex tolerance, just add it, we're done */
2346  if (nvertices <= maxvertices)
2347  {
2349  return 1;
2350  }
2351 
2352  split_ordinate = (width > height) ? 0 : 1;
2353  if (split_ordinate == 0)
2354  center = (clip.xmin + clip.xmax) / 2;
2355  else
2356  center = (clip.ymin + clip.ymax) / 2;
2357 
2358  if (geom->type == POLYGONTYPE)
2359  {
2360  uint32_t ring_to_trim = 0;
2361  double ring_area = 0;
2362  double pivot_eps = DBL_MAX;
2363  double pt_eps = DBL_MAX;
2364  POINTARRAY *pa;
2365  pivot = DBL_MAX;
2366  lwpoly = (LWPOLY *)geom;
2367 
2368  /* if there are more points in holes than in outer ring */
2369  if (nvertices >= 2 * lwpoly->rings[0]->npoints)
2370  {
2371  /* trim holes starting from biggest */
2372  for (i = 1; i < lwpoly->nrings; i++)
2373  {
2374  double current_ring_area = fabs(ptarray_signed_area(lwpoly->rings[i]));
2375  if (current_ring_area >= ring_area)
2376  {
2377  ring_area = current_ring_area;
2378  ring_to_trim = i;
2379  }
2380  }
2381  }
2382 
2383  pa = lwpoly->rings[ring_to_trim];
2384 
2385  /* find most central point in chosen ring */
2386  for (i = 0; i < pa->npoints; i++)
2387  {
2388  double pt;
2389  if (split_ordinate == 0)
2390  pt = getPoint2d_cp(pa, i)->x;
2391  else
2392  pt = getPoint2d_cp(pa, i)->y;
2393  pt_eps = fabs(pt - center);
2394  if (pivot_eps > pt_eps)
2395  {
2396  pivot = pt;
2397  pivot_eps = pt_eps;
2398  }
2399  }
2400  }
2401 
2402  gbox_duplicate(&clip, &subbox1);
2403  gbox_duplicate(&clip, &subbox2);
2404 
2405  if (pivot == DBL_MAX) pivot = center;
2406 
2407  if (split_ordinate == 0)
2408  {
2409  if (FP_NEQUALS(subbox1.xmax, pivot) && FP_NEQUALS(subbox1.xmin, pivot))
2410  subbox1.xmax = subbox2.xmin = pivot;
2411  else
2412  subbox1.xmax = subbox2.xmin = center;
2413  }
2414  else
2415  {
2416  if (FP_NEQUALS(subbox1.ymax, pivot) && FP_NEQUALS(subbox1.ymin, pivot))
2417  subbox1.ymax = subbox2.ymin = pivot;
2418  else
2419  subbox1.ymax = subbox2.ymin = center;
2420  }
2421 
2422  ++depth;
2423 
2424  LWGEOM* subbox = (LWGEOM*) lwpoly_construct_envelope(geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
2425  clipped = lwgeom_intersection(geom, subbox);
2426  lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
2427  lwgeom_free(subbox);
2428  if (clipped && !lwgeom_is_empty(clipped))
2429  {
2430  n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
2431  lwgeom_free(clipped);
2432  }
2433 
2434  subbox = (LWGEOM*) lwpoly_construct_envelope(geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
2435  clipped = lwgeom_intersection(geom, subbox);
2436  lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
2437  lwgeom_free(subbox);
2438  if (clipped && !lwgeom_is_empty(clipped))
2439  {
2440  n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
2441  lwgeom_free(clipped);
2442  }
2443 
2444  return n;
2445 }
void gbox_duplicate(const GBOX *original, GBOX *duplicate)
Copy the values of original GBOX into duplicate.
Definition: g_box.c:440
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
#define MULTIPOINTTYPE
Definition: liblwgeom.h:88
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:85
#define TINTYPE
Definition: liblwgeom.h:99
#define POLYGONTYPE
Definition: liblwgeom.h:87
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:97
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
LWPOLY * lwpoly_construct_envelope(int srid, double x1, double y1, double x2, double y2)
Definition: lwpoly.c:98
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from.
Definition: lwgeom_api.c:374
double ptarray_signed_area(const POINTARRAY *pa)
Returns the area in cartesian units.
Definition: ptarray.c:997
#define FP_TOLERANCE
Floating point comparators.
#define FP_NEQUALS(A, B)
void lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
Definition: lwgeom.c:1750
int lwgeom_is_collection(const LWGEOM *geom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep-clone an LWGEOM object.
Definition: lwgeom.c:520
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count points in an LWGEOM.
Definition: lwgeom.c:1235
int lwgeom_dimension(const LWGEOM *geom)
For an LWGEOM, returns 0 for points, 1 for lines, 2 for polygons, 3 for volume, and the max dimension...
Definition: lwgeom.c:1287
const GBOX * lwgeom_get_bbox(const LWGEOM *lwg)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:734
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:1393
void lwgeom_free(LWGEOM *lwgeom)
Definition: lwgeom.c:1144
static int lwgeom_subdivide_recursive(const LWGEOM *geom, uint8_t dimension, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col)
Definition: lwgeom.c:2265
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static double pivot(double *left, double *right)
Definition: rt_statistics.c:40
double ymax
Definition: liblwgeom.h:298
double xmax
Definition: liblwgeom.h:296
double ymin
Definition: liblwgeom.h:297
double xmin
Definition: liblwgeom.h:295
uint32_t ngeoms
Definition: liblwgeom.h:510
LWGEOM ** geoms
Definition: liblwgeom.h:512
uint8_t type
Definition: liblwgeom.h:399
int32_t srid
Definition: liblwgeom.h:402
POINTARRAY ** rings
Definition: liblwgeom.h:460
uint32_t nrings
Definition: liblwgeom.h:458
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
uint32_t npoints
Definition: liblwgeom.h:374
unsigned int uint32_t
Definition: uthash.h:78

References FP_NEQUALS, FP_TOLERANCE, gbox_duplicate(), LWCOLLECTION::geoms, getPoint2d_cp(), LW_TRUE, lwcollection_add_lwgeom(), lwerror(), lwgeom_clone_deep(), lwgeom_count_vertices(), lwgeom_dimension(), lwgeom_free(), lwgeom_get_bbox(), lwgeom_intersection(), lwgeom_is_collection(), lwgeom_is_empty(), lwgeom_simplify_in_place(), lwpoly_construct_envelope(), lwtype_name(), MULTIPOINTTYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, LWPOLY::nrings, pivot(), POINTTYPE, POLYGONTYPE, POLYHEDRALSURFACETYPE, ptarray_signed_area(), LWPOLY::rings, LWGEOM::srid, TINTYPE, LWGEOM::type, POINT2D::x, GBOX::xmax, GBOX::xmin, POINT2D::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_subdivide().

Here is the call graph for this function:
Here is the caller graph for this function: