1 /******************************************************************************
2 * $Id: shptree.c,v 1.17 2012-01-27 21:09:26 fwarmerdam Exp $
5 * Purpose: Implementation of quadtree building and searching functions.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
8 ******************************************************************************
9 * Copyright (c) 1999, Frank Warmerdam
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.
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:
24 * The above copyright notice and this permission notice shall be included
25 * in all copies or substantial portions of the Software.
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 ******************************************************************************
37 * Revision 1.17 2012-01-27 21:09:26 fwarmerdam
38 * optimize .qix output (gdal #4472)
40 * Revision 1.16 2011-12-11 22:26:46 fwarmerdam
41 * upgrade .qix access code to use SAHooks (gdal #3365)
43 * Revision 1.15 2011-07-24 05:59:25 fwarmerdam
44 * minimize use of CPLError in favor of SAHooks.Error()
46 * Revision 1.14 2010-08-27 23:43:27 fwarmerdam
47 * add SHPAPI_CALL attribute in code
49 * Revision 1.13 2010-06-29 05:50:15 fwarmerdam
50 * fix sign of Z/M comparisons in SHPCheckObjectContained (#2223)
52 * Revision 1.12 2008-11-12 15:39:50 fwarmerdam
53 * improve safety in face of buggy .shp file.
55 * Revision 1.11 2007/10/27 03:31:14 fwarmerdam
56 * limit default depth of tree to 12 levels (gdal ticket #1594)
58 * Revision 1.10 2005/01/03 22:30:13 fwarmerdam
59 * added support for saved quadtrees
61 * Revision 1.9 2003/01/28 15:53:41 warmerda
62 * Avoid build warnings.
64 * Revision 1.8 2002/05/07 13:07:45 warmerda
65 * use qsort() - patch from Bernhard Herzog
67 * Revision 1.7 2002/01/15 14:36:07 warmerda
68 * updated email address
70 * Revision 1.6 2001/05/23 13:36:52 warmerda
71 * added use of SHPAPI_CALL
73 * Revision 1.5 1999/11/05 14:12:05 warmerda
74 * updated license terms
76 * Revision 1.4 1999/06/02 18:24:21 warmerda
79 * Revision 1.3 1999/06/02 17:56:12 warmerda
80 * added quad'' subnode support for trees
82 * Revision 1.2 1999/05/18 19:11:11 warmerda
83 * Added example searching capability
85 * Revision 1.1 1999/05/18 17:49:20 warmerda
98 #include "cpl_error.h"
101 SHP_CVSID("$Id: shptree.c,v 1.17 2012-01-27 21:09:26 fwarmerdam Exp $")
108 static int bBigEndian = 0;
111 /* -------------------------------------------------------------------- */
112 /* If the following is 0.5, nodes will be split in half. If it */
113 /* is 0.6 then each subnode will contain 60% of the parent */
114 /* node, with 20% representing overlap. This can be help to */
115 /* prevent small objects on a boundary from shifting too high */
117 /* -------------------------------------------------------------------- */
119 #define SHP_SPLIT_RATIO 0.55
121 /************************************************************************/
124 /* A realloc cover function that will access a NULL pointer as */
126 /************************************************************************/
128 static void * SfRealloc( void * pMem, int nNewSize )
132 return( (void *) malloc(nNewSize) );
134 return( (void *) realloc(pMem,nNewSize) );
137 /************************************************************************/
138 /* SHPTreeNodeInit() */
140 /* Initialize a tree node. */
141 /************************************************************************/
143 static SHPTreeNode *SHPTreeNodeCreate( double * padfBoundsMin,
144 double * padfBoundsMax )
147 SHPTreeNode *psTreeNode;
149 psTreeNode = (SHPTreeNode *) malloc(sizeof(SHPTreeNode));
150 if( NULL == psTreeNode )
153 psTreeNode->nShapeCount = 0;
154 psTreeNode->panShapeIds = NULL;
155 psTreeNode->papsShapeObj = NULL;
157 psTreeNode->nSubNodes = 0;
159 if( padfBoundsMin != NULL )
160 memcpy( psTreeNode->adfBoundsMin, padfBoundsMin, sizeof(double) * 4 );
162 if( padfBoundsMax != NULL )
163 memcpy( psTreeNode->adfBoundsMax, padfBoundsMax, sizeof(double) * 4 );
169 /************************************************************************/
170 /* SHPCreateTree() */
171 /************************************************************************/
173 SHPTree SHPAPI_CALL1(*)
174 SHPCreateTree( SHPHandle hSHP, int nDimension, int nMaxDepth,
175 double *padfBoundsMin, double *padfBoundsMax )
180 if( padfBoundsMin == NULL && hSHP == NULL )
183 /* -------------------------------------------------------------------- */
184 /* Allocate the tree object */
185 /* -------------------------------------------------------------------- */
186 psTree = (SHPTree *) malloc(sizeof(SHPTree));
193 psTree->nMaxDepth = nMaxDepth;
194 psTree->nDimension = nDimension;
195 psTree->nTotalCount = 0;
197 /* -------------------------------------------------------------------- */
198 /* If no max depth was defined, try to select a reasonable one */
199 /* that implies approximately 8 shapes per node. */
200 /* -------------------------------------------------------------------- */
201 if( psTree->nMaxDepth == 0 && hSHP != NULL )
203 int nMaxNodeCount = 1;
206 SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
207 while( nMaxNodeCount*4 < nShapeCount )
209 psTree->nMaxDepth += 1;
210 nMaxNodeCount = nMaxNodeCount * 2;
215 "Estimated spatial index tree depth: %d",
219 /* NOTE: Due to problems with memory allocation for deep trees,
220 * automatically estimated depth is limited up to 12 levels.
221 * See Ticket #1594 for detailed discussion.
223 if( psTree->nMaxDepth > MAX_DEFAULT_TREE_DEPTH )
225 psTree->nMaxDepth = MAX_DEFAULT_TREE_DEPTH;
229 "Falling back to max number of allowed index tree levels (%d).",
230 MAX_DEFAULT_TREE_DEPTH );
235 /* -------------------------------------------------------------------- */
236 /* Allocate the root node. */
237 /* -------------------------------------------------------------------- */
238 psTree->psRoot = SHPTreeNodeCreate( padfBoundsMin, padfBoundsMax );
239 if( NULL == psTree->psRoot )
244 /* -------------------------------------------------------------------- */
245 /* Assign the bounds to the root node. If none are passed in, */
246 /* use the bounds of the provided file otherwise the create */
247 /* function will have already set the bounds. */
248 /* -------------------------------------------------------------------- */
249 assert( NULL != psTree );
250 assert( NULL != psTree->psRoot );
252 if( padfBoundsMin == NULL )
254 SHPGetInfo( hSHP, NULL, NULL,
255 psTree->psRoot->adfBoundsMin,
256 psTree->psRoot->adfBoundsMax );
259 /* -------------------------------------------------------------------- */
260 /* If we have a file, insert all it's shapes into the tree. */
261 /* -------------------------------------------------------------------- */
264 int iShape, nShapeCount;
266 SHPGetInfo( hSHP, &nShapeCount, NULL, NULL, NULL );
268 for( iShape = 0; iShape < nShapeCount; iShape++ )
272 psShape = SHPReadObject( hSHP, iShape );
273 if( psShape != NULL )
275 SHPTreeAddShapeId( psTree, psShape );
276 SHPDestroyObject( psShape );
284 /************************************************************************/
285 /* SHPDestroyTreeNode() */
286 /************************************************************************/
288 static void SHPDestroyTreeNode( SHPTreeNode * psTreeNode )
293 assert( NULL != psTreeNode );
295 for( i = 0; i < psTreeNode->nSubNodes; i++ )
297 if( psTreeNode->apsSubNode[i] != NULL )
298 SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
301 if( psTreeNode->panShapeIds != NULL )
302 free( psTreeNode->panShapeIds );
304 if( psTreeNode->papsShapeObj != NULL )
306 for( i = 0; i < psTreeNode->nShapeCount; i++ )
308 if( psTreeNode->papsShapeObj[i] != NULL )
309 SHPDestroyObject( psTreeNode->papsShapeObj[i] );
312 free( psTreeNode->papsShapeObj );
318 /************************************************************************/
319 /* SHPDestroyTree() */
320 /************************************************************************/
323 SHPDestroyTree( SHPTree * psTree )
326 SHPDestroyTreeNode( psTree->psRoot );
330 /************************************************************************/
331 /* SHPCheckBoundsOverlap() */
333 /* Do the given boxes overlap at all? */
334 /************************************************************************/
337 SHPCheckBoundsOverlap( double * padfBox1Min, double * padfBox1Max,
338 double * padfBox2Min, double * padfBox2Max,
344 for( iDim = 0; iDim < nDimension; iDim++ )
346 if( padfBox2Max[iDim] < padfBox1Min[iDim] )
349 if( padfBox1Max[iDim] < padfBox2Min[iDim] )
356 /************************************************************************/
357 /* SHPCheckObjectContained() */
359 /* Does the given shape fit within the indicated extents? */
360 /************************************************************************/
362 static int SHPCheckObjectContained( SHPObject * psObject, int nDimension,
363 double * padfBoundsMin, double * padfBoundsMax )
366 if( psObject->dfXMin < padfBoundsMin[0]
367 || psObject->dfXMax > padfBoundsMax[0] )
370 if( psObject->dfYMin < padfBoundsMin[1]
371 || psObject->dfYMax > padfBoundsMax[1] )
374 if( nDimension == 2 )
377 if( psObject->dfZMin < padfBoundsMin[2]
378 || psObject->dfZMax > padfBoundsMax[2] )
381 if( nDimension == 3 )
384 if( psObject->dfMMin < padfBoundsMin[3]
385 || psObject->dfMMax > padfBoundsMax[3] )
391 /************************************************************************/
392 /* SHPTreeSplitBounds() */
394 /* Split a region into two subregion evenly, cutting along the */
395 /* longest dimension. */
396 /************************************************************************/
399 SHPTreeSplitBounds( double *padfBoundsMinIn, double *padfBoundsMaxIn,
400 double *padfBoundsMin1, double * padfBoundsMax1,
401 double *padfBoundsMin2, double * padfBoundsMax2 )
404 /* -------------------------------------------------------------------- */
405 /* The output bounds will be very similar to the input bounds, */
406 /* so just copy over to start. */
407 /* -------------------------------------------------------------------- */
408 memcpy( padfBoundsMin1, padfBoundsMinIn, sizeof(double) * 4 );
409 memcpy( padfBoundsMax1, padfBoundsMaxIn, sizeof(double) * 4 );
410 memcpy( padfBoundsMin2, padfBoundsMinIn, sizeof(double) * 4 );
411 memcpy( padfBoundsMax2, padfBoundsMaxIn, sizeof(double) * 4 );
413 /* -------------------------------------------------------------------- */
414 /* Split in X direction. */
415 /* -------------------------------------------------------------------- */
416 if( (padfBoundsMaxIn[0] - padfBoundsMinIn[0])
417 > (padfBoundsMaxIn[1] - padfBoundsMinIn[1]) )
419 double dfRange = padfBoundsMaxIn[0] - padfBoundsMinIn[0];
421 padfBoundsMax1[0] = padfBoundsMinIn[0] + dfRange * SHP_SPLIT_RATIO;
422 padfBoundsMin2[0] = padfBoundsMaxIn[0] - dfRange * SHP_SPLIT_RATIO;
425 /* -------------------------------------------------------------------- */
426 /* Otherwise split in Y direction. */
427 /* -------------------------------------------------------------------- */
430 double dfRange = padfBoundsMaxIn[1] - padfBoundsMinIn[1];
432 padfBoundsMax1[1] = padfBoundsMinIn[1] + dfRange * SHP_SPLIT_RATIO;
433 padfBoundsMin2[1] = padfBoundsMaxIn[1] - dfRange * SHP_SPLIT_RATIO;
437 /************************************************************************/
438 /* SHPTreeNodeAddShapeId() */
439 /************************************************************************/
442 SHPTreeNodeAddShapeId( SHPTreeNode * psTreeNode, SHPObject * psObject,
443 int nMaxDepth, int nDimension )
448 /* -------------------------------------------------------------------- */
449 /* If there are subnodes, then consider wiether this object */
450 /* will fit in them. */
451 /* -------------------------------------------------------------------- */
452 if( nMaxDepth > 1 && psTreeNode->nSubNodes > 0 )
454 for( i = 0; i < psTreeNode->nSubNodes; i++ )
456 if( SHPCheckObjectContained(psObject, nDimension,
457 psTreeNode->apsSubNode[i]->adfBoundsMin,
458 psTreeNode->apsSubNode[i]->adfBoundsMax))
460 return SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[i],
461 psObject, nMaxDepth-1,
467 /* -------------------------------------------------------------------- */
468 /* Otherwise, consider creating four subnodes if could fit into */
469 /* them, and adding to the appropriate subnode. */
470 /* -------------------------------------------------------------------- */
472 else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
474 double adfBoundsMinH1[4], adfBoundsMaxH1[4];
475 double adfBoundsMinH2[4], adfBoundsMaxH2[4];
476 double adfBoundsMin1[4], adfBoundsMax1[4];
477 double adfBoundsMin2[4], adfBoundsMax2[4];
478 double adfBoundsMin3[4], adfBoundsMax3[4];
479 double adfBoundsMin4[4], adfBoundsMax4[4];
481 SHPTreeSplitBounds( psTreeNode->adfBoundsMin,
482 psTreeNode->adfBoundsMax,
483 adfBoundsMinH1, adfBoundsMaxH1,
484 adfBoundsMinH2, adfBoundsMaxH2 );
486 SHPTreeSplitBounds( adfBoundsMinH1, adfBoundsMaxH1,
487 adfBoundsMin1, adfBoundsMax1,
488 adfBoundsMin2, adfBoundsMax2 );
490 SHPTreeSplitBounds( adfBoundsMinH2, adfBoundsMaxH2,
491 adfBoundsMin3, adfBoundsMax3,
492 adfBoundsMin4, adfBoundsMax4 );
494 if( SHPCheckObjectContained(psObject, nDimension,
495 adfBoundsMin1, adfBoundsMax1)
496 || SHPCheckObjectContained(psObject, nDimension,
497 adfBoundsMin2, adfBoundsMax2)
498 || SHPCheckObjectContained(psObject, nDimension,
499 adfBoundsMin3, adfBoundsMax3)
500 || SHPCheckObjectContained(psObject, nDimension,
501 adfBoundsMin4, adfBoundsMax4) )
503 psTreeNode->nSubNodes = 4;
504 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
506 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
508 psTreeNode->apsSubNode[2] = SHPTreeNodeCreate( adfBoundsMin3,
510 psTreeNode->apsSubNode[3] = SHPTreeNodeCreate( adfBoundsMin4,
513 /* recurse back on this node now that it has subnodes */
514 return( SHPTreeNodeAddShapeId( psTreeNode, psObject,
515 nMaxDepth, nDimension ) );
518 #endif /* MAX_SUBNODE == 4 */
520 /* -------------------------------------------------------------------- */
521 /* Otherwise, consider creating two subnodes if could fit into */
522 /* them, and adding to the appropriate subnode. */
523 /* -------------------------------------------------------------------- */
525 else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
527 double adfBoundsMin1[4], adfBoundsMax1[4];
528 double adfBoundsMin2[4], adfBoundsMax2[4];
530 SHPTreeSplitBounds( psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
531 adfBoundsMin1, adfBoundsMax1,
532 adfBoundsMin2, adfBoundsMax2 );
534 if( SHPCheckObjectContained(psObject, nDimension,
535 adfBoundsMin1, adfBoundsMax1))
537 psTreeNode->nSubNodes = 2;
538 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
540 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
543 return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[0], psObject,
544 nMaxDepth - 1, nDimension ) );
546 else if( SHPCheckObjectContained(psObject, nDimension,
547 adfBoundsMin2, adfBoundsMax2) )
549 psTreeNode->nSubNodes = 2;
550 psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
552 psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
555 return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[1], psObject,
556 nMaxDepth - 1, nDimension ) );
559 #endif /* MAX_SUBNODE == 2 */
561 /* -------------------------------------------------------------------- */
562 /* If none of that worked, just add it to this nodes list. */
563 /* -------------------------------------------------------------------- */
564 psTreeNode->nShapeCount++;
566 psTreeNode->panShapeIds = (int *)
567 SfRealloc( psTreeNode->panShapeIds,
568 sizeof(int) * psTreeNode->nShapeCount );
569 psTreeNode->panShapeIds[psTreeNode->nShapeCount-1] = psObject->nShapeId;
571 if( psTreeNode->papsShapeObj != NULL )
573 psTreeNode->papsShapeObj = (SHPObject **)
574 SfRealloc( psTreeNode->papsShapeObj,
575 sizeof(void *) * psTreeNode->nShapeCount );
576 psTreeNode->papsShapeObj[psTreeNode->nShapeCount-1] = NULL;
582 /************************************************************************/
583 /* SHPTreeAddShapeId() */
585 /* Add a shape to the tree, but don't keep a pointer to the */
586 /* object data, just keep the shapeid. */
587 /************************************************************************/
590 SHPTreeAddShapeId( SHPTree * psTree, SHPObject * psObject )
593 psTree->nTotalCount++;
595 return( SHPTreeNodeAddShapeId( psTree->psRoot, psObject,
596 psTree->nMaxDepth, psTree->nDimension ) );
599 /************************************************************************/
600 /* SHPTreeCollectShapesIds() */
602 /* Work function implementing SHPTreeFindLikelyShapes() on a */
603 /* tree node by tree node basis. */
604 /************************************************************************/
607 SHPTreeCollectShapeIds( SHPTree *hTree, SHPTreeNode * psTreeNode,
608 double * padfBoundsMin, double * padfBoundsMax,
609 int * pnShapeCount, int * pnMaxShapes,
610 int ** ppanShapeList )
615 /* -------------------------------------------------------------------- */
616 /* Does this node overlap the area of interest at all? If not, */
617 /* return without adding to the list at all. */
618 /* -------------------------------------------------------------------- */
619 if( !SHPCheckBoundsOverlap( psTreeNode->adfBoundsMin,
620 psTreeNode->adfBoundsMax,
623 hTree->nDimension ) )
626 /* -------------------------------------------------------------------- */
627 /* Grow the list to hold the shapes on this node. */
628 /* -------------------------------------------------------------------- */
629 if( *pnShapeCount + psTreeNode->nShapeCount > *pnMaxShapes )
631 *pnMaxShapes = (*pnShapeCount + psTreeNode->nShapeCount) * 2 + 20;
632 *ppanShapeList = (int *)
633 SfRealloc(*ppanShapeList,sizeof(int) * *pnMaxShapes);
636 /* -------------------------------------------------------------------- */
637 /* Add the local nodes shapeids to the list. */
638 /* -------------------------------------------------------------------- */
639 for( i = 0; i < psTreeNode->nShapeCount; i++ )
641 (*ppanShapeList)[(*pnShapeCount)++] = psTreeNode->panShapeIds[i];
644 /* -------------------------------------------------------------------- */
645 /* Recurse to subnodes if they exist. */
646 /* -------------------------------------------------------------------- */
647 for( i = 0; i < psTreeNode->nSubNodes; i++ )
649 if( psTreeNode->apsSubNode[i] != NULL )
650 SHPTreeCollectShapeIds( hTree, psTreeNode->apsSubNode[i],
651 padfBoundsMin, padfBoundsMax,
652 pnShapeCount, pnMaxShapes,
657 /************************************************************************/
658 /* SHPTreeFindLikelyShapes() */
660 /* Find all shapes within tree nodes for which the tree node */
661 /* bounding box overlaps the search box. The return value is */
662 /* an array of shapeids terminated by a -1. The shapeids will */
663 /* be in order, as hopefully this will result in faster (more */
664 /* sequential) reading from the file. */
665 /************************************************************************/
667 /* helper for qsort */
669 compare_ints( const void * a, const void * b)
671 return (*(int*)a) - (*(int*)b);
675 SHPTreeFindLikelyShapes( SHPTree * hTree,
676 double * padfBoundsMin, double * padfBoundsMax,
680 int *panShapeList=NULL, nMaxShapes = 0;
682 /* -------------------------------------------------------------------- */
683 /* Perform the search by recursive descent. */
684 /* -------------------------------------------------------------------- */
687 SHPTreeCollectShapeIds( hTree, hTree->psRoot,
688 padfBoundsMin, padfBoundsMax,
689 pnShapeCount, &nMaxShapes,
692 /* -------------------------------------------------------------------- */
693 /* Sort the id array */
694 /* -------------------------------------------------------------------- */
696 qsort(panShapeList, *pnShapeCount, sizeof(int), compare_ints);
701 /************************************************************************/
702 /* SHPTreeNodeTrim() */
704 /* This is the recurve version of SHPTreeTrimExtraNodes() that */
705 /* walks the tree cleaning it up. */
706 /************************************************************************/
708 static int SHPTreeNodeTrim( SHPTreeNode * psTreeNode )
713 /* -------------------------------------------------------------------- */
714 /* Trim subtrees, and free subnodes that come back empty. */
715 /* -------------------------------------------------------------------- */
716 for( i = 0; i < psTreeNode->nSubNodes; i++ )
718 if( SHPTreeNodeTrim( psTreeNode->apsSubNode[i] ) )
720 SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
722 psTreeNode->apsSubNode[i] =
723 psTreeNode->apsSubNode[psTreeNode->nSubNodes-1];
725 psTreeNode->nSubNodes--;
727 i--; /* process the new occupant of this subnode entry */
731 /* -------------------------------------------------------------------- */
732 /* If the current node has 1 subnode and no shapes, promote that */
733 /* subnode to the current node position. */
734 /* -------------------------------------------------------------------- */
735 if( psTreeNode->nSubNodes == 1 && psTreeNode->nShapeCount == 0)
737 SHPTreeNode* psSubNode = psTreeNode->apsSubNode[0];
739 memcpy(psTreeNode->adfBoundsMin, psSubNode->adfBoundsMin,
740 sizeof(psSubNode->adfBoundsMin));
741 memcpy(psTreeNode->adfBoundsMax, psSubNode->adfBoundsMax,
742 sizeof(psSubNode->adfBoundsMax));
743 psTreeNode->nShapeCount = psSubNode->nShapeCount;
744 assert(psTreeNode->panShapeIds == NULL);
745 psTreeNode->panShapeIds = psSubNode->panShapeIds;
746 assert(psTreeNode->papsShapeObj == NULL);
747 psTreeNode->papsShapeObj = psSubNode->papsShapeObj;
748 psTreeNode->nSubNodes = psSubNode->nSubNodes;
749 for( i = 0; i < psSubNode->nSubNodes; i++ )
750 psTreeNode->apsSubNode[i] = psSubNode->apsSubNode[i];
754 /* -------------------------------------------------------------------- */
755 /* We should be trimmed if we have no subnodes, and no shapes. */
756 /* -------------------------------------------------------------------- */
757 return( psTreeNode->nSubNodes == 0 && psTreeNode->nShapeCount == 0 );
760 /************************************************************************/
761 /* SHPTreeTrimExtraNodes() */
763 /* Trim empty nodes from the tree. Note that we never trim an */
764 /* empty root node. */
765 /************************************************************************/
768 SHPTreeTrimExtraNodes( SHPTree * hTree )
771 SHPTreeNodeTrim( hTree->psRoot );
774 /************************************************************************/
777 /* Swap a 2, 4 or 8 byte word. */
778 /************************************************************************/
780 static void SwapWord( int length, void * wordP )
786 for( i=0; i < length/2; i++ )
788 temp = ((unsigned char *) wordP)[i];
789 ((unsigned char *)wordP)[i] = ((unsigned char *) wordP)[length-i-1];
790 ((unsigned char *) wordP)[length-i-1] = temp;
795 struct SHPDiskTreeInfo
801 /************************************************************************/
802 /* SHPOpenDiskTree() */
803 /************************************************************************/
805 SHPTreeDiskHandle SHPOpenDiskTree( const char* pszQIXFilename,
808 SHPTreeDiskHandle hDiskTree;
810 hDiskTree = (SHPTreeDiskHandle) calloc(sizeof(struct SHPDiskTreeInfo),1);
813 SASetupDefaultHooks( &(hDiskTree->sHooks) );
815 memcpy( &(hDiskTree->sHooks), psHooks, sizeof(SAHooks) );
817 hDiskTree->fpQIX = hDiskTree->sHooks.FOpen(pszQIXFilename, "rb");
818 if (hDiskTree->fpQIX == NULL)
827 /***********************************************************************/
828 /* SHPCloseDiskTree() */
829 /************************************************************************/
831 void SHPCloseDiskTree( SHPTreeDiskHandle hDiskTree )
833 if (hDiskTree == NULL)
836 hDiskTree->sHooks.FClose(hDiskTree->fpQIX);
840 /************************************************************************/
841 /* SHPSearchDiskTreeNode() */
842 /************************************************************************/
845 SHPSearchDiskTreeNode( SHPTreeDiskHandle hDiskTree, double *padfBoundsMin, double *padfBoundsMax,
846 int **ppanResultBuffer, int *pnBufferMax,
847 int *pnResultCount, int bNeedSwap )
852 int numshapes, numsubnodes;
853 double adfNodeBoundsMin[2], adfNodeBoundsMax[2];
855 /* -------------------------------------------------------------------- */
856 /* Read and unswap first part of node info. */
857 /* -------------------------------------------------------------------- */
858 hDiskTree->sHooks.FRead( &offset, 4, 1, hDiskTree->fpQIX );
859 if ( bNeedSwap ) SwapWord ( 4, &offset );
861 hDiskTree->sHooks.FRead( adfNodeBoundsMin, sizeof(double), 2, hDiskTree->fpQIX );
862 hDiskTree->sHooks.FRead( adfNodeBoundsMax, sizeof(double), 2, hDiskTree->fpQIX );
865 SwapWord( 8, adfNodeBoundsMin + 0 );
866 SwapWord( 8, adfNodeBoundsMin + 1 );
867 SwapWord( 8, adfNodeBoundsMax + 0 );
868 SwapWord( 8, adfNodeBoundsMax + 1 );
871 hDiskTree->sHooks.FRead( &numshapes, 4, 1, hDiskTree->fpQIX );
872 if ( bNeedSwap ) SwapWord ( 4, &numshapes );
874 /* -------------------------------------------------------------------- */
875 /* If we don't overlap this node at all, we can just fseek() */
876 /* pass this node info and all subnodes. */
877 /* -------------------------------------------------------------------- */
878 if( !SHPCheckBoundsOverlap( adfNodeBoundsMin, adfNodeBoundsMax,
879 padfBoundsMin, padfBoundsMax, 2 ) )
881 offset += numshapes*sizeof(int) + sizeof(int);
882 hDiskTree->sHooks.FSeek(hDiskTree->fpQIX, offset, SEEK_CUR);
886 /* -------------------------------------------------------------------- */
887 /* Add all the shapeids at this node to our list. */
888 /* -------------------------------------------------------------------- */
891 if( *pnResultCount + numshapes > *pnBufferMax )
893 *pnBufferMax = (int) ((*pnResultCount + numshapes + 100) * 1.25);
894 *ppanResultBuffer = (int *)
895 SfRealloc( *ppanResultBuffer, *pnBufferMax * sizeof(int) );
898 hDiskTree->sHooks.FRead( *ppanResultBuffer + *pnResultCount,
899 sizeof(int), numshapes, hDiskTree->fpQIX );
903 for( i=0; i<numshapes; i++ )
904 SwapWord( 4, *ppanResultBuffer + *pnResultCount + i );
907 *pnResultCount += numshapes;
910 /* -------------------------------------------------------------------- */
911 /* Process the subnodes. */
912 /* -------------------------------------------------------------------- */
913 hDiskTree->sHooks.FRead( &numsubnodes, 4, 1, hDiskTree->fpQIX );
914 if ( bNeedSwap ) SwapWord ( 4, &numsubnodes );
916 for(i=0; i<numsubnodes; i++)
918 if( !SHPSearchDiskTreeNode( hDiskTree, padfBoundsMin, padfBoundsMax,
919 ppanResultBuffer, pnBufferMax,
920 pnResultCount, bNeedSwap ) )
927 /************************************************************************/
928 /* SHPTreeReadLibc() */
929 /************************************************************************/
932 SAOffset SHPTreeReadLibc( void *p, SAOffset size, SAOffset nmemb, SAFile file )
935 return (SAOffset) fread( p, (size_t) size, (size_t) nmemb,
939 /************************************************************************/
940 /* SHPTreeSeekLibc() */
941 /************************************************************************/
944 SAOffset SHPTreeSeekLibc( SAFile file, SAOffset offset, int whence )
947 return (SAOffset) fseek( (FILE *) file, (long) offset, whence );
950 /************************************************************************/
951 /* SHPSearchDiskTree() */
952 /************************************************************************/
955 SHPSearchDiskTree( FILE *fp,
956 double *padfBoundsMin, double *padfBoundsMax,
959 struct SHPDiskTreeInfo sDiskTree;
960 memset(&sDiskTree.sHooks, 0, sizeof(sDiskTree.sHooks));
962 /* We do not use SASetupDefaultHooks() because the FILE* */
963 /* is a libc FILE* */
964 sDiskTree.sHooks.FSeek = SHPTreeSeekLibc;
965 sDiskTree.sHooks.FRead = SHPTreeReadLibc;
967 sDiskTree.fpQIX = (SAFile)fp;
969 return SHPSearchDiskTreeEx( &sDiskTree, padfBoundsMin, padfBoundsMax,
973 /***********************************************************************/
974 /* SHPSearchDiskTreeEx() */
975 /************************************************************************/
977 int* SHPSearchDiskTreeEx( SHPTreeDiskHandle hDiskTree,
978 double *padfBoundsMin, double *padfBoundsMax,
982 int i, bNeedSwap, nBufferMax = 0;
983 unsigned char abyBuf[16];
984 int *panResultBuffer = NULL;
988 /* -------------------------------------------------------------------- */
989 /* Establish the byte order on this machine. */
990 /* -------------------------------------------------------------------- */
992 if( *((unsigned char *) &i) == 1 )
997 /* -------------------------------------------------------------------- */
998 /* Read the header. */
999 /* -------------------------------------------------------------------- */
1000 hDiskTree->sHooks.FSeek( hDiskTree->fpQIX, 0, SEEK_SET );
1001 hDiskTree->sHooks.FRead( abyBuf, 16, 1, hDiskTree->fpQIX );
1003 if( memcmp( abyBuf, "SQT", 3 ) != 0 )
1006 if( (abyBuf[3] == 2 && bBigEndian)
1007 || (abyBuf[3] == 1 && !bBigEndian) )
1012 /* -------------------------------------------------------------------- */
1013 /* Search through root node and it's decendents. */
1014 /* -------------------------------------------------------------------- */
1015 if( !SHPSearchDiskTreeNode( hDiskTree, padfBoundsMin, padfBoundsMax,
1016 &panResultBuffer, &nBufferMax,
1017 pnShapeCount, bNeedSwap ) )
1019 if( panResultBuffer != NULL )
1020 free( panResultBuffer );
1024 /* -------------------------------------------------------------------- */
1025 /* Sort the id array */
1026 /* -------------------------------------------------------------------- */
1027 qsort(panResultBuffer, *pnShapeCount, sizeof(int), compare_ints);
1029 return panResultBuffer;
1032 /************************************************************************/
1033 /* SHPGetSubNodeOffset() */
1035 /* Determine how big all the subnodes of this node (and their */
1036 /* children) will be. This will allow disk based searchers to */
1037 /* seek past them all efficiently. */
1038 /************************************************************************/
1040 static int SHPGetSubNodeOffset( SHPTreeNode *node)
1045 for(i=0; i<node->nSubNodes; i++ )
1047 if(node->apsSubNode[i])
1049 offset += 4*sizeof(double)
1050 + (node->apsSubNode[i]->nShapeCount+3)*sizeof(int);
1051 offset += SHPGetSubNodeOffset(node->apsSubNode[i]);
1058 /************************************************************************/
1059 /* SHPWriteTreeNode() */
1060 /************************************************************************/
1062 static void SHPWriteTreeNode( SAFile fp, SHPTreeNode *node, SAHooks* psHooks)
1066 unsigned char *pabyRec = NULL;
1067 assert( NULL != node );
1069 offset = SHPGetSubNodeOffset(node);
1071 pabyRec = (unsigned char *)
1072 malloc(sizeof(double) * 4
1073 + (3 * sizeof(int)) + (node->nShapeCount * sizeof(int)) );
1074 if( NULL == pabyRec )
1077 CPLError( CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
1083 memcpy( pabyRec, &offset, 4);
1085 /* minx, miny, maxx, maxy */
1086 memcpy( pabyRec+ 4, node->adfBoundsMin+0, sizeof(double) );
1087 memcpy( pabyRec+12, node->adfBoundsMin+1, sizeof(double) );
1088 memcpy( pabyRec+20, node->adfBoundsMax+0, sizeof(double) );
1089 memcpy( pabyRec+28, node->adfBoundsMax+1, sizeof(double) );
1091 memcpy( pabyRec+36, &node->nShapeCount, 4);
1092 j = node->nShapeCount * sizeof(int);
1093 memcpy( pabyRec+40, node->panShapeIds, j);
1094 memcpy( pabyRec+j+40, &node->nSubNodes, 4);
1096 psHooks->FWrite( pabyRec, 44+j, 1, fp );
1099 for(i=0; i<node->nSubNodes; i++ )
1101 if(node->apsSubNode[i])
1102 SHPWriteTreeNode( fp, node->apsSubNode[i], psHooks);
1106 /************************************************************************/
1107 /* SHPWriteTree() */
1108 /************************************************************************/
1110 int SHPAPI_CALL SHPWriteTree(SHPTree *tree, const char *filename )
1114 SASetupDefaultHooks( &sHooks );
1116 return SHPWriteTreeLL(tree, filename, &sHooks);
1119 /************************************************************************/
1120 /* SHPWriteTreeLL() */
1121 /************************************************************************/
1123 int SHPWriteTreeLL(SHPTree *tree, const char *filename, SAHooks* psHooks )
1125 char signature[4] = "SQT";
1131 if (psHooks == NULL)
1133 SASetupDefaultHooks( &sHooks );
1137 /* -------------------------------------------------------------------- */
1138 /* Open the output file. */
1139 /* -------------------------------------------------------------------- */
1140 fp = psHooks->FOpen(filename, "wb");
1146 /* -------------------------------------------------------------------- */
1147 /* Establish the byte order on this machine. */
1148 /* -------------------------------------------------------------------- */
1150 if( *((unsigned char *) &i) == 1 )
1155 /* -------------------------------------------------------------------- */
1156 /* Write the header. */
1157 /* -------------------------------------------------------------------- */
1158 memcpy( abyBuf+0, signature, 3 );
1161 abyBuf[3] = 2; /* New MSB */
1163 abyBuf[3] = 1; /* New LSB */
1165 abyBuf[4] = 1; /* version */
1166 abyBuf[5] = 0; /* next 3 reserved */
1170 psHooks->FWrite( abyBuf, 8, 1, fp );
1172 psHooks->FWrite( &(tree->nTotalCount), 4, 1, fp );
1174 /* write maxdepth */
1176 psHooks->FWrite( &(tree->nMaxDepth), 4, 1, fp );
1178 /* -------------------------------------------------------------------- */
1179 /* Write all the nodes "in order". */
1180 /* -------------------------------------------------------------------- */
1182 SHPWriteTreeNode( fp, tree->psRoot, psHooks );
1184 psHooks->FClose( fp );