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