1 /******************************************************************************
2  *
3  * Project:  OpenGIS Simple Features Reference Implementation
4  * Purpose:  The OGRLinearRing geometry class.
5  * Author:   Frank Warmerdam, warmerdam@pobox.com
6  *
7  ******************************************************************************
8  * Copyright (c) 1999, Frank Warmerdam
9  * Copyright (c) 2008-2014, Even Rouault <even dot rouault at spatialys.com>
10  *
11  * Permission is hereby granted, free of charge, to any person obtaining a
12  * copy of this software and associated documentation files (the "Software"),
13  * to deal in the Software without restriction, including without limitation
14  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15  * and/or sell copies of the Software, and to permit persons to whom the
16  * Software is furnished to do so, subject to the following conditions:
17  *
18  * The above copyright notice and this permission notice shall be included
19  * in all copies or substantial portions of the Software.
20  *
21  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27  * DEALINGS IN THE SOFTWARE.
28  ****************************************************************************/
29 
30 #include "cpl_port.h"
31 #include "ogr_geometry.h"
32 
33 #include <climits>
34 #include <cmath>
35 #include <cstring>
36 
37 #include "cpl_error.h"
38 #include "ogr_core.h"
39 #include "ogr_geometry.h"
40 #include "ogr_p.h"
41 
42 CPL_CVSID("$Id: ogrlinearring.cpp 3798cbe48457b7127606931896549f26507469db 2021-04-09 15:04:16 +0200 Even Rouault $")
43 
44 /************************************************************************/
45 /*                           OGRLinearRing()                            */
46 /************************************************************************/
47 
48 /** Constructor */
49 OGRLinearRing::OGRLinearRing() = default;
50 
51 /************************************************************************/
52 /*                  OGRLinearRing( const OGRLinearRing& )               */
53 /************************************************************************/
54 
55 /**
56  * \brief Copy constructor.
57  *
58  * Note: before GDAL 2.1, only the default implementation of the constructor
59  * existed, which could be unsafe to use.
60  *
61  * @since GDAL 2.1
62  */
63 
64 OGRLinearRing::OGRLinearRing( const OGRLinearRing& ) = default;
65 
66 /************************************************************************/
67 /*                          ~OGRLinearRing()                            */
68 /************************************************************************/
69 
70 OGRLinearRing::~OGRLinearRing() = default;
71 
72 /************************************************************************/
73 /*                           OGRLinearRing()                            */
74 /************************************************************************/
75 
76 /** Constructor
77  * @param poSrcRing source ring.
78  */
OGRLinearRing(OGRLinearRing * poSrcRing)79 OGRLinearRing::OGRLinearRing( OGRLinearRing * poSrcRing )
80 
81 {
82     if( poSrcRing == nullptr )
83     {
84         CPLDebug( "OGR",
85                   "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - "
86                   "passed in ring is NULL!" );
87         return;
88     }
89 
90     setNumPoints( poSrcRing->getNumPoints(), FALSE );
91 
92     memcpy( paoPoints, poSrcRing->paoPoints,
93             sizeof(OGRRawPoint) * getNumPoints() );
94 
95     if( poSrcRing->padfZ )
96     {
97         Make3D();
98         memcpy( padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints() );
99     }
100 }
101 
102 /************************************************************************/
103 /*                    operator=( const OGRLinearRing& )                 */
104 /************************************************************************/
105 
106 /**
107  * \brief Assignment operator.
108  *
109  * Note: before GDAL 2.1, only the default implementation of the operator
110  * existed, which could be unsafe to use.
111  *
112  * @since GDAL 2.1
113  */
114 
operator =(const OGRLinearRing & other)115 OGRLinearRing& OGRLinearRing::operator=( const OGRLinearRing& other )
116 {
117     if( this != &other)
118     {
119         OGRLineString::operator=( other );
120     }
121     return *this;
122 }
123 
124 /************************************************************************/
125 /*                          getGeometryName()                           */
126 /************************************************************************/
127 
getGeometryName() const128 const char * OGRLinearRing::getGeometryName() const
129 
130 {
131     return "LINEARRING";
132 }
133 
134 /************************************************************************/
135 /*                              WkbSize()                               */
136 /*                                                                      */
137 /*      Disable this method.                                            */
138 /************************************************************************/
139 
WkbSize() const140 size_t OGRLinearRing::WkbSize() const
141 
142 {
143     return 0;
144 }
145 
146 /************************************************************************/
147 /*                           importFromWkb()                            */
148 /*                                                                      */
149 /*      Disable method for this class.                                  */
150 /************************************************************************/
151 
importFromWkb(const unsigned char *,size_t,OGRwkbVariant,size_t &)152 OGRErr OGRLinearRing::importFromWkb( const unsigned char * /*pabyData*/,
153                                      size_t /*nSize*/,
154                                      OGRwkbVariant /*eWkbVariant*/,
155                                      size_t& /* nBytesConsumedOut */ )
156 
157 {
158     return OGRERR_UNSUPPORTED_OPERATION;
159 }
160 
161 /************************************************************************/
162 /*                            exportToWkb()                             */
163 /*                                                                      */
164 /*      Disable method for this class.                                  */
165 /************************************************************************/
166 
exportToWkb(CPL_UNUSED OGRwkbByteOrder eByteOrder,CPL_UNUSED unsigned char * pabyData,CPL_UNUSED OGRwkbVariant eWkbVariant) const167 OGRErr OGRLinearRing::exportToWkb( CPL_UNUSED OGRwkbByteOrder eByteOrder,
168                                    CPL_UNUSED unsigned char * pabyData,
169                                    CPL_UNUSED OGRwkbVariant eWkbVariant ) const
170 
171 {
172     return OGRERR_UNSUPPORTED_OPERATION;
173 }
174 
175 /************************************************************************/
176 /*                           _importFromWkb()                           */
177 /*                                                                      */
178 /*      Helper method for OGRPolygon.  NOT A NORMAL importFromWkb()     */
179 /*      method.                                                         */
180 /************************************************************************/
181 
182 //! @cond Doxygen_Suppress
_importFromWkb(OGRwkbByteOrder eByteOrder,int _flags,const unsigned char * pabyData,size_t nBytesAvailable,size_t & nBytesConsumedOut)183 OGRErr OGRLinearRing::_importFromWkb( OGRwkbByteOrder eByteOrder, int _flags,
184                                       const unsigned char * pabyData,
185                                       size_t nBytesAvailable,
186                                       size_t& nBytesConsumedOut )
187 
188 {
189     nBytesConsumedOut = 0;
190     if( nBytesAvailable < 4 && nBytesAvailable != static_cast<size_t>(-1) )
191         return OGRERR_NOT_ENOUGH_DATA;
192 
193 /* -------------------------------------------------------------------- */
194 /*      Get the vertex count.                                           */
195 /* -------------------------------------------------------------------- */
196     int nNewNumPoints = 0;
197 
198     memcpy( &nNewNumPoints, pabyData, 4 );
199 
200     if( OGR_SWAP( eByteOrder ) )
201         nNewNumPoints = CPL_SWAP32(nNewNumPoints);
202 
203     // Check if the wkb stream buffer is big enough to store
204     // fetched number of points.
205     // 16, 24, or 32 - size of point structure.
206     size_t nPointSize = 0;
207     if( (_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED) )
208         nPointSize = 32;
209     else if( (_flags & OGR_G_3D) || (_flags & OGR_G_MEASURED) )
210         nPointSize = 24;
211     else
212         nPointSize = 16;
213 
214     if( nNewNumPoints < 0 ||
215         static_cast<size_t>(nNewNumPoints) > std::numeric_limits<size_t>::max() / nPointSize )
216     {
217         return OGRERR_CORRUPT_DATA;
218     }
219     const size_t nBufferMinSize = nPointSize * nNewNumPoints;
220     if( nBytesAvailable != static_cast<size_t>(-1) && nBufferMinSize > nBytesAvailable - 4 )
221     {
222         CPLError( CE_Failure, CPLE_AppDefined,
223                   "Length of input WKB is too small" );
224         return OGRERR_NOT_ENOUGH_DATA;
225     }
226 
227     // (Re)Allocation of paoPoints buffer.
228     setNumPoints( nNewNumPoints, FALSE );
229 
230     if( _flags & OGR_G_3D )
231         Make3D();
232     else
233         Make2D();
234 
235     if( _flags & OGR_G_MEASURED )
236         AddM();
237     else
238         RemoveM();
239 
240 
241     nBytesConsumedOut =  4 + nPointCount * nPointSize;
242 
243 /* -------------------------------------------------------------------- */
244 /*      Get the vertices                                                */
245 /* -------------------------------------------------------------------- */
246     if( (flags & OGR_G_3D) && (flags & OGR_G_MEASURED) )
247     {
248         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
249         {
250             memcpy( &(paoPoints[i].x), pabyData + 4 + 32 * i, 8 );
251             memcpy( &(paoPoints[i].y), pabyData + 4 + 32 * i + 8, 8 );
252             memcpy( padfZ + i, pabyData + 4 + 32 * i + 16, 8 );
253             memcpy( padfM + i, pabyData + 4 + 32 * i + 24, 8 );
254         }
255     }
256     else if( flags & OGR_G_MEASURED )
257     {
258         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
259         {
260             memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
261             memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
262             memcpy( padfM + i, pabyData + 4 + 24 * i + 16, 8 );
263         }
264     }
265     else if( flags & OGR_G_3D )
266     {
267         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
268         {
269             memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
270             memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
271             memcpy( padfZ + i, pabyData + 4 + 24 * i + 16, 8 );
272         }
273     }
274     else
275     {
276         memcpy( paoPoints, pabyData + 4, 16 * static_cast<size_t>(nPointCount) );
277     }
278 
279 /* -------------------------------------------------------------------- */
280 /*      Byte swap if needed.                                            */
281 /* -------------------------------------------------------------------- */
282     if( OGR_SWAP( eByteOrder ) )
283     {
284         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
285         {
286             CPL_SWAPDOUBLE( &(paoPoints[i].x) );
287             CPL_SWAPDOUBLE( &(paoPoints[i].y) );
288 
289             if( flags & OGR_G_3D )
290             {
291                 CPL_SWAPDOUBLE( padfZ + i );
292             }
293             if( flags & OGR_G_MEASURED )
294             {
295                 CPL_SWAPDOUBLE( padfM + i );
296             }
297         }
298     }
299 
300     return OGRERR_NONE;
301 }
302 
303 /************************************************************************/
304 /*                            _exportToWkb()                            */
305 /*                                                                      */
306 /*      Helper method for OGRPolygon.  THIS IS NOT THE NORMAL           */
307 /*      exportToWkb() METHOD.                                           */
308 /************************************************************************/
309 
_exportToWkb(OGRwkbByteOrder eByteOrder,int _flags,unsigned char * pabyData) const310 OGRErr OGRLinearRing::_exportToWkb( OGRwkbByteOrder eByteOrder, int _flags,
311                                     unsigned char * pabyData ) const
312 
313 {
314 
315 /* -------------------------------------------------------------------- */
316 /*      Copy in the raw data.                                           */
317 /* -------------------------------------------------------------------- */
318     memcpy( pabyData, &nPointCount, 4 );
319 
320 /* -------------------------------------------------------------------- */
321 /*      Copy in the raw data.                                           */
322 /* -------------------------------------------------------------------- */
323     size_t nWords = 0;
324     if( (_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED) )
325     {
326         nWords = 4 * static_cast<size_t>(nPointCount);
327         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
328         {
329             memcpy( pabyData+4+i*32, &(paoPoints[i].x), 8 );
330             memcpy( pabyData+4+i*32+8, &(paoPoints[i].y), 8 );
331             if( padfZ == nullptr )
332                 memset( pabyData+4+i*32+16, 0, 8 );
333             else
334                 memcpy( pabyData+4+i*32+16, padfZ + i, 8 );
335             if( padfM == nullptr )
336                 memset( pabyData+4+i*32+24, 0, 8 );
337             else
338                 memcpy( pabyData+4+i*32+24, padfM + i, 8 );
339         }
340     }
341     else if( _flags & OGR_G_MEASURED )
342     {
343         nWords = 3 * static_cast<size_t>(nPointCount);
344         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
345         {
346             memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
347             memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
348             if( padfM == nullptr )
349                 memset( pabyData+4+i*24+16, 0, 8 );
350             else
351                 memcpy( pabyData+4+i*24+16, padfM + i, 8 );
352         }
353     }
354     else if( _flags & OGR_G_3D )
355     {
356         nWords = 3 * static_cast<size_t>(nPointCount);
357         for( size_t i = 0; i < static_cast<size_t>(nPointCount); i++ )
358         {
359             memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
360             memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
361             if( padfZ == nullptr )
362                 memset( pabyData+4+i*24+16, 0, 8 );
363             else
364                 memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
365         }
366     }
367     else
368     {
369         nWords = 2 * static_cast<size_t>(nPointCount);
370         memcpy( pabyData+4, paoPoints, 16 * static_cast<size_t>(nPointCount) );
371     }
372 
373 /* -------------------------------------------------------------------- */
374 /*      Swap if needed.                                                 */
375 /* -------------------------------------------------------------------- */
376     if( OGR_SWAP( eByteOrder ) )
377     {
378         const int nCount = CPL_SWAP32( nPointCount );
379         memcpy( pabyData, &nCount, 4 );
380 
381         for( size_t i = 0; i < nWords; i++ )
382         {
383             CPL_SWAPDOUBLE( pabyData + 4 + 8 * i );
384         }
385     }
386 
387     return OGRERR_NONE;
388 }
389 
390 /************************************************************************/
391 /*                              _WkbSize()                              */
392 /*                                                                      */
393 /*      Helper method for OGRPolygon.  NOT THE NORMAL WkbSize() METHOD. */
394 /************************************************************************/
395 
_WkbSize(int _flags) const396 size_t OGRLinearRing::_WkbSize( int _flags ) const
397 
398 {
399     if( (_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED) )
400         return 4 + 32 * static_cast<size_t>(nPointCount);
401     else if( (_flags & OGR_G_3D) || (_flags & OGR_G_MEASURED) )
402         return 4 + 24 * static_cast<size_t>(nPointCount);
403     else
404         return 4 + 16 * static_cast<size_t>(nPointCount);
405 }
406 //! @endcond
407 
408 /************************************************************************/
409 /*                               clone()                                */
410 /*                                                                      */
411 /*      We override the OGRCurve clone() to ensure that we get the      */
412 /*      correct virtual table.                                          */
413 /************************************************************************/
414 
clone() const415 OGRLinearRing *OGRLinearRing::clone() const
416 
417 {
418     OGRLinearRing *poNewLinearRing = new OGRLinearRing();
419     poNewLinearRing->assignSpatialReference( getSpatialReference() );
420 
421     poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ, padfM );
422     poNewLinearRing->flags = flags;
423 
424     return poNewLinearRing;
425 }
426 
427 /************************************************************************/
428 /*                            epsilonEqual()                            */
429 /************************************************************************/
430 
431 constexpr double EPSILON = 1.0E-5;
432 
epsilonEqual(double a,double b,double eps)433 static inline bool epsilonEqual(double a, double b, double eps)
434 {
435     return ::fabs(a - b) < eps;
436 }
437 
438 /************************************************************************/
439 /*                            isClockwise()                             */
440 /************************************************************************/
441 
442 /**
443  * \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
444  *
445  * @return TRUE if clockwise otherwise FALSE.
446  */
447 
isClockwise() const448 int OGRLinearRing::isClockwise() const
449 
450 {
451     if( nPointCount < 2 )
452         return TRUE;
453 
454     bool bUseFallback = false;
455 
456     // Find the lowest rightmost vertex.
457     int v = 0;  // Used after for.
458     for( int i = 1; i < nPointCount - 1; i++ )
459     {
460         // => v < end.
461         if( paoPoints[i].y< paoPoints[v].y ||
462             ( paoPoints[i].y== paoPoints[v].y &&
463               paoPoints[i].x > paoPoints[v].x ) )
464         {
465             v = i;
466             bUseFallback = false;
467         }
468         else if( paoPoints[i].y == paoPoints[v].y &&
469                  paoPoints[i].x == paoPoints[v].x )
470         {
471             // Two vertex with same coordinates are the lowest rightmost
472             // vertex.  Cannot use that point as the pivot (#5342).
473             bUseFallback = true;
474         }
475     }
476 
477     // Previous.
478     int next = v - 1;
479     if( next < 0 )
480     {
481         next = nPointCount - 1 - 1;
482     }
483 
484     if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
485         epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
486     {
487         // Don't try to be too clever by retrying with a next point.
488         // This can lead to false results as in the case of #3356.
489         bUseFallback = true;
490     }
491 
492     const double dx0 = paoPoints[next].x - paoPoints[v].x;
493     const double dy0 = paoPoints[next].y - paoPoints[v].y;
494 
495     // Following.
496     next = v + 1;
497     if( next >= nPointCount - 1 )
498     {
499         next = 0;
500     }
501 
502     if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
503         epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
504     {
505         // Don't try to be too clever by retrying with a next point.
506         // This can lead to false results as in the case of #3356.
507         bUseFallback = true;
508     }
509 
510     const double dx1 = paoPoints[next].x - paoPoints[v].x;
511     const double dy1 = paoPoints[next].y - paoPoints[v].y;
512 
513     const double crossproduct = dx1 * dy0 - dx0 * dy1;
514 
515     if( !bUseFallback )
516     {
517         if( crossproduct > 0 )       // CCW
518             return FALSE;
519         else if( crossproduct < 0 )  // CW
520             return TRUE;
521     }
522 
523     // This is a degenerate case: the extent of the polygon is less than EPSILON
524     // or 2 nearly identical points were found.
525     // Try with Green Formula as a fallback, but this is not a guarantee
526     // as we'll probably be affected by numerical instabilities.
527 
528     double dfSum =
529         paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
530 
531     for( int i = 1; i < nPointCount-1; i++ )
532     {
533         dfSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
534     }
535 
536     dfSum +=
537         paoPoints[nPointCount-1].x *
538         (paoPoints[0].y - paoPoints[nPointCount-2].y);
539 
540     return dfSum < 0;
541 }
542 
543 /************************************************************************/
544 /*                             reverseWindingOrder()                    */
545 /************************************************************************/
546 
547 /** Reverse order of points.
548  */
reverseWindingOrder()549 void OGRLinearRing::reverseWindingOrder()
550 
551 {
552     OGRPoint pointA;
553     OGRPoint pointB;
554 
555     for( int i = 0; i < nPointCount / 2; i++ )
556     {
557         getPoint( i, &pointA );
558         const int pos = nPointCount - i - 1;
559         getPoint( pos, &pointB );
560         setPoint( i, &pointB );
561         setPoint( pos, &pointA );
562     }
563 }
564 
565 /************************************************************************/
566 /*                             closeRing()                              */
567 /************************************************************************/
568 
closeRings()569 void OGRLinearRing::closeRings()
570 
571 {
572     if( nPointCount < 2 )
573         return;
574 
575     if( getX(0) != getX(nPointCount-1)
576         || getY(0) != getY(nPointCount-1)
577         || getZ(0) != getZ(nPointCount-1) )
578     {
579         OGRPoint oFirstPoint;
580         getPoint( 0, &oFirstPoint );
581         addPoint( &oFirstPoint );
582     }
583 }
584 
585 /************************************************************************/
586 /*                              isPointInRing()                         */
587 /************************************************************************/
588 
589 /** Returns whether the point is inside the ring.
590  * @param poPoint point
591  * @param bTestEnvelope set to TRUE if the presence of the point inside the
592  *                      ring envelope must be checked first.
593  * @return TRUE or FALSE.
594  */
isPointInRing(const OGRPoint * poPoint,int bTestEnvelope) const595 OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint,
596                                         int bTestEnvelope) const
597 {
598     if( nullptr == poPoint )
599     {
600         CPLDebug( "OGR",
601                   "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - "
602                   "passed point is NULL!" );
603         return FALSE;
604     }
605     if( poPoint->IsEmpty() )
606     {
607         return FALSE;
608     }
609 
610     const int iNumPoints = getNumPoints();
611 
612     // Simple validation
613     if( iNumPoints < 4 )
614         return FALSE;
615 
616     const double dfTestX = poPoint->getX();
617     const double dfTestY = poPoint->getY();
618 
619     // Fast test if point is inside extent of the ring.
620     if( bTestEnvelope )
621     {
622         OGREnvelope extent;
623         getEnvelope(&extent);
624         if( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
625                && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
626         {
627             return FALSE;
628         }
629     }
630 
631     // For every point p in ring,
632     // test if ray starting from given point crosses segment (p - 1, p)
633     int iNumCrossings = 0;
634 
635     double prev_diff_x = getX(0) - dfTestX;
636     double prev_diff_y = getY(0) - dfTestY;
637 
638     for( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
639     {
640         const double x1 = getX(iPoint) - dfTestX;
641         const double y1 = getY(iPoint) - dfTestY;
642 
643         const double x2 = prev_diff_x;
644         const double y2 = prev_diff_y;
645 
646         if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) )
647         {
648             // Check if ray intersects with segment of the ring
649             const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
650             if( 0.0 < dfIntersection )
651             {
652                 // Count intersections
653                 iNumCrossings++;
654             }
655         }
656 
657         prev_diff_x = x1;
658         prev_diff_y = y1;
659     }
660 
661     // If iNumCrossings number is even, given point is outside the ring,
662     // when the crossings number is odd, the point is inside the ring.
663     return iNumCrossings % 2;  // OGRBoolean
664 }
665 
666 /************************************************************************/
667 /*                       isPointOnRingBoundary()                        */
668 /************************************************************************/
669 
670 /** Returns whether the point is on the ring boundary.
671  * @param poPoint point
672  * @param bTestEnvelope set to TRUE if the presence of the point inside the
673  *                      ring envelope must be checked first.
674  * @return TRUE or FALSE.
675  */
isPointOnRingBoundary(const OGRPoint * poPoint,int bTestEnvelope) const676 OGRBoolean OGRLinearRing::isPointOnRingBoundary( const OGRPoint* poPoint,
677                                                  int bTestEnvelope ) const
678 {
679     if( nullptr == poPoint )
680     {
681         CPLDebug( "OGR",
682                   "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* "
683                   "poPoint) - passed point is NULL!" );
684         return 0;
685     }
686 
687     const int iNumPoints = getNumPoints();
688 
689     // Simple validation.
690     if( iNumPoints < 4 )
691         return 0;
692 
693     const double dfTestX = poPoint->getX();
694     const double dfTestY = poPoint->getY();
695 
696     // Fast test if point is inside extent of the ring
697     if( bTestEnvelope )
698     {
699         OGREnvelope extent;
700         getEnvelope(&extent);
701         if( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
702                && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
703         {
704             return 0;
705         }
706     }
707 
708     double prev_diff_x = dfTestX - getX(0);
709     double prev_diff_y = dfTestY - getY(0);
710 
711     for( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
712     {
713         const double dx1 = dfTestX - getX(iPoint);
714         const double dy1 = dfTestY - getY(iPoint);
715 
716         const double dx2 = prev_diff_x;
717         const double dy2 = prev_diff_y;
718 
719         // If the point is on the segment, return immediately.
720         // FIXME? If the test point is not exactly identical to one of
721         // the vertices of the ring, but somewhere on a segment, there's
722         // little chance that we get 0. So that should be tested against some
723         // epsilon.
724 
725         if( dx1 * dy2 - dx2 * dy1 == 0 )
726         {
727             // If iPoint and iPointPrev are the same, go on.
728             if( !(dx1 == dx2 && dy1 == dy2) )
729             {
730                 const double dx_segment = getX(iPoint) - getX(iPoint-1);
731                 const double dy_segment = getY(iPoint) - getY(iPoint-1);
732                 const double crossproduct =
733                     dx2 * dx_segment + dy2 * dy_segment;
734                 if( crossproduct >= 0 )
735                 {
736                     const double sq_length_seg = dx_segment * dx_segment +
737                                                  dy_segment * dy_segment;
738                     if( crossproduct <= sq_length_seg )
739                     {
740                         return 1;
741                     }
742                 }
743             }
744         }
745 
746         prev_diff_x = dx1;
747         prev_diff_y = dy1;
748     }
749 
750     return 0;
751 }
752 
753 /************************************************************************/
754 /*                             transform()                              */
755 /************************************************************************/
756 
transform(OGRCoordinateTransformation * poCT)757 OGRErr OGRLinearRing::transform( OGRCoordinateTransformation *poCT )
758 
759 {
760     const bool bIsClosed = getNumPoints() > 2 && CPL_TO_BOOL(get_IsClosed());
761     OGRErr eErr = OGRLineString::transform(poCT);
762     if( bIsClosed && eErr == OGRERR_NONE && !get_IsClosed() )
763     {
764         CPLDebug("OGR", "Linearring is not closed after coordinate "
765                   "transformation. Forcing last point to be identical to "
766                   "first one");
767         // Force last point to be identical to first point.
768         // This is a safety belt in case the reprojection of the same coordinate
769         // isn't perfectly stable. This can for example happen in very rare cases
770         // when reprojecting a cutline with a RPC transform with a DEM that
771         // is a VRT whose sources are resampled...
772         OGRPoint oStartPoint;
773         StartPoint( &oStartPoint );
774 
775         setPoint( getNumPoints()-1, &oStartPoint);
776     }
777     return eErr;
778 }
779 
780 /************************************************************************/
781 /*                          CastToLineString()                          */
782 /************************************************************************/
783 
784 /**
785  * \brief Cast to line string.
786  *
787  * The passed in geometry is consumed and a new one returned .
788  *
789  * @param poLR the input geometry - ownership is passed to the method.
790  * @return new geometry.
791  */
792 
CastToLineString(OGRLinearRing * poLR)793 OGRLineString* OGRLinearRing::CastToLineString( OGRLinearRing* poLR )
794 {
795     return TransferMembersAndDestroy(poLR, new OGRLineString());
796 }
797 
798 //! @cond Doxygen_Suppress
799 /************************************************************************/
800 /*                     GetCasterToLineString()                          */
801 /************************************************************************/
802 
CasterToLineString(OGRCurve * poCurve)803 OGRLineString* OGRLinearRing::CasterToLineString( OGRCurve* poCurve )
804 {
805     return OGRLinearRing::CastToLineString(poCurve->toLinearRing());
806 }
807 
GetCasterToLineString() const808 OGRCurveCasterToLineString OGRLinearRing::GetCasterToLineString() const
809 {
810     return OGRLinearRing::CasterToLineString;
811 }
812 
813 /************************************************************************/
814 /*                        GetCasterToLinearRing()                       */
815 /************************************************************************/
816 
CasterToLinearRing(OGRCurve * poCurve)817 static OGRLinearRing* CasterToLinearRing(OGRCurve* poCurve)
818 {
819     return poCurve->toLinearRing();
820 }
821 
GetCasterToLinearRing() const822 OGRCurveCasterToLinearRing OGRLinearRing::GetCasterToLinearRing() const
823 {
824     return ::CasterToLinearRing;
825 }
826 //! @endcond
827