Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: CPL - Common Portability Library
4 : * Purpose: Implementation of quadtree building and searching functions.
5 : * Derived from shapelib and mapserver implementations
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : * Even Rouault, <even dot rouault at spatialys.com>
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 1999-2008, Frank Warmerdam
11 : * Copyright (c) 2008-2014, Even Rouault <even dot rouault at spatialys.com>
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ******************************************************************************
31 : */
32 :
33 : #include "cpl_port.h"
34 : #include "cpl_quad_tree.h"
35 :
36 : #include <algorithm>
37 : #include <cstdio>
38 : #include <cstring>
39 :
40 : #include "cpl_conv.h"
41 : #include "cpl_error.h"
42 :
43 : constexpr int MAX_DEFAULT_TREE_DEPTH = 12;
44 : constexpr int MAX_SUBNODES = 4;
45 :
46 : typedef struct _QuadTreeNode QuadTreeNode;
47 :
48 : struct _QuadTreeNode
49 : {
50 : /* area covered by this psNode */
51 : CPLRectObj rect;
52 :
53 : int nFeatures; /* number of shapes stored at this psNode. */
54 :
55 : int nNumSubNodes; /* number of active subnodes */
56 :
57 : void **pahFeatures; /* list of shapes stored at this psNode. */
58 : CPLRectObj *pasBounds;
59 :
60 : QuadTreeNode *apSubNode[MAX_SUBNODES];
61 : };
62 :
63 : struct _CPLQuadTree
64 : {
65 : QuadTreeNode *psRoot;
66 : CPLQuadTreeGetBoundsFunc pfnGetBounds;
67 : CPLQuadTreeGetBoundsExFunc pfnGetBoundsEx;
68 : void *pUserData;
69 : int nFeatures;
70 : int nMaxDepth;
71 : int nBucketCapacity;
72 : double dfSplitRatio;
73 : bool bForceUseOfSubNodes;
74 : };
75 :
76 : static void CPLQuadTreeAddFeatureInternal(CPLQuadTree *hQuadTree,
77 : void *hFeature,
78 : const CPLRectObj *pRect);
79 : static void CPLQuadTreeNodeDestroy(QuadTreeNode *psNode);
80 :
81 : /* -------------------------------------------------------------------- */
82 : /* If the following is 0.5, psNodes will be split in half. If it */
83 : /* is 0.6 then each apSubNode will contain 60% of the parent */
84 : /* psNode, with 20% representing overlap. This can be help to */
85 : /* prevent small objects on a boundary from shifting too high */
86 : /* up the hQuadTree. */
87 : /* -------------------------------------------------------------------- */
88 : constexpr double DEFAULT_SPLIT_RATIO = 0.55;
89 :
90 : /*
91 : ** Returns TRUE if rectangle a is contained in rectangle b
92 : */
93 2193100 : static CPL_INLINE bool CPL_RectContained(const CPLRectObj *a,
94 : const CPLRectObj *b)
95 : {
96 3424870 : return a->minx >= b->minx && a->maxx <= b->maxx && a->miny >= b->miny &&
97 3424870 : a->maxy <= b->maxy;
98 : }
99 :
100 : /*
101 : ** Returns TRUE if rectangles a and b overlap
102 : */
103 22282700 : static CPL_INLINE bool CPL_RectOverlap(const CPLRectObj *a, const CPLRectObj *b)
104 : {
105 22282700 : if (a->minx > b->maxx)
106 7156720 : return false;
107 15126000 : if (a->maxx < b->minx)
108 6500010 : return false;
109 8625980 : if (a->miny > b->maxy)
110 1877920 : return false;
111 6748060 : if (a->maxy < b->miny)
112 1840210 : return false;
113 4907850 : return true;
114 : }
115 :
116 : /************************************************************************/
117 : /* CPLQuadTreeNodeCreate() */
118 : /************************************************************************/
119 :
120 42691 : static QuadTreeNode *CPLQuadTreeNodeCreate(const CPLRectObj *pRect)
121 : {
122 : QuadTreeNode *psNode =
123 42691 : static_cast<QuadTreeNode *>(CPLMalloc(sizeof(QuadTreeNode)));
124 :
125 42691 : psNode->nFeatures = 0;
126 42691 : psNode->pahFeatures = nullptr;
127 42691 : psNode->pasBounds = nullptr;
128 :
129 42691 : psNode->nNumSubNodes = 0;
130 :
131 42691 : memcpy(&(psNode->rect), pRect, sizeof(CPLRectObj));
132 :
133 42691 : return psNode;
134 : }
135 :
136 : /************************************************************************/
137 : /* CPLQuadTreeCreate() */
138 : /************************************************************************/
139 :
140 : /**
141 : * Create a new quadtree
142 : *
143 : * @param pGlobalBounds a pointer to the global extent of all
144 : * the elements that will be inserted
145 : * @param pfnGetBounds a user provided function to get the bounding box of
146 : * the inserted elements. If it is set to NULL, then
147 : * CPLQuadTreeInsertWithBounds() must be used, and
148 : * extra memory will be used to keep features bounds in the
149 : * quad tree.
150 : *
151 : * @return a newly allocated quadtree
152 : */
153 :
154 239 : CPLQuadTree *CPLQuadTreeCreate(const CPLRectObj *pGlobalBounds,
155 : CPLQuadTreeGetBoundsFunc pfnGetBounds)
156 : {
157 239 : CPLAssert(pGlobalBounds);
158 :
159 : /* -------------------------------------------------------------------- */
160 : /* Allocate the hQuadTree object */
161 : /* -------------------------------------------------------------------- */
162 : CPLQuadTree *hQuadTree =
163 239 : static_cast<CPLQuadTree *>(CPLMalloc(sizeof(CPLQuadTree)));
164 :
165 239 : hQuadTree->nFeatures = 0;
166 239 : hQuadTree->pfnGetBounds = pfnGetBounds;
167 239 : hQuadTree->pfnGetBoundsEx = nullptr;
168 239 : hQuadTree->nMaxDepth = 0;
169 239 : hQuadTree->nBucketCapacity = 8;
170 :
171 239 : hQuadTree->dfSplitRatio = DEFAULT_SPLIT_RATIO;
172 239 : hQuadTree->bForceUseOfSubNodes = false;
173 :
174 : /* -------------------------------------------------------------------- */
175 : /* Allocate the psRoot psNode. */
176 : /* -------------------------------------------------------------------- */
177 239 : hQuadTree->psRoot = CPLQuadTreeNodeCreate(pGlobalBounds);
178 :
179 239 : hQuadTree->pUserData = nullptr;
180 :
181 239 : return hQuadTree;
182 : }
183 :
184 : /************************************************************************/
185 : /* CPLQuadTreeCreateEx() */
186 : /************************************************************************/
187 :
188 : /**
189 : * Create a new quadtree
190 : *
191 : * @param pGlobalBounds a pointer to the global extent of all
192 : * the elements that will be inserted
193 : * @param pfnGetBoundsEx a user provided function to get the bounding box of
194 : * the inserted elements. If it is set to NULL, then
195 : * CPLQuadTreeInsertWithBounds() must be used, and
196 : * extra memory will be used to keep features bounds in the
197 : * quad tree.
198 : * @param pUserData user data passed to pfnGetBoundsEx
199 : *
200 : * @return a newly allocated quadtree
201 : */
202 :
203 4 : CPLQuadTree *CPLQuadTreeCreateEx(const CPLRectObj *pGlobalBounds,
204 : CPLQuadTreeGetBoundsExFunc pfnGetBoundsEx,
205 : void *pUserData)
206 : {
207 4 : CPLAssert(pGlobalBounds);
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Allocate the hQuadTree object */
211 : /* -------------------------------------------------------------------- */
212 : CPLQuadTree *hQuadTree =
213 4 : static_cast<CPLQuadTree *>(CPLMalloc(sizeof(CPLQuadTree)));
214 :
215 4 : hQuadTree->nFeatures = 0;
216 4 : hQuadTree->pfnGetBounds = nullptr;
217 4 : hQuadTree->pfnGetBoundsEx = pfnGetBoundsEx;
218 4 : hQuadTree->nMaxDepth = 0;
219 4 : hQuadTree->nBucketCapacity = 8;
220 :
221 4 : hQuadTree->dfSplitRatio = DEFAULT_SPLIT_RATIO;
222 4 : hQuadTree->bForceUseOfSubNodes = false;
223 :
224 : /* -------------------------------------------------------------------- */
225 : /* Allocate the psRoot psNode. */
226 : /* -------------------------------------------------------------------- */
227 4 : hQuadTree->psRoot = CPLQuadTreeNodeCreate(pGlobalBounds);
228 :
229 4 : hQuadTree->pUserData = pUserData;
230 :
231 4 : return hQuadTree;
232 : }
233 :
234 : /************************************************************************/
235 : /* CPLQuadTreeGetAdvisedMaxDepth() */
236 : /************************************************************************/
237 :
238 : /**
239 : * Returns the optimal depth of a quadtree to hold nExpectedFeatures
240 : *
241 : * @param nExpectedFeatures the expected maximum number of elements to be
242 : * inserted.
243 : *
244 : * @return the optimal depth of a quadtree to hold nExpectedFeatures
245 : */
246 :
247 91 : int CPLQuadTreeGetAdvisedMaxDepth(int nExpectedFeatures)
248 : {
249 : /* -------------------------------------------------------------------- */
250 : /* Try to select a reasonable one */
251 : /* that implies approximately 8 shapes per node. */
252 : /* -------------------------------------------------------------------- */
253 91 : int nMaxDepth = 0;
254 91 : int nMaxNodeCount = 1;
255 :
256 106 : while (nMaxNodeCount < nExpectedFeatures / 4)
257 : {
258 15 : nMaxDepth += 1;
259 15 : nMaxNodeCount = nMaxNodeCount * 2;
260 : }
261 :
262 91 : CPLDebug("CPLQuadTree", "Estimated spatial index tree depth: %d",
263 : nMaxDepth);
264 :
265 : /* NOTE: Due to problems with memory allocation for deep trees,
266 : * automatically estimated depth is limited up to 12 levels.
267 : * See Ticket #1594 for detailed discussion.
268 : */
269 91 : if (nMaxDepth > MAX_DEFAULT_TREE_DEPTH)
270 : {
271 0 : nMaxDepth = MAX_DEFAULT_TREE_DEPTH;
272 :
273 0 : CPLDebug("CPLQuadTree",
274 : "Falling back to max number of allowed index tree "
275 : "levels (%d).",
276 : MAX_DEFAULT_TREE_DEPTH);
277 : }
278 :
279 91 : return nMaxDepth;
280 : }
281 :
282 : /************************************************************************/
283 : /* CPLQuadTreeSetMaxDepth() */
284 : /************************************************************************/
285 :
286 : /**
287 : * Set the maximum depth of a quadtree. By default, quad trees have
288 : * no maximum depth, but a maximum bucket capacity.
289 : *
290 : * @param hQuadTree the quad tree
291 : * @param nMaxDepth the maximum depth allowed
292 : */
293 :
294 91 : void CPLQuadTreeSetMaxDepth(CPLQuadTree *hQuadTree, int nMaxDepth)
295 : {
296 91 : hQuadTree->nMaxDepth = nMaxDepth;
297 91 : }
298 :
299 : /************************************************************************/
300 : /* CPLQuadTreeSetBucketCapacity() */
301 : /************************************************************************/
302 :
303 : /**
304 : * Set the maximum capacity of a node of a quadtree. The default value is 8.
305 : * Note that the maximum capacity will only be honoured if the features
306 : * inserted have a point geometry. Otherwise it may be exceeded.
307 : *
308 : * @param hQuadTree the quad tree
309 : * @param nBucketCapacity the maximum capacity of a node of a quadtree
310 : */
311 :
312 1 : void CPLQuadTreeSetBucketCapacity(CPLQuadTree *hQuadTree, int nBucketCapacity)
313 : {
314 1 : if (nBucketCapacity > 0)
315 1 : hQuadTree->nBucketCapacity = nBucketCapacity;
316 1 : }
317 :
318 : /***********************************************************************/
319 : /* CPLQuadTreeForceUseOfSubNodes() */
320 : /************************************************************************/
321 :
322 : /**
323 : * Force the quadtree to insert as much as possible a feature whose bbox
324 : * spread over multiple subnodes into those subnodes, rather than in the
325 : * list of features attached to the node.
326 : *
327 : * @param hQuadTree the quad tree
328 : */
329 :
330 4 : void CPLQuadTreeForceUseOfSubNodes(CPLQuadTree *hQuadTree)
331 : {
332 4 : hQuadTree->bForceUseOfSubNodes = true;
333 4 : }
334 :
335 : /************************************************************************/
336 : /* CPLQuadTreeInsert() */
337 : /************************************************************************/
338 :
339 : /**
340 : * Insert a feature into a quadtree
341 : *
342 : * @param hQuadTree the quad tree
343 : * @param hFeature the feature to insert
344 : */
345 :
346 29118 : void CPLQuadTreeInsert(CPLQuadTree *hQuadTree, void *hFeature)
347 : {
348 29118 : if (hQuadTree->pfnGetBounds == nullptr &&
349 1062 : hQuadTree->pfnGetBoundsEx == nullptr)
350 : {
351 0 : CPLError(CE_Failure, CPLE_AppDefined,
352 : "hQuadTree->pfnGetBounds == NULL");
353 0 : return;
354 : }
355 29118 : hQuadTree->nFeatures++;
356 : CPLRectObj bounds;
357 29118 : if (hQuadTree->pfnGetBoundsEx)
358 1062 : hQuadTree->pfnGetBoundsEx(hFeature, hQuadTree->pUserData, &bounds);
359 : else
360 28056 : hQuadTree->pfnGetBounds(hFeature, &bounds);
361 29118 : CPLQuadTreeAddFeatureInternal(hQuadTree, hFeature, &bounds);
362 : }
363 :
364 : /************************************************************************/
365 : /* CPLQuadTreeInsertWithBounds() */
366 : /************************************************************************/
367 :
368 : /**
369 : * Insert a feature into a quadtree
370 : *
371 : * @param hQuadTree the quad tree
372 : * @param hFeature the feature to insert
373 : * @param psBounds bounds of the feature
374 : */
375 94583 : void CPLQuadTreeInsertWithBounds(CPLQuadTree *hQuadTree, void *hFeature,
376 : const CPLRectObj *psBounds)
377 : {
378 94583 : hQuadTree->nFeatures++;
379 94583 : CPLQuadTreeAddFeatureInternal(hQuadTree, hFeature, psBounds);
380 94583 : }
381 :
382 : /************************************************************************/
383 : /* CPLQuadTreeRemove() */
384 : /************************************************************************/
385 :
386 23998 : static bool CPLQuadTreeRemoveInternal(QuadTreeNode *psNode, void *hFeature,
387 : const CPLRectObj *psBounds)
388 : {
389 23998 : bool bRemoved = false;
390 :
391 480592 : for (int i = 0; i < psNode->nFeatures; i++)
392 : {
393 458604 : if (psNode->pahFeatures[i] == hFeature)
394 : {
395 2010 : if (i < psNode->nFeatures - 1)
396 : {
397 1894 : memmove(psNode->pahFeatures + i, psNode->pahFeatures + i + 1,
398 1894 : (psNode->nFeatures - 1 - i) * sizeof(void *));
399 1894 : if (psNode->pasBounds)
400 : {
401 1894 : memmove(psNode->pasBounds + i, psNode->pasBounds + i + 1,
402 1894 : (psNode->nFeatures - 1 - i) * sizeof(CPLRectObj));
403 : }
404 : }
405 2010 : bRemoved = true;
406 2010 : psNode->nFeatures--;
407 2010 : break;
408 : }
409 : }
410 23998 : if (psNode->nFeatures == 0 && psNode->pahFeatures != nullptr)
411 : {
412 116 : CPLFree(psNode->pahFeatures);
413 116 : CPLFree(psNode->pasBounds);
414 116 : psNode->pahFeatures = nullptr;
415 116 : psNode->pasBounds = nullptr;
416 : }
417 :
418 : /* -------------------------------------------------------------------- */
419 : /* Recurse to subnodes if they exist. */
420 : /* -------------------------------------------------------------------- */
421 62784 : for (int i = 0; i < psNode->nNumSubNodes; i++)
422 : {
423 77572 : if (psNode->apSubNode[i] &&
424 38786 : CPL_RectOverlap(&(psNode->apSubNode[i]->rect), psBounds))
425 : {
426 21988 : bRemoved |= CPLQuadTreeRemoveInternal(psNode->apSubNode[i],
427 : hFeature, psBounds);
428 :
429 21988 : if (psNode->apSubNode[i]->nFeatures == 0 &&
430 646 : psNode->apSubNode[i]->nNumSubNodes == 0)
431 : {
432 120 : CPLQuadTreeNodeDestroy(psNode->apSubNode[i]);
433 120 : if (i < psNode->nNumSubNodes - 1)
434 : {
435 77 : memmove(psNode->apSubNode + i, psNode->apSubNode + i + 1,
436 77 : (psNode->nNumSubNodes - 1 - i) *
437 : sizeof(QuadTreeNode *));
438 : }
439 120 : i--;
440 120 : psNode->nNumSubNodes--;
441 : }
442 : }
443 : }
444 :
445 23998 : return bRemoved;
446 : }
447 :
448 : /**
449 : * Remove a feature from a quadtree.
450 : *
451 : * Currently the quadtree is not re-balanced.
452 : *
453 : * @param hQuadTree the quad tree
454 : * @param hFeature the feature to remove
455 : * @param psBounds bounds of the feature (or NULL if pfnGetBounds has been
456 : * filled)
457 : */
458 2010 : void CPLQuadTreeRemove(CPLQuadTree *hQuadTree, void *hFeature,
459 : const CPLRectObj *psBounds)
460 : {
461 2010 : if (psBounds == nullptr && hQuadTree->pfnGetBounds == nullptr &&
462 0 : hQuadTree->pfnGetBoundsEx == nullptr)
463 : {
464 0 : CPLError(CE_Failure, CPLE_AppDefined,
465 : "hQuadTree->pfnGetBounds == NULL");
466 0 : return;
467 : }
468 : CPLRectObj bounds; // keep variable in this outer scope
469 2010 : if (psBounds == nullptr)
470 : {
471 0 : if (hQuadTree->pfnGetBoundsEx)
472 0 : hQuadTree->pfnGetBoundsEx(hFeature, hQuadTree->pUserData, &bounds);
473 : else
474 0 : hQuadTree->pfnGetBounds(hFeature, &bounds);
475 0 : psBounds = &bounds;
476 : }
477 2010 : if (CPLQuadTreeRemoveInternal(hQuadTree->psRoot, hFeature, psBounds))
478 : {
479 2010 : hQuadTree->nFeatures--;
480 : }
481 : }
482 :
483 : /************************************************************************/
484 : /* CPLQuadTreeNodeDestroy() */
485 : /************************************************************************/
486 :
487 42691 : static void CPLQuadTreeNodeDestroy(QuadTreeNode *psNode)
488 : {
489 85019 : for (int i = 0; i < psNode->nNumSubNodes; i++)
490 : {
491 42328 : if (psNode->apSubNode[i])
492 42328 : CPLQuadTreeNodeDestroy(psNode->apSubNode[i]);
493 : }
494 :
495 42691 : if (psNode->pahFeatures)
496 : {
497 29824 : CPLFree(psNode->pahFeatures);
498 29824 : CPLFree(psNode->pasBounds);
499 : }
500 :
501 42691 : CPLFree(psNode);
502 42691 : }
503 :
504 : /************************************************************************/
505 : /* CPLQuadTreeDestroy() */
506 : /************************************************************************/
507 :
508 : /**
509 : * Destroy a quadtree
510 : *
511 : * @param hQuadTree the quad tree to destroy
512 : */
513 :
514 243 : void CPLQuadTreeDestroy(CPLQuadTree *hQuadTree)
515 : {
516 243 : CPLAssert(hQuadTree);
517 243 : CPLQuadTreeNodeDestroy(hQuadTree->psRoot);
518 243 : CPLFree(hQuadTree);
519 243 : }
520 :
521 : /************************************************************************/
522 : /* CPLQuadTreeSplitBounds() */
523 : /************************************************************************/
524 :
525 46584 : static void CPLQuadTreeSplitBounds(double dfSplitRatio, const CPLRectObj *in,
526 : CPLRectObj *out1, CPLRectObj *out2)
527 : {
528 : /* -------------------------------------------------------------------- */
529 : /* The output bounds will be very similar to the input bounds, */
530 : /* so just copy over to start. */
531 : /* -------------------------------------------------------------------- */
532 46584 : memcpy(out1, in, sizeof(CPLRectObj));
533 46584 : memcpy(out2, in, sizeof(CPLRectObj));
534 :
535 : /* -------------------------------------------------------------------- */
536 : /* Split in X direction. */
537 : /* -------------------------------------------------------------------- */
538 46584 : if ((in->maxx - in->minx) > (in->maxy - in->miny))
539 : {
540 20084 : const double range = in->maxx - in->minx;
541 :
542 20084 : out1->maxx = in->minx + range * dfSplitRatio;
543 20084 : out2->minx = in->maxx - range * dfSplitRatio;
544 : }
545 :
546 : /* -------------------------------------------------------------------- */
547 : /* Otherwise split in Y direction. */
548 : /* -------------------------------------------------------------------- */
549 : else
550 : {
551 26500 : const double range = in->maxy - in->miny;
552 :
553 26500 : out1->maxy = in->miny + range * dfSplitRatio;
554 26500 : out2->miny = in->maxy - range * dfSplitRatio;
555 : }
556 46584 : }
557 :
558 : /************************************************************************/
559 : /* CPLQuadTreeNodeAddFeatureAlg1() */
560 : /************************************************************************/
561 :
562 985517 : static void CPLQuadTreeNodeAddFeatureAlg1(CPLQuadTree *hQuadTree,
563 : QuadTreeNode *psNode, void *hFeature,
564 : const CPLRectObj *pRect)
565 : {
566 985517 : if (psNode->nNumSubNodes == 0)
567 : {
568 : // If we have reached the max bucket capacity, try to insert
569 : // in a subnode if possible.
570 195229 : if (psNode->nFeatures >= hQuadTree->nBucketCapacity)
571 : {
572 15526 : CPLRectObj half1 = {0.0, 0.0, 0.0, 0.0};
573 15526 : CPLRectObj half2 = {0.0, 0.0, 0.0, 0.0};
574 15526 : CPLRectObj quad1 = {0.0, 0.0, 0.0, 0.0};
575 15526 : CPLRectObj quad2 = {0.0, 0.0, 0.0, 0.0};
576 15526 : CPLRectObj quad3 = {0.0, 0.0, 0.0, 0.0};
577 15526 : CPLRectObj quad4 = {0.0, 0.0, 0.0, 0.0};
578 :
579 15526 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &psNode->rect,
580 : &half1, &half2);
581 15526 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &half1, &quad1,
582 : &quad2);
583 15526 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &half2, &quad3,
584 : &quad4);
585 :
586 46578 : if (memcmp(&psNode->rect, &quad1, sizeof(CPLRectObj)) != 0 &&
587 15526 : memcmp(&psNode->rect, &quad2, sizeof(CPLRectObj)) != 0 &&
588 15526 : memcmp(&psNode->rect, &quad3, sizeof(CPLRectObj)) != 0 &&
589 46578 : memcmp(&psNode->rect, &quad4, sizeof(CPLRectObj)) != 0 &&
590 30944 : (hQuadTree->bForceUseOfSubNodes ||
591 29757 : CPL_RectContained(pRect, &quad1) ||
592 26459 : CPL_RectContained(pRect, &quad2) ||
593 21347 : CPL_RectContained(pRect, &quad3) ||
594 9227 : CPL_RectContained(pRect, &quad4)))
595 : {
596 10610 : psNode->nNumSubNodes = 4;
597 10610 : psNode->apSubNode[0] = CPLQuadTreeNodeCreate(&quad1);
598 10610 : psNode->apSubNode[1] = CPLQuadTreeNodeCreate(&quad2);
599 10610 : psNode->apSubNode[2] = CPLQuadTreeNodeCreate(&quad3);
600 10610 : psNode->apSubNode[3] = CPLQuadTreeNodeCreate(&quad4);
601 :
602 10610 : const int oldNumFeatures = psNode->nFeatures;
603 10610 : void **oldFeatures = psNode->pahFeatures;
604 10610 : CPLRectObj *pasOldBounds = psNode->pasBounds;
605 10610 : psNode->nFeatures = 0;
606 10610 : psNode->pahFeatures = nullptr;
607 10610 : psNode->pasBounds = nullptr;
608 :
609 : // Redispatch existing pahFeatures in apSubNodes.
610 97783 : for (int i = 0; i < oldNumFeatures; i++)
611 : {
612 87173 : if (hQuadTree->pfnGetBounds == nullptr &&
613 64017 : hQuadTree->pfnGetBoundsEx == nullptr)
614 63153 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree, psNode,
615 63153 : oldFeatures[i],
616 63153 : &pasOldBounds[i]);
617 : else
618 : {
619 : CPLRectObj bounds;
620 24020 : if (hQuadTree->pfnGetBoundsEx)
621 864 : hQuadTree->pfnGetBoundsEx(
622 864 : oldFeatures[i], hQuadTree->pUserData, &bounds);
623 : else
624 23156 : hQuadTree->pfnGetBounds(oldFeatures[i], &bounds);
625 24020 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree, psNode,
626 24020 : oldFeatures[i], &bounds);
627 : }
628 : }
629 :
630 10610 : CPLFree(oldFeatures);
631 10610 : CPLFree(pasOldBounds);
632 :
633 : /* recurse back on this psNode now that it has apSubNodes */
634 10610 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree, psNode, hFeature,
635 : pRect);
636 10610 : return;
637 : }
638 : }
639 : }
640 : else
641 : {
642 : /* --------------------------------------------------------------------
643 : */
644 : /* If there are apSubNodes, then consider whether this object */
645 : /* will fit in them. */
646 : /* --------------------------------------------------------------------
647 : */
648 2169020 : for (int i = 0; i < psNode->nNumSubNodes; i++)
649 : {
650 2141990 : if (CPL_RectContained(pRect, &psNode->apSubNode[i]->rect))
651 : {
652 763251 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree, psNode->apSubNode[i],
653 : hFeature, pRect);
654 763251 : return;
655 : }
656 : }
657 27037 : if (hQuadTree->bForceUseOfSubNodes)
658 : {
659 : bool overlaps[4];
660 576 : bool overlapAll = true;
661 2880 : for (int i = 0; i < psNode->nNumSubNodes; i++)
662 : {
663 2304 : overlaps[i] =
664 2304 : CPL_RectOverlap(pRect, &psNode->apSubNode[i]->rect);
665 2304 : if (!overlaps[i])
666 1142 : overlapAll = false;
667 : }
668 576 : if (!overlapAll)
669 : {
670 2440 : for (int i = 0; i < psNode->nNumSubNodes; i++)
671 : {
672 1952 : if (overlaps[i])
673 : {
674 : CPLRectObj intersection;
675 810 : intersection.minx = std::max(
676 810 : pRect->minx, psNode->apSubNode[i]->rect.minx);
677 810 : intersection.miny = std::max(
678 810 : pRect->miny, psNode->apSubNode[i]->rect.miny);
679 810 : intersection.maxx = std::min(
680 810 : pRect->maxx, psNode->apSubNode[i]->rect.maxx);
681 810 : intersection.maxy = std::min(
682 810 : pRect->maxy, psNode->apSubNode[i]->rect.maxy);
683 810 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree,
684 : psNode->apSubNode[i],
685 : hFeature, &intersection);
686 : }
687 : }
688 488 : return;
689 : }
690 : }
691 : }
692 :
693 : /* -------------------------------------------------------------------- */
694 : /* If none of that worked, just add it to this psNodes list. */
695 : /* -------------------------------------------------------------------- */
696 211168 : psNode->nFeatures++;
697 :
698 211168 : if (psNode->nFeatures == 1)
699 : {
700 40546 : CPLAssert(psNode->pahFeatures == nullptr);
701 40546 : psNode->pahFeatures = static_cast<void **>(
702 40546 : CPLMalloc(hQuadTree->nBucketCapacity * sizeof(void *)));
703 40546 : if (hQuadTree->pfnGetBounds == nullptr &&
704 29558 : hQuadTree->pfnGetBoundsEx == nullptr)
705 29052 : psNode->pasBounds = static_cast<CPLRectObj *>(
706 29052 : CPLMalloc(hQuadTree->nBucketCapacity * sizeof(CPLRectObj)));
707 : }
708 170622 : else if (psNode->nFeatures > hQuadTree->nBucketCapacity)
709 : {
710 16454 : psNode->pahFeatures = static_cast<void **>(CPLRealloc(
711 8227 : psNode->pahFeatures, sizeof(void *) * psNode->nFeatures));
712 8227 : if (hQuadTree->pfnGetBounds == nullptr &&
713 8227 : hQuadTree->pfnGetBoundsEx == nullptr)
714 8227 : psNode->pasBounds = static_cast<CPLRectObj *>(CPLRealloc(
715 8227 : psNode->pasBounds, sizeof(CPLRectObj) * psNode->nFeatures));
716 : }
717 211168 : psNode->pahFeatures[psNode->nFeatures - 1] = hFeature;
718 211168 : if (hQuadTree->pfnGetBounds == nullptr &&
719 159956 : hQuadTree->pfnGetBoundsEx == nullptr)
720 157708 : psNode->pasBounds[psNode->nFeatures - 1] = *pRect;
721 :
722 211168 : return;
723 : }
724 :
725 : /************************************************************************/
726 : /* CPLQuadTreeNodeAddFeatureAlg2() */
727 : /************************************************************************/
728 :
729 32 : static void CPLQuadTreeNodeAddFeatureAlg2(CPLQuadTree *hQuadTree,
730 : QuadTreeNode *psNode, void *hFeature,
731 : const CPLRectObj *pRect,
732 : int nMaxDepth)
733 : {
734 : /* -------------------------------------------------------------------- */
735 : /* If there are apSubNodes, then consider whether this object */
736 : /* will fit in them. */
737 : /* -------------------------------------------------------------------- */
738 32 : if (nMaxDepth > 1 && psNode->nNumSubNodes > 0)
739 : {
740 2 : for (int i = 0; i < psNode->nNumSubNodes; i++)
741 : {
742 2 : if (CPL_RectContained(pRect, &psNode->apSubNode[i]->rect))
743 : {
744 2 : CPLQuadTreeNodeAddFeatureAlg2(hQuadTree, psNode->apSubNode[i],
745 : hFeature, pRect, nMaxDepth - 1);
746 2 : return;
747 : }
748 0 : }
749 : }
750 :
751 : /* -------------------------------------------------------------------- */
752 : /* Otherwise, consider creating four apSubNodes if could fit into */
753 : /* them, and adding to the appropriate apSubNode. */
754 : /* -------------------------------------------------------------------- */
755 30 : else if (nMaxDepth > 1 && psNode->nNumSubNodes == 0)
756 : {
757 : CPLRectObj half1, half2, quad1, quad2, quad3, quad4;
758 :
759 2 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &psNode->rect, &half1,
760 : &half2);
761 2 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &half1, &quad1, &quad2);
762 2 : CPLQuadTreeSplitBounds(hQuadTree->dfSplitRatio, &half2, &quad3, &quad4);
763 :
764 6 : if (memcmp(&psNode->rect, &quad1, sizeof(CPLRectObj)) != 0 &&
765 2 : memcmp(&psNode->rect, &quad2, sizeof(CPLRectObj)) != 0 &&
766 2 : memcmp(&psNode->rect, &quad3, sizeof(CPLRectObj)) != 0 &&
767 6 : memcmp(&psNode->rect, &quad4, sizeof(CPLRectObj)) != 0 &&
768 2 : (CPL_RectContained(pRect, &quad1) ||
769 0 : CPL_RectContained(pRect, &quad2) ||
770 0 : CPL_RectContained(pRect, &quad3) ||
771 0 : CPL_RectContained(pRect, &quad4)))
772 : {
773 2 : psNode->nNumSubNodes = 4;
774 2 : psNode->apSubNode[0] = CPLQuadTreeNodeCreate(&quad1);
775 2 : psNode->apSubNode[1] = CPLQuadTreeNodeCreate(&quad2);
776 2 : psNode->apSubNode[2] = CPLQuadTreeNodeCreate(&quad3);
777 2 : psNode->apSubNode[3] = CPLQuadTreeNodeCreate(&quad4);
778 :
779 : /* recurse back on this psNode now that it has apSubNodes */
780 2 : CPLQuadTreeNodeAddFeatureAlg2(hQuadTree, psNode, hFeature, pRect,
781 : nMaxDepth);
782 2 : return;
783 : }
784 : }
785 :
786 : /* -------------------------------------------------------------------- */
787 : /* If none of that worked, just add it to this psNodes list. */
788 : /* -------------------------------------------------------------------- */
789 28 : psNode->nFeatures++;
790 :
791 28 : psNode->pahFeatures = static_cast<void **>(
792 28 : CPLRealloc(psNode->pahFeatures, sizeof(void *) * psNode->nFeatures));
793 28 : if (hQuadTree->pfnGetBounds == nullptr &&
794 28 : hQuadTree->pfnGetBoundsEx == nullptr)
795 : {
796 28 : psNode->pasBounds = static_cast<CPLRectObj *>(CPLRealloc(
797 28 : psNode->pasBounds, sizeof(CPLRectObj) * psNode->nFeatures));
798 : }
799 28 : psNode->pahFeatures[psNode->nFeatures - 1] = hFeature;
800 28 : if (hQuadTree->pfnGetBounds == nullptr &&
801 28 : hQuadTree->pfnGetBoundsEx == nullptr)
802 : {
803 28 : psNode->pasBounds[psNode->nFeatures - 1] = *pRect;
804 : }
805 : }
806 :
807 : /************************************************************************/
808 : /* CPLQuadTreeAddFeatureInternal() */
809 : /************************************************************************/
810 :
811 123701 : static void CPLQuadTreeAddFeatureInternal(CPLQuadTree *hQuadTree,
812 : void *hFeature,
813 : const CPLRectObj *pRect)
814 : {
815 123701 : if (hQuadTree->nMaxDepth == 0)
816 : {
817 123673 : CPLQuadTreeNodeAddFeatureAlg1(hQuadTree, hQuadTree->psRoot, hFeature,
818 : pRect);
819 : }
820 : else
821 : {
822 28 : CPLQuadTreeNodeAddFeatureAlg2(hQuadTree, hQuadTree->psRoot, hFeature,
823 : pRect, hQuadTree->nMaxDepth);
824 : }
825 123701 : }
826 :
827 : /************************************************************************/
828 : /* CPLQuadTreeCollectFeatures() */
829 : /************************************************************************/
830 :
831 10521900 : static void CPLQuadTreeCollectFeatures(const CPLQuadTree *hQuadTree,
832 : const QuadTreeNode *psNode,
833 : const CPLRectObj *pAoi,
834 : int *pnFeatureCount, int *pnMaxFeatures,
835 : void ***pppFeatureList)
836 : {
837 : /* -------------------------------------------------------------------- */
838 : /* Does this psNode overlap the area of interest at all? If not, */
839 : /* return without adding to the list at all. */
840 : /* -------------------------------------------------------------------- */
841 10521900 : if (!CPL_RectOverlap(&psNode->rect, pAoi))
842 6813970 : return;
843 :
844 : /* -------------------------------------------------------------------- */
845 : /* Grow the list to hold the features on this psNode. */
846 : /* -------------------------------------------------------------------- */
847 3707540 : if (*pnFeatureCount + psNode->nFeatures > *pnMaxFeatures)
848 : {
849 : // TODO(schwehr): Symbolic constant.
850 251315 : *pnMaxFeatures = (*pnFeatureCount + psNode->nFeatures) * 2 + 20;
851 251313 : *pppFeatureList = static_cast<void **>(
852 251315 : CPLRealloc(*pppFeatureList, sizeof(void *) * *pnMaxFeatures));
853 : }
854 :
855 : /* -------------------------------------------------------------------- */
856 : /* Add the local features to the list. */
857 : /* -------------------------------------------------------------------- */
858 15430900 : for (int i = 0; i < psNode->nFeatures; i++)
859 : {
860 11722200 : if (hQuadTree->pfnGetBounds == nullptr &&
861 10991100 : hQuadTree->pfnGetBoundsEx == nullptr)
862 : {
863 10990900 : if (CPL_RectOverlap(&psNode->pasBounds[i], pAoi))
864 803632 : (*pppFeatureList)[(*pnFeatureCount)++] = psNode->pahFeatures[i];
865 : }
866 : else
867 : {
868 : CPLRectObj bounds;
869 731245 : if (hQuadTree->pfnGetBoundsEx)
870 188 : hQuadTree->pfnGetBoundsEx(psNode->pahFeatures[i],
871 188 : hQuadTree->pUserData, &bounds);
872 : else
873 731057 : hQuadTree->pfnGetBounds(psNode->pahFeatures[i], &bounds);
874 :
875 731189 : if (CPL_RectOverlap(&bounds, pAoi))
876 376559 : (*pppFeatureList)[(*pnFeatureCount)++] = psNode->pahFeatures[i];
877 : }
878 : }
879 :
880 : /* -------------------------------------------------------------------- */
881 : /* Recurse to subnodes if they exist. */
882 : /* -------------------------------------------------------------------- */
883 13986000 : for (int i = 0; i < psNode->nNumSubNodes; i++)
884 : {
885 10276900 : if (psNode->apSubNode[i])
886 10276900 : CPLQuadTreeCollectFeatures(hQuadTree, psNode->apSubNode[i], pAoi,
887 : pnFeatureCount, pnMaxFeatures,
888 : pppFeatureList);
889 : }
890 : }
891 :
892 : /************************************************************************/
893 : /* CPLQuadTreeSearch() */
894 : /************************************************************************/
895 :
896 : /**
897 : * Returns all the elements inserted whose bounding box intersects the
898 : * provided area of interest
899 : *
900 : * @param hQuadTree the quad tree
901 : * @param pAoi the pointer to the area of interest
902 : * @param pnFeatureCount the user data provided to the function.
903 : *
904 : * @return an array of features that must be freed with CPLFree
905 : */
906 :
907 245352 : void **CPLQuadTreeSearch(const CPLQuadTree *hQuadTree, const CPLRectObj *pAoi,
908 : int *pnFeatureCount)
909 : {
910 245352 : CPLAssert(hQuadTree);
911 245352 : CPLAssert(pAoi);
912 :
913 245352 : int nFeatureCount = 0;
914 245352 : if (pnFeatureCount == nullptr)
915 0 : pnFeatureCount = &nFeatureCount;
916 :
917 245352 : *pnFeatureCount = 0;
918 :
919 245352 : int nMaxFeatures = 0;
920 245352 : void **ppFeatureList = nullptr;
921 245352 : CPLQuadTreeCollectFeatures(hQuadTree, hQuadTree->psRoot, pAoi,
922 : pnFeatureCount, &nMaxFeatures, &ppFeatureList);
923 :
924 245358 : return ppFeatureList;
925 : }
926 :
927 : /************************************************************************/
928 : /* CPLQuadTreeNodeForeach() */
929 : /************************************************************************/
930 :
931 25 : static bool CPLQuadTreeNodeForeach(const QuadTreeNode *psNode,
932 : CPLQuadTreeForeachFunc pfnForeach,
933 : void *pUserData)
934 : {
935 49 : for (int i = 0; i < psNode->nNumSubNodes; i++)
936 : {
937 24 : if (!CPLQuadTreeNodeForeach(psNode->apSubNode[i], pfnForeach,
938 : pUserData))
939 0 : return false;
940 : }
941 :
942 50 : for (int i = 0; i < psNode->nFeatures; i++)
943 : {
944 25 : if (pfnForeach(psNode->pahFeatures[i], pUserData) == FALSE)
945 0 : return false;
946 : }
947 :
948 25 : return true;
949 : }
950 :
951 : /************************************************************************/
952 : /* CPLQuadTreeForeach() */
953 : /************************************************************************/
954 :
955 : /**
956 : * Walk through the quadtree and runs the provided function on all the
957 : * elements
958 : *
959 : * This function is provided with the user_data argument of pfnForeach.
960 : * It must return TRUE to go on the walk through the hash set, or FALSE to
961 : * make it stop.
962 : *
963 : * Note : the structure of the quadtree must *NOT* be modified during the
964 : * walk.
965 : *
966 : * @param hQuadTree the quad tree
967 : * @param pfnForeach the function called on each element.
968 : * @param pUserData the user data provided to the function.
969 : */
970 :
971 1 : void CPLQuadTreeForeach(const CPLQuadTree *hQuadTree,
972 : CPLQuadTreeForeachFunc pfnForeach, void *pUserData)
973 : {
974 1 : CPLAssert(hQuadTree);
975 1 : CPLAssert(pfnForeach);
976 1 : CPLQuadTreeNodeForeach(hQuadTree->psRoot, pfnForeach, pUserData);
977 1 : }
978 :
979 : /************************************************************************/
980 : /* CPLQuadTreeDumpNode() */
981 : /************************************************************************/
982 :
983 0 : static void CPLQuadTreeDumpNode(const QuadTreeNode *psNode, int nIndentLevel,
984 : CPLQuadTreeDumpFeatureFunc pfnDumpFeatureFunc,
985 : void *pUserData)
986 : {
987 0 : if (psNode->nNumSubNodes)
988 : {
989 0 : for (int count = nIndentLevel; --count >= 0;)
990 : {
991 0 : printf(" "); /*ok*/
992 : }
993 0 : printf("SubhQuadTrees :\n"); /*ok*/
994 0 : for (int i = 0; i < psNode->nNumSubNodes; i++)
995 : {
996 0 : for (int count = nIndentLevel + 1; --count >= 0;)
997 : {
998 0 : printf(" "); /*ok*/
999 : }
1000 0 : printf("SubhQuadTree %d :\n", i + 1); /*ok*/
1001 0 : CPLQuadTreeDumpNode(psNode->apSubNode[i], nIndentLevel + 2,
1002 : pfnDumpFeatureFunc, pUserData);
1003 : }
1004 : }
1005 0 : if (psNode->nFeatures)
1006 : {
1007 0 : for (int count = nIndentLevel; --count >= 0;)
1008 0 : printf(" "); /*ok*/
1009 0 : printf("Leaves (%d):\n", psNode->nFeatures); /*ok*/
1010 0 : for (int i = 0; i < psNode->nFeatures; i++)
1011 : {
1012 0 : if (pfnDumpFeatureFunc)
1013 : {
1014 0 : pfnDumpFeatureFunc(psNode->pahFeatures[i], nIndentLevel + 2,
1015 : pUserData);
1016 : }
1017 : else
1018 : {
1019 0 : for (int count = nIndentLevel + 1; --count >= 0;)
1020 : {
1021 0 : printf(" "); /*ok*/
1022 : }
1023 0 : printf("%p\n", psNode->pahFeatures[i]); /*ok*/
1024 : }
1025 : }
1026 : }
1027 0 : }
1028 :
1029 : /************************************************************************/
1030 : /* CPLQuadTreeDump() */
1031 : /************************************************************************/
1032 :
1033 : /** Dump quad tree */
1034 0 : void CPLQuadTreeDump(const CPLQuadTree *hQuadTree,
1035 : CPLQuadTreeDumpFeatureFunc pfnDumpFeatureFunc,
1036 : void *pUserData)
1037 : {
1038 0 : CPLQuadTreeDumpNode(hQuadTree->psRoot, 0, pfnDumpFeatureFunc, pUserData);
1039 0 : }
1040 :
1041 : /************************************************************************/
1042 : /* CPLQuadTreeGetStatsNode() */
1043 : /************************************************************************/
1044 :
1045 0 : static void CPLQuadTreeGetStatsNode(const QuadTreeNode *psNode, int nDepthLevel,
1046 : int *pnNodeCount, int *pnMaxDepth,
1047 : int *pnMaxBucketCapacity)
1048 : {
1049 0 : (*pnNodeCount)++;
1050 0 : if (nDepthLevel > *pnMaxDepth)
1051 0 : *pnMaxDepth = nDepthLevel;
1052 0 : if (psNode->nFeatures > *pnMaxBucketCapacity)
1053 0 : *pnMaxBucketCapacity = psNode->nFeatures;
1054 :
1055 0 : for (int i = 0; i < psNode->nNumSubNodes; i++)
1056 : {
1057 0 : CPLQuadTreeGetStatsNode(psNode->apSubNode[i], nDepthLevel + 1,
1058 : pnNodeCount, pnMaxDepth, pnMaxBucketCapacity);
1059 : }
1060 0 : }
1061 :
1062 : /************************************************************************/
1063 : /* CPLQuadTreeGetStats() */
1064 : /************************************************************************/
1065 :
1066 : /** Get stats */
1067 0 : void CPLQuadTreeGetStats(const CPLQuadTree *hQuadTree, int *pnFeatureCount,
1068 : int *pnNodeCount, int *pnMaxDepth,
1069 : int *pnMaxBucketCapacity)
1070 : {
1071 0 : CPLAssert(hQuadTree);
1072 :
1073 0 : int nFeatureCount = 0;
1074 0 : if (pnFeatureCount == nullptr)
1075 0 : pnFeatureCount = &nFeatureCount;
1076 0 : int nNodeCount = 0;
1077 0 : if (pnNodeCount == nullptr)
1078 0 : pnNodeCount = &nNodeCount;
1079 0 : int nMaxDepth = 0;
1080 0 : if (pnMaxDepth == nullptr)
1081 0 : pnMaxDepth = &nMaxDepth;
1082 0 : int nMaxBucketCapacity = 0;
1083 0 : if (pnMaxBucketCapacity == nullptr)
1084 0 : pnMaxBucketCapacity = &nMaxBucketCapacity;
1085 :
1086 0 : *pnFeatureCount = hQuadTree->nFeatures;
1087 0 : *pnNodeCount = 0;
1088 0 : *pnMaxDepth = 1;
1089 0 : *pnMaxBucketCapacity = 0;
1090 :
1091 0 : CPLQuadTreeGetStatsNode(hQuadTree->psRoot, 0, pnNodeCount, pnMaxDepth,
1092 : pnMaxBucketCapacity);
1093 :
1094 : // TODO(schwehr): If any of the pointers were set to local vars,
1095 : // do they need to be reset to a nullptr?
1096 0 : }
|
