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