Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Shapelib
4 : * Purpose: Implementation of quadtree building and searching functions.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1999, Frank Warmerdam
9 : * Copyright (c) 2012-2024, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT OR LGPL-2.0-or-later
12 : ******************************************************************************
13 : *
14 : */
15 :
16 : #include "shapefil_private.h"
17 :
18 : #include <math.h>
19 : #include <assert.h>
20 : #include <stdbool.h>
21 : #include <stdlib.h>
22 : #include <string.h>
23 : #include <limits.h>
24 :
25 : #ifdef USE_CPL
26 : #include "cpl_error.h"
27 : #endif
28 :
29 : #ifndef TRUE
30 : #define TRUE 1
31 : #define FALSE 0
32 : #endif
33 :
34 : /* -------------------------------------------------------------------- */
35 : /* If the following is 0.5, nodes will be split in half. If it */
36 : /* is 0.6 then each subnode will contain 60% of the parent */
37 : /* node, with 20% representing overlap. This can be help to */
38 : /* prevent small objects on a boundary from shifting too high */
39 : /* up the tree. */
40 : /* -------------------------------------------------------------------- */
41 :
42 : #define SHP_SPLIT_RATIO 0.55
43 :
44 : /************************************************************************/
45 : /* SHPTreeNodeInit() */
46 : /* */
47 : /* Initialize a tree node. */
48 : /************************************************************************/
49 :
50 740917 : static SHPTreeNode *SHPTreeNodeCreate(const double *padfBoundsMin,
51 : const double *padfBoundsMax)
52 :
53 : {
54 : SHPTreeNode *psTreeNode;
55 :
56 740917 : psTreeNode = STATIC_CAST(SHPTreeNode *, malloc(sizeof(SHPTreeNode)));
57 740917 : if (SHPLIB_NULLPTR == psTreeNode)
58 0 : return SHPLIB_NULLPTR;
59 :
60 740917 : psTreeNode->nShapeCount = 0;
61 740917 : psTreeNode->panShapeIds = SHPLIB_NULLPTR;
62 740917 : psTreeNode->papsShapeObj = SHPLIB_NULLPTR;
63 :
64 740917 : psTreeNode->nSubNodes = 0;
65 :
66 740917 : if (padfBoundsMin != SHPLIB_NULLPTR)
67 740900 : memcpy(psTreeNode->adfBoundsMin, padfBoundsMin, sizeof(double) * 4);
68 :
69 740917 : if (padfBoundsMax != SHPLIB_NULLPTR)
70 740900 : memcpy(psTreeNode->adfBoundsMax, padfBoundsMax, sizeof(double) * 4);
71 :
72 740917 : return psTreeNode;
73 : }
74 :
75 : /************************************************************************/
76 : /* SHPCreateTree() */
77 : /************************************************************************/
78 :
79 : SHPTree SHPAPI_CALL1(*)
80 17 : SHPCreateTree(SHPHandle hSHP, int nDimension, int nMaxDepth,
81 : const double *padfBoundsMin, const double *padfBoundsMax)
82 :
83 : {
84 : SHPTree *psTree;
85 :
86 17 : if (padfBoundsMin == SHPLIB_NULLPTR && hSHP == SHPLIB_NULLPTR)
87 0 : return SHPLIB_NULLPTR;
88 :
89 : /* -------------------------------------------------------------------- */
90 : /* Allocate the tree object */
91 : /* -------------------------------------------------------------------- */
92 17 : psTree = STATIC_CAST(SHPTree *, malloc(sizeof(SHPTree)));
93 17 : if (SHPLIB_NULLPTR == psTree)
94 : {
95 0 : return SHPLIB_NULLPTR;
96 : }
97 :
98 17 : psTree->hSHP = hSHP;
99 17 : psTree->nMaxDepth = nMaxDepth;
100 17 : psTree->nDimension = nDimension;
101 17 : psTree->nTotalCount = 0;
102 :
103 : /* -------------------------------------------------------------------- */
104 : /* If no max depth was defined, try to select a reasonable one */
105 : /* that implies approximately 8 shapes per node. */
106 : /* -------------------------------------------------------------------- */
107 17 : if (psTree->nMaxDepth == 0 && hSHP != SHPLIB_NULLPTR)
108 : {
109 17 : int nMaxNodeCount = 1;
110 : int nShapeCount;
111 :
112 17 : SHPGetInfo(hSHP, &nShapeCount, SHPLIB_NULLPTR, SHPLIB_NULLPTR,
113 : SHPLIB_NULLPTR);
114 95 : while (nMaxNodeCount * 4 < nShapeCount)
115 : {
116 78 : psTree->nMaxDepth += 1;
117 78 : nMaxNodeCount = nMaxNodeCount * 2;
118 : }
119 :
120 : #ifdef USE_CPL
121 17 : CPLDebug("Shape", "Estimated spatial index tree depth: %d",
122 : psTree->nMaxDepth);
123 : #endif
124 :
125 : /* NOTE: Due to problems with memory allocation for deep trees,
126 : * automatically estimated depth is limited up to 12 levels.
127 : * See Ticket #1594 for detailed discussion.
128 : */
129 17 : if (psTree->nMaxDepth > MAX_DEFAULT_TREE_DEPTH)
130 : {
131 0 : psTree->nMaxDepth = MAX_DEFAULT_TREE_DEPTH;
132 :
133 : #ifdef USE_CPL
134 0 : CPLDebug(
135 : "Shape",
136 : "Falling back to max number of allowed index tree levels (%d).",
137 : MAX_DEFAULT_TREE_DEPTH);
138 : #endif
139 : }
140 : }
141 :
142 : /* -------------------------------------------------------------------- */
143 : /* Allocate the root node. */
144 : /* -------------------------------------------------------------------- */
145 17 : psTree->psRoot = SHPTreeNodeCreate(padfBoundsMin, padfBoundsMax);
146 17 : if (SHPLIB_NULLPTR == psTree->psRoot)
147 : {
148 0 : free(psTree);
149 0 : return SHPLIB_NULLPTR;
150 : }
151 :
152 : /* -------------------------------------------------------------------- */
153 : /* Assign the bounds to the root node. If none are passed in, */
154 : /* use the bounds of the provided file otherwise the create */
155 : /* function will have already set the bounds. */
156 : /* -------------------------------------------------------------------- */
157 17 : if (padfBoundsMin == SHPLIB_NULLPTR)
158 : {
159 17 : SHPGetInfo(hSHP, SHPLIB_NULLPTR, SHPLIB_NULLPTR,
160 17 : psTree->psRoot->adfBoundsMin, psTree->psRoot->adfBoundsMax);
161 : }
162 :
163 : /* -------------------------------------------------------------------- */
164 : /* If we have a file, insert all its shapes into the tree. */
165 : /* -------------------------------------------------------------------- */
166 17 : if (hSHP != SHPLIB_NULLPTR)
167 : {
168 : int iShape, nShapeCount;
169 :
170 17 : SHPGetInfo(hSHP, &nShapeCount, SHPLIB_NULLPTR, SHPLIB_NULLPTR,
171 : SHPLIB_NULLPTR);
172 :
173 56322 : for (iShape = 0; iShape < nShapeCount; iShape++)
174 : {
175 : SHPObject *psShape;
176 :
177 56305 : psShape = SHPReadObject(hSHP, iShape);
178 56305 : if (psShape != SHPLIB_NULLPTR)
179 : {
180 56305 : SHPTreeAddShapeId(psTree, psShape);
181 56305 : SHPDestroyObject(psShape);
182 : }
183 : }
184 : }
185 :
186 17 : return psTree;
187 : }
188 :
189 : /************************************************************************/
190 : /* SHPDestroyTreeNode() */
191 : /************************************************************************/
192 :
193 586664 : static void SHPDestroyTreeNode(SHPTreeNode *psTreeNode)
194 :
195 : {
196 : int i;
197 :
198 586664 : assert(SHPLIB_NULLPTR != psTreeNode);
199 :
200 672600 : for (i = 0; i < psTreeNode->nSubNodes; i++)
201 : {
202 85936 : if (psTreeNode->apsSubNode[i] != SHPLIB_NULLPTR)
203 85936 : SHPDestroyTreeNode(psTreeNode->apsSubNode[i]);
204 : }
205 :
206 586664 : if (psTreeNode->panShapeIds != SHPLIB_NULLPTR)
207 54983 : free(psTreeNode->panShapeIds);
208 :
209 586664 : if (psTreeNode->papsShapeObj != SHPLIB_NULLPTR)
210 : {
211 0 : for (i = 0; i < psTreeNode->nShapeCount; i++)
212 : {
213 0 : if (psTreeNode->papsShapeObj[i] != SHPLIB_NULLPTR)
214 0 : SHPDestroyObject(psTreeNode->papsShapeObj[i]);
215 : }
216 :
217 0 : free(psTreeNode->papsShapeObj);
218 : }
219 :
220 586664 : free(psTreeNode);
221 586664 : }
222 :
223 : /************************************************************************/
224 : /* SHPDestroyTree() */
225 : /************************************************************************/
226 :
227 17 : void SHPAPI_CALL SHPDestroyTree(SHPTree *psTree)
228 :
229 : {
230 17 : SHPDestroyTreeNode(psTree->psRoot);
231 17 : free(psTree);
232 17 : }
233 :
234 : /************************************************************************/
235 : /* SHPCheckBoundsOverlap() */
236 : /* */
237 : /* Do the given boxes overlap at all? */
238 : /************************************************************************/
239 :
240 525416 : int SHPAPI_CALL SHPCheckBoundsOverlap(const double *padfBox1Min,
241 : const double *padfBox1Max,
242 : const double *padfBox2Min,
243 : const double *padfBox2Max, int nDimension)
244 : {
245 988737 : for (int iDim = 0; iDim < nDimension; iDim++)
246 : {
247 817216 : if (padfBox2Max[iDim] < padfBox1Min[iDim])
248 199983 : return FALSE;
249 :
250 617233 : if (padfBox1Max[iDim] < padfBox2Min[iDim])
251 153912 : return FALSE;
252 : }
253 :
254 171521 : return TRUE;
255 : }
256 :
257 : /************************************************************************/
258 : /* SHPCheckObjectContained() */
259 : /* */
260 : /* Does the given shape fit within the indicated extents? */
261 : /************************************************************************/
262 :
263 2226510 : static bool SHPCheckObjectContained(const SHPObject *psObject, int nDimension,
264 : const double *padfBoundsMin,
265 : const double *padfBoundsMax)
266 :
267 : {
268 2226510 : if (psObject->dfXMin < padfBoundsMin[0] ||
269 2088720 : psObject->dfXMax > padfBoundsMax[0])
270 922534 : return false;
271 :
272 1303980 : if (psObject->dfYMin < padfBoundsMin[1] ||
273 1300830 : psObject->dfYMax > padfBoundsMax[1])
274 531627 : return false;
275 :
276 772349 : if (nDimension == 2)
277 772349 : return true;
278 :
279 0 : if (psObject->dfZMin < padfBoundsMin[2] ||
280 0 : psObject->dfZMax > padfBoundsMax[2])
281 0 : return false;
282 :
283 0 : if (nDimension == 3)
284 0 : return true;
285 :
286 0 : if (psObject->dfMMin < padfBoundsMin[3] ||
287 0 : psObject->dfMMax > padfBoundsMax[3])
288 0 : return false;
289 :
290 0 : return true;
291 : }
292 :
293 : /************************************************************************/
294 : /* SHPTreeSplitBounds() */
295 : /* */
296 : /* Split a region into two subregion evenly, cutting along the */
297 : /* longest dimension. */
298 : /************************************************************************/
299 :
300 585681 : static void SHPTreeSplitBounds(const double *padfBoundsMinIn,
301 : const double *padfBoundsMaxIn,
302 : double *padfBoundsMin1, double *padfBoundsMax1,
303 : double *padfBoundsMin2, double *padfBoundsMax2)
304 :
305 : {
306 : /* -------------------------------------------------------------------- */
307 : /* The output bounds will be very similar to the input bounds, */
308 : /* so just copy over to start. */
309 : /* -------------------------------------------------------------------- */
310 585681 : memcpy(padfBoundsMin1, padfBoundsMinIn, sizeof(double) * 4);
311 585681 : memcpy(padfBoundsMax1, padfBoundsMaxIn, sizeof(double) * 4);
312 585681 : memcpy(padfBoundsMin2, padfBoundsMinIn, sizeof(double) * 4);
313 585681 : memcpy(padfBoundsMax2, padfBoundsMaxIn, sizeof(double) * 4);
314 :
315 : /* -------------------------------------------------------------------- */
316 : /* Split in X direction. */
317 : /* -------------------------------------------------------------------- */
318 585681 : if ((padfBoundsMaxIn[0] - padfBoundsMinIn[0]) >
319 585681 : (padfBoundsMaxIn[1] - padfBoundsMinIn[1]))
320 : {
321 383299 : double dfRange = padfBoundsMaxIn[0] - padfBoundsMinIn[0];
322 :
323 383299 : padfBoundsMax1[0] = padfBoundsMinIn[0] + dfRange * SHP_SPLIT_RATIO;
324 383299 : padfBoundsMin2[0] = padfBoundsMaxIn[0] - dfRange * SHP_SPLIT_RATIO;
325 : }
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* Otherwise split in Y direction. */
329 : /* -------------------------------------------------------------------- */
330 : else
331 : {
332 202382 : double dfRange = padfBoundsMaxIn[1] - padfBoundsMinIn[1];
333 :
334 202382 : padfBoundsMax1[1] = padfBoundsMinIn[1] + dfRange * SHP_SPLIT_RATIO;
335 202382 : padfBoundsMin2[1] = padfBoundsMaxIn[1] - dfRange * SHP_SPLIT_RATIO;
336 : }
337 585681 : }
338 :
339 : /************************************************************************/
340 : /* SHPTreeNodeAddShapeId() */
341 : /************************************************************************/
342 :
343 828654 : static bool SHPTreeNodeAddShapeId(SHPTreeNode *psTreeNode, SHPObject *psObject,
344 : int nMaxDepth, int nDimension)
345 :
346 : {
347 : int i;
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* If there are subnodes, then consider whether this object */
351 : /* will fit in them. */
352 : /* -------------------------------------------------------------------- */
353 828654 : if (nMaxDepth > 1 && psTreeNode->nSubNodes > 0)
354 : {
355 1626980 : for (i = 0; i < psTreeNode->nSubNodes; i++)
356 : {
357 1626980 : if (SHPCheckObjectContained(
358 : psObject, nDimension,
359 1626980 : psTreeNode->apsSubNode[i]->adfBoundsMin,
360 1626980 : psTreeNode->apsSubNode[i]->adfBoundsMax))
361 : {
362 587124 : return SHPTreeNodeAddShapeId(psTreeNode->apsSubNode[i],
363 : psObject, nMaxDepth - 1,
364 587124 : nDimension);
365 : }
366 : }
367 : }
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Otherwise, consider creating four subnodes if could fit into */
371 : /* them, and adding to the appropriate subnode. */
372 : /* -------------------------------------------------------------------- */
373 : #if MAX_SUBNODE == 4
374 241526 : else if (nMaxDepth > 1 && psTreeNode->nSubNodes == 0)
375 : {
376 : double adfBoundsMinH1[4], adfBoundsMaxH1[4];
377 : double adfBoundsMinH2[4], adfBoundsMaxH2[4];
378 : double adfBoundsMin1[4], adfBoundsMax1[4];
379 : double adfBoundsMin2[4], adfBoundsMax2[4];
380 : double adfBoundsMin3[4], adfBoundsMax3[4];
381 : double adfBoundsMin4[4], adfBoundsMax4[4];
382 :
383 195227 : SHPTreeSplitBounds(psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
384 : adfBoundsMinH1, adfBoundsMaxH1, adfBoundsMinH2,
385 : adfBoundsMaxH2);
386 :
387 195227 : SHPTreeSplitBounds(adfBoundsMinH1, adfBoundsMaxH1, adfBoundsMin1,
388 : adfBoundsMax1, adfBoundsMin2, adfBoundsMax2);
389 :
390 195227 : SHPTreeSplitBounds(adfBoundsMinH2, adfBoundsMaxH2, adfBoundsMin3,
391 : adfBoundsMax3, adfBoundsMin4, adfBoundsMax4);
392 :
393 195227 : if (SHPCheckObjectContained(psObject, nDimension, adfBoundsMin1,
394 174771 : adfBoundsMax1) ||
395 174771 : SHPCheckObjectContained(psObject, nDimension, adfBoundsMin2,
396 140387 : adfBoundsMax2) ||
397 140387 : SHPCheckObjectContained(psObject, nDimension, adfBoundsMin3,
398 369998 : adfBoundsMax3) ||
399 89146 : SHPCheckObjectContained(psObject, nDimension, adfBoundsMin4,
400 : adfBoundsMax4))
401 : {
402 185225 : psTreeNode->nSubNodes = 4;
403 185225 : psTreeNode->apsSubNode[0] =
404 185225 : SHPTreeNodeCreate(adfBoundsMin1, adfBoundsMax1);
405 185225 : psTreeNode->apsSubNode[1] =
406 185225 : SHPTreeNodeCreate(adfBoundsMin2, adfBoundsMax2);
407 185225 : psTreeNode->apsSubNode[2] =
408 185225 : SHPTreeNodeCreate(adfBoundsMin3, adfBoundsMax3);
409 185225 : psTreeNode->apsSubNode[3] =
410 185225 : SHPTreeNodeCreate(adfBoundsMin4, adfBoundsMax4);
411 :
412 : /* recurse back on this node now that it has subnodes */
413 185225 : return (SHPTreeNodeAddShapeId(psTreeNode, psObject, nMaxDepth,
414 185225 : nDimension));
415 : }
416 : }
417 : #endif /* MAX_SUBNODE == 4 */
418 :
419 : /* -------------------------------------------------------------------- */
420 : /* Otherwise, consider creating two subnodes if could fit into */
421 : /* them, and adding to the appropriate subnode. */
422 : /* -------------------------------------------------------------------- */
423 : #if MAX_SUBNODE == 2
424 : else if (nMaxDepth > 1 && psTreeNode->nSubNodes == 0)
425 : {
426 : double adfBoundsMin1[4], adfBoundsMax1[4];
427 : double adfBoundsMin2[4], adfBoundsMax2[4];
428 :
429 : SHPTreeSplitBounds(psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
430 : adfBoundsMin1, adfBoundsMax1, adfBoundsMin2,
431 : adfBoundsMax2);
432 :
433 : if (SHPCheckObjectContained(psObject, nDimension, adfBoundsMin1,
434 : adfBoundsMax1))
435 : {
436 : psTreeNode->nSubNodes = 2;
437 : psTreeNode->apsSubNode[0] =
438 : SHPTreeNodeCreate(adfBoundsMin1, adfBoundsMax1);
439 : psTreeNode->apsSubNode[1] =
440 : SHPTreeNodeCreate(adfBoundsMin2, adfBoundsMax2);
441 :
442 : return (SHPTreeNodeAddShapeId(psTreeNode->apsSubNode[0], psObject,
443 : nMaxDepth - 1, nDimension));
444 : }
445 : else if (SHPCheckObjectContained(psObject, nDimension, adfBoundsMin2,
446 : adfBoundsMax2))
447 : {
448 : psTreeNode->nSubNodes = 2;
449 : psTreeNode->apsSubNode[0] =
450 : SHPTreeNodeCreate(adfBoundsMin1, adfBoundsMax1);
451 : psTreeNode->apsSubNode[1] =
452 : SHPTreeNodeCreate(adfBoundsMin2, adfBoundsMax2);
453 :
454 : return (SHPTreeNodeAddShapeId(psTreeNode->apsSubNode[1], psObject,
455 : nMaxDepth - 1, nDimension));
456 : }
457 : }
458 : #endif /* MAX_SUBNODE == 2 */
459 :
460 : /* -------------------------------------------------------------------- */
461 : /* If none of that worked, just add it to this nodes list. */
462 : /* -------------------------------------------------------------------- */
463 56305 : psTreeNode->nShapeCount++;
464 :
465 56305 : psTreeNode->panShapeIds =
466 56305 : STATIC_CAST(int *, realloc(psTreeNode->panShapeIds,
467 : sizeof(int) * psTreeNode->nShapeCount));
468 56305 : psTreeNode->panShapeIds[psTreeNode->nShapeCount - 1] = psObject->nShapeId;
469 :
470 56305 : if (psTreeNode->papsShapeObj != SHPLIB_NULLPTR)
471 : {
472 0 : psTreeNode->papsShapeObj = STATIC_CAST(
473 : SHPObject **, realloc(psTreeNode->papsShapeObj,
474 : sizeof(void *) * psTreeNode->nShapeCount));
475 0 : psTreeNode->papsShapeObj[psTreeNode->nShapeCount - 1] = SHPLIB_NULLPTR;
476 : }
477 :
478 56305 : return true;
479 : }
480 :
481 : /************************************************************************/
482 : /* SHPTreeAddShapeId() */
483 : /* */
484 : /* Add a shape to the tree, but don't keep a pointer to the */
485 : /* object data, just keep the shapeid. */
486 : /************************************************************************/
487 :
488 56305 : int SHPAPI_CALL SHPTreeAddShapeId(SHPTree *psTree, SHPObject *psObject)
489 :
490 : {
491 56305 : psTree->nTotalCount++;
492 :
493 56305 : return (SHPTreeNodeAddShapeId(psTree->psRoot, psObject, psTree->nMaxDepth,
494 56305 : psTree->nDimension));
495 : }
496 :
497 : /************************************************************************/
498 : /* SHPTreeCollectShapesIds() */
499 : /* */
500 : /* Work function implementing SHPTreeFindLikelyShapes() on a */
501 : /* tree node by tree node basis. */
502 : /************************************************************************/
503 :
504 0 : static void SHPTreeCollectShapeIds(const SHPTree *hTree,
505 : const SHPTreeNode *psTreeNode,
506 : double *padfBoundsMin, double *padfBoundsMax,
507 : int *pnShapeCount, int *pnMaxShapes,
508 : int **ppanShapeList)
509 :
510 : {
511 : int i;
512 :
513 : /* -------------------------------------------------------------------- */
514 : /* Does this node overlap the area of interest at all? If not, */
515 : /* return without adding to the list at all. */
516 : /* -------------------------------------------------------------------- */
517 0 : if (!SHPCheckBoundsOverlap(psTreeNode->adfBoundsMin,
518 0 : psTreeNode->adfBoundsMax, padfBoundsMin,
519 0 : padfBoundsMax, hTree->nDimension))
520 0 : return;
521 :
522 : /* -------------------------------------------------------------------- */
523 : /* Grow the list to hold the shapes on this node. */
524 : /* -------------------------------------------------------------------- */
525 0 : if (*pnShapeCount + psTreeNode->nShapeCount > *pnMaxShapes)
526 : {
527 0 : *pnMaxShapes = (*pnShapeCount + psTreeNode->nShapeCount) * 2 + 20;
528 0 : *ppanShapeList = STATIC_CAST(
529 : int *, realloc(*ppanShapeList, sizeof(int) * *pnMaxShapes));
530 : }
531 :
532 : /* -------------------------------------------------------------------- */
533 : /* Add the local nodes shapeids to the list. */
534 : /* -------------------------------------------------------------------- */
535 0 : for (i = 0; i < psTreeNode->nShapeCount; i++)
536 : {
537 0 : (*ppanShapeList)[(*pnShapeCount)++] = psTreeNode->panShapeIds[i];
538 : }
539 :
540 : /* -------------------------------------------------------------------- */
541 : /* Recurse to subnodes if they exist. */
542 : /* -------------------------------------------------------------------- */
543 0 : for (i = 0; i < psTreeNode->nSubNodes; i++)
544 : {
545 0 : if (psTreeNode->apsSubNode[i] != SHPLIB_NULLPTR)
546 0 : SHPTreeCollectShapeIds(hTree, psTreeNode->apsSubNode[i],
547 : padfBoundsMin, padfBoundsMax, pnShapeCount,
548 : pnMaxShapes, ppanShapeList);
549 : }
550 : }
551 :
552 : /************************************************************************/
553 : /* SHPTreeFindLikelyShapes() */
554 : /* */
555 : /* Find all shapes within tree nodes for which the tree node */
556 : /* bounding box overlaps the search box. The return value is */
557 : /* an array of shapeids terminated by a -1. The shapeids will */
558 : /* be in order, as hopefully this will result in faster (more */
559 : /* sequential) reading from the file. */
560 : /************************************************************************/
561 :
562 : /* helper for qsort */
563 919395 : static int SHPTreeCompareInts(const void *a, const void *b)
564 : {
565 919395 : return *REINTERPRET_CAST(const int *, a) -
566 919395 : *REINTERPRET_CAST(const int *, b);
567 : }
568 :
569 : int SHPAPI_CALL1(*)
570 0 : SHPTreeFindLikelyShapes(const SHPTree *hTree, double *padfBoundsMin,
571 : double *padfBoundsMax, int *pnShapeCount)
572 :
573 : {
574 0 : int *panShapeList = SHPLIB_NULLPTR, nMaxShapes = 0;
575 :
576 : /* -------------------------------------------------------------------- */
577 : /* Perform the search by recursive descent. */
578 : /* -------------------------------------------------------------------- */
579 0 : *pnShapeCount = 0;
580 :
581 0 : SHPTreeCollectShapeIds(hTree, hTree->psRoot, padfBoundsMin, padfBoundsMax,
582 : pnShapeCount, &nMaxShapes, &panShapeList);
583 :
584 : /* -------------------------------------------------------------------- */
585 : /* Sort the id array */
586 : /* -------------------------------------------------------------------- */
587 :
588 0 : if (panShapeList != SHPLIB_NULLPTR)
589 0 : qsort(panShapeList, *pnShapeCount, sizeof(int), SHPTreeCompareInts);
590 :
591 0 : return panShapeList;
592 : }
593 :
594 : /************************************************************************/
595 : /* SHPTreeNodeTrim() */
596 : /* */
597 : /* This is the recursive version of SHPTreeTrimExtraNodes() that */
598 : /* walks the tree cleaning it up. */
599 : /************************************************************************/
600 :
601 740917 : static int SHPTreeNodeTrim(SHPTreeNode *psTreeNode)
602 : {
603 : int i;
604 :
605 : /* -------------------------------------------------------------------- */
606 : /* Trim subtrees, and free subnodes that come back empty. */
607 : /* -------------------------------------------------------------------- */
608 1481820 : for (i = 0; i < psTreeNode->nSubNodes; i++)
609 : {
610 740900 : if (SHPTreeNodeTrim(psTreeNode->apsSubNode[i]))
611 : {
612 500711 : SHPDestroyTreeNode(psTreeNode->apsSubNode[i]);
613 :
614 500711 : psTreeNode->apsSubNode[i] =
615 500711 : psTreeNode->apsSubNode[psTreeNode->nSubNodes - 1];
616 :
617 500711 : psTreeNode->nSubNodes--;
618 :
619 500711 : i--; /* process the new occupant of this subnode entry */
620 : }
621 : }
622 :
623 : /* -------------------------------------------------------------------- */
624 : /* If the current node has 1 subnode and no shapes, promote that */
625 : /* subnode to the current node position. */
626 : /* -------------------------------------------------------------------- */
627 740917 : if (psTreeNode->nSubNodes == 1 && psTreeNode->nShapeCount == 0)
628 : {
629 154253 : SHPTreeNode *psSubNode = psTreeNode->apsSubNode[0];
630 :
631 154253 : memcpy(psTreeNode->adfBoundsMin, psSubNode->adfBoundsMin,
632 : sizeof(psSubNode->adfBoundsMin));
633 154253 : memcpy(psTreeNode->adfBoundsMax, psSubNode->adfBoundsMax,
634 : sizeof(psSubNode->adfBoundsMax));
635 154253 : psTreeNode->nShapeCount = psSubNode->nShapeCount;
636 154253 : assert(psTreeNode->panShapeIds == SHPLIB_NULLPTR);
637 154253 : psTreeNode->panShapeIds = psSubNode->panShapeIds;
638 154253 : assert(psTreeNode->papsShapeObj == SHPLIB_NULLPTR);
639 154253 : psTreeNode->papsShapeObj = psSubNode->papsShapeObj;
640 154253 : psTreeNode->nSubNodes = psSubNode->nSubNodes;
641 157789 : for (i = 0; i < psSubNode->nSubNodes; i++)
642 3536 : psTreeNode->apsSubNode[i] = psSubNode->apsSubNode[i];
643 154253 : free(psSubNode);
644 : }
645 :
646 : /* -------------------------------------------------------------------- */
647 : /* We should be trimmed if we have no subnodes, and no shapes. */
648 : /* -------------------------------------------------------------------- */
649 740917 : return (psTreeNode->nSubNodes == 0 && psTreeNode->nShapeCount == 0);
650 : }
651 :
652 : /************************************************************************/
653 : /* SHPTreeTrimExtraNodes() */
654 : /* */
655 : /* Trim empty nodes from the tree. Note that we never trim an */
656 : /* empty root node. */
657 : /************************************************************************/
658 :
659 17 : void SHPAPI_CALL SHPTreeTrimExtraNodes(SHPTree *hTree)
660 : {
661 17 : SHPTreeNodeTrim(hTree->psRoot);
662 17 : }
663 :
664 : struct SHPDiskTreeInfo
665 : {
666 : SAHooks sHooks;
667 : SAFile fpQIX;
668 : };
669 :
670 : /************************************************************************/
671 : /* SHPOpenDiskTree() */
672 : /************************************************************************/
673 :
674 1673 : SHPTreeDiskHandle SHPOpenDiskTree(const char *pszQIXFilename,
675 : const SAHooks *psHooks)
676 : {
677 : SHPTreeDiskHandle hDiskTree;
678 :
679 1673 : hDiskTree = STATIC_CAST(SHPTreeDiskHandle,
680 : calloc(1, sizeof(struct SHPDiskTreeInfo)));
681 1673 : if (!hDiskTree)
682 0 : return SHPLIB_NULLPTR;
683 :
684 1673 : if (psHooks == SHPLIB_NULLPTR)
685 1673 : SASetupDefaultHooks(&(hDiskTree->sHooks));
686 : else
687 0 : memcpy(&(hDiskTree->sHooks), psHooks, sizeof(SAHooks));
688 :
689 1673 : hDiskTree->fpQIX = hDiskTree->sHooks.FOpen(pszQIXFilename, "rb",
690 : hDiskTree->sHooks.pvUserData);
691 1673 : if (hDiskTree->fpQIX == SHPLIB_NULLPTR)
692 : {
693 1630 : free(hDiskTree);
694 1630 : return SHPLIB_NULLPTR;
695 : }
696 :
697 43 : return hDiskTree;
698 : }
699 :
700 : /***********************************************************************/
701 : /* SHPCloseDiskTree() */
702 : /************************************************************************/
703 :
704 44 : void SHPCloseDiskTree(SHPTreeDiskHandle hDiskTree)
705 : {
706 44 : if (hDiskTree == SHPLIB_NULLPTR)
707 1 : return;
708 :
709 43 : hDiskTree->sHooks.FClose(hDiskTree->fpQIX);
710 43 : free(hDiskTree);
711 : }
712 :
713 : /************************************************************************/
714 : /* SHPSearchDiskTreeNode() */
715 : /************************************************************************/
716 :
717 525416 : static bool SHPSearchDiskTreeNode(const SHPTreeDiskHandle hDiskTree,
718 : double *padfBoundsMin, double *padfBoundsMax,
719 : int **ppanResultBuffer, int *pnBufferMax,
720 : int *pnResultCount, int bNeedSwap,
721 : int nRecLevel)
722 :
723 : {
724 : unsigned int i;
725 : unsigned int offset;
726 : unsigned int numshapes, numsubnodes;
727 : double adfNodeBoundsMin[2], adfNodeBoundsMax[2];
728 : int nFReadAcc;
729 :
730 : /* -------------------------------------------------------------------- */
731 : /* Read and unswap first part of node info. */
732 : /* -------------------------------------------------------------------- */
733 525416 : nFReadAcc = STATIC_CAST(
734 : int, hDiskTree->sHooks.FRead(&offset, 4, 1, hDiskTree->fpQIX));
735 525416 : if (bNeedSwap)
736 0 : SHP_SWAP32(&offset);
737 :
738 525416 : nFReadAcc += STATIC_CAST(int, hDiskTree->sHooks.FRead(adfNodeBoundsMin,
739 : sizeof(double), 2,
740 : hDiskTree->fpQIX));
741 525416 : nFReadAcc += STATIC_CAST(int, hDiskTree->sHooks.FRead(adfNodeBoundsMax,
742 : sizeof(double), 2,
743 : hDiskTree->fpQIX));
744 525416 : if (bNeedSwap)
745 : {
746 0 : SHP_SWAPDOUBLE(adfNodeBoundsMin + 0);
747 0 : SHP_SWAPDOUBLE(adfNodeBoundsMin + 1);
748 0 : SHP_SWAPDOUBLE(adfNodeBoundsMax + 0);
749 0 : SHP_SWAPDOUBLE(adfNodeBoundsMax + 1);
750 : }
751 :
752 525416 : nFReadAcc += STATIC_CAST(
753 : int, hDiskTree->sHooks.FRead(&numshapes, 4, 1, hDiskTree->fpQIX));
754 525416 : if (bNeedSwap)
755 0 : SHP_SWAP32(&numshapes);
756 :
757 : /* Check that we could read all previous values */
758 525416 : if (nFReadAcc != 1 + 2 + 2 + 1)
759 : {
760 0 : hDiskTree->sHooks.Error("I/O error");
761 0 : return false;
762 : }
763 :
764 : /* Sanity checks to avoid int overflows in later computation */
765 525416 : if (offset > INT_MAX - sizeof(int))
766 : {
767 0 : hDiskTree->sHooks.Error("Invalid value for offset");
768 0 : return false;
769 : }
770 :
771 525416 : if (numshapes > (INT_MAX - offset - sizeof(int)) / sizeof(int) ||
772 525416 : numshapes > INT_MAX / sizeof(int) - *pnResultCount)
773 : {
774 0 : hDiskTree->sHooks.Error("Invalid value for numshapes");
775 0 : return false;
776 : }
777 :
778 : /* -------------------------------------------------------------------- */
779 : /* If we don't overlap this node at all, we can just fseek() */
780 : /* pass this node info and all subnodes. */
781 : /* -------------------------------------------------------------------- */
782 525416 : if (!SHPCheckBoundsOverlap(adfNodeBoundsMin, adfNodeBoundsMax,
783 : padfBoundsMin, padfBoundsMax, 2))
784 : {
785 353895 : offset += numshapes * sizeof(int) + sizeof(int);
786 353895 : hDiskTree->sHooks.FSeek(hDiskTree->fpQIX, offset, SEEK_CUR);
787 353895 : return true;
788 : }
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* Add all the shapeids at this node to our list. */
792 : /* -------------------------------------------------------------------- */
793 171521 : if (numshapes > 0)
794 : {
795 24613 : if (*pnResultCount + numshapes >
796 24613 : STATIC_CAST(unsigned int, *pnBufferMax))
797 : {
798 : int *pNewBuffer;
799 :
800 12421 : *pnBufferMax = (*pnResultCount + numshapes + 100) * 5 / 4;
801 :
802 12421 : if (STATIC_CAST(size_t, *pnBufferMax) > INT_MAX / sizeof(int))
803 0 : *pnBufferMax = *pnResultCount + numshapes;
804 :
805 12421 : pNewBuffer = STATIC_CAST(
806 : int *, realloc(*ppanResultBuffer, *pnBufferMax * sizeof(int)));
807 :
808 12421 : if (pNewBuffer == SHPLIB_NULLPTR)
809 : {
810 0 : hDiskTree->sHooks.Error("Out of memory error");
811 0 : return false;
812 : }
813 :
814 12421 : *ppanResultBuffer = pNewBuffer;
815 : }
816 :
817 24613 : if (hDiskTree->sHooks.FRead(*ppanResultBuffer + *pnResultCount,
818 : sizeof(int), numshapes,
819 24613 : hDiskTree->fpQIX) != numshapes)
820 : {
821 0 : hDiskTree->sHooks.Error("I/O error");
822 0 : return false;
823 : }
824 :
825 24613 : if (bNeedSwap)
826 : {
827 0 : for (i = 0; i < numshapes; i++)
828 0 : SHP_SWAP32(*ppanResultBuffer + *pnResultCount + i);
829 : }
830 :
831 24613 : *pnResultCount += numshapes;
832 : }
833 :
834 : /* -------------------------------------------------------------------- */
835 : /* Process the subnodes. */
836 : /* -------------------------------------------------------------------- */
837 171521 : if (hDiskTree->sHooks.FRead(&numsubnodes, 4, 1, hDiskTree->fpQIX) != 1)
838 : {
839 0 : hDiskTree->sHooks.Error("I/O error");
840 0 : return false;
841 : }
842 171521 : if (bNeedSwap)
843 0 : SHP_SWAP32(&numsubnodes);
844 171521 : if (numsubnodes > 0 && nRecLevel == 32)
845 : {
846 0 : hDiskTree->sHooks.Error("Shape tree is too deep");
847 0 : return false;
848 : }
849 :
850 684613 : for (i = 0; i < numsubnodes; i++)
851 : {
852 513092 : if (!SHPSearchDiskTreeNode(hDiskTree, padfBoundsMin, padfBoundsMax,
853 : ppanResultBuffer, pnBufferMax, pnResultCount,
854 : bNeedSwap, nRecLevel + 1))
855 0 : return false;
856 : }
857 :
858 171521 : return true;
859 : }
860 :
861 : /************************************************************************/
862 : /* SHPTreeReadLibc() */
863 : /************************************************************************/
864 :
865 0 : static SAOffset SHPTreeReadLibc(void *p, SAOffset size, SAOffset nmemb,
866 : SAFile file)
867 :
868 : {
869 0 : return STATIC_CAST(SAOffset, fread(p, STATIC_CAST(size_t, size),
870 : STATIC_CAST(size_t, nmemb),
871 : REINTERPRET_CAST(FILE *, file)));
872 : }
873 :
874 : /************************************************************************/
875 : /* SHPTreeSeekLibc() */
876 : /************************************************************************/
877 :
878 0 : static SAOffset SHPTreeSeekLibc(SAFile file, SAOffset offset, int whence)
879 :
880 : {
881 0 : return STATIC_CAST(SAOffset, fseek(REINTERPRET_CAST(FILE *, file),
882 : STATIC_CAST(long, offset), whence));
883 : }
884 :
885 : /************************************************************************/
886 : /* SHPSearchDiskTree() */
887 : /************************************************************************/
888 :
889 0 : int SHPAPI_CALL1(*) SHPSearchDiskTree(FILE *fp, double *padfBoundsMin,
890 : double *padfBoundsMax, int *pnShapeCount)
891 : {
892 : struct SHPDiskTreeInfo sDiskTree;
893 0 : memset(&sDiskTree.sHooks, 0, sizeof(sDiskTree.sHooks));
894 :
895 : /* We do not use SASetupDefaultHooks() because the FILE* */
896 : /* is a libc FILE* */
897 0 : sDiskTree.sHooks.FSeek = SHPTreeSeekLibc;
898 0 : sDiskTree.sHooks.FRead = SHPTreeReadLibc;
899 :
900 0 : sDiskTree.fpQIX = REINTERPRET_CAST(SAFile, fp);
901 :
902 0 : return SHPSearchDiskTreeEx(&sDiskTree, padfBoundsMin, padfBoundsMax,
903 0 : pnShapeCount);
904 : }
905 :
906 : /***********************************************************************/
907 : /* SHPSearchDiskTreeEx() */
908 : /************************************************************************/
909 :
910 : int SHPAPI_CALL1(*)
911 12324 : SHPSearchDiskTreeEx(const SHPTreeDiskHandle hDiskTree,
912 : double *padfBoundsMin, double *padfBoundsMax,
913 : int *pnShapeCount)
914 :
915 : {
916 12324 : int nBufferMax = 0;
917 : unsigned char abyBuf[16];
918 12324 : int *panResultBuffer = SHPLIB_NULLPTR;
919 :
920 12324 : *pnShapeCount = 0;
921 :
922 : /* -------------------------------------------------------------------- */
923 : /* Read the header. */
924 : /* -------------------------------------------------------------------- */
925 12324 : hDiskTree->sHooks.FSeek(hDiskTree->fpQIX, 0, SEEK_SET);
926 12324 : hDiskTree->sHooks.FRead(abyBuf, 16, 1, hDiskTree->fpQIX);
927 :
928 12324 : if (memcmp(abyBuf, "SQT", 3) != 0)
929 0 : return SHPLIB_NULLPTR;
930 :
931 : #if defined(SHP_BIG_ENDIAN)
932 : bool bNeedSwap = abyBuf[3] != 2;
933 : #else
934 12324 : bool bNeedSwap = abyBuf[3] != 1;
935 : #endif
936 :
937 : /* -------------------------------------------------------------------- */
938 : /* Search through root node and its descendants. */
939 : /* -------------------------------------------------------------------- */
940 12324 : if (!SHPSearchDiskTreeNode(hDiskTree, padfBoundsMin, padfBoundsMax,
941 : &panResultBuffer, &nBufferMax, pnShapeCount,
942 : bNeedSwap, 0))
943 : {
944 0 : if (panResultBuffer != SHPLIB_NULLPTR)
945 0 : free(panResultBuffer);
946 0 : *pnShapeCount = 0;
947 0 : return SHPLIB_NULLPTR;
948 : }
949 : /* -------------------------------------------------------------------- */
950 : /* Sort the id array */
951 : /* -------------------------------------------------------------------- */
952 :
953 : /* To distinguish between empty intersection from error case */
954 12324 : if (panResultBuffer == SHPLIB_NULLPTR)
955 0 : panResultBuffer = STATIC_CAST(int *, calloc(1, sizeof(int)));
956 : else
957 12324 : qsort(panResultBuffer, *pnShapeCount, sizeof(int), SHPTreeCompareInts);
958 :
959 12324 : return panResultBuffer;
960 : }
961 :
962 : /************************************************************************/
963 : /* SHPGetSubNodeOffset() */
964 : /* */
965 : /* Determine how big all the subnodes of this node (and their */
966 : /* children) will be. This will allow disk based searchers to */
967 : /* seek past them all efficiently. */
968 : /************************************************************************/
969 :
970 700492 : static int SHPGetSubNodeOffset(SHPTreeNode *node)
971 : {
972 : int i;
973 700492 : int offset = 0;
974 :
975 1315030 : for (i = 0; i < node->nSubNodes; i++)
976 : {
977 614539 : if (node->apsSubNode[i])
978 : {
979 614539 : offset += 4 * sizeof(double) +
980 614539 : (node->apsSubNode[i]->nShapeCount + 3) * sizeof(int);
981 614539 : offset += SHPGetSubNodeOffset(node->apsSubNode[i]);
982 : }
983 : }
984 :
985 700492 : return (offset);
986 : }
987 :
988 : /************************************************************************/
989 : /* SHPWriteTreeNode() */
990 : /************************************************************************/
991 :
992 85953 : static void SHPWriteTreeNode(SAFile fp, SHPTreeNode *node,
993 : const SAHooks *psHooks)
994 : {
995 : int i, j;
996 : int offset;
997 : unsigned char *pabyRec;
998 85953 : assert(SHPLIB_NULLPTR != node);
999 :
1000 85953 : offset = SHPGetSubNodeOffset(node);
1001 :
1002 85953 : pabyRec = STATIC_CAST(unsigned char *,
1003 : malloc(sizeof(double) * 4 + (3 * sizeof(int)) +
1004 : (node->nShapeCount * sizeof(int))));
1005 85953 : if (SHPLIB_NULLPTR == pabyRec)
1006 : {
1007 : #ifdef USE_CPL
1008 0 : CPLError(CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
1009 : #endif
1010 0 : assert(0);
1011 : return;
1012 : }
1013 :
1014 85953 : memcpy(pabyRec, &offset, 4);
1015 :
1016 : /* minx, miny, maxx, maxy */
1017 85953 : memcpy(pabyRec + 4, node->adfBoundsMin + 0, sizeof(double));
1018 85953 : memcpy(pabyRec + 12, node->adfBoundsMin + 1, sizeof(double));
1019 85953 : memcpy(pabyRec + 20, node->adfBoundsMax + 0, sizeof(double));
1020 85953 : memcpy(pabyRec + 28, node->adfBoundsMax + 1, sizeof(double));
1021 :
1022 85953 : memcpy(pabyRec + 36, &node->nShapeCount, 4);
1023 85953 : j = node->nShapeCount * sizeof(int);
1024 85953 : if (j)
1025 54983 : memcpy(pabyRec + 40, node->panShapeIds, j);
1026 85953 : memcpy(pabyRec + j + 40, &node->nSubNodes, 4);
1027 :
1028 85953 : psHooks->FWrite(pabyRec, 44 + j, 1, fp);
1029 85953 : free(pabyRec);
1030 :
1031 171889 : for (i = 0; i < node->nSubNodes; i++)
1032 : {
1033 85936 : if (node->apsSubNode[i])
1034 85936 : SHPWriteTreeNode(fp, node->apsSubNode[i], psHooks);
1035 : }
1036 : }
1037 :
1038 : /************************************************************************/
1039 : /* SHPWriteTree() */
1040 : /************************************************************************/
1041 :
1042 17 : int SHPAPI_CALL SHPWriteTree(SHPTree *tree, const char *filename)
1043 : {
1044 : SAHooks sHooks;
1045 :
1046 17 : SASetupDefaultHooks(&sHooks);
1047 :
1048 34 : return SHPWriteTreeLL(tree, filename, &sHooks);
1049 : }
1050 :
1051 : /************************************************************************/
1052 : /* SHPWriteTreeLL() */
1053 : /************************************************************************/
1054 :
1055 17 : int SHPWriteTreeLL(SHPTree *tree, const char *filename, const SAHooks *psHooks)
1056 : {
1057 17 : const char signature[4] = "SQT";
1058 : char abyBuf[32];
1059 : SAFile fp;
1060 :
1061 : SAHooks sHooks;
1062 17 : if (psHooks == SHPLIB_NULLPTR)
1063 : {
1064 0 : SASetupDefaultHooks(&sHooks);
1065 0 : psHooks = &sHooks;
1066 : }
1067 :
1068 : /* -------------------------------------------------------------------- */
1069 : /* Open the output file. */
1070 : /* -------------------------------------------------------------------- */
1071 17 : fp = psHooks->FOpen(filename, "wb", psHooks->pvUserData);
1072 17 : if (fp == SHPLIB_NULLPTR)
1073 : {
1074 0 : return FALSE;
1075 : }
1076 :
1077 : /* -------------------------------------------------------------------- */
1078 : /* Write the header. */
1079 : /* -------------------------------------------------------------------- */
1080 17 : memcpy(abyBuf + 0, signature, 3);
1081 :
1082 : #if defined(SHP_BIG_ENDIAN)
1083 : abyBuf[3] = 2; /* New MSB */
1084 : #else
1085 17 : abyBuf[3] = 1; /* New LSB */
1086 : #endif
1087 :
1088 17 : abyBuf[4] = 1; /* version */
1089 17 : abyBuf[5] = 0; /* next 3 reserved */
1090 17 : abyBuf[6] = 0;
1091 17 : abyBuf[7] = 0;
1092 :
1093 17 : psHooks->FWrite(abyBuf, 8, 1, fp);
1094 :
1095 17 : psHooks->FWrite(&(tree->nTotalCount), 4, 1, fp);
1096 :
1097 : /* write maxdepth */
1098 :
1099 17 : psHooks->FWrite(&(tree->nMaxDepth), 4, 1, fp);
1100 :
1101 : /* -------------------------------------------------------------------- */
1102 : /* Write all the nodes "in order". */
1103 : /* -------------------------------------------------------------------- */
1104 :
1105 17 : SHPWriteTreeNode(fp, tree->psRoot, psHooks);
1106 :
1107 17 : psHooks->FClose(fp);
1108 :
1109 17 : return TRUE;
1110 : }
|
