Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OGR
4 : * Purpose: WKB geometry related methods
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2022, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_error.h"
14 : #include "ogr_wkb.h"
15 : #include "ogr_core.h"
16 : #include "ogr_geometry.h"
17 : #include "ogr_p.h"
18 :
19 : #include <algorithm>
20 : #include <cmath>
21 : #include <climits>
22 : #include <limits>
23 :
24 : #include <algorithm>
25 : #include <limits>
26 :
27 : #define USE_FAST_FLOAT
28 : #ifdef USE_FAST_FLOAT
29 : #include "include_fast_float.h"
30 : #endif
31 :
32 : /************************************************************************/
33 : /* OGRWKBNeedSwap() */
34 : /************************************************************************/
35 :
36 200 : static inline bool OGRWKBNeedSwap(GByte b)
37 : {
38 : #if CPL_IS_LSB
39 200 : const bool bNeedSwap = b == 0;
40 : #else
41 : const bool bNeedSwap = b == 1;
42 : #endif
43 200 : return bNeedSwap;
44 : }
45 :
46 : /************************************************************************/
47 : /* OGRWKBReadUInt32() */
48 : /************************************************************************/
49 :
50 223 : static inline uint32_t OGRWKBReadUInt32(const GByte *pabyWkb, bool bNeedSwap)
51 : {
52 : uint32_t nVal;
53 223 : memcpy(&nVal, pabyWkb, sizeof(nVal));
54 223 : if (bNeedSwap)
55 0 : CPL_SWAP32PTR(&nVal);
56 223 : return nVal;
57 : }
58 :
59 : /************************************************************************/
60 : /* OGRWKBReadFloat64() */
61 : /************************************************************************/
62 :
63 100 : static inline double OGRWKBReadFloat64(const GByte *pabyWkb, bool bNeedSwap)
64 : {
65 : double dfVal;
66 100 : memcpy(&dfVal, pabyWkb, sizeof(dfVal));
67 100 : if (bNeedSwap)
68 0 : CPL_SWAP64PTR(&dfVal);
69 100 : return dfVal;
70 : }
71 :
72 : /************************************************************************/
73 : /* OGRWKBRingGetArea() */
74 : /************************************************************************/
75 :
76 12 : static bool OGRWKBRingGetArea(const GByte *&pabyWkb, size_t &nWKBSize, int nDim,
77 : bool bNeedSwap, double &dfArea)
78 : {
79 12 : const uint32_t nPoints = OGRWKBReadUInt32(pabyWkb, bNeedSwap);
80 12 : if (nPoints >= 4 &&
81 12 : (nWKBSize - sizeof(uint32_t)) / (nDim * sizeof(double)) >= nPoints)
82 : {
83 12 : nWKBSize -= sizeof(uint32_t) + nDim * sizeof(double);
84 12 : pabyWkb += sizeof(uint32_t);
85 : // Computation according to Green's Theorem
86 : // Cf OGRSimpleCurve::get_LinearArea()
87 12 : double x_m1 = OGRWKBReadFloat64(pabyWkb, bNeedSwap);
88 12 : double y_m1 = OGRWKBReadFloat64(pabyWkb + sizeof(double), bNeedSwap);
89 12 : double y_m2 = y_m1;
90 12 : dfArea = 0;
91 12 : pabyWkb += nDim * sizeof(double);
92 50 : for (uint32_t i = 1; i < nPoints; ++i)
93 : {
94 38 : const double x = OGRWKBReadFloat64(pabyWkb, bNeedSwap);
95 : const double y =
96 38 : OGRWKBReadFloat64(pabyWkb + sizeof(double), bNeedSwap);
97 38 : pabyWkb += nDim * sizeof(double);
98 38 : dfArea += x_m1 * (y - y_m2);
99 38 : y_m2 = y_m1;
100 38 : x_m1 = x;
101 38 : y_m1 = y;
102 : }
103 12 : dfArea += x_m1 * (y_m1 - y_m2);
104 12 : dfArea = 0.5 * std::fabs(dfArea);
105 12 : return true;
106 : }
107 0 : return false;
108 : }
109 :
110 : /************************************************************************/
111 : /* OGRWKBGetGeomType() */
112 : /************************************************************************/
113 :
114 195 : bool OGRWKBGetGeomType(const GByte *pabyWkb, size_t nWKBSize, bool &bNeedSwap,
115 : uint32_t &nType)
116 : {
117 195 : if (nWKBSize >= 5)
118 : {
119 195 : bNeedSwap = OGRWKBNeedSwap(pabyWkb[0]);
120 195 : nType = OGRWKBReadUInt32(pabyWkb + 1, bNeedSwap);
121 195 : return true;
122 : }
123 0 : return false;
124 : }
125 :
126 : /************************************************************************/
127 : /* OGRWKBPolygonGetArea() */
128 : /************************************************************************/
129 :
130 11 : bool OGRWKBPolygonGetArea(const GByte *&pabyWkb, size_t &nWKBSize,
131 : double &dfArea)
132 : {
133 : bool bNeedSwap;
134 : uint32_t nType;
135 11 : if (nWKBSize < 9 || !OGRWKBGetGeomType(pabyWkb, nWKBSize, bNeedSwap, nType))
136 0 : return false;
137 :
138 11 : int nDims = 2;
139 11 : if (nType == wkbPolygon)
140 : {
141 : // do nothing
142 : }
143 6 : else if (nType == wkbPolygon + 1000 || // wkbPolygonZ
144 4 : nType == wkbPolygon25D || nType == wkbPolygonM)
145 : {
146 4 : nDims = 3;
147 : }
148 2 : else if (nType == wkbPolygonZM)
149 : {
150 2 : nDims = 4;
151 : }
152 : else
153 : {
154 0 : return false;
155 : }
156 :
157 11 : const uint32_t nRings = OGRWKBReadUInt32(pabyWkb + 5, bNeedSwap);
158 11 : if ((nWKBSize - 9) / sizeof(uint32_t) >= nRings)
159 : {
160 11 : pabyWkb += 9;
161 11 : nWKBSize -= 9;
162 11 : dfArea = 0;
163 11 : if (nRings > 0)
164 : {
165 11 : if (!OGRWKBRingGetArea(pabyWkb, nWKBSize, nDims, bNeedSwap, dfArea))
166 0 : return false;
167 12 : for (uint32_t i = 1; i < nRings; ++i)
168 : {
169 : double dfRingArea;
170 1 : if (!OGRWKBRingGetArea(pabyWkb, nWKBSize, nDims, bNeedSwap,
171 : dfRingArea))
172 0 : return false;
173 1 : dfArea -= dfRingArea;
174 : }
175 : }
176 11 : return true;
177 : }
178 0 : return false;
179 : }
180 :
181 : /************************************************************************/
182 : /* OGRWKBMultiPolygonGetArea() */
183 : /************************************************************************/
184 :
185 5 : bool OGRWKBMultiPolygonGetArea(const GByte *&pabyWkb, size_t &nWKBSize,
186 : double &dfArea)
187 : {
188 5 : if (nWKBSize < 9)
189 0 : return false;
190 :
191 5 : const bool bNeedSwap = OGRWKBNeedSwap(pabyWkb[0]);
192 5 : const uint32_t nPolys = OGRWKBReadUInt32(pabyWkb + 5, bNeedSwap);
193 5 : if ((nWKBSize - 9) / 9 >= nPolys)
194 : {
195 5 : pabyWkb += 9;
196 5 : nWKBSize -= 9;
197 5 : dfArea = 0;
198 11 : for (uint32_t i = 0; i < nPolys; ++i)
199 : {
200 : double dfPolyArea;
201 6 : if (!OGRWKBPolygonGetArea(pabyWkb, nWKBSize, dfPolyArea))
202 0 : return false;
203 6 : dfArea += dfPolyArea;
204 : }
205 5 : return true;
206 : }
207 0 : return false;
208 : }
209 :
210 : /************************************************************************/
211 : /* WKBFromEWKB() */
212 : /************************************************************************/
213 :
214 1429 : const GByte *WKBFromEWKB(GByte *pabyEWKB, size_t nEWKBSize, size_t &nWKBSizeOut,
215 : int *pnSRIDOut)
216 : {
217 1429 : if (nEWKBSize < 5U)
218 : {
219 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid EWKB content : %u bytes",
220 : static_cast<unsigned>(nEWKBSize));
221 0 : return nullptr;
222 : }
223 :
224 1429 : const GByte *pabyWKB = pabyEWKB;
225 :
226 : /* -------------------------------------------------------------------- */
227 : /* PostGIS EWKB format includes an SRID, but this won't be */
228 : /* understood by OGR, so if the SRID flag is set, we remove the */
229 : /* SRID (bytes at offset 5 to 8). */
230 : /* -------------------------------------------------------------------- */
231 1429 : if (nEWKBSize > 9 &&
232 1422 : ((pabyEWKB[0] == 0 /* big endian */ && (pabyEWKB[1] & 0x20)) ||
233 1422 : (pabyEWKB[0] != 0 /* little endian */ && (pabyEWKB[4] & 0x20))))
234 : {
235 514 : if (pnSRIDOut)
236 : {
237 0 : memcpy(pnSRIDOut, pabyEWKB + 5, 4);
238 0 : const OGRwkbByteOrder eByteOrder =
239 0 : (pabyEWKB[0] == 0 ? wkbXDR : wkbNDR);
240 0 : if (OGR_SWAP(eByteOrder))
241 0 : *pnSRIDOut = CPL_SWAP32(*pnSRIDOut);
242 : }
243 :
244 : // Drop the SRID flag
245 514 : if (pabyEWKB[0] == 0)
246 0 : pabyEWKB[1] &= (~0x20);
247 : else
248 514 : pabyEWKB[4] &= (~0x20);
249 :
250 : // Move 5 first bytes of EWKB 4 bytes later to create regular WKB
251 514 : memmove(pabyEWKB + 4, pabyEWKB, 5);
252 514 : memset(pabyEWKB, 0, 4);
253 : // and make pabyWKB point to that
254 514 : pabyWKB += 4;
255 514 : nWKBSizeOut = nEWKBSize - 4;
256 : }
257 : else
258 : {
259 915 : if (pnSRIDOut)
260 : {
261 0 : *pnSRIDOut = INT_MIN;
262 : }
263 915 : nWKBSizeOut = nEWKBSize;
264 : }
265 :
266 1429 : return pabyWKB;
267 : }
268 :
269 : /************************************************************************/
270 : /* OGRWKBReadUInt32AtOffset() */
271 : /************************************************************************/
272 :
273 8435 : static uint32_t OGRWKBReadUInt32AtOffset(const uint8_t *data,
274 : OGRwkbByteOrder eByteOrder,
275 : size_t &iOffset)
276 : {
277 : uint32_t v;
278 8435 : memcpy(&v, data + iOffset, sizeof(v));
279 8435 : iOffset += sizeof(v);
280 8435 : if (OGR_SWAP(eByteOrder))
281 : {
282 1666 : CPL_SWAP32PTR(&v);
283 : }
284 8435 : return v;
285 : }
286 :
287 : /************************************************************************/
288 : /* ReadWKBPointSequence() */
289 : /************************************************************************/
290 :
291 : template <bool INCLUDE_Z, typename EnvelopeType>
292 844 : static bool ReadWKBPointSequence(const uint8_t *data, size_t size,
293 : OGRwkbByteOrder eByteOrder, int nDim,
294 : bool bHasZ, size_t &iOffset,
295 : EnvelopeType &sEnvelope)
296 : {
297 : const uint32_t nPoints =
298 844 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
299 844 : if (nPoints > (size - iOffset) / (nDim * sizeof(double)))
300 0 : return false;
301 844 : double dfX = 0;
302 844 : double dfY = 0;
303 844 : [[maybe_unused]] double dfZ = 0;
304 14492 : for (uint32_t j = 0; j < nPoints; j++)
305 : {
306 13648 : memcpy(&dfX, data + iOffset, sizeof(double));
307 13648 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
308 : if constexpr (INCLUDE_Z)
309 : {
310 52 : if (bHasZ)
311 34 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
312 : sizeof(double));
313 : }
314 13648 : iOffset += nDim * sizeof(double);
315 13648 : if (OGR_SWAP(eByteOrder))
316 : {
317 560 : CPL_SWAP64PTR(&dfX);
318 560 : CPL_SWAP64PTR(&dfY);
319 : if constexpr (INCLUDE_Z)
320 : {
321 0 : CPL_SWAP64PTR(&dfZ);
322 : }
323 : }
324 13648 : sEnvelope.MinX = std::min(sEnvelope.MinX, dfX);
325 13648 : sEnvelope.MinY = std::min(sEnvelope.MinY, dfY);
326 13648 : sEnvelope.MaxX = std::max(sEnvelope.MaxX, dfX);
327 13648 : sEnvelope.MaxY = std::max(sEnvelope.MaxY, dfY);
328 : if constexpr (INCLUDE_Z)
329 : {
330 52 : if (bHasZ)
331 : {
332 34 : sEnvelope.MinZ = std::min(sEnvelope.MinZ, dfZ);
333 34 : sEnvelope.MaxZ = std::max(sEnvelope.MaxZ, dfZ);
334 : }
335 : }
336 : }
337 844 : return true;
338 : }
339 :
340 : /************************************************************************/
341 : /* ReadWKBRingSequence() */
342 : /************************************************************************/
343 :
344 : template <bool INCLUDE_Z, typename EnvelopeType>
345 592 : static bool ReadWKBRingSequence(const uint8_t *data, size_t size,
346 : OGRwkbByteOrder eByteOrder, int nDim,
347 : bool bHasZ, size_t &iOffset,
348 : EnvelopeType &sEnvelope)
349 : {
350 592 : const uint32_t nRings = OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
351 592 : if (nRings > (size - iOffset) / sizeof(uint32_t))
352 0 : return false;
353 1242 : for (uint32_t i = 0; i < nRings; i++)
354 : {
355 650 : if (iOffset + sizeof(uint32_t) > size)
356 0 : return false;
357 650 : if (!ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
358 : bHasZ, iOffset, sEnvelope))
359 0 : return false;
360 : }
361 592 : return true;
362 : }
363 :
364 : /************************************************************************/
365 : /* OGRWKBGetBoundingBox() */
366 : /************************************************************************/
367 :
368 : constexpr uint32_t WKB_PREFIX_SIZE = 1 + sizeof(uint32_t);
369 : constexpr uint32_t MIN_WKB_SIZE = WKB_PREFIX_SIZE + sizeof(uint32_t);
370 :
371 : template <bool INCLUDE_Z, typename EnvelopeType>
372 958550 : static bool OGRWKBGetBoundingBox(const uint8_t *data, size_t size,
373 : size_t &iOffset, EnvelopeType &sEnvelope,
374 : int nRec)
375 : {
376 958550 : if (size - iOffset < MIN_WKB_SIZE)
377 0 : return false;
378 958550 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffset]);
379 958550 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
380 0 : return false;
381 958550 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
382 :
383 958550 : OGRwkbGeometryType eGeometryType = wkbUnknown;
384 958550 : OGRReadWKBGeometryType(data + iOffset, wkbVariantIso, &eGeometryType);
385 958550 : iOffset += 5;
386 958550 : const auto eFlatType = wkbFlatten(eGeometryType);
387 958550 : const bool bHasZ = CPL_TO_BOOL(OGR_GT_HasZ(eGeometryType));
388 958550 : const int nDim = 2 + (bHasZ ? 1 : 0) + (OGR_GT_HasM(eGeometryType) ? 1 : 0);
389 :
390 958550 : if (eFlatType == wkbPoint)
391 : {
392 957986 : if (size - iOffset < nDim * sizeof(double))
393 0 : return false;
394 957986 : double dfX = 0;
395 957986 : double dfY = 0;
396 957986 : [[maybe_unused]] double dfZ = 0;
397 957986 : memcpy(&dfX, data + iOffset, sizeof(double));
398 957986 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
399 : if constexpr (INCLUDE_Z)
400 : {
401 11 : if (bHasZ)
402 4 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
403 : sizeof(double));
404 : }
405 957986 : iOffset += nDim * sizeof(double);
406 957986 : if (OGR_SWAP(eByteOrder))
407 : {
408 40 : CPL_SWAP64PTR(&dfX);
409 40 : CPL_SWAP64PTR(&dfY);
410 : if constexpr (INCLUDE_Z)
411 : {
412 0 : CPL_SWAP64PTR(&dfZ);
413 : }
414 : }
415 957986 : if (std::isnan(dfX))
416 : {
417 : // Point empty
418 6 : sEnvelope = EnvelopeType();
419 : }
420 : else
421 : {
422 957980 : sEnvelope.MinX = dfX;
423 957980 : sEnvelope.MinY = dfY;
424 957980 : sEnvelope.MaxX = dfX;
425 957980 : sEnvelope.MaxY = dfY;
426 : if constexpr (INCLUDE_Z)
427 : {
428 10 : if (bHasZ)
429 : {
430 4 : sEnvelope.MinZ = dfZ;
431 4 : sEnvelope.MaxZ = dfZ;
432 : }
433 : }
434 : }
435 957986 : return true;
436 : }
437 :
438 564 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
439 : {
440 106 : sEnvelope = EnvelopeType();
441 :
442 106 : return ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
443 106 : bHasZ, iOffset, sEnvelope);
444 : }
445 :
446 458 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
447 : {
448 182 : sEnvelope = EnvelopeType();
449 :
450 182 : return ReadWKBRingSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
451 182 : bHasZ, iOffset, sEnvelope);
452 : }
453 :
454 276 : if (eFlatType == wkbMultiPoint)
455 : {
456 57 : sEnvelope = EnvelopeType();
457 :
458 57 : uint32_t nParts = OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
459 57 : if (nParts >
460 57 : (size - iOffset) / (WKB_PREFIX_SIZE + nDim * sizeof(double)))
461 0 : return false;
462 57 : double dfX = 0;
463 57 : double dfY = 0;
464 57 : [[maybe_unused]] double dfZ = 0;
465 170 : for (uint32_t k = 0; k < nParts; k++)
466 : {
467 113 : iOffset += WKB_PREFIX_SIZE;
468 113 : memcpy(&dfX, data + iOffset, sizeof(double));
469 113 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
470 : if constexpr (INCLUDE_Z)
471 : {
472 5 : if (bHasZ)
473 3 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
474 : sizeof(double));
475 : }
476 113 : iOffset += nDim * sizeof(double);
477 113 : if (OGR_SWAP(eByteOrder))
478 : {
479 40 : CPL_SWAP64PTR(&dfX);
480 40 : CPL_SWAP64PTR(&dfY);
481 : if constexpr (INCLUDE_Z)
482 : {
483 0 : CPL_SWAP64PTR(&dfZ);
484 : }
485 : }
486 113 : sEnvelope.MinX = std::min(sEnvelope.MinX, dfX);
487 113 : sEnvelope.MinY = std::min(sEnvelope.MinY, dfY);
488 113 : sEnvelope.MaxX = std::max(sEnvelope.MaxX, dfX);
489 113 : sEnvelope.MaxY = std::max(sEnvelope.MaxY, dfY);
490 : if constexpr (INCLUDE_Z)
491 : {
492 5 : if (bHasZ)
493 : {
494 3 : sEnvelope.MinZ = std::min(sEnvelope.MinZ, dfZ);
495 3 : sEnvelope.MaxZ = std::max(sEnvelope.MaxZ, dfZ);
496 : }
497 : }
498 : }
499 57 : return true;
500 : }
501 :
502 219 : if (eFlatType == wkbMultiLineString)
503 : {
504 59 : sEnvelope = EnvelopeType();
505 :
506 : const uint32_t nParts =
507 59 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
508 59 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
509 0 : return false;
510 147 : for (uint32_t k = 0; k < nParts; k++)
511 : {
512 88 : if (iOffset + MIN_WKB_SIZE > size)
513 0 : return false;
514 88 : iOffset += WKB_PREFIX_SIZE;
515 88 : if (!ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
516 : bHasZ, iOffset, sEnvelope))
517 0 : return false;
518 : }
519 59 : return true;
520 : }
521 :
522 160 : if (eFlatType == wkbMultiPolygon)
523 : {
524 82 : sEnvelope = EnvelopeType();
525 :
526 : const uint32_t nParts =
527 82 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
528 82 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
529 0 : return false;
530 492 : for (uint32_t k = 0; k < nParts; k++)
531 : {
532 410 : if (iOffset + MIN_WKB_SIZE > size)
533 0 : return false;
534 410 : CPLAssert(data[iOffset] == eByteOrder);
535 410 : iOffset += WKB_PREFIX_SIZE;
536 410 : if (!ReadWKBRingSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
537 : bHasZ, iOffset, sEnvelope))
538 0 : return false;
539 : }
540 82 : return true;
541 : }
542 :
543 78 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve ||
544 4 : eFlatType == wkbCurvePolygon || eFlatType == wkbMultiCurve ||
545 2 : eFlatType == wkbMultiSurface || eFlatType == wkbPolyhedralSurface ||
546 : eFlatType == wkbTIN)
547 : {
548 78 : if (nRec == 128)
549 0 : return false;
550 78 : sEnvelope = EnvelopeType();
551 :
552 : const uint32_t nParts =
553 78 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
554 78 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
555 0 : return false;
556 78 : EnvelopeType sEnvelopeSubGeom;
557 176 : for (uint32_t k = 0; k < nParts; k++)
558 : {
559 98 : if (!OGRWKBGetBoundingBox<INCLUDE_Z>(data, size, iOffset,
560 : sEnvelopeSubGeom, nRec + 1))
561 0 : return false;
562 98 : sEnvelope.Merge(sEnvelopeSubGeom);
563 : }
564 78 : return true;
565 : }
566 :
567 0 : return false;
568 : }
569 :
570 : /************************************************************************/
571 : /* OGRWKBGetBoundingBox() */
572 : /************************************************************************/
573 :
574 958417 : bool OGRWKBGetBoundingBox(const GByte *pabyWkb, size_t nWKBSize,
575 : OGREnvelope &sEnvelope)
576 : {
577 958417 : size_t iOffset = 0;
578 958417 : return OGRWKBGetBoundingBox<false>(pabyWkb, nWKBSize, iOffset, sEnvelope,
579 1916830 : 0);
580 : }
581 :
582 : /************************************************************************/
583 : /* OGRWKBGetBoundingBox() */
584 : /************************************************************************/
585 :
586 35 : bool OGRWKBGetBoundingBox(const GByte *pabyWkb, size_t nWKBSize,
587 : OGREnvelope3D &sEnvelope)
588 : {
589 35 : size_t iOffset = 0;
590 70 : return OGRWKBGetBoundingBox<true>(pabyWkb, nWKBSize, iOffset, sEnvelope, 0);
591 : }
592 :
593 : /************************************************************************/
594 : /* OGRWKBIntersectsPointSequencePessimistic() */
595 : /************************************************************************/
596 :
597 190 : static bool OGRWKBIntersectsPointSequencePessimistic(
598 : const uint8_t *data, const size_t size, const OGRwkbByteOrder eByteOrder,
599 : const int nDim, size_t &iOffsetInOut, const OGREnvelope &sEnvelope,
600 : bool &bErrorOut)
601 : {
602 : const uint32_t nPoints =
603 190 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
604 190 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
605 : {
606 12 : bErrorOut = true;
607 12 : return false;
608 : }
609 :
610 178 : double dfX = 0;
611 178 : double dfY = 0;
612 1211 : for (uint32_t j = 0; j < nPoints; j++)
613 : {
614 1149 : memcpy(&dfX, data + iOffsetInOut, sizeof(double));
615 1149 : memcpy(&dfY, data + iOffsetInOut + sizeof(double), sizeof(double));
616 1149 : iOffsetInOut += nDim * sizeof(double);
617 1149 : if (OGR_SWAP(eByteOrder))
618 : {
619 3 : CPL_SWAP64PTR(&dfX);
620 3 : CPL_SWAP64PTR(&dfY);
621 : }
622 1149 : if (dfX >= sEnvelope.MinX && dfY >= sEnvelope.MinY &&
623 709 : dfX <= sEnvelope.MaxX && dfY <= sEnvelope.MaxY)
624 : {
625 116 : return true;
626 : }
627 : }
628 :
629 62 : return false;
630 : }
631 :
632 : /************************************************************************/
633 : /* OGRWKBIntersectsRingSequencePessimistic() */
634 : /************************************************************************/
635 :
636 140 : static bool OGRWKBIntersectsRingSequencePessimistic(
637 : const uint8_t *data, const size_t size, const OGRwkbByteOrder eByteOrder,
638 : const int nDim, size_t &iOffsetInOut, const OGREnvelope &sEnvelope,
639 : bool &bErrorOut)
640 : {
641 : const uint32_t nRings =
642 140 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
643 140 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
644 : {
645 16 : bErrorOut = true;
646 16 : return false;
647 : }
648 124 : if (nRings == 0)
649 2 : return false;
650 122 : if (iOffsetInOut + sizeof(uint32_t) > size)
651 : {
652 0 : bErrorOut = true;
653 0 : return false;
654 : }
655 122 : if (OGRWKBIntersectsPointSequencePessimistic(
656 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut))
657 : {
658 89 : return true;
659 : }
660 33 : if (bErrorOut)
661 0 : return false;
662 :
663 : // skip inner rings
664 39 : for (uint32_t i = 1; i < nRings; ++i)
665 : {
666 6 : if (iOffsetInOut + sizeof(uint32_t) > size)
667 : {
668 0 : bErrorOut = true;
669 0 : return false;
670 : }
671 : const uint32_t nPoints =
672 6 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
673 6 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
674 : {
675 0 : bErrorOut = true;
676 0 : return false;
677 : }
678 6 : iOffsetInOut += sizeof(double) * nPoints * nDim;
679 : }
680 33 : return false;
681 : }
682 :
683 : /************************************************************************/
684 : /* OGRWKBIntersectsPessimistic() */
685 : /************************************************************************/
686 :
687 23428 : static bool OGRWKBIntersectsPessimistic(const GByte *data, const size_t size,
688 : size_t &iOffsetInOut,
689 : const OGREnvelope &sEnvelope,
690 : const int nRec, bool &bErrorOut)
691 : {
692 23428 : if (size - iOffsetInOut < MIN_WKB_SIZE)
693 : {
694 0 : bErrorOut = true;
695 0 : return false;
696 : }
697 23428 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
698 23428 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
699 : {
700 0 : bErrorOut = true;
701 0 : return false;
702 : }
703 23428 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
704 :
705 23428 : OGRwkbGeometryType eGeometryType = wkbUnknown;
706 23428 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
707 23428 : iOffsetInOut += 5;
708 23428 : const auto eFlatType = wkbFlatten(eGeometryType);
709 23428 : const int nDim = 2 + (OGR_GT_HasZ(eGeometryType) ? 1 : 0) +
710 23428 : (OGR_GT_HasM(eGeometryType) ? 1 : 0);
711 :
712 23428 : if (eFlatType == wkbPoint)
713 : {
714 23070 : if (size - iOffsetInOut < nDim * sizeof(double))
715 8 : return false;
716 23062 : double dfX = 0;
717 23062 : double dfY = 0;
718 23062 : memcpy(&dfX, data + iOffsetInOut, sizeof(double));
719 23062 : memcpy(&dfY, data + iOffsetInOut + sizeof(double), sizeof(double));
720 23062 : iOffsetInOut += nDim * sizeof(double);
721 23062 : if (OGR_SWAP(eByteOrder))
722 : {
723 0 : CPL_SWAP64PTR(&dfX);
724 0 : CPL_SWAP64PTR(&dfY);
725 : }
726 23062 : if (std::isnan(dfX))
727 : {
728 1 : return false;
729 : }
730 : else
731 : {
732 23054 : return dfX >= sEnvelope.MinX && dfX <= sEnvelope.MaxX &&
733 46115 : dfY >= sEnvelope.MinY && dfY <= sEnvelope.MaxY;
734 : }
735 : }
736 :
737 358 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
738 : {
739 68 : return OGRWKBIntersectsPointSequencePessimistic(
740 68 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut);
741 : }
742 :
743 290 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
744 : {
745 140 : return OGRWKBIntersectsRingSequencePessimistic(
746 140 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut);
747 : }
748 :
749 150 : if (eFlatType == wkbMultiPoint || eFlatType == wkbMultiLineString ||
750 81 : eFlatType == wkbMultiPolygon || eFlatType == wkbGeometryCollection ||
751 65 : eFlatType == wkbCompoundCurve || eFlatType == wkbCurvePolygon ||
752 39 : eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface ||
753 13 : eFlatType == wkbPolyhedralSurface || eFlatType == wkbTIN)
754 : {
755 150 : if (nRec == 128)
756 : {
757 0 : bErrorOut = true;
758 0 : return false;
759 : }
760 : const uint32_t nParts =
761 150 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
762 150 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
763 : {
764 80 : bErrorOut = true;
765 80 : return false;
766 : }
767 115 : for (uint32_t k = 0; k < nParts; k++)
768 : {
769 79 : if (OGRWKBIntersectsPessimistic(data, size, iOffsetInOut, sEnvelope,
770 : nRec + 1, bErrorOut))
771 : {
772 34 : return true;
773 : }
774 45 : else if (bErrorOut)
775 : {
776 0 : return false;
777 : }
778 : }
779 36 : return false;
780 : }
781 :
782 0 : bErrorOut = true;
783 0 : return false;
784 : }
785 :
786 : /************************************************************************/
787 : /* OGRWKBIntersectsPessimistic() */
788 : /************************************************************************/
789 :
790 : /* Returns whether the geometry (pabyWkb, nWKBSize) intersects, for sure,
791 : * the passed envelope.
792 : * When it returns true, the geometry intersects the envelope.
793 : * When it returns false, the geometry may or may not intersect the envelope.
794 : */
795 23349 : bool OGRWKBIntersectsPessimistic(const GByte *pabyWkb, size_t nWKBSize,
796 : const OGREnvelope &sEnvelope)
797 : {
798 23349 : size_t iOffsetInOut = 0;
799 23349 : bool bErrorOut = false;
800 23349 : bool bRet = OGRWKBIntersectsPessimistic(pabyWkb, nWKBSize, iOffsetInOut,
801 : sEnvelope, 0, bErrorOut);
802 23349 : if (!bRet && !bErrorOut)
803 : {
804 : // The following assert only holds if there is no trailing data
805 : // after the WKB
806 : // CPLAssert(iOffsetInOut == nWKBSize);
807 : }
808 23349 : return bRet;
809 : }
810 :
811 : /************************************************************************/
812 : /* epsilonEqual() */
813 : /************************************************************************/
814 :
815 68 : static inline bool epsilonEqual(double a, double b, double eps)
816 : {
817 68 : return ::fabs(a - b) < eps;
818 : }
819 :
820 : /************************************************************************/
821 : /* OGRWKBIsClockwiseRing() */
822 : /************************************************************************/
823 :
824 318 : static inline double GetX(const GByte *data, uint32_t i, int nDim,
825 : bool bNeedSwap)
826 : {
827 : double dfX;
828 318 : memcpy(&dfX, data + static_cast<size_t>(i) * nDim * sizeof(double),
829 : sizeof(double));
830 318 : if (bNeedSwap)
831 0 : CPL_SWAP64PTR(&dfX);
832 318 : return dfX;
833 : }
834 :
835 597 : static inline double GetY(const GByte *data, uint32_t i, int nDim,
836 : bool bNeedSwap)
837 : {
838 : double dfY;
839 597 : memcpy(&dfY, data + (static_cast<size_t>(i) * nDim + 1) * sizeof(double),
840 : sizeof(double));
841 597 : if (bNeedSwap)
842 0 : CPL_SWAP64PTR(&dfY);
843 597 : return dfY;
844 : }
845 :
846 29 : static bool OGRWKBIsClockwiseRing(const GByte *data, const uint32_t nPoints,
847 : const int nDim, const bool bNeedSwap)
848 : {
849 29 : constexpr double EPSILON = 1.0E-5;
850 :
851 : // WARNING: keep in sync OGRLineString::isClockwise(),
852 : // OGRCurve::isClockwise() and OGRWKBIsClockwiseRing()
853 :
854 29 : bool bUseFallback = false;
855 :
856 : // Find the lowest rightmost vertex.
857 29 : uint32_t v = 0; // Used after for.
858 29 : double vX = GetX(data, v, nDim, bNeedSwap);
859 29 : double vY = GetY(data, v, nDim, bNeedSwap);
860 503 : for (uint32_t i = 1; i < nPoints - 1; i++)
861 : {
862 : // => v < end.
863 474 : const double y = GetY(data, i, nDim, bNeedSwap);
864 474 : if (y < vY)
865 : {
866 155 : v = i;
867 155 : vX = GetX(data, i, nDim, bNeedSwap);
868 155 : vY = y;
869 155 : bUseFallback = false;
870 : }
871 319 : else if (y == vY)
872 : {
873 2 : const double x = GetX(data, i, nDim, bNeedSwap);
874 2 : if (x > vX)
875 : {
876 0 : v = i;
877 0 : vX = x;
878 : // vY = y;
879 0 : bUseFallback = false;
880 : }
881 2 : else if (x == vX)
882 : {
883 : // Two vertex with same coordinates are the lowest rightmost
884 : // vertex. Cannot use that point as the pivot (#5342).
885 2 : bUseFallback = true;
886 : }
887 : }
888 : }
889 :
890 : // Previous.
891 29 : uint32_t next = (v == 0) ? nPoints - 2 : v - 1;
892 33 : if (epsilonEqual(GetX(data, next, nDim, bNeedSwap), vX, EPSILON) &&
893 4 : epsilonEqual(GetY(data, next, nDim, bNeedSwap), vY, EPSILON))
894 : {
895 : // Don't try to be too clever by retrying with a next point.
896 : // This can lead to false results as in the case of #3356.
897 0 : bUseFallback = true;
898 : }
899 :
900 29 : const double dx0 = GetX(data, next, nDim, bNeedSwap) - vX;
901 29 : const double dy0 = GetY(data, next, nDim, bNeedSwap) - vY;
902 :
903 : // Following.
904 29 : next = v + 1;
905 29 : if (next >= nPoints - 1)
906 : {
907 2 : next = 0;
908 : }
909 :
910 35 : if (epsilonEqual(GetX(data, next, nDim, bNeedSwap), vX, EPSILON) &&
911 6 : epsilonEqual(GetY(data, next, nDim, bNeedSwap), vY, EPSILON))
912 : {
913 : // Don't try to be too clever by retrying with a next point.
914 : // This can lead to false results as in the case of #3356.
915 2 : bUseFallback = true;
916 : }
917 :
918 29 : const double dx1 = GetX(data, next, nDim, bNeedSwap) - vX;
919 29 : const double dy1 = GetY(data, next, nDim, bNeedSwap) - vY;
920 :
921 29 : const double crossproduct = dx1 * dy0 - dx0 * dy1;
922 :
923 29 : if (!bUseFallback)
924 : {
925 27 : if (crossproduct > 0) // CCW
926 23 : return false;
927 4 : else if (crossproduct < 0) // CW
928 4 : return true;
929 : }
930 :
931 : // This is a degenerate case: the extent of the polygon is less than EPSILON
932 : // or 2 nearly identical points were found.
933 : // Try with Green Formula as a fallback, but this is not a guarantee
934 : // as we'll probably be affected by numerical instabilities.
935 :
936 2 : double dfSum = GetX(data, 0, nDim, bNeedSwap) *
937 2 : (GetY(data, 1, nDim, bNeedSwap) -
938 2 : GetY(data, nPoints - 1, nDim, bNeedSwap));
939 :
940 12 : for (uint32_t i = 1; i < nPoints - 1; i++)
941 : {
942 10 : dfSum += GetX(data, i, nDim, bNeedSwap) *
943 20 : (GetY(data, i + 1, nDim, bNeedSwap) -
944 10 : GetY(data, i - 1, nDim, bNeedSwap));
945 : }
946 :
947 2 : dfSum += GetX(data, nPoints - 1, nDim, bNeedSwap) *
948 2 : (GetY(data, 0, nDim, bNeedSwap) -
949 2 : GetX(data, nPoints - 2, nDim, bNeedSwap));
950 :
951 2 : return dfSum < 0;
952 : }
953 :
954 : /************************************************************************/
955 : /* OGRWKBFixupCounterClockWiseExternalRing() */
956 : /************************************************************************/
957 :
958 55 : static bool OGRWKBFixupCounterClockWiseExternalRingInternal(
959 : GByte *data, size_t size, size_t &iOffsetInOut, const int nRec)
960 : {
961 55 : if (size - iOffsetInOut < MIN_WKB_SIZE)
962 : {
963 0 : return false;
964 : }
965 55 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
966 55 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
967 : {
968 0 : return false;
969 : }
970 55 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
971 :
972 55 : OGRwkbGeometryType eGeometryType = wkbUnknown;
973 55 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
974 55 : iOffsetInOut += 5;
975 55 : const auto eFlatType = wkbFlatten(eGeometryType);
976 55 : const int nDim = 2 + (OGR_GT_HasZ(eGeometryType) ? 1 : 0) +
977 55 : (OGR_GT_HasM(eGeometryType) ? 1 : 0);
978 :
979 55 : if (eFlatType == wkbPolygon)
980 : {
981 : const uint32_t nRings =
982 26 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
983 26 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
984 : {
985 0 : return false;
986 : }
987 55 : for (uint32_t iRing = 0; iRing < nRings; ++iRing)
988 : {
989 29 : if (iOffsetInOut + sizeof(uint32_t) > size)
990 0 : return false;
991 : const uint32_t nPoints =
992 29 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
993 29 : const size_t sizeOfPoint = nDim * sizeof(double);
994 29 : if (nPoints > (size - iOffsetInOut) / sizeOfPoint)
995 : {
996 0 : return false;
997 : }
998 :
999 29 : if (nPoints >= 4)
1000 : {
1001 29 : const bool bIsClockwiseRing = OGRWKBIsClockwiseRing(
1002 29 : data + iOffsetInOut, nPoints, nDim, OGR_SWAP(eByteOrder));
1003 29 : if ((bIsClockwiseRing && iRing == 0) ||
1004 25 : (!bIsClockwiseRing && iRing > 0))
1005 : {
1006 : GByte abyTmp[4 * sizeof(double)];
1007 19 : for (uint32_t i = 0; i < nPoints / 2; ++i)
1008 : {
1009 13 : GByte *pBegin = data + iOffsetInOut + i * sizeOfPoint;
1010 13 : GByte *pEnd = data + iOffsetInOut +
1011 13 : (nPoints - 1 - i) * sizeOfPoint;
1012 13 : memcpy(abyTmp, pBegin, sizeOfPoint);
1013 13 : memcpy(pBegin, pEnd, sizeOfPoint);
1014 13 : memcpy(pEnd, abyTmp, sizeOfPoint);
1015 : }
1016 : }
1017 : }
1018 :
1019 29 : iOffsetInOut += nPoints * sizeOfPoint;
1020 : }
1021 : }
1022 :
1023 55 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbMultiPolygon ||
1024 : eFlatType == wkbMultiSurface)
1025 : {
1026 6 : if (nRec == 128)
1027 : {
1028 0 : return false;
1029 : }
1030 : const uint32_t nParts =
1031 6 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1032 6 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
1033 : {
1034 0 : return false;
1035 : }
1036 11 : for (uint32_t k = 0; k < nParts; k++)
1037 : {
1038 5 : if (!OGRWKBFixupCounterClockWiseExternalRingInternal(
1039 : data, size, iOffsetInOut, nRec))
1040 : {
1041 0 : return false;
1042 : }
1043 : }
1044 : }
1045 :
1046 55 : return true;
1047 : }
1048 :
1049 : /** Modifies the geometry such that exterior rings of polygons are
1050 : * counter-clockwise oriented and inner rings clockwise oriented.
1051 : */
1052 50 : void OGRWKBFixupCounterClockWiseExternalRing(GByte *pabyWkb, size_t nWKBSize)
1053 : {
1054 50 : size_t iOffsetInOut = 0;
1055 50 : OGRWKBFixupCounterClockWiseExternalRingInternal(
1056 : pabyWkb, nWKBSize, iOffsetInOut, /* nRec = */ 0);
1057 50 : }
1058 :
1059 : /************************************************************************/
1060 : /* OGRWKBPointUpdater() */
1061 : /************************************************************************/
1062 :
1063 : OGRWKBPointUpdater::OGRWKBPointUpdater() = default;
1064 :
1065 : /************************************************************************/
1066 : /* OGRWKBIntersectsPointSequencePessimistic() */
1067 : /************************************************************************/
1068 :
1069 3019 : static bool OGRWKBUpdatePointsSequence(uint8_t *data, const size_t size,
1070 : OGRWKBPointUpdater &oUpdater,
1071 : const OGRwkbByteOrder eByteOrder,
1072 : const int nDim, const bool bHasZ,
1073 : const bool bHasM, size_t &iOffsetInOut)
1074 : {
1075 : const uint32_t nPoints =
1076 3019 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1077 3019 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
1078 : {
1079 1497 : return false;
1080 : }
1081 1522 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1082 5812 : for (uint32_t j = 0; j < nPoints; j++)
1083 : {
1084 4292 : void *pdfX = data + iOffsetInOut;
1085 4292 : void *pdfY = data + iOffsetInOut + sizeof(double);
1086 4292 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1087 4292 : void *pdfM =
1088 4292 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1089 : : nullptr;
1090 4292 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1091 2 : return false;
1092 :
1093 4290 : iOffsetInOut += nDim * sizeof(double);
1094 : }
1095 :
1096 1520 : return true;
1097 : }
1098 :
1099 : /************************************************************************/
1100 : /* OGRWKBVisitRingSequence() */
1101 : /************************************************************************/
1102 :
1103 1766 : static bool OGRWKBVisitRingSequence(uint8_t *data, const size_t size,
1104 : OGRWKBPointUpdater &oUpdater,
1105 : const OGRwkbByteOrder eByteOrder,
1106 : const int nDim, const bool bHasZ,
1107 : const bool bHasM, size_t &iOffsetInOut)
1108 : {
1109 : const uint32_t nRings =
1110 1766 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1111 1766 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
1112 : {
1113 150 : return false;
1114 : }
1115 :
1116 2514 : for (uint32_t i = 0; i < nRings; ++i)
1117 : {
1118 1717 : if (iOffsetInOut + sizeof(uint32_t) > size)
1119 : {
1120 4 : return false;
1121 : }
1122 1713 : if (!OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder, nDim,
1123 : bHasZ, bHasM, iOffsetInOut))
1124 : {
1125 815 : return false;
1126 : }
1127 : }
1128 797 : return true;
1129 : }
1130 :
1131 : /************************************************************************/
1132 : /* OGRWKBUpdatePoints() */
1133 : /************************************************************************/
1134 :
1135 16254 : static bool OGRWKBUpdatePoints(uint8_t *data, const size_t size,
1136 : OGRWKBPointUpdater &oUpdater,
1137 : size_t &iOffsetInOut, const int nRec)
1138 : {
1139 16254 : if (size - iOffsetInOut < MIN_WKB_SIZE)
1140 : {
1141 514 : return false;
1142 : }
1143 15740 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
1144 15740 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
1145 : {
1146 68 : return false;
1147 : }
1148 15672 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
1149 :
1150 15672 : OGRwkbGeometryType eGeometryType = wkbUnknown;
1151 15672 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
1152 15689 : iOffsetInOut += 5;
1153 15689 : const auto eFlatType = wkbFlatten(eGeometryType);
1154 :
1155 15689 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve ||
1156 15218 : eFlatType == wkbCurvePolygon || eFlatType == wkbMultiPoint ||
1157 14908 : eFlatType == wkbMultiLineString || eFlatType == wkbMultiPolygon ||
1158 14721 : eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface ||
1159 14460 : eFlatType == wkbPolyhedralSurface || eFlatType == wkbTIN)
1160 : {
1161 1392 : if (nRec == 128)
1162 1 : return false;
1163 :
1164 : const uint32_t nParts =
1165 1391 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1166 1391 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
1167 : {
1168 176 : return false;
1169 : }
1170 1856 : for (uint32_t k = 0; k < nParts; k++)
1171 : {
1172 1358 : if (!OGRWKBUpdatePoints(data, size, oUpdater, iOffsetInOut,
1173 : nRec + 1))
1174 717 : return false;
1175 : }
1176 498 : return true;
1177 : }
1178 :
1179 14297 : const bool bHasZ = OGR_GT_HasZ(eGeometryType);
1180 14295 : const bool bHasM = OGR_GT_HasM(eGeometryType);
1181 14288 : const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0);
1182 :
1183 14288 : if (eFlatType == wkbPoint)
1184 : {
1185 11014 : if (size - iOffsetInOut < nDim * sizeof(double))
1186 376 : return false;
1187 10638 : void *pdfX = data + iOffsetInOut;
1188 10638 : void *pdfY = data + iOffsetInOut + sizeof(double);
1189 10638 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1190 10638 : void *pdfM =
1191 10638 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1192 : : nullptr;
1193 10638 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1194 10638 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1195 1 : return false;
1196 10637 : iOffsetInOut += nDim * sizeof(double);
1197 10637 : return true;
1198 : }
1199 :
1200 3274 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
1201 : {
1202 1305 : return OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder,
1203 1306 : nDim, bHasZ, bHasM, iOffsetInOut);
1204 : }
1205 :
1206 1969 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
1207 : {
1208 1766 : return OGRWKBVisitRingSequence(data, size, oUpdater, eByteOrder, nDim,
1209 1766 : bHasZ, bHasM, iOffsetInOut);
1210 : }
1211 :
1212 203 : CPLDebug("OGR", "Unknown WKB geometry type");
1213 203 : return false;
1214 : }
1215 :
1216 : /** Visit all points of a WKB geometry to update them.
1217 : */
1218 14896 : bool OGRWKBUpdatePoints(GByte *pabyWkb, size_t nWKBSize,
1219 : OGRWKBPointUpdater &oUpdater)
1220 : {
1221 14896 : size_t iOffsetInOut = 0;
1222 14896 : return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater, iOffsetInOut,
1223 29806 : /* nRec = */ 0);
1224 : }
1225 :
1226 : /************************************************************************/
1227 : /* OGRWKBTransformCache::clear() */
1228 : /************************************************************************/
1229 :
1230 : #ifdef OGR_WKB_TRANSFORM_ALL_AT_ONCE
1231 : void OGRWKBTransformCache::clear()
1232 : {
1233 : abNeedSwap.clear();
1234 : abIsEmpty.clear();
1235 : apdfX.clear();
1236 : apdfY.clear();
1237 : apdfZ.clear();
1238 : apdfM.clear();
1239 : adfX.clear();
1240 : adfY.clear();
1241 : adfZ.clear();
1242 : adfM.clear();
1243 : anErrorCodes.clear();
1244 : }
1245 : #endif
1246 :
1247 : /************************************************************************/
1248 : /* OGRWKBTransform() */
1249 : /************************************************************************/
1250 :
1251 : /** Visit all points of a WKB geometry to transform them.
1252 : */
1253 14903 : bool OGRWKBTransform(GByte *pabyWkb, size_t nWKBSize,
1254 : OGRCoordinateTransformation *poCT,
1255 : [[maybe_unused]] OGRWKBTransformCache &oCache,
1256 : OGREnvelope3D &sEnvelope)
1257 : {
1258 : #ifndef OGR_WKB_TRANSFORM_ALL_AT_ONCE
1259 : struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater
1260 : {
1261 : OGRCoordinateTransformation *m_poCT;
1262 : OGREnvelope3D &m_sEnvelope;
1263 :
1264 14897 : explicit OGRWKBPointUpdaterReproj(OGRCoordinateTransformation *poCTIn,
1265 : OGREnvelope3D &sEnvelopeIn)
1266 14897 : : m_poCT(poCTIn), m_sEnvelope(sEnvelopeIn)
1267 : {
1268 14896 : }
1269 :
1270 14931 : bool update(bool bNeedSwap, void *x, void *y, void *z,
1271 : void * /* m */) override
1272 : {
1273 : double dfX, dfY, dfZ;
1274 14931 : memcpy(&dfX, x, sizeof(double));
1275 14931 : memcpy(&dfY, y, sizeof(double));
1276 14931 : if (bNeedSwap)
1277 : {
1278 1116 : CPL_SWAP64PTR(&dfX);
1279 1116 : CPL_SWAP64PTR(&dfY);
1280 : }
1281 14931 : if (!(std::isnan(dfX) && std::isnan(dfY)))
1282 : {
1283 14729 : if (z)
1284 : {
1285 1166 : memcpy(&dfZ, z, sizeof(double));
1286 1166 : if (bNeedSwap)
1287 : {
1288 646 : CPL_SWAP64PTR(&dfZ);
1289 : }
1290 : }
1291 : else
1292 13563 : dfZ = 0;
1293 14729 : int nErrorCode = 0;
1294 14729 : m_poCT->TransformWithErrorCodes(1, &dfX, &dfY, &dfZ, nullptr,
1295 14729 : &nErrorCode);
1296 14734 : if (nErrorCode)
1297 3 : return false;
1298 14731 : m_sEnvelope.Merge(dfX, dfY, dfZ);
1299 14729 : if (bNeedSwap)
1300 : {
1301 1016 : CPL_SWAP64PTR(&dfX);
1302 1016 : CPL_SWAP64PTR(&dfY);
1303 1016 : CPL_SWAP64PTR(&dfZ);
1304 : }
1305 14727 : memcpy(x, &dfX, sizeof(double));
1306 14727 : memcpy(y, &dfY, sizeof(double));
1307 14727 : if (z)
1308 1166 : memcpy(z, &dfZ, sizeof(double));
1309 : }
1310 14927 : return true;
1311 : }
1312 :
1313 : private:
1314 : OGRWKBPointUpdaterReproj(const OGRWKBPointUpdaterReproj &) = delete;
1315 : OGRWKBPointUpdaterReproj &
1316 : operator=(const OGRWKBPointUpdaterReproj &) = delete;
1317 : };
1318 :
1319 14903 : sEnvelope = OGREnvelope3D();
1320 27098 : OGRWKBPointUpdaterReproj oUpdater(poCT, sEnvelope);
1321 29793 : return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater);
1322 :
1323 : #else
1324 : struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater
1325 : {
1326 : OGRWKBTransformCache &m_oCache;
1327 :
1328 : explicit OGRWKBPointUpdaterReproj(OGRWKBTransformCache &oCacheIn)
1329 : : m_oCache(oCacheIn)
1330 : {
1331 : }
1332 :
1333 : bool update(bool bNeedSwap, void *x, void *y, void *z,
1334 : void * /* m */) override
1335 : {
1336 : m_oCache.abNeedSwap.push_back(bNeedSwap);
1337 : m_oCache.apdfX.push_back(x);
1338 : m_oCache.apdfY.push_back(y);
1339 : m_oCache.apdfZ.push_back(z);
1340 : return true;
1341 : }
1342 : };
1343 :
1344 : oCache.clear();
1345 : OGRWKBPointUpdaterReproj oUpdater(oCache);
1346 : if (!OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater))
1347 : return false;
1348 :
1349 : oCache.adfX.resize(oCache.apdfX.size());
1350 : oCache.adfY.resize(oCache.apdfX.size());
1351 : oCache.adfZ.resize(oCache.apdfX.size());
1352 :
1353 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1354 : {
1355 : memcpy(&oCache.adfX[i], oCache.apdfX[i], sizeof(double));
1356 : memcpy(&oCache.adfY[i], oCache.apdfY[i], sizeof(double));
1357 : if (oCache.apdfZ[i])
1358 : memcpy(&oCache.adfZ[i], oCache.apdfZ[i], sizeof(double));
1359 : if (oCache.abNeedSwap[i])
1360 : {
1361 : CPL_SWAP64PTR(&oCache.adfX[i]);
1362 : CPL_SWAP64PTR(&oCache.adfY[i]);
1363 : CPL_SWAP64PTR(&oCache.adfZ[i]);
1364 : }
1365 : oCache.abIsEmpty.push_back(std::isnan(oCache.adfX[i]) &&
1366 : std::isnan(oCache.adfY[i]));
1367 : }
1368 :
1369 : oCache.anErrorCodes.resize(oCache.apdfX.size());
1370 : poCT->TransformWithErrorCodes(static_cast<int>(oCache.apdfX.size()),
1371 : oCache.adfX.data(), oCache.adfY.data(),
1372 : oCache.adfZ.data(), nullptr,
1373 : oCache.anErrorCodes.data());
1374 :
1375 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1376 : {
1377 : if (!oCache.abIsEmpty[i] && oCache.anErrorCodes[i])
1378 : return false;
1379 : }
1380 :
1381 : sEnvelope = OGREnvelope3D();
1382 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1383 : {
1384 : if (oCache.abIsEmpty[i])
1385 : {
1386 : oCache.adfX[i] = std::numeric_limits<double>::quiet_NaN();
1387 : oCache.adfY[i] = std::numeric_limits<double>::quiet_NaN();
1388 : oCache.adfZ[i] = std::numeric_limits<double>::quiet_NaN();
1389 : }
1390 : else
1391 : {
1392 : sEnvelope.Merge(oCache.adfX[i], oCache.adfY[i], oCache.adfZ[i]);
1393 : }
1394 : if (oCache.abNeedSwap[i])
1395 : {
1396 : CPL_SWAP64PTR(&oCache.adfX[i]);
1397 : CPL_SWAP64PTR(&oCache.adfY[i]);
1398 : CPL_SWAP64PTR(&oCache.adfZ[i]);
1399 : }
1400 : memcpy(oCache.apdfX[i], &oCache.adfX[i], sizeof(double));
1401 : memcpy(oCache.apdfY[i], &oCache.adfY[i], sizeof(double));
1402 : if (oCache.apdfZ[i])
1403 : memcpy(oCache.apdfZ[i], &oCache.adfZ[i], sizeof(double));
1404 : }
1405 :
1406 : return true;
1407 : #endif
1408 : }
1409 :
1410 : /************************************************************************/
1411 : /* OGRAppendBuffer() */
1412 : /************************************************************************/
1413 :
1414 : OGRAppendBuffer::OGRAppendBuffer() = default;
1415 :
1416 : /************************************************************************/
1417 : /* ~OGRAppendBuffer() */
1418 : /************************************************************************/
1419 :
1420 : OGRAppendBuffer::~OGRAppendBuffer() = default;
1421 :
1422 : /************************************************************************/
1423 : /* OGRWKTToWKBTranslator() */
1424 : /************************************************************************/
1425 :
1426 18 : OGRWKTToWKBTranslator::OGRWKTToWKBTranslator(OGRAppendBuffer &oAppendBuffer)
1427 18 : : m_oAppendBuffer(oAppendBuffer)
1428 : {
1429 : #ifndef USE_FAST_FLOAT
1430 : // Test if current locale decimal separator is decimal point
1431 : char szTest[10];
1432 : snprintf(szTest, sizeof(szTest), "%f", 1.5);
1433 : m_bCanUseStrtod = strchr(szTest, '.') != nullptr;
1434 : #endif
1435 18 : CPL_IGNORE_RET_VAL(m_bCanUseStrtod);
1436 18 : }
1437 :
1438 : /************************************************************************/
1439 : /* TranslateWKT() */
1440 : /************************************************************************/
1441 :
1442 138 : size_t OGRWKTToWKBTranslator::TranslateWKT(void *pabyWKTStart, size_t nLength,
1443 : bool bCanAlterByteAfter)
1444 : {
1445 138 : const char *pszPtrStart = static_cast<const char *>(pabyWKTStart);
1446 : // Optimize single-part single-ring multipolygon WKT->WKB translation
1447 138 : if (bCanAlterByteAfter && nLength > strlen("MULTIPOLYGON") &&
1448 28 : EQUALN(pszPtrStart, "MULTIPOLYGON", strlen("MULTIPOLYGON")))
1449 : {
1450 28 : int nCountOpenPar = 0;
1451 28 : size_t nCountComma = 0;
1452 28 : bool bHasZ = false;
1453 28 : bool bHasM = false;
1454 :
1455 28 : char *pszEnd = static_cast<char *>(pabyWKTStart) + nLength;
1456 28 : const char chBackup = *pszEnd;
1457 28 : *pszEnd = 0;
1458 :
1459 : // Checks that the multipolygon consists of a single part with
1460 : // only an exterior ring.
1461 852 : for (const char *pszPtr = pszPtrStart + strlen("MULTIPOLYGON"); *pszPtr;
1462 : ++pszPtr)
1463 : {
1464 832 : const char ch = *pszPtr;
1465 832 : if (ch == 'Z')
1466 8 : bHasZ = true;
1467 824 : else if (ch == 'M')
1468 8 : bHasM = true;
1469 832 : if (ch == '(')
1470 : {
1471 84 : nCountOpenPar++;
1472 84 : if (nCountOpenPar == 4)
1473 0 : break;
1474 : }
1475 748 : else if (ch == ')')
1476 : {
1477 72 : nCountOpenPar--;
1478 72 : if (nCountOpenPar < 0)
1479 0 : break;
1480 : }
1481 676 : else if (ch == ',')
1482 : {
1483 92 : if (nCountOpenPar < 3)
1484 : {
1485 : // multipart / multi-ring
1486 8 : break;
1487 : }
1488 84 : nCountComma++;
1489 : }
1490 : }
1491 28 : const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0);
1492 48 : if (nCountOpenPar == 0 && nCountComma > 0 &&
1493 20 : nCountComma < std::numeric_limits<uint32_t>::max())
1494 : {
1495 20 : const uint32_t nVerticesCount =
1496 20 : static_cast<uint32_t>(nCountComma + 1);
1497 20 : const size_t nWKBSize =
1498 : sizeof(GByte) + // Endianness
1499 : sizeof(uint32_t) + // multipolygon WKB geometry type
1500 : sizeof(uint32_t) + // number of parts
1501 : sizeof(GByte) + // Endianness
1502 : sizeof(uint32_t) + // polygon WKB geometry type
1503 : sizeof(uint32_t) + // number of rings
1504 : sizeof(uint32_t) + // number of vertices
1505 20 : nDim * sizeof(double) * nVerticesCount;
1506 : GByte *const pabyCurStart = static_cast<GByte *>(
1507 20 : m_oAppendBuffer.GetPtrForNewBytes(nWKBSize));
1508 20 : if (!pabyCurStart)
1509 : {
1510 0 : return static_cast<size_t>(-1);
1511 : }
1512 20 : GByte *pabyCur = pabyCurStart;
1513 : // Multipolygon byte order
1514 : {
1515 20 : *pabyCur = wkbNDR;
1516 20 : pabyCur++;
1517 : }
1518 : // Multipolygon geometry type
1519 : {
1520 20 : uint32_t nWKBGeomType =
1521 20 : wkbMultiPolygon + (bHasZ ? 1000 : 0) + (bHasM ? 2000 : 0);
1522 20 : CPL_LSBPTR32(&nWKBGeomType);
1523 20 : memcpy(pabyCur, &nWKBGeomType, sizeof(uint32_t));
1524 20 : pabyCur += sizeof(uint32_t);
1525 : }
1526 : // Number of parts
1527 : {
1528 20 : uint32_t nOne = 1;
1529 20 : CPL_LSBPTR32(&nOne);
1530 20 : memcpy(pabyCur, &nOne, sizeof(uint32_t));
1531 20 : pabyCur += sizeof(uint32_t);
1532 : }
1533 : // Polygon byte order
1534 : {
1535 20 : *pabyCur = wkbNDR;
1536 20 : pabyCur++;
1537 : }
1538 : // Polygon geometry type
1539 : {
1540 20 : uint32_t nWKBGeomType =
1541 20 : wkbPolygon + (bHasZ ? 1000 : 0) + (bHasM ? 2000 : 0);
1542 20 : CPL_LSBPTR32(&nWKBGeomType);
1543 20 : memcpy(pabyCur, &nWKBGeomType, sizeof(uint32_t));
1544 20 : pabyCur += sizeof(uint32_t);
1545 : }
1546 : // Number of rings
1547 : {
1548 20 : uint32_t nOne = 1;
1549 20 : CPL_LSBPTR32(&nOne);
1550 20 : memcpy(pabyCur, &nOne, sizeof(uint32_t));
1551 20 : pabyCur += sizeof(uint32_t);
1552 : }
1553 : // Number of vertices
1554 : {
1555 20 : uint32_t nVerticesCountToWrite = nVerticesCount;
1556 20 : CPL_LSBPTR32(&nVerticesCountToWrite);
1557 20 : memcpy(pabyCur, &nVerticesCountToWrite, sizeof(uint32_t));
1558 20 : pabyCur += sizeof(uint32_t);
1559 : }
1560 20 : uint32_t nDoubleCount = 0;
1561 20 : const uint32_t nExpectedDoubleCount = nVerticesCount * nDim;
1562 616 : for (const char *pszPtr = pszPtrStart + strlen("MULTIPOLYGON");
1563 616 : *pszPtr;
1564 : /* nothing */)
1565 : {
1566 596 : const char ch = *pszPtr;
1567 596 : if (ch == '-' || ch == '.' || (ch >= '0' && ch <= '9'))
1568 : {
1569 224 : nDoubleCount++;
1570 224 : if (nDoubleCount > nExpectedDoubleCount)
1571 : {
1572 0 : break;
1573 : }
1574 : #ifdef USE_FAST_FLOAT
1575 : double dfVal;
1576 224 : auto answer = fast_float::from_chars(pszPtr, pszEnd, dfVal);
1577 224 : if (answer.ec != std::errc())
1578 : {
1579 0 : nDoubleCount = 0;
1580 0 : break;
1581 : }
1582 224 : pszPtr = answer.ptr;
1583 : #else
1584 : char *endptr = nullptr;
1585 : const double dfVal =
1586 : m_bCanUseStrtod ? strtod(pszPtr, &endptr)
1587 : : CPLStrtodDelim(pszPtr, &endptr, '.');
1588 : pszPtr = endptr;
1589 : #endif
1590 224 : CPL_LSBPTR64(&dfVal);
1591 224 : memcpy(pabyCur, &dfVal, sizeof(double));
1592 224 : pabyCur += sizeof(double);
1593 : }
1594 : else
1595 : {
1596 372 : ++pszPtr;
1597 : }
1598 : }
1599 20 : if (nDoubleCount == nExpectedDoubleCount)
1600 : {
1601 20 : CPLAssert(static_cast<size_t>(pabyCur - pabyCurStart) ==
1602 : nWKBSize);
1603 : // cppcheck-suppress selfAssignment
1604 20 : *pszEnd = chBackup;
1605 20 : return nWKBSize;
1606 : }
1607 : else
1608 : {
1609 0 : CPLError(CE_Failure, CPLE_AppDefined,
1610 : "Invalid WKT geometry: %s", pszPtrStart);
1611 : // cppcheck-suppress selfAssignment
1612 0 : *pszEnd = chBackup;
1613 0 : return static_cast<size_t>(-1);
1614 : }
1615 : }
1616 : // cppcheck-suppress selfAssignment
1617 8 : *pszEnd = chBackup;
1618 : }
1619 :
1620 : // General case going through a OGRGeometry
1621 118 : OGRGeometry *poGeometry = nullptr;
1622 118 : if (bCanAlterByteAfter)
1623 : {
1624 : // Slight optimization for all geometries but the final one, to
1625 : // avoid creating a new string each time.
1626 : // We set the ending byte to '\0' and restore it back after parsing
1627 : // the WKT
1628 100 : char *pszEnd = static_cast<char *>(pabyWKTStart) + nLength;
1629 100 : const char chBackup = *pszEnd;
1630 100 : *pszEnd = 0;
1631 100 : OGRGeometryFactory::createFromWkt(pszPtrStart, nullptr, &poGeometry);
1632 : // cppcheck-suppress selfAssignment
1633 100 : *pszEnd = chBackup;
1634 : }
1635 : else
1636 : {
1637 36 : std::string osTmp;
1638 18 : osTmp.assign(pszPtrStart, nLength);
1639 18 : OGRGeometryFactory::createFromWkt(osTmp.c_str(), nullptr, &poGeometry);
1640 : }
1641 118 : if (!poGeometry)
1642 : {
1643 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid WKT geometry");
1644 0 : return static_cast<size_t>(-1);
1645 : }
1646 118 : const size_t nWKBSize = poGeometry->WkbSize();
1647 : GByte *pabyWKB =
1648 118 : static_cast<GByte *>(m_oAppendBuffer.GetPtrForNewBytes(nWKBSize));
1649 118 : if (!pabyWKB)
1650 : {
1651 0 : return static_cast<size_t>(-1);
1652 : }
1653 118 : poGeometry->exportToWkb(wkbNDR, pabyWKB, wkbVariantIso);
1654 118 : delete poGeometry;
1655 118 : return nWKBSize;
1656 : }
|
