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