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