Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: The OGRCurve geometry class.
5 : * Author: Frank Warmerdam, warmerda@home.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1999, Frank Warmerdam
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "ogr_geometry.h"
14 : #include "ogr_p.h"
15 :
16 : //! @cond Doxygen_Suppress
17 :
18 : /************************************************************************/
19 : /* operator=( const OGRCurve& ) */
20 : /************************************************************************/
21 :
22 16 : OGRCurve &OGRCurve::operator=(const OGRCurve &other)
23 : {
24 16 : if (this != &other)
25 : {
26 16 : OGRGeometry::operator=(other);
27 : }
28 16 : return *this;
29 : }
30 :
31 : //! @endcond
32 :
33 : /************************************************************************/
34 : /* getDimension() */
35 : /************************************************************************/
36 :
37 278 : int OGRCurve::getDimension() const
38 :
39 : {
40 278 : return 1;
41 : }
42 :
43 : /************************************************************************/
44 : /* get_IsClosed() */
45 : /************************************************************************/
46 :
47 : /**
48 : * \brief Return TRUE if curve is closed.
49 : *
50 : * Tests if a curve is closed. A curve is closed if its start point is
51 : * equal to its end point.
52 : *
53 : * For equality tests, the M dimension is ignored.
54 : *
55 : * This method relates to the SFCOM ICurve::get_IsClosed() method.
56 : *
57 : * @return TRUE if closed, else FALSE.
58 : */
59 :
60 228046 : int OGRCurve::get_IsClosed() const
61 :
62 : {
63 456092 : OGRPoint oStartPoint;
64 228046 : StartPoint(&oStartPoint);
65 :
66 456091 : OGRPoint oEndPoint;
67 228046 : EndPoint(&oEndPoint);
68 :
69 228046 : if (oStartPoint.Is3D() && oEndPoint.Is3D())
70 : {
71 : // XYZ type
72 40788 : if (oStartPoint.getX() == oEndPoint.getX() &&
73 40788 : oStartPoint.getY() == oEndPoint.getY() &&
74 19522 : oStartPoint.getZ() == oEndPoint.getZ())
75 : {
76 18923 : return TRUE;
77 : }
78 : else
79 1769 : return FALSE;
80 : }
81 :
82 : // one of the points is 3D
83 414706 : else if (((oStartPoint.Is3D() & oEndPoint.Is3D()) == 0) &&
84 207354 : ((oStartPoint.Is3D() | oEndPoint.Is3D()) == 1))
85 : {
86 0 : return FALSE;
87 : }
88 :
89 : else
90 : {
91 : // XY type
92 411670 : if (oStartPoint.getX() == oEndPoint.getX() &&
93 204318 : oStartPoint.getY() == oEndPoint.getY())
94 203440 : return TRUE;
95 : else
96 3913 : return FALSE;
97 : }
98 : }
99 :
100 : /**
101 : * \fn double OGRCurve::get_Length() const;
102 : *
103 : * \brief Returns the length of the curve.
104 : *
105 : * This method relates to the SFCOM ICurve::get_Length() method.
106 : *
107 : * @return the length of the curve, zero if the curve hasn't been
108 : * initialized.
109 : *
110 : * @see get_GeodesicLength() for an alternative method returning lengths
111 : * computed on the ellipsoid, and in meters.
112 : */
113 :
114 : /**
115 : * \fn double OGRCurve::get_GeodesicLength(const OGRSpatialReference* poSRSOverride = nullptr) const;
116 : *
117 : * \brief Get the length of the curve, considered as a geodesic line on the
118 : * underlying ellipsoid of the SRS attached to the geometry.
119 : *
120 : * The returned length will always be in meters.
121 : *
122 : * <a href="https://geographiclib.sourceforge.io/html/python/geodesics.html">Geodesics</a>
123 : * follow the shortest route on the surface of the ellipsoid.
124 : *
125 : * If the geometry' SRS is not a geographic one, geometries are reprojected to
126 : * the underlying geographic SRS of the geometry' SRS.
127 : * OGRSpatialReference::GetDataAxisToSRSAxisMapping() is honored.
128 : *
129 : * Note that geometries with circular arcs will be linearized in their original
130 : * coordinate space first, so the resulting geodesic length will be an
131 : * approximation.
132 : *
133 : * @param poSRSOverride If not null, overrides OGRGeometry::getSpatialReference()
134 : * @return the length of the geometry in meters, or a negative value in case
135 : * of error.
136 : *
137 : * @see get_Length() for an alternative method returning areas computed in
138 : * 2D Cartesian space.
139 : *
140 : * @since GDAL 3.10
141 : */
142 :
143 : /**
144 : * \fn void OGRCurve::StartPoint( OGRPoint * poPoint ) const;
145 : *
146 : * \brief Return the curve start point.
147 : *
148 : * This method relates to the SF COM ICurve::get_StartPoint() method.
149 : *
150 : * @param poPoint the point to be assigned the start location.
151 : */
152 :
153 : /**
154 : * \fn void OGRCurve::EndPoint( OGRPoint * poPoint ) const;
155 : *
156 : * \brief Return the curve end point.
157 : *
158 : * This method relates to the SF COM ICurve::get_EndPoint() method.
159 : *
160 : * @param poPoint the point to be assigned the end location.
161 : */
162 :
163 : /**
164 : * \fn void OGRCurve::Value( double dfDistance, OGRPoint * poPoint ) const;
165 : *
166 : * \brief Fetch point at given distance along curve.
167 : *
168 : * This method relates to the SF COM ICurve::get_Value() method.
169 : *
170 : * This function is the same as the C function OGR_G_Value().
171 : *
172 : * @param dfDistance distance along the curve at which to sample position.
173 : * This distance should be between zero and get_Length()
174 : * for this curve.
175 : * @param poPoint the point to be assigned the curve position.
176 : */
177 :
178 : /**
179 : * \fn OGRLineString* OGRCurve::CurveToLine( double dfMaxAngleStepSizeDegrees,
180 : * const char* const* papszOptions ) const;
181 : *
182 : * \brief Return a linestring from a curve geometry.
183 : *
184 : * The returned geometry is a new instance whose ownership belongs to the
185 : * caller.
186 : *
187 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
188 : * used. This is currently 4 degrees unless the user has overridden the
189 : * value with the OGR_ARC_STEPSIZE configuration variable.
190 : *
191 : * This method relates to the ISO SQL/MM Part 3 ICurve::CurveToLine() method.
192 : *
193 : * This function is the same as C function OGR_G_CurveToLine().
194 : *
195 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
196 : * arc, zero to use the default setting.
197 : * @param papszOptions options as a null-terminated list of strings or NULL.
198 : * See OGRGeometryFactory::curveToLineString() for valid
199 : * options.
200 : *
201 : * @return a line string approximating the curve
202 : *
203 : */
204 :
205 : /**
206 : * \fn int OGRCurve::getNumPoints() const;
207 : *
208 : * \brief Return the number of points of a curve geometry.
209 : *
210 : *
211 : * This method, as a method of OGRCurve, does not relate to a standard.
212 : * For circular strings or linestrings, it returns the number of points,
213 : * conforming to SF COM NumPoints().
214 : * For compound curves, it returns the sum of the number of points of each
215 : * of its components (non including intermediate starting/ending points of
216 : * the different parts).
217 : *
218 : * @return the number of points of the curve.
219 : *
220 : */
221 :
222 : /**
223 : * \fn OGRPointIterator* OGRCurve::getPointIterator() const;
224 : *
225 : * \brief Returns a point iterator over the curve.
226 : *
227 : * The curve must not be modified while an iterator exists on it.
228 : *
229 : * The iterator must be destroyed with OGRPointIterator::destroy().
230 : *
231 : * @return a point iterator over the curve.
232 : *
233 : */
234 :
235 : /**
236 : * \fn double OGRCurve::get_Area() const;
237 : *
238 : * \brief Get the area of the (closed) curve.
239 : *
240 : * This method is designed to be used by OGRCurvePolygon::get_Area().
241 : *
242 : * @return the area of the geometry in square units of the spatial reference
243 : * system in use.
244 : *
245 : * @see get_GeodesicArea() for an alternative method returning areas
246 : * computed on the ellipsoid, and in square meters.
247 : *
248 : */
249 :
250 : /**
251 : * \fn double OGRCurve::get_GeodesicArea(const OGRSpatialReference* poSRSOverride = nullptr) const;
252 : *
253 : * \brief Get the area of the (closed) curve, considered as a surface on the
254 : * underlying ellipsoid of the SRS attached to the geometry.
255 : *
256 : * This method is designed to be used by OGRCurvePolygon::get_GeodesicArea().
257 : *
258 : * The returned area will always be in square meters, and assumes that
259 : * polygon edges describe geodesic lines on the ellipsoid.
260 : *
261 : * <a href="https://geographiclib.sourceforge.io/html/python/geodesics.html">Geodesics</a>
262 : * follow the shortest route on the surface of the ellipsoid.
263 : *
264 : * If the geometry' SRS is not a geographic one, geometries are reprojected to
265 : * the underlying geographic SRS of the geometry' SRS.
266 : * OGRSpatialReference::GetDataAxisToSRSAxisMapping() is honored.
267 : *
268 : * Note that geometries with circular arcs will be linearized in their original
269 : * coordinate space first, so the resulting geodesic area will be an
270 : * approximation.
271 : *
272 : * @param poSRSOverride If not null, overrides OGRGeometry::getSpatialReference()
273 : * @return the area of the geometry in square meters, or a negative value in case
274 : * of error.
275 : *
276 : * @see get_Area() for an alternative method returning areas computed in
277 : * 2D Cartesian space.
278 : *
279 : * @since GDAL 3.9
280 : */
281 :
282 : /**
283 : * \fn double OGRCurve::get_AreaOfCurveSegments() const;
284 : *
285 : * \brief Get the area of the purely curve portions of a (closed) curve.
286 : *
287 : * This method is designed to be used on a closed convex curve.
288 : *
289 : * @return the area of the feature in square units of the spatial reference
290 : * system in use.
291 : *
292 : */
293 :
294 : /************************************************************************/
295 : /* IsConvex() */
296 : /************************************************************************/
297 :
298 : /**
299 : * \brief Returns if a (closed) curve forms a convex shape.
300 : *
301 : * @return TRUE if the curve forms a convex shape.
302 : *
303 : */
304 :
305 42 : OGRBoolean OGRCurve::IsConvex() const
306 : {
307 42 : bool bRet = true;
308 42 : OGRPointIterator *poPointIter = getPointIterator();
309 84 : OGRPoint p1;
310 42 : OGRPoint p2;
311 42 : if (poPointIter->getNextPoint(&p1) && poPointIter->getNextPoint(&p2))
312 : {
313 84 : OGRPoint p3;
314 168 : while (poPointIter->getNextPoint(&p3))
315 : {
316 : const double crossproduct =
317 140 : (p2.getX() - p1.getX()) * (p3.getY() - p2.getY()) -
318 140 : (p2.getY() - p1.getY()) * (p3.getX() - p2.getX());
319 140 : if (crossproduct > 0)
320 : {
321 14 : bRet = false;
322 14 : break;
323 : }
324 126 : p1.setX(p2.getX());
325 126 : p1.setY(p2.getY());
326 126 : p2.setX(p3.getX());
327 126 : p2.setY(p3.getY());
328 : }
329 : }
330 42 : delete poPointIter;
331 84 : return bRet;
332 : }
333 :
334 : /************************************************************************/
335 : /* CastToCompoundCurve() */
336 : /************************************************************************/
337 :
338 : /**
339 : * \brief Cast to compound curve
340 : *
341 : * The passed in geometry is consumed and a new one returned (or NULL in case
342 : * of failure)
343 : *
344 : * @param poCurve the input geometry - ownership is passed to the method.
345 : * @return new geometry
346 : *
347 : */
348 :
349 32 : OGRCompoundCurve *OGRCurve::CastToCompoundCurve(OGRCurve *poCurve)
350 : {
351 32 : OGRCompoundCurve *poCC = new OGRCompoundCurve();
352 32 : if (wkbFlatten(poCurve->getGeometryType()) == wkbLineString)
353 23 : poCurve = CastToLineString(poCurve);
354 32 : if (!poCurve->IsEmpty() && poCC->addCurveDirectly(poCurve) != OGRERR_NONE)
355 : {
356 2 : delete poCC;
357 2 : delete poCurve;
358 2 : return nullptr;
359 : }
360 30 : poCC->assignSpatialReference(poCurve->getSpatialReference());
361 30 : return poCC;
362 : }
363 :
364 : /************************************************************************/
365 : /* CastToLineString() */
366 : /************************************************************************/
367 :
368 : /**
369 : * \brief Cast to linestring
370 : *
371 : * The passed in geometry is consumed and a new one returned (or NULL in case
372 : * of failure)
373 : *
374 : * @param poCurve the input geometry - ownership is passed to the method.
375 : * @return new geometry.
376 : *
377 : */
378 :
379 224 : OGRLineString *OGRCurve::CastToLineString(OGRCurve *poCurve)
380 : {
381 224 : OGRCurveCasterToLineString pfn = poCurve->GetCasterToLineString();
382 224 : return pfn(poCurve);
383 : }
384 :
385 : /************************************************************************/
386 : /* CastToLinearRing() */
387 : /************************************************************************/
388 :
389 : /**
390 : * \brief Cast to linear ring
391 : *
392 : * The passed in geometry is consumed and a new one returned (or NULL in case
393 : * of failure)
394 : *
395 : * @param poCurve the input geometry - ownership is passed to the method.
396 : * @return new geometry.
397 : *
398 : */
399 :
400 551 : OGRLinearRing *OGRCurve::CastToLinearRing(OGRCurve *poCurve)
401 : {
402 551 : OGRCurveCasterToLinearRing pfn = poCurve->GetCasterToLinearRing();
403 551 : return pfn(poCurve);
404 : }
405 :
406 : /************************************************************************/
407 : /* ContainsPoint() */
408 : /************************************************************************/
409 :
410 : /**
411 : * \brief Returns if a point is contained in a (closed) curve.
412 : *
413 : * Final users should use OGRGeometry::Contains() instead.
414 : *
415 : * @param p the point to test
416 : * @return TRUE if it is inside the curve, FALSE otherwise or -1 if unknown.
417 : *
418 : */
419 :
420 17 : int OGRCurve::ContainsPoint(CPL_UNUSED const OGRPoint *p) const
421 : {
422 17 : return -1;
423 : }
424 :
425 : /************************************************************************/
426 : /* IntersectsPoint() */
427 : /************************************************************************/
428 :
429 : /**
430 : * \brief Returns if a point intersects a (closed) curve.
431 : *
432 : * Final users should use OGRGeometry::Intersects() instead.
433 : *
434 : * @param p the point to test
435 : * @return TRUE if it intersects the curve, FALSE otherwise or -1 if unknown.
436 : *
437 : */
438 :
439 1 : int OGRCurve::IntersectsPoint(CPL_UNUSED const OGRPoint *p) const
440 : {
441 1 : return -1;
442 : }
443 :
444 : /************************************************************************/
445 : /* ~OGRPointIterator() */
446 : /************************************************************************/
447 :
448 : OGRPointIterator::~OGRPointIterator() = default;
449 :
450 : /**
451 : * \fn OGRBoolean OGRPointIterator::getNextPoint(OGRPoint* p);
452 : *
453 : * \brief Returns the next point followed by the iterator.
454 : *
455 : * @param p point to fill.
456 : *
457 : * @return TRUE in case of success, or FALSE if the end of the curve is reached.
458 : *
459 : */
460 :
461 : /************************************************************************/
462 : /* destroy() */
463 : /************************************************************************/
464 :
465 : /**
466 : * \brief Destroys a point iterator.
467 : *
468 : */
469 0 : void OGRPointIterator::destroy(OGRPointIterator *poIter)
470 : {
471 0 : delete poIter;
472 0 : }
473 :
474 : /************************************************************************/
475 : /* OGRSimpleCurve::Iterator */
476 : /************************************************************************/
477 :
478 : OGRIteratedPoint::~OGRIteratedPoint() = default;
479 :
480 2 : void OGRIteratedPoint::setX(double xIn)
481 : {
482 2 : OGRPoint::setX(xIn);
483 2 : m_poCurve->setPoint(m_nPos, xIn, getY());
484 2 : }
485 :
486 2 : void OGRIteratedPoint::setY(double yIn)
487 : {
488 2 : OGRPoint::setY(yIn);
489 2 : m_poCurve->setPoint(m_nPos, getX(), yIn);
490 2 : }
491 :
492 2 : void OGRIteratedPoint::setZ(double zIn)
493 : {
494 2 : OGRPoint::setZ(zIn);
495 2 : m_poCurve->setZ(m_nPos, zIn);
496 2 : }
497 :
498 2 : void OGRIteratedPoint::setM(double mIn)
499 : {
500 2 : OGRPoint::setM(mIn);
501 2 : m_poCurve->setM(m_nPos, mIn);
502 2 : }
503 :
504 : struct OGRSimpleCurve::Iterator::Private
505 : {
506 : CPL_DISALLOW_COPY_ASSIGN(Private)
507 228 : Private() = default;
508 :
509 : bool m_bUpdateChecked = true;
510 : OGRIteratedPoint m_oPoint{};
511 : };
512 :
513 2528 : void OGRSimpleCurve::Iterator::update()
514 : {
515 2528 : if (!m_poPrivate->m_bUpdateChecked)
516 : {
517 1149 : OGRPoint oPointBefore;
518 1149 : m_poPrivate->m_oPoint.m_poCurve->getPoint(m_poPrivate->m_oPoint.m_nPos,
519 : &oPointBefore);
520 1149 : if (oPointBefore != m_poPrivate->m_oPoint)
521 : {
522 899 : if (m_poPrivate->m_oPoint.Is3D())
523 739 : m_poPrivate->m_oPoint.m_poCurve->set3D(true);
524 899 : if (m_poPrivate->m_oPoint.IsMeasured())
525 4 : m_poPrivate->m_oPoint.m_poCurve->setMeasured(true);
526 1798 : m_poPrivate->m_oPoint.m_poCurve->setPoint(
527 899 : m_poPrivate->m_oPoint.m_nPos, &m_poPrivate->m_oPoint);
528 : }
529 1149 : m_poPrivate->m_bUpdateChecked = true;
530 : }
531 2528 : }
532 :
533 228 : OGRSimpleCurve::Iterator::Iterator(OGRSimpleCurve *poSelf, int nPos)
534 228 : : m_poPrivate(new Private())
535 : {
536 228 : m_poPrivate->m_oPoint.m_poCurve = poSelf;
537 228 : m_poPrivate->m_oPoint.m_nPos = nPos;
538 228 : }
539 :
540 228 : OGRSimpleCurve::Iterator::~Iterator()
541 : {
542 228 : update();
543 228 : }
544 :
545 1149 : OGRIteratedPoint &OGRSimpleCurve::Iterator::operator*()
546 : {
547 1149 : update();
548 2298 : m_poPrivate->m_oPoint.m_poCurve->getPoint(m_poPrivate->m_oPoint.m_nPos,
549 1149 : &m_poPrivate->m_oPoint);
550 1149 : m_poPrivate->m_bUpdateChecked = false;
551 1149 : return m_poPrivate->m_oPoint;
552 : }
553 :
554 1151 : OGRSimpleCurve::Iterator &OGRSimpleCurve::Iterator::operator++()
555 : {
556 1151 : update();
557 1151 : ++m_poPrivate->m_oPoint.m_nPos;
558 1151 : return *this;
559 : }
560 :
561 1266 : bool OGRSimpleCurve::Iterator::operator!=(const Iterator &it) const
562 : {
563 1266 : return m_poPrivate->m_oPoint.m_nPos != it.m_poPrivate->m_oPoint.m_nPos;
564 : }
565 :
566 115 : OGRSimpleCurve::Iterator OGRSimpleCurve::begin()
567 : {
568 115 : return {this, 0};
569 : }
570 :
571 113 : OGRSimpleCurve::Iterator OGRSimpleCurve::end()
572 : {
573 113 : return {this, nPointCount};
574 : }
575 :
576 : /************************************************************************/
577 : /* OGRSimpleCurve::ConstIterator */
578 : /************************************************************************/
579 :
580 : struct OGRSimpleCurve::ConstIterator::Private
581 : {
582 : CPL_DISALLOW_COPY_ASSIGN(Private)
583 250 : Private() = default;
584 :
585 : mutable OGRIteratedPoint m_oPoint{};
586 : const OGRSimpleCurve *m_poSelf = nullptr;
587 : int m_nPos = 0;
588 : };
589 :
590 250 : OGRSimpleCurve::ConstIterator::ConstIterator(const OGRSimpleCurve *poSelf,
591 250 : int nPos)
592 250 : : m_poPrivate(new Private())
593 : {
594 250 : m_poPrivate->m_poSelf = poSelf;
595 250 : m_poPrivate->m_nPos = nPos;
596 250 : }
597 :
598 : OGRSimpleCurve::ConstIterator::~ConstIterator() = default;
599 :
600 2623 : const OGRPoint &OGRSimpleCurve::ConstIterator::operator*() const
601 : {
602 5246 : m_poPrivate->m_poSelf->getPoint(m_poPrivate->m_nPos,
603 2623 : &m_poPrivate->m_oPoint);
604 2623 : return m_poPrivate->m_oPoint;
605 : }
606 :
607 2623 : OGRSimpleCurve::ConstIterator &OGRSimpleCurve::ConstIterator::operator++()
608 : {
609 2623 : ++m_poPrivate->m_nPos;
610 2623 : return *this;
611 : }
612 :
613 2748 : bool OGRSimpleCurve::ConstIterator::operator!=(const ConstIterator &it) const
614 : {
615 2748 : return m_poPrivate->m_nPos != it.m_poPrivate->m_nPos;
616 : }
617 :
618 125 : OGRSimpleCurve::ConstIterator OGRSimpleCurve::begin() const
619 : {
620 125 : return {this, 0};
621 : }
622 :
623 125 : OGRSimpleCurve::ConstIterator OGRSimpleCurve::end() const
624 : {
625 125 : return {this, nPointCount};
626 : }
627 :
628 : /************************************************************************/
629 : /* OGRCurve::ConstIterator */
630 : /************************************************************************/
631 :
632 : struct OGRCurve::ConstIterator::Private
633 : {
634 : CPL_DISALLOW_COPY_ASSIGN(Private)
635 38 : Private() = default;
636 : Private(Private &&) = delete;
637 : Private &operator=(Private &&) = default;
638 :
639 : OGRPoint m_oPoint{};
640 : const OGRCurve *m_poCurve{};
641 : int m_nStep = 0;
642 : std::unique_ptr<OGRPointIterator> m_poIterator{};
643 : };
644 :
645 38 : OGRCurve::ConstIterator::ConstIterator(const OGRCurve *poSelf, bool bStart)
646 38 : : m_poPrivate(new Private())
647 : {
648 38 : m_poPrivate->m_poCurve = poSelf;
649 38 : if (bStart)
650 : {
651 30 : m_poPrivate->m_poIterator.reset(poSelf->getPointIterator());
652 30 : if (!m_poPrivate->m_poIterator->getNextPoint(&m_poPrivate->m_oPoint))
653 : {
654 2 : m_poPrivate->m_nStep = -1;
655 2 : m_poPrivate->m_poIterator.reset();
656 : }
657 : }
658 : else
659 : {
660 8 : m_poPrivate->m_nStep = -1;
661 : }
662 38 : }
663 :
664 0 : OGRCurve::ConstIterator::ConstIterator(ConstIterator &&oOther) noexcept
665 0 : : m_poPrivate(std::move(oOther.m_poPrivate))
666 : {
667 0 : }
668 :
669 : OGRCurve::ConstIterator &
670 4 : OGRCurve::ConstIterator::operator=(ConstIterator &&oOther)
671 : {
672 4 : m_poPrivate = std::move(oOther.m_poPrivate);
673 4 : return *this;
674 : }
675 :
676 : OGRCurve::ConstIterator::~ConstIterator() = default;
677 :
678 119 : const OGRPoint &OGRCurve::ConstIterator::operator*() const
679 : {
680 119 : return m_poPrivate->m_oPoint;
681 : }
682 :
683 83 : OGRCurve::ConstIterator &OGRCurve::ConstIterator::operator++()
684 : {
685 83 : CPLAssert(m_poPrivate->m_nStep >= 0);
686 83 : ++m_poPrivate->m_nStep;
687 83 : if (!m_poPrivate->m_poIterator->getNextPoint(&m_poPrivate->m_oPoint))
688 : {
689 6 : m_poPrivate->m_nStep = -1;
690 6 : m_poPrivate->m_poIterator.reset();
691 : }
692 83 : return *this;
693 : }
694 :
695 24 : bool OGRCurve::ConstIterator::operator!=(const ConstIterator &it) const
696 : {
697 48 : return m_poPrivate->m_poCurve != it.m_poPrivate->m_poCurve ||
698 48 : m_poPrivate->m_nStep != it.m_poPrivate->m_nStep;
699 : }
700 :
701 30 : OGRCurve::ConstIterator OGRCurve::begin() const
702 : {
703 30 : return {this, true};
704 : }
705 :
706 8 : OGRCurve::ConstIterator OGRCurve::end() const
707 : {
708 8 : return {this, false};
709 : }
710 :
711 : /************************************************************************/
712 : /* isClockwise() */
713 : /************************************************************************/
714 :
715 : /**
716 : * \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
717 : *
718 : * Assumes that the line is closed.
719 : *
720 : * @return TRUE if clockwise otherwise FALSE.
721 : */
722 :
723 18 : int OGRCurve::isClockwise() const
724 :
725 : {
726 : // WARNING: keep in sync OGRLineString::isClockwise(),
727 : // OGRCurve::isClockwise() and OGRWKBIsClockwiseRing()
728 :
729 18 : const int nPointCount = getNumPoints();
730 18 : if (nPointCount < 3)
731 0 : return TRUE;
732 :
733 18 : bool bUseFallback = false;
734 :
735 : // Find the lowest rightmost vertex.
736 36 : auto oIter = begin();
737 36 : const OGRPoint oStartPoint = *oIter;
738 36 : OGRPoint oPointBefore = oStartPoint;
739 36 : OGRPoint oPointBeforeSel;
740 36 : OGRPoint oPointSel = oStartPoint;
741 36 : OGRPoint oPointNextSel;
742 18 : bool bNextPointIsNextSel = true;
743 18 : int v = 0;
744 :
745 69 : for (int i = 1; i < nPointCount - 1; i++)
746 : {
747 51 : ++oIter;
748 102 : const OGRPoint oPointCur = *oIter;
749 51 : if (bNextPointIsNextSel)
750 : {
751 28 : oPointNextSel = oPointCur;
752 28 : bNextPointIsNextSel = false;
753 : }
754 102 : if (oPointCur.getY() < oPointSel.getY() ||
755 51 : (oPointCur.getY() == oPointSel.getY() &&
756 11 : oPointCur.getX() > oPointSel.getX()))
757 : {
758 18 : v = i;
759 18 : oPointBeforeSel = oPointBefore;
760 18 : oPointSel = oPointCur;
761 18 : bUseFallback = false;
762 18 : bNextPointIsNextSel = true;
763 : }
764 37 : else if (oPointCur.getY() == oPointSel.getY() &&
765 4 : oPointCur.getX() == oPointSel.getX())
766 : {
767 : // Two vertex with same coordinates are the lowest rightmost
768 : // vertex. Cannot use that point as the pivot (#5342).
769 4 : bUseFallback = true;
770 : }
771 51 : oPointBefore = oPointCur;
772 : }
773 36 : const OGRPoint oPointN_m2 = *oIter;
774 :
775 18 : if (bNextPointIsNextSel)
776 : {
777 8 : oPointNextSel = oPointN_m2;
778 : }
779 :
780 : // Previous.
781 18 : if (v == 0)
782 : {
783 5 : oPointBeforeSel = oPointN_m2;
784 : }
785 :
786 18 : constexpr double EPSILON = 1.0E-5;
787 40 : const auto epsilonEqual = [](double a, double b, double eps)
788 40 : { return ::fabs(a - b) < eps; };
789 :
790 20 : if (epsilonEqual(oPointBeforeSel.getX(), oPointSel.getX(), EPSILON) &&
791 2 : epsilonEqual(oPointBeforeSel.getY(), oPointSel.getY(), EPSILON))
792 : {
793 : // Don't try to be too clever by retrying with a next point.
794 : // This can lead to false results as in the case of #3356.
795 2 : bUseFallback = true;
796 : }
797 :
798 18 : const double dx0 = oPointBeforeSel.getX() - oPointSel.getX();
799 18 : const double dy0 = oPointBeforeSel.getY() - oPointSel.getY();
800 :
801 : // Following.
802 18 : if (v + 1 >= nPointCount - 1)
803 : {
804 8 : oPointNextSel = oStartPoint;
805 : }
806 :
807 20 : if (epsilonEqual(oPointNextSel.getX(), oPointSel.getX(), EPSILON) &&
808 2 : epsilonEqual(oPointNextSel.getY(), oPointSel.getY(), EPSILON))
809 : {
810 : // Don't try to be too clever by retrying with a next point.
811 : // This can lead to false results as in the case of #3356.
812 2 : bUseFallback = true;
813 : }
814 :
815 18 : const double dx1 = oPointNextSel.getX() - oPointSel.getX();
816 18 : const double dy1 = oPointNextSel.getY() - oPointSel.getY();
817 :
818 18 : const double crossproduct = dx1 * dy0 - dx0 * dy1;
819 :
820 18 : if (!bUseFallback)
821 : {
822 14 : if (crossproduct > 0) // CCW
823 5 : return FALSE;
824 9 : else if (crossproduct < 0) // CW
825 9 : return TRUE;
826 : }
827 :
828 : // This is a degenerate case: the extent of the polygon is less than EPSILON
829 : // or 2 nearly identical points were found.
830 : // Try with Green Formula as a fallback, but this is not a guarantee
831 : // as we'll probably be affected by numerical instabilities.
832 4 : oIter = begin();
833 4 : oPointBefore = oStartPoint;
834 4 : ++oIter;
835 4 : auto oPointCur = *oIter;
836 4 : double dfSum = oStartPoint.getX() * (oPointCur.getY() - oStartPoint.getY());
837 :
838 16 : for (int i = 1; i < nPointCount - 1; i++)
839 : {
840 12 : ++oIter;
841 12 : const auto &oPointNext = *oIter;
842 12 : dfSum += oPointCur.getX() * (oPointNext.getY() - oPointBefore.getY());
843 12 : oPointBefore = oPointCur;
844 12 : oPointCur = oPointNext;
845 : }
846 :
847 4 : dfSum += oPointCur.getX() * (oStartPoint.getY() - oPointBefore.getY());
848 :
849 4 : return dfSum < 0;
850 : }
851 :
852 : /**
853 : * \fn void OGRCurve::reversePoints();
854 : *
855 : * \brief Reverse point order.
856 : *
857 : * This method updates the points in this curve in place
858 : * reversing the point ordering (first for last, etc) and component ordering
859 : * for a compound curve.
860 : *
861 : * @since 3.10
862 : */
|
