src/shapelib/shptree.c

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