/******************************************************************************
* $Id: shptree.c,v 1.20 2018-08-16 15:39:07 erouault Exp $
*
* Project: Shapelib
* Purpose: Implementation of quadtree building and searching functions.
* Author: Frank Warmerdam, warmerdam@pobox.com
*
******************************************************************************
* Copyright (c) 1999, Frank Warmerdam
* Copyright (c) 2012, Even Rouault <even dot rouault at mines-paris dot org>
*
* This software is available under the following "MIT Style" license,
* or at the option of the licensee under the LGPL (see COPYING). This
* option is discussed in more detail in shapelib.html.
*
* --
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
******************************************************************************
*
* $Log: shptree.c,v $
* Revision 1.20 2018-08-16 15:39:07 erouault
* * shpopen.c, dbfopen.c, shptree.c, sbnsearch.c: resyc with GDAL
* internal shapelib. Mostly to allow building those files as C++
* without warning. Also add FTDate entry in DBFFieldType
* (see https://github.com/OSGeo/gdal/pull/308). And some other
* code cleanups
*
* Revision 1.19 2016-12-05 12:44:06 erouault
* * Major overhaul of Makefile build system to use autoconf/automake.
*
* * Warning fixes in contrib/
*
* Revision 1.18 2016-12-04 15:30:15 erouault
* * shpopen.c, dbfopen.c, shptree.c, shapefil.h: resync with
* GDAL Shapefile driver. Mostly cleanups. SHPObject and DBFInfo
* structures extended with new members. New functions:
* DBFSetLastModifiedDate, SHPOpenLLEx, SHPRestoreSHX,
* SHPSetFastModeReadObject
*
* * sbnsearch.c: new file to implement original ESRI .sbn spatial
* index reading. (no write support). New functions:
* SBNOpenDiskTree, SBNCloseDiskTree, SBNSearchDiskTree,
* SBNSearchDiskTreeInteger, SBNSearchFreeIds
*
* * Makefile, makefile.vc, CMakeLists.txt, shapelib.def: updates
* with new file and symbols.
*
* * commit: helper script to cvs commit
*
* Revision 1.17 2012-01-27 21:09:26 fwarmerdam
* optimize .qix output (gdal #4472)
*
* Revision 1.16 2011-12-11 22:26:46 fwarmerdam
* upgrade .qix access code to use SAHooks (gdal #3365)
*
* Revision 1.15 2011-07-24 05:59:25 fwarmerdam
* minimize use of CPLError in favor of SAHooks.Error()
*
* Revision 1.14 2010-08-27 23:43:27 fwarmerdam
* add SHPAPI_CALL attribute in code
*
* Revision 1.13 2010-06-29 05:50:15 fwarmerdam
* fix sign of Z/M comparisons in SHPCheckObjectContained (#2223)
*
* Revision 1.12 2008-11-12 15:39:50 fwarmerdam
* improve safety in face of buggy .shp file.
*
* Revision 1.11 2007/10/27 03:31:14 fwarmerdam
* limit default depth of tree to 12 levels (gdal ticket #1594)
*
* Revision 1.10 2005/01/03 22:30:13 fwarmerdam
* added support for saved quadtrees
*
* Revision 1.9 2003/01/28 15:53:41 warmerda
* Avoid build warnings.
*
* Revision 1.8 2002/05/07 13:07:45 warmerda
* use qsort() - patch from Bernhard Herzog
*
* Revision 1.7 2002/01/15 14:36:07 warmerda
* updated email address
*
* Revision 1.6 2001/05/23 13:36:52 warmerda
* added use of SHPAPI_CALL
*
* Revision 1.5 1999/11/05 14:12:05 warmerda
* updated license terms
*
* Revision 1.4 1999/06/02 18:24:21 warmerda
* added trimming code
*
* Revision 1.3 1999/06/02 17:56:12 warmerda
* added quad'' subnode support for trees
*
* Revision 1.2 1999/05/18 19:11:11 warmerda
* Added example searching capability
*
* Revision 1.1 1999/05/18 17:49:20 warmerda
* New
*
*/
#include "shapefil.h"
#include <math.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#ifdef USE_CPL
#include "cpl_error.h"
#endif
SHP_CVSID("$Id: shptree.c,v 1.20 2018-08-16 15:39:07 erouault Exp $")
#ifndef TRUE
# define TRUE 1
# define FALSE 0
#endif
static int bBigEndian = 0;
/* -------------------------------------------------------------------- */
/* If the following is 0.5, nodes will be split in half. If it */
/* is 0.6 then each subnode will contain 60% of the parent */
/* node, with 20% representing overlap. This can be help to */
/* prevent small objects on a boundary from shifting too high */
/* up the tree. */
/* -------------------------------------------------------------------- */
#define SHP_SPLIT_RATIO 0.55
#ifdef __cplusplus
#define STATIC_CAST(type,x) static_cast<type>(x)
#define REINTERPRET_CAST(type,x) reinterpret_cast<type>(x)
#define CONST_CAST(type,x) const_cast<type>(x)
#define SHPLIB_NULLPTR nullptr
#else
#define STATIC_CAST(type,x) ((type)(x))
#define REINTERPRET_CAST(type,x) ((type)(x))
#define CONST_CAST(type,x) ((type)(x))
#define SHPLIB_NULLPTR NULL
#endif
/************************************************************************/
/* SfRealloc() */
/* */
/* A realloc cover function that will access a NULL pointer as */
/* a valid input. */
/************************************************************************/
static void * SfRealloc( void * pMem, int nNewSize )
{
if( pMem == SHPLIB_NULLPTR )
return malloc(nNewSize);
else
return realloc(pMem,nNewSize);
}
/************************************************************************/
/* SHPTreeNodeInit() */
/* */
/* Initialize a tree node. */
/************************************************************************/
static SHPTreeNode *SHPTreeNodeCreate( double * padfBoundsMin,
double * padfBoundsMax )
{
SHPTreeNode *psTreeNode;
psTreeNode = STATIC_CAST(SHPTreeNode *, malloc(sizeof(SHPTreeNode)));
if( SHPLIB_NULLPTR == psTreeNode )
return SHPLIB_NULLPTR;
psTreeNode->nShapeCount = 0;
psTreeNode->panShapeIds = SHPLIB_NULLPTR;
psTreeNode->papsShapeObj = SHPLIB_NULLPTR;
psTreeNode->nSubNodes = 0;
if( padfBoundsMin != SHPLIB_NULLPTR )
memcpy( psTreeNode->adfBoundsMin, padfBoundsMin, sizeof(double) * 4 );
if( padfBoundsMax != SHPLIB_NULLPTR )
memcpy( psTreeNode->adfBoundsMax, padfBoundsMax, sizeof(double) * 4 );
return psTreeNode;
}
/************************************************************************/
/* SHPCreateTree() */
/************************************************************************/
SHPTree SHPAPI_CALL1(*)
SHPCreateTree( SHPHandle hSHP, int nDimension, int nMaxDepth,
double *padfBoundsMin, double *padfBoundsMax )
{
SHPTree *psTree;
if( padfBoundsMin == SHPLIB_NULLPTR && hSHP == SHPLIB_NULLPTR )
return SHPLIB_NULLPTR;
/* -------------------------------------------------------------------- */
/* Allocate the tree object */
/* -------------------------------------------------------------------- */
psTree = STATIC_CAST(SHPTree *, malloc(sizeof(SHPTree)));
if( SHPLIB_NULLPTR == psTree )
{
return SHPLIB_NULLPTR;
}
psTree->hSHP = hSHP;
psTree->nMaxDepth = nMaxDepth;
psTree->nDimension = nDimension;
psTree->nTotalCount = 0;
/* -------------------------------------------------------------------- */
/* If no max depth was defined, try to select a reasonable one */
/* that implies approximately 8 shapes per node. */
/* -------------------------------------------------------------------- */
if( psTree->nMaxDepth == 0 && hSHP != SHPLIB_NULLPTR )
{
int nMaxNodeCount = 1;
int nShapeCount;
SHPGetInfo( hSHP, &nShapeCount, SHPLIB_NULLPTR, SHPLIB_NULLPTR, SHPLIB_NULLPTR );
while( nMaxNodeCount*4 < nShapeCount )
{
psTree->nMaxDepth += 1;
nMaxNodeCount = nMaxNodeCount * 2;
}
#ifdef USE_CPL
CPLDebug( "Shape",
"Estimated spatial index tree depth: %d",
psTree->nMaxDepth );
#endif
/* NOTE: Due to problems with memory allocation for deep trees,
* automatically estimated depth is limited up to 12 levels.
* See Ticket #1594 for detailed discussion.
*/
if( psTree->nMaxDepth > MAX_DEFAULT_TREE_DEPTH )
{
psTree->nMaxDepth = MAX_DEFAULT_TREE_DEPTH;
#ifdef USE_CPL
CPLDebug( "Shape",
"Falling back to max number of allowed index tree levels (%d).",
MAX_DEFAULT_TREE_DEPTH );
#endif
}
}
/* -------------------------------------------------------------------- */
/* Allocate the root node. */
/* -------------------------------------------------------------------- */
psTree->psRoot = SHPTreeNodeCreate( padfBoundsMin, padfBoundsMax );
if( SHPLIB_NULLPTR == psTree->psRoot )
{
free( psTree );
return SHPLIB_NULLPTR;
}
/* -------------------------------------------------------------------- */
/* Assign the bounds to the root node. If none are passed in, */
/* use the bounds of the provided file otherwise the create */
/* function will have already set the bounds. */
/* -------------------------------------------------------------------- */
if( padfBoundsMin == SHPLIB_NULLPTR )
{
SHPGetInfo( hSHP, SHPLIB_NULLPTR, SHPLIB_NULLPTR,
psTree->psRoot->adfBoundsMin,
psTree->psRoot->adfBoundsMax );
}
/* -------------------------------------------------------------------- */
/* If we have a file, insert all its shapes into the tree. */
/* -------------------------------------------------------------------- */
if( hSHP != SHPLIB_NULLPTR )
{
int iShape, nShapeCount;
SHPGetInfo( hSHP, &nShapeCount, SHPLIB_NULLPTR, SHPLIB_NULLPTR, SHPLIB_NULLPTR );
for( iShape = 0; iShape < nShapeCount; iShape++ )
{
SHPObject *psShape;
psShape = SHPReadObject( hSHP, iShape );
if( psShape != SHPLIB_NULLPTR )
{
SHPTreeAddShapeId( psTree, psShape );
SHPDestroyObject( psShape );
}
}
}
return psTree;
}
/************************************************************************/
/* SHPDestroyTreeNode() */
/************************************************************************/
static void SHPDestroyTreeNode( SHPTreeNode * psTreeNode )
{
int i;
assert( SHPLIB_NULLPTR != psTreeNode );
for( i = 0; i < psTreeNode->nSubNodes; i++ )
{
if( psTreeNode->apsSubNode[i] != SHPLIB_NULLPTR )
SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
}
if( psTreeNode->panShapeIds != SHPLIB_NULLPTR )
free( psTreeNode->panShapeIds );
if( psTreeNode->papsShapeObj != SHPLIB_NULLPTR )
{
for( i = 0; i < psTreeNode->nShapeCount; i++ )
{
if( psTreeNode->papsShapeObj[i] != SHPLIB_NULLPTR )
SHPDestroyObject( psTreeNode->papsShapeObj[i] );
}
free( psTreeNode->papsShapeObj );
}
free( psTreeNode );
}
/************************************************************************/
/* SHPDestroyTree() */
/************************************************************************/
void SHPAPI_CALL
SHPDestroyTree( SHPTree * psTree )
{
SHPDestroyTreeNode( psTree->psRoot );
free( psTree );
}
/************************************************************************/
/* SHPCheckBoundsOverlap() */
/* */
/* Do the given boxes overlap at all? */
/************************************************************************/
int SHPAPI_CALL
SHPCheckBoundsOverlap( double * padfBox1Min, double * padfBox1Max,
double * padfBox2Min, double * padfBox2Max,
int nDimension )
{
int iDim;
for( iDim = 0; iDim < nDimension; iDim++ )
{
if( padfBox2Max[iDim] < padfBox1Min[iDim] )
return FALSE;
if( padfBox1Max[iDim] < padfBox2Min[iDim] )
return FALSE;
}
return TRUE;
}
/************************************************************************/
/* SHPCheckObjectContained() */
/* */
/* Does the given shape fit within the indicated extents? */
/************************************************************************/
static int SHPCheckObjectContained( SHPObject * psObject, int nDimension,
double * padfBoundsMin, double * padfBoundsMax )
{
if( psObject->dfXMin < padfBoundsMin[0]
|| psObject->dfXMax > padfBoundsMax[0] )
return FALSE;
if( psObject->dfYMin < padfBoundsMin[1]
|| psObject->dfYMax > padfBoundsMax[1] )
return FALSE;
if( nDimension == 2 )
return TRUE;
if( psObject->dfZMin < padfBoundsMin[2]
|| psObject->dfZMax > padfBoundsMax[2] )
return FALSE;
if( nDimension == 3 )
return TRUE;
if( psObject->dfMMin < padfBoundsMin[3]
|| psObject->dfMMax > padfBoundsMax[3] )
return FALSE;
return TRUE;
}
/************************************************************************/
/* SHPTreeSplitBounds() */
/* */
/* Split a region into two subregion evenly, cutting along the */
/* longest dimension. */
/************************************************************************/
static void
SHPTreeSplitBounds( double *padfBoundsMinIn, double *padfBoundsMaxIn,
double *padfBoundsMin1, double * padfBoundsMax1,
double *padfBoundsMin2, double * padfBoundsMax2 )
{
/* -------------------------------------------------------------------- */
/* The output bounds will be very similar to the input bounds, */
/* so just copy over to start. */
/* -------------------------------------------------------------------- */
memcpy( padfBoundsMin1, padfBoundsMinIn, sizeof(double) * 4 );
memcpy( padfBoundsMax1, padfBoundsMaxIn, sizeof(double) * 4 );
memcpy( padfBoundsMin2, padfBoundsMinIn, sizeof(double) * 4 );
memcpy( padfBoundsMax2, padfBoundsMaxIn, sizeof(double) * 4 );
/* -------------------------------------------------------------------- */
/* Split in X direction. */
/* -------------------------------------------------------------------- */
if( (padfBoundsMaxIn[0] - padfBoundsMinIn[0])
> (padfBoundsMaxIn[1] - padfBoundsMinIn[1]) )
{
double dfRange = padfBoundsMaxIn[0] - padfBoundsMinIn[0];
padfBoundsMax1[0] = padfBoundsMinIn[0] + dfRange * SHP_SPLIT_RATIO;
padfBoundsMin2[0] = padfBoundsMaxIn[0] - dfRange * SHP_SPLIT_RATIO;
}
/* -------------------------------------------------------------------- */
/* Otherwise split in Y direction. */
/* -------------------------------------------------------------------- */
else
{
double dfRange = padfBoundsMaxIn[1] - padfBoundsMinIn[1];
padfBoundsMax1[1] = padfBoundsMinIn[1] + dfRange * SHP_SPLIT_RATIO;
padfBoundsMin2[1] = padfBoundsMaxIn[1] - dfRange * SHP_SPLIT_RATIO;
}
}
/************************************************************************/
/* SHPTreeNodeAddShapeId() */
/************************************************************************/
static int
SHPTreeNodeAddShapeId( SHPTreeNode * psTreeNode, SHPObject * psObject,
int nMaxDepth, int nDimension )
{
int i;
/* -------------------------------------------------------------------- */
/* If there are subnodes, then consider whether this object */
/* will fit in them. */
/* -------------------------------------------------------------------- */
if( nMaxDepth > 1 && psTreeNode->nSubNodes > 0 )
{
for( i = 0; i < psTreeNode->nSubNodes; i++ )
{
if( SHPCheckObjectContained(psObject, nDimension,
psTreeNode->apsSubNode[i]->adfBoundsMin,
psTreeNode->apsSubNode[i]->adfBoundsMax))
{
return SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[i],
psObject, nMaxDepth-1,
nDimension );
}
}
}
/* -------------------------------------------------------------------- */
/* Otherwise, consider creating four subnodes if could fit into */
/* them, and adding to the appropriate subnode. */
/* -------------------------------------------------------------------- */
#if MAX_SUBNODE == 4
else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
{
double adfBoundsMinH1[4], adfBoundsMaxH1[4];
double adfBoundsMinH2[4], adfBoundsMaxH2[4];
double adfBoundsMin1[4], adfBoundsMax1[4];
double adfBoundsMin2[4], adfBoundsMax2[4];
double adfBoundsMin3[4], adfBoundsMax3[4];
double adfBoundsMin4[4], adfBoundsMax4[4];
SHPTreeSplitBounds( psTreeNode->adfBoundsMin,
psTreeNode->adfBoundsMax,
adfBoundsMinH1, adfBoundsMaxH1,
adfBoundsMinH2, adfBoundsMaxH2 );
SHPTreeSplitBounds( adfBoundsMinH1, adfBoundsMaxH1,
adfBoundsMin1, adfBoundsMax1,
adfBoundsMin2, adfBoundsMax2 );
SHPTreeSplitBounds( adfBoundsMinH2, adfBoundsMaxH2,
adfBoundsMin3, adfBoundsMax3,
adfBoundsMin4, adfBoundsMax4 );
if( SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin1, adfBoundsMax1)
|| SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin2, adfBoundsMax2)
|| SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin3, adfBoundsMax3)
|| SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin4, adfBoundsMax4) )
{
psTreeNode->nSubNodes = 4;
psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
adfBoundsMax1 );
psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
adfBoundsMax2 );
psTreeNode->apsSubNode[2] = SHPTreeNodeCreate( adfBoundsMin3,
adfBoundsMax3 );
psTreeNode->apsSubNode[3] = SHPTreeNodeCreate( adfBoundsMin4,
adfBoundsMax4 );
/* recurse back on this node now that it has subnodes */
return( SHPTreeNodeAddShapeId( psTreeNode, psObject,
nMaxDepth, nDimension ) );
}
}
#endif /* MAX_SUBNODE == 4 */
/* -------------------------------------------------------------------- */
/* Otherwise, consider creating two subnodes if could fit into */
/* them, and adding to the appropriate subnode. */
/* -------------------------------------------------------------------- */
#if MAX_SUBNODE == 2
else if( nMaxDepth > 1 && psTreeNode->nSubNodes == 0 )
{
double adfBoundsMin1[4], adfBoundsMax1[4];
double adfBoundsMin2[4], adfBoundsMax2[4];
SHPTreeSplitBounds( psTreeNode->adfBoundsMin, psTreeNode->adfBoundsMax,
adfBoundsMin1, adfBoundsMax1,
adfBoundsMin2, adfBoundsMax2 );
if( SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin1, adfBoundsMax1))
{
psTreeNode->nSubNodes = 2;
psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
adfBoundsMax1 );
psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
adfBoundsMax2 );
return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[0], psObject,
nMaxDepth - 1, nDimension ) );
}
else if( SHPCheckObjectContained(psObject, nDimension,
adfBoundsMin2, adfBoundsMax2) )
{
psTreeNode->nSubNodes = 2;
psTreeNode->apsSubNode[0] = SHPTreeNodeCreate( adfBoundsMin1,
adfBoundsMax1 );
psTreeNode->apsSubNode[1] = SHPTreeNodeCreate( adfBoundsMin2,
adfBoundsMax2 );
return( SHPTreeNodeAddShapeId( psTreeNode->apsSubNode[1], psObject,
nMaxDepth - 1, nDimension ) );
}
}
#endif /* MAX_SUBNODE == 2 */
/* -------------------------------------------------------------------- */
/* If none of that worked, just add it to this nodes list. */
/* -------------------------------------------------------------------- */
psTreeNode->nShapeCount++;
psTreeNode->panShapeIds = STATIC_CAST(int *,
SfRealloc( psTreeNode->panShapeIds,
sizeof(int) * psTreeNode->nShapeCount ));
psTreeNode->panShapeIds[psTreeNode->nShapeCount-1] = psObject->nShapeId;
if( psTreeNode->papsShapeObj != SHPLIB_NULLPTR )
{
psTreeNode->papsShapeObj = STATIC_CAST(SHPObject **,
SfRealloc( psTreeNode->papsShapeObj,
sizeof(void *) * psTreeNode->nShapeCount ));
psTreeNode->papsShapeObj[psTreeNode->nShapeCount-1] = SHPLIB_NULLPTR;
}
return TRUE;
}
/************************************************************************/
/* SHPTreeAddShapeId() */
/* */
/* Add a shape to the tree, but don't keep a pointer to the */
/* object data, just keep the shapeid. */
/************************************************************************/
int SHPAPI_CALL
SHPTreeAddShapeId( SHPTree * psTree, SHPObject * psObject )
{
psTree->nTotalCount++;
return( SHPTreeNodeAddShapeId( psTree->psRoot, psObject,
psTree->nMaxDepth, psTree->nDimension ) );
}
/************************************************************************/
/* SHPTreeCollectShapesIds() */
/* */
/* Work function implementing SHPTreeFindLikelyShapes() on a */
/* tree node by tree node basis. */
/************************************************************************/
static void
SHPTreeCollectShapeIds( SHPTree *hTree, SHPTreeNode * psTreeNode,
double * padfBoundsMin, double * padfBoundsMax,
int * pnShapeCount, int * pnMaxShapes,
int ** ppanShapeList )
{
int i;
/* -------------------------------------------------------------------- */
/* Does this node overlap the area of interest at all? If not, */
/* return without adding to the list at all. */
/* -------------------------------------------------------------------- */
if( !SHPCheckBoundsOverlap( psTreeNode->adfBoundsMin,
psTreeNode->adfBoundsMax,
padfBoundsMin,
padfBoundsMax,
hTree->nDimension ) )
return;
/* -------------------------------------------------------------------- */
/* Grow the list to hold the shapes on this node. */
/* -------------------------------------------------------------------- */
if( *pnShapeCount + psTreeNode->nShapeCount > *pnMaxShapes )
{
*pnMaxShapes = (*pnShapeCount + psTreeNode->nShapeCount) * 2 + 20;
*ppanShapeList = STATIC_CAST(int *,
SfRealloc(*ppanShapeList,sizeof(int) * *pnMaxShapes));
}
/* -------------------------------------------------------------------- */
/* Add the local nodes shapeids to the list. */
/* -------------------------------------------------------------------- */
for( i = 0; i < psTreeNode->nShapeCount; i++ )
{
(*ppanShapeList)[(*pnShapeCount)++] = psTreeNode->panShapeIds[i];
}
/* -------------------------------------------------------------------- */
/* Recurse to subnodes if they exist. */
/* -------------------------------------------------------------------- */
for( i = 0; i < psTreeNode->nSubNodes; i++ )
{
if( psTreeNode->apsSubNode[i] != SHPLIB_NULLPTR )
SHPTreeCollectShapeIds( hTree, psTreeNode->apsSubNode[i],
padfBoundsMin, padfBoundsMax,
pnShapeCount, pnMaxShapes,
ppanShapeList );
}
}
/************************************************************************/
/* SHPTreeFindLikelyShapes() */
/* */
/* Find all shapes within tree nodes for which the tree node */
/* bounding box overlaps the search box. The return value is */
/* an array of shapeids terminated by a -1. The shapeids will */
/* be in order, as hopefully this will result in faster (more */
/* sequential) reading from the file. */
/************************************************************************/
/* helper for qsort */
static int
compare_ints( const void * a, const void * b)
{
return *REINTERPRET_CAST(const int*, a) - *REINTERPRET_CAST(const int*, b);
}
int SHPAPI_CALL1(*)
SHPTreeFindLikelyShapes( SHPTree * hTree,
double * padfBoundsMin, double * padfBoundsMax,
int * pnShapeCount )
{
int *panShapeList=SHPLIB_NULLPTR, nMaxShapes = 0;
/* -------------------------------------------------------------------- */
/* Perform the search by recursive descent. */
/* -------------------------------------------------------------------- */
*pnShapeCount = 0;
SHPTreeCollectShapeIds( hTree, hTree->psRoot,
padfBoundsMin, padfBoundsMax,
pnShapeCount, &nMaxShapes,
&panShapeList );
/* -------------------------------------------------------------------- */
/* Sort the id array */
/* -------------------------------------------------------------------- */
if( panShapeList != SHPLIB_NULLPTR )
qsort(panShapeList, *pnShapeCount, sizeof(int), compare_ints);
return panShapeList;
}
/************************************************************************/
/* SHPTreeNodeTrim() */
/* */
/* This is the recursive version of SHPTreeTrimExtraNodes() that */
/* walks the tree cleaning it up. */
/************************************************************************/
static int SHPTreeNodeTrim( SHPTreeNode * psTreeNode )
{
int i;
/* -------------------------------------------------------------------- */
/* Trim subtrees, and free subnodes that come back empty. */
/* -------------------------------------------------------------------- */
for( i = 0; i < psTreeNode->nSubNodes; i++ )
{
if( SHPTreeNodeTrim( psTreeNode->apsSubNode[i] ) )
{
SHPDestroyTreeNode( psTreeNode->apsSubNode[i] );
psTreeNode->apsSubNode[i] =
psTreeNode->apsSubNode[psTreeNode->nSubNodes-1];
psTreeNode->nSubNodes--;
i--; /* process the new occupant of this subnode entry */
}
}
/* -------------------------------------------------------------------- */
/* If the current node has 1 subnode and no shapes, promote that */
/* subnode to the current node position. */
/* -------------------------------------------------------------------- */
if( psTreeNode->nSubNodes == 1 && psTreeNode->nShapeCount == 0)
{
SHPTreeNode* psSubNode = psTreeNode->apsSubNode[0];
memcpy(psTreeNode->adfBoundsMin, psSubNode->adfBoundsMin,
sizeof(psSubNode->adfBoundsMin));
memcpy(psTreeNode->adfBoundsMax, psSubNode->adfBoundsMax,
sizeof(psSubNode->adfBoundsMax));
psTreeNode->nShapeCount = psSubNode->nShapeCount;
assert(psTreeNode->panShapeIds == SHPLIB_NULLPTR);
psTreeNode->panShapeIds = psSubNode->panShapeIds;
assert(psTreeNode->papsShapeObj == SHPLIB_NULLPTR);
psTreeNode->papsShapeObj = psSubNode->papsShapeObj;
psTreeNode->nSubNodes = psSubNode->nSubNodes;
for( i = 0; i < psSubNode->nSubNodes; i++ )
psTreeNode->apsSubNode[i] = psSubNode->apsSubNode[i];
free(psSubNode);
}
/* -------------------------------------------------------------------- */
/* We should be trimmed if we have no subnodes, and no shapes. */
/* -------------------------------------------------------------------- */
return( psTreeNode->nSubNodes == 0 && psTreeNode->nShapeCount == 0 );
}
/************************************************************************/
/* SHPTreeTrimExtraNodes() */
/* */
/* Trim empty nodes from the tree. Note that we never trim an */
/* empty root node. */
/************************************************************************/
void SHPAPI_CALL
SHPTreeTrimExtraNodes( SHPTree * hTree )
{
SHPTreeNodeTrim( hTree->psRoot );
}
/************************************************************************/
/* SwapWord() */
/* */
/* Swap a 2, 4 or 8 byte word. */
/************************************************************************/
static void SwapWord( int length, void * wordP )
{
int i;
unsigned char temp;
for( i=0; i < length/2; i++ )
{
temp = STATIC_CAST(unsigned char*, wordP)[i];
STATIC_CAST(unsigned char*, wordP)[i] = STATIC_CAST(unsigned char*, wordP)[length-i-1];
STATIC_CAST(unsigned char*, wordP)[length-i-1] = temp;
}
}
struct SHPDiskTreeInfo
{
SAHooks sHooks;
SAFile fpQIX;
};
/************************************************************************/
/* SHPOpenDiskTree() */
/************************************************************************/
SHPTreeDiskHandle SHPOpenDiskTree( const char* pszQIXFilename,
SAHooks *psHooks )
{
SHPTreeDiskHandle hDiskTree;
hDiskTree = STATIC_CAST(SHPTreeDiskHandle, calloc(sizeof(struct SHPDiskTreeInfo),1));
if (psHooks == SHPLIB_NULLPTR)
SASetupDefaultHooks( &(hDiskTree->sHooks) );
else
memcpy( &(hDiskTree->sHooks), psHooks, sizeof(SAHooks) );
hDiskTree->fpQIX = hDiskTree->sHooks.FOpen(pszQIXFilename, "rb");
if (hDiskTree->fpQIX == SHPLIB_NULLPTR)
{
free(hDiskTree);
return SHPLIB_NULLPTR;
}
return hDiskTree;
}
/***********************************************************************/
/* SHPCloseDiskTree() */
/************************************************************************/
void SHPCloseDiskTree( SHPTreeDiskHandle hDiskTree )
{
if (hDiskTree == SHPLIB_NULLPTR)
return;
hDiskTree->sHooks.FClose(hDiskTree->fpQIX);
free(hDiskTree);
}
/************************************************************************/
/* SHPSearchDiskTreeNode() */
/************************************************************************/
static int
SHPSearchDiskTreeNode( SHPTreeDiskHandle hDiskTree, double *padfBoundsMin, double *padfBoundsMax,
int **ppanResultBuffer, int *pnBufferMax,
int *pnResultCount, int bNeedSwap, int nRecLevel )
{
unsigned int i;
unsigned int offset;
unsigned int numshapes, numsubnodes;
double adfNodeBoundsMin[2], adfNodeBoundsMax[2];
int nFReadAcc;
/* -------------------------------------------------------------------- */
/* Read and unswap first part of node info. */
/* -------------------------------------------------------------------- */
nFReadAcc = STATIC_CAST(int, hDiskTree->sHooks.FRead( &offset, 4, 1, hDiskTree->fpQIX ));
if ( bNeedSwap ) SwapWord ( 4, &offset );
nFReadAcc += STATIC_CAST(int, hDiskTree->sHooks.FRead( adfNodeBoundsMin, sizeof(double), 2, hDiskTree->fpQIX ));
nFReadAcc += STATIC_CAST(int, hDiskTree->sHooks.FRead( adfNodeBoundsMax, sizeof(double), 2, hDiskTree->fpQIX ));
if ( bNeedSwap )
{
SwapWord( 8, adfNodeBoundsMin + 0 );
SwapWord( 8, adfNodeBoundsMin + 1 );
SwapWord( 8, adfNodeBoundsMax + 0 );
SwapWord( 8, adfNodeBoundsMax + 1 );
}
nFReadAcc += STATIC_CAST(int, hDiskTree->sHooks.FRead( &numshapes, 4, 1, hDiskTree->fpQIX ));
if ( bNeedSwap ) SwapWord ( 4, &numshapes );
/* Check that we could read all previous values */
if( nFReadAcc != 1 + 2 + 2 + 1 )
{
hDiskTree->sHooks.Error("I/O error");
return FALSE;
}
/* Sanity checks to avoid int overflows in later computation */
if( offset > INT_MAX - sizeof(int) )
{
hDiskTree->sHooks.Error("Invalid value for offset");
return FALSE;
}
if( numshapes > (INT_MAX - offset - sizeof(int)) / sizeof(int) ||
numshapes > INT_MAX / sizeof(int) - *pnResultCount )
{
hDiskTree->sHooks.Error("Invalid value for numshapes");
return FALSE;
}
/* -------------------------------------------------------------------- */
/* If we don't overlap this node at all, we can just fseek() */
/* pass this node info and all subnodes. */
/* -------------------------------------------------------------------- */
if( !SHPCheckBoundsOverlap( adfNodeBoundsMin, adfNodeBoundsMax,
padfBoundsMin, padfBoundsMax, 2 ) )
{
offset += numshapes*sizeof(int) + sizeof(int);
hDiskTree->sHooks.FSeek(hDiskTree->fpQIX, offset, SEEK_CUR);
return TRUE;
}
/* -------------------------------------------------------------------- */
/* Add all the shapeids at this node to our list. */
/* -------------------------------------------------------------------- */
if(numshapes > 0)
{
if( *pnResultCount + numshapes > STATIC_CAST(unsigned int, *pnBufferMax) )
{
int* pNewBuffer;
*pnBufferMax = (*pnResultCount + numshapes + 100) * 5 / 4;
if( STATIC_CAST(size_t, *pnBufferMax) > INT_MAX / sizeof(int) )
*pnBufferMax = *pnResultCount + numshapes;
pNewBuffer = STATIC_CAST(int *,
SfRealloc( *ppanResultBuffer, *pnBufferMax * sizeof(int) ));
if( pNewBuffer == SHPLIB_NULLPTR )
{
hDiskTree->sHooks.Error("Out of memory error");
return FALSE;
}
*ppanResultBuffer = pNewBuffer;
}
if( hDiskTree->sHooks.FRead( *ppanResultBuffer + *pnResultCount,
sizeof(int), numshapes, hDiskTree->fpQIX ) != numshapes )
{
hDiskTree->sHooks.Error("I/O error");
return FALSE;
}
if (bNeedSwap )
{
for( i=0; i<numshapes; i++ )
SwapWord( 4, *ppanResultBuffer + *pnResultCount + i );
}
*pnResultCount += numshapes;
}
/* -------------------------------------------------------------------- */
/* Process the subnodes. */
/* -------------------------------------------------------------------- */
if( hDiskTree->sHooks.FRead( &numsubnodes, 4, 1, hDiskTree->fpQIX ) != 1 )
{
hDiskTree->sHooks.Error("I/O error");
return FALSE;
}
if ( bNeedSwap ) SwapWord ( 4, &numsubnodes );
if( numsubnodes > 0 && nRecLevel == 32 )
{
hDiskTree->sHooks.Error("Shape tree is too deep");
return FALSE;
}
for(i=0; i<numsubnodes; i++)
{
if( !SHPSearchDiskTreeNode( hDiskTree, padfBoundsMin, padfBoundsMax,
ppanResultBuffer, pnBufferMax,
pnResultCount, bNeedSwap, nRecLevel + 1 ) )
return FALSE;
}
return TRUE;
}
/************************************************************************/
/* SHPTreeReadLibc() */
/************************************************************************/
static
SAOffset SHPTreeReadLibc( void *p, SAOffset size, SAOffset nmemb, SAFile file )
{
return STATIC_CAST(SAOffset, fread( p, STATIC_CAST(size_t, size),
STATIC_CAST(size_t, nmemb),
REINTERPRET_CAST(FILE*, file) ));
}
/************************************************************************/
/* SHPTreeSeekLibc() */
/************************************************************************/
static
SAOffset SHPTreeSeekLibc( SAFile file, SAOffset offset, int whence )
{
return STATIC_CAST(SAOffset, fseek( REINTERPRET_CAST(FILE*, file),
STATIC_CAST(long, offset), whence ));
}
/************************************************************************/
/* SHPSearchDiskTree() */
/************************************************************************/
int SHPAPI_CALL1(*)
SHPSearchDiskTree( FILE *fp,
double *padfBoundsMin, double *padfBoundsMax,
int *pnShapeCount )
{
struct SHPDiskTreeInfo sDiskTree;
memset(&sDiskTree.sHooks, 0, sizeof(sDiskTree.sHooks));
/* We do not use SASetupDefaultHooks() because the FILE* */
/* is a libc FILE* */
sDiskTree.sHooks.FSeek = SHPTreeSeekLibc;
sDiskTree.sHooks.FRead = SHPTreeReadLibc;
sDiskTree.fpQIX = REINTERPRET_CAST(SAFile, fp);
return SHPSearchDiskTreeEx( &sDiskTree, padfBoundsMin, padfBoundsMax,
pnShapeCount );
}
/***********************************************************************/
/* SHPSearchDiskTreeEx() */
/************************************************************************/
int* SHPSearchDiskTreeEx( SHPTreeDiskHandle hDiskTree,
double *padfBoundsMin, double *padfBoundsMax,
int *pnShapeCount )
{
int i, bNeedSwap, nBufferMax = 0;
unsigned char abyBuf[16];
int *panResultBuffer = SHPLIB_NULLPTR;
*pnShapeCount = 0;
/* -------------------------------------------------------------------- */
/* Establish the byte order on this machine. */
/* -------------------------------------------------------------------- */
i = 1;
if( *REINTERPRET_CAST(unsigned char *, &i) == 1 )
bBigEndian = FALSE;
else
bBigEndian = TRUE;
/* -------------------------------------------------------------------- */
/* Read the header. */
/* -------------------------------------------------------------------- */
hDiskTree->sHooks.FSeek( hDiskTree->fpQIX, 0, SEEK_SET );
hDiskTree->sHooks.FRead( abyBuf, 16, 1, hDiskTree->fpQIX );
if( memcmp( abyBuf, "SQT", 3 ) != 0 )
return SHPLIB_NULLPTR;
if( (abyBuf[3] == 2 && bBigEndian)
|| (abyBuf[3] == 1 && !bBigEndian) )
bNeedSwap = FALSE;
else
bNeedSwap = TRUE;
/* -------------------------------------------------------------------- */
/* Search through root node and its descendants. */
/* -------------------------------------------------------------------- */
if( !SHPSearchDiskTreeNode( hDiskTree, padfBoundsMin, padfBoundsMax,
&panResultBuffer, &nBufferMax,
pnShapeCount, bNeedSwap, 0 ) )
{
if( panResultBuffer != SHPLIB_NULLPTR )
free( panResultBuffer );
*pnShapeCount = 0;
return SHPLIB_NULLPTR;
}
/* -------------------------------------------------------------------- */
/* Sort the id array */
/* -------------------------------------------------------------------- */
/* To distinguish between empty intersection from error case */
if( panResultBuffer == SHPLIB_NULLPTR )
panResultBuffer = STATIC_CAST(int*, calloc(1, sizeof(int)));
else
qsort(panResultBuffer, *pnShapeCount, sizeof(int), compare_ints);
return panResultBuffer;
}
/************************************************************************/
/* SHPGetSubNodeOffset() */
/* */
/* Determine how big all the subnodes of this node (and their */
/* children) will be. This will allow disk based searchers to */
/* seek past them all efficiently. */
/************************************************************************/
static int SHPGetSubNodeOffset( SHPTreeNode *node)
{
int i;
int offset=0;
for(i=0; i<node->nSubNodes; i++ )
{
if(node->apsSubNode[i])
{
offset += 4*sizeof(double)
+ (node->apsSubNode[i]->nShapeCount+3)*sizeof(int);
offset += SHPGetSubNodeOffset(node->apsSubNode[i]);
}
}
return(offset);
}
/************************************************************************/
/* SHPWriteTreeNode() */
/************************************************************************/
static void SHPWriteTreeNode( SAFile fp, SHPTreeNode *node, SAHooks* psHooks)
{
int i,j;
int offset;
unsigned char *pabyRec;
assert( SHPLIB_NULLPTR != node );
offset = SHPGetSubNodeOffset(node);
pabyRec = STATIC_CAST(unsigned char *,
malloc(sizeof(double) * 4
+ (3 * sizeof(int)) + (node->nShapeCount * sizeof(int)) ));
if( SHPLIB_NULLPTR == pabyRec )
{
#ifdef USE_CPL
CPLError( CE_Fatal, CPLE_OutOfMemory, "Memory allocation failure");
#endif
assert( 0 );
return;
}
memcpy( pabyRec, &offset, 4);
/* minx, miny, maxx, maxy */
memcpy( pabyRec+ 4, node->adfBoundsMin+0, sizeof(double) );
memcpy( pabyRec+12, node->adfBoundsMin+1, sizeof(double) );
memcpy( pabyRec+20, node->adfBoundsMax+0, sizeof(double) );
memcpy( pabyRec+28, node->adfBoundsMax+1, sizeof(double) );
memcpy( pabyRec+36, &node->nShapeCount, 4);
j = node->nShapeCount * sizeof(int);
if( j )
memcpy( pabyRec+40, node->panShapeIds, j);
memcpy( pabyRec+j+40, &node->nSubNodes, 4);
psHooks->FWrite( pabyRec, 44+j, 1, fp );
free (pabyRec);
for(i=0; i<node->nSubNodes; i++ )
{
if(node->apsSubNode[i])
SHPWriteTreeNode( fp, node->apsSubNode[i], psHooks);
}
}
/************************************************************************/
/* SHPWriteTree() */
/************************************************************************/
int SHPAPI_CALL SHPWriteTree(SHPTree *tree, const char *filename )
{
SAHooks sHooks;
SASetupDefaultHooks( &sHooks );
return SHPWriteTreeLL(tree, filename, &sHooks);
}
/************************************************************************/
/* SHPWriteTreeLL() */
/************************************************************************/
int SHPWriteTreeLL(SHPTree *tree, const char *filename, SAHooks* psHooks )
{
char signature[4] = "SQT";
int i;
char abyBuf[32];
SAFile fp;
SAHooks sHooks;
if (psHooks == SHPLIB_NULLPTR)
{
SASetupDefaultHooks( &sHooks );
psHooks = &sHooks;
}
/* -------------------------------------------------------------------- */
/* Open the output file. */
/* -------------------------------------------------------------------- */
fp = psHooks->FOpen(filename, "wb");
if( fp == SHPLIB_NULLPTR )
{
return FALSE;
}
/* -------------------------------------------------------------------- */
/* Establish the byte order on this machine. */
/* -------------------------------------------------------------------- */
i = 1;
if( *REINTERPRET_CAST(unsigned char *, &i) == 1 )
bBigEndian = FALSE;
else
bBigEndian = TRUE;
/* -------------------------------------------------------------------- */
/* Write the header. */
/* -------------------------------------------------------------------- */
memcpy( abyBuf+0, signature, 3 );
if( bBigEndian )
abyBuf[3] = 2; /* New MSB */
else
abyBuf[3] = 1; /* New LSB */
abyBuf[4] = 1; /* version */
abyBuf[5] = 0; /* next 3 reserved */
abyBuf[6] = 0;
abyBuf[7] = 0;
psHooks->FWrite( abyBuf, 8, 1, fp );
psHooks->FWrite( &(tree->nTotalCount), 4, 1, fp );
/* write maxdepth */
psHooks->FWrite( &(tree->nMaxDepth), 4, 1, fp );
/* -------------------------------------------------------------------- */
/* Write all the nodes "in order". */
/* -------------------------------------------------------------------- */
SHPWriteTreeNode( fp, tree->psRoot, psHooks );
psHooks->FClose( fp );
return TRUE;
}