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