1 : /******************************************************************************
2 : * $Id: sbnsearch.c 24638 2012-07-01 16:31:40Z rouault $
3 : *
4 : * Project: Shapelib
5 : * Purpose: Implementation of search in ESRI SBN spatial index.
6 : * Author: Even Rouault, even dot rouault at mines dash paris dot org
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2012, Even Rouault
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 : #include "shapefil.h"
37 :
38 : #include <math.h>
39 : #include <assert.h>
40 : #include <stdlib.h>
41 : #include <string.h>
42 :
43 : SHP_CVSID("$Id: sbnsearch.c 24638 2012-07-01 16:31:40Z rouault $")
44 :
45 : #ifndef TRUE
46 : # define TRUE 1
47 : # define FALSE 0
48 : #endif
49 :
50 : #define READ_MSB_INT(ptr) \
51 : (((ptr)[0] << 24) | ((ptr)[1] << 16) | ((ptr)[2] << 8) | (ptr)[3])
52 :
53 : #define CACHED_DEPTH_LIMIT 8
54 :
55 : typedef unsigned char uchar;
56 :
57 : typedef int coord;
58 : /*typedef uchar coord;*/
59 :
60 : typedef struct
61 : {
62 : uchar *pabyShapeDesc; /* Cache of (nShapeCount * 8) bytes of the bins. May be NULL. */
63 : int nBinStart; /* Index of first bin for this node. */
64 : int nShapeCount; /* Number of shapes attached to this node. */
65 : int nBinCount; /* Number of bins for this node. May be 0 if node is empty. */
66 : int nBinOffset; /* Offset in file of the start of the first bin. May be 0 if node is empty. */
67 :
68 : int bBBoxInit; /* TRUE if the following bounding box has been computed. */
69 : coord bMinX; /* Bounding box of the shapes directly attached to this node. */
70 : coord bMinY; /* This is *not* the theoretical footprint of the node. */
71 : coord bMaxX;
72 : coord bMaxY;
73 : } SBNNodeDescriptor;
74 :
75 : struct SBNSearchInfo
76 : {
77 : SAHooks sHooks;
78 : SAFile fpSBN;
79 : SBNNodeDescriptor *pasNodeDescriptor;
80 : int nShapeCount; /* Total number of shapes */
81 : int nMaxDepth; /* Tree depth */
82 : double dfMinX; /* Bounding box of all shapes */
83 : double dfMaxX;
84 : double dfMinY;
85 : double dfMaxY;
86 :
87 : #ifdef DEBUG_IO
88 : int nTotalBytesRead;
89 : #endif
90 : };
91 :
92 : typedef struct
93 : {
94 : SBNSearchHandle hSBN;
95 :
96 : coord bMinX; /* Search bounding box */
97 : coord bMinY;
98 : coord bMaxX;
99 : coord bMaxY;
100 :
101 : int nShapeCount;
102 : int nShapeAlloc;
103 : int *panShapeId; /* 0 based */
104 :
105 : uchar abyBinShape[8 * 100];
106 :
107 : #ifdef DEBUG_IO
108 : int nBytesRead;
109 : #endif
110 : } SearchStruct;
111 :
112 : /************************************************************************/
113 : /* SwapWord() */
114 : /* */
115 : /* Swap a 2, 4 or 8 byte word. */
116 : /************************************************************************/
117 :
118 32 : static void SwapWord( int length, void * wordP )
119 :
120 : {
121 : int i;
122 : uchar temp;
123 :
124 160 : for( i=0; i < length/2; i++ )
125 : {
126 128 : temp = ((uchar *) wordP)[i];
127 128 : ((uchar *)wordP)[i] = ((uchar *) wordP)[length-i-1];
128 128 : ((uchar *) wordP)[length-i-1] = temp;
129 : }
130 32 : }
131 :
132 : /************************************************************************/
133 : /* SBNOpenDiskTree() */
134 : /************************************************************************/
135 :
136 1255 : SBNSearchHandle SBNOpenDiskTree( const char* pszSBNFilename,
137 : SAHooks *psHooks )
138 : {
139 : int i;
140 : SBNSearchHandle hSBN;
141 : uchar abyHeader[108];
142 : int nShapeCount;
143 : int nMaxDepth;
144 : int nMaxNodes;
145 : int nNodeDescSize;
146 : int nNodeDescCount;
147 1255 : uchar* pabyData = NULL;
148 1255 : SBNNodeDescriptor* pasNodeDescriptor = NULL;
149 : uchar abyBinHeader[8];
150 : int nCurNode;
151 : int nNextNonEmptyNode;
152 : int nExpectedBinId;
153 : int bBigEndian;
154 :
155 : /* -------------------------------------------------------------------- */
156 : /* Establish the byte order on this machine. */
157 : /* -------------------------------------------------------------------- */
158 1255 : i = 1;
159 1255 : if( *((unsigned char *) &i) == 1 )
160 1255 : bBigEndian = FALSE;
161 : else
162 0 : bBigEndian = TRUE;
163 :
164 : /* -------------------------------------------------------------------- */
165 : /* Initialize the handle structure. */
166 : /* -------------------------------------------------------------------- */
167 1255 : hSBN = (SBNSearchHandle)
168 : calloc(sizeof(struct SBNSearchInfo),1);
169 :
170 1255 : if (psHooks == NULL)
171 1255 : SASetupDefaultHooks( &(hSBN->sHooks) );
172 : else
173 0 : memcpy( &(hSBN->sHooks), psHooks, sizeof(SAHooks) );
174 :
175 1255 : hSBN->fpSBN = hSBN->sHooks.FOpen(pszSBNFilename, "rb");
176 1255 : if (hSBN->fpSBN == NULL)
177 : {
178 1247 : free(hSBN);
179 1247 : return NULL;
180 : }
181 :
182 : /* -------------------------------------------------------------------- */
183 : /* Check file header signature. */
184 : /* -------------------------------------------------------------------- */
185 72 : if (hSBN->sHooks.FRead(abyHeader, 108, 1, hSBN->fpSBN) != 1 ||
186 8 : abyHeader[0] != 0 ||
187 8 : abyHeader[1] != 0 ||
188 8 : abyHeader[2] != 0x27 ||
189 8 : (abyHeader[3] != 0x0A && abyHeader[3] != 0x0D) ||
190 8 : abyHeader[4] != 0xFF ||
191 8 : abyHeader[5] != 0xFF ||
192 8 : abyHeader[6] != 0xFE ||
193 8 : abyHeader[7] != 0x70)
194 : {
195 0 : hSBN->sHooks.Error( ".sbn file is unreadable, or corrupt." );
196 0 : SBNCloseDiskTree(hSBN);
197 0 : return NULL;
198 : }
199 :
200 : /* -------------------------------------------------------------------- */
201 : /* Read shapes bounding box. */
202 : /* -------------------------------------------------------------------- */
203 :
204 8 : memcpy(&hSBN->dfMinX, abyHeader + 32, 8);
205 8 : memcpy(&hSBN->dfMinY, abyHeader + 40, 8);
206 8 : memcpy(&hSBN->dfMaxX, abyHeader + 48, 8);
207 8 : memcpy(&hSBN->dfMaxY, abyHeader + 56, 8);
208 :
209 8 : if( !bBigEndian )
210 : {
211 8 : SwapWord(8, &hSBN->dfMinX);
212 8 : SwapWord(8, &hSBN->dfMinY);
213 8 : SwapWord(8, &hSBN->dfMaxX);
214 8 : SwapWord(8, &hSBN->dfMaxY);
215 : }
216 :
217 16 : if( hSBN->dfMinX > hSBN->dfMaxX ||
218 8 : hSBN->dfMinY > hSBN->dfMaxY )
219 : {
220 0 : hSBN->sHooks.Error( "Invalid extent in .sbn file." );
221 0 : SBNCloseDiskTree(hSBN);
222 0 : return NULL;
223 : }
224 :
225 : /* -------------------------------------------------------------------- */
226 : /* Read and check number of shapes. */
227 : /* -------------------------------------------------------------------- */
228 8 : nShapeCount = READ_MSB_INT(abyHeader + 28);
229 8 : hSBN->nShapeCount = nShapeCount;
230 8 : if (nShapeCount < 0 || nShapeCount > 256000000 )
231 : {
232 : char szErrorMsg[64];
233 0 : sprintf(szErrorMsg, "Invalid shape count in .sbn : %d", nShapeCount );
234 0 : hSBN->sHooks.Error( szErrorMsg );
235 0 : SBNCloseDiskTree(hSBN);
236 0 : return NULL;
237 : }
238 :
239 : /* -------------------------------------------------------------------- */
240 : /* Compute tree depth. */
241 : /* It is computed such as in average there are not more than 8 */
242 : /* shapes per node. With a minimum depth of 2, and a maximum of 15 */
243 : /* -------------------------------------------------------------------- */
244 8 : nMaxDepth = 2;
245 35 : while( nMaxDepth < 15 && nShapeCount > ((1 << nMaxDepth) - 1) * 8 )
246 19 : nMaxDepth ++;
247 8 : hSBN->nMaxDepth = nMaxDepth;
248 8 : nMaxNodes = (1 << nMaxDepth) - 1;
249 :
250 : /* -------------------------------------------------------------------- */
251 : /* Check that the first bin id is 1. */
252 : /* -------------------------------------------------------------------- */
253 :
254 8 : if( READ_MSB_INT(abyHeader + 100) != 1 )
255 : {
256 0 : hSBN->sHooks.Error( "Unexpected bin id" );
257 0 : SBNCloseDiskTree(hSBN);
258 0 : return NULL;
259 : }
260 :
261 : /* -------------------------------------------------------------------- */
262 : /* Read and check number of node descriptors to be read. */
263 : /* There are at most (2^nMaxDepth) - 1, but all are not necessary */
264 : /* described. Non described nodes are empty. */
265 : /* -------------------------------------------------------------------- */
266 8 : nNodeDescSize = READ_MSB_INT(abyHeader + 104);
267 8 : nNodeDescSize *= 2; /* 16-bit words */
268 :
269 : /* each bin descriptor is made of 2 ints */
270 8 : nNodeDescCount = nNodeDescSize / 8;
271 :
272 8 : if ((nNodeDescSize % 8) != 0 ||
273 : nNodeDescCount < 0 || nNodeDescCount > nMaxNodes )
274 : {
275 : char szErrorMsg[64];
276 0 : sprintf(szErrorMsg,
277 : "Invalid node descriptor size in .sbn : %d", nNodeDescSize );
278 0 : hSBN->sHooks.Error( szErrorMsg );
279 0 : SBNCloseDiskTree(hSBN);
280 0 : return NULL;
281 : }
282 :
283 8 : pabyData = (uchar*) malloc( nNodeDescSize );
284 8 : pasNodeDescriptor = (SBNNodeDescriptor*)
285 : calloc ( nMaxNodes, sizeof(SBNNodeDescriptor) );
286 8 : if (pabyData == NULL || pasNodeDescriptor == NULL)
287 : {
288 0 : free(pabyData);
289 0 : free(pasNodeDescriptor);
290 0 : hSBN->sHooks.Error( "Out of memory error" );
291 0 : SBNCloseDiskTree(hSBN);
292 0 : return NULL;
293 : }
294 :
295 : /* -------------------------------------------------------------------- */
296 : /* Read node descriptors. */
297 : /* -------------------------------------------------------------------- */
298 8 : if (hSBN->sHooks.FRead(pabyData, nNodeDescSize, 1,
299 : hSBN->fpSBN) != 1)
300 : {
301 0 : free(pabyData);
302 0 : free(pasNodeDescriptor);
303 0 : hSBN->sHooks.Error( "Cannot read node descriptors" );
304 0 : SBNCloseDiskTree(hSBN);
305 0 : return NULL;
306 : }
307 :
308 8 : hSBN->pasNodeDescriptor = pasNodeDescriptor;
309 :
310 1042 : for(i = 0; i < nNodeDescCount; i++)
311 : {
312 : /* -------------------------------------------------------------------- */
313 : /* Each node descriptor contains the index of the first bin that */
314 : /* described it, and the number of shapes in this first bin and */
315 : /* the following bins (in the relevant case). */
316 : /* -------------------------------------------------------------------- */
317 1034 : int nBinStart = READ_MSB_INT(pabyData + 8 * i);
318 1034 : int nNodeShapeCount = READ_MSB_INT(pabyData + 8 * i + 4);
319 1034 : pasNodeDescriptor[i].nBinStart = nBinStart > 0 ? nBinStart : 0;
320 1034 : pasNodeDescriptor[i].nShapeCount = nNodeShapeCount;
321 :
322 1034 : if ((nBinStart > 0 && nNodeShapeCount == 0) ||
323 : nNodeShapeCount < 0 || nNodeShapeCount > nShapeCount)
324 : {
325 0 : hSBN->sHooks.Error( "Inconsistant shape count in bin" );
326 0 : SBNCloseDiskTree(hSBN);
327 0 : return NULL;
328 : }
329 : }
330 :
331 8 : free(pabyData);
332 8 : pabyData = NULL;
333 :
334 : /* Locate first non-empty node */
335 8 : nCurNode = 0;
336 20 : while(nCurNode < nMaxNodes && pasNodeDescriptor[nCurNode].nBinStart <= 0)
337 4 : nCurNode ++;
338 :
339 8 : if( nCurNode >= nMaxNodes)
340 : {
341 0 : hSBN->sHooks.Error( "All nodes are empty" );
342 0 : SBNCloseDiskTree(hSBN);
343 0 : return NULL;
344 : }
345 :
346 16 : pasNodeDescriptor[nCurNode].nBinOffset =
347 8 : hSBN->sHooks.FTell(hSBN->fpSBN);
348 :
349 : /* Compute the index of the next non empty node. */
350 8 : nNextNonEmptyNode = nCurNode + 1;
351 27 : while(nNextNonEmptyNode < nMaxNodes &&
352 9 : pasNodeDescriptor[nNextNonEmptyNode].nBinStart <= 0)
353 2 : nNextNonEmptyNode ++;
354 :
355 8 : nExpectedBinId = 1;
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Traverse bins to compute the offset of the first bin of each */
359 : /* node. */
360 : /* Note: we could use the .sbx file to compute the offsets instead.*/
361 : /* -------------------------------------------------------------------- */
362 676 : while( hSBN->sHooks.FRead(abyBinHeader, 8, 1,
363 : hSBN->fpSBN) == 1 )
364 : {
365 : int nBinId;
366 : int nBinSize;
367 :
368 660 : nExpectedBinId ++;
369 :
370 660 : nBinId = READ_MSB_INT(abyBinHeader);
371 660 : nBinSize = READ_MSB_INT(abyBinHeader + 4);
372 660 : nBinSize *= 2; /* 16-bit words */
373 :
374 660 : if( nBinId != nExpectedBinId )
375 : {
376 0 : hSBN->sHooks.Error( "Unexpected bin id" );
377 0 : SBNCloseDiskTree(hSBN);
378 0 : return NULL;
379 : }
380 :
381 : /* Bins are always limited to 100 features */
382 : /* If there are more, then they are located in continuous bins */
383 660 : if( (nBinSize % 8) != 0 || nBinSize <= 0 || nBinSize > 100 * 8)
384 : {
385 0 : hSBN->sHooks.Error( "Unexpected bin size" );
386 0 : SBNCloseDiskTree(hSBN);
387 0 : return NULL;
388 : }
389 :
390 1319 : if( nNextNonEmptyNode < nMaxNodes &&
391 659 : nBinId == pasNodeDescriptor[nNextNonEmptyNode].nBinStart )
392 : {
393 652 : nCurNode = nNextNonEmptyNode;
394 652 : pasNodeDescriptor[nCurNode].nBinOffset =
395 : hSBN->sHooks.FTell(hSBN->fpSBN) - 8;
396 :
397 : /* Compute the index of the next non empty node. */
398 652 : nNextNonEmptyNode = nCurNode + 1;
399 2689 : while(nNextNonEmptyNode < nMaxNodes &&
400 1015 : pasNodeDescriptor[nNextNonEmptyNode].nBinStart <= 0)
401 370 : nNextNonEmptyNode ++;
402 : }
403 :
404 660 : pasNodeDescriptor[nCurNode].nBinCount ++;
405 :
406 : /* Skip shape description */
407 660 : hSBN->sHooks.FSeek(hSBN->fpSBN, nBinSize, SEEK_CUR);
408 : }
409 :
410 8 : return hSBN;
411 : }
412 :
413 : /***********************************************************************/
414 : /* SBNCloseDiskTree() */
415 : /************************************************************************/
416 :
417 14 : void SBNCloseDiskTree( SBNSearchHandle hSBN )
418 : {
419 : int i;
420 : int nMaxNodes;
421 :
422 14 : if (hSBN == NULL)
423 6 : return;
424 :
425 8 : nMaxNodes = (1 << hSBN->nMaxDepth) - 1;
426 1044 : for(i = 0; i < nMaxNodes; i++)
427 : {
428 1036 : if( hSBN->pasNodeDescriptor[i].pabyShapeDesc != NULL )
429 517 : free(hSBN->pasNodeDescriptor[i].pabyShapeDesc);
430 : }
431 :
432 : /* printf("hSBN->nTotalBytesRead = %d\n", hSBN->nTotalBytesRead); */
433 :
434 8 : hSBN->sHooks.FClose(hSBN->fpSBN);
435 8 : free(hSBN->pasNodeDescriptor);
436 8 : free(hSBN);
437 : }
438 :
439 :
440 : /************************************************************************/
441 : /* SfRealloc() */
442 : /* */
443 : /* A realloc cover function that will access a NULL pointer as */
444 : /* a valid input. */
445 : /************************************************************************/
446 :
447 5186 : static void * SfRealloc( void * pMem, int nNewSize )
448 :
449 : {
450 5186 : if( pMem == NULL )
451 5164 : return( (void *) malloc(nNewSize) );
452 : else
453 22 : return( (void *) realloc(pMem,nNewSize) );
454 : }
455 :
456 : /************************************************************************/
457 : /* SBNAddShapeId() */
458 : /************************************************************************/
459 :
460 50297 : static int SBNAddShapeId( SearchStruct* psSearch,
461 : int nShapeId )
462 : {
463 50297 : if (psSearch->nShapeCount == psSearch->nShapeAlloc)
464 : {
465 : int* pNewPtr;
466 :
467 5186 : psSearch->nShapeAlloc =
468 5186 : (int) (((psSearch->nShapeCount + 100) * 5) / 4);
469 5186 : pNewPtr =
470 5186 : (int *) SfRealloc( psSearch->panShapeId,
471 : psSearch->nShapeAlloc * sizeof(int) );
472 5186 : if( pNewPtr == NULL )
473 : {
474 0 : psSearch->hSBN->sHooks.Error( "Out of memory error" );
475 0 : return FALSE;
476 : }
477 5186 : psSearch->panShapeId = pNewPtr;
478 : }
479 :
480 50297 : psSearch->panShapeId[psSearch->nShapeCount] = nShapeId;
481 50297 : psSearch->nShapeCount ++;
482 50297 : return TRUE;
483 : }
484 :
485 : /************************************************************************/
486 : /* SBNSearchDiskInternal() */
487 : /************************************************************************/
488 :
489 : /* Due to the way integer coordinates are rounded, */
490 : /* we can use a strict intersection test, except when the node */
491 : /* bounding box or the search bounding box is degenerated. */
492 : #define SEARCH_BB_INTERSECTS(_bMinX, _bMinY, _bMaxX, _bMaxY) \
493 : (((bSearchMinX < _bMaxX && bSearchMaxX > _bMinX) || \
494 : ((_bMinX == _bMaxX || bSearchMinX == bSearchMaxX) && \
495 : bSearchMinX <= _bMaxX && bSearchMaxX >= _bMinX)) && \
496 : ((bSearchMinY < _bMaxY && bSearchMaxY > _bMinY) || \
497 : ((_bMinY == _bMaxY || bSearchMinY == bSearchMaxY ) && \
498 : bSearchMinY <= _bMaxY && bSearchMaxY >= _bMinY)))
499 :
500 :
501 54834 : static int SBNSearchDiskInternal( SearchStruct* psSearch,
502 : int nDepth,
503 : int nNodeId,
504 : coord bNodeMinX,
505 : coord bNodeMinY,
506 : coord bNodeMaxX,
507 : coord bNodeMaxY )
508 : {
509 : SBNSearchHandle hSBN;
510 : SBNNodeDescriptor* psNode;
511 54834 : coord bSearchMinX = psSearch->bMinX;
512 54834 : coord bSearchMinY = psSearch->bMinY;
513 54834 : coord bSearchMaxX = psSearch->bMaxX;
514 54834 : coord bSearchMaxY = psSearch->bMaxY;
515 :
516 54834 : hSBN = psSearch->hSBN;
517 :
518 54834 : psNode = &(hSBN->pasNodeDescriptor[nNodeId]);
519 :
520 : /* -------------------------------------------------------------------- */
521 : /* Check if this node contains shapes that intersect the search */
522 : /* bounding box. */
523 : /* -------------------------------------------------------------------- */
524 54834 : if ( psNode->bBBoxInit &&
525 : !(SEARCH_BB_INTERSECTS(psNode->bMinX, psNode->bMinY,
526 : psNode->bMaxX, psNode->bMaxY)) )
527 :
528 : {
529 : /* No intersection, then don't try to read the shapes attached */
530 : /* to this node */
531 : }
532 :
533 : /* -------------------------------------------------------------------- */
534 : /* If this node contains shapes that are cached, then read them. */
535 : /* -------------------------------------------------------------------- */
536 27470 : else if (psNode->pabyShapeDesc != NULL)
537 : {
538 : int j;
539 17907 : uchar* pabyShapeDesc = psNode->pabyShapeDesc;
540 :
541 : /* printf("nNodeId = %d, nDepth = %d\n", nNodeId, nDepth); */
542 :
543 468690 : for(j = 0; j < psNode->nShapeCount; j++)
544 : {
545 450783 : coord bMinX = pabyShapeDesc[0];
546 450783 : coord bMinY = pabyShapeDesc[1];
547 450783 : coord bMaxX = pabyShapeDesc[2];
548 450783 : coord bMaxY = pabyShapeDesc[3];
549 :
550 450783 : if( SEARCH_BB_INTERSECTS(bMinX, bMinY, bMaxX, bMaxY) )
551 : {
552 : int nShapeId;
553 :
554 35771 : nShapeId = READ_MSB_INT(pabyShapeDesc + 4);
555 :
556 : /* Caution : we count shape id starting from 0, and not 1 */
557 35771 : nShapeId --;
558 :
559 : /*printf("shape=%d, minx=%d, miny=%d, maxx=%d, maxy=%d\n",
560 : nShapeId, bMinX, bMinY, bMaxX, bMaxY);*/
561 :
562 35771 : if( !SBNAddShapeId( psSearch, nShapeId ) )
563 0 : return FALSE;
564 : }
565 :
566 450783 : pabyShapeDesc += 8;
567 : }
568 : }
569 :
570 : /* -------------------------------------------------------------------- */
571 : /* If the node has attached shapes (that are not (yet) cached), */
572 : /* then retrieve them from disk. */
573 : /* -------------------------------------------------------------------- */
574 :
575 9563 : else if (psNode->nBinCount > 0)
576 : {
577 : uchar abyBinHeader[8];
578 : int nBinSize, nShapes;
579 3947 : int nShapeCountAcc = 0;
580 : int i, j;
581 :
582 : /* printf("nNodeId = %d, nDepth = %d\n", nNodeId, nDepth); */
583 :
584 3947 : hSBN->sHooks.FSeek(hSBN->fpSBN, psNode->nBinOffset, SEEK_SET);
585 :
586 3947 : if (nDepth < CACHED_DEPTH_LIMIT)
587 517 : psNode->pabyShapeDesc = (uchar*) malloc(psNode->nShapeCount * 8);
588 :
589 7894 : for(i = 0; i < psNode->nBinCount; i++)
590 : {
591 : uchar* pabyBinShape;
592 :
593 : #ifdef DEBUG_IO
594 : psSearch->nBytesRead += 8;
595 : #endif
596 3947 : if( hSBN->sHooks.FRead(abyBinHeader, 8, 1,
597 : hSBN->fpSBN) != 1)
598 : {
599 0 : hSBN->sHooks.Error( "I/O error" );
600 0 : free(psNode->pabyShapeDesc);
601 0 : psNode->pabyShapeDesc = NULL;
602 0 : return FALSE;
603 : }
604 :
605 3947 : if ( READ_MSB_INT(abyBinHeader + 0) != psNode->nBinStart + i )
606 : {
607 0 : hSBN->sHooks.Error( "Unexpected bin id" );
608 0 : free(psNode->pabyShapeDesc);
609 0 : psNode->pabyShapeDesc = NULL;
610 0 : return FALSE;
611 : }
612 :
613 3947 : nBinSize = READ_MSB_INT(abyBinHeader + 4);
614 3947 : nBinSize *= 2; /* 16-bit words */
615 :
616 3947 : nShapes = nBinSize / 8;
617 :
618 : /* Bins are always limited to 100 features */
619 3947 : if( (nBinSize % 8) != 0 || nShapes <= 0 || nShapes > 100)
620 : {
621 0 : hSBN->sHooks.Error( "Unexpected bin size" );
622 0 : free(psNode->pabyShapeDesc);
623 0 : psNode->pabyShapeDesc = NULL;
624 0 : return FALSE;
625 : }
626 :
627 3947 : if( nShapeCountAcc + nShapes > psNode->nShapeCount)
628 : {
629 0 : free(psNode->pabyShapeDesc);
630 0 : psNode->pabyShapeDesc = NULL;
631 0 : hSBN->sHooks.Error( "Inconsistant shape count for bin" );
632 0 : return FALSE;
633 : }
634 :
635 4464 : if (nDepth < CACHED_DEPTH_LIMIT && psNode->pabyShapeDesc != NULL)
636 : {
637 517 : pabyBinShape = psNode->pabyShapeDesc + nShapeCountAcc * 8;
638 : }
639 : else
640 : {
641 3430 : pabyBinShape = psSearch->abyBinShape;
642 : }
643 :
644 : #ifdef DEBUG_IO
645 : psSearch->nBytesRead += nBinSize;
646 : #endif
647 3947 : if (hSBN->sHooks.FRead(pabyBinShape, nBinSize, 1,
648 : hSBN->fpSBN) != 1)
649 : {
650 0 : hSBN->sHooks.Error( "I/O error" );
651 0 : free(psNode->pabyShapeDesc);
652 0 : psNode->pabyShapeDesc = NULL;
653 0 : return FALSE;
654 : }
655 :
656 3947 : nShapeCountAcc += nShapes;
657 :
658 3947 : if (i == 0 && !psNode->bBBoxInit)
659 : {
660 660 : psNode->bMinX = pabyBinShape[0];
661 660 : psNode->bMinY = pabyBinShape[1];
662 660 : psNode->bMaxX = pabyBinShape[2];
663 660 : psNode->bMaxY = pabyBinShape[3];
664 : }
665 :
666 63333 : for(j = 0; j < nShapes; j++)
667 : {
668 59386 : coord bMinX = pabyBinShape[0];
669 59386 : coord bMinY = pabyBinShape[1];
670 59386 : coord bMaxX = pabyBinShape[2];
671 59386 : coord bMaxY = pabyBinShape[3];
672 :
673 59386 : if( !psNode->bBBoxInit )
674 : {
675 : #ifdef sanity_checks
676 : /* -------------------------------------------------------------------- */
677 : /* Those tests only check that the shape bounding box in the bin */
678 : /* are consistant (self-consistant and consistant with the node */
679 : /* they are attached to). They are optional however (as far as */
680 : /* the safety of runtime is concerned at least). */
681 : /* -------------------------------------------------------------------- */
682 :
683 : if( !(((bMinX < bMaxX ||
684 : (bMinX == 0 && bMaxX == 0) ||
685 : (bMinX == 255 && bMaxX == 255))) &&
686 : ((bMinY < bMaxY ||
687 : (bMinY == 0 && bMaxY == 0) ||
688 : (bMinY == 255 && bMaxY == 255)))) ||
689 : bMaxX < bNodeMinX || bMaxY < bNodeMinY ||
690 : bMinX > bNodeMaxX || bMinY > bNodeMaxY )
691 : {
692 : /*printf("shape %d %d %d %d\n", bMinX, bMinY, bMaxX, bMaxY);
693 : printf("node %d %d %d %d\n", bNodeMinX, bNodeMinY, bNodeMaxX, bNodeMaxY);*/
694 : hSBN->sHooks.Error(
695 : "Invalid shape bounding box in bin" );
696 : free(psNode->pabyShapeDesc);
697 : psNode->pabyShapeDesc = NULL;
698 : return FALSE;
699 : }
700 : #endif
701 5146 : if (bMinX < psNode->bMinX) psNode->bMinX = bMinX;
702 5146 : if (bMinY < psNode->bMinY) psNode->bMinY = bMinY;
703 5146 : if (bMaxX > psNode->bMaxX) psNode->bMaxX = bMaxX;
704 5146 : if (bMaxY > psNode->bMaxY) psNode->bMaxY = bMaxY;
705 : }
706 :
707 59386 : if( SEARCH_BB_INTERSECTS(bMinX, bMinY, bMaxX, bMaxY) )
708 : {
709 : int nShapeId;
710 :
711 14526 : nShapeId = READ_MSB_INT(pabyBinShape + 4);
712 :
713 : /* Caution : we count shape id starting from 0, and not 1 */
714 14526 : nShapeId --;
715 :
716 : /*printf("shape=%d, minx=%d, miny=%d, maxx=%d, maxy=%d\n",
717 : nShapeId, bMinX, bMinY, bMaxX, bMaxY);*/
718 :
719 14526 : if( !SBNAddShapeId( psSearch, nShapeId ) )
720 0 : return FALSE;
721 : }
722 :
723 59386 : pabyBinShape += 8;
724 : }
725 : }
726 :
727 3947 : if( nShapeCountAcc != psNode->nShapeCount)
728 : {
729 0 : free(psNode->pabyShapeDesc);
730 0 : psNode->pabyShapeDesc = NULL;
731 0 : hSBN->sHooks.Error( "Inconsistant shape count for bin" );
732 0 : return FALSE;
733 : }
734 :
735 3947 : psNode->bBBoxInit = TRUE;
736 : }
737 :
738 : /* -------------------------------------------------------------------- */
739 : /* Look up in child nodes. */
740 : /* -------------------------------------------------------------------- */
741 54834 : if( nDepth + 1 < hSBN->nMaxDepth )
742 : {
743 45563 : nNodeId = nNodeId * 2 + 1;
744 :
745 45563 : if( (nDepth % 2) == 0 ) /* x split */
746 : {
747 23608 : coord bMid = (coord) (1 + ((int)bNodeMinX + bNodeMaxX) / 2);
748 36836 : if( bSearchMinX <= bMid - 1 &&
749 13228 : !SBNSearchDiskInternal( psSearch, nDepth + 1, nNodeId + 1,
750 : bNodeMinX, bNodeMinY,
751 : bMid - 1, bNodeMaxY ) )
752 : {
753 0 : return FALSE;
754 : }
755 36252 : if( bSearchMaxX >= bMid &&
756 12644 : !SBNSearchDiskInternal( psSearch, nDepth + 1, nNodeId,
757 : bMid, bNodeMinY,
758 : bNodeMaxX, bNodeMaxY ) )
759 : {
760 0 : return FALSE;
761 : }
762 : }
763 : else /* y split */
764 : {
765 21955 : coord bMid = (coord) (1 + ((int)bNodeMinY + bNodeMaxY) / 2);
766 34515 : if( bSearchMinY <= bMid - 1 &&
767 12560 : !SBNSearchDiskInternal( psSearch, nDepth + 1, nNodeId + 1,
768 : bNodeMinX, bNodeMinY,
769 : bNodeMaxX, bMid - 1 ) )
770 : {
771 0 : return FALSE;
772 : }
773 33193 : if( bSearchMaxY >= bMid &&
774 11238 : !SBNSearchDiskInternal( psSearch, nDepth + 1, nNodeId,
775 : bNodeMinX, bMid,
776 : bNodeMaxX, bNodeMaxY ) )
777 : {
778 0 : return FALSE;
779 : }
780 : }
781 : }
782 :
783 54834 : return TRUE;
784 : }
785 :
786 : /************************************************************************/
787 : /* compare_ints() */
788 : /************************************************************************/
789 :
790 : /* helper for qsort */
791 : static int
792 152961 : compare_ints( const void * a, const void * b)
793 : {
794 152961 : return (*(int*)a) - (*(int*)b);
795 : }
796 :
797 : /************************************************************************/
798 : /* SBNSearchDiskTree() */
799 : /************************************************************************/
800 :
801 5164 : int* SBNSearchDiskTree( SBNSearchHandle hSBN,
802 : double *padfBoundsMin, double *padfBoundsMax,
803 : int *pnShapeCount )
804 : {
805 : double dfMinX, dfMinY, dfMaxX, dfMaxY;
806 : double dfDiskXExtent, dfDiskYExtent;
807 : int bMinX, bMinY, bMaxX, bMaxY;
808 :
809 5164 : *pnShapeCount = 0;
810 :
811 5164 : dfMinX = padfBoundsMin[0];
812 5164 : dfMinY = padfBoundsMin[1];
813 5164 : dfMaxX = padfBoundsMax[0];
814 5164 : dfMaxY = padfBoundsMax[1];
815 :
816 5164 : if( dfMinX > dfMaxX || dfMinY > dfMaxY )
817 0 : return NULL;
818 :
819 15492 : if( dfMaxX < hSBN->dfMinX || dfMaxY < hSBN->dfMinY ||
820 10328 : dfMinX > hSBN->dfMaxX || dfMinY > hSBN->dfMaxY )
821 0 : return NULL;
822 :
823 : /* -------------------------------------------------------------------- */
824 : /* Compute the search coordinates in [0,255]x[0,255] coord. space */
825 : /* -------------------------------------------------------------------- */
826 5164 : dfDiskXExtent = hSBN->dfMaxX - hSBN->dfMinX;
827 5164 : dfDiskYExtent = hSBN->dfMaxY - hSBN->dfMinY;
828 :
829 5164 : if ( dfDiskXExtent == 0.0 )
830 : {
831 0 : bMinX = 0;
832 0 : bMaxX = 255;
833 : }
834 : else
835 : {
836 5164 : if( dfMinX < hSBN->dfMinX )
837 0 : bMinX = 0;
838 : else
839 : {
840 5164 : double dfMinX_255 = (dfMinX - hSBN->dfMinX)
841 5164 : / dfDiskXExtent * 255.0;
842 5164 : bMinX = (int)floor(dfMinX_255 - 0.005);
843 5164 : if( bMinX < 0 ) bMinX = 0;
844 : }
845 :
846 5164 : if( dfMaxX > hSBN->dfMaxX )
847 3 : bMaxX = 255;
848 : else
849 : {
850 5161 : double dfMaxX_255 = (dfMaxX - hSBN->dfMinX)
851 5161 : / dfDiskXExtent * 255.0;
852 5161 : bMaxX = (int)ceil(dfMaxX_255 + 0.005);
853 5161 : if( bMaxX > 255 ) bMaxX = 255;
854 : }
855 : }
856 :
857 5164 : if ( dfDiskYExtent == 0.0 )
858 : {
859 0 : bMinY = 0;
860 0 : bMaxY = 255;
861 : }
862 : else
863 : {
864 5164 : if( dfMinY < hSBN->dfMinY )
865 0 : bMinY = 0;
866 : else
867 : {
868 5164 : double dfMinY_255 = (dfMinY - hSBN->dfMinY)
869 5164 : / dfDiskYExtent * 255.0;
870 5164 : bMinY = (int)floor(dfMinY_255 - 0.005);
871 5164 : if( bMinY < 0 ) bMinY = 0;
872 : }
873 :
874 5164 : if( dfMaxY > hSBN->dfMaxY )
875 0 : bMaxY = 255;
876 : else
877 : {
878 5164 : double dfMaxY_255 = (dfMaxY - hSBN->dfMinY)
879 5164 : / dfDiskYExtent * 255.0;
880 5164 : bMaxY = (int)ceil(dfMaxY_255 + 0.005);
881 5164 : if( bMaxY > 255 ) bMaxY = 255;
882 : }
883 : }
884 :
885 : /* -------------------------------------------------------------------- */
886 : /* Run the search. */
887 : /* -------------------------------------------------------------------- */
888 :
889 5164 : return SBNSearchDiskTreeInteger(hSBN,
890 : bMinX, bMinY, bMaxX, bMaxY,
891 : pnShapeCount);
892 : }
893 :
894 : /************************************************************************/
895 : /* SBNSearchDiskTreeInteger() */
896 : /************************************************************************/
897 :
898 5164 : int* SBNSearchDiskTreeInteger( SBNSearchHandle hSBN,
899 : int bMinX, int bMinY, int bMaxX, int bMaxY,
900 : int *pnShapeCount )
901 : {
902 : SearchStruct sSearch;
903 : int bRet;
904 :
905 5164 : *pnShapeCount = 0;
906 :
907 5164 : if( bMinX > bMaxX || bMinY > bMaxY )
908 0 : return NULL;
909 :
910 5164 : if( bMaxX < 0 || bMaxY < 0 || bMinX > 255 || bMinX > 255 )
911 0 : return NULL;
912 :
913 : /* -------------------------------------------------------------------- */
914 : /* Run the search. */
915 : /* -------------------------------------------------------------------- */
916 5164 : sSearch.hSBN = hSBN;
917 5164 : sSearch.bMinX = (coord) (bMinX >= 0 ? bMinX : 0);
918 5164 : sSearch.bMinY = (coord) (bMinY >= 0 ? bMinY : 0);
919 5164 : sSearch.bMaxX = (coord) (bMaxX <= 255 ? bMaxX : 255);
920 5164 : sSearch.bMaxY = (coord) (bMaxY <= 255 ? bMaxY : 255);
921 5164 : sSearch.nShapeCount = 0;
922 5164 : sSearch.nShapeAlloc = 0;
923 5164 : sSearch.panShapeId = NULL;
924 : #ifdef DEBUG_IO
925 : sSearch.nBytesRead = 0;
926 : #endif
927 :
928 5164 : bRet = SBNSearchDiskInternal(&sSearch, 0, 0, 0, 0, 255, 255);
929 :
930 : #ifdef DEBUG_IO
931 : hSBN->nTotalBytesRead += sSearch.nBytesRead;
932 : /* printf("nBytesRead = %d\n", sSearch.nBytesRead); */
933 : #endif
934 :
935 5164 : if( !bRet )
936 : {
937 0 : if( sSearch.panShapeId != NULL )
938 0 : free( sSearch.panShapeId );
939 0 : *pnShapeCount = 0;
940 0 : return NULL;
941 : }
942 :
943 5164 : *pnShapeCount = sSearch.nShapeCount;
944 :
945 : /* -------------------------------------------------------------------- */
946 : /* Sort the id array */
947 : /* -------------------------------------------------------------------- */
948 5164 : qsort(sSearch.panShapeId, *pnShapeCount, sizeof(int), compare_ints);
949 :
950 : /* To distinguish between empty intersection from error case */
951 5164 : if( sSearch.panShapeId == NULL )
952 0 : sSearch.panShapeId = (int*) calloc(1, sizeof(int));
953 :
954 5164 : return sSearch.panShapeId;
955 : }
956 :
957 : /************************************************************************/
958 : /* SBNSearchFreeIds() */
959 : /************************************************************************/
960 :
961 0 : void SBNSearchFreeIds( int* panShapeId )
962 : {
963 0 : free( panShapeId );
964 0 : }
|
