1 /******************************************************************************
2 * $Id: shptree.c 361 2019-10-03 12:20:24Z rsbivand $
3 *
4 * Project: Shapelib
5 * Purpose: Implementation of quadtree building and searching functions.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 1999, Frank Warmerdam
10 *
11 * This software is available under the following "MIT Style" license,
12 * or at the option of the licensee under the LGPL (see LICENSE.LGPL). This
13 * option is discussed in more detail in shapelib.html.
14 *
15 * --
16 *
17 * Permission is hereby granted, free of charge, to any person obtaining a
18 * copy of this software and associated documentation files (the "Software"),
19 * to deal in the Software without restriction, including without limitation
20 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
21 * and/or sell copies of the Software, and to permit persons to whom the
22 * Software is furnished to do so, subject to the following conditions:
23 *
24 * The above copyright notice and this permission notice shall be included
25 * in all copies or substantial portions of the Software.
26 *
27 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
28 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
30 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
32 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
33 * DEALINGS IN THE SOFTWARE.
34 ******************************************************************************
35 *
36 * $Log: shptree.c,v $
37 * Revision 1.2 2007/11/10 13:17:42 rsbivand
38 * assert
39 *
40 * Revision 1.1.1.1 2005/09/01 18:22:22 rsbivand
41 * Initial import.
42 *
43 * Revision 1.10 2005/01/03 22:30:13 fwarmerdam
44 * added support for saved quadtrees
45 *
46 * Revision 1.9 2003/01/28 15:53:41 warmerda
47 * Avoid build warnings.
48 *
49 * Revision 1.8 2002/05/07 13:07:45 warmerda
50 * use qsort() - patch from Bernhard Herzog
51 *
52 * Revision 1.7 2002/01/15 14:36:07 warmerda
53 * updated email address
54 *
55 * Revision 1.6 2001/05/23 13:36:52 warmerda
56 * added use of SHPAPI_CALL
57 *
58 * Revision 1.5 1999/11/05 14:12:05 warmerda
59 * updated license terms
60 *
61 * Revision 1.4 1999/06/02 18:24:21 warmerda
62 * added trimming code
63 *
64 * Revision 1.3 1999/06/02 17:56:12 warmerda
65 * added quad'' subnode support for trees
66 *
67 * Revision 1.2 1999/05/18 19:11:11 warmerda
68 * Added example searching capability
69 *
70 * Revision 1.1 1999/05/18 17:49:20 warmerda
71 * New
72 *
73 */
74
75 #include "shapefil.h"
76
77 #include <math.h>
78 /* #include <assert.h> RSB 071110 */
79 #include <stdlib.h>
80 #include <string.h>
81
82 //SHP_CVSID("$Id: shptree.c 361 2019-10-03 12:20:24Z rsbivand $")
83
84 #ifndef TRUE
85 # define TRUE 1
86 # define FALSE 0
87 #endif
88
89 static int bBigEndian = 0;
90
91
92 /* -------------------------------------------------------------------- */
93 /* If the following is 0.5, nodes will be split in half. If it */
94 /* is 0.6 then each subnode will contain 60% of the parent */
95 /* node, with 20% representing overlap. This can be help to */
96 /* prevent small objects on a boundary from shifting too high */
97 /* up the tree. */
98 /* -------------------------------------------------------------------- */
99
100 #define SHP_SPLIT_RATIO 0.55
101
102 /************************************************************************/
103 /* SfRealloc() */
104 /* */
105 /* A realloc cover function that will access a NULL pointer as */
106 /* a valid input. */
107 /************************************************************************/
108
SfRealloc(void * pMem,size_t nNewSize)109 static void * SfRealloc( void * pMem, size_t nNewSize )
110
111 {
112 if( pMem == NULL )
113 return( (void *) malloc((size_t) nNewSize) );
114 else
115 return( (void *) realloc(pMem, (size_t) nNewSize) );
116 }
117
118 /************************************************************************/
119 /* SHPTreeNodeInit() */
120 /* */
121 /* Initialize a tree node. */
122 /************************************************************************/
123
SHPTreeNodeCreate(double * padfBoundsMin,double * padfBoundsMax)124 static SHPTreeNode *SHPTreeNodeCreate( double * padfBoundsMin,
125 double * padfBoundsMax )
126
127 {
128 SHPTreeNode *psTreeNode;
129
130 psTreeNode = (SHPTreeNode *) malloc(sizeof(SHPTreeNode));
131
132 psTreeNode->nShapeCount = 0;
133 psTreeNode->panShapeIds = NULL;
134 psTreeNode->papsShapeObj = NULL;
135
136 psTreeNode->nSubNodes = 0;
137
138 if( padfBoundsMin != NULL )
139 memcpy( psTreeNode->adfBoundsMin, padfBoundsMin, sizeof(double) * 4 );
140
141 if( padfBoundsMax != NULL )
142 memcpy( psTreeNode->adfBoundsMax, padfBoundsMax, sizeof(double) * 4 );
143
144 return psTreeNode;
145 }
146
147
148 /************************************************************************/
149 /* SHPCreateTree() */
150 /************************************************************************/
151
152 SHPTree SHPAPI_CALL1(*)
SHPCreateTree(SHPHandle hSHP,int nDimension,int nMaxDepth,double * padfBoundsMin,double * padfBoundsMax)153 SHPCreateTree( SHPHandle hSHP, int nDimension, int nMaxDepth,
154 double *padfBoundsMin, double *padfBoundsMax )
155
156 {
157 SHPTree *psTree;
158
159 if( padfBoundsMin == NULL && hSHP == NULL )
160 return NULL;
161
162 /* -------------------------------------------------------------------- */
163 /* Allocate the tree object */
164 /* -------------------------------------------------------------------- */
165 psTree = (SHPTree *) malloc(sizeof(SHPTree));
166
167 psTree->hSHP = hSHP;
168 psTree->nMaxDepth = nMaxDepth;
169 psTree->nDimension = nDimension;
170 psTree->nTotalCount = 0;
171
172 /* -------------------------------------------------------------------- */
173 /* If no max depth was defined, try to select a reasonable one */
174 /* that implies approximately 8 shapes per node. */
175 /* -------------------------------------------------------------------- */
176 if( psTree->nMaxDepth == 0 && hSHP != NULL )
177 {
178 int nMaxNodeCount = 1;
179 int nShapeCount;
180
181 SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
182 while( nMaxNodeCount*4 < nShapeCount )
183 {
184 psTree->nMaxDepth += 1;
185 nMaxNodeCount = nMaxNodeCount * 2;
186 }
187 }
188
189 /* -------------------------------------------------------------------- */
190 /* Allocate the root node. */
191 /* -------------------------------------------------------------------- */
192 psTree->psRoot = SHPTreeNodeCreate( padfBoundsMin, padfBoundsMax );
193
194 /* -------------------------------------------------------------------- */
195 /* Assign the bounds to the root node. If none are passed in, */
196 /* use the bounds of the provided file otherwise the create */
197 /* function will have already set the bounds. */
198 /* -------------------------------------------------------------------- */
199 if( padfBoundsMin == NULL )
200 {
201 SHPGetInfo( hSHP, NULL, NULL,
202 psTree->psRoot->adfBoundsMin,
203 psTree->psRoot->adfBoundsMax );
204 }
205
206 /* -------------------------------------------------------------------- */
207 /* If we have a file, insert all it's shapes into the tree. */
208 /* -------------------------------------------------------------------- */
209 if( hSHP != NULL )
210 {
211 int iShape, nShapeCount;
212
213 SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
214
215 for( iShape = 0; iShape < nShapeCount; iShape++ )
216 {
217 SHPObject *psShape;
218
219 psShape = SHPReadObject( hSHP, iShape );
220 SHPTreeAddShapeId( psTree, psShape );
221 SHPDestroyObject( psShape );
222 }
223 }
224
225 return psTree;
226 }
227
228 /************************************************************************/
229 /* SHPDestroyTreeNode() */
230 /************************************************************************/
231
SHPDestroyTreeNode(SHPTreeNode * psTreeNode)232 static void SHPDestroyTreeNode( SHPTreeNode * psTreeNode )
233
234 {
235 int i;
236
237 for( i = 0; i < psTreeNode->nSubNodes; i++ )
238 {
239 if( psTreeNode->apsSubNode[i] != NULL )
240 SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
241 }
242
243 if( psTreeNode->panShapeIds != NULL )
244 free( psTreeNode->panShapeIds );
245
246 if( psTreeNode->papsShapeObj != NULL )
247 {
248 for( i = 0; i < psTreeNode->nShapeCount; i++ )
249 {
250 if( psTreeNode->papsShapeObj[i] != NULL )
251 SHPDestroyObject( psTreeNode->papsShapeObj[i] );
252 }
253
254 free( psTreeNode->papsShapeObj );
255 }
256
257 free( psTreeNode );
258 }
259
260 /************************************************************************/
261 /* SHPDestroyTree() */
262 /************************************************************************/
263
264 void SHPAPI_CALL
SHPDestroyTree(SHPTree * psTree)265 SHPDestroyTree( SHPTree * psTree )
266
267 {
268 SHPDestroyTreeNode( psTree->psRoot );
269 free( psTree );
270 }
271
272 /************************************************************************/
273 /* SHPCheckBoundsOverlap() */
274 /* */
275 /* Do the given boxes overlap at all? */
276 /************************************************************************/
277
278 int SHPAPI_CALL
SHPCheckBoundsOverlap(double * padfBox1Min,double * padfBox1Max,double * padfBox2Min,double * padfBox2Max,int nDimension)279 SHPCheckBoundsOverlap( double * padfBox1Min, double * padfBox1Max,
280 double * padfBox2Min, double * padfBox2Max,
281 int nDimension )
282
283 {
284 int iDim;
285
286 for( iDim = 0; iDim < nDimension; iDim++ )
287 {
288 if( padfBox2Max[iDim] < padfBox1Min[iDim] )
289 return FALSE;
290
291 if( padfBox1Max[iDim] < padfBox2Min[iDim] )
292 return FALSE;
293 }
294
295 return TRUE;
296 }
297
298 /************************************************************************/
299 /* SHPCheckObjectContained() */
300 /* */
301 /* Does the given shape fit within the indicated extents? */
302 /************************************************************************/
303
SHPCheckObjectContained(SHPObject * psObject,int nDimension,double * padfBoundsMin,double * padfBoundsMax)304 static int SHPCheckObjectContained( SHPObject * psObject, int nDimension,
305 double * padfBoundsMin, double * padfBoundsMax )
306
307 {
308 if( psObject->dfXMin < padfBoundsMin[0]
309 || psObject->dfXMax > padfBoundsMax[0] )
310 return FALSE;
311
312 if( psObject->dfYMin < padfBoundsMin[1]
313 || psObject->dfYMax > padfBoundsMax[1] )
314 return FALSE;
315
316 if( nDimension == 2 )
317 return TRUE;
318
319 if( psObject->dfZMin < padfBoundsMin[2]
320 || psObject->dfZMax < padfBoundsMax[2] )
321 return FALSE;
322
323 if( nDimension == 3 )
324 return TRUE;
325
326 if( psObject->dfMMin < padfBoundsMin[3]
327 || psObject->dfMMax < padfBoundsMax[3] )
328 return FALSE;
329
330 return TRUE;
331 }
332
333 /************************************************************************/
334 /* SHPTreeSplitBounds() */
335 /* */
336 /* Split a region into two subregion evenly, cutting along the */
337 /* longest dimension. */
338 /************************************************************************/
339
340 void SHPAPI_CALL
SHPTreeSplitBounds(double * padfBoundsMinIn,double * padfBoundsMaxIn,double * padfBoundsMin1,double * padfBoundsMax1,double * padfBoundsMin2,double * padfBoundsMax2)341 SHPTreeSplitBounds( double *padfBoundsMinIn, double *padfBoundsMaxIn,
342 double *padfBoundsMin1, double * padfBoundsMax1,
343 double *padfBoundsMin2, double * padfBoundsMax2 )
344
345 {
346 /* -------------------------------------------------------------------- */
347 /* The output bounds will be very similar to the input bounds, */
348 /* so just copy over to start. */
349 /* -------------------------------------------------------------------- */
350 memcpy( padfBoundsMin1, padfBoundsMinIn, sizeof(double) * 4 );
351 memcpy( padfBoundsMax1, padfBoundsMaxIn, sizeof(double) * 4 );
352 memcpy( padfBoundsMin2, padfBoundsMinIn, sizeof(double) * 4 );
353 memcpy( padfBoundsMax2, padfBoundsMaxIn, sizeof(double) * 4 );
354
355 /* -------------------------------------------------------------------- */
356 /* Split in X direction. */
357 /* -------------------------------------------------------------------- */
358 if( (padfBoundsMaxIn[0] - padfBoundsMinIn[0])
359 > (padfBoundsMaxIn[1] - padfBoundsMinIn[1]) )
360 {
361 double dfRange = padfBoundsMaxIn[0] - padfBoundsMinIn[0];
362
363 padfBoundsMax1[0] = padfBoundsMinIn[0] + dfRange * SHP_SPLIT_RATIO;
364 padfBoundsMin2[0] = padfBoundsMaxIn[0] - dfRange * SHP_SPLIT_RATIO;
365 }
366
367 /* -------------------------------------------------------------------- */
368 /* Otherwise split in Y direction. */
369 /* -------------------------------------------------------------------- */
370 else
371 {
372 double dfRange = padfBoundsMaxIn[1] - padfBoundsMinIn[1];
373
374 padfBoundsMax1[1] = padfBoundsMinIn[1] + dfRange * SHP_SPLIT_RATIO;
375 padfBoundsMin2[1] = padfBoundsMaxIn[1] - dfRange * SHP_SPLIT_RATIO;
376 }
377 }
378
379 /************************************************************************/
380 /* SHPTreeNodeAddShapeId() */
381 /************************************************************************/
382
383 static int
SHPTreeNodeAddShapeId(SHPTreeNode * psTreeNode,SHPObject * psObject,int nMaxDepth,int nDimension)384 SHPTreeNodeAddShapeId( SHPTreeNode * psTreeNode, SHPObject * psObject,
385 int nMaxDepth, int nDimension )
386
387 {
388 int i;
389
390 /* -------------------------------------------------------------------- */
391 /* If there are subnodes, then consider wiether this object */
392 /* will fit in them. */
393 /* -------------------------------------------------------------------- */
394 if( nMaxDepth > 1 && psTreeNode->nSubNodes > 0 )
395 {
396 for( i = 0; i < psTreeNode->nSubNodes; i++ )
397 {
398 if( SHPCheckObjectContained(psObject, nDimension,
399 psTreeNode->apsSubNode[i]->adfBoundsMin,
400 psTreeNode->apsSubNode[i]->adfBoundsMax))
401 {
402 return SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[i],
403 psObject, nMaxDepth-1,
404 nDimension );
405 }
406 }
407 }
408
409 /* -------------------------------------------------------------------- */
410 /* Otherwise, consider creating four subnodes if could fit into */
411 /* them, and adding to the appropriate subnode. */
412 /* -------------------------------------------------------------------- */
413 #if MAX_SUBNODE == 4
414 else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
415 {
416 double adfBoundsMinH1[4], adfBoundsMaxH1[4];
417 double adfBoundsMinH2[4], adfBoundsMaxH2[4];
418 double adfBoundsMin1[4], adfBoundsMax1[4];
419 double adfBoundsMin2[4], adfBoundsMax2[4];
420 double adfBoundsMin3[4], adfBoundsMax3[4];
421 double adfBoundsMin4[4], adfBoundsMax4[4];
422
423 SHPTreeSplitBounds( psTreeNode->adfBoundsMin,
424 psTreeNode->adfBoundsMax,
425 adfBoundsMinH1, adfBoundsMaxH1,
426 adfBoundsMinH2, adfBoundsMaxH2 );
427
428 SHPTreeSplitBounds( adfBoundsMinH1, adfBoundsMaxH1,
429 adfBoundsMin1, adfBoundsMax1,
430 adfBoundsMin2, adfBoundsMax2 );
431
432 SHPTreeSplitBounds( adfBoundsMinH2, adfBoundsMaxH2,
433 adfBoundsMin3, adfBoundsMax3,
434 adfBoundsMin4, adfBoundsMax4 );
435
436 if( SHPCheckObjectContained(psObject, nDimension,
437 adfBoundsMin1, adfBoundsMax1)
438 || SHPCheckObjectContained(psObject, nDimension,
439 adfBoundsMin2, adfBoundsMax2)
440 || SHPCheckObjectContained(psObject, nDimension,
441 adfBoundsMin3, adfBoundsMax3)
442 || SHPCheckObjectContained(psObject, nDimension,
443 adfBoundsMin4, adfBoundsMax4) )
444 {
445 psTreeNode->nSubNodes = 4;
446 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
447 adfBoundsMax1 );
448 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
449 adfBoundsMax2 );
450 psTreeNode->apsSubNode[2] = SHPTreeNodeCreate( adfBoundsMin3,
451 adfBoundsMax3 );
452 psTreeNode->apsSubNode[3] = SHPTreeNodeCreate( adfBoundsMin4,
453 adfBoundsMax4 );
454
455 /* recurse back on this node now that it has subnodes */
456 return( SHPTreeNodeAddShapeId( psTreeNode, psObject,
457 nMaxDepth, nDimension ) );
458 }
459 }
460 #endif /* MAX_SUBNODE == 4 */
461
462 /* -------------------------------------------------------------------- */
463 /* Otherwise, consider creating two subnodes if could fit into */
464 /* them, and adding to the appropriate subnode. */
465 /* -------------------------------------------------------------------- */
466 #if MAX_SUBNODE == 2
467 else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
468 {
469 double adfBoundsMin1[4], adfBoundsMax1[4];
470 double adfBoundsMin2[4], adfBoundsMax2[4];
471
472 SHPTreeSplitBounds( psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
473 adfBoundsMin1, adfBoundsMax1,
474 adfBoundsMin2, adfBoundsMax2 );
475
476 if( SHPCheckObjectContained(psObject, nDimension,
477 adfBoundsMin1, adfBoundsMax1))
478 {
479 psTreeNode->nSubNodes = 2;
480 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
481 adfBoundsMax1 );
482 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
483 adfBoundsMax2 );
484
485 return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[0], psObject,
486 nMaxDepth - 1, nDimension ) );
487 }
488 else if( SHPCheckObjectContained(psObject, nDimension,
489 adfBoundsMin2, adfBoundsMax2) )
490 {
491 psTreeNode->nSubNodes = 2;
492 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
493 adfBoundsMax1 );
494 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
495 adfBoundsMax2 );
496
497 return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[1], psObject,
498 nMaxDepth - 1, nDimension ) );
499 }
500 }
501 #endif /* MAX_SUBNODE == 2 */
502
503 /* -------------------------------------------------------------------- */
504 /* If none of that worked, just add it to this nodes list. */
505 /* -------------------------------------------------------------------- */
506 psTreeNode->nShapeCount++;
507
508 psTreeNode->panShapeIds = (int *)
509 SfRealloc( psTreeNode->panShapeIds,
510 sizeof(int) * (size_t) psTreeNode->nShapeCount );
511 psTreeNode->panShapeIds[psTreeNode->nShapeCount-1] = psObject->nShapeId;
512
513 if( psTreeNode->papsShapeObj != NULL )
514 {
515 psTreeNode->papsShapeObj = (SHPObject **)
516 SfRealloc( psTreeNode->papsShapeObj,
517 sizeof(void *) * (size_t) psTreeNode->nShapeCount );
518 psTreeNode->papsShapeObj[psTreeNode->nShapeCount-1] = NULL;
519 }
520
521 return TRUE;
522 }
523
524 /************************************************************************/
525 /* SHPTreeAddShapeId() */
526 /* */
527 /* Add a shape to the tree, but don't keep a pointer to the */
528 /* object data, just keep the shapeid. */
529 /************************************************************************/
530
531 int SHPAPI_CALL
SHPTreeAddShapeId(SHPTree * psTree,SHPObject * psObject)532 SHPTreeAddShapeId( SHPTree * psTree, SHPObject * psObject )
533
534 {
535 psTree->nTotalCount++;
536
537 return( SHPTreeNodeAddShapeId( psTree->psRoot, psObject,
538 psTree->nMaxDepth, psTree->nDimension ) );
539 }
540
541 /************************************************************************/
542 /* SHPTreeCollectShapesIds() */
543 /* */
544 /* Work function implementing SHPTreeFindLikelyShapes() on a */
545 /* tree node by tree node basis. */
546 /************************************************************************/
547
548 void SHPAPI_CALL
SHPTreeCollectShapeIds(SHPTree * hTree,SHPTreeNode * psTreeNode,double * padfBoundsMin,double * padfBoundsMax,int * pnShapeCount,int * pnMaxShapes,int ** ppanShapeList)549 SHPTreeCollectShapeIds( SHPTree *hTree, SHPTreeNode * psTreeNode,
550 double * padfBoundsMin, double * padfBoundsMax,
551 int * pnShapeCount, int * pnMaxShapes,
552 int ** ppanShapeList )
553
554 {
555 int i;
556
557 /* -------------------------------------------------------------------- */
558 /* Does this node overlap the area of interest at all? If not, */
559 /* return without adding to the list at all. */
560 /* -------------------------------------------------------------------- */
561 if( !SHPCheckBoundsOverlap( psTreeNode->adfBoundsMin,
562 psTreeNode->adfBoundsMax,
563 padfBoundsMin,
564 padfBoundsMax,
565 hTree->nDimension ) )
566 return;
567
568 /* -------------------------------------------------------------------- */
569 /* Grow the list to hold the shapes on this node. */
570 /* -------------------------------------------------------------------- */
571 if( *pnShapeCount + psTreeNode->nShapeCount > *pnMaxShapes )
572 {
573 *pnMaxShapes = (*pnShapeCount + psTreeNode->nShapeCount) * 2 + 20;
574 *ppanShapeList = (int *)
575 SfRealloc(*ppanShapeList,sizeof(int) * (size_t) *pnMaxShapes);
576 }
577
578 /* -------------------------------------------------------------------- */
579 /* Add the local nodes shapeids to the list. */
580 /* -------------------------------------------------------------------- */
581 for( i = 0; i < psTreeNode->nShapeCount; i++ )
582 {
583 (*ppanShapeList)[(*pnShapeCount)++] = psTreeNode->panShapeIds[i];
584 }
585
586 /* -------------------------------------------------------------------- */
587 /* Recurse to subnodes if they exist. */
588 /* -------------------------------------------------------------------- */
589 for( i = 0; i < psTreeNode->nSubNodes; i++ )
590 {
591 if( psTreeNode->apsSubNode[i] != NULL )
592 SHPTreeCollectShapeIds( hTree, psTreeNode->apsSubNode[i],
593 padfBoundsMin, padfBoundsMax,
594 pnShapeCount, pnMaxShapes,
595 ppanShapeList );
596 }
597 }
598
599 /************************************************************************/
600 /* SHPTreeFindLikelyShapes() */
601 /* */
602 /* Find all shapes within tree nodes for which the tree node */
603 /* bounding box overlaps the search box. The return value is */
604 /* an array of shapeids terminated by a -1. The shapeids will */
605 /* be in order, as hopefully this will result in faster (more */
606 /* sequential) reading from the file. */
607 /************************************************************************/
608
609 /* helper for qsort */
610 static int
compare_ints(const void * a,const void * b)611 compare_ints( const void * a, const void * b)
612 {
613 return (*(int*)a) - (*(int*)b);
614 }
615
616 int SHPAPI_CALL1(*)
SHPTreeFindLikelyShapes(SHPTree * hTree,double * padfBoundsMin,double * padfBoundsMax,int * pnShapeCount)617 SHPTreeFindLikelyShapes( SHPTree * hTree,
618 double * padfBoundsMin, double * padfBoundsMax,
619 int * pnShapeCount )
620
621 {
622 int *panShapeList=NULL, nMaxShapes = 0;
623
624 /* -------------------------------------------------------------------- */
625 /* Perform the search by recursive descent. */
626 /* -------------------------------------------------------------------- */
627 *pnShapeCount = 0;
628
629 SHPTreeCollectShapeIds( hTree, hTree->psRoot,
630 padfBoundsMin, padfBoundsMax,
631 pnShapeCount, &nMaxShapes,
632 &panShapeList );
633
634 /* -------------------------------------------------------------------- */
635 /* Sort the id array */
636 /* -------------------------------------------------------------------- */
637
638 qsort(panShapeList, (size_t) *pnShapeCount, sizeof(int), compare_ints);
639
640 return panShapeList;
641 }
642
643 /************************************************************************/
644 /* SHPTreeNodeTrim() */
645 /* */
646 /* This is the recurve version of SHPTreeTrimExtraNodes() that */
647 /* walks the tree cleaning it up. */
648 /************************************************************************/
649
SHPTreeNodeTrim(SHPTreeNode * psTreeNode)650 static int SHPTreeNodeTrim( SHPTreeNode * psTreeNode )
651
652 {
653 int i;
654
655 /* -------------------------------------------------------------------- */
656 /* Trim subtrees, and free subnodes that come back empty. */
657 /* -------------------------------------------------------------------- */
658 for( i = 0; i < psTreeNode->nSubNodes; i++ )
659 {
660 if( SHPTreeNodeTrim( psTreeNode->apsSubNode[i] ) )
661 {
662 SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
663
664 psTreeNode->apsSubNode[i] =
665 psTreeNode->apsSubNode[psTreeNode->nSubNodes-1];
666
667 psTreeNode->nSubNodes--;
668
669 i--; /* process the new occupant of this subnode entry */
670 }
671 }
672
673 /* -------------------------------------------------------------------- */
674 /* We should be trimmed if we have no subnodes, and no shapes. */
675 /* -------------------------------------------------------------------- */
676 return( psTreeNode->nSubNodes == 0 && psTreeNode->nShapeCount == 0 );
677 }
678
679 /************************************************************************/
680 /* SHPTreeTrimExtraNodes() */
681 /* */
682 /* Trim empty nodes from the tree. Note that we never trim an */
683 /* empty root node. */
684 /************************************************************************/
685
686 void SHPAPI_CALL
SHPTreeTrimExtraNodes(SHPTree * hTree)687 SHPTreeTrimExtraNodes( SHPTree * hTree )
688
689 {
690 SHPTreeNodeTrim( hTree->psRoot );
691 }
692
693 /************************************************************************/
694 /* SwapWord() */
695 /* */
696 /* Swap a 2, 4 or 8 byte word. */
697 /************************************************************************/
698
SwapWord(int length,void * wordP)699 static void SwapWord( int length, void * wordP )
700
701 {
702 int i;
703 unsigned char temp;
704
705 for( i=0; i < length/2; i++ )
706 {
707 temp = ((unsigned char *) wordP)[i];
708 ((unsigned char *)wordP)[i] = ((unsigned char *) wordP)[length-i-1];
709 ((unsigned char *) wordP)[length-i-1] = temp;
710 }
711 }
712
713 /************************************************************************/
714 /* SHPSearchDiskTreeNode() */
715 /************************************************************************/
716
717 static int
SHPSearchDiskTreeNode(FILE * fp,double * padfBoundsMin,double * padfBoundsMax,int ** ppanResultBuffer,int * pnBufferMax,int * pnResultCount,int bNeedSwap)718 SHPSearchDiskTreeNode( FILE *fp, double *padfBoundsMin, double *padfBoundsMax,
719 int **ppanResultBuffer, int *pnBufferMax,
720 int *pnResultCount, int bNeedSwap )
721
722 {
723 int i;
724 size_t offset, out;
725 int numshapes, numsubnodes;
726 double adfNodeBoundsMin[2], adfNodeBoundsMax[2];
727
728 /* -------------------------------------------------------------------- */
729 /* Read and unswap first part of node info. */
730 /* -------------------------------------------------------------------- */
731 out = fread( &offset, 4, 1, fp );
732 if ( bNeedSwap ) SwapWord ( 4, &offset );
733
734 out = fread( adfNodeBoundsMin, sizeof(double), 2, fp );
735 out = fread( adfNodeBoundsMax, sizeof(double), 2, fp );
736 if ( bNeedSwap )
737 {
738 SwapWord( 8, adfNodeBoundsMin + 0 );
739 SwapWord( 8, adfNodeBoundsMin + 1 );
740 SwapWord( 8, adfNodeBoundsMax + 0 );
741 SwapWord( 8, adfNodeBoundsMax + 1 );
742 }
743
744 out = fread( &numshapes, 4, 1, fp );
745 if ( bNeedSwap ) SwapWord ( 4, &numshapes );
746
747 /* -------------------------------------------------------------------- */
748 /* If we don't overlap this node at all, we can just fseek() */
749 /* pass this node info and all subnodes. */
750 /* -------------------------------------------------------------------- */
751 if( !SHPCheckBoundsOverlap( adfNodeBoundsMin, adfNodeBoundsMax,
752 padfBoundsMin, padfBoundsMax, 2 ) )
753 {
754 offset += (size_t) numshapes*sizeof(int) + sizeof(int);
755 fseek(fp, (long) offset, SEEK_CUR);
756 return TRUE;
757 }
758
759 /* -------------------------------------------------------------------- */
760 /* Add all the shapeids at this node to our list. */
761 /* -------------------------------------------------------------------- */
762 if(numshapes > 0)
763 {
764 if( *pnResultCount + numshapes > *pnBufferMax )
765 {
766 *pnBufferMax = (int) ((*pnResultCount + numshapes + 100) * 1.25);
767 *ppanResultBuffer = (int *)
768 SfRealloc( *ppanResultBuffer, (size_t) *pnBufferMax * sizeof(int) );
769 }
770
771 out = fread( *ppanResultBuffer + *pnResultCount,
772 sizeof(int), (size_t) numshapes, fp );
773
774 if (bNeedSwap )
775 {
776 for( i=0; i<numshapes; i++ )
777 SwapWord( 4, *ppanResultBuffer + *pnResultCount + i );
778 }
779
780 *pnResultCount += numshapes;
781 }
782
783 /* -------------------------------------------------------------------- */
784 /* Process the subnodes. */
785 /* -------------------------------------------------------------------- */
786 out = fread( &numsubnodes, 4, 1, fp );
787 if ( bNeedSwap ) SwapWord ( 4, &numsubnodes );
788
789 for(i=0; i<numsubnodes; i++)
790 {
791 if( !SHPSearchDiskTreeNode( fp, padfBoundsMin, padfBoundsMax,
792 ppanResultBuffer, pnBufferMax,
793 pnResultCount, bNeedSwap ) )
794 return FALSE;
795 }
796
797 return TRUE;
798 }
799
800 /************************************************************************/
801 /* SHPSearchDiskTree() */
802 /************************************************************************/
803
804 int SHPAPI_CALL1(*)
SHPSearchDiskTree(FILE * fp,double * padfBoundsMin,double * padfBoundsMax,int * pnShapeCount)805 SHPSearchDiskTree( FILE *fp,
806 double *padfBoundsMin, double *padfBoundsMax,
807 int *pnShapeCount )
808
809 {
810 int i, bNeedSwap, nBufferMax = 0;
811 unsigned char abyBuf[16];
812 int *panResultBuffer = NULL;
813 size_t out;
814
815 *pnShapeCount = 0;
816
817 /* -------------------------------------------------------------------- */
818 /* Establish the byte order on this machine. */
819 /* -------------------------------------------------------------------- */
820 i = 1;
821 if( *((unsigned char *) &i) == 1 )
822 bBigEndian = FALSE;
823 else
824 bBigEndian = TRUE;
825
826 /* -------------------------------------------------------------------- */
827 /* Read the header. */
828 /* -------------------------------------------------------------------- */
829 fseek( fp, 0, SEEK_SET );
830 out = fread( abyBuf, 16, 1, fp );
831
832 if( memcmp( abyBuf, "SQT", 3 ) != 0 )
833 return NULL;
834
835 if( (abyBuf[3] == 2 && bBigEndian)
836 || (abyBuf[3] == 1 && !bBigEndian) )
837 bNeedSwap = FALSE;
838 else
839 bNeedSwap = TRUE;
840
841 /* -------------------------------------------------------------------- */
842 /* Search through root node and it's decendents. */
843 /* -------------------------------------------------------------------- */
844 if( !SHPSearchDiskTreeNode( fp, padfBoundsMin, padfBoundsMax,
845 &panResultBuffer, &nBufferMax,
846 pnShapeCount, bNeedSwap ) )
847 {
848 if( panResultBuffer != NULL )
849 free( panResultBuffer );
850 *pnShapeCount = 0;
851 return NULL;
852 }
853 /* -------------------------------------------------------------------- */
854 /* Sort the id array */
855 /* -------------------------------------------------------------------- */
856 qsort(panResultBuffer, (size_t) *pnShapeCount, sizeof(int), compare_ints);
857
858 return panResultBuffer;
859 }
860
861 /************************************************************************/
862 /* SHPGetSubNodeOffset() */
863 /* */
864 /* Determine how big all the subnodes of this node (and their */
865 /* children) will be. This will allow disk based searchers to */
866 /* seek past them all efficiently. */
867 /************************************************************************/
868
SHPGetSubNodeOffset(SHPTreeNode * node)869 static int SHPGetSubNodeOffset( SHPTreeNode *node)
870 {
871 int i;
872 size_t offset=0;
873
874 for(i=0; i<node->nSubNodes; i++ )
875 {
876 if(node->apsSubNode[i])
877 {
878 offset += 4*sizeof(double)
879 + (size_t) (node->apsSubNode[i]->nShapeCount+3)*sizeof(int);
880 offset += (size_t) SHPGetSubNodeOffset(node->apsSubNode[i]);
881 }
882 }
883
884 return((int) offset);
885 }
886
887 /************************************************************************/
888 /* SHPWriteTreeNode() */
889 /************************************************************************/
890
SHPWriteTreeNode(FILE * fp,SHPTreeNode * node)891 static void SHPWriteTreeNode( FILE *fp, SHPTreeNode *node)
892 {
893 int i,j;
894 int offset;
895 unsigned char *pabyRec = NULL;
896
897 offset = SHPGetSubNodeOffset(node);
898
899 pabyRec = (unsigned char *)
900 malloc(sizeof(double) * 4
901 + (3 * sizeof(int)) + ((size_t) node->nShapeCount * sizeof(int)) );
902
903 memcpy( pabyRec, &offset, 4);
904
905 /* minx, miny, maxx, maxy */
906 memcpy( pabyRec+ 4, node->adfBoundsMin+0, sizeof(double) );
907 memcpy( pabyRec+12, node->adfBoundsMin+1, sizeof(double) );
908 memcpy( pabyRec+20, node->adfBoundsMax+0, sizeof(double) );
909 memcpy( pabyRec+28, node->adfBoundsMax+1, sizeof(double) );
910
911 memcpy( pabyRec+36, &node->nShapeCount, 4);
912 j = node->nShapeCount * (int) sizeof(int);
913 memcpy( pabyRec+40, node->panShapeIds, (size_t) j);
914 memcpy( pabyRec+j+40, &node->nSubNodes, 4);
915
916 fwrite( pabyRec, (size_t) (44+j), 1, fp );
917 free (pabyRec);
918
919 for(i=0; i<node->nSubNodes; i++ )
920 {
921 if(node->apsSubNode[i])
922 SHPWriteTreeNode( fp, node->apsSubNode[i]);
923 }
924 }
925
926 /************************************************************************/
927 /* SHPWriteTree() */
928 /************************************************************************/
929
SHPWriteTree(SHPTree * tree,const char * filename)930 int SHPWriteTree(SHPTree *tree, const char *filename )
931 {
932 char signature[4] = "SQT";
933 int i;
934 char abyBuf[32];
935 FILE *fp;
936
937 /* -------------------------------------------------------------------- */
938 /* Open the output file. */
939 /* -------------------------------------------------------------------- */
940 fp = fopen(filename, "wb");
941 if( fp == NULL )
942 {
943 return FALSE;
944 }
945
946 /* -------------------------------------------------------------------- */
947 /* Establish the byte order on this machine. */
948 /* -------------------------------------------------------------------- */
949 i = 1;
950 if( *((unsigned char *) &i) == 1 )
951 bBigEndian = FALSE;
952 else
953 bBigEndian = TRUE;
954
955 /* -------------------------------------------------------------------- */
956 /* Write the header. */
957 /* -------------------------------------------------------------------- */
958 memcpy( abyBuf+0, signature, 3 );
959
960 if( bBigEndian )
961 abyBuf[3] = 2; /* New MSB */
962 else
963 abyBuf[3] = 1; /* New LSB */
964
965 abyBuf[4] = 1; /* version */
966 abyBuf[5] = 0; /* next 3 reserved */
967 abyBuf[6] = 0;
968 abyBuf[7] = 0;
969
970 fwrite( abyBuf, 8, 1, fp );
971
972 fwrite( &(tree->nTotalCount), 4, 1, fp );
973
974 /* write maxdepth */
975
976 fwrite( &(tree->nMaxDepth), 4, 1, fp );
977
978 /* -------------------------------------------------------------------- */
979 /* Write all the nodes "in order". */
980 /* -------------------------------------------------------------------- */
981
982 SHPWriteTreeNode( fp, tree->psRoot );
983
984 fclose( fp );
985
986 return TRUE;
987 }
988