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 8473 : static uint32_t OGRWKBReadUInt32AtOffset(const uint8_t *data,
274 : OGRwkbByteOrder eByteOrder,
275 : size_t &iOffset)
276 : {
277 : uint32_t v;
278 8473 : memcpy(&v, data + iOffset, sizeof(v));
279 8473 : iOffset += sizeof(v);
280 8473 : if (OGR_SWAP(eByteOrder))
281 : {
282 1666 : CPL_SWAP32PTR(&v);
283 : }
284 8473 : 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 958554 : static bool OGRWKBGetBoundingBox(const uint8_t *data, size_t size,
373 : size_t &iOffset, EnvelopeType &sEnvelope,
374 : int nRec)
375 : {
376 958554 : if (size - iOffset < MIN_WKB_SIZE)
377 0 : return false;
378 958554 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffset]);
379 958554 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
380 0 : return false;
381 958554 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
382 :
383 958554 : OGRwkbGeometryType eGeometryType = wkbUnknown;
384 958554 : OGRReadWKBGeometryType(data + iOffset, wkbVariantIso, &eGeometryType);
385 958554 : iOffset += 5;
386 958554 : const auto eFlatType = wkbFlatten(eGeometryType);
387 958554 : const bool bHasZ = CPL_TO_BOOL(OGR_GT_HasZ(eGeometryType));
388 958554 : const int nDim = 2 + (bHasZ ? 1 : 0) + (OGR_GT_HasM(eGeometryType) ? 1 : 0);
389 :
390 958554 : if (eFlatType == wkbPoint)
391 : {
392 957990 : if (size - iOffset < nDim * sizeof(double))
393 0 : return false;
394 957990 : double dfX = 0;
395 957990 : double dfY = 0;
396 957990 : [[maybe_unused]] double dfZ = 0;
397 957990 : memcpy(&dfX, data + iOffset, sizeof(double));
398 957990 : 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 957990 : iOffset += nDim * sizeof(double);
406 957990 : 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 957990 : if (std::isnan(dfX))
416 : {
417 : // Point empty
418 6 : sEnvelope = EnvelopeType();
419 : }
420 : else
421 : {
422 957984 : sEnvelope.MinX = dfX;
423 957984 : sEnvelope.MinY = dfY;
424 957984 : sEnvelope.MaxX = dfX;
425 957984 : 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 957990 : 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 958421 : bool OGRWKBGetBoundingBox(const GByte *pabyWkb, size_t nWKBSize,
575 : OGREnvelope &sEnvelope)
576 : {
577 958421 : size_t iOffset = 0;
578 958421 : return OGRWKBGetBoundingBox<false>(pabyWkb, nWKBSize, iOffset, sEnvelope,
579 1916840 : 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 199 : 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 199 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
604 199 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
605 : {
606 12 : bErrorOut = true;
607 12 : return false;
608 : }
609 :
610 187 : double dfX = 0;
611 187 : double dfY = 0;
612 1318 : for (uint32_t j = 0; j < nPoints; j++)
613 : {
614 1256 : memcpy(&dfX, data + iOffsetInOut, sizeof(double));
615 1256 : memcpy(&dfY, data + iOffsetInOut + sizeof(double), sizeof(double));
616 1256 : iOffsetInOut += nDim * sizeof(double);
617 1256 : if (OGR_SWAP(eByteOrder))
618 : {
619 3 : CPL_SWAP64PTR(&dfX);
620 3 : CPL_SWAP64PTR(&dfY);
621 : }
622 1256 : if (dfX >= sEnvelope.MinX && dfY >= sEnvelope.MinY &&
623 775 : dfX <= sEnvelope.MaxX && dfY <= sEnvelope.MaxY)
624 : {
625 125 : return true;
626 : }
627 : }
628 :
629 62 : return false;
630 : }
631 :
632 : /************************************************************************/
633 : /* OGRWKBIntersectsRingSequencePessimistic() */
634 : /************************************************************************/
635 :
636 149 : 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 149 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
643 149 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
644 : {
645 16 : bErrorOut = true;
646 16 : return false;
647 : }
648 133 : if (nRings == 0)
649 2 : return false;
650 131 : if (iOffsetInOut + sizeof(uint32_t) > size)
651 : {
652 0 : bErrorOut = true;
653 0 : return false;
654 : }
655 131 : if (OGRWKBIntersectsPointSequencePessimistic(
656 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut))
657 : {
658 98 : 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 23437 : 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 23437 : if (size - iOffsetInOut < MIN_WKB_SIZE)
693 : {
694 0 : bErrorOut = true;
695 0 : return false;
696 : }
697 23437 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
698 23437 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
699 : {
700 0 : bErrorOut = true;
701 0 : return false;
702 : }
703 23437 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
704 :
705 23437 : OGRwkbGeometryType eGeometryType = wkbUnknown;
706 23437 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
707 23437 : iOffsetInOut += 5;
708 23437 : const auto eFlatType = wkbFlatten(eGeometryType);
709 23437 : const int nDim = 2 + (OGR_GT_HasZ(eGeometryType) ? 1 : 0) +
710 23437 : (OGR_GT_HasM(eGeometryType) ? 1 : 0);
711 :
712 23437 : 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 367 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
738 : {
739 68 : return OGRWKBIntersectsPointSequencePessimistic(
740 68 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut);
741 : }
742 :
743 299 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
744 : {
745 149 : return OGRWKBIntersectsRingSequencePessimistic(
746 149 : 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 23358 : bool OGRWKBIntersectsPessimistic(const GByte *pabyWkb, size_t nWKBSize,
796 : const OGREnvelope &sEnvelope)
797 : {
798 23358 : size_t iOffsetInOut = 0;
799 23358 : bool bErrorOut = false;
800 23358 : bool bRet = OGRWKBIntersectsPessimistic(pabyWkb, nWKBSize, iOffsetInOut,
801 : sEnvelope, 0, bErrorOut);
802 23358 : 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 23358 : 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 : OGRWKBPointUpdater::~OGRWKBPointUpdater() = default;
1066 :
1067 : /************************************************************************/
1068 : /* OGRWKBIntersectsPointSequencePessimistic() */
1069 : /************************************************************************/
1070 :
1071 3029 : static bool OGRWKBUpdatePointsSequence(uint8_t *data, const size_t size,
1072 : OGRWKBPointUpdater &oUpdater,
1073 : const OGRwkbByteOrder eByteOrder,
1074 : const int nDim, const bool bHasZ,
1075 : const bool bHasM, size_t &iOffsetInOut)
1076 : {
1077 : const uint32_t nPoints =
1078 3029 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1079 3029 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
1080 : {
1081 1497 : return false;
1082 : }
1083 1532 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1084 6067 : for (uint32_t j = 0; j < nPoints; j++)
1085 : {
1086 4537 : void *pdfX = data + iOffsetInOut;
1087 4537 : void *pdfY = data + iOffsetInOut + sizeof(double);
1088 4537 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1089 4537 : void *pdfM =
1090 4537 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1091 : : nullptr;
1092 4537 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1093 2 : return false;
1094 :
1095 4535 : iOffsetInOut += nDim * sizeof(double);
1096 : }
1097 :
1098 1530 : return true;
1099 : }
1100 :
1101 : /************************************************************************/
1102 : /* OGRWKBVisitRingSequence() */
1103 : /************************************************************************/
1104 :
1105 1776 : static bool OGRWKBVisitRingSequence(uint8_t *data, const size_t size,
1106 : OGRWKBPointUpdater &oUpdater,
1107 : const OGRwkbByteOrder eByteOrder,
1108 : const int nDim, const bool bHasZ,
1109 : const bool bHasM, size_t &iOffsetInOut)
1110 : {
1111 : const uint32_t nRings =
1112 1776 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1113 1776 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
1114 : {
1115 150 : return false;
1116 : }
1117 :
1118 2534 : for (uint32_t i = 0; i < nRings; ++i)
1119 : {
1120 1727 : if (iOffsetInOut + sizeof(uint32_t) > size)
1121 : {
1122 4 : return false;
1123 : }
1124 1723 : if (!OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder, nDim,
1125 : bHasZ, bHasM, iOffsetInOut))
1126 : {
1127 815 : return false;
1128 : }
1129 : }
1130 807 : return true;
1131 : }
1132 :
1133 : /************************************************************************/
1134 : /* OGRWKBUpdatePoints() */
1135 : /************************************************************************/
1136 :
1137 25838 : static bool OGRWKBUpdatePoints(uint8_t *data, const size_t size,
1138 : OGRWKBPointUpdater &oUpdater,
1139 : size_t &iOffsetInOut, const int nRec)
1140 : {
1141 25838 : if (size - iOffsetInOut < MIN_WKB_SIZE)
1142 : {
1143 514 : return false;
1144 : }
1145 25324 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
1146 25324 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
1147 : {
1148 68 : return false;
1149 : }
1150 25256 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
1151 :
1152 25256 : OGRwkbGeometryType eGeometryType = wkbUnknown;
1153 25256 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
1154 25133 : iOffsetInOut += 5;
1155 25133 : const auto eFlatType = wkbFlatten(eGeometryType);
1156 :
1157 24988 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve ||
1158 24560 : eFlatType == wkbCurvePolygon || eFlatType == wkbMultiPoint ||
1159 24172 : eFlatType == wkbMultiLineString || eFlatType == wkbMultiPolygon ||
1160 23944 : eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface ||
1161 23679 : eFlatType == wkbPolyhedralSurface || eFlatType == wkbTIN)
1162 : {
1163 1400 : if (nRec == 128)
1164 1 : return false;
1165 :
1166 : const uint32_t nParts =
1167 1399 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1168 1391 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
1169 : {
1170 176 : return false;
1171 : }
1172 1856 : for (uint32_t k = 0; k < nParts; k++)
1173 : {
1174 1358 : if (!OGRWKBUpdatePoints(data, size, oUpdater, iOffsetInOut,
1175 : nRec + 1))
1176 717 : return false;
1177 : }
1178 498 : return true;
1179 : }
1180 :
1181 23588 : const bool bHasZ = OGR_GT_HasZ(eGeometryType);
1182 23820 : const bool bHasM = OGR_GT_HasM(eGeometryType);
1183 23685 : const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0);
1184 :
1185 23685 : if (eFlatType == wkbPoint)
1186 : {
1187 20122 : if (size - iOffsetInOut < nDim * sizeof(double))
1188 376 : return false;
1189 19746 : void *pdfX = data + iOffsetInOut;
1190 19746 : void *pdfY = data + iOffsetInOut + sizeof(double);
1191 19746 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1192 19746 : void *pdfM =
1193 19746 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1194 : : nullptr;
1195 19746 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1196 19746 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1197 1 : return false;
1198 20471 : iOffsetInOut += nDim * sizeof(double);
1199 20471 : return true;
1200 : }
1201 :
1202 3563 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
1203 : {
1204 1584 : return OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder,
1205 1306 : nDim, bHasZ, bHasM, iOffsetInOut);
1206 : }
1207 :
1208 1979 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
1209 : {
1210 1776 : return OGRWKBVisitRingSequence(data, size, oUpdater, eByteOrder, nDim,
1211 1776 : bHasZ, bHasM, iOffsetInOut);
1212 : }
1213 :
1214 203 : CPLDebug("OGR", "Unknown WKB geometry type");
1215 203 : return false;
1216 : }
1217 :
1218 : /** Visit all points of a WKB geometry to update them.
1219 : */
1220 24552 : bool OGRWKBUpdatePoints(GByte *pabyWkb, size_t nWKBSize,
1221 : OGRWKBPointUpdater &oUpdater)
1222 : {
1223 24552 : size_t iOffsetInOut = 0;
1224 24552 : return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater, iOffsetInOut,
1225 49502 : /* nRec = */ 0);
1226 : }
1227 :
1228 : /************************************************************************/
1229 : /* OGRWKBTransformCache::clear() */
1230 : /************************************************************************/
1231 :
1232 : #ifdef OGR_WKB_TRANSFORM_ALL_AT_ONCE
1233 : void OGRWKBTransformCache::clear()
1234 : {
1235 : abNeedSwap.clear();
1236 : abIsEmpty.clear();
1237 : apdfX.clear();
1238 : apdfY.clear();
1239 : apdfZ.clear();
1240 : apdfM.clear();
1241 : adfX.clear();
1242 : adfY.clear();
1243 : adfZ.clear();
1244 : adfM.clear();
1245 : anErrorCodes.clear();
1246 : }
1247 : #endif
1248 :
1249 : /************************************************************************/
1250 : /* OGRWKBTransform() */
1251 : /************************************************************************/
1252 :
1253 : /** Visit all points of a WKB geometry to transform them.
1254 : */
1255 14507 : bool OGRWKBTransform(GByte *pabyWkb, size_t nWKBSize,
1256 : OGRCoordinateTransformation *poCT,
1257 : [[maybe_unused]] OGRWKBTransformCache &oCache,
1258 : OGREnvelope3D &sEnvelope)
1259 : {
1260 : #ifndef OGR_WKB_TRANSFORM_ALL_AT_ONCE
1261 : struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater
1262 : {
1263 : OGRCoordinateTransformation *m_poCT;
1264 : OGREnvelope3D &m_sEnvelope;
1265 :
1266 14525 : explicit OGRWKBPointUpdaterReproj(OGRCoordinateTransformation *poCTIn,
1267 : OGREnvelope3D &sEnvelopeIn)
1268 14525 : : m_poCT(poCTIn), m_sEnvelope(sEnvelopeIn)
1269 : {
1270 14576 : }
1271 :
1272 14301 : bool update(bool bNeedSwap, void *x, void *y, void *z,
1273 : void * /* m */) override
1274 : {
1275 : double dfX, dfY, dfZ;
1276 14301 : memcpy(&dfX, x, sizeof(double));
1277 14301 : memcpy(&dfY, y, sizeof(double));
1278 14301 : if (bNeedSwap)
1279 : {
1280 1116 : CPL_SWAP64PTR(&dfX);
1281 1116 : CPL_SWAP64PTR(&dfY);
1282 : }
1283 14301 : if (!(std::isnan(dfX) && std::isnan(dfY)))
1284 : {
1285 13895 : if (z)
1286 : {
1287 1166 : memcpy(&dfZ, z, sizeof(double));
1288 1166 : if (bNeedSwap)
1289 : {
1290 646 : CPL_SWAP64PTR(&dfZ);
1291 : }
1292 : }
1293 : else
1294 12729 : dfZ = 0;
1295 13895 : int nErrorCode = 0;
1296 13895 : m_poCT->TransformWithErrorCodes(1, &dfX, &dfY, &dfZ, nullptr,
1297 13895 : &nErrorCode);
1298 14534 : if (nErrorCode)
1299 3 : return false;
1300 14531 : m_sEnvelope.Merge(dfX, dfY, dfZ);
1301 14592 : if (bNeedSwap)
1302 : {
1303 1016 : CPL_SWAP64PTR(&dfX);
1304 1016 : CPL_SWAP64PTR(&dfY);
1305 1016 : CPL_SWAP64PTR(&dfZ);
1306 : }
1307 14526 : memcpy(x, &dfX, sizeof(double));
1308 14526 : memcpy(y, &dfY, sizeof(double));
1309 14526 : if (z)
1310 1166 : memcpy(z, &dfZ, sizeof(double));
1311 : }
1312 14536 : return true;
1313 : }
1314 :
1315 : private:
1316 : OGRWKBPointUpdaterReproj(const OGRWKBPointUpdaterReproj &) = delete;
1317 : OGRWKBPointUpdaterReproj &
1318 : operator=(const OGRWKBPointUpdaterReproj &) = delete;
1319 : };
1320 :
1321 14507 : sEnvelope = OGREnvelope3D();
1322 29237 : OGRWKBPointUpdaterReproj oUpdater(poCT, sEnvelope);
1323 29104 : return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater);
1324 :
1325 : #else
1326 : struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater
1327 : {
1328 : OGRWKBTransformCache &m_oCache;
1329 :
1330 : explicit OGRWKBPointUpdaterReproj(OGRWKBTransformCache &oCacheIn)
1331 : : m_oCache(oCacheIn)
1332 : {
1333 : }
1334 :
1335 : bool update(bool bNeedSwap, void *x, void *y, void *z,
1336 : void * /* m */) override
1337 : {
1338 : m_oCache.abNeedSwap.push_back(bNeedSwap);
1339 : m_oCache.apdfX.push_back(x);
1340 : m_oCache.apdfY.push_back(y);
1341 : m_oCache.apdfZ.push_back(z);
1342 : return true;
1343 : }
1344 : };
1345 :
1346 : oCache.clear();
1347 : OGRWKBPointUpdaterReproj oUpdater(oCache);
1348 : if (!OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater))
1349 : return false;
1350 :
1351 : oCache.adfX.resize(oCache.apdfX.size());
1352 : oCache.adfY.resize(oCache.apdfX.size());
1353 : oCache.adfZ.resize(oCache.apdfX.size());
1354 :
1355 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1356 : {
1357 : memcpy(&oCache.adfX[i], oCache.apdfX[i], sizeof(double));
1358 : memcpy(&oCache.adfY[i], oCache.apdfY[i], sizeof(double));
1359 : if (oCache.apdfZ[i])
1360 : memcpy(&oCache.adfZ[i], oCache.apdfZ[i], sizeof(double));
1361 : if (oCache.abNeedSwap[i])
1362 : {
1363 : CPL_SWAP64PTR(&oCache.adfX[i]);
1364 : CPL_SWAP64PTR(&oCache.adfY[i]);
1365 : CPL_SWAP64PTR(&oCache.adfZ[i]);
1366 : }
1367 : oCache.abIsEmpty.push_back(std::isnan(oCache.adfX[i]) &&
1368 : std::isnan(oCache.adfY[i]));
1369 : }
1370 :
1371 : oCache.anErrorCodes.resize(oCache.apdfX.size());
1372 : poCT->TransformWithErrorCodes(static_cast<int>(oCache.apdfX.size()),
1373 : oCache.adfX.data(), oCache.adfY.data(),
1374 : oCache.adfZ.data(), nullptr,
1375 : oCache.anErrorCodes.data());
1376 :
1377 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1378 : {
1379 : if (!oCache.abIsEmpty[i] && oCache.anErrorCodes[i])
1380 : return false;
1381 : }
1382 :
1383 : sEnvelope = OGREnvelope3D();
1384 : for (size_t i = 0; i < oCache.apdfX.size(); ++i)
1385 : {
1386 : if (oCache.abIsEmpty[i])
1387 : {
1388 : oCache.adfX[i] = std::numeric_limits<double>::quiet_NaN();
1389 : oCache.adfY[i] = std::numeric_limits<double>::quiet_NaN();
1390 : oCache.adfZ[i] = std::numeric_limits<double>::quiet_NaN();
1391 : }
1392 : else
1393 : {
1394 : sEnvelope.Merge(oCache.adfX[i], oCache.adfY[i], oCache.adfZ[i]);
1395 : }
1396 : if (oCache.abNeedSwap[i])
1397 : {
1398 : CPL_SWAP64PTR(&oCache.adfX[i]);
1399 : CPL_SWAP64PTR(&oCache.adfY[i]);
1400 : CPL_SWAP64PTR(&oCache.adfZ[i]);
1401 : }
1402 : memcpy(oCache.apdfX[i], &oCache.adfX[i], sizeof(double));
1403 : memcpy(oCache.apdfY[i], &oCache.adfY[i], sizeof(double));
1404 : if (oCache.apdfZ[i])
1405 : memcpy(oCache.apdfZ[i], &oCache.adfZ[i], sizeof(double));
1406 : }
1407 :
1408 : return true;
1409 : #endif
1410 : }
1411 :
1412 : /************************************************************************/
1413 : /* OGRAppendBuffer() */
1414 : /************************************************************************/
1415 :
1416 : OGRAppendBuffer::OGRAppendBuffer() = default;
1417 :
1418 : /************************************************************************/
1419 : /* ~OGRAppendBuffer() */
1420 : /************************************************************************/
1421 :
1422 : OGRAppendBuffer::~OGRAppendBuffer() = default;
1423 :
1424 : /************************************************************************/
1425 : /* OGRWKTToWKBTranslator() */
1426 : /************************************************************************/
1427 :
1428 18 : OGRWKTToWKBTranslator::OGRWKTToWKBTranslator(OGRAppendBuffer &oAppendBuffer)
1429 18 : : m_oAppendBuffer(oAppendBuffer)
1430 : {
1431 : #ifndef USE_FAST_FLOAT
1432 : // Test if current locale decimal separator is decimal point
1433 : char szTest[10];
1434 : snprintf(szTest, sizeof(szTest), "%f", 1.5);
1435 : m_bCanUseStrtod = strchr(szTest, '.') != nullptr;
1436 : #endif
1437 18 : CPL_IGNORE_RET_VAL(m_bCanUseStrtod);
1438 18 : }
1439 :
1440 : /************************************************************************/
1441 : /* TranslateWKT() */
1442 : /************************************************************************/
1443 :
1444 138 : size_t OGRWKTToWKBTranslator::TranslateWKT(void *pabyWKTStart, size_t nLength,
1445 : bool bCanAlterByteAfter)
1446 : {
1447 138 : const char *pszPtrStart = static_cast<const char *>(pabyWKTStart);
1448 : // Optimize single-part single-ring multipolygon WKT->WKB translation
1449 138 : if (bCanAlterByteAfter && nLength > strlen("MULTIPOLYGON") &&
1450 28 : EQUALN(pszPtrStart, "MULTIPOLYGON", strlen("MULTIPOLYGON")))
1451 : {
1452 28 : int nCountOpenPar = 0;
1453 28 : size_t nCountComma = 0;
1454 28 : bool bHasZ = false;
1455 28 : bool bHasM = false;
1456 :
1457 28 : char *pszEnd = static_cast<char *>(pabyWKTStart) + nLength;
1458 28 : const char chBackup = *pszEnd;
1459 28 : *pszEnd = 0;
1460 :
1461 : // Checks that the multipolygon consists of a single part with
1462 : // only an exterior ring.
1463 852 : for (const char *pszPtr = pszPtrStart + strlen("MULTIPOLYGON"); *pszPtr;
1464 : ++pszPtr)
1465 : {
1466 832 : const char ch = *pszPtr;
1467 832 : if (ch == 'Z')
1468 8 : bHasZ = true;
1469 824 : else if (ch == 'M')
1470 8 : bHasM = true;
1471 832 : if (ch == '(')
1472 : {
1473 84 : nCountOpenPar++;
1474 84 : if (nCountOpenPar == 4)
1475 0 : break;
1476 : }
1477 748 : else if (ch == ')')
1478 : {
1479 72 : nCountOpenPar--;
1480 72 : if (nCountOpenPar < 0)
1481 0 : break;
1482 : }
1483 676 : else if (ch == ',')
1484 : {
1485 92 : if (nCountOpenPar < 3)
1486 : {
1487 : // multipart / multi-ring
1488 8 : break;
1489 : }
1490 84 : nCountComma++;
1491 : }
1492 : }
1493 28 : const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0);
1494 48 : if (nCountOpenPar == 0 && nCountComma > 0 &&
1495 20 : nCountComma < std::numeric_limits<uint32_t>::max())
1496 : {
1497 20 : const uint32_t nVerticesCount =
1498 20 : static_cast<uint32_t>(nCountComma + 1);
1499 20 : const size_t nWKBSize =
1500 : sizeof(GByte) + // Endianness
1501 : sizeof(uint32_t) + // multipolygon WKB geometry type
1502 : sizeof(uint32_t) + // number of parts
1503 : sizeof(GByte) + // Endianness
1504 : sizeof(uint32_t) + // polygon WKB geometry type
1505 : sizeof(uint32_t) + // number of rings
1506 : sizeof(uint32_t) + // number of vertices
1507 20 : nDim * sizeof(double) * nVerticesCount;
1508 : GByte *const pabyCurStart = static_cast<GByte *>(
1509 20 : m_oAppendBuffer.GetPtrForNewBytes(nWKBSize));
1510 20 : if (!pabyCurStart)
1511 : {
1512 0 : return static_cast<size_t>(-1);
1513 : }
1514 20 : GByte *pabyCur = pabyCurStart;
1515 : // Multipolygon byte order
1516 : {
1517 20 : *pabyCur = wkbNDR;
1518 20 : pabyCur++;
1519 : }
1520 : // Multipolygon geometry type
1521 : {
1522 20 : uint32_t nWKBGeomType =
1523 20 : wkbMultiPolygon + (bHasZ ? 1000 : 0) + (bHasM ? 2000 : 0);
1524 20 : CPL_LSBPTR32(&nWKBGeomType);
1525 20 : memcpy(pabyCur, &nWKBGeomType, sizeof(uint32_t));
1526 20 : pabyCur += sizeof(uint32_t);
1527 : }
1528 : // Number of parts
1529 : {
1530 20 : uint32_t nOne = 1;
1531 20 : CPL_LSBPTR32(&nOne);
1532 20 : memcpy(pabyCur, &nOne, sizeof(uint32_t));
1533 20 : pabyCur += sizeof(uint32_t);
1534 : }
1535 : // Polygon byte order
1536 : {
1537 20 : *pabyCur = wkbNDR;
1538 20 : pabyCur++;
1539 : }
1540 : // Polygon geometry type
1541 : {
1542 20 : uint32_t nWKBGeomType =
1543 20 : wkbPolygon + (bHasZ ? 1000 : 0) + (bHasM ? 2000 : 0);
1544 20 : CPL_LSBPTR32(&nWKBGeomType);
1545 20 : memcpy(pabyCur, &nWKBGeomType, sizeof(uint32_t));
1546 20 : pabyCur += sizeof(uint32_t);
1547 : }
1548 : // Number of rings
1549 : {
1550 20 : uint32_t nOne = 1;
1551 20 : CPL_LSBPTR32(&nOne);
1552 20 : memcpy(pabyCur, &nOne, sizeof(uint32_t));
1553 20 : pabyCur += sizeof(uint32_t);
1554 : }
1555 : // Number of vertices
1556 : {
1557 20 : uint32_t nVerticesCountToWrite = nVerticesCount;
1558 20 : CPL_LSBPTR32(&nVerticesCountToWrite);
1559 20 : memcpy(pabyCur, &nVerticesCountToWrite, sizeof(uint32_t));
1560 20 : pabyCur += sizeof(uint32_t);
1561 : }
1562 20 : uint32_t nDoubleCount = 0;
1563 20 : const uint32_t nExpectedDoubleCount = nVerticesCount * nDim;
1564 616 : for (const char *pszPtr = pszPtrStart + strlen("MULTIPOLYGON");
1565 616 : *pszPtr;
1566 : /* nothing */)
1567 : {
1568 596 : const char ch = *pszPtr;
1569 596 : if (ch == '-' || ch == '.' || (ch >= '0' && ch <= '9'))
1570 : {
1571 224 : nDoubleCount++;
1572 224 : if (nDoubleCount > nExpectedDoubleCount)
1573 : {
1574 0 : break;
1575 : }
1576 : #ifdef USE_FAST_FLOAT
1577 : double dfVal;
1578 224 : auto answer = fast_float::from_chars(pszPtr, pszEnd, dfVal);
1579 224 : if (answer.ec != std::errc())
1580 : {
1581 0 : nDoubleCount = 0;
1582 0 : break;
1583 : }
1584 224 : pszPtr = answer.ptr;
1585 : #else
1586 : char *endptr = nullptr;
1587 : const double dfVal =
1588 : m_bCanUseStrtod ? strtod(pszPtr, &endptr)
1589 : : CPLStrtodDelim(pszPtr, &endptr, '.');
1590 : pszPtr = endptr;
1591 : #endif
1592 224 : CPL_LSBPTR64(&dfVal);
1593 224 : memcpy(pabyCur, &dfVal, sizeof(double));
1594 224 : pabyCur += sizeof(double);
1595 : }
1596 : else
1597 : {
1598 372 : ++pszPtr;
1599 : }
1600 : }
1601 20 : if (nDoubleCount == nExpectedDoubleCount)
1602 : {
1603 20 : CPLAssert(static_cast<size_t>(pabyCur - pabyCurStart) ==
1604 : nWKBSize);
1605 : // cppcheck-suppress selfAssignment
1606 20 : *pszEnd = chBackup;
1607 20 : return nWKBSize;
1608 : }
1609 : else
1610 : {
1611 0 : CPLError(CE_Failure, CPLE_AppDefined,
1612 : "Invalid WKT geometry: %s", pszPtrStart);
1613 : // cppcheck-suppress selfAssignment
1614 0 : *pszEnd = chBackup;
1615 0 : return static_cast<size_t>(-1);
1616 : }
1617 : }
1618 : // cppcheck-suppress selfAssignment
1619 8 : *pszEnd = chBackup;
1620 : }
1621 :
1622 : // General case going through a OGRGeometry
1623 236 : std::unique_ptr<OGRGeometry> poGeometry = nullptr;
1624 118 : if (bCanAlterByteAfter)
1625 : {
1626 : // Slight optimization for all geometries but the final one, to
1627 : // avoid creating a new string each time.
1628 : // We set the ending byte to '\0' and restore it back after parsing
1629 : // the WKT
1630 100 : char *pszEnd = static_cast<char *>(pabyWKTStart) + nLength;
1631 100 : const char chBackup = *pszEnd;
1632 100 : *pszEnd = 0;
1633 100 : poGeometry = OGRGeometryFactory::createFromWkt(pszPtrStart).first;
1634 : // cppcheck-suppress selfAssignment
1635 100 : *pszEnd = chBackup;
1636 : }
1637 : else
1638 : {
1639 18 : std::string osTmp;
1640 18 : osTmp.assign(pszPtrStart, nLength);
1641 18 : poGeometry = OGRGeometryFactory::createFromWkt(osTmp.c_str()).first;
1642 : }
1643 118 : if (!poGeometry)
1644 : {
1645 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid WKT geometry");
1646 0 : return static_cast<size_t>(-1);
1647 : }
1648 118 : const size_t nWKBSize = poGeometry->WkbSize();
1649 : GByte *pabyWKB =
1650 118 : static_cast<GByte *>(m_oAppendBuffer.GetPtrForNewBytes(nWKBSize));
1651 118 : if (!pabyWKB)
1652 : {
1653 0 : return static_cast<size_t>(-1);
1654 : }
1655 118 : poGeometry->exportToWkb(wkbNDR, pabyWKB, wkbVariantIso);
1656 118 : return nWKBSize;
1657 : }
|
