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 208 : static inline bool OGRWKBNeedSwap(GByte b)
37 : {
38 : if constexpr (CPL_IS_LSB)
39 208 : return b != 1;
40 : else
41 : return b != 0;
42 : }
43 :
44 : /************************************************************************/
45 : /* OGRWKBReadUInt32() */
46 : /************************************************************************/
47 :
48 231 : static inline uint32_t OGRWKBReadUInt32(const GByte *pabyWkb, bool bNeedSwap)
49 : {
50 : uint32_t nVal;
51 231 : memcpy(&nVal, pabyWkb, sizeof(nVal));
52 231 : if (bNeedSwap)
53 0 : CPL_SWAP32PTR(&nVal);
54 231 : return nVal;
55 : }
56 :
57 : /************************************************************************/
58 : /* OGRWKBReadFloat64() */
59 : /************************************************************************/
60 :
61 100 : static inline double OGRWKBReadFloat64(const GByte *pabyWkb, bool bNeedSwap)
62 : {
63 : double dfVal;
64 100 : memcpy(&dfVal, pabyWkb, sizeof(dfVal));
65 100 : if (bNeedSwap)
66 0 : CPL_SWAP64PTR(&dfVal);
67 100 : return dfVal;
68 : }
69 :
70 : /************************************************************************/
71 : /* OGRWKBRingGetArea() */
72 : /************************************************************************/
73 :
74 12 : static bool OGRWKBRingGetArea(const GByte *&pabyWkb, size_t &nWKBSize, int nDim,
75 : bool bNeedSwap, double &dfArea)
76 : {
77 12 : const uint32_t nPoints = OGRWKBReadUInt32(pabyWkb, bNeedSwap);
78 12 : if (nPoints >= 4 &&
79 12 : (nWKBSize - sizeof(uint32_t)) / (nDim * sizeof(double)) >= nPoints)
80 : {
81 12 : nWKBSize -= sizeof(uint32_t) + nDim * sizeof(double);
82 12 : pabyWkb += sizeof(uint32_t);
83 : // Computation according to Green's Theorem
84 : // Cf OGRSimpleCurve::get_LinearArea()
85 12 : double x_m1 = OGRWKBReadFloat64(pabyWkb, bNeedSwap);
86 12 : double y_m1 = OGRWKBReadFloat64(pabyWkb + sizeof(double), bNeedSwap);
87 12 : double y_m2 = y_m1;
88 12 : dfArea = 0;
89 12 : pabyWkb += nDim * sizeof(double);
90 50 : for (uint32_t i = 1; i < nPoints; ++i)
91 : {
92 38 : const double x = OGRWKBReadFloat64(pabyWkb, bNeedSwap);
93 : const double y =
94 38 : OGRWKBReadFloat64(pabyWkb + sizeof(double), bNeedSwap);
95 38 : pabyWkb += nDim * sizeof(double);
96 38 : dfArea += x_m1 * (y - y_m2);
97 38 : y_m2 = y_m1;
98 38 : x_m1 = x;
99 38 : y_m1 = y;
100 : }
101 12 : dfArea += x_m1 * (y_m1 - y_m2);
102 12 : dfArea = 0.5 * std::fabs(dfArea);
103 12 : return true;
104 : }
105 0 : return false;
106 : }
107 :
108 : /************************************************************************/
109 : /* OGRWKBGetGeomType() */
110 : /************************************************************************/
111 :
112 203 : bool OGRWKBGetGeomType(const GByte *pabyWkb, size_t nWKBSize, bool &bNeedSwap,
113 : uint32_t &nType)
114 : {
115 203 : if (nWKBSize >= 5)
116 : {
117 203 : bNeedSwap = OGRWKBNeedSwap(pabyWkb[0]);
118 203 : nType = OGRWKBReadUInt32(pabyWkb + 1, bNeedSwap);
119 203 : return true;
120 : }
121 0 : return false;
122 : }
123 :
124 : /************************************************************************/
125 : /* OGRWKBPolygonGetArea() */
126 : /************************************************************************/
127 :
128 11 : bool OGRWKBPolygonGetArea(const GByte *&pabyWkb, size_t &nWKBSize,
129 : double &dfArea)
130 : {
131 : bool bNeedSwap;
132 : uint32_t nType;
133 11 : if (nWKBSize < 9 || !OGRWKBGetGeomType(pabyWkb, nWKBSize, bNeedSwap, nType))
134 0 : return false;
135 :
136 11 : int nDims = 2;
137 11 : if (nType == wkbPolygon)
138 : {
139 : // do nothing
140 : }
141 6 : else if (nType == wkbPolygon + 1000 || // wkbPolygonZ
142 4 : nType == wkbPolygon25D || nType == wkbPolygonM)
143 : {
144 4 : nDims = 3;
145 : }
146 2 : else if (nType == wkbPolygonZM)
147 : {
148 2 : nDims = 4;
149 : }
150 : else
151 : {
152 0 : return false;
153 : }
154 :
155 11 : const uint32_t nRings = OGRWKBReadUInt32(pabyWkb + 5, bNeedSwap);
156 11 : if ((nWKBSize - 9) / sizeof(uint32_t) >= nRings)
157 : {
158 11 : pabyWkb += 9;
159 11 : nWKBSize -= 9;
160 11 : dfArea = 0;
161 11 : if (nRings > 0)
162 : {
163 11 : if (!OGRWKBRingGetArea(pabyWkb, nWKBSize, nDims, bNeedSwap, dfArea))
164 0 : return false;
165 12 : for (uint32_t i = 1; i < nRings; ++i)
166 : {
167 : double dfRingArea;
168 1 : if (!OGRWKBRingGetArea(pabyWkb, nWKBSize, nDims, bNeedSwap,
169 : dfRingArea))
170 0 : return false;
171 1 : dfArea -= dfRingArea;
172 : }
173 : }
174 11 : return true;
175 : }
176 0 : return false;
177 : }
178 :
179 : /************************************************************************/
180 : /* OGRWKBMultiPolygonGetArea() */
181 : /************************************************************************/
182 :
183 5 : bool OGRWKBMultiPolygonGetArea(const GByte *&pabyWkb, size_t &nWKBSize,
184 : double &dfArea)
185 : {
186 5 : if (nWKBSize < 9)
187 0 : return false;
188 :
189 5 : const bool bNeedSwap = OGRWKBNeedSwap(pabyWkb[0]);
190 5 : const uint32_t nPolys = OGRWKBReadUInt32(pabyWkb + 5, bNeedSwap);
191 5 : if ((nWKBSize - 9) / 9 >= nPolys)
192 : {
193 5 : pabyWkb += 9;
194 5 : nWKBSize -= 9;
195 5 : dfArea = 0;
196 11 : for (uint32_t i = 0; i < nPolys; ++i)
197 : {
198 : double dfPolyArea;
199 6 : if (!OGRWKBPolygonGetArea(pabyWkb, nWKBSize, dfPolyArea))
200 0 : return false;
201 6 : dfArea += dfPolyArea;
202 : }
203 5 : return true;
204 : }
205 0 : return false;
206 : }
207 :
208 : /************************************************************************/
209 : /* WKBFromEWKB() */
210 : /************************************************************************/
211 :
212 1445 : const GByte *WKBFromEWKB(GByte *pabyEWKB, size_t nEWKBSize, size_t &nWKBSizeOut,
213 : int *pnSRIDOut)
214 : {
215 1445 : if (nEWKBSize < 5U)
216 : {
217 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid EWKB content : %u bytes",
218 : static_cast<unsigned>(nEWKBSize));
219 0 : return nullptr;
220 : }
221 :
222 1445 : const GByte *pabyWKB = pabyEWKB;
223 :
224 : /* -------------------------------------------------------------------- */
225 : /* PostGIS EWKB format includes an SRID, but this won't be */
226 : /* understood by OGR, so if the SRID flag is set, we remove the */
227 : /* SRID (bytes at offset 5 to 8). */
228 : /* -------------------------------------------------------------------- */
229 1445 : if (nEWKBSize > 9 &&
230 1438 : ((pabyEWKB[0] == 0 /* big endian */ && (pabyEWKB[1] & 0x20)) ||
231 1438 : (pabyEWKB[0] != 0 /* little endian */ && (pabyEWKB[4] & 0x20))))
232 : {
233 526 : if (pnSRIDOut)
234 : {
235 0 : memcpy(pnSRIDOut, pabyEWKB + 5, 4);
236 0 : const OGRwkbByteOrder eByteOrder =
237 0 : (pabyEWKB[0] == 0 ? wkbXDR : wkbNDR);
238 0 : if (OGR_SWAP(eByteOrder))
239 0 : *pnSRIDOut = CPL_SWAP32(*pnSRIDOut);
240 : }
241 :
242 : // Drop the SRID flag
243 526 : if (pabyEWKB[0] == 0)
244 0 : pabyEWKB[1] &= (~0x20);
245 : else
246 526 : pabyEWKB[4] &= (~0x20);
247 :
248 : // Move 5 first bytes of EWKB 4 bytes later to create regular WKB
249 526 : memmove(pabyEWKB + 4, pabyEWKB, 5);
250 526 : memset(pabyEWKB, 0, 4);
251 : // and make pabyWKB point to that
252 526 : pabyWKB += 4;
253 526 : nWKBSizeOut = nEWKBSize - 4;
254 : }
255 : else
256 : {
257 919 : if (pnSRIDOut)
258 : {
259 0 : *pnSRIDOut = INT_MIN;
260 : }
261 919 : nWKBSizeOut = nEWKBSize;
262 : }
263 :
264 1445 : return pabyWKB;
265 : }
266 :
267 : /************************************************************************/
268 : /* OGRWKBReadUInt32AtOffset() */
269 : /************************************************************************/
270 :
271 8505 : static uint32_t OGRWKBReadUInt32AtOffset(const uint8_t *data,
272 : OGRwkbByteOrder eByteOrder,
273 : size_t &iOffset)
274 : {
275 : uint32_t v;
276 8505 : memcpy(&v, data + iOffset, sizeof(v));
277 8505 : iOffset += sizeof(v);
278 8505 : if (OGR_SWAP(eByteOrder))
279 : {
280 1666 : CPL_SWAP32PTR(&v);
281 : }
282 8505 : return v;
283 : }
284 :
285 : /************************************************************************/
286 : /* ReadWKBPointSequence() */
287 : /************************************************************************/
288 :
289 : template <bool INCLUDE_Z, typename EnvelopeType>
290 852 : static bool ReadWKBPointSequence(const uint8_t *data, size_t size,
291 : OGRwkbByteOrder eByteOrder, int nDim,
292 : bool bHasZ, size_t &iOffset,
293 : EnvelopeType &sEnvelope)
294 : {
295 : const uint32_t nPoints =
296 852 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
297 852 : if (nPoints > (size - iOffset) / (nDim * sizeof(double)))
298 0 : return false;
299 852 : double dfX = 0;
300 852 : double dfY = 0;
301 852 : [[maybe_unused]] double dfZ = 0;
302 14540 : for (uint32_t j = 0; j < nPoints; j++)
303 : {
304 13688 : memcpy(&dfX, data + iOffset, sizeof(double));
305 13688 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
306 : if constexpr (INCLUDE_Z)
307 : {
308 52 : if (bHasZ)
309 34 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
310 : sizeof(double));
311 : }
312 13688 : iOffset += nDim * sizeof(double);
313 13688 : if (OGR_SWAP(eByteOrder))
314 : {
315 560 : CPL_SWAP64PTR(&dfX);
316 560 : CPL_SWAP64PTR(&dfY);
317 : if constexpr (INCLUDE_Z)
318 : {
319 0 : CPL_SWAP64PTR(&dfZ);
320 : }
321 : }
322 13688 : sEnvelope.MinX = std::min(sEnvelope.MinX, dfX);
323 13688 : sEnvelope.MinY = std::min(sEnvelope.MinY, dfY);
324 13688 : sEnvelope.MaxX = std::max(sEnvelope.MaxX, dfX);
325 13688 : sEnvelope.MaxY = std::max(sEnvelope.MaxY, dfY);
326 : if constexpr (INCLUDE_Z)
327 : {
328 52 : if (bHasZ)
329 : {
330 34 : sEnvelope.MinZ = std::min(sEnvelope.MinZ, dfZ);
331 34 : sEnvelope.MaxZ = std::max(sEnvelope.MaxZ, dfZ);
332 : }
333 : }
334 : }
335 852 : return true;
336 : }
337 :
338 : /************************************************************************/
339 : /* ReadWKBRingSequence() */
340 : /************************************************************************/
341 :
342 : template <bool INCLUDE_Z, typename EnvelopeType>
343 600 : static bool ReadWKBRingSequence(const uint8_t *data, size_t size,
344 : OGRwkbByteOrder eByteOrder, int nDim,
345 : bool bHasZ, size_t &iOffset,
346 : EnvelopeType &sEnvelope)
347 : {
348 600 : const uint32_t nRings = OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
349 600 : if (nRings > (size - iOffset) / sizeof(uint32_t))
350 0 : return false;
351 1258 : for (uint32_t i = 0; i < nRings; i++)
352 : {
353 658 : if (iOffset + sizeof(uint32_t) > size)
354 0 : return false;
355 658 : if (!ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
356 : bHasZ, iOffset, sEnvelope))
357 0 : return false;
358 : }
359 600 : return true;
360 : }
361 :
362 : /************************************************************************/
363 : /* OGRWKBGetBoundingBox() */
364 : /************************************************************************/
365 :
366 : constexpr uint32_t WKB_PREFIX_SIZE = 1 + sizeof(uint32_t);
367 : constexpr uint32_t MIN_WKB_SIZE = WKB_PREFIX_SIZE + sizeof(uint32_t);
368 :
369 : template <bool INCLUDE_Z, typename EnvelopeType>
370 990738 : static bool OGRWKBGetBoundingBox(const uint8_t *data, size_t size,
371 : size_t &iOffset, EnvelopeType &sEnvelope,
372 : int nRec)
373 : {
374 990738 : if (size - iOffset < MIN_WKB_SIZE)
375 0 : return false;
376 990738 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffset]);
377 990738 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
378 0 : return false;
379 990738 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
380 :
381 990738 : OGRwkbGeometryType eGeometryType = wkbUnknown;
382 990738 : OGRReadWKBGeometryType(data + iOffset, wkbVariantIso, &eGeometryType);
383 990738 : iOffset += 5;
384 990738 : const auto eFlatType = wkbFlatten(eGeometryType);
385 990738 : const bool bHasZ = CPL_TO_BOOL(OGR_GT_HasZ(eGeometryType));
386 990738 : const int nDim = 2 + (bHasZ ? 1 : 0) + (OGR_GT_HasM(eGeometryType) ? 1 : 0);
387 :
388 990738 : if (eFlatType == wkbPoint)
389 : {
390 990166 : if (size - iOffset < nDim * sizeof(double))
391 0 : return false;
392 990166 : double dfX = 0;
393 990166 : double dfY = 0;
394 990166 : [[maybe_unused]] double dfZ = 0;
395 990166 : memcpy(&dfX, data + iOffset, sizeof(double));
396 990166 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
397 : if constexpr (INCLUDE_Z)
398 : {
399 11 : if (bHasZ)
400 4 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
401 : sizeof(double));
402 : }
403 990166 : iOffset += nDim * sizeof(double);
404 990166 : if (OGR_SWAP(eByteOrder))
405 : {
406 40 : CPL_SWAP64PTR(&dfX);
407 40 : CPL_SWAP64PTR(&dfY);
408 : if constexpr (INCLUDE_Z)
409 : {
410 0 : CPL_SWAP64PTR(&dfZ);
411 : }
412 : }
413 990166 : if (std::isnan(dfX))
414 : {
415 : // Point empty
416 6 : sEnvelope = EnvelopeType();
417 : }
418 : else
419 : {
420 990160 : sEnvelope.MinX = dfX;
421 990160 : sEnvelope.MinY = dfY;
422 990160 : sEnvelope.MaxX = dfX;
423 990160 : sEnvelope.MaxY = dfY;
424 : if constexpr (INCLUDE_Z)
425 : {
426 10 : if (bHasZ)
427 : {
428 4 : sEnvelope.MinZ = dfZ;
429 4 : sEnvelope.MaxZ = dfZ;
430 : }
431 : }
432 : }
433 990166 : return true;
434 : }
435 :
436 572 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
437 : {
438 106 : sEnvelope = EnvelopeType();
439 :
440 106 : return ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
441 106 : bHasZ, iOffset, sEnvelope);
442 : }
443 :
444 466 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
445 : {
446 190 : sEnvelope = EnvelopeType();
447 :
448 190 : return ReadWKBRingSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
449 190 : bHasZ, iOffset, sEnvelope);
450 : }
451 :
452 276 : if (eFlatType == wkbMultiPoint)
453 : {
454 57 : sEnvelope = EnvelopeType();
455 :
456 57 : uint32_t nParts = OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
457 57 : if (nParts >
458 57 : (size - iOffset) / (WKB_PREFIX_SIZE + nDim * sizeof(double)))
459 0 : return false;
460 57 : double dfX = 0;
461 57 : double dfY = 0;
462 57 : [[maybe_unused]] double dfZ = 0;
463 170 : for (uint32_t k = 0; k < nParts; k++)
464 : {
465 113 : iOffset += WKB_PREFIX_SIZE;
466 113 : memcpy(&dfX, data + iOffset, sizeof(double));
467 113 : memcpy(&dfY, data + iOffset + sizeof(double), sizeof(double));
468 : if constexpr (INCLUDE_Z)
469 : {
470 5 : if (bHasZ)
471 3 : memcpy(&dfZ, data + iOffset + 2 * sizeof(double),
472 : sizeof(double));
473 : }
474 113 : iOffset += nDim * sizeof(double);
475 113 : if (OGR_SWAP(eByteOrder))
476 : {
477 40 : CPL_SWAP64PTR(&dfX);
478 40 : CPL_SWAP64PTR(&dfY);
479 : if constexpr (INCLUDE_Z)
480 : {
481 0 : CPL_SWAP64PTR(&dfZ);
482 : }
483 : }
484 113 : sEnvelope.MinX = std::min(sEnvelope.MinX, dfX);
485 113 : sEnvelope.MinY = std::min(sEnvelope.MinY, dfY);
486 113 : sEnvelope.MaxX = std::max(sEnvelope.MaxX, dfX);
487 113 : sEnvelope.MaxY = std::max(sEnvelope.MaxY, dfY);
488 : if constexpr (INCLUDE_Z)
489 : {
490 5 : if (bHasZ)
491 : {
492 3 : sEnvelope.MinZ = std::min(sEnvelope.MinZ, dfZ);
493 3 : sEnvelope.MaxZ = std::max(sEnvelope.MaxZ, dfZ);
494 : }
495 : }
496 : }
497 57 : return true;
498 : }
499 :
500 219 : if (eFlatType == wkbMultiLineString)
501 : {
502 59 : sEnvelope = EnvelopeType();
503 :
504 : const uint32_t nParts =
505 59 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
506 59 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
507 0 : return false;
508 147 : for (uint32_t k = 0; k < nParts; k++)
509 : {
510 88 : if (iOffset + MIN_WKB_SIZE > size)
511 0 : return false;
512 88 : iOffset += WKB_PREFIX_SIZE;
513 88 : if (!ReadWKBPointSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
514 : bHasZ, iOffset, sEnvelope))
515 0 : return false;
516 : }
517 59 : return true;
518 : }
519 :
520 160 : if (eFlatType == wkbMultiPolygon)
521 : {
522 82 : sEnvelope = EnvelopeType();
523 :
524 : const uint32_t nParts =
525 82 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
526 82 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
527 0 : return false;
528 492 : for (uint32_t k = 0; k < nParts; k++)
529 : {
530 410 : if (iOffset + MIN_WKB_SIZE > size)
531 0 : return false;
532 410 : CPLAssert(data[iOffset] == eByteOrder);
533 410 : iOffset += WKB_PREFIX_SIZE;
534 410 : if (!ReadWKBRingSequence<INCLUDE_Z>(data, size, eByteOrder, nDim,
535 : bHasZ, iOffset, sEnvelope))
536 0 : return false;
537 : }
538 82 : return true;
539 : }
540 :
541 78 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve ||
542 4 : eFlatType == wkbCurvePolygon || eFlatType == wkbMultiCurve ||
543 2 : eFlatType == wkbMultiSurface || eFlatType == wkbPolyhedralSurface ||
544 : eFlatType == wkbTIN)
545 : {
546 78 : if (nRec == 128)
547 0 : return false;
548 78 : sEnvelope = EnvelopeType();
549 :
550 : const uint32_t nParts =
551 78 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffset);
552 78 : if (nParts > (size - iOffset) / MIN_WKB_SIZE)
553 0 : return false;
554 78 : EnvelopeType sEnvelopeSubGeom;
555 176 : for (uint32_t k = 0; k < nParts; k++)
556 : {
557 98 : if (!OGRWKBGetBoundingBox<INCLUDE_Z>(data, size, iOffset,
558 : sEnvelopeSubGeom, nRec + 1))
559 0 : return false;
560 98 : sEnvelope.Merge(sEnvelopeSubGeom);
561 : }
562 78 : return true;
563 : }
564 :
565 0 : return false;
566 : }
567 :
568 : /************************************************************************/
569 : /* OGRWKBGetBoundingBox() */
570 : /************************************************************************/
571 :
572 990605 : bool OGRWKBGetBoundingBox(const GByte *pabyWkb, size_t nWKBSize,
573 : OGREnvelope &sEnvelope)
574 : {
575 990605 : size_t iOffset = 0;
576 990605 : return OGRWKBGetBoundingBox<false>(pabyWkb, nWKBSize, iOffset, sEnvelope,
577 1981210 : 0);
578 : }
579 :
580 : /************************************************************************/
581 : /* OGRWKBGetBoundingBox() */
582 : /************************************************************************/
583 :
584 35 : bool OGRWKBGetBoundingBox(const GByte *pabyWkb, size_t nWKBSize,
585 : OGREnvelope3D &sEnvelope)
586 : {
587 35 : size_t iOffset = 0;
588 70 : return OGRWKBGetBoundingBox<true>(pabyWkb, nWKBSize, iOffset, sEnvelope, 0);
589 : }
590 :
591 : /************************************************************************/
592 : /* OGRWKBIntersectsPointSequencePessimistic() */
593 : /************************************************************************/
594 :
595 199 : static bool OGRWKBIntersectsPointSequencePessimistic(
596 : const uint8_t *data, const size_t size, const OGRwkbByteOrder eByteOrder,
597 : const int nDim, size_t &iOffsetInOut, const OGREnvelope &sEnvelope,
598 : bool &bErrorOut)
599 : {
600 : const uint32_t nPoints =
601 199 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
602 199 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
603 : {
604 12 : bErrorOut = true;
605 12 : return false;
606 : }
607 :
608 187 : double dfX = 0;
609 187 : double dfY = 0;
610 1318 : for (uint32_t j = 0; j < nPoints; j++)
611 : {
612 1256 : memcpy(&dfX, data + iOffsetInOut, sizeof(double));
613 1256 : memcpy(&dfY, data + iOffsetInOut + sizeof(double), sizeof(double));
614 1256 : iOffsetInOut += nDim * sizeof(double);
615 1256 : if (OGR_SWAP(eByteOrder))
616 : {
617 3 : CPL_SWAP64PTR(&dfX);
618 3 : CPL_SWAP64PTR(&dfY);
619 : }
620 1256 : if (dfX >= sEnvelope.MinX && dfY >= sEnvelope.MinY &&
621 775 : dfX <= sEnvelope.MaxX && dfY <= sEnvelope.MaxY)
622 : {
623 125 : return true;
624 : }
625 : }
626 :
627 62 : return false;
628 : }
629 :
630 : /************************************************************************/
631 : /* OGRWKBIntersectsRingSequencePessimistic() */
632 : /************************************************************************/
633 :
634 149 : static bool OGRWKBIntersectsRingSequencePessimistic(
635 : const uint8_t *data, const size_t size, const OGRwkbByteOrder eByteOrder,
636 : const int nDim, size_t &iOffsetInOut, const OGREnvelope &sEnvelope,
637 : bool &bErrorOut)
638 : {
639 : const uint32_t nRings =
640 149 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
641 149 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
642 : {
643 16 : bErrorOut = true;
644 16 : return false;
645 : }
646 133 : if (nRings == 0)
647 2 : return false;
648 131 : if (iOffsetInOut + sizeof(uint32_t) > size)
649 : {
650 0 : bErrorOut = true;
651 0 : return false;
652 : }
653 131 : if (OGRWKBIntersectsPointSequencePessimistic(
654 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut))
655 : {
656 98 : return true;
657 : }
658 33 : if (bErrorOut)
659 0 : return false;
660 :
661 : // skip inner rings
662 39 : for (uint32_t i = 1; i < nRings; ++i)
663 : {
664 6 : if (iOffsetInOut + sizeof(uint32_t) > size)
665 : {
666 0 : bErrorOut = true;
667 0 : return false;
668 : }
669 : const uint32_t nPoints =
670 6 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
671 6 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
672 : {
673 0 : bErrorOut = true;
674 0 : return false;
675 : }
676 6 : iOffsetInOut += sizeof(double) * nPoints * nDim;
677 : }
678 33 : return false;
679 : }
680 :
681 : /************************************************************************/
682 : /* OGRWKBIntersectsPessimistic() */
683 : /************************************************************************/
684 :
685 23538 : static bool OGRWKBIntersectsPessimistic(const GByte *data, const size_t size,
686 : size_t &iOffsetInOut,
687 : const OGREnvelope &sEnvelope,
688 : const int nRec, bool &bErrorOut)
689 : {
690 23538 : if (size - iOffsetInOut < MIN_WKB_SIZE)
691 : {
692 0 : bErrorOut = true;
693 0 : return false;
694 : }
695 23538 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
696 23538 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
697 : {
698 0 : bErrorOut = true;
699 0 : return false;
700 : }
701 23538 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
702 :
703 23538 : OGRwkbGeometryType eGeometryType = wkbUnknown;
704 23538 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
705 23538 : iOffsetInOut += 5;
706 23538 : const auto eFlatType = wkbFlatten(eGeometryType);
707 23538 : const int nDim = 2 + (OGR_GT_HasZ(eGeometryType) ? 1 : 0) +
708 23538 : (OGR_GT_HasM(eGeometryType) ? 1 : 0);
709 :
710 23538 : if (eFlatType == wkbPoint)
711 : {
712 23171 : if (size - iOffsetInOut < nDim * sizeof(double))
713 8 : return false;
714 23163 : double dfX = 0;
715 23163 : double dfY = 0;
716 23163 : memcpy(&dfX, data + iOffsetInOut, sizeof(double));
717 23163 : memcpy(&dfY, data + iOffsetInOut + sizeof(double), sizeof(double));
718 23163 : iOffsetInOut += nDim * sizeof(double);
719 23163 : if (OGR_SWAP(eByteOrder))
720 : {
721 0 : CPL_SWAP64PTR(&dfX);
722 0 : CPL_SWAP64PTR(&dfY);
723 : }
724 23163 : if (std::isnan(dfX))
725 : {
726 1 : return false;
727 : }
728 : else
729 : {
730 23155 : return dfX >= sEnvelope.MinX && dfX <= sEnvelope.MaxX &&
731 46317 : dfY >= sEnvelope.MinY && dfY <= sEnvelope.MaxY;
732 : }
733 : }
734 :
735 367 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
736 : {
737 68 : return OGRWKBIntersectsPointSequencePessimistic(
738 68 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut);
739 : }
740 :
741 299 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
742 : {
743 149 : return OGRWKBIntersectsRingSequencePessimistic(
744 149 : data, size, eByteOrder, nDim, iOffsetInOut, sEnvelope, bErrorOut);
745 : }
746 :
747 150 : if (eFlatType == wkbMultiPoint || eFlatType == wkbMultiLineString ||
748 81 : eFlatType == wkbMultiPolygon || eFlatType == wkbGeometryCollection ||
749 65 : eFlatType == wkbCompoundCurve || eFlatType == wkbCurvePolygon ||
750 39 : eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface ||
751 13 : eFlatType == wkbPolyhedralSurface || eFlatType == wkbTIN)
752 : {
753 150 : if (nRec == 128)
754 : {
755 0 : bErrorOut = true;
756 0 : return false;
757 : }
758 : const uint32_t nParts =
759 150 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
760 150 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
761 : {
762 80 : bErrorOut = true;
763 80 : return false;
764 : }
765 115 : for (uint32_t k = 0; k < nParts; k++)
766 : {
767 79 : if (OGRWKBIntersectsPessimistic(data, size, iOffsetInOut, sEnvelope,
768 : nRec + 1, bErrorOut))
769 : {
770 34 : return true;
771 : }
772 45 : else if (bErrorOut)
773 : {
774 0 : return false;
775 : }
776 : }
777 36 : return false;
778 : }
779 :
780 0 : bErrorOut = true;
781 0 : return false;
782 : }
783 :
784 : /************************************************************************/
785 : /* OGRWKBIntersectsPessimistic() */
786 : /************************************************************************/
787 :
788 : /* Returns whether the geometry (pabyWkb, nWKBSize) intersects, for sure,
789 : * the passed envelope.
790 : * When it returns true, the geometry intersects the envelope.
791 : * When it returns false, the geometry may or may not intersect the envelope.
792 : */
793 23459 : bool OGRWKBIntersectsPessimistic(const GByte *pabyWkb, size_t nWKBSize,
794 : const OGREnvelope &sEnvelope)
795 : {
796 23459 : size_t iOffsetInOut = 0;
797 23459 : bool bErrorOut = false;
798 23459 : bool bRet = OGRWKBIntersectsPessimistic(pabyWkb, nWKBSize, iOffsetInOut,
799 : sEnvelope, 0, bErrorOut);
800 23459 : if (!bRet && !bErrorOut)
801 : {
802 : // The following assert only holds if there is no trailing data
803 : // after the WKB
804 : // CPLAssert(iOffsetInOut == nWKBSize);
805 : }
806 23459 : return bRet;
807 : }
808 :
809 : /************************************************************************/
810 : /* epsilonEqual() */
811 : /************************************************************************/
812 :
813 92 : static inline bool epsilonEqual(double a, double b, double eps)
814 : {
815 92 : return ::fabs(a - b) < eps;
816 : }
817 :
818 : /************************************************************************/
819 : /* OGRWKBIsClockwiseRing() */
820 : /************************************************************************/
821 :
822 382 : static inline double GetX(const GByte *data, uint32_t i, int nDim,
823 : bool bNeedSwap)
824 : {
825 : double dfX;
826 382 : memcpy(&dfX, data + static_cast<size_t>(i) * nDim * sizeof(double),
827 : sizeof(double));
828 382 : if (bNeedSwap)
829 0 : CPL_SWAP64PTR(&dfX);
830 382 : return dfX;
831 : }
832 :
833 653 : static inline double GetY(const GByte *data, uint32_t i, int nDim,
834 : bool bNeedSwap)
835 : {
836 : double dfY;
837 653 : memcpy(&dfY, data + (static_cast<size_t>(i) * nDim + 1) * sizeof(double),
838 : sizeof(double));
839 653 : if (bNeedSwap)
840 0 : CPL_SWAP64PTR(&dfY);
841 653 : return dfY;
842 : }
843 :
844 37 : static bool OGRWKBIsClockwiseRing(const GByte *data, const uint32_t nPoints,
845 : const int nDim, const bool bNeedSwap)
846 : {
847 37 : constexpr double EPSILON = 1.0E-5;
848 :
849 : // WARNING: keep in sync OGRLineString::isClockwise(),
850 : // OGRCurve::isClockwise() and OGRWKBIsClockwiseRing()
851 :
852 37 : bool bUseFallback = false;
853 :
854 : // Find the lowest rightmost vertex.
855 37 : uint32_t v = 0; // Used after for.
856 37 : double vX = GetX(data, v, nDim, bNeedSwap);
857 37 : double vY = GetY(data, v, nDim, bNeedSwap);
858 535 : for (uint32_t i = 1; i < nPoints - 1; i++)
859 : {
860 : // => v < end.
861 498 : const double y = GetY(data, i, nDim, bNeedSwap);
862 498 : if (y < vY)
863 : {
864 163 : v = i;
865 163 : vX = GetX(data, i, nDim, bNeedSwap);
866 163 : vY = y;
867 163 : bUseFallback = false;
868 : }
869 335 : else if (y == vY)
870 : {
871 18 : const double x = GetX(data, i, nDim, bNeedSwap);
872 18 : if (x > vX)
873 : {
874 8 : v = i;
875 8 : vX = x;
876 : // vY = y;
877 8 : bUseFallback = false;
878 : }
879 10 : else if (x == vX)
880 : {
881 : // Two vertex with same coordinates are the lowest rightmost
882 : // vertex. Cannot use that point as the pivot (#5342).
883 2 : bUseFallback = true;
884 : }
885 : }
886 : }
887 :
888 : // Previous.
889 37 : uint32_t next = (v == 0) ? nPoints - 2 : v - 1;
890 49 : if (epsilonEqual(GetX(data, next, nDim, bNeedSwap), vX, EPSILON) &&
891 12 : epsilonEqual(GetY(data, next, nDim, bNeedSwap), vY, EPSILON))
892 : {
893 : // Don't try to be too clever by retrying with a next point.
894 : // This can lead to false results as in the case of #3356.
895 0 : bUseFallback = true;
896 : }
897 :
898 37 : const double dx0 = GetX(data, next, nDim, bNeedSwap) - vX;
899 37 : const double dy0 = GetY(data, next, nDim, bNeedSwap) - vY;
900 :
901 : // Following.
902 37 : next = v + 1;
903 37 : if (next >= nPoints - 1)
904 : {
905 2 : next = 0;
906 : }
907 :
908 43 : if (epsilonEqual(GetX(data, next, nDim, bNeedSwap), vX, EPSILON) &&
909 6 : epsilonEqual(GetY(data, next, nDim, bNeedSwap), vY, EPSILON))
910 : {
911 : // Don't try to be too clever by retrying with a next point.
912 : // This can lead to false results as in the case of #3356.
913 2 : bUseFallback = true;
914 : }
915 :
916 37 : const double dx1 = GetX(data, next, nDim, bNeedSwap) - vX;
917 37 : const double dy1 = GetY(data, next, nDim, bNeedSwap) - vY;
918 :
919 37 : const double crossproduct = dx1 * dy0 - dx0 * dy1;
920 :
921 37 : if (!bUseFallback)
922 : {
923 35 : if (crossproduct > 0) // CCW
924 23 : return false;
925 12 : else if (crossproduct < 0) // CW
926 12 : return true;
927 : }
928 :
929 : // This is a degenerate case: the extent of the polygon is less than EPSILON
930 : // or 2 nearly identical points were found.
931 : // Try with Green Formula as a fallback, but this is not a guarantee
932 : // as we'll probably be affected by numerical instabilities.
933 :
934 2 : double dfSum = GetX(data, 0, nDim, bNeedSwap) *
935 2 : (GetY(data, 1, nDim, bNeedSwap) -
936 2 : GetY(data, nPoints - 1, nDim, bNeedSwap));
937 :
938 12 : for (uint32_t i = 1; i < nPoints - 1; i++)
939 : {
940 10 : dfSum += GetX(data, i, nDim, bNeedSwap) *
941 20 : (GetY(data, i + 1, nDim, bNeedSwap) -
942 10 : GetY(data, i - 1, nDim, bNeedSwap));
943 : }
944 :
945 2 : dfSum += GetX(data, nPoints - 1, nDim, bNeedSwap) *
946 2 : (GetY(data, 0, nDim, bNeedSwap) -
947 2 : GetX(data, nPoints - 2, nDim, bNeedSwap));
948 :
949 2 : return dfSum < 0;
950 : }
951 :
952 : /************************************************************************/
953 : /* OGRWKBFixupCounterClockWiseExternalRing() */
954 : /************************************************************************/
955 :
956 63 : static bool OGRWKBFixupCounterClockWiseExternalRingInternal(
957 : GByte *data, size_t size, size_t &iOffsetInOut, const int nRec)
958 : {
959 63 : if (size - iOffsetInOut < MIN_WKB_SIZE)
960 : {
961 0 : return false;
962 : }
963 63 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
964 63 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
965 : {
966 0 : return false;
967 : }
968 63 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
969 :
970 63 : OGRwkbGeometryType eGeometryType = wkbUnknown;
971 63 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
972 63 : iOffsetInOut += 5;
973 63 : const auto eFlatType = wkbFlatten(eGeometryType);
974 63 : const int nDim = 2 + (OGR_GT_HasZ(eGeometryType) ? 1 : 0) +
975 63 : (OGR_GT_HasM(eGeometryType) ? 1 : 0);
976 :
977 63 : if (eFlatType == wkbPolygon)
978 : {
979 : const uint32_t nRings =
980 34 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
981 34 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
982 : {
983 0 : return false;
984 : }
985 71 : for (uint32_t iRing = 0; iRing < nRings; ++iRing)
986 : {
987 37 : if (iOffsetInOut + sizeof(uint32_t) > size)
988 0 : return false;
989 : const uint32_t nPoints =
990 37 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
991 37 : const size_t sizeOfPoint = nDim * sizeof(double);
992 37 : if (nPoints > (size - iOffsetInOut) / sizeOfPoint)
993 : {
994 0 : return false;
995 : }
996 :
997 37 : if (nPoints >= 4)
998 : {
999 37 : const bool bIsClockwiseRing = OGRWKBIsClockwiseRing(
1000 37 : data + iOffsetInOut, nPoints, nDim, OGR_SWAP(eByteOrder));
1001 37 : if ((bIsClockwiseRing && iRing == 0) ||
1002 25 : (!bIsClockwiseRing && iRing > 0))
1003 : {
1004 : GByte abyTmp[4 * sizeof(double)];
1005 43 : for (uint32_t i = 0; i < nPoints / 2; ++i)
1006 : {
1007 29 : GByte *pBegin = data + iOffsetInOut + i * sizeOfPoint;
1008 29 : GByte *pEnd = data + iOffsetInOut +
1009 29 : (nPoints - 1 - i) * sizeOfPoint;
1010 29 : memcpy(abyTmp, pBegin, sizeOfPoint);
1011 29 : memcpy(pBegin, pEnd, sizeOfPoint);
1012 29 : memcpy(pEnd, abyTmp, sizeOfPoint);
1013 : }
1014 : }
1015 : }
1016 :
1017 37 : iOffsetInOut += nPoints * sizeOfPoint;
1018 : }
1019 : }
1020 :
1021 63 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbMultiPolygon ||
1022 : eFlatType == wkbMultiSurface)
1023 : {
1024 6 : if (nRec == 128)
1025 : {
1026 0 : return false;
1027 : }
1028 : const uint32_t nParts =
1029 6 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1030 6 : if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE)
1031 : {
1032 0 : return false;
1033 : }
1034 11 : for (uint32_t k = 0; k < nParts; k++)
1035 : {
1036 5 : if (!OGRWKBFixupCounterClockWiseExternalRingInternal(
1037 : data, size, iOffsetInOut, nRec))
1038 : {
1039 0 : return false;
1040 : }
1041 : }
1042 : }
1043 :
1044 63 : return true;
1045 : }
1046 :
1047 : /** Modifies the geometry such that exterior rings of polygons are
1048 : * counter-clockwise oriented and inner rings clockwise oriented.
1049 : */
1050 58 : void OGRWKBFixupCounterClockWiseExternalRing(GByte *pabyWkb, size_t nWKBSize)
1051 : {
1052 58 : size_t iOffsetInOut = 0;
1053 58 : OGRWKBFixupCounterClockWiseExternalRingInternal(
1054 : pabyWkb, nWKBSize, iOffsetInOut, /* nRec = */ 0);
1055 58 : }
1056 :
1057 : /************************************************************************/
1058 : /* OGRWKBPointUpdater() */
1059 : /************************************************************************/
1060 :
1061 : OGRWKBPointUpdater::OGRWKBPointUpdater() = default;
1062 :
1063 : OGRWKBPointUpdater::~OGRWKBPointUpdater() = default;
1064 :
1065 : /************************************************************************/
1066 : /* OGRWKBIntersectsPointSequencePessimistic() */
1067 : /************************************************************************/
1068 :
1069 3029 : 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 3029 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1077 3029 : if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double)))
1078 : {
1079 1497 : return false;
1080 : }
1081 1532 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1082 6067 : for (uint32_t j = 0; j < nPoints; j++)
1083 : {
1084 4537 : void *pdfX = data + iOffsetInOut;
1085 4537 : void *pdfY = data + iOffsetInOut + sizeof(double);
1086 4537 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1087 4537 : void *pdfM =
1088 4537 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1089 : : nullptr;
1090 4537 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1091 2 : return false;
1092 :
1093 4535 : iOffsetInOut += nDim * sizeof(double);
1094 : }
1095 :
1096 1530 : return true;
1097 : }
1098 :
1099 : /************************************************************************/
1100 : /* OGRWKBVisitRingSequence() */
1101 : /************************************************************************/
1102 :
1103 1776 : 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 1776 : OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut);
1111 1776 : if (nRings > (size - iOffsetInOut) / sizeof(uint32_t))
1112 : {
1113 150 : return false;
1114 : }
1115 :
1116 2534 : for (uint32_t i = 0; i < nRings; ++i)
1117 : {
1118 1727 : if (iOffsetInOut + sizeof(uint32_t) > size)
1119 : {
1120 4 : return false;
1121 : }
1122 1723 : if (!OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder, nDim,
1123 : bHasZ, bHasM, iOffsetInOut))
1124 : {
1125 815 : return false;
1126 : }
1127 : }
1128 807 : return true;
1129 : }
1130 :
1131 : /************************************************************************/
1132 : /* OGRWKBUpdatePoints() */
1133 : /************************************************************************/
1134 :
1135 26303 : static bool OGRWKBUpdatePoints(uint8_t *data, const size_t size,
1136 : OGRWKBPointUpdater &oUpdater,
1137 : size_t &iOffsetInOut, const int nRec)
1138 : {
1139 26303 : if (size - iOffsetInOut < MIN_WKB_SIZE)
1140 : {
1141 514 : return false;
1142 : }
1143 25789 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]);
1144 25789 : if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR))
1145 : {
1146 68 : return false;
1147 : }
1148 25721 : const OGRwkbByteOrder eByteOrder = static_cast<OGRwkbByteOrder>(nByteOrder);
1149 :
1150 25721 : OGRwkbGeometryType eGeometryType = wkbUnknown;
1151 25721 : OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType);
1152 25721 : iOffsetInOut += 5;
1153 25721 : const auto eFlatType = wkbFlatten(eGeometryType);
1154 :
1155 25721 : if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve ||
1156 25250 : eFlatType == wkbCurvePolygon || eFlatType == wkbMultiPoint ||
1157 24939 : eFlatType == wkbMultiLineString || eFlatType == wkbMultiPolygon ||
1158 24751 : eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface ||
1159 24491 : 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 24329 : const bool bHasZ = CPL_TO_BOOL(OGR_GT_HasZ(eGeometryType));
1180 24329 : const bool bHasM = CPL_TO_BOOL(OGR_GT_HasM(eGeometryType));
1181 24329 : const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0);
1182 :
1183 24329 : if (eFlatType == wkbPoint)
1184 : {
1185 21044 : if (size - iOffsetInOut < nDim * sizeof(double))
1186 376 : return false;
1187 20668 : void *pdfX = data + iOffsetInOut;
1188 20668 : void *pdfY = data + iOffsetInOut + sizeof(double);
1189 20668 : void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr;
1190 20668 : void *pdfM =
1191 20668 : bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double)
1192 : : nullptr;
1193 20668 : const bool bNeedSwap = OGR_SWAP(eByteOrder);
1194 20668 : if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM))
1195 1 : return false;
1196 20667 : iOffsetInOut += nDim * sizeof(double);
1197 20667 : return true;
1198 : }
1199 :
1200 3285 : if (eFlatType == wkbLineString || eFlatType == wkbCircularString)
1201 : {
1202 1306 : return OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder,
1203 1306 : nDim, bHasZ, bHasM, iOffsetInOut);
1204 : }
1205 :
1206 1979 : if (eFlatType == wkbPolygon || eFlatType == wkbTriangle)
1207 : {
1208 1776 : return OGRWKBVisitRingSequence(data, size, oUpdater, eByteOrder, nDim,
1209 1776 : 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 24945 : bool OGRWKBUpdatePoints(GByte *pabyWkb, size_t nWKBSize,
1219 : OGRWKBPointUpdater &oUpdater)
1220 : {
1221 24945 : size_t iOffsetInOut = 0;
1222 24945 : return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater, iOffsetInOut,
1223 49890 : /* 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 14918 : 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 14918 : explicit OGRWKBPointUpdaterReproj(OGRCoordinateTransformation *poCTIn,
1265 : OGREnvelope3D &sEnvelopeIn)
1266 14918 : : m_poCT(poCTIn), m_sEnvelope(sEnvelopeIn)
1267 : {
1268 14918 : }
1269 :
1270 14943 : bool update(bool bNeedSwap, void *x, void *y, void *z,
1271 : void * /* m */) override
1272 : {
1273 : double dfX, dfY, dfZ;
1274 14943 : memcpy(&dfX, x, sizeof(double));
1275 14943 : memcpy(&dfY, y, sizeof(double));
1276 14943 : if (bNeedSwap)
1277 : {
1278 1116 : CPL_SWAP64PTR(&dfX);
1279 1116 : CPL_SWAP64PTR(&dfY);
1280 : }
1281 14943 : if (!(std::isnan(dfX) && std::isnan(dfY)))
1282 : {
1283 14743 : 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 13577 : dfZ = 0;
1293 14743 : int nErrorCode = 0;
1294 14743 : m_poCT->TransformWithErrorCodes(1, &dfX, &dfY, &dfZ, nullptr,
1295 14743 : &nErrorCode);
1296 14743 : if (nErrorCode)
1297 3 : return false;
1298 14740 : m_sEnvelope.Merge(dfX, dfY, dfZ);
1299 14740 : if (bNeedSwap)
1300 : {
1301 1016 : CPL_SWAP64PTR(&dfX);
1302 1016 : CPL_SWAP64PTR(&dfY);
1303 1016 : CPL_SWAP64PTR(&dfZ);
1304 : }
1305 14740 : memcpy(x, &dfX, sizeof(double));
1306 14740 : memcpy(y, &dfY, sizeof(double));
1307 14740 : if (z)
1308 1166 : memcpy(z, &dfZ, sizeof(double));
1309 : }
1310 14940 : return true;
1311 : }
1312 :
1313 : private:
1314 : OGRWKBPointUpdaterReproj(const OGRWKBPointUpdaterReproj &) = delete;
1315 : OGRWKBPointUpdaterReproj &
1316 : operator=(const OGRWKBPointUpdaterReproj &) = delete;
1317 : };
1318 :
1319 14918 : sEnvelope = OGREnvelope3D();
1320 29836 : OGRWKBPointUpdaterReproj oUpdater(poCT, sEnvelope);
1321 29836 : 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 236 : std::unique_ptr<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 : poGeometry = OGRGeometryFactory::createFromWkt(pszPtrStart).first;
1632 : // cppcheck-suppress selfAssignment
1633 100 : *pszEnd = chBackup;
1634 : }
1635 : else
1636 : {
1637 18 : std::string osTmp;
1638 18 : osTmp.assign(pszPtrStart, nLength);
1639 18 : poGeometry = OGRGeometryFactory::createFromWkt(osTmp.c_str()).first;
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 : return nWKBSize;
1655 : }
|
