Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: Implements decoder of shapebin geometry for PGeo
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : * Paul Ramsey, pramsey at cleverelephant.ca
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
10 : * Copyright (c) 2011, Paul Ramsey <pramsey at cleverelephant.ca>
11 : * Copyright (c) 2011-2014, Even Rouault <even dot rouault at spatialys.com>
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : // PGeo == ESRI Personal GeoDatabase.
17 :
18 : #include "cpl_port.h"
19 : #include "ogrpgeogeometry.h"
20 :
21 : #include <algorithm>
22 : #include <cassert>
23 : #include <cmath>
24 : #include <cstddef>
25 : #include <cstring>
26 : #include <limits>
27 : #include <map>
28 : #include <memory>
29 : #include <set>
30 : #include <utility>
31 : #include <vector>
32 :
33 : #include "cpl_conv.h"
34 : #include "cpl_error.h"
35 : #include "cpl_string.h"
36 : #include "cpl_vsi.h"
37 : #include "ogr_api.h"
38 : #include "ogr_core.h"
39 : #include "ogr_p.h"
40 :
41 : #ifdef HAVE_WFLAG_UNREACHABLE_CODE_AGGRESSIVE
42 : #if defined(__clang__)
43 : #pragma clang diagnostic push
44 : #pragma clang diagnostic ignored "-Wunreachable-code"
45 : #endif
46 : #endif
47 :
48 : constexpr int SHPP_TRISTRIP = 0;
49 : constexpr int SHPP_TRIFAN = 1;
50 : constexpr int SHPP_OUTERRING = 2;
51 : constexpr int SHPP_INNERRING = 3;
52 : constexpr int SHPP_FIRSTRING = 4;
53 : constexpr int SHPP_RING = 5;
54 : constexpr int SHPP_TRIANGLES = 6; // Multipatch 9.0 specific.
55 :
56 : enum CurveType
57 : {
58 : CURVE_ARC_INTERIOR_POINT,
59 : CURVE_ARC_CENTER_POINT,
60 : CURVE_BEZIER,
61 : CURVE_ELLIPSE_BY_CENTER
62 : };
63 :
64 : namespace
65 : {
66 : struct CurveSegment
67 : {
68 : int nStartPointIdx;
69 : CurveType eType;
70 :
71 : union
72 : {
73 : // Arc defined by an intermediate point.
74 : struct
75 : {
76 : double dfX;
77 : double dfY;
78 : } ArcByIntermediatePoint;
79 :
80 : // Deprecated way of defining circular arc by its center and
81 : // winding order.
82 : struct
83 : {
84 : double dfX;
85 : double dfY;
86 : bool bIsCCW;
87 : } ArcByCenterPoint;
88 :
89 : struct
90 : {
91 : double dfX1;
92 : double dfY1;
93 : double dfX2;
94 : double dfY2;
95 : } Bezier;
96 :
97 : struct
98 : {
99 : double dfX;
100 : double dfY;
101 : double dfRotationDeg;
102 : double dfSemiMajor;
103 : double dfRatioSemiMinor;
104 : bool bIsMinor;
105 : bool bIsComplete;
106 : } EllipseByCenter;
107 : } u;
108 : };
109 : } /* namespace */
110 :
111 : constexpr int EXT_SHAPE_SEGMENT_ARC = 1;
112 : constexpr int EXT_SHAPE_SEGMENT_BEZIER = 4;
113 : constexpr int EXT_SHAPE_SEGMENT_ELLIPSE = 5;
114 :
115 : constexpr int EXT_SHAPE_ARC_EMPTY = 0x1;
116 : constexpr int EXT_SHAPE_ARC_CCW = 0x8;
117 : #ifdef DEBUG_VERBOSE
118 : constexpr int EXT_SHAPE_ARC_MINOR = 0x10;
119 : #endif
120 : constexpr int EXT_SHAPE_ARC_LINE = 0x20;
121 : constexpr int EXT_SHAPE_ARC_POINT = 0x40;
122 : constexpr int EXT_SHAPE_ARC_IP = 0x80;
123 :
124 : #ifdef DEBUG_VERBOSE
125 : constexpr int EXT_SHAPE_ELLIPSE_EMPTY = 0x1;
126 : constexpr int EXT_SHAPE_ELLIPSE_LINE = 0x40;
127 : constexpr int EXT_SHAPE_ELLIPSE_POINT = 0x80;
128 : constexpr int EXT_SHAPE_ELLIPSE_CIRCULAR = 0x100;
129 : constexpr int EXT_SHAPE_ELLIPSE_CCW = 0x800;
130 : #endif
131 :
132 : constexpr int EXT_SHAPE_ELLIPSE_CENTER_TO = 0x200;
133 : constexpr int EXT_SHAPE_ELLIPSE_CENTER_FROM = 0x400;
134 : constexpr int EXT_SHAPE_ELLIPSE_MINOR = 0x1000;
135 : constexpr int EXT_SHAPE_ELLIPSE_COMPLETE = 0x2000;
136 :
137 : /************************************************************************/
138 : /* OGRCreateFromMultiPatchPart() */
139 : /************************************************************************/
140 :
141 3536 : static void OGRCreateFromMultiPatchPart(OGRGeometryCollection *poGC,
142 : OGRMultiPolygon *&poMP,
143 : OGRPolygon *&poLastPoly, int nPartType,
144 : int nPartPoints, const double *padfX,
145 : const double *padfY,
146 : const double *padfZ)
147 : {
148 3536 : nPartType &= 0xf;
149 :
150 3536 : if (nPartType == SHPP_TRISTRIP)
151 : {
152 800 : if (poMP != nullptr && poLastPoly != nullptr)
153 : {
154 0 : poMP->addGeometryDirectly(poLastPoly);
155 0 : poLastPoly = nullptr;
156 : }
157 :
158 800 : OGRTriangulatedSurface *poTIN = new OGRTriangulatedSurface();
159 2404 : for (int iBaseVert = 0; iBaseVert < nPartPoints - 2; iBaseVert++)
160 : {
161 1604 : const int iSrcVert = iBaseVert;
162 :
163 3208 : OGRPoint oPoint1(padfX[iSrcVert], padfY[iSrcVert], padfZ[iSrcVert]);
164 :
165 1604 : OGRPoint oPoint2(padfX[iSrcVert + 1], padfY[iSrcVert + 1],
166 3208 : padfZ[iSrcVert + 1]);
167 :
168 1604 : OGRPoint oPoint3(padfX[iSrcVert + 2], padfY[iSrcVert + 2],
169 3208 : padfZ[iSrcVert + 2]);
170 :
171 : OGRTriangle *poTriangle =
172 1604 : new OGRTriangle(oPoint1, oPoint2, oPoint3);
173 :
174 1604 : poTIN->addGeometryDirectly(poTriangle);
175 : }
176 800 : poGC->addGeometryDirectly(poTIN);
177 : }
178 2736 : else if (nPartType == SHPP_TRIFAN)
179 : {
180 691 : if (poMP != nullptr && poLastPoly != nullptr)
181 : {
182 0 : poMP->addGeometryDirectly(poLastPoly);
183 0 : poLastPoly = nullptr;
184 : }
185 :
186 691 : OGRTriangulatedSurface *poTIN = new OGRTriangulatedSurface();
187 2068 : for (int iBaseVert = 0; iBaseVert < nPartPoints - 2; iBaseVert++)
188 : {
189 1377 : const int iSrcVert = iBaseVert;
190 :
191 2754 : OGRPoint oPoint1(padfX[0], padfY[0], padfZ[0]);
192 :
193 1377 : OGRPoint oPoint2(padfX[iSrcVert + 1], padfY[iSrcVert + 1],
194 2754 : padfZ[iSrcVert + 1]);
195 :
196 1377 : OGRPoint oPoint3(padfX[iSrcVert + 2], padfY[iSrcVert + 2],
197 2754 : padfZ[iSrcVert + 2]);
198 :
199 : OGRTriangle *poTriangle =
200 1377 : new OGRTriangle(oPoint1, oPoint2, oPoint3);
201 :
202 1377 : poTIN->addGeometryDirectly(poTriangle);
203 : }
204 691 : poGC->addGeometryDirectly(poTIN);
205 : }
206 2045 : else if (nPartType == SHPP_OUTERRING || nPartType == SHPP_INNERRING ||
207 682 : nPartType == SHPP_FIRSTRING || nPartType == SHPP_RING)
208 : {
209 1363 : if (poMP == nullptr)
210 : {
211 681 : poMP = new OGRMultiPolygon();
212 : }
213 :
214 1363 : if (poMP != nullptr && poLastPoly != nullptr &&
215 682 : (nPartType == SHPP_OUTERRING || nPartType == SHPP_FIRSTRING))
216 : {
217 0 : poMP->addGeometryDirectly(poLastPoly);
218 0 : poLastPoly = nullptr;
219 : }
220 :
221 1363 : if (poLastPoly == nullptr)
222 681 : poLastPoly = new OGRPolygon();
223 :
224 1363 : OGRLinearRing *poRing = new OGRLinearRing;
225 :
226 1363 : poRing->setPoints(nPartPoints, const_cast<double *>(padfX),
227 : const_cast<double *>(padfY),
228 : const_cast<double *>(padfZ));
229 :
230 1363 : poRing->closeRings();
231 :
232 1363 : poLastPoly->addRingDirectly(poRing);
233 : }
234 682 : else if (nPartType == SHPP_TRIANGLES)
235 : {
236 682 : if (poMP != nullptr && poLastPoly != nullptr)
237 : {
238 0 : poMP->addGeometryDirectly(poLastPoly);
239 0 : poLastPoly = nullptr;
240 : }
241 :
242 682 : OGRTriangulatedSurface *poTIN = new OGRTriangulatedSurface();
243 1364 : for (int iBaseVert = 0; iBaseVert < nPartPoints - 2; iBaseVert += 3)
244 : {
245 682 : const int iSrcVert = iBaseVert;
246 :
247 1364 : OGRPoint oPoint1(padfX[iSrcVert], padfY[iSrcVert], padfZ[iSrcVert]);
248 :
249 682 : OGRPoint oPoint2(padfX[iSrcVert + 1], padfY[iSrcVert + 1],
250 1364 : padfZ[iSrcVert + 1]);
251 :
252 682 : OGRPoint oPoint3(padfX[iSrcVert + 2], padfY[iSrcVert + 2],
253 1364 : padfZ[iSrcVert + 2]);
254 :
255 : OGRTriangle *poTriangle =
256 682 : new OGRTriangle(oPoint1, oPoint2, oPoint3);
257 :
258 682 : poTIN->addGeometryDirectly(poTriangle);
259 : }
260 682 : poGC->addGeometryDirectly(poTIN);
261 : }
262 : else
263 0 : CPLDebug("OGR", "Unrecognized parttype %d, ignored.", nPartType);
264 3536 : }
265 :
266 : static bool
267 12 : RegisterEdge(const double *padfX, const double *padfY, const double *padfZ,
268 : int nPart,
269 : std::map<std::vector<double>, std::pair<int, int>> &oMapEdges)
270 : {
271 12 : int idx = 0;
272 12 : if (padfX[0] > padfX[1])
273 : {
274 6 : idx = 1;
275 : }
276 6 : else if (padfX[0] == padfX[1])
277 : {
278 2 : if (padfY[0] > padfY[1])
279 : {
280 0 : idx = 1;
281 : }
282 2 : else if (padfY[0] == padfY[1])
283 : {
284 0 : if (padfZ[0] > padfZ[1])
285 : {
286 0 : idx = 1;
287 : }
288 : }
289 : }
290 24 : std::vector<double> oVector;
291 12 : oVector.push_back(padfX[idx]);
292 12 : oVector.push_back(padfY[idx]);
293 12 : oVector.push_back(padfZ[idx]);
294 12 : oVector.push_back(padfX[1 - idx]);
295 12 : oVector.push_back(padfY[1 - idx]);
296 12 : oVector.push_back(padfZ[1 - idx]);
297 12 : const auto oIter = oMapEdges.find(oVector);
298 12 : if (oIter == oMapEdges.end())
299 : {
300 10 : oMapEdges[oVector] = std::pair(nPart, -1);
301 : }
302 : else
303 : {
304 2 : CPLAssert(oIter->second.first >= 0);
305 2 : if (oIter->second.second < 0)
306 2 : oIter->second.second = nPart;
307 : else
308 0 : return false;
309 : }
310 12 : return true;
311 : }
312 :
313 12 : static const std::pair<int, int> &GetEdgeOwners(
314 : const double *padfX, const double *padfY, const double *padfZ,
315 : const std::map<std::vector<double>, std::pair<int, int>> &oMapEdges)
316 : {
317 12 : int idx = 0;
318 12 : if (padfX[0] > padfX[1])
319 : {
320 6 : idx = 1;
321 : }
322 6 : else if (padfX[0] == padfX[1])
323 : {
324 2 : if (padfY[0] > padfY[1])
325 : {
326 0 : idx = 1;
327 : }
328 2 : else if (padfY[0] == padfY[1])
329 : {
330 0 : if (padfZ[0] > padfZ[1])
331 : {
332 0 : idx = 1;
333 : }
334 : }
335 : }
336 12 : std::vector<double> oVector;
337 12 : oVector.push_back(padfX[idx]);
338 12 : oVector.push_back(padfY[idx]);
339 12 : oVector.push_back(padfZ[idx]);
340 12 : oVector.push_back(padfX[1 - idx]);
341 12 : oVector.push_back(padfY[1 - idx]);
342 12 : oVector.push_back(padfZ[1 - idx]);
343 12 : const auto oIter = oMapEdges.find(oVector);
344 12 : CPLAssert(oIter != oMapEdges.end());
345 24 : return oIter->second;
346 : }
347 :
348 : /************************************************************************/
349 : /* OGRCreateFromMultiPatch() */
350 : /* */
351 : /* Translate a multipatch representation to an OGR geometry */
352 : /************************************************************************/
353 :
354 812 : OGRGeometry *OGRCreateFromMultiPatch(int nParts, const GInt32 *panPartStart,
355 : const GInt32 *panPartType, int nPoints,
356 : const double *padfX, const double *padfY,
357 : const double *padfZ)
358 : {
359 : // Deal with particular case of a patch of OuterRing of 4 points
360 : // that form a TIN. And be robust to consecutive duplicated triangles !
361 1624 : std::map<std::vector<double>, std::pair<int, int>> oMapEdges;
362 812 : bool bTINCandidate = nParts >= 2;
363 1624 : std::set<int> oSetDuplicated;
364 817 : for (int iPart = 0; iPart < nParts && panPartStart != nullptr; iPart++)
365 : {
366 815 : int nPartPoints = 0;
367 :
368 : // Figure out details about this part's vertex list.
369 815 : if (iPart == nParts - 1)
370 129 : nPartPoints = nPoints - panPartStart[iPart];
371 : else
372 686 : nPartPoints = panPartStart[iPart + 1] - panPartStart[iPart];
373 815 : const int nPartStart = panPartStart[iPart];
374 :
375 5 : if (panPartType[iPart] == SHPP_OUTERRING && nPartPoints == 4 &&
376 5 : padfX[nPartStart] == padfX[nPartStart + 3] &&
377 5 : padfY[nPartStart] == padfY[nPartStart + 3] &&
378 5 : padfZ[nPartStart] == padfZ[nPartStart + 3] &&
379 5 : !std::isnan(padfX[nPartStart]) &&
380 5 : !std::isnan(padfX[nPartStart + 1]) &&
381 5 : !std::isnan(padfX[nPartStart + 2]) &&
382 5 : !std::isnan(padfY[nPartStart]) &&
383 5 : !std::isnan(padfY[nPartStart + 1]) &&
384 5 : !std::isnan(padfY[nPartStart + 2]) &&
385 5 : !std::isnan(padfZ[nPartStart]) &&
386 825 : !std::isnan(padfZ[nPartStart + 1]) &&
387 5 : !std::isnan(padfZ[nPartStart + 2]))
388 : {
389 5 : bool bDuplicate = false;
390 5 : if (iPart > 0)
391 : {
392 3 : bDuplicate = true;
393 3 : const int nPrevPartStart = panPartStart[iPart - 1];
394 6 : for (int j = 0; j < 3; j++)
395 : {
396 5 : if (padfX[nPartStart + j] == padfX[nPrevPartStart + j] &&
397 3 : padfY[nPartStart + j] == padfY[nPrevPartStart + j] &&
398 3 : padfZ[nPartStart + j] == padfZ[nPrevPartStart + j])
399 : {
400 : }
401 : else
402 : {
403 2 : bDuplicate = false;
404 2 : break;
405 : }
406 : }
407 : }
408 5 : if (bDuplicate)
409 : {
410 1 : oSetDuplicated.insert(iPart);
411 : }
412 12 : else if (RegisterEdge(padfX + nPartStart, padfY + nPartStart,
413 4 : padfZ + nPartStart, iPart, oMapEdges) &&
414 4 : RegisterEdge(padfX + nPartStart + 1,
415 4 : padfY + nPartStart + 1,
416 8 : padfZ + nPartStart + 1, iPart, oMapEdges) &&
417 4 : RegisterEdge(padfX + nPartStart + 2,
418 4 : padfY + nPartStart + 2,
419 4 : padfZ + nPartStart + 2, iPart, oMapEdges))
420 : {
421 : // ok
422 : }
423 : else
424 : {
425 0 : bTINCandidate = false;
426 0 : break;
427 : }
428 : }
429 : else
430 : {
431 810 : bTINCandidate = false;
432 810 : break;
433 : }
434 : }
435 812 : if (bTINCandidate && panPartStart != nullptr)
436 : {
437 2 : std::set<int> oVisitedParts;
438 2 : std::set<int> oToBeVisitedParts;
439 2 : oToBeVisitedParts.insert(0);
440 6 : while (!oToBeVisitedParts.empty())
441 : {
442 4 : const int iPart = *(oToBeVisitedParts.begin());
443 4 : oVisitedParts.insert(iPart);
444 4 : oToBeVisitedParts.erase(iPart);
445 :
446 4 : const int nPartStart = panPartStart[iPart];
447 16 : for (int j = 0; j < 3; j++)
448 : {
449 : const auto &oPair = GetEdgeOwners(
450 12 : padfX + nPartStart + j, padfY + nPartStart + j,
451 12 : padfZ + nPartStart + j, oMapEdges);
452 12 : const int iOtherPart =
453 12 : (oPair.first == iPart) ? oPair.second : oPair.first;
454 16 : if (iOtherPart >= 0 &&
455 16 : oVisitedParts.find(iOtherPart) == oVisitedParts.end())
456 : {
457 2 : oToBeVisitedParts.insert(iOtherPart);
458 : }
459 : }
460 : }
461 2 : if (static_cast<int>(oVisitedParts.size()) ==
462 2 : nParts - static_cast<int>(oSetDuplicated.size()))
463 : {
464 2 : OGRTriangulatedSurface *poTIN = new OGRTriangulatedSurface();
465 7 : for (int iPart = 0; iPart < nParts; iPart++)
466 : {
467 5 : if (oSetDuplicated.find(iPart) != oSetDuplicated.end())
468 1 : continue;
469 :
470 4 : const int nPartStart = panPartStart[iPart];
471 4 : OGRPoint oPoint1(padfX[nPartStart], padfY[nPartStart],
472 8 : padfZ[nPartStart]);
473 :
474 4 : OGRPoint oPoint2(padfX[nPartStart + 1], padfY[nPartStart + 1],
475 8 : padfZ[nPartStart + 1]);
476 :
477 4 : OGRPoint oPoint3(padfX[nPartStart + 2], padfY[nPartStart + 2],
478 8 : padfZ[nPartStart + 2]);
479 :
480 : OGRTriangle *poTriangle =
481 4 : new OGRTriangle(oPoint1, oPoint2, oPoint3);
482 :
483 4 : poTIN->addGeometryDirectly(poTriangle);
484 : }
485 2 : return poTIN;
486 : }
487 : }
488 :
489 810 : OGRGeometryCollection *poGC = new OGRGeometryCollection();
490 810 : OGRMultiPolygon *poMP = nullptr;
491 810 : OGRPolygon *poLastPoly = nullptr;
492 4346 : for (int iPart = 0; iPart < nParts; iPart++)
493 : {
494 3536 : int nPartPoints = 0;
495 3536 : int nPartStart = 0;
496 :
497 : // Figure out details about this part's vertex list.
498 3536 : if (panPartStart == nullptr)
499 : {
500 0 : nPartPoints = nPoints;
501 : }
502 : else
503 : {
504 :
505 3536 : if (iPart == nParts - 1)
506 810 : nPartPoints = nPoints - panPartStart[iPart];
507 : else
508 2726 : nPartPoints = panPartStart[iPart + 1] - panPartStart[iPart];
509 3536 : nPartStart = panPartStart[iPart];
510 : }
511 :
512 3536 : OGRCreateFromMultiPatchPart(poGC, poMP, poLastPoly, panPartType[iPart],
513 3536 : nPartPoints, padfX + nPartStart,
514 3536 : padfY + nPartStart, padfZ + nPartStart);
515 : }
516 :
517 810 : if (poMP != nullptr && poLastPoly != nullptr)
518 : {
519 681 : poMP->addGeometryDirectly(poLastPoly);
520 : // poLastPoly = NULL;
521 : }
522 810 : if (poMP != nullptr)
523 681 : poGC->addGeometryDirectly(poMP);
524 :
525 810 : if (poGC->getNumGeometries() == 1)
526 : {
527 127 : OGRGeometry *poResultGeom = poGC->getGeometryRef(0);
528 127 : poGC->removeGeometry(0, FALSE);
529 127 : delete poGC;
530 127 : return poResultGeom;
531 : }
532 :
533 : else
534 683 : return poGC;
535 : }
536 :
537 : /************************************************************************/
538 : /* OGRWriteToShapeBin() */
539 : /* */
540 : /* Translate OGR geometry to a shapefile binary representation */
541 : /************************************************************************/
542 :
543 0 : OGRErr OGRWriteToShapeBin(const OGRGeometry *poGeom, GByte **ppabyShape,
544 : int *pnBytes)
545 : {
546 0 : int nShpSize = 4; // All types start with integer type number.
547 :
548 : /* -------------------------------------------------------------------- */
549 : /* Null or Empty input maps to SHPT_NULL. */
550 : /* -------------------------------------------------------------------- */
551 0 : if (!poGeom || poGeom->IsEmpty())
552 : {
553 0 : *ppabyShape = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nShpSize));
554 0 : if (*ppabyShape == nullptr)
555 0 : return OGRERR_FAILURE;
556 0 : GUInt32 zero = SHPT_NULL;
557 0 : memcpy(*ppabyShape, &zero, nShpSize);
558 0 : *pnBytes = nShpSize;
559 0 : return OGRERR_NONE;
560 : }
561 :
562 0 : const OGRwkbGeometryType nOGRType = wkbFlatten(poGeom->getGeometryType());
563 0 : const bool b3d = wkbHasZ(poGeom->getGeometryType());
564 0 : const bool bHasM = wkbHasM(poGeom->getGeometryType());
565 0 : const int nCoordDims = poGeom->CoordinateDimension();
566 :
567 0 : int nShpZSize = 0; // Z gets tacked onto the end.
568 0 : GUInt32 nPoints = 0;
569 0 : GUInt32 nParts = 0;
570 :
571 : /* -------------------------------------------------------------------- */
572 : /* Calculate the shape buffer size */
573 : /* -------------------------------------------------------------------- */
574 0 : if (nOGRType == wkbPoint)
575 : {
576 0 : nShpSize += 8 * nCoordDims;
577 : }
578 0 : else if (nOGRType == wkbLineString)
579 : {
580 0 : const OGRLineString *poLine = poGeom->toLineString();
581 0 : nPoints = poLine->getNumPoints();
582 0 : nParts = 1;
583 0 : nShpSize += 16 * nCoordDims; // xy(z)(m) box.
584 0 : nShpSize += 4; // nparts.
585 0 : nShpSize += 4; // npoints.
586 0 : nShpSize += 4; // Parts[1].
587 0 : nShpSize += 8 * nCoordDims * nPoints; // Points.
588 0 : nShpZSize = 16 + 8 * nPoints;
589 : }
590 0 : else if (nOGRType == wkbPolygon)
591 : {
592 0 : std::unique_ptr<OGRPolygon> poPoly(poGeom->toPolygon()->clone());
593 0 : poPoly->closeRings();
594 0 : nParts = poPoly->getNumInteriorRings() + 1;
595 0 : for (const auto poRing : *poPoly)
596 : {
597 0 : nPoints += poRing->getNumPoints();
598 : }
599 0 : nShpSize += 16 * nCoordDims; // xy(z)(m) box.
600 0 : nShpSize += 4; // nparts.
601 0 : nShpSize += 4; // npoints.
602 0 : nShpSize += 4 * nParts; // parts[nparts]
603 0 : nShpSize += 8 * nCoordDims * nPoints; // Points.
604 0 : nShpZSize = 16 + 8 * nPoints;
605 : }
606 0 : else if (nOGRType == wkbMultiPoint)
607 : {
608 0 : const OGRMultiPoint *poMPoint = poGeom->toMultiPoint();
609 0 : for (const auto poPoint : *poMPoint)
610 : {
611 0 : if (poPoint->IsEmpty())
612 0 : continue;
613 0 : nPoints++;
614 : }
615 0 : nShpSize += 16 * nCoordDims; // xy(z)(m) box.
616 0 : nShpSize += 4; // npoints.
617 0 : nShpSize += 8 * nCoordDims * nPoints; // Points.
618 0 : nShpZSize = 16 + 8 * nPoints;
619 : }
620 0 : else if (nOGRType == wkbMultiLineString)
621 : {
622 0 : const OGRMultiLineString *poMLine = poGeom->toMultiLineString();
623 0 : for (const auto poLine : *poMLine)
624 : {
625 : // Skip empties.
626 0 : if (poLine->IsEmpty())
627 0 : continue;
628 0 : nParts++;
629 0 : nPoints += poLine->getNumPoints();
630 : }
631 0 : nShpSize += 16 * nCoordDims; //* xy(z)(m) box.
632 0 : nShpSize += 4; // nparts.
633 0 : nShpSize += 4; // npoints.
634 0 : nShpSize += 4 * nParts; // parts[nparts].
635 0 : nShpSize += 8 * nCoordDims * nPoints; // Points.
636 0 : nShpZSize = 16 + 8 * nPoints;
637 : }
638 0 : else if (nOGRType == wkbMultiPolygon)
639 : {
640 : std::unique_ptr<OGRMultiPolygon> poMPoly(
641 0 : poGeom->toMultiPolygon()->clone());
642 0 : poMPoly->closeRings();
643 0 : for (const auto poPoly : *poMPoly)
644 : {
645 : // Skip empties.
646 0 : if (poPoly->IsEmpty())
647 0 : continue;
648 :
649 0 : const int nRings = poPoly->getNumInteriorRings() + 1;
650 0 : nParts += nRings;
651 0 : for (const auto poRing : *poPoly)
652 : {
653 0 : nPoints += poRing->getNumPoints();
654 : }
655 : }
656 0 : nShpSize += 16 * nCoordDims; // xy(z)(m) box.
657 0 : nShpSize += 4; // nparts.
658 0 : nShpSize += 4; // npoints.
659 0 : nShpSize += 4 * nParts; // parts[nparts].
660 0 : nShpSize += 8 * nCoordDims * nPoints; // Points.
661 0 : nShpZSize = 16 + 8 * nPoints;
662 : }
663 : else
664 : {
665 0 : return OGRERR_UNSUPPORTED_OPERATION;
666 : }
667 :
668 : // #define WRITE_ARC_HACK
669 : #ifdef WRITE_ARC_HACK
670 : int nShpSizeBeforeCurve = nShpSize;
671 : nShpSize += 4 + 4 + 4 + 20;
672 : #endif
673 : // Allocate our shape buffer.
674 0 : *ppabyShape = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nShpSize));
675 0 : if (!*ppabyShape)
676 0 : return OGRERR_FAILURE;
677 :
678 : #ifdef WRITE_ARC_HACK
679 : /* To be used with:
680 : id,WKT
681 : 1,"LINESTRING (1 0,0 1)"
682 : 2,"LINESTRING (5 1,6 0)"
683 : 3,"LINESTRING (10 1,11 0)"
684 : 4,"LINESTRING (16 0,15 1)"
685 : 5,"LINESTRING (21 0,20 1)"
686 : 6,"LINESTRING (31 0,30 2)" <-- not constant radius
687 : */
688 :
689 : GUInt32 nTmp = 1;
690 : memcpy((*ppabyShape) + nShpSizeBeforeCurve, &nTmp, 4);
691 : nTmp = 0;
692 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 4, &nTmp, 4);
693 : nTmp = EXT_SHAPE_SEGMENT_ARC;
694 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8, &nTmp, 4);
695 : static int nCounter = 0;
696 : nCounter++;
697 : if (nCounter == 1)
698 : {
699 : double dfVal = 0;
700 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
701 : dfVal = 0;
702 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
703 : nTmp = 0;
704 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
705 : }
706 : else if (nCounter == 2)
707 : {
708 : double dfVal = 5;
709 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
710 : dfVal = 0;
711 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
712 : nTmp = EXT_SHAPE_ARC_MINOR;
713 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
714 : }
715 : else if (nCounter == 3)
716 : {
717 : double dfVal = 10;
718 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
719 : dfVal = 0;
720 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
721 : nTmp = EXT_SHAPE_ARC_CCW;
722 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
723 : }
724 : else if (nCounter == 4)
725 : {
726 : double dfVal = 15;
727 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
728 : dfVal = 0;
729 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
730 : nTmp = EXT_SHAPE_ARC_CCW | EXT_SHAPE_ARC_MINOR;
731 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
732 : }
733 : else if (nCounter == 5)
734 : {
735 : double dfVal = 20;
736 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
737 : dfVal = 0;
738 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
739 : // Inconsistent with SP and EP. Only the CCW/not CCW is taken into
740 : // account by ArcGIS.
741 : nTmp = EXT_SHAPE_ARC_MINOR;
742 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
743 : }
744 : else if (nCounter == 6)
745 : {
746 : double dfVal = 30; // Radius inconsistent with SP and EP.
747 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4, &dfVal, 8);
748 : dfVal = 0;
749 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 8, &dfVal, 8);
750 : nTmp = EXT_SHAPE_ARC_MINOR;
751 : memcpy((*ppabyShape) + nShpSizeBeforeCurve + 8 + 4 + 16, &nTmp, 4);
752 : }
753 : #endif
754 :
755 : // Fill in the output size.
756 0 : *pnBytes = nShpSize;
757 :
758 : // Set up write pointers.
759 0 : unsigned char *pabyPtr = *ppabyShape;
760 0 : unsigned char *pabyPtrM = bHasM ? pabyPtr + nShpSize - nShpZSize : nullptr;
761 :
762 0 : unsigned char *pabyPtrZ = nullptr;
763 0 : if (b3d)
764 : {
765 0 : if (bHasM)
766 0 : pabyPtrZ = pabyPtrM - nShpZSize;
767 : else
768 0 : pabyPtrZ = pabyPtr + nShpSize - nShpZSize;
769 : }
770 :
771 : /* -------------------------------------------------------------------- */
772 : /* Write in the Shape type number now */
773 : /* -------------------------------------------------------------------- */
774 0 : GUInt32 nGType = SHPT_NULL;
775 :
776 0 : switch (nOGRType)
777 : {
778 0 : case wkbPoint:
779 : {
780 0 : nGType = (b3d && bHasM) ? SHPT_POINTZM
781 : : (b3d) ? SHPT_POINTZ
782 : : (bHasM) ? SHPT_POINTM
783 : : SHPT_POINT;
784 0 : break;
785 : }
786 0 : case wkbMultiPoint:
787 : {
788 0 : nGType = (b3d && bHasM) ? SHPT_MULTIPOINTZM
789 : : (b3d) ? SHPT_MULTIPOINTZ
790 : : (bHasM) ? SHPT_MULTIPOINTM
791 : : SHPT_MULTIPOINT;
792 0 : break;
793 : }
794 0 : case wkbLineString:
795 : case wkbMultiLineString:
796 : {
797 0 : nGType = (b3d && bHasM) ? SHPT_ARCZM
798 : : (b3d) ? SHPT_ARCZ
799 : : (bHasM) ? SHPT_ARCM
800 : : SHPT_ARC;
801 0 : break;
802 : }
803 0 : case wkbPolygon:
804 : case wkbMultiPolygon:
805 : {
806 0 : nGType = (b3d && bHasM) ? SHPT_POLYGONZM
807 : : (b3d) ? SHPT_POLYGONZ
808 : : (bHasM) ? SHPT_POLYGONM
809 : : SHPT_POLYGON;
810 0 : break;
811 : }
812 0 : default:
813 : {
814 0 : return OGRERR_UNSUPPORTED_OPERATION;
815 : }
816 : }
817 : // Write in the type number and advance the pointer.
818 : #ifdef WRITE_ARC_HACK
819 : nGType = SHPT_GENERALPOLYLINE | 0x20000000;
820 : #endif
821 :
822 0 : CPL_LSBPTR32(&nGType);
823 0 : memcpy(pabyPtr, &nGType, 4);
824 0 : pabyPtr += 4;
825 :
826 : /* -------------------------------------------------------------------- */
827 : /* POINT and POINTZ */
828 : /* -------------------------------------------------------------------- */
829 0 : if (nOGRType == wkbPoint)
830 : {
831 0 : auto poPoint = poGeom->toPoint();
832 0 : const double x = poPoint->getX();
833 0 : const double y = poPoint->getY();
834 :
835 : // Copy in the raw data.
836 0 : memcpy(pabyPtr, &x, 8);
837 0 : memcpy(pabyPtr + 8, &y, 8);
838 0 : if (b3d)
839 : {
840 0 : double z = poPoint->getZ();
841 0 : memcpy(pabyPtr + 8 + 8, &z, 8);
842 : }
843 0 : if (bHasM)
844 : {
845 0 : double m = poPoint->getM();
846 0 : memcpy(pabyPtr + 8 + ((b3d) ? 16 : 8), &m, 8);
847 : }
848 :
849 : // Swap if needed. Shape doubles always LSB.
850 : if constexpr (OGR_SWAP(wkbNDR))
851 : {
852 : CPL_SWAPDOUBLE(pabyPtr);
853 : CPL_SWAPDOUBLE(pabyPtr + 8);
854 : if (b3d)
855 : CPL_SWAPDOUBLE(pabyPtr + 8 + 8);
856 : if (bHasM)
857 : CPL_SWAPDOUBLE(pabyPtr + 8 + ((b3d) ? 16 : 8));
858 : }
859 :
860 0 : return OGRERR_NONE;
861 : }
862 :
863 : /* -------------------------------------------------------------------- */
864 : /* All the non-POINT types require an envelope next */
865 : /* -------------------------------------------------------------------- */
866 0 : OGREnvelope3D envelope;
867 0 : poGeom->getEnvelope(&envelope);
868 0 : memcpy(pabyPtr, &(envelope.MinX), 8);
869 0 : memcpy(pabyPtr + 8, &(envelope.MinY), 8);
870 0 : memcpy(pabyPtr + 8 + 8, &(envelope.MaxX), 8);
871 0 : memcpy(pabyPtr + 8 + 8 + 8, &(envelope.MaxY), 8);
872 :
873 : // Swap box if needed. Shape doubles are always LSB.
874 : if constexpr (OGR_SWAP(wkbNDR))
875 : {
876 : for (int i = 0; i < 4; i++)
877 : CPL_SWAPDOUBLE(pabyPtr + 8 * i);
878 : }
879 0 : pabyPtr += 32;
880 :
881 : // Write in the Z bounds at the end of the XY buffer.
882 0 : if (b3d)
883 : {
884 0 : memcpy(pabyPtrZ, &(envelope.MinZ), 8);
885 0 : memcpy(pabyPtrZ + 8, &(envelope.MaxZ), 8);
886 :
887 : // Swap Z bounds if necessary.
888 : if constexpr (OGR_SWAP(wkbNDR))
889 : {
890 : for (int i = 0; i < 2; i++)
891 : CPL_SWAPDOUBLE(pabyPtrZ + 8 * i);
892 : }
893 0 : pabyPtrZ += 16;
894 : }
895 :
896 : // Reserve space for the M bounds at the end of the XY buffer.
897 0 : GByte *pabyPtrMBounds = nullptr;
898 0 : double dfMinM = std::numeric_limits<double>::max();
899 0 : double dfMaxM = -dfMinM;
900 0 : if (bHasM)
901 : {
902 0 : pabyPtrMBounds = pabyPtrM;
903 0 : pabyPtrM += 16;
904 : }
905 :
906 : /* -------------------------------------------------------------------- */
907 : /* LINESTRING and LINESTRINGZ */
908 : /* -------------------------------------------------------------------- */
909 0 : if (nOGRType == wkbLineString)
910 : {
911 0 : auto poLine = poGeom->toLineString();
912 :
913 : // Write in the nparts (1).
914 0 : const GUInt32 nPartsLsb = CPL_LSBWORD32(nParts);
915 0 : memcpy(pabyPtr, &nPartsLsb, 4);
916 0 : pabyPtr += 4;
917 :
918 : // Write in the npoints.
919 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
920 0 : memcpy(pabyPtr, &nPointsLsb, 4);
921 0 : pabyPtr += 4;
922 :
923 : // Write in the part index (0).
924 0 : GUInt32 nPartIndex = 0;
925 0 : memcpy(pabyPtr, &nPartIndex, 4);
926 0 : pabyPtr += 4;
927 :
928 : // Write in the point data.
929 0 : poLine->getPoints(reinterpret_cast<OGRRawPoint *>(pabyPtr),
930 : reinterpret_cast<double *>(pabyPtrZ));
931 0 : if (bHasM)
932 : {
933 0 : for (GUInt32 k = 0; k < nPoints; k++)
934 : {
935 0 : double dfM = poLine->getM(k);
936 0 : memcpy(pabyPtrM + 8 * k, &dfM, 8);
937 0 : if (dfM < dfMinM)
938 0 : dfMinM = dfM;
939 0 : if (dfM > dfMaxM)
940 0 : dfMaxM = dfM;
941 : }
942 : }
943 :
944 : // Swap if necessary.
945 : if constexpr (OGR_SWAP(wkbNDR))
946 : {
947 : for (GUInt32 k = 0; k < nPoints; k++)
948 : {
949 : CPL_SWAPDOUBLE(pabyPtr + 16 * k);
950 : CPL_SWAPDOUBLE(pabyPtr + 16 * k + 8);
951 : if (b3d)
952 : CPL_SWAPDOUBLE(pabyPtrZ + 8 * k);
953 : if (bHasM)
954 : CPL_SWAPDOUBLE(pabyPtrM + 8 * k);
955 : }
956 : }
957 : }
958 : /* -------------------------------------------------------------------- */
959 : /* POLYGON and POLYGONZ */
960 : /* -------------------------------------------------------------------- */
961 0 : else if (nOGRType == wkbPolygon)
962 : {
963 0 : auto poPoly = poGeom->toPolygon();
964 :
965 : // Write in the part count.
966 0 : GUInt32 nPartsLsb = CPL_LSBWORD32(nParts);
967 0 : memcpy(pabyPtr, &nPartsLsb, 4);
968 0 : pabyPtr += 4;
969 :
970 : // Write in the total point count.
971 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
972 0 : memcpy(pabyPtr, &nPointsLsb, 4);
973 0 : pabyPtr += 4;
974 :
975 : /* --------------------------------------------------------------------
976 : */
977 : /* Now we have to visit each ring and write an index number into */
978 : /* the parts list, and the coordinates into the points list. */
979 : /* to do it in one pass, we will use three write pointers. */
980 : /* pabyPtr writes the part indexes */
981 : /* pabyPoints writes the xy coordinates */
982 : /* pabyPtrZ writes the z coordinates */
983 : /* --------------------------------------------------------------------
984 : */
985 :
986 : // Just past the partindex[nparts] array.
987 0 : unsigned char *pabyPoints = pabyPtr + 4 * nParts;
988 :
989 0 : int nPointIndexCount = 0;
990 :
991 0 : for (GUInt32 i = 0; i < nParts; i++)
992 : {
993 : // Check our Ring and condition it.
994 0 : std::unique_ptr<OGRLinearRing> poRing;
995 0 : if (i == 0)
996 : {
997 0 : poRing.reset(poPoly->getExteriorRing()->clone());
998 0 : assert(poRing);
999 : // Outer ring must be clockwise.
1000 0 : if (!poRing->isClockwise())
1001 0 : poRing->reversePoints();
1002 : }
1003 : else
1004 : {
1005 0 : poRing.reset(poPoly->getInteriorRing(i - 1)->clone());
1006 0 : assert(poRing);
1007 : // Inner rings should be anti-clockwise.
1008 0 : if (poRing->isClockwise())
1009 0 : poRing->reversePoints();
1010 : }
1011 :
1012 0 : int nRingNumPoints = poRing->getNumPoints();
1013 :
1014 : #ifndef WRITE_ARC_HACK
1015 : // Cannot write un-closed rings to shape.
1016 0 : if (nRingNumPoints <= 2 || !poRing->get_IsClosed())
1017 0 : return OGRERR_FAILURE;
1018 : #endif
1019 :
1020 : // Write in the part index.
1021 0 : GUInt32 nPartIndex = CPL_LSBWORD32(nPointIndexCount);
1022 0 : memcpy(pabyPtr, &nPartIndex, 4);
1023 :
1024 : // Write in the point data.
1025 0 : poRing->getPoints(reinterpret_cast<OGRRawPoint *>(pabyPoints),
1026 : reinterpret_cast<double *>(pabyPtrZ));
1027 0 : if (bHasM)
1028 : {
1029 0 : for (int k = 0; k < nRingNumPoints; k++)
1030 : {
1031 0 : double dfM = poRing->getM(k);
1032 0 : memcpy(pabyPtrM + 8 * k, &dfM, 8);
1033 0 : if (dfM < dfMinM)
1034 0 : dfMinM = dfM;
1035 0 : if (dfM > dfMaxM)
1036 0 : dfMaxM = dfM;
1037 : }
1038 : }
1039 :
1040 : // Swap if necessary.
1041 : if constexpr (OGR_SWAP(wkbNDR))
1042 : {
1043 : for (int k = 0; k < nRingNumPoints; k++)
1044 : {
1045 : CPL_SWAPDOUBLE(pabyPoints + 16 * k);
1046 : CPL_SWAPDOUBLE(pabyPoints + 16 * k + 8);
1047 : if (b3d)
1048 : CPL_SWAPDOUBLE(pabyPtrZ + 8 * k);
1049 : if (bHasM)
1050 : CPL_SWAPDOUBLE(pabyPtrM + 8 * k);
1051 : }
1052 : }
1053 :
1054 0 : nPointIndexCount += nRingNumPoints;
1055 : // Advance the write pointers.
1056 0 : pabyPtr += 4;
1057 0 : pabyPoints += 16 * nRingNumPoints;
1058 0 : if (b3d)
1059 0 : pabyPtrZ += 8 * nRingNumPoints;
1060 0 : if (bHasM)
1061 0 : pabyPtrM += 8 * nRingNumPoints;
1062 : }
1063 : }
1064 : /* -------------------------------------------------------------------- */
1065 : /* MULTIPOINT and MULTIPOINTZ */
1066 : /* -------------------------------------------------------------------- */
1067 0 : else if (nOGRType == wkbMultiPoint)
1068 : {
1069 0 : auto poMPoint = poGeom->toMultiPoint();
1070 :
1071 : // Write in the total point count.
1072 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
1073 0 : memcpy(pabyPtr, &nPointsLsb, 4);
1074 0 : pabyPtr += 4;
1075 :
1076 : /* --------------------------------------------------------------------
1077 : */
1078 : /* Now we have to visit each point write it into the points list */
1079 : /* We will use two write pointers. */
1080 : /* pabyPtr writes the xy coordinates */
1081 : /* pabyPtrZ writes the z coordinates */
1082 : /* --------------------------------------------------------------------
1083 : */
1084 :
1085 0 : for (auto &&poPt : poMPoint)
1086 : {
1087 : // Skip empties.
1088 0 : if (poPt->IsEmpty())
1089 0 : continue;
1090 :
1091 : // Write the coordinates.
1092 0 : double x = poPt->getX();
1093 0 : double y = poPt->getY();
1094 0 : memcpy(pabyPtr, &x, 8);
1095 0 : memcpy(pabyPtr + 8, &y, 8);
1096 0 : if (b3d)
1097 : {
1098 0 : double z = poPt->getZ();
1099 0 : memcpy(pabyPtrZ, &z, 8);
1100 : }
1101 0 : if (bHasM)
1102 : {
1103 0 : double dfM = poPt->getM();
1104 0 : memcpy(pabyPtrM, &dfM, 8);
1105 0 : if (dfM < dfMinM)
1106 0 : dfMinM = dfM;
1107 0 : if (dfM > dfMaxM)
1108 0 : dfMaxM = dfM;
1109 : }
1110 :
1111 : // Swap if necessary.
1112 : if constexpr (OGR_SWAP(wkbNDR))
1113 : {
1114 : CPL_SWAPDOUBLE(pabyPtr);
1115 : CPL_SWAPDOUBLE(pabyPtr + 8);
1116 : if (b3d)
1117 : CPL_SWAPDOUBLE(pabyPtrZ);
1118 : if (bHasM)
1119 : CPL_SWAPDOUBLE(pabyPtrM);
1120 : }
1121 :
1122 : // Advance the write pointers.
1123 0 : pabyPtr += 16;
1124 0 : if (b3d)
1125 0 : pabyPtrZ += 8;
1126 0 : if (bHasM)
1127 0 : pabyPtrM += 8;
1128 : }
1129 : }
1130 :
1131 : /* -------------------------------------------------------------------- */
1132 : /* MULTILINESTRING and MULTILINESTRINGZ */
1133 : /* -------------------------------------------------------------------- */
1134 0 : else if (nOGRType == wkbMultiLineString)
1135 : {
1136 0 : auto poMLine = poGeom->toMultiLineString();
1137 :
1138 : // Write in the part count.
1139 0 : GUInt32 nPartsLsb = CPL_LSBWORD32(nParts);
1140 0 : memcpy(pabyPtr, &nPartsLsb, 4);
1141 0 : pabyPtr += 4;
1142 :
1143 : // Write in the total point count.
1144 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
1145 0 : memcpy(pabyPtr, &nPointsLsb, 4);
1146 0 : pabyPtr += 4;
1147 :
1148 : // Just past the partindex[nparts] array.
1149 0 : unsigned char *pabyPoints = pabyPtr + 4 * nParts;
1150 :
1151 0 : int nPointIndexCount = 0;
1152 :
1153 0 : for (auto &&poLine : poMLine)
1154 : {
1155 : // Skip empties.
1156 0 : if (poLine->IsEmpty())
1157 0 : continue;
1158 :
1159 0 : int nLineNumPoints = poLine->getNumPoints();
1160 :
1161 : // Write in the part index.
1162 0 : GUInt32 nPartIndex = CPL_LSBWORD32(nPointIndexCount);
1163 0 : memcpy(pabyPtr, &nPartIndex, 4);
1164 :
1165 : // Write in the point data.
1166 0 : poLine->getPoints(reinterpret_cast<OGRRawPoint *>(pabyPoints),
1167 : reinterpret_cast<double *>(pabyPtrZ));
1168 0 : if (bHasM)
1169 : {
1170 0 : for (int k = 0; k < nLineNumPoints; k++)
1171 : {
1172 0 : double dfM = poLine->getM(k);
1173 0 : memcpy(pabyPtrM + 8 * k, &dfM, 8);
1174 0 : if (dfM < dfMinM)
1175 0 : dfMinM = dfM;
1176 0 : if (dfM > dfMaxM)
1177 0 : dfMaxM = dfM;
1178 : }
1179 : }
1180 :
1181 : // Swap if necessary.
1182 : if constexpr (OGR_SWAP(wkbNDR))
1183 : {
1184 : for (int k = 0; k < nLineNumPoints; k++)
1185 : {
1186 : CPL_SWAPDOUBLE(pabyPoints + 16 * k);
1187 : CPL_SWAPDOUBLE(pabyPoints + 16 * k + 8);
1188 : if (b3d)
1189 : CPL_SWAPDOUBLE(pabyPtrZ + 8 * k);
1190 : if (bHasM)
1191 : CPL_SWAPDOUBLE(pabyPtrM + 8 * k);
1192 : }
1193 : }
1194 :
1195 0 : nPointIndexCount += nLineNumPoints;
1196 :
1197 : // Advance the write pointers.
1198 0 : pabyPtr += 4;
1199 0 : pabyPoints += 16 * nLineNumPoints;
1200 0 : if (b3d)
1201 0 : pabyPtrZ += 8 * nLineNumPoints;
1202 0 : if (bHasM)
1203 0 : pabyPtrM += 8 * nLineNumPoints;
1204 : }
1205 : }
1206 : /* -------------------------------------------------------------------- */
1207 : /* MULTIPOLYGON and MULTIPOLYGONZ */
1208 : /* -------------------------------------------------------------------- */
1209 : else // if( nOGRType == wkbMultiPolygon )
1210 : {
1211 0 : auto poMPoly = poGeom->toMultiPolygon();
1212 :
1213 : // Write in the part count.
1214 0 : GUInt32 nPartsLsb = CPL_LSBWORD32(nParts);
1215 0 : memcpy(pabyPtr, &nPartsLsb, 4);
1216 0 : pabyPtr += 4;
1217 :
1218 : // Write in the total point count.
1219 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
1220 0 : memcpy(pabyPtr, &nPointsLsb, 4);
1221 0 : pabyPtr += 4;
1222 :
1223 : /* --------------------------------------------------------------------
1224 : */
1225 : /* Now we have to visit each ring and write an index number into */
1226 : /* the parts list, and the coordinates into the points list. */
1227 : /* to do it in one pass, we will use three write pointers. */
1228 : /* pabyPtr writes the part indexes */
1229 : /* pabyPoints writes the xy coordinates */
1230 : /* pabyPtrZ writes the z coordinates */
1231 : /* --------------------------------------------------------------------
1232 : */
1233 :
1234 : // Just past the partindex[nparts] array.
1235 0 : unsigned char *pabyPoints = pabyPtr + 4 * nParts;
1236 :
1237 0 : int nPointIndexCount = 0;
1238 :
1239 0 : for (auto &&poPoly : poMPoly)
1240 : {
1241 : // Skip empties.
1242 0 : if (poPoly->IsEmpty())
1243 0 : continue;
1244 :
1245 0 : int nRings = 1 + poPoly->getNumInteriorRings();
1246 :
1247 0 : for (int j = 0; j < nRings; j++)
1248 : {
1249 : // Check our Ring and condition it.
1250 0 : std::unique_ptr<OGRLinearRing> poRing;
1251 0 : if (j == 0)
1252 : {
1253 0 : poRing.reset(poPoly->getExteriorRing()->clone());
1254 0 : assert(poRing != nullptr);
1255 : // Outer ring must be clockwise.
1256 0 : if (!poRing->isClockwise())
1257 0 : poRing->reversePoints();
1258 : }
1259 : else
1260 : {
1261 0 : poRing.reset(poPoly->getInteriorRing(j - 1)->clone());
1262 0 : assert(poRing != nullptr);
1263 : // Inner rings should be anti-clockwise.
1264 0 : if (poRing->isClockwise())
1265 0 : poRing->reversePoints();
1266 : }
1267 :
1268 0 : int nRingNumPoints = poRing->getNumPoints();
1269 :
1270 : // Cannot write closed rings to shape.
1271 0 : if (nRingNumPoints <= 2 || !poRing->get_IsClosed())
1272 0 : return OGRERR_FAILURE;
1273 :
1274 : // Write in the part index.
1275 0 : GUInt32 nPartIndex = CPL_LSBWORD32(nPointIndexCount);
1276 0 : memcpy(pabyPtr, &nPartIndex, 4);
1277 :
1278 : // Write in the point data.
1279 0 : poRing->getPoints(reinterpret_cast<OGRRawPoint *>(pabyPoints),
1280 : reinterpret_cast<double *>(pabyPtrZ));
1281 0 : if (bHasM)
1282 : {
1283 0 : for (int k = 0; k < nRingNumPoints; k++)
1284 : {
1285 0 : double dfM = poRing->getM(k);
1286 0 : memcpy(pabyPtrM + 8 * k, &dfM, 8);
1287 0 : if (dfM < dfMinM)
1288 0 : dfMinM = dfM;
1289 0 : if (dfM > dfMaxM)
1290 0 : dfMaxM = dfM;
1291 : }
1292 : }
1293 :
1294 : // Swap if necessary.
1295 : if constexpr (OGR_SWAP(wkbNDR))
1296 : {
1297 : for (int k = 0; k < nRingNumPoints; k++)
1298 : {
1299 : CPL_SWAPDOUBLE(pabyPoints + 16 * k);
1300 : CPL_SWAPDOUBLE(pabyPoints + 16 * k + 8);
1301 : if (b3d)
1302 : CPL_SWAPDOUBLE(pabyPtrZ + 8 * k);
1303 : if (bHasM)
1304 : CPL_SWAPDOUBLE(pabyPtrM + 8 * k);
1305 : }
1306 : }
1307 :
1308 0 : nPointIndexCount += nRingNumPoints;
1309 : // Advance the write pointers.
1310 0 : pabyPtr += 4;
1311 0 : pabyPoints += 16 * nRingNumPoints;
1312 0 : if (b3d)
1313 0 : pabyPtrZ += 8 * nRingNumPoints;
1314 0 : if (bHasM)
1315 0 : pabyPtrM += 8 * nRingNumPoints;
1316 : }
1317 : }
1318 : }
1319 :
1320 0 : if (bHasM)
1321 : {
1322 0 : if (dfMinM > dfMaxM)
1323 : {
1324 0 : dfMinM = 0.0;
1325 0 : dfMaxM = 0.0;
1326 : }
1327 0 : memcpy(pabyPtrMBounds, &(dfMinM), 8);
1328 0 : memcpy(pabyPtrMBounds + 8, &(dfMaxM), 8);
1329 :
1330 : // Swap M bounds if necessary.
1331 : if constexpr (OGR_SWAP(wkbNDR))
1332 : {
1333 : for (int i = 0; i < 2; i++)
1334 : CPL_SWAPDOUBLE(pabyPtrMBounds + 8 * i);
1335 : }
1336 : }
1337 :
1338 0 : return OGRERR_NONE;
1339 : }
1340 :
1341 : /************************************************************************/
1342 : /* OGRCreateMultiPatch() */
1343 : /************************************************************************/
1344 :
1345 36 : OGRErr OGRCreateMultiPatch(const OGRGeometry *poGeomConst,
1346 : int bAllowSHPTTriangle, int &nParts,
1347 : std::vector<int> &anPartStart,
1348 : std::vector<int> &anPartType, int &nPoints,
1349 : std::vector<OGRRawPoint> &aoPoints,
1350 : std::vector<double> &adfZ)
1351 : {
1352 36 : const OGRwkbGeometryType eType = wkbFlatten(poGeomConst->getGeometryType());
1353 36 : if (eType != wkbPolygon && eType != wkbTriangle &&
1354 32 : eType != wkbMultiPolygon && eType != wkbMultiSurface &&
1355 16 : eType != wkbTIN && eType != wkbPolyhedralSurface &&
1356 : eType != wkbGeometryCollection)
1357 : {
1358 0 : return OGRERR_UNSUPPORTED_OPERATION;
1359 : }
1360 :
1361 72 : std::unique_ptr<OGRGeometry> poGeom(poGeomConst->clone());
1362 36 : poGeom->closeRings();
1363 :
1364 36 : OGRMultiPolygon *poMPoly = nullptr;
1365 36 : std::unique_ptr<OGRGeometry> poGeomToDelete;
1366 36 : if (eType == wkbMultiPolygon)
1367 1 : poMPoly = poGeom->toMultiPolygon();
1368 : else
1369 : {
1370 70 : poGeomToDelete = std::unique_ptr<OGRGeometry>(
1371 70 : OGRGeometryFactory::forceToMultiPolygon(poGeom->clone()));
1372 70 : if (poGeomToDelete.get() &&
1373 35 : wkbFlatten(poGeomToDelete->getGeometryType()) == wkbMultiPolygon)
1374 : {
1375 32 : poMPoly = poGeomToDelete->toMultiPolygon();
1376 : }
1377 : }
1378 36 : if (poMPoly == nullptr)
1379 : {
1380 3 : return OGRERR_UNSUPPORTED_OPERATION;
1381 : }
1382 :
1383 33 : nParts = 0;
1384 33 : anPartStart.clear();
1385 33 : anPartType.clear();
1386 33 : nPoints = 0;
1387 33 : aoPoints.clear();
1388 33 : adfZ.clear();
1389 33 : int nBeginLastPart = 0;
1390 133 : for (const auto poPoly : *poMPoly)
1391 : {
1392 : // Skip empties.
1393 100 : if (poPoly->IsEmpty())
1394 0 : continue;
1395 :
1396 100 : const int nRings = poPoly->getNumInteriorRings() + 1;
1397 100 : const OGRLinearRing *poRing = poPoly->getExteriorRing();
1398 100 : if (nRings == 1 && poRing->getNumPoints() == 4)
1399 : {
1400 90 : int nCorrectedPoints = nPoints;
1401 95 : if (nParts > 0 && anPartType[nParts - 1] == SHPP_OUTERRING &&
1402 5 : nPoints - anPartStart[nParts - 1] == 4)
1403 : {
1404 5 : nCorrectedPoints--;
1405 : }
1406 :
1407 237 : if (nParts > 0 &&
1408 57 : ((anPartType[nParts - 1] == SHPP_TRIANGLES &&
1409 31 : nPoints - anPartStart[nParts - 1] == 3) ||
1410 26 : (anPartType[nParts - 1] == SHPP_OUTERRING &&
1411 5 : nPoints - anPartStart[nParts - 1] == 4) ||
1412 21 : anPartType[nParts - 1] == SHPP_TRIFAN) &&
1413 47 : poRing->getX(0) == aoPoints[nBeginLastPart].x &&
1414 34 : poRing->getY(0) == aoPoints[nBeginLastPart].y &&
1415 13 : poRing->getZ(0) == adfZ[nBeginLastPart] &&
1416 13 : poRing->getX(1) == aoPoints[nCorrectedPoints - 1].x &&
1417 159 : poRing->getY(1) == aoPoints[nCorrectedPoints - 1].y &&
1418 12 : poRing->getZ(1) == adfZ[nCorrectedPoints - 1])
1419 : {
1420 12 : nPoints = nCorrectedPoints;
1421 12 : anPartType[nParts - 1] = SHPP_TRIFAN;
1422 :
1423 12 : aoPoints.resize(nPoints + 1);
1424 12 : adfZ.resize(nPoints + 1);
1425 12 : aoPoints[nPoints].x = poRing->getX(2);
1426 12 : aoPoints[nPoints].y = poRing->getY(2);
1427 12 : adfZ[nPoints] = poRing->getZ(2);
1428 12 : nPoints++;
1429 : }
1430 201 : else if (nParts > 0 &&
1431 45 : ((anPartType[nParts - 1] == SHPP_TRIANGLES &&
1432 20 : nPoints - anPartStart[nParts - 1] == 3) ||
1433 25 : (anPartType[nParts - 1] == SHPP_OUTERRING &&
1434 4 : nPoints - anPartStart[nParts - 1] == 4) ||
1435 21 : anPartType[nParts - 1] == SHPP_TRISTRIP) &&
1436 34 : poRing->getX(0) == aoPoints[nCorrectedPoints - 2].x &&
1437 22 : poRing->getY(0) == aoPoints[nCorrectedPoints - 2].y &&
1438 21 : poRing->getZ(0) == adfZ[nCorrectedPoints - 2] &&
1439 21 : poRing->getX(1) == aoPoints[nCorrectedPoints - 1].x &&
1440 144 : poRing->getY(1) == aoPoints[nCorrectedPoints - 1].y &&
1441 21 : poRing->getZ(1) == adfZ[nCorrectedPoints - 1])
1442 : {
1443 21 : nPoints = nCorrectedPoints;
1444 21 : anPartType[nParts - 1] = SHPP_TRISTRIP;
1445 :
1446 21 : aoPoints.resize(nPoints + 1);
1447 21 : adfZ.resize(nPoints + 1);
1448 21 : aoPoints[nPoints].x = poRing->getX(2);
1449 21 : aoPoints[nPoints].y = poRing->getY(2);
1450 21 : adfZ[nPoints] = poRing->getZ(2);
1451 21 : nPoints++;
1452 : }
1453 : else
1454 : {
1455 57 : if (nParts == 0 || anPartType[nParts - 1] != SHPP_TRIANGLES ||
1456 : !bAllowSHPTTriangle)
1457 : {
1458 57 : nBeginLastPart = nPoints;
1459 :
1460 57 : anPartStart.resize(nParts + 1);
1461 57 : anPartType.resize(nParts + 1);
1462 57 : anPartStart[nParts] = nPoints;
1463 57 : anPartType[nParts] =
1464 57 : bAllowSHPTTriangle ? SHPP_TRIANGLES : SHPP_OUTERRING;
1465 57 : nParts++;
1466 : }
1467 :
1468 57 : aoPoints.resize(nPoints + 4);
1469 57 : adfZ.resize(nPoints + 4);
1470 285 : for (int i = 0; i < 4; i++)
1471 : {
1472 228 : aoPoints[nPoints + i].x = poRing->getX(i);
1473 228 : aoPoints[nPoints + i].y = poRing->getY(i);
1474 228 : adfZ[nPoints + i] = poRing->getZ(i);
1475 : }
1476 57 : nPoints += bAllowSHPTTriangle ? 3 : 4;
1477 : }
1478 : }
1479 : else
1480 : {
1481 10 : anPartStart.resize(nParts + nRings);
1482 10 : anPartType.resize(nParts + nRings);
1483 :
1484 30 : for (int i = 0; i < nRings; i++)
1485 : {
1486 20 : anPartStart[nParts + i] = nPoints;
1487 20 : if (i == 0)
1488 : {
1489 10 : poRing = poPoly->getExteriorRing();
1490 10 : anPartType[nParts + i] = SHPP_OUTERRING;
1491 : }
1492 : else
1493 : {
1494 10 : poRing = poPoly->getInteriorRing(i - 1);
1495 10 : anPartType[nParts + i] = SHPP_INNERRING;
1496 : }
1497 20 : aoPoints.resize(nPoints + poRing->getNumPoints());
1498 20 : adfZ.resize(nPoints + poRing->getNumPoints());
1499 120 : for (int k = 0; k < poRing->getNumPoints(); k++)
1500 : {
1501 100 : aoPoints[nPoints + k].x = poRing->getX(k);
1502 100 : aoPoints[nPoints + k].y = poRing->getY(k);
1503 100 : adfZ[nPoints + k] = poRing->getZ(k);
1504 : }
1505 20 : nPoints += poRing->getNumPoints();
1506 : }
1507 :
1508 10 : nParts += nRings;
1509 : }
1510 : }
1511 :
1512 33 : if (nParts == 1 && anPartType[0] == SHPP_OUTERRING && nPoints == 4)
1513 : {
1514 8 : anPartType[0] = SHPP_TRIFAN;
1515 8 : nPoints = 3;
1516 : }
1517 :
1518 33 : return OGRERR_NONE;
1519 : }
1520 :
1521 : /************************************************************************/
1522 : /* OGRWriteMultiPatchToShapeBin() */
1523 : /************************************************************************/
1524 :
1525 0 : OGRErr OGRWriteMultiPatchToShapeBin(const OGRGeometry *poGeom,
1526 : GByte **ppabyShape, int *pnBytes)
1527 : {
1528 0 : int nParts = 0;
1529 0 : std::vector<int> anPartStart;
1530 0 : std::vector<int> anPartType;
1531 0 : int nPoints = 0;
1532 0 : std::vector<OGRRawPoint> aoPoints;
1533 0 : std::vector<double> adfZ;
1534 0 : OGRErr eErr = OGRCreateMultiPatch(poGeom, TRUE, nParts, anPartStart,
1535 : anPartType, nPoints, aoPoints, adfZ);
1536 0 : if (eErr != OGRERR_NONE)
1537 0 : return eErr;
1538 :
1539 : const bool bOmitZ =
1540 0 : !poGeom->Is3D() &&
1541 0 : CPLTestBool(CPLGetConfigOption("OGR_MULTIPATCH_OMIT_Z", "NO"));
1542 :
1543 0 : int nShpSize = 4; // All types start with integer type number.
1544 0 : nShpSize += 16 * 2; // xy bbox.
1545 0 : nShpSize += 4; // nparts.
1546 0 : nShpSize += 4; // npoints.
1547 0 : nShpSize += 4 * nParts; // panPartStart[nparts].
1548 0 : nShpSize += 4 * nParts; // panPartType[nparts].
1549 0 : nShpSize += 8 * 2 * nPoints; // xy points.
1550 0 : if (!bOmitZ)
1551 : {
1552 0 : nShpSize += 16; // z bbox.
1553 0 : nShpSize += 8 * nPoints; // z points.
1554 : }
1555 :
1556 0 : *pnBytes = nShpSize;
1557 0 : *ppabyShape = static_cast<GByte *>(CPLMalloc(nShpSize));
1558 :
1559 0 : GByte *pabyPtr = *ppabyShape;
1560 :
1561 : // Write in the type number and advance the pointer.
1562 0 : GUInt32 nGType = bOmitZ ? CPL_LSBWORD32(SHPT_GENERALMULTIPATCH)
1563 : : CPL_LSBWORD32(SHPT_MULTIPATCH);
1564 0 : memcpy(pabyPtr, &nGType, 4);
1565 0 : pabyPtr += 4;
1566 :
1567 0 : OGREnvelope3D envelope;
1568 0 : poGeom->getEnvelope(&envelope);
1569 0 : memcpy(pabyPtr, &(envelope.MinX), 8);
1570 0 : memcpy(pabyPtr + 8, &(envelope.MinY), 8);
1571 0 : memcpy(pabyPtr + 8 + 8, &(envelope.MaxX), 8);
1572 0 : memcpy(pabyPtr + 8 + 8 + 8, &(envelope.MaxY), 8);
1573 :
1574 : // Swap box if needed. Shape doubles are always LSB.
1575 : if constexpr (OGR_SWAP(wkbNDR))
1576 : {
1577 : for (int i = 0; i < 4; i++)
1578 : CPL_SWAPDOUBLE(pabyPtr + 8 * i);
1579 : }
1580 0 : pabyPtr += 32;
1581 :
1582 : // Write in the part count.
1583 0 : GUInt32 nPartsLsb = CPL_LSBWORD32(nParts);
1584 0 : memcpy(pabyPtr, &nPartsLsb, 4);
1585 0 : pabyPtr += 4;
1586 :
1587 : // Write in the total point count.
1588 0 : GUInt32 nPointsLsb = CPL_LSBWORD32(nPoints);
1589 0 : memcpy(pabyPtr, &nPointsLsb, 4);
1590 0 : pabyPtr += 4;
1591 :
1592 0 : for (int i = 0; i < nParts; i++)
1593 : {
1594 0 : int nPartStart = CPL_LSBWORD32(anPartStart[i]);
1595 0 : memcpy(pabyPtr, &nPartStart, 4);
1596 0 : pabyPtr += 4;
1597 : }
1598 0 : for (int i = 0; i < nParts; i++)
1599 : {
1600 0 : int nPartType = CPL_LSBWORD32(anPartType[i]);
1601 0 : memcpy(pabyPtr, &nPartType, 4);
1602 0 : pabyPtr += 4;
1603 : }
1604 :
1605 0 : if (!aoPoints.empty())
1606 0 : memcpy(pabyPtr, aoPoints.data(), 2 * 8 * nPoints);
1607 :
1608 : // Swap box if needed. Shape doubles are always LSB.
1609 : if constexpr (OGR_SWAP(wkbNDR))
1610 : {
1611 : for (int i = 0; i < 2 * nPoints; i++)
1612 : CPL_SWAPDOUBLE(pabyPtr + 8 * i);
1613 : }
1614 0 : pabyPtr += 2 * 8 * nPoints;
1615 :
1616 0 : if (!bOmitZ)
1617 : {
1618 0 : memcpy(pabyPtr, &(envelope.MinZ), 8);
1619 0 : memcpy(pabyPtr + 8, &(envelope.MaxZ), 8);
1620 : if constexpr (OGR_SWAP(wkbNDR))
1621 : {
1622 : for (int i = 0; i < 2; i++)
1623 : CPL_SWAPDOUBLE(pabyPtr + 8 * i);
1624 : }
1625 0 : pabyPtr += 16;
1626 :
1627 0 : if (!adfZ.empty())
1628 0 : memcpy(pabyPtr, adfZ.data(), 8 * nPoints);
1629 : // Swap box if needed. Shape doubles are always LSB.
1630 : if constexpr (OGR_SWAP(wkbNDR))
1631 : {
1632 : for (int i = 0; i < nPoints; i++)
1633 : CPL_SWAPDOUBLE(pabyPtr + 8 * i);
1634 : }
1635 : // pabyPtr += 8 * nPoints;
1636 : }
1637 :
1638 0 : return OGRERR_NONE;
1639 : }
1640 :
1641 : /************************************************************************/
1642 : /* GetAngleOnEllipse() */
1643 : /************************************************************************/
1644 :
1645 : // Return the angle in deg [0, 360] of dfArcX,dfArcY regarding the
1646 : // ellipse semi-major axis.
1647 26 : static double GetAngleOnEllipse(double dfPointOnArcX, double dfPointOnArcY,
1648 : double dfCenterX, double dfCenterY,
1649 : double dfRotationDeg, // Ellipse rotation.
1650 : double dfSemiMajor, double dfSemiMinor)
1651 : {
1652 : // Invert the following equation where cosA, sinA are unknown:
1653 : // dfPointOnArcX-dfCenterX = cosA*M*cosRot + sinA*m*sinRot
1654 : // dfPointOnArcY-dfCenterY = -cosA*M*sinRot + sinA*m*cosRot
1655 :
1656 26 : if (dfSemiMajor == 0.0 || dfSemiMinor == 0.0)
1657 0 : return 0.0;
1658 26 : const double dfRotationRadians = dfRotationDeg * M_PI / 180.0;
1659 26 : const double dfCosRot = cos(dfRotationRadians);
1660 26 : const double dfSinRot = sin(dfRotationRadians);
1661 26 : const double dfDeltaX = dfPointOnArcX - dfCenterX;
1662 26 : const double dfDeltaY = dfPointOnArcY - dfCenterY;
1663 26 : const double dfCosA =
1664 26 : (dfCosRot * dfDeltaX - dfSinRot * dfDeltaY) / dfSemiMajor;
1665 26 : const double dfSinA =
1666 26 : (dfSinRot * dfDeltaX + dfCosRot * dfDeltaY) / dfSemiMinor;
1667 : // We could check that dfCosA^2 + dfSinA^2 ~= 1 to verify that the point
1668 : // is on the ellipse.
1669 26 : const double dfAngle = atan2(dfSinA, dfCosA) / M_PI * 180;
1670 26 : if (dfAngle < -180)
1671 0 : return dfAngle + 360;
1672 26 : return dfAngle;
1673 : }
1674 :
1675 : /************************************************************************/
1676 : /* OGRShapeCreateCompoundCurve() */
1677 : /************************************************************************/
1678 :
1679 122 : static OGRCurve *OGRShapeCreateCompoundCurve(int nPartStartIdx, int nPartPoints,
1680 : const CurveSegment *pasCurves,
1681 : int nCurves, int nFirstCurveIdx,
1682 : /* const */ double *padfX,
1683 : /* const */ double *padfY,
1684 : /* const */ double *padfZ,
1685 : /* const */ double *padfM,
1686 : int *pnLastCurveIdx)
1687 : {
1688 244 : auto poCC = std::make_unique<OGRCompoundCurve>();
1689 122 : int nLastPointIdx = nPartStartIdx;
1690 122 : bool bHasCircularArcs = false;
1691 122 : int i = nFirstCurveIdx; // Used after for.
1692 296 : for (; i < nCurves; i++)
1693 : {
1694 200 : const int nStartPointIdx = pasCurves[i].nStartPointIdx;
1695 :
1696 200 : if (nStartPointIdx < nPartStartIdx)
1697 : {
1698 : // Shouldn't happen normally, but who knows.
1699 0 : continue;
1700 : }
1701 :
1702 : // For a multi-part geometry, stop when the start index of the curve
1703 : // exceeds the last point index of the part
1704 200 : if (nStartPointIdx >= nPartStartIdx + nPartPoints)
1705 : {
1706 26 : if (pnLastCurveIdx)
1707 26 : *pnLastCurveIdx = i;
1708 26 : break;
1709 : }
1710 :
1711 : // Add linear segments between end of last curve portion (or beginning
1712 : // of the part) and start of current curve.
1713 174 : if (nStartPointIdx > nLastPointIdx)
1714 : {
1715 39 : OGRLineString *poLine = new OGRLineString();
1716 43 : poLine->setPoints(
1717 39 : nStartPointIdx - nLastPointIdx + 1, padfX + nLastPointIdx,
1718 39 : padfY + nLastPointIdx,
1719 8 : padfZ != nullptr ? padfZ + nLastPointIdx : nullptr,
1720 4 : padfM != nullptr ? padfM + nLastPointIdx : nullptr);
1721 39 : poCC->addCurveDirectly(poLine);
1722 : }
1723 :
1724 174 : if (pasCurves[i].eType == CURVE_ARC_INTERIOR_POINT &&
1725 113 : nStartPointIdx + 1 < nPartStartIdx + nPartPoints)
1726 : {
1727 113 : OGRPoint p1(padfX[nStartPointIdx], padfY[nStartPointIdx],
1728 26 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0,
1729 113 : padfM != nullptr ? padfM[nStartPointIdx] : 0.0);
1730 113 : OGRPoint p2(pasCurves[i].u.ArcByIntermediatePoint.dfX,
1731 113 : pasCurves[i].u.ArcByIntermediatePoint.dfY,
1732 113 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0);
1733 113 : OGRPoint p3(padfX[nStartPointIdx + 1], padfY[nStartPointIdx + 1],
1734 26 : padfZ != nullptr ? padfZ[nStartPointIdx + 1] : 0.0,
1735 113 : padfM != nullptr ? padfM[nStartPointIdx + 1] : 0.0);
1736 :
1737 : // Some software (like QGIS, see https://hub.qgis.org/issues/15116)
1738 : // do not like 3-point circles, so use a 5 point variant.
1739 113 : if (p1.getX() == p3.getX() && p1.getY() == p3.getY())
1740 : {
1741 8 : if (p1.getX() != p2.getX() || p1.getY() != p2.getY())
1742 : {
1743 8 : bHasCircularArcs = true;
1744 8 : OGRCircularString *poCS = new OGRCircularString();
1745 8 : const double dfCenterX = (p1.getX() + p2.getX()) / 2;
1746 8 : const double dfCenterY = (p1.getY() + p2.getY()) / 2;
1747 8 : poCS->addPoint(&p1);
1748 16 : OGRPoint pInterm1(dfCenterX - (p2.getY() - dfCenterY),
1749 8 : dfCenterY + (p1.getX() - dfCenterX),
1750 2 : padfZ != nullptr ? padfZ[nStartPointIdx]
1751 24 : : 0.0);
1752 8 : poCS->addPoint(&pInterm1);
1753 8 : poCS->addPoint(&p2);
1754 16 : OGRPoint pInterm2(dfCenterX + (p2.getY() - dfCenterY),
1755 8 : dfCenterY - (p1.getX() - dfCenterX),
1756 2 : padfZ != nullptr ? padfZ[nStartPointIdx]
1757 24 : : 0.0);
1758 8 : poCS->addPoint(&pInterm2);
1759 8 : poCS->addPoint(&p3);
1760 8 : poCS->set3D(padfZ != nullptr);
1761 8 : poCS->setMeasured(padfM != nullptr);
1762 8 : poCC->addCurveDirectly(poCS);
1763 : }
1764 : }
1765 : else
1766 : {
1767 105 : bHasCircularArcs = true;
1768 105 : OGRCircularString *poCS = new OGRCircularString();
1769 105 : poCS->addPoint(&p1);
1770 105 : poCS->addPoint(&p2);
1771 105 : poCS->addPoint(&p3);
1772 105 : poCS->set3D(padfZ != nullptr);
1773 105 : poCS->setMeasured(padfM != nullptr);
1774 105 : if (poCC->addCurveDirectly(poCS) != OGRERR_NONE)
1775 : {
1776 0 : delete poCS;
1777 0 : return nullptr;
1778 : }
1779 113 : }
1780 : }
1781 :
1782 61 : else if (pasCurves[i].eType == CURVE_ARC_CENTER_POINT &&
1783 12 : nStartPointIdx + 1 < nPartStartIdx + nPartPoints)
1784 : {
1785 12 : const double dfCenterX = pasCurves[i].u.ArcByCenterPoint.dfX;
1786 12 : const double dfCenterY = pasCurves[i].u.ArcByCenterPoint.dfY;
1787 12 : double dfDeltaY = padfY[nStartPointIdx] - dfCenterY;
1788 12 : double dfDeltaX = padfX[nStartPointIdx] - dfCenterX;
1789 12 : double dfAngleStart = atan2(dfDeltaY, dfDeltaX);
1790 12 : dfDeltaY = padfY[nStartPointIdx + 1] - dfCenterY;
1791 12 : dfDeltaX = padfX[nStartPointIdx + 1] - dfCenterX;
1792 12 : double dfAngleEnd = atan2(dfDeltaY, dfDeltaX);
1793 : // Note: This definition from center and 2 points may be
1794 : // not a circle.
1795 12 : double dfRadius = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
1796 12 : if (pasCurves[i].u.ArcByCenterPoint.bIsCCW)
1797 : {
1798 4 : if (dfAngleStart >= dfAngleEnd)
1799 2 : dfAngleEnd += 2 * M_PI;
1800 : }
1801 : else
1802 : {
1803 8 : if (dfAngleStart <= dfAngleEnd)
1804 6 : dfAngleEnd -= 2 * M_PI;
1805 : }
1806 12 : const double dfMidAngle = (dfAngleStart + dfAngleEnd) / 2;
1807 12 : OGRPoint p1(padfX[nStartPointIdx], padfY[nStartPointIdx],
1808 0 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0,
1809 24 : padfM != nullptr ? padfM[nStartPointIdx] : 0.0);
1810 12 : OGRPoint p2(dfCenterX + dfRadius * cos(dfMidAngle),
1811 12 : dfCenterY + dfRadius * sin(dfMidAngle),
1812 24 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0);
1813 12 : OGRPoint p3(padfX[nStartPointIdx + 1], padfY[nStartPointIdx + 1],
1814 0 : padfZ != nullptr ? padfZ[nStartPointIdx + 1] : 0.0,
1815 24 : padfM != nullptr ? padfM[nStartPointIdx + 1] : 0.0);
1816 :
1817 12 : bHasCircularArcs = true;
1818 12 : OGRCircularString *poCS = new OGRCircularString();
1819 12 : poCS->addPoint(&p1);
1820 12 : poCS->addPoint(&p2);
1821 12 : poCS->addPoint(&p3);
1822 12 : poCS->set3D(padfZ != nullptr);
1823 12 : poCS->setMeasured(padfM != nullptr);
1824 24 : poCC->addCurveDirectly(poCS);
1825 : }
1826 :
1827 49 : else if (pasCurves[i].eType == CURVE_BEZIER &&
1828 36 : nStartPointIdx + 1 < nPartStartIdx + nPartPoints)
1829 : {
1830 36 : OGRLineString *poLine = new OGRLineString();
1831 36 : const double dfX0 = padfX[nStartPointIdx];
1832 36 : const double dfY0 = padfY[nStartPointIdx];
1833 36 : const double dfX1 = pasCurves[i].u.Bezier.dfX1;
1834 36 : const double dfY1 = pasCurves[i].u.Bezier.dfY1;
1835 36 : const double dfX2 = pasCurves[i].u.Bezier.dfX2;
1836 36 : const double dfY2 = pasCurves[i].u.Bezier.dfY2;
1837 36 : const double dfX3 = padfX[nStartPointIdx + 1];
1838 36 : const double dfY3 = padfY[nStartPointIdx + 1];
1839 36 : double dfStartAngle = atan2(dfY1 - dfY0, dfX1 - dfX0);
1840 36 : double dfEndAngle = atan2(dfY3 - dfY2, dfX3 - dfX2);
1841 36 : if (dfStartAngle + M_PI < dfEndAngle)
1842 2 : dfStartAngle += 2 * M_PI;
1843 34 : else if (dfEndAngle + M_PI < dfStartAngle)
1844 0 : dfEndAngle += 2 * M_PI;
1845 : const double dfStepSizeRad =
1846 36 : OGRGeometryFactory::GetDefaultArcStepSize() / 180.0 * M_PI;
1847 36 : const double dfLengthTangentStart =
1848 36 : (dfX1 - dfX0) * (dfX1 - dfX0) + (dfY1 - dfY0) * (dfY1 - dfY0);
1849 36 : const double dfLengthTangentEnd =
1850 36 : (dfX3 - dfX2) * (dfX3 - dfX2) + (dfY3 - dfY2) * (dfY3 - dfY2);
1851 36 : const double dfLength =
1852 36 : (dfX3 - dfX0) * (dfX3 - dfX0) + (dfY3 - dfY0) * (dfY3 - dfY0);
1853 : // Heuristics to compute number of steps: we take into account the
1854 : // angular difference between the start and end tangent. And we
1855 : // also take into account the relative length of the tangent vs
1856 : // the length of the straight segment
1857 : const int nSteps =
1858 : (dfLength < 1e-9)
1859 36 : ? 1
1860 36 : : static_cast<int>(std::min(
1861 72 : 1000.0,
1862 144 : ceil(std::max(2.0, fabs(dfEndAngle - dfStartAngle) /
1863 36 : dfStepSizeRad) *
1864 108 : std::max(1.0, 5.0 *
1865 36 : (dfLengthTangentStart +
1866 72 : dfLengthTangentEnd) /
1867 36 : dfLength))));
1868 36 : poLine->setNumPoints(nSteps + 1);
1869 84 : poLine->setPoint(0, dfX0, dfY0,
1870 24 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0,
1871 24 : padfM != nullptr ? padfM[nStartPointIdx] : 0.0);
1872 1317 : for (int j = 1; j < nSteps; j++)
1873 : {
1874 1281 : const double t = static_cast<double>(j) / nSteps;
1875 : // Third-order Bezier interpolation.
1876 1281 : poLine->setPoint(
1877 : j,
1878 1281 : (1 - t) * (1 - t) * (1 - t) * dfX0 +
1879 1281 : 3 * (1 - t) * (1 - t) * t * dfX1 +
1880 1281 : 3 * (1 - t) * t * t * dfX2 + t * t * t * dfX3,
1881 1281 : (1 - t) * (1 - t) * (1 - t) * dfY0 +
1882 1281 : 3 * (1 - t) * (1 - t) * t * dfY1 +
1883 1281 : 3 * (1 - t) * t * t * dfY2 + t * t * t * dfY3);
1884 : }
1885 84 : poLine->setPoint(nSteps, dfX3, dfY3,
1886 24 : padfZ != nullptr ? padfZ[nStartPointIdx + 1] : 0.0,
1887 24 : padfM != nullptr ? padfM[nStartPointIdx + 1]
1888 : : 0.0);
1889 36 : poLine->set3D(padfZ != nullptr);
1890 36 : poLine->setMeasured(padfM != nullptr);
1891 36 : if (poCC->addCurveDirectly(poLine) != OGRERR_NONE)
1892 : {
1893 0 : delete poLine;
1894 0 : return nullptr;
1895 36 : }
1896 : }
1897 :
1898 13 : else if (pasCurves[i].eType == CURVE_ELLIPSE_BY_CENTER &&
1899 13 : nStartPointIdx + 1 < nPartStartIdx + nPartPoints)
1900 : {
1901 13 : const double dfSemiMinor =
1902 13 : pasCurves[i].u.EllipseByCenter.dfSemiMajor *
1903 13 : pasCurves[i].u.EllipseByCenter.dfRatioSemiMinor;
1904 : // Different sign conventions between ext shape
1905 : // (trigonometric, CCW) and approximateArcAngles (CW).
1906 13 : const double dfRotationDeg =
1907 13 : -pasCurves[i].u.EllipseByCenter.dfRotationDeg;
1908 13 : const double dfAngleStart = GetAngleOnEllipse(
1909 13 : padfX[nStartPointIdx], padfY[nStartPointIdx],
1910 13 : pasCurves[i].u.EllipseByCenter.dfX,
1911 13 : pasCurves[i].u.EllipseByCenter.dfY, dfRotationDeg,
1912 13 : pasCurves[i].u.EllipseByCenter.dfSemiMajor, dfSemiMinor);
1913 13 : const double dfAngleEnd = GetAngleOnEllipse(
1914 13 : padfX[nStartPointIdx + 1], padfY[nStartPointIdx + 1],
1915 13 : pasCurves[i].u.EllipseByCenter.dfX,
1916 13 : pasCurves[i].u.EllipseByCenter.dfY, dfRotationDeg,
1917 13 : pasCurves[i].u.EllipseByCenter.dfSemiMajor, dfSemiMinor);
1918 : // CPLDebug("OGR", "Start angle=%f, End angle=%f",
1919 : // dfAngleStart, dfAngleEnd);
1920 : // Approximatearcangles() use CW.
1921 13 : double dfAngleStartForApprox = -dfAngleStart;
1922 13 : double dfAngleEndForApprox = -dfAngleEnd;
1923 13 : if (pasCurves[i].u.EllipseByCenter.bIsComplete)
1924 : {
1925 9 : dfAngleEndForApprox = dfAngleStartForApprox + 360;
1926 : }
1927 4 : else if (pasCurves[i].u.EllipseByCenter.bIsMinor)
1928 : {
1929 4 : if (dfAngleEndForApprox > dfAngleStartForApprox + 180)
1930 : {
1931 0 : dfAngleEndForApprox -= 360;
1932 : }
1933 4 : else if (dfAngleEndForApprox < dfAngleStartForApprox - 180)
1934 : {
1935 0 : dfAngleEndForApprox += 360;
1936 : }
1937 : }
1938 : else
1939 : {
1940 0 : if (dfAngleEndForApprox > dfAngleStartForApprox &&
1941 0 : dfAngleEndForApprox < dfAngleStartForApprox + 180)
1942 : {
1943 0 : dfAngleEndForApprox -= 360;
1944 : }
1945 0 : else if (dfAngleEndForApprox < dfAngleStartForApprox &&
1946 0 : dfAngleEndForApprox > dfAngleStartForApprox - 180)
1947 : {
1948 0 : dfAngleEndForApprox += 360;
1949 : }
1950 : }
1951 : OGRLineString *poLine =
1952 : OGRGeometryFactory::approximateArcAngles(
1953 13 : pasCurves[i].u.EllipseByCenter.dfX,
1954 13 : pasCurves[i].u.EllipseByCenter.dfY,
1955 2 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0,
1956 13 : pasCurves[i].u.EllipseByCenter.dfSemiMajor, dfSemiMinor,
1957 : dfRotationDeg, dfAngleStartForApprox, dfAngleEndForApprox,
1958 : 0)
1959 13 : ->toLineString();
1960 13 : if (poLine->getNumPoints() >= 2)
1961 : {
1962 15 : poLine->setPoint(
1963 13 : 0, padfX[nStartPointIdx], padfY[nStartPointIdx],
1964 2 : padfZ != nullptr ? padfZ[nStartPointIdx] : 0.0,
1965 2 : padfM != nullptr ? padfM[nStartPointIdx] : 0.0);
1966 15 : poLine->setPoint(
1967 13 : poLine->getNumPoints() - 1, padfX[nStartPointIdx + 1],
1968 13 : padfY[nStartPointIdx + 1],
1969 2 : padfZ != nullptr ? padfZ[nStartPointIdx + 1] : 0.0,
1970 2 : padfM != nullptr ? padfM[nStartPointIdx + 1] : 0.0);
1971 : }
1972 13 : poLine->set3D(padfZ != nullptr);
1973 13 : poLine->setMeasured(padfM != nullptr);
1974 13 : poCC->addCurveDirectly(poLine);
1975 : }
1976 :
1977 : // Should not happen normally.
1978 0 : else if (nStartPointIdx + 1 < nPartStartIdx + nPartPoints)
1979 : {
1980 0 : OGRLineString *poLine = new OGRLineString();
1981 0 : poLine->setPoints(
1982 0 : 2, padfX + nStartPointIdx, padfY + nStartPointIdx,
1983 0 : padfZ != nullptr ? padfZ + nStartPointIdx : nullptr,
1984 0 : padfM != nullptr ? padfM + nStartPointIdx : nullptr);
1985 0 : poCC->addCurveDirectly(poLine);
1986 : }
1987 174 : nLastPointIdx = nStartPointIdx + 1;
1988 : }
1989 122 : if (i == nCurves)
1990 : {
1991 96 : if (pnLastCurveIdx)
1992 17 : *pnLastCurveIdx = i;
1993 : }
1994 :
1995 : // Add terminating linear segments
1996 122 : if (nLastPointIdx < nPartStartIdx + nPartPoints - 1)
1997 : {
1998 49 : OGRLineString *poLine = new OGRLineString();
1999 59 : poLine->setPoints(nPartStartIdx + nPartPoints - 1 - nLastPointIdx + 1,
2000 49 : padfX + nLastPointIdx, padfY + nLastPointIdx,
2001 14 : padfZ != nullptr ? padfZ + nLastPointIdx : nullptr,
2002 10 : padfM != nullptr ? padfM + nLastPointIdx : nullptr);
2003 49 : if (poCC->addCurveDirectly(poLine) != OGRERR_NONE)
2004 : {
2005 0 : delete poLine;
2006 0 : return nullptr;
2007 : }
2008 : }
2009 :
2010 122 : if (!bHasCircularArcs)
2011 : {
2012 29 : OGRwkbGeometryType eLSType = wkbLineString;
2013 29 : if (OGR_GT_HasZ(poCC->getGeometryType()))
2014 7 : eLSType = OGR_GT_SetZ(eLSType);
2015 29 : if (OGR_GT_HasM(poCC->getGeometryType()))
2016 7 : eLSType = OGR_GT_SetM(eLSType);
2017 29 : return reinterpret_cast<OGRCurve *>(OGR_G_ForceTo(
2018 58 : OGRGeometry::ToHandle(poCC.release()), eLSType, nullptr));
2019 : }
2020 : else
2021 93 : return poCC.release();
2022 : }
2023 :
2024 : /************************************************************************/
2025 : /* OGRCreateFromShapeBin() */
2026 : /* */
2027 : /* Translate shapefile binary representation to an OGR */
2028 : /* geometry. */
2029 : /************************************************************************/
2030 :
2031 1767 : OGRErr OGRCreateFromShapeBin(GByte *pabyShape, OGRGeometry **ppoGeom,
2032 : int nBytes)
2033 :
2034 : {
2035 1767 : *ppoGeom = nullptr;
2036 :
2037 1767 : if (nBytes < 4)
2038 : {
2039 0 : CPLError(CE_Failure, CPLE_AppDefined,
2040 : "Shape buffer size (%d) too small", nBytes);
2041 0 : return OGRERR_FAILURE;
2042 : }
2043 :
2044 : /* -------------------------------------------------------------------- */
2045 : /* Detect zlib compressed shapes and uncompress buffer if necessary */
2046 : /* NOTE: this seems to be an undocumented feature, even in the */
2047 : /* extended_shapefile_format.pdf found in the FileGDB API */
2048 : /* documentation. */
2049 : /* -------------------------------------------------------------------- */
2050 1767 : if (nBytes >= 14 && pabyShape[12] == 0x78 &&
2051 0 : pabyShape[13] == 0xDA /* zlib marker */)
2052 : {
2053 0 : GInt32 nUncompressedSize = 0;
2054 0 : GInt32 nCompressedSize = 0;
2055 0 : memcpy(&nUncompressedSize, pabyShape + 4, 4);
2056 0 : memcpy(&nCompressedSize, pabyShape + 8, 4);
2057 0 : CPL_LSBPTR32(&nUncompressedSize);
2058 0 : CPL_LSBPTR32(&nCompressedSize);
2059 0 : if (nCompressedSize + 12 == nBytes && nUncompressedSize > 0)
2060 : {
2061 : GByte *pabyUncompressedBuffer =
2062 0 : static_cast<GByte *>(VSI_MALLOC_VERBOSE(nUncompressedSize));
2063 0 : if (pabyUncompressedBuffer == nullptr)
2064 : {
2065 0 : return OGRERR_FAILURE;
2066 : }
2067 :
2068 0 : size_t nRealUncompressedSize = 0;
2069 0 : if (CPLZLibInflate(pabyShape + 12, nCompressedSize,
2070 : pabyUncompressedBuffer, nUncompressedSize,
2071 0 : &nRealUncompressedSize) == nullptr)
2072 : {
2073 0 : CPLError(CE_Failure, CPLE_AppDefined,
2074 : "CPLZLibInflate() failed");
2075 0 : VSIFree(pabyUncompressedBuffer);
2076 0 : return OGRERR_FAILURE;
2077 : }
2078 :
2079 : const OGRErr eErr =
2080 0 : OGRCreateFromShapeBin(pabyUncompressedBuffer, ppoGeom,
2081 : static_cast<int>(nRealUncompressedSize));
2082 :
2083 0 : VSIFree(pabyUncompressedBuffer);
2084 :
2085 0 : return eErr;
2086 : }
2087 : }
2088 :
2089 1767 : int nSHPType = pabyShape[0];
2090 :
2091 : /* -------------------------------------------------------------------- */
2092 : /* Return a NULL geometry when SHPT_NULL is encountered. */
2093 : /* Watch out, null return does not mean "bad data" it means */
2094 : /* "no geometry here". Watch the OGRErr for the error status */
2095 : /* -------------------------------------------------------------------- */
2096 1767 : if (nSHPType == SHPT_NULL)
2097 : {
2098 0 : *ppoGeom = nullptr;
2099 0 : return OGRERR_NONE;
2100 : }
2101 :
2102 : #if DEBUG_VERBOSE
2103 : CPLDebug("PGeo", "Shape type read from PGeo data is nSHPType = %d",
2104 : nSHPType);
2105 : #endif
2106 :
2107 1767 : const bool bIsExtended =
2108 1767 : nSHPType >= SHPT_GENERALPOLYLINE && nSHPType <= SHPT_GENERALMULTIPATCH;
2109 :
2110 1767 : const bool bHasZ =
2111 1649 : (nSHPType == SHPT_POINTZ || nSHPType == SHPT_POINTZM ||
2112 1531 : nSHPType == SHPT_MULTIPOINTZ || nSHPType == SHPT_MULTIPOINTZM ||
2113 1295 : nSHPType == SHPT_POLYGONZ || nSHPType == SHPT_POLYGONZM ||
2114 1055 : nSHPType == SHPT_ARCZ || nSHPType == SHPT_ARCZM ||
2115 3513 : nSHPType == SHPT_MULTIPATCH || nSHPType == SHPT_MULTIPATCHM ||
2116 97 : (bIsExtended && (pabyShape[3] & 0x80) != 0));
2117 :
2118 1767 : const bool bHasM =
2119 1767 : (nSHPType == SHPT_POINTM || nSHPType == SHPT_POINTZM ||
2120 1767 : nSHPType == SHPT_MULTIPOINTM || nSHPType == SHPT_MULTIPOINTZM ||
2121 1767 : nSHPType == SHPT_POLYGONM || nSHPType == SHPT_POLYGONZM ||
2122 1763 : nSHPType == SHPT_ARCM || nSHPType == SHPT_ARCZM ||
2123 3631 : nSHPType == SHPT_MULTIPATCHM ||
2124 97 : (bIsExtended && (pabyShape[3] & 0x40) != 0));
2125 :
2126 1767 : const bool bHasCurves = (bIsExtended && (pabyShape[3] & 0x20) != 0);
2127 :
2128 1767 : switch (nSHPType)
2129 : {
2130 44 : case SHPT_GENERALPOLYLINE:
2131 44 : nSHPType = SHPT_ARC;
2132 44 : break;
2133 53 : case SHPT_GENERALPOLYGON:
2134 53 : nSHPType = SHPT_POLYGON;
2135 53 : break;
2136 0 : case SHPT_GENERALPOINT:
2137 0 : nSHPType = SHPT_POINT;
2138 0 : break;
2139 0 : case SHPT_GENERALMULTIPOINT:
2140 0 : nSHPType = SHPT_MULTIPOINT;
2141 0 : break;
2142 0 : case SHPT_GENERALMULTIPATCH:
2143 0 : nSHPType = SHPT_MULTIPATCH;
2144 : }
2145 :
2146 : /* ==================================================================== */
2147 : /* Extract vertices for a Polygon or Arc. */
2148 : /* ==================================================================== */
2149 1767 : if (nSHPType == SHPT_ARC || nSHPType == SHPT_ARCZ ||
2150 1248 : nSHPType == SHPT_ARCM || nSHPType == SHPT_ARCZM ||
2151 954 : nSHPType == SHPT_POLYGON || nSHPType == SHPT_POLYGONZ ||
2152 718 : nSHPType == SHPT_POLYGONM || nSHPType == SHPT_POLYGONZM ||
2153 478 : nSHPType == SHPT_MULTIPATCH || nSHPType == SHPT_MULTIPATCHM)
2154 : {
2155 1289 : if (nBytes < 44)
2156 : {
2157 0 : CPLError(CE_Failure, CPLE_AppDefined,
2158 : "Corrupted Shape : nBytes=%d, nSHPType=%d", nBytes,
2159 : nSHPType);
2160 0 : return OGRERR_FAILURE;
2161 : }
2162 :
2163 : /* --------------------------------------------------------------------
2164 : */
2165 : /* Extract part/point count, and build vertex and part arrays */
2166 : /* to proper size. */
2167 : /* --------------------------------------------------------------------
2168 : */
2169 1289 : GInt32 nPoints = 0;
2170 1289 : memcpy(&nPoints, pabyShape + 40, 4);
2171 1289 : GInt32 nParts = 0;
2172 1289 : memcpy(&nParts, pabyShape + 36, 4);
2173 :
2174 1289 : CPL_LSBPTR32(&nPoints);
2175 1289 : CPL_LSBPTR32(&nParts);
2176 :
2177 1289 : if (nPoints < 0 || nParts < 0 || nPoints > 50 * 1000 * 1000 ||
2178 1289 : nParts > 10 * 1000 * 1000)
2179 : {
2180 0 : CPLError(CE_Failure, CPLE_AppDefined,
2181 : "Corrupted Shape : nPoints=%d, nParts=%d.", nPoints,
2182 : nParts);
2183 0 : return OGRERR_FAILURE;
2184 : }
2185 :
2186 1289 : const bool bIsMultiPatch =
2187 1289 : nSHPType == SHPT_MULTIPATCH || nSHPType == SHPT_MULTIPATCHM;
2188 :
2189 : // With the previous checks on nPoints and nParts,
2190 : // we should not overflow here and after
2191 : // since 50 M * (16 + 8 + 8) = 1 600 MB.
2192 1289 : int nRequiredSize = 44 + 4 * nParts + 16 * nPoints;
2193 1289 : if (bHasZ)
2194 : {
2195 735 : nRequiredSize += 16 + 8 * nPoints;
2196 : }
2197 1289 : if (bHasM)
2198 : {
2199 25 : nRequiredSize += 16 + 8 * nPoints;
2200 : }
2201 1289 : if (bIsMultiPatch)
2202 : {
2203 236 : nRequiredSize += 4 * nParts;
2204 : }
2205 1289 : if (nRequiredSize > nBytes)
2206 : {
2207 0 : CPLError(CE_Failure, CPLE_AppDefined,
2208 : "Corrupted Shape : nPoints=%d, nParts=%d, nBytes=%d, "
2209 : "nSHPType=%d, nRequiredSize=%d",
2210 : nPoints, nParts, nBytes, nSHPType, nRequiredSize);
2211 0 : return OGRERR_FAILURE;
2212 : }
2213 :
2214 : GInt32 *panPartStart =
2215 1289 : static_cast<GInt32 *>(VSI_MALLOC2_VERBOSE(nParts, sizeof(GInt32)));
2216 1289 : if (nParts != 0 && panPartStart == nullptr)
2217 : {
2218 0 : return OGRERR_FAILURE;
2219 : }
2220 :
2221 : /* --------------------------------------------------------------------
2222 : */
2223 : /* Copy out the part array from the record. */
2224 : /* --------------------------------------------------------------------
2225 : */
2226 1289 : memcpy(panPartStart, pabyShape + 44, 4 * nParts);
2227 3446 : for (int i = 0; i < nParts; i++)
2228 : {
2229 2157 : CPL_LSBPTR32(panPartStart + i);
2230 :
2231 : // Check that the offset is inside the vertex array.
2232 2157 : if (panPartStart[i] < 0 || panPartStart[i] >= nPoints)
2233 : {
2234 0 : CPLError(
2235 : CE_Failure, CPLE_AppDefined,
2236 : "Corrupted Shape : panPartStart[%d] = %d, nPoints = %d", i,
2237 0 : panPartStart[i], nPoints);
2238 0 : CPLFree(panPartStart);
2239 0 : return OGRERR_FAILURE;
2240 : }
2241 2157 : if (i > 0 && panPartStart[i] <= panPartStart[i - 1])
2242 : {
2243 0 : CPLError(CE_Failure, CPLE_AppDefined,
2244 : "Corrupted Shape : panPartStart[%d] = %d, "
2245 : "panPartStart[%d] = %d",
2246 0 : i, panPartStart[i], i - 1, panPartStart[i - 1]);
2247 0 : CPLFree(panPartStart);
2248 0 : return OGRERR_FAILURE;
2249 : }
2250 : }
2251 :
2252 1289 : int nOffset = 44 + 4 * nParts;
2253 :
2254 : /* --------------------------------------------------------------------
2255 : */
2256 : /* If this is a multipatch, we will also have parts types. */
2257 : /* --------------------------------------------------------------------
2258 : */
2259 1289 : GInt32 *panPartType = nullptr;
2260 :
2261 1289 : if (bIsMultiPatch)
2262 : {
2263 : panPartType = static_cast<GInt32 *>(
2264 236 : VSI_MALLOC2_VERBOSE(nParts, sizeof(GInt32)));
2265 236 : if (panPartType == nullptr)
2266 : {
2267 0 : CPLFree(panPartStart);
2268 0 : return OGRERR_FAILURE;
2269 : }
2270 :
2271 236 : memcpy(panPartType, pabyShape + nOffset, 4 * nParts);
2272 944 : for (int i = 0; i < nParts; i++)
2273 : {
2274 708 : CPL_LSBPTR32(panPartType + i);
2275 : }
2276 236 : nOffset += 4 * nParts;
2277 : }
2278 :
2279 : /* --------------------------------------------------------------------
2280 : */
2281 : /* Copy out the vertices from the record. */
2282 : /* --------------------------------------------------------------------
2283 : */
2284 : double *padfX =
2285 1289 : static_cast<double *>(VSI_MALLOC_VERBOSE(sizeof(double) * nPoints));
2286 : double *padfY =
2287 1289 : static_cast<double *>(VSI_MALLOC_VERBOSE(sizeof(double) * nPoints));
2288 : double *padfZ =
2289 1289 : static_cast<double *>(VSI_CALLOC_VERBOSE(sizeof(double), nPoints));
2290 1289 : double *padfM = static_cast<double *>(
2291 25 : bHasM ? VSI_CALLOC_VERBOSE(sizeof(double), nPoints) : nullptr);
2292 :
2293 1289 : if (nPoints != 0 && (padfX == nullptr || padfY == nullptr ||
2294 1288 : padfZ == nullptr || (bHasM && padfM == nullptr)))
2295 : {
2296 0 : CPLFree(panPartStart);
2297 0 : CPLFree(panPartType);
2298 0 : CPLFree(padfX);
2299 0 : CPLFree(padfY);
2300 0 : CPLFree(padfZ);
2301 0 : CPLFree(padfM);
2302 0 : return OGRERR_FAILURE;
2303 : }
2304 :
2305 10974 : for (int i = 0; i < nPoints; i++)
2306 : {
2307 9685 : memcpy(padfX + i, pabyShape + nOffset + i * 16, 8);
2308 9685 : memcpy(padfY + i, pabyShape + nOffset + i * 16 + 8, 8);
2309 9685 : CPL_LSBPTR64(padfX + i);
2310 9685 : CPL_LSBPTR64(padfY + i);
2311 : }
2312 :
2313 1289 : nOffset += 16 * nPoints;
2314 :
2315 : /* --------------------------------------------------------------------
2316 : */
2317 : /* If we have a Z coordinate, collect that now. */
2318 : /* --------------------------------------------------------------------
2319 : */
2320 1289 : if (bHasZ)
2321 : {
2322 5477 : for (int i = 0; i < nPoints; i++)
2323 : {
2324 4742 : memcpy(padfZ + i, pabyShape + nOffset + 16 + i * 8, 8);
2325 4742 : CPL_LSBPTR64(padfZ + i);
2326 : }
2327 :
2328 735 : nOffset += 16 + 8 * nPoints;
2329 : }
2330 :
2331 : /* --------------------------------------------------------------------
2332 : */
2333 : /* If we have a M coordinate, collect that now. */
2334 : /* --------------------------------------------------------------------
2335 : */
2336 1289 : if (bHasM)
2337 : {
2338 25 : bool bIsAllNAN = nPoints > 0;
2339 155 : for (int i = 0; i < nPoints; i++)
2340 : {
2341 130 : memcpy(padfM + i, pabyShape + nOffset + 16 + i * 8, 8);
2342 130 : CPL_LSBPTR64(padfM + i);
2343 130 : bIsAllNAN &= std::isnan(padfM[i]);
2344 : }
2345 25 : if (bIsAllNAN)
2346 : {
2347 4 : CPLFree(padfM);
2348 4 : padfM = nullptr;
2349 : }
2350 :
2351 25 : nOffset += 16 + 8 * nPoints;
2352 : }
2353 :
2354 : /* --------------------------------------------------------------------
2355 : */
2356 : /* If we have curves, collect them now. */
2357 : /* --------------------------------------------------------------------
2358 : */
2359 1289 : int nCurves = 0;
2360 1289 : CurveSegment *pasCurves = nullptr;
2361 1289 : if (bHasCurves && nOffset + 4 <= nBytes)
2362 : {
2363 97 : memcpy(&nCurves, pabyShape + nOffset, 4);
2364 97 : CPL_LSBPTR32(&nCurves);
2365 97 : nOffset += 4;
2366 : #ifdef DEBUG_VERBOSE
2367 : CPLDebug("OGR", "nCurves = %d", nCurves);
2368 : #endif
2369 97 : if (nCurves < 0 || nCurves > (nBytes - nOffset) / (8 + 20))
2370 : {
2371 0 : CPLDebug("OGR", "Invalid nCurves = %d", nCurves);
2372 0 : nCurves = 0;
2373 : }
2374 : pasCurves = static_cast<CurveSegment *>(
2375 97 : VSI_MALLOC2_VERBOSE(sizeof(CurveSegment), nCurves));
2376 97 : if (pasCurves == nullptr)
2377 : {
2378 0 : nCurves = 0;
2379 : }
2380 97 : int iCurve = 0;
2381 272 : for (int i = 0; i < nCurves; i++)
2382 : {
2383 175 : if (nOffset + 8 > nBytes)
2384 : {
2385 0 : CPLDebug("OGR", "Not enough bytes");
2386 0 : break;
2387 : }
2388 175 : int nStartPointIdx = 0;
2389 175 : memcpy(&nStartPointIdx, pabyShape + nOffset, 4);
2390 175 : CPL_LSBPTR32(&nStartPointIdx);
2391 175 : nOffset += 4;
2392 175 : int nSegmentType = 0;
2393 175 : memcpy(&nSegmentType, pabyShape + nOffset, 4);
2394 175 : CPL_LSBPTR32(&nSegmentType);
2395 175 : nOffset += 4;
2396 : #ifdef DEBUG_VERBOSE
2397 : CPLDebug("OGR", "[%d] nStartPointIdx = %d, segmentType = %d", i,
2398 : nSegmentType, nSegmentType);
2399 : #endif
2400 175 : if (nStartPointIdx < 0 || nStartPointIdx >= nPoints ||
2401 78 : (iCurve > 0 &&
2402 78 : nStartPointIdx <= pasCurves[iCurve - 1].nStartPointIdx))
2403 : {
2404 0 : CPLDebug("OGR", "Invalid nStartPointIdx = %d",
2405 : nStartPointIdx);
2406 0 : break;
2407 : }
2408 175 : pasCurves[iCurve].nStartPointIdx = nStartPointIdx;
2409 175 : if (nSegmentType == EXT_SHAPE_SEGMENT_ARC)
2410 : {
2411 126 : if (nOffset + 20 > nBytes)
2412 : {
2413 0 : CPLDebug("OGR", "Not enough bytes");
2414 0 : break;
2415 : }
2416 126 : double dfVal1 = 0.0;
2417 126 : double dfVal2 = 0.0;
2418 126 : memcpy(&dfVal1, pabyShape + nOffset + 0, 8);
2419 126 : CPL_LSBPTR64(&dfVal1);
2420 126 : memcpy(&dfVal2, pabyShape + nOffset + 8, 8);
2421 126 : CPL_LSBPTR64(&dfVal2);
2422 126 : int nBits = 0;
2423 126 : memcpy(&nBits, pabyShape + nOffset + 16, 4);
2424 126 : CPL_LSBPTR32(&nBits);
2425 :
2426 : #ifdef DEBUG_VERBOSE
2427 : CPLDebug("OGR", "Arc: ");
2428 : CPLDebug("OGR", " dfVal1 = %f, dfVal2 = %f, nBits=%X",
2429 : dfVal1, dfVal2, nBits);
2430 : if (nBits & EXT_SHAPE_ARC_EMPTY)
2431 : CPLDebug("OGR", " IsEmpty");
2432 : if (nBits & EXT_SHAPE_ARC_CCW)
2433 : CPLDebug("OGR", " IsCCW");
2434 : if (nBits & EXT_SHAPE_ARC_MINOR)
2435 : CPLDebug("OGR", " IsMinor");
2436 : if (nBits & EXT_SHAPE_ARC_LINE)
2437 : CPLDebug("OGR", " IsLine");
2438 : if (nBits & EXT_SHAPE_ARC_POINT)
2439 : CPLDebug("OGR", " IsPoint");
2440 : if (nBits & EXT_SHAPE_ARC_IP)
2441 : CPLDebug("OGR", " DefinedIP");
2442 : #endif
2443 126 : if ((nBits & EXT_SHAPE_ARC_IP) != 0 &&
2444 114 : (nBits & EXT_SHAPE_ARC_LINE) == 0)
2445 : {
2446 113 : pasCurves[iCurve].eType = CURVE_ARC_INTERIOR_POINT;
2447 113 : pasCurves[iCurve].u.ArcByIntermediatePoint.dfX = dfVal1;
2448 113 : pasCurves[iCurve].u.ArcByIntermediatePoint.dfY = dfVal2;
2449 113 : iCurve++;
2450 : }
2451 13 : else if ((nBits & EXT_SHAPE_ARC_EMPTY) == 0 &&
2452 13 : (nBits & EXT_SHAPE_ARC_LINE) == 0 &&
2453 12 : (nBits & EXT_SHAPE_ARC_POINT) == 0)
2454 : {
2455 : // This is the old deprecated way
2456 12 : pasCurves[iCurve].eType = CURVE_ARC_CENTER_POINT;
2457 12 : pasCurves[iCurve].u.ArcByCenterPoint.dfX = dfVal1;
2458 12 : pasCurves[iCurve].u.ArcByCenterPoint.dfY = dfVal2;
2459 12 : pasCurves[iCurve].u.ArcByCenterPoint.bIsCCW =
2460 12 : (nBits & EXT_SHAPE_ARC_CCW) != 0;
2461 12 : iCurve++;
2462 : }
2463 126 : nOffset += 16 + 4;
2464 : }
2465 49 : else if (nSegmentType == EXT_SHAPE_SEGMENT_BEZIER)
2466 : {
2467 36 : if (nOffset + 32 > nBytes)
2468 : {
2469 0 : CPLDebug("OGR", "Not enough bytes");
2470 0 : break;
2471 : }
2472 36 : double dfX1 = 0.0;
2473 36 : double dfY1 = 0.0;
2474 36 : memcpy(&dfX1, pabyShape + nOffset + 0, 8);
2475 36 : CPL_LSBPTR64(&dfX1);
2476 36 : memcpy(&dfY1, pabyShape + nOffset + 8, 8);
2477 36 : CPL_LSBPTR64(&dfY1);
2478 36 : double dfX2 = 0.0;
2479 36 : double dfY2 = 0.0;
2480 36 : memcpy(&dfX2, pabyShape + nOffset + 16, 8);
2481 36 : CPL_LSBPTR64(&dfX2);
2482 36 : memcpy(&dfY2, pabyShape + nOffset + 24, 8);
2483 36 : CPL_LSBPTR64(&dfY2);
2484 : #ifdef DEBUG_VERBOSE
2485 : CPLDebug("OGR", "Bezier:");
2486 : CPLDebug("OGR", " dfX1 = %f, dfY1 = %f", dfX1, dfY1);
2487 : CPLDebug("OGR", " dfX2 = %f, dfY2 = %f", dfX2, dfY2);
2488 : #endif
2489 36 : pasCurves[iCurve].eType = CURVE_BEZIER;
2490 36 : pasCurves[iCurve].u.Bezier.dfX1 = dfX1;
2491 36 : pasCurves[iCurve].u.Bezier.dfY1 = dfY1;
2492 36 : pasCurves[iCurve].u.Bezier.dfX2 = dfX2;
2493 36 : pasCurves[iCurve].u.Bezier.dfY2 = dfY2;
2494 36 : iCurve++;
2495 36 : nOffset += 32;
2496 : }
2497 13 : else if (nSegmentType == EXT_SHAPE_SEGMENT_ELLIPSE)
2498 : {
2499 13 : if (nOffset + 44 > nBytes)
2500 : {
2501 0 : CPLDebug("OGR", "Not enough bytes");
2502 0 : break;
2503 : }
2504 13 : double dfVS0 = 0.0;
2505 13 : memcpy(&dfVS0, pabyShape + nOffset, 8);
2506 13 : nOffset += 8;
2507 13 : CPL_LSBPTR64(&dfVS0);
2508 :
2509 13 : double dfVS1 = 0.0;
2510 13 : memcpy(&dfVS1, pabyShape + nOffset, 8);
2511 13 : nOffset += 8;
2512 13 : CPL_LSBPTR64(&dfVS1);
2513 :
2514 13 : double dfRotationOrFromV = 0.0;
2515 13 : memcpy(&dfRotationOrFromV, pabyShape + nOffset, 8);
2516 13 : nOffset += 8;
2517 13 : CPL_LSBPTR64(&dfRotationOrFromV);
2518 :
2519 13 : double dfSemiMajor = 0.0;
2520 13 : memcpy(&dfSemiMajor, pabyShape + nOffset, 8);
2521 13 : nOffset += 8;
2522 13 : CPL_LSBPTR64(&dfSemiMajor);
2523 :
2524 13 : double dfMinorMajorRatioOrDeltaV = 0.0;
2525 13 : memcpy(&dfMinorMajorRatioOrDeltaV, pabyShape + nOffset, 8);
2526 13 : nOffset += 8;
2527 13 : CPL_LSBPTR64(&dfMinorMajorRatioOrDeltaV);
2528 :
2529 13 : int nBits = 0;
2530 13 : memcpy(&nBits, pabyShape + nOffset, 4);
2531 13 : CPL_LSBPTR32(&nBits);
2532 13 : nOffset += 4;
2533 :
2534 : #ifdef DEBUG_VERBOSE
2535 : CPLDebug("OGR", "Ellipse:");
2536 : CPLDebug("OGR", " dfVS0 = %f", dfVS0);
2537 : CPLDebug("OGR", " dfVS1 = %f", dfVS1);
2538 : CPLDebug("OGR", " dfRotationOrFromV = %f",
2539 : dfRotationOrFromV);
2540 : CPLDebug("OGR", " dfSemiMajor = %f", dfSemiMajor);
2541 : CPLDebug("OGR", " dfMinorMajorRatioOrDeltaV = %f",
2542 : dfMinorMajorRatioOrDeltaV);
2543 : CPLDebug("OGR", " nBits=%X", nBits);
2544 :
2545 : if (nBits & EXT_SHAPE_ELLIPSE_EMPTY)
2546 : CPLDebug("OGR", " IsEmpty");
2547 : if (nBits & EXT_SHAPE_ELLIPSE_LINE)
2548 : CPLDebug("OGR", " IsLine");
2549 : if (nBits & EXT_SHAPE_ELLIPSE_POINT)
2550 : CPLDebug("OGR", " IsPoint");
2551 : if (nBits & EXT_SHAPE_ELLIPSE_CIRCULAR)
2552 : CPLDebug("OGR", " IsCircular");
2553 : if (nBits & EXT_SHAPE_ELLIPSE_CENTER_TO)
2554 : CPLDebug("OGR", " CenterTo");
2555 : if (nBits & EXT_SHAPE_ELLIPSE_CENTER_FROM)
2556 : CPLDebug("OGR", " CenterFrom");
2557 : if (nBits & EXT_SHAPE_ELLIPSE_CCW)
2558 : CPLDebug("OGR", " IsCCW");
2559 : if (nBits & EXT_SHAPE_ELLIPSE_MINOR)
2560 : CPLDebug("OGR", " IsMinor");
2561 : if (nBits & EXT_SHAPE_ELLIPSE_COMPLETE)
2562 : CPLDebug("OGR", " IsComplete");
2563 : #endif
2564 :
2565 13 : if ((nBits & EXT_SHAPE_ELLIPSE_CENTER_TO) == 0 &&
2566 13 : (nBits & EXT_SHAPE_ELLIPSE_CENTER_FROM) == 0)
2567 : {
2568 13 : pasCurves[iCurve].eType = CURVE_ELLIPSE_BY_CENTER;
2569 13 : pasCurves[iCurve].u.EllipseByCenter.dfX = dfVS0;
2570 13 : pasCurves[iCurve].u.EllipseByCenter.dfY = dfVS1;
2571 13 : pasCurves[iCurve].u.EllipseByCenter.dfRotationDeg =
2572 13 : dfRotationOrFromV / M_PI * 180;
2573 13 : pasCurves[iCurve].u.EllipseByCenter.dfSemiMajor =
2574 : dfSemiMajor;
2575 13 : pasCurves[iCurve].u.EllipseByCenter.dfRatioSemiMinor =
2576 : dfMinorMajorRatioOrDeltaV;
2577 13 : pasCurves[iCurve].u.EllipseByCenter.bIsMinor =
2578 13 : ((nBits & EXT_SHAPE_ELLIPSE_MINOR) != 0);
2579 13 : pasCurves[iCurve].u.EllipseByCenter.bIsComplete =
2580 13 : ((nBits & EXT_SHAPE_ELLIPSE_COMPLETE) != 0);
2581 13 : iCurve++;
2582 : }
2583 : }
2584 : else
2585 : {
2586 0 : CPLDebug("OGR", "unhandled segmentType = %d", nSegmentType);
2587 : }
2588 : }
2589 :
2590 97 : nCurves = iCurve;
2591 : }
2592 :
2593 : /* --------------------------------------------------------------------
2594 : */
2595 : /* Build corresponding OGR objects. */
2596 : /* --------------------------------------------------------------------
2597 : */
2598 1289 : if (nSHPType == SHPT_ARC || nSHPType == SHPT_ARCZ ||
2599 770 : nSHPType == SHPT_ARCM || nSHPType == SHPT_ARCZM)
2600 : {
2601 : /* --------------------------------------------------------------------
2602 : */
2603 : /* Arc - As LineString */
2604 : /* --------------------------------------------------------------------
2605 : */
2606 519 : if (nParts == 1)
2607 : {
2608 396 : if (nCurves > 0)
2609 : {
2610 38 : *ppoGeom = OGRShapeCreateCompoundCurve(
2611 : 0, nPoints, pasCurves, nCurves, 0, padfX, padfY,
2612 : bHasZ ? padfZ : nullptr, padfM, nullptr);
2613 : }
2614 : else
2615 : {
2616 358 : OGRLineString *poLine = new OGRLineString();
2617 358 : *ppoGeom = poLine;
2618 :
2619 358 : poLine->setPoints(nPoints, padfX, padfY, padfZ, padfM);
2620 : }
2621 : }
2622 :
2623 : /* --------------------------------------------------------------------
2624 : */
2625 : /* Arc - As MultiLineString */
2626 : /* --------------------------------------------------------------------
2627 : */
2628 : else
2629 : {
2630 123 : if (nCurves > 0)
2631 : {
2632 5 : OGRMultiCurve *poMulti = new OGRMultiCurve;
2633 5 : *ppoGeom = poMulti;
2634 :
2635 5 : int iCurveIdx = 0;
2636 18 : for (int i = 0; i < nParts; i++)
2637 : {
2638 13 : const int nVerticesInThisPart =
2639 13 : i == nParts - 1
2640 13 : ? nPoints - panPartStart[i]
2641 8 : : panPartStart[i + 1] - panPartStart[i];
2642 :
2643 13 : OGRCurve *poCurve = OGRShapeCreateCompoundCurve(
2644 13 : panPartStart[i], nVerticesInThisPart, pasCurves,
2645 : nCurves, iCurveIdx, padfX, padfY,
2646 : bHasZ ? padfZ : nullptr, padfM, &iCurveIdx);
2647 26 : if (poCurve == nullptr || poMulti->addGeometryDirectly(
2648 13 : poCurve) != OGRERR_NONE)
2649 : {
2650 0 : delete poCurve;
2651 0 : delete poMulti;
2652 0 : *ppoGeom = nullptr;
2653 0 : break;
2654 : }
2655 : }
2656 : }
2657 : else
2658 : {
2659 118 : OGRMultiLineString *poMulti = new OGRMultiLineString;
2660 118 : *ppoGeom = poMulti;
2661 :
2662 354 : for (int i = 0; i < nParts; i++)
2663 : {
2664 236 : OGRLineString *poLine = new OGRLineString;
2665 236 : const int nVerticesInThisPart =
2666 236 : i == nParts - 1
2667 236 : ? nPoints - panPartStart[i]
2668 118 : : panPartStart[i + 1] - panPartStart[i];
2669 :
2670 236 : poLine->setPoints(
2671 236 : nVerticesInThisPart, padfX + panPartStart[i],
2672 236 : padfY + panPartStart[i], padfZ + panPartStart[i],
2673 0 : padfM != nullptr ? padfM + panPartStart[i]
2674 : : nullptr);
2675 :
2676 236 : poMulti->addGeometryDirectly(poLine);
2677 : }
2678 : }
2679 519 : }
2680 : } // ARC.
2681 :
2682 : /* --------------------------------------------------------------------
2683 : */
2684 : /* Polygon */
2685 : /* --------------------------------------------------------------------
2686 : */
2687 770 : else if (nSHPType == SHPT_POLYGON || nSHPType == SHPT_POLYGONZ ||
2688 240 : nSHPType == SHPT_POLYGONM || nSHPType == SHPT_POLYGONZM)
2689 : {
2690 534 : if (nCurves > 0 && nParts != 0)
2691 : {
2692 53 : if (nParts == 1)
2693 : {
2694 41 : OGRCurvePolygon *poOGRPoly = new OGRCurvePolygon;
2695 41 : *ppoGeom = poOGRPoly;
2696 41 : const int nVerticesInThisPart = nPoints - panPartStart[0];
2697 :
2698 41 : OGRCurve *poRing = OGRShapeCreateCompoundCurve(
2699 : panPartStart[0], nVerticesInThisPart, pasCurves,
2700 : nCurves, 0, padfX, padfY, bHasZ ? padfZ : nullptr,
2701 : padfM, nullptr);
2702 82 : if (poRing == nullptr ||
2703 41 : poOGRPoly->addRingDirectly(poRing) != OGRERR_NONE)
2704 : {
2705 0 : delete poRing;
2706 0 : delete poOGRPoly;
2707 0 : *ppoGeom = nullptr;
2708 : }
2709 : }
2710 : else
2711 : {
2712 12 : OGRGeometry *poOGR = nullptr;
2713 : OGRCurvePolygon **tabPolygons =
2714 12 : new OGRCurvePolygon *[nParts];
2715 :
2716 12 : int iCurveIdx = 0;
2717 42 : for (int i = 0; i < nParts; i++)
2718 : {
2719 30 : tabPolygons[i] = new OGRCurvePolygon();
2720 30 : const int nVerticesInThisPart =
2721 30 : i == nParts - 1
2722 30 : ? nPoints - panPartStart[i]
2723 18 : : panPartStart[i + 1] - panPartStart[i];
2724 :
2725 30 : OGRCurve *poRing = OGRShapeCreateCompoundCurve(
2726 30 : panPartStart[i], nVerticesInThisPart, pasCurves,
2727 : nCurves, iCurveIdx, padfX, padfY,
2728 : bHasZ ? padfZ : nullptr, padfM, &iCurveIdx);
2729 60 : if (poRing == nullptr ||
2730 30 : tabPolygons[i]->addRingDirectly(poRing) !=
2731 : OGRERR_NONE)
2732 : {
2733 0 : delete poRing;
2734 0 : for (; i >= 0; --i)
2735 0 : delete tabPolygons[i];
2736 0 : delete[] tabPolygons;
2737 0 : tabPolygons = nullptr;
2738 0 : *ppoGeom = nullptr;
2739 0 : break;
2740 : }
2741 : }
2742 :
2743 12 : if (tabPolygons != nullptr)
2744 : {
2745 12 : int isValidGeometry = FALSE;
2746 12 : const char *papszOptions[] = {"METHOD=ONLY_CCW",
2747 : nullptr};
2748 12 : poOGR = OGRGeometryFactory::organizePolygons(
2749 : reinterpret_cast<OGRGeometry **>(tabPolygons),
2750 : nParts, &isValidGeometry, papszOptions);
2751 :
2752 12 : if (!isValidGeometry)
2753 : {
2754 0 : CPLError(CE_Warning, CPLE_AppDefined,
2755 : "Geometry of polygon cannot be translated "
2756 : "to Simple Geometry. All polygons will "
2757 : "be contained in a multipolygon.");
2758 : }
2759 :
2760 12 : *ppoGeom = poOGR;
2761 12 : delete[] tabPolygons;
2762 : }
2763 53 : }
2764 : }
2765 481 : else if (nParts != 0)
2766 : {
2767 480 : if (nParts == 1)
2768 : {
2769 360 : OGRPolygon *poOGRPoly = new OGRPolygon;
2770 360 : *ppoGeom = poOGRPoly;
2771 360 : OGRLinearRing *poRing = new OGRLinearRing;
2772 360 : int nVerticesInThisPart = nPoints - panPartStart[0];
2773 :
2774 360 : poRing->setPoints(
2775 360 : nVerticesInThisPart, padfX + panPartStart[0],
2776 360 : padfY + panPartStart[0], padfZ + panPartStart[0],
2777 2 : padfM != nullptr ? padfM + panPartStart[0] : nullptr);
2778 :
2779 360 : if (poOGRPoly->addRingDirectly(poRing) != OGRERR_NONE)
2780 : {
2781 0 : delete poRing;
2782 0 : delete poOGRPoly;
2783 0 : *ppoGeom = nullptr;
2784 : }
2785 : }
2786 : else
2787 : {
2788 120 : OGRGeometry *poOGR = nullptr;
2789 120 : OGRPolygon **tabPolygons = new OGRPolygon *[nParts];
2790 :
2791 493 : for (int i = 0; i < nParts; i++)
2792 : {
2793 373 : tabPolygons[i] = new OGRPolygon();
2794 373 : OGRLinearRing *poRing = new OGRLinearRing;
2795 373 : const int nVerticesInThisPart =
2796 373 : i == nParts - 1
2797 373 : ? nPoints - panPartStart[i]
2798 253 : : panPartStart[i + 1] - panPartStart[i];
2799 :
2800 373 : poRing->setPoints(
2801 373 : nVerticesInThisPart, padfX + panPartStart[i],
2802 373 : padfY + panPartStart[i], padfZ + panPartStart[i],
2803 0 : padfM != nullptr ? padfM + panPartStart[i]
2804 : : nullptr);
2805 373 : if (tabPolygons[i]->addRingDirectly(poRing) !=
2806 : OGRERR_NONE)
2807 : {
2808 0 : delete poRing;
2809 0 : for (; i >= 0; --i)
2810 0 : delete tabPolygons[i];
2811 0 : delete[] tabPolygons;
2812 0 : tabPolygons = nullptr;
2813 0 : *ppoGeom = nullptr;
2814 0 : break;
2815 : }
2816 : }
2817 :
2818 120 : if (tabPolygons != nullptr)
2819 : {
2820 120 : int isValidGeometry = FALSE;
2821 : // The outer ring is supposed to be clockwise oriented
2822 : // If it is not, then use the default/slow method.
2823 120 : const char *papszOptions[] = {
2824 120 : !(tabPolygons[0]->getExteriorRing()->isClockwise())
2825 120 : ? "METHOD=DEFAULT"
2826 : : "METHOD=ONLY_CCW",
2827 120 : nullptr};
2828 120 : poOGR = OGRGeometryFactory::organizePolygons(
2829 : reinterpret_cast<OGRGeometry **>(tabPolygons),
2830 : nParts, &isValidGeometry, papszOptions);
2831 :
2832 120 : if (!isValidGeometry)
2833 : {
2834 0 : CPLError(CE_Warning, CPLE_AppDefined,
2835 : "Geometry of polygon cannot be translated "
2836 : "to Simple Geometry. All polygons will be "
2837 : "contained in a multipolygon.");
2838 : }
2839 :
2840 120 : *ppoGeom = poOGR;
2841 120 : delete[] tabPolygons;
2842 : }
2843 : }
2844 : }
2845 : else
2846 : {
2847 1 : *ppoGeom = new OGRPolygon();
2848 534 : }
2849 : } // Polygon.
2850 :
2851 : /* --------------------------------------------------------------------
2852 : */
2853 : /* Multipatch */
2854 : /* --------------------------------------------------------------------
2855 : */
2856 236 : else if (bIsMultiPatch)
2857 : {
2858 236 : *ppoGeom =
2859 236 : OGRCreateFromMultiPatch(nParts, panPartStart, panPartType,
2860 : nPoints, padfX, padfY, padfZ);
2861 : }
2862 :
2863 1289 : CPLFree(panPartStart);
2864 1289 : CPLFree(panPartType);
2865 1289 : CPLFree(padfX);
2866 1289 : CPLFree(padfY);
2867 1289 : CPLFree(padfZ);
2868 1289 : CPLFree(padfM);
2869 1289 : CPLFree(pasCurves);
2870 :
2871 1289 : if (*ppoGeom != nullptr)
2872 : {
2873 1289 : if (!bHasZ)
2874 554 : (*ppoGeom)->set3D(FALSE);
2875 : }
2876 :
2877 1289 : return *ppoGeom != nullptr ? OGRERR_NONE : OGRERR_FAILURE;
2878 : }
2879 :
2880 : /* ==================================================================== */
2881 : /* Extract vertices for a MultiPoint. */
2882 : /* ==================================================================== */
2883 478 : else if (nSHPType == SHPT_MULTIPOINT || nSHPType == SHPT_MULTIPOINTM ||
2884 242 : nSHPType == SHPT_MULTIPOINTZ || nSHPType == SHPT_MULTIPOINTZM)
2885 : {
2886 236 : GInt32 nPoints = 0;
2887 236 : memcpy(&nPoints, pabyShape + 36, 4);
2888 236 : CPL_LSBPTR32(&nPoints);
2889 :
2890 236 : if (nPoints < 0 || nPoints > 50 * 1000 * 1000)
2891 : {
2892 0 : CPLError(CE_Failure, CPLE_AppDefined,
2893 : "Corrupted Shape : nPoints=%d.", nPoints);
2894 0 : return OGRERR_FAILURE;
2895 : }
2896 :
2897 236 : const GInt32 nOffsetZ = 40 + 2 * 8 * nPoints + 2 * 8;
2898 236 : GInt32 nOffsetM = 0;
2899 236 : if (bHasM)
2900 0 : nOffsetM = bHasZ ? nOffsetZ + 2 * 8 * 8 * nPoints : nOffsetZ;
2901 :
2902 236 : OGRMultiPoint *poMultiPt = new OGRMultiPoint;
2903 236 : *ppoGeom = poMultiPt;
2904 :
2905 708 : for (int i = 0; i < nPoints; i++)
2906 : {
2907 472 : OGRPoint *poPt = new OGRPoint;
2908 :
2909 : // Copy X.
2910 472 : double x = 0.0;
2911 472 : memcpy(&x, pabyShape + 40 + i * 16, 8);
2912 472 : CPL_LSBPTR64(&x);
2913 472 : poPt->setX(x);
2914 :
2915 : // Copy Y.
2916 472 : double y = 0.0;
2917 472 : memcpy(&y, pabyShape + 40 + i * 16 + 8, 8);
2918 472 : CPL_LSBPTR64(&y);
2919 472 : poPt->setY(y);
2920 :
2921 : // Copy Z.
2922 472 : if (bHasZ)
2923 : {
2924 236 : double z = 0.0;
2925 236 : memcpy(&z, pabyShape + nOffsetZ + i * 8, 8);
2926 236 : CPL_LSBPTR64(&z);
2927 236 : poPt->setZ(z);
2928 : }
2929 :
2930 : // Copy M.
2931 472 : if (bHasM)
2932 : {
2933 0 : double m = 0.0;
2934 0 : memcpy(&m, pabyShape + nOffsetM + i * 8, 8);
2935 0 : CPL_LSBPTR64(&m);
2936 0 : poPt->setM(m);
2937 : }
2938 :
2939 472 : poMultiPt->addGeometryDirectly(poPt);
2940 : }
2941 :
2942 236 : poMultiPt->set3D(bHasZ);
2943 236 : poMultiPt->setMeasured(bHasM);
2944 :
2945 236 : return OGRERR_NONE;
2946 : }
2947 :
2948 : /* ==================================================================== */
2949 : /* Extract vertices for a point. */
2950 : /* ==================================================================== */
2951 242 : else if (nSHPType == SHPT_POINT || nSHPType == SHPT_POINTM ||
2952 0 : nSHPType == SHPT_POINTZ || nSHPType == SHPT_POINTZM)
2953 : {
2954 242 : if (nBytes < 4 + 8 + 8 + ((bHasZ) ? 8 : 0) + ((bHasM) ? 8 : 0))
2955 : {
2956 0 : CPLError(CE_Failure, CPLE_AppDefined,
2957 : "Corrupted Shape : nBytes=%d, nSHPType=%d", nBytes,
2958 : nSHPType);
2959 0 : return OGRERR_FAILURE;
2960 : }
2961 :
2962 242 : double dfX = 0.0;
2963 242 : double dfY = 0.0;
2964 :
2965 242 : memcpy(&dfX, pabyShape + 4, 8);
2966 242 : memcpy(&dfY, pabyShape + 4 + 8, 8);
2967 :
2968 242 : CPL_LSBPTR64(&dfX);
2969 242 : CPL_LSBPTR64(&dfY);
2970 : // int nOffset = 20 + 8;
2971 :
2972 242 : double dfZ = 0.0;
2973 242 : if (bHasZ)
2974 : {
2975 118 : memcpy(&dfZ, pabyShape + 4 + 16, 8);
2976 118 : CPL_LSBPTR64(&dfZ);
2977 : }
2978 :
2979 242 : double dfM = 0.0;
2980 242 : if (bHasM)
2981 : {
2982 0 : memcpy(&dfM, pabyShape + 4 + 16 + ((bHasZ) ? 8 : 0), 8);
2983 0 : CPL_LSBPTR64(&dfM);
2984 : }
2985 :
2986 242 : if (bHasZ && bHasM)
2987 : {
2988 0 : *ppoGeom = new OGRPoint(dfX, dfY, dfZ, dfM);
2989 : }
2990 242 : else if (bHasZ)
2991 : {
2992 118 : *ppoGeom = new OGRPoint(dfX, dfY, dfZ);
2993 : }
2994 124 : else if (bHasM)
2995 : {
2996 0 : OGRPoint *poPoint = new OGRPoint(dfX, dfY);
2997 0 : poPoint->setM(dfM);
2998 0 : *ppoGeom = poPoint;
2999 : }
3000 : else
3001 : {
3002 124 : *ppoGeom = new OGRPoint(dfX, dfY);
3003 : }
3004 :
3005 242 : return OGRERR_NONE;
3006 : }
3007 :
3008 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unsupported geometry type: %d",
3009 : nSHPType);
3010 :
3011 0 : return OGRERR_FAILURE;
3012 : }
3013 :
3014 : #ifdef HAVE_WFLAG_UNREACHABLE_CODE_AGGRESSIVE
3015 : #if defined(__clang__)
3016 : #pragma clang diagnostic pop
3017 : #endif
3018 : #endif
|
