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