PostGIS  3.1.6dev-r@@SVN_REVISION@@

◆ lwgeom_subdivide_recursive()

static void lwgeom_subdivide_recursive ( const LWGEOM geom,
uint8_t  dimension,
uint32_t  maxvertices,
uint32_t  depth,
LWCOLLECTION col,
double  gridSize 
)
static

Definition at line 2261 of file lwgeom.c.

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

References FP_NEQUALS, FP_TOLERANCE, gbox_duplicate(), LWCOLLECTION::geoms, getPoint2d_cp(), LW_ON_INTERRUPT, LW_TRUE, lwcollection_add_lwgeom(), lwerror(), lwgeom_clone_deep(), lwgeom_count_vertices(), lwgeom_dimension(), lwgeom_free(), lwgeom_get_bbox(), lwgeom_intersection_prec(), 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_prec().

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