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