1 : /******************************************************************************
2 : * $Id: shptree.c,v 1.12 2008/11/12 15:39:50 fwarmerdam Exp $
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.12 2008/11/12 15:39:50 fwarmerdam
38 : * improve safety in face of buggy .shp file.
39 : *
40 : * Revision 1.11 2007/10/27 03:31:14 fwarmerdam
41 : * limit default depth of tree to 12 levels (gdal ticket #1594)
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>
79 : #include <stdlib.h>
80 : #include <string.h>
81 : #ifdef USE_CPL
82 : #include <cpl_error.h>
83 : #endif
84 :
85 0 : SHP_CVSID("$Id: shptree.c,v 1.12 2008/11/12 15:39:50 fwarmerdam Exp $")
86 :
87 : #ifndef TRUE
88 : # define TRUE 1
89 : # define FALSE 0
90 : #endif
91 :
92 : static int bBigEndian = 0;
93 :
94 :
95 : /* -------------------------------------------------------------------- */
96 : /* If the following is 0.5, nodes will be split in half. If it */
97 : /* is 0.6 then each subnode will contain 60% of the parent */
98 : /* node, with 20% representing overlap. This can be help to */
99 : /* prevent small objects on a boundary from shifting too high */
100 : /* up the tree. */
101 : /* -------------------------------------------------------------------- */
102 :
103 : #define SHP_SPLIT_RATIO 0.55
104 :
105 : /************************************************************************/
106 : /* SfRealloc() */
107 : /* */
108 : /* A realloc cover function that will access a NULL pointer as */
109 : /* a valid input. */
110 : /************************************************************************/
111 :
112 14665 : static void * SfRealloc( void * pMem, int nNewSize )
113 :
114 : {
115 14665 : if( pMem == NULL )
116 14654 : return( (void *) malloc(nNewSize) );
117 : else
118 11 : return( (void *) realloc(pMem,nNewSize) );
119 : }
120 :
121 : /************************************************************************/
122 : /* SHPTreeNodeInit() */
123 : /* */
124 : /* Initialize a tree node. */
125 : /************************************************************************/
126 :
127 219099 : static SHPTreeNode *SHPTreeNodeCreate( double * padfBoundsMin,
128 : double * padfBoundsMax )
129 :
130 : {
131 : SHPTreeNode *psTreeNode;
132 :
133 219099 : psTreeNode = (SHPTreeNode *) malloc(sizeof(SHPTreeNode));
134 219099 : if( NULL == psTreeNode )
135 : {
136 : #ifdef USE_CPL
137 0 : CPLError( CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
138 : #endif
139 0 : return NULL;
140 : }
141 :
142 219099 : psTreeNode->nShapeCount = 0;
143 219099 : psTreeNode->panShapeIds = NULL;
144 219099 : psTreeNode->papsShapeObj = NULL;
145 :
146 219099 : psTreeNode->nSubNodes = 0;
147 :
148 219099 : if( padfBoundsMin != NULL )
149 219096 : memcpy( psTreeNode->adfBoundsMin, padfBoundsMin, sizeof(double) * 4 );
150 :
151 219099 : if( padfBoundsMax != NULL )
152 219096 : memcpy( psTreeNode->adfBoundsMax, padfBoundsMax, sizeof(double) * 4 );
153 :
154 219099 : return psTreeNode;
155 : }
156 :
157 :
158 : /************************************************************************/
159 : /* SHPCreateTree() */
160 : /************************************************************************/
161 :
162 : SHPTree SHPAPI_CALL1(*)
163 3 : SHPCreateTree( SHPHandle hSHP, int nDimension, int nMaxDepth,
164 : double *padfBoundsMin, double *padfBoundsMax )
165 :
166 : {
167 : SHPTree *psTree;
168 :
169 3 : if( padfBoundsMin == NULL && hSHP == NULL )
170 0 : return NULL;
171 :
172 : /* -------------------------------------------------------------------- */
173 : /* Allocate the tree object */
174 : /* -------------------------------------------------------------------- */
175 3 : psTree = (SHPTree *) malloc(sizeof(SHPTree));
176 3 : if( NULL == psTree )
177 : {
178 : #ifdef USE_CPL
179 0 : CPLError( CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
180 : #endif
181 0 : return NULL;
182 : }
183 :
184 3 : psTree->hSHP = hSHP;
185 3 : psTree->nMaxDepth = nMaxDepth;
186 3 : psTree->nDimension = nDimension;
187 3 : psTree->nTotalCount = 0;
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* If no max depth was defined, try to select a reasonable one */
191 : /* that implies approximately 8 shapes per node. */
192 : /* -------------------------------------------------------------------- */
193 3 : if( psTree->nMaxDepth == 0 && hSHP != NULL )
194 : {
195 3 : int nMaxNodeCount = 1;
196 : int nShapeCount;
197 :
198 3 : SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
199 20 : while( nMaxNodeCount*4 < nShapeCount )
200 : {
201 14 : psTree->nMaxDepth += 1;
202 14 : nMaxNodeCount = nMaxNodeCount * 2;
203 : }
204 :
205 : #ifdef USE_CPL
206 3 : CPLDebug( "Shape",
207 : "Estimated spatial index tree depth: %d",
208 : psTree->nMaxDepth );
209 : #endif
210 :
211 : /* NOTE: Due to problems with memory allocation for deep trees,
212 : * automatically estimated depth is limited up to 12 levels.
213 : * See Ticket #1594 for detailed discussion.
214 : */
215 3 : if( psTree->nMaxDepth > MAX_DEFAULT_TREE_DEPTH )
216 : {
217 0 : psTree->nMaxDepth = MAX_DEFAULT_TREE_DEPTH;
218 :
219 : #ifdef USE_CPL
220 0 : CPLDebug( "Shape",
221 : "Falling back to max number of allowed index tree levels (%d).",
222 : MAX_DEFAULT_TREE_DEPTH );
223 : #endif
224 : }
225 : }
226 :
227 : /* -------------------------------------------------------------------- */
228 : /* Allocate the root node. */
229 : /* -------------------------------------------------------------------- */
230 3 : psTree->psRoot = SHPTreeNodeCreate( padfBoundsMin, padfBoundsMax );
231 3 : if( NULL == psTree->psRoot )
232 : {
233 0 : return NULL;
234 : }
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Assign the bounds to the root node. If none are passed in, */
238 : /* use the bounds of the provided file otherwise the create */
239 : /* function will have already set the bounds. */
240 : /* -------------------------------------------------------------------- */
241 3 : assert( NULL != psTree );
242 3 : assert( NULL != psTree->psRoot );
243 :
244 3 : if( padfBoundsMin == NULL )
245 : {
246 6 : SHPGetInfo( hSHP, NULL, NULL,
247 3 : psTree->psRoot->adfBoundsMin,
248 3 : psTree->psRoot->adfBoundsMax );
249 : }
250 :
251 : /* -------------------------------------------------------------------- */
252 : /* If we have a file, insert all it's shapes into the tree. */
253 : /* -------------------------------------------------------------------- */
254 3 : if( hSHP != NULL )
255 : {
256 : int iShape, nShapeCount;
257 :
258 3 : SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
259 :
260 14659 : for( iShape = 0; iShape < nShapeCount; iShape++ )
261 : {
262 : SHPObject *psShape;
263 :
264 14656 : psShape = SHPReadObject( hSHP, iShape );
265 14656 : if( psShape != NULL )
266 : {
267 14656 : SHPTreeAddShapeId( psTree, psShape );
268 14656 : SHPDestroyObject( psShape );
269 : }
270 : }
271 : }
272 :
273 3 : return psTree;
274 : }
275 :
276 : /************************************************************************/
277 : /* SHPDestroyTreeNode() */
278 : /************************************************************************/
279 :
280 219099 : static void SHPDestroyTreeNode( SHPTreeNode * psTreeNode )
281 :
282 : {
283 : int i;
284 :
285 219099 : assert( NULL != psTreeNode );
286 :
287 288514 : for( i = 0; i < psTreeNode->nSubNodes; i++ )
288 : {
289 69415 : if( psTreeNode->apsSubNode[i] != NULL )
290 69415 : SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
291 : }
292 :
293 219099 : if( psTreeNode->panShapeIds != NULL )
294 14645 : free( psTreeNode->panShapeIds );
295 :
296 219099 : if( psTreeNode->papsShapeObj != NULL )
297 : {
298 0 : for( i = 0; i < psTreeNode->nShapeCount; i++ )
299 : {
300 0 : if( psTreeNode->papsShapeObj[i] != NULL )
301 0 : SHPDestroyObject( psTreeNode->papsShapeObj[i] );
302 : }
303 :
304 0 : free( psTreeNode->papsShapeObj );
305 : }
306 :
307 219099 : free( psTreeNode );
308 219099 : }
309 :
310 : /************************************************************************/
311 : /* SHPDestroyTree() */
312 : /************************************************************************/
313 :
314 : void SHPAPI_CALL
315 3 : SHPDestroyTree( SHPTree * psTree )
316 :
317 : {
318 3 : SHPDestroyTreeNode( psTree->psRoot );
319 3 : free( psTree );
320 3 : }
321 :
322 : /************************************************************************/
323 : /* SHPCheckBoundsOverlap() */
324 : /* */
325 : /* Do the given boxes overlap at all? */
326 : /************************************************************************/
327 :
328 : int SHPAPI_CALL
329 28 : SHPCheckBoundsOverlap( double * padfBox1Min, double * padfBox1Max,
330 : double * padfBox2Min, double * padfBox2Max,
331 : int nDimension )
332 :
333 : {
334 : int iDim;
335 :
336 67 : for( iDim = 0; iDim < nDimension; iDim++ )
337 : {
338 50 : if( padfBox2Max[iDim] < padfBox1Min[iDim] )
339 11 : return FALSE;
340 :
341 39 : if( padfBox1Max[iDim] < padfBox2Min[iDim] )
342 0 : return FALSE;
343 : }
344 :
345 17 : return TRUE;
346 : }
347 :
348 : /************************************************************************/
349 : /* SHPCheckObjectContained() */
350 : /* */
351 : /* Does the given shape fit within the indicated extents? */
352 : /************************************************************************/
353 :
354 621218 : static int SHPCheckObjectContained( SHPObject * psObject, int nDimension,
355 : double * padfBoundsMin, double * padfBoundsMax )
356 :
357 : {
358 1826822 : if( psObject->dfXMin < padfBoundsMin[0]
359 1826822 : || psObject->dfXMax > padfBoundsMax[0] )
360 252118 : return FALSE;
361 :
362 1107300 : if( psObject->dfYMin < padfBoundsMin[1]
363 1107300 : || psObject->dfYMax > padfBoundsMax[1] )
364 153268 : return FALSE;
365 :
366 215832 : if( nDimension == 2 )
367 215832 : return TRUE;
368 :
369 0 : if( psObject->dfZMin < padfBoundsMin[2]
370 0 : || psObject->dfZMax < padfBoundsMax[2] )
371 0 : return FALSE;
372 :
373 0 : if( nDimension == 3 )
374 0 : return TRUE;
375 :
376 0 : if( psObject->dfMMin < padfBoundsMin[3]
377 0 : || psObject->dfMMax < padfBoundsMax[3] )
378 0 : return FALSE;
379 :
380 0 : return TRUE;
381 : }
382 :
383 : /************************************************************************/
384 : /* SHPTreeSplitBounds() */
385 : /* */
386 : /* Split a region into two subregion evenly, cutting along the */
387 : /* longest dimension. */
388 : /************************************************************************/
389 :
390 : void SHPAPI_CALL
391 164325 : SHPTreeSplitBounds( double *padfBoundsMinIn, double *padfBoundsMaxIn,
392 : double *padfBoundsMin1, double * padfBoundsMax1,
393 : double *padfBoundsMin2, double * padfBoundsMax2 )
394 :
395 : {
396 : /* -------------------------------------------------------------------- */
397 : /* The output bounds will be very similar to the input bounds, */
398 : /* so just copy over to start. */
399 : /* -------------------------------------------------------------------- */
400 164325 : memcpy( padfBoundsMin1, padfBoundsMinIn, sizeof(double) * 4 );
401 164325 : memcpy( padfBoundsMax1, padfBoundsMaxIn, sizeof(double) * 4 );
402 164325 : memcpy( padfBoundsMin2, padfBoundsMinIn, sizeof(double) * 4 );
403 164325 : memcpy( padfBoundsMax2, padfBoundsMaxIn, sizeof(double) * 4 );
404 :
405 : /* -------------------------------------------------------------------- */
406 : /* Split in X direction. */
407 : /* -------------------------------------------------------------------- */
408 328650 : if( (padfBoundsMaxIn[0] - padfBoundsMinIn[0])
409 164325 : > (padfBoundsMaxIn[1] - padfBoundsMinIn[1]) )
410 : {
411 109548 : double dfRange = padfBoundsMaxIn[0] - padfBoundsMinIn[0];
412 :
413 109548 : padfBoundsMax1[0] = padfBoundsMinIn[0] + dfRange * SHP_SPLIT_RATIO;
414 109548 : padfBoundsMin2[0] = padfBoundsMaxIn[0] - dfRange * SHP_SPLIT_RATIO;
415 : }
416 :
417 : /* -------------------------------------------------------------------- */
418 : /* Otherwise split in Y direction. */
419 : /* -------------------------------------------------------------------- */
420 : else
421 : {
422 54777 : double dfRange = padfBoundsMaxIn[1] - padfBoundsMinIn[1];
423 :
424 54777 : padfBoundsMax1[1] = padfBoundsMinIn[1] + dfRange * SHP_SPLIT_RATIO;
425 54777 : padfBoundsMin2[1] = padfBoundsMaxIn[1] - dfRange * SHP_SPLIT_RATIO;
426 : }
427 164325 : }
428 :
429 : /************************************************************************/
430 : /* SHPTreeNodeAddShapeId() */
431 : /************************************************************************/
432 :
433 : static int
434 230488 : SHPTreeNodeAddShapeId( SHPTreeNode * psTreeNode, SHPObject * psObject,
435 : int nMaxDepth, int nDimension )
436 :
437 : {
438 : int i;
439 :
440 : /* -------------------------------------------------------------------- */
441 : /* If there are subnodes, then consider wiether this object */
442 : /* will fit in them. */
443 : /* -------------------------------------------------------------------- */
444 230491 : if( nMaxDepth > 1 && psTreeNode->nSubNodes > 0 )
445 : {
446 452571 : for( i = 0; i < psTreeNode->nSubNodes; i++ )
447 : {
448 905136 : if( SHPCheckObjectContained(psObject, nDimension,
449 452568 : psTreeNode->apsSubNode[i]->adfBoundsMin,
450 452568 : psTreeNode->apsSubNode[i]->adfBoundsMax))
451 : {
452 161058 : return SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[i],
453 : psObject, nMaxDepth-1,
454 : nDimension );
455 : }
456 : }
457 : }
458 :
459 : /* -------------------------------------------------------------------- */
460 : /* Otherwise, consider creating four subnodes if could fit into */
461 : /* them, and adding to the appropriate subnode. */
462 : /* -------------------------------------------------------------------- */
463 : #if MAX_SUBNODE == 4
464 69427 : else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
465 : {
466 : double adfBoundsMinH1[4], adfBoundsMaxH1[4];
467 : double adfBoundsMinH2[4], adfBoundsMaxH2[4];
468 : double adfBoundsMin1[4], adfBoundsMax1[4];
469 : double adfBoundsMin2[4], adfBoundsMax2[4];
470 : double adfBoundsMin3[4], adfBoundsMax3[4];
471 : double adfBoundsMin4[4], adfBoundsMax4[4];
472 :
473 54775 : SHPTreeSplitBounds( psTreeNode->adfBoundsMin,
474 : psTreeNode->adfBoundsMax,
475 : adfBoundsMinH1, adfBoundsMaxH1,
476 : adfBoundsMinH2, adfBoundsMaxH2 );
477 :
478 54775 : SHPTreeSplitBounds( adfBoundsMinH1, adfBoundsMaxH1,
479 : adfBoundsMin1, adfBoundsMax1,
480 : adfBoundsMin2, adfBoundsMax2 );
481 :
482 54775 : SHPTreeSplitBounds( adfBoundsMinH2, adfBoundsMaxH2,
483 : adfBoundsMin3, adfBoundsMax3,
484 : adfBoundsMin4, adfBoundsMax4 );
485 :
486 168650 : if( SHPCheckObjectContained(psObject, nDimension,
487 : adfBoundsMin1, adfBoundsMax1)
488 54775 : || SHPCheckObjectContained(psObject, nDimension,
489 : adfBoundsMin2, adfBoundsMax2)
490 49676 : || SHPCheckObjectContained(psObject, nDimension,
491 : adfBoundsMin3, adfBoundsMax3)
492 39873 : || SHPCheckObjectContained(psObject, nDimension,
493 24326 : adfBoundsMin4, adfBoundsMax4) )
494 : {
495 54774 : psTreeNode->nSubNodes = 4;
496 54774 : psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
497 : adfBoundsMax1 );
498 54774 : psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
499 : adfBoundsMax2 );
500 54774 : psTreeNode->apsSubNode[2] = SHPTreeNodeCreate( adfBoundsMin3,
501 : adfBoundsMax3 );
502 54774 : psTreeNode->apsSubNode[3] = SHPTreeNodeCreate( adfBoundsMin4,
503 : adfBoundsMax4 );
504 :
505 : /* recurse back on this node now that it has subnodes */
506 54774 : return( SHPTreeNodeAddShapeId( psTreeNode, psObject,
507 : nMaxDepth, nDimension ) );
508 : }
509 : }
510 : #endif /* MAX_SUBNODE == 4 */
511 :
512 : /* -------------------------------------------------------------------- */
513 : /* Otherwise, consider creating two subnodes if could fit into */
514 : /* them, and adding to the appropriate subnode. */
515 : /* -------------------------------------------------------------------- */
516 : #if MAX_SUBNODE == 2
517 : else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
518 : {
519 : double adfBoundsMin1[4], adfBoundsMax1[4];
520 : double adfBoundsMin2[4], adfBoundsMax2[4];
521 :
522 : SHPTreeSplitBounds( psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
523 : adfBoundsMin1, adfBoundsMax1,
524 : adfBoundsMin2, adfBoundsMax2 );
525 :
526 : if( SHPCheckObjectContained(psObject, nDimension,
527 : adfBoundsMin1, adfBoundsMax1))
528 : {
529 : psTreeNode->nSubNodes = 2;
530 : psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
531 : adfBoundsMax1 );
532 : psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
533 : adfBoundsMax2 );
534 :
535 : return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[0], psObject,
536 : nMaxDepth - 1, nDimension ) );
537 : }
538 : else if( SHPCheckObjectContained(psObject, nDimension,
539 : adfBoundsMin2, adfBoundsMax2) )
540 : {
541 : psTreeNode->nSubNodes = 2;
542 : psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
543 : adfBoundsMax1 );
544 : psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
545 : adfBoundsMax2 );
546 :
547 : return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[1], psObject,
548 : nMaxDepth - 1, nDimension ) );
549 : }
550 : }
551 : #endif /* MAX_SUBNODE == 2 */
552 :
553 : /* -------------------------------------------------------------------- */
554 : /* If none of that worked, just add it to this nodes list. */
555 : /* -------------------------------------------------------------------- */
556 14656 : psTreeNode->nShapeCount++;
557 :
558 29312 : psTreeNode->panShapeIds = (int *)
559 14656 : SfRealloc( psTreeNode->panShapeIds,
560 14656 : sizeof(int) * psTreeNode->nShapeCount );
561 14656 : psTreeNode->panShapeIds[psTreeNode->nShapeCount-1] = psObject->nShapeId;
562 :
563 14656 : if( psTreeNode->papsShapeObj != NULL )
564 : {
565 0 : psTreeNode->papsShapeObj = (SHPObject **)
566 0 : SfRealloc( psTreeNode->papsShapeObj,
567 0 : sizeof(void *) * psTreeNode->nShapeCount );
568 0 : psTreeNode->papsShapeObj[psTreeNode->nShapeCount-1] = NULL;
569 : }
570 :
571 14656 : return TRUE;
572 : }
573 :
574 : /************************************************************************/
575 : /* SHPTreeAddShapeId() */
576 : /* */
577 : /* Add a shape to the tree, but don't keep a pointer to the */
578 : /* object data, just keep the shapeid. */
579 : /************************************************************************/
580 :
581 : int SHPAPI_CALL
582 14656 : SHPTreeAddShapeId( SHPTree * psTree, SHPObject * psObject )
583 :
584 : {
585 14656 : psTree->nTotalCount++;
586 :
587 14656 : return( SHPTreeNodeAddShapeId( psTree->psRoot, psObject,
588 : psTree->nMaxDepth, psTree->nDimension ) );
589 : }
590 :
591 : /************************************************************************/
592 : /* SHPTreeCollectShapesIds() */
593 : /* */
594 : /* Work function implementing SHPTreeFindLikelyShapes() on a */
595 : /* tree node by tree node basis. */
596 : /************************************************************************/
597 :
598 : void SHPAPI_CALL
599 0 : SHPTreeCollectShapeIds( SHPTree *hTree, SHPTreeNode * psTreeNode,
600 : double * padfBoundsMin, double * padfBoundsMax,
601 : int * pnShapeCount, int * pnMaxShapes,
602 : int ** ppanShapeList )
603 :
604 : {
605 : int i;
606 :
607 : /* -------------------------------------------------------------------- */
608 : /* Does this node overlap the area of interest at all? If not, */
609 : /* return without adding to the list at all. */
610 : /* -------------------------------------------------------------------- */
611 0 : if( !SHPCheckBoundsOverlap( psTreeNode->adfBoundsMin,
612 : psTreeNode->adfBoundsMax,
613 : padfBoundsMin,
614 : padfBoundsMax,
615 : hTree->nDimension ) )
616 0 : return;
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* Grow the list to hold the shapes on this node. */
620 : /* -------------------------------------------------------------------- */
621 0 : if( *pnShapeCount + psTreeNode->nShapeCount > *pnMaxShapes )
622 : {
623 0 : *pnMaxShapes = (*pnShapeCount + psTreeNode->nShapeCount) * 2 + 20;
624 0 : *ppanShapeList = (int *)
625 0 : SfRealloc(*ppanShapeList,sizeof(int) * *pnMaxShapes);
626 : }
627 :
628 : /* -------------------------------------------------------------------- */
629 : /* Add the local nodes shapeids to the list. */
630 : /* -------------------------------------------------------------------- */
631 0 : for( i = 0; i < psTreeNode->nShapeCount; i++ )
632 : {
633 0 : (*ppanShapeList)[(*pnShapeCount)++] = psTreeNode->panShapeIds[i];
634 : }
635 :
636 : /* -------------------------------------------------------------------- */
637 : /* Recurse to subnodes if they exist. */
638 : /* -------------------------------------------------------------------- */
639 0 : for( i = 0; i < psTreeNode->nSubNodes; i++ )
640 : {
641 0 : if( psTreeNode->apsSubNode[i] != NULL )
642 0 : SHPTreeCollectShapeIds( hTree, psTreeNode->apsSubNode[i],
643 : padfBoundsMin, padfBoundsMax,
644 : pnShapeCount, pnMaxShapes,
645 : ppanShapeList );
646 : }
647 : }
648 :
649 : /************************************************************************/
650 : /* SHPTreeFindLikelyShapes() */
651 : /* */
652 : /* Find all shapes within tree nodes for which the tree node */
653 : /* bounding box overlaps the search box. The return value is */
654 : /* an array of shapeids terminated by a -1. The shapeids will */
655 : /* be in order, as hopefully this will result in faster (more */
656 : /* sequential) reading from the file. */
657 : /************************************************************************/
658 :
659 : /* helper for qsort */
660 : static int
661 72 : compare_ints( const void * a, const void * b)
662 : {
663 72 : return (*(int*)a) - (*(int*)b);
664 : }
665 :
666 : int SHPAPI_CALL1(*)
667 0 : SHPTreeFindLikelyShapes( SHPTree * hTree,
668 : double * padfBoundsMin, double * padfBoundsMax,
669 : int * pnShapeCount )
670 :
671 : {
672 0 : int *panShapeList=NULL, nMaxShapes = 0;
673 :
674 : /* -------------------------------------------------------------------- */
675 : /* Perform the search by recursive descent. */
676 : /* -------------------------------------------------------------------- */
677 0 : *pnShapeCount = 0;
678 :
679 0 : SHPTreeCollectShapeIds( hTree, hTree->psRoot,
680 : padfBoundsMin, padfBoundsMax,
681 : pnShapeCount, &nMaxShapes,
682 : &panShapeList );
683 :
684 : /* -------------------------------------------------------------------- */
685 : /* Sort the id array */
686 : /* -------------------------------------------------------------------- */
687 :
688 0 : qsort(panShapeList, *pnShapeCount, sizeof(int), compare_ints);
689 :
690 0 : return panShapeList;
691 : }
692 :
693 : /************************************************************************/
694 : /* SHPTreeNodeTrim() */
695 : /* */
696 : /* This is the recurve version of SHPTreeTrimExtraNodes() that */
697 : /* walks the tree cleaning it up. */
698 : /************************************************************************/
699 :
700 219099 : static int SHPTreeNodeTrim( SHPTreeNode * psTreeNode )
701 :
702 : {
703 : int i;
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* Trim subtrees, and free subnodes that come back empty. */
707 : /* -------------------------------------------------------------------- */
708 438195 : for( i = 0; i < psTreeNode->nSubNodes; i++ )
709 : {
710 219096 : if( SHPTreeNodeTrim( psTreeNode->apsSubNode[i] ) )
711 : {
712 149681 : SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
713 :
714 299362 : psTreeNode->apsSubNode[i] =
715 149681 : psTreeNode->apsSubNode[psTreeNode->nSubNodes-1];
716 :
717 149681 : psTreeNode->nSubNodes--;
718 :
719 149681 : i--; /* process the new occupant of this subnode entry */
720 : }
721 : }
722 :
723 : /* -------------------------------------------------------------------- */
724 : /* We should be trimmed if we have no subnodes, and no shapes. */
725 : /* -------------------------------------------------------------------- */
726 219099 : return( psTreeNode->nSubNodes == 0 && psTreeNode->nShapeCount == 0 );
727 : }
728 :
729 : /************************************************************************/
730 : /* SHPTreeTrimExtraNodes() */
731 : /* */
732 : /* Trim empty nodes from the tree. Note that we never trim an */
733 : /* empty root node. */
734 : /************************************************************************/
735 :
736 : void SHPAPI_CALL
737 3 : SHPTreeTrimExtraNodes( SHPTree * hTree )
738 :
739 : {
740 3 : SHPTreeNodeTrim( hTree->psRoot );
741 3 : }
742 :
743 : /************************************************************************/
744 : /* SwapWord() */
745 : /* */
746 : /* Swap a 2, 4 or 8 byte word. */
747 : /************************************************************************/
748 :
749 0 : static void SwapWord( int length, void * wordP )
750 :
751 : {
752 : int i;
753 : unsigned char temp;
754 :
755 0 : for( i=0; i < length/2; i++ )
756 : {
757 0 : temp = ((unsigned char *) wordP)[i];
758 0 : ((unsigned char *)wordP)[i] = ((unsigned char *) wordP)[length-i-1];
759 0 : ((unsigned char *) wordP)[length-i-1] = temp;
760 : }
761 0 : }
762 :
763 : /************************************************************************/
764 : /* SHPSearchDiskTreeNode() */
765 : /************************************************************************/
766 :
767 : static int
768 28 : SHPSearchDiskTreeNode( FILE *fp, double *padfBoundsMin, double *padfBoundsMax,
769 : int **ppanResultBuffer, int *pnBufferMax,
770 : int *pnResultCount, int bNeedSwap )
771 :
772 : {
773 : int i;
774 : int offset;
775 : int numshapes, numsubnodes;
776 : double adfNodeBoundsMin[2], adfNodeBoundsMax[2];
777 :
778 : /* -------------------------------------------------------------------- */
779 : /* Read and unswap first part of node info. */
780 : /* -------------------------------------------------------------------- */
781 28 : fread( &offset, 4, 1, fp );
782 28 : if ( bNeedSwap ) SwapWord ( 4, &offset );
783 :
784 28 : fread( adfNodeBoundsMin, sizeof(double), 2, fp );
785 28 : fread( adfNodeBoundsMax, sizeof(double), 2, fp );
786 28 : if ( bNeedSwap )
787 : {
788 0 : SwapWord( 8, adfNodeBoundsMin + 0 );
789 0 : SwapWord( 8, adfNodeBoundsMin + 1 );
790 0 : SwapWord( 8, adfNodeBoundsMax + 0 );
791 0 : SwapWord( 8, adfNodeBoundsMax + 1 );
792 : }
793 :
794 28 : fread( &numshapes, 4, 1, fp );
795 28 : if ( bNeedSwap ) SwapWord ( 4, &numshapes );
796 :
797 : /* -------------------------------------------------------------------- */
798 : /* If we don't overlap this node at all, we can just fseek() */
799 : /* pass this node info and all subnodes. */
800 : /* -------------------------------------------------------------------- */
801 28 : if( !SHPCheckBoundsOverlap( adfNodeBoundsMin, adfNodeBoundsMax,
802 : padfBoundsMin, padfBoundsMax, 2 ) )
803 : {
804 11 : offset += numshapes*sizeof(int) + sizeof(int);
805 11 : fseek(fp, offset, SEEK_CUR);
806 11 : return TRUE;
807 : }
808 :
809 : /* -------------------------------------------------------------------- */
810 : /* Add all the shapeids at this node to our list. */
811 : /* -------------------------------------------------------------------- */
812 17 : if(numshapes > 0)
813 : {
814 17 : if( *pnResultCount + numshapes > *pnBufferMax )
815 : {
816 9 : *pnBufferMax = (int) ((*pnResultCount + numshapes + 100) * 1.25);
817 18 : *ppanResultBuffer = (int *)
818 18 : SfRealloc( *ppanResultBuffer, *pnBufferMax * sizeof(int) );
819 : }
820 :
821 17 : fread( *ppanResultBuffer + *pnResultCount,
822 : sizeof(int), numshapes, fp );
823 :
824 17 : if (bNeedSwap )
825 : {
826 0 : for( i=0; i<numshapes; i++ )
827 0 : SwapWord( 4, *ppanResultBuffer + *pnResultCount + i );
828 : }
829 :
830 17 : *pnResultCount += numshapes;
831 : }
832 :
833 : /* -------------------------------------------------------------------- */
834 : /* Process the subnodes. */
835 : /* -------------------------------------------------------------------- */
836 17 : fread( &numsubnodes, 4, 1, fp );
837 17 : if ( bNeedSwap ) SwapWord ( 4, &numsubnodes );
838 :
839 35 : for(i=0; i<numsubnodes; i++)
840 : {
841 18 : if( !SHPSearchDiskTreeNode( fp, padfBoundsMin, padfBoundsMax,
842 : ppanResultBuffer, pnBufferMax,
843 : pnResultCount, bNeedSwap ) )
844 0 : return FALSE;
845 : }
846 :
847 17 : return TRUE;
848 : }
849 :
850 : /************************************************************************/
851 : /* SHPSearchDiskTree() */
852 : /************************************************************************/
853 :
854 : int SHPAPI_CALL1(*)
855 10 : SHPSearchDiskTree( FILE *fp,
856 : double *padfBoundsMin, double *padfBoundsMax,
857 : int *pnShapeCount )
858 :
859 : {
860 10 : int i, bNeedSwap, nBufferMax = 0;
861 : unsigned char abyBuf[16];
862 10 : int *panResultBuffer = NULL;
863 :
864 10 : *pnShapeCount = 0;
865 :
866 : /* -------------------------------------------------------------------- */
867 : /* Establish the byte order on this machine. */
868 : /* -------------------------------------------------------------------- */
869 10 : i = 1;
870 10 : if( *((unsigned char *) &i) == 1 )
871 10 : bBigEndian = FALSE;
872 : else
873 0 : bBigEndian = TRUE;
874 :
875 : /* -------------------------------------------------------------------- */
876 : /* Read the header. */
877 : /* -------------------------------------------------------------------- */
878 10 : fseek( fp, 0, SEEK_SET );
879 10 : fread( abyBuf, 16, 1, fp );
880 :
881 10 : if( memcmp( abyBuf, "SQT", 3 ) != 0 )
882 0 : return NULL;
883 :
884 40 : if( (abyBuf[3] == 2 && bBigEndian)
885 20 : || (abyBuf[3] == 1 && !bBigEndian) )
886 10 : bNeedSwap = FALSE;
887 : else
888 0 : bNeedSwap = TRUE;
889 :
890 : /* -------------------------------------------------------------------- */
891 : /* Search through root node and it's decendents. */
892 : /* -------------------------------------------------------------------- */
893 10 : if( !SHPSearchDiskTreeNode( fp, padfBoundsMin, padfBoundsMax,
894 : &panResultBuffer, &nBufferMax,
895 : pnShapeCount, bNeedSwap ) )
896 : {
897 0 : if( panResultBuffer != NULL )
898 0 : free( panResultBuffer );
899 0 : *pnShapeCount = 0;
900 0 : return NULL;
901 : }
902 : /* -------------------------------------------------------------------- */
903 : /* Sort the id array */
904 : /* -------------------------------------------------------------------- */
905 10 : qsort(panResultBuffer, *pnShapeCount, sizeof(int), compare_ints);
906 :
907 10 : return panResultBuffer;
908 : }
909 :
910 : /************************************************************************/
911 : /* SHPGetSubNodeOffset() */
912 : /* */
913 : /* Determine how big all the subnodes of this node (and their */
914 : /* children) will be. This will allow disk based searchers to */
915 : /* seek past them all efficiently. */
916 : /************************************************************************/
917 :
918 696276 : static int SHPGetSubNodeOffset( SHPTreeNode *node)
919 : {
920 : int i;
921 696276 : long offset=0;
922 :
923 1323134 : for(i=0; i<node->nSubNodes; i++ )
924 : {
925 626858 : if(node->apsSubNode[i])
926 : {
927 626858 : offset += 4*sizeof(double)
928 626858 : + (node->apsSubNode[i]->nShapeCount+3)*sizeof(int);
929 626858 : offset += SHPGetSubNodeOffset(node->apsSubNode[i]);
930 : }
931 : }
932 :
933 696276 : return(offset);
934 : }
935 :
936 : /************************************************************************/
937 : /* SHPWriteTreeNode() */
938 : /************************************************************************/
939 :
940 69418 : static void SHPWriteTreeNode( FILE *fp, SHPTreeNode *node)
941 : {
942 : int i,j;
943 : int offset;
944 69418 : unsigned char *pabyRec = NULL;
945 69418 : assert( NULL != node );
946 :
947 69418 : offset = SHPGetSubNodeOffset(node);
948 :
949 69418 : pabyRec = (unsigned char *)
950 : malloc(sizeof(double) * 4
951 : + (3 * sizeof(int)) + (node->nShapeCount * sizeof(int)) );
952 69418 : if( NULL == pabyRec )
953 : {
954 : #ifdef USE_CPL
955 0 : CPLError( CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
956 : #endif
957 0 : assert( 0 );
958 : }
959 69418 : assert( NULL != pabyRec );
960 :
961 69418 : memcpy( pabyRec, &offset, 4);
962 :
963 : /* minx, miny, maxx, maxy */
964 69418 : memcpy( pabyRec+ 4, node->adfBoundsMin+0, sizeof(double) );
965 69418 : memcpy( pabyRec+12, node->adfBoundsMin+1, sizeof(double) );
966 69418 : memcpy( pabyRec+20, node->adfBoundsMax+0, sizeof(double) );
967 69418 : memcpy( pabyRec+28, node->adfBoundsMax+1, sizeof(double) );
968 :
969 69418 : memcpy( pabyRec+36, &node->nShapeCount, 4);
970 69418 : j = node->nShapeCount * sizeof(int);
971 69418 : memcpy( pabyRec+40, node->panShapeIds, j);
972 69418 : memcpy( pabyRec+j+40, &node->nSubNodes, 4);
973 :
974 69418 : fwrite( pabyRec, 44+j, 1, fp );
975 69418 : free (pabyRec);
976 :
977 138833 : for(i=0; i<node->nSubNodes; i++ )
978 : {
979 69415 : if(node->apsSubNode[i])
980 69415 : SHPWriteTreeNode( fp, node->apsSubNode[i]);
981 : }
982 69418 : }
983 :
984 : /************************************************************************/
985 : /* SHPWriteTree() */
986 : /************************************************************************/
987 :
988 3 : int SHPWriteTree(SHPTree *tree, const char *filename )
989 : {
990 3 : char signature[4] = "SQT";
991 : int i;
992 : char abyBuf[32];
993 : FILE *fp;
994 :
995 : /* -------------------------------------------------------------------- */
996 : /* Open the output file. */
997 : /* -------------------------------------------------------------------- */
998 3 : fp = fopen(filename, "wb");
999 3 : if( fp == NULL )
1000 : {
1001 0 : return FALSE;
1002 : }
1003 :
1004 : /* -------------------------------------------------------------------- */
1005 : /* Establish the byte order on this machine. */
1006 : /* -------------------------------------------------------------------- */
1007 3 : i = 1;
1008 3 : if( *((unsigned char *) &i) == 1 )
1009 3 : bBigEndian = FALSE;
1010 : else
1011 0 : bBigEndian = TRUE;
1012 :
1013 : /* -------------------------------------------------------------------- */
1014 : /* Write the header. */
1015 : /* -------------------------------------------------------------------- */
1016 3 : memcpy( abyBuf+0, signature, 3 );
1017 :
1018 3 : if( bBigEndian )
1019 0 : abyBuf[3] = 2; /* New MSB */
1020 : else
1021 3 : abyBuf[3] = 1; /* New LSB */
1022 :
1023 3 : abyBuf[4] = 1; /* version */
1024 3 : abyBuf[5] = 0; /* next 3 reserved */
1025 3 : abyBuf[6] = 0;
1026 3 : abyBuf[7] = 0;
1027 :
1028 3 : fwrite( abyBuf, 8, 1, fp );
1029 :
1030 3 : fwrite( &(tree->nTotalCount), 4, 1, fp );
1031 :
1032 : /* write maxdepth */
1033 :
1034 3 : fwrite( &(tree->nMaxDepth), 4, 1, fp );
1035 :
1036 : /* -------------------------------------------------------------------- */
1037 : /* Write all the nodes "in order". */
1038 : /* -------------------------------------------------------------------- */
1039 :
1040 3 : SHPWriteTreeNode( fp, tree->psRoot );
1041 :
1042 3 : fclose( fp );
1043 :
1044 3 : return TRUE;
1045 : }
|
