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