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