1 /******************************************************************************
2 *
3 * Project: OpenGIS Simple Features Reference Implementation
4 * Purpose: The Point geometry class.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 1999, Frank Warmerdam
9 * Copyright (c) 2008-2011, 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 <cmath>
34 #include <cstdio>
35 #include <cstring>
36 #include <algorithm>
37 #include <limits>
38 #include <new>
39
40 #include "cpl_conv.h"
41 #include "ogr_core.h"
42 #include "ogr_p.h"
43 #include "ogr_spatialref.h"
44
45 CPL_CVSID("$Id: ogrpoint.cpp 3798cbe48457b7127606931896549f26507469db 2021-04-09 15:04:16 +0200 Even Rouault $")
46
47 /************************************************************************/
48 /* GetEmptyNonEmptyFlag() */
49 /************************************************************************/
50
GetEmptyNonEmptyFlag(double x,double y)51 static int GetEmptyNonEmptyFlag(double x, double y)
52 {
53 if( std::isnan(x) || std::isnan(y) )
54 return 0;
55 return OGRGeometry::OGR_G_NOT_EMPTY_POINT;
56 }
57
58 /************************************************************************/
59 /* OGRPoint() */
60 /************************************************************************/
61
62 /**
63 * \brief Create an empty point.
64 */
65
OGRPoint()66 OGRPoint::OGRPoint(): x(0.0), y(0.0), z(0.0), m(0.0)
67
68 {
69 flags = 0;
70 }
71
72 /************************************************************************/
73 /* OGRPoint() */
74 /************************************************************************/
75
76 /**
77 * \brief Create a point.
78 * @param xIn x
79 * @param yIn y
80 * @param zIn z
81 */
82
OGRPoint(double xIn,double yIn,double zIn)83 OGRPoint::OGRPoint( double xIn, double yIn, double zIn ) :
84 x(xIn),
85 y(yIn),
86 z(zIn),
87 m(0.0)
88 {
89 flags = GetEmptyNonEmptyFlag(xIn, yIn) | OGR_G_3D;
90 }
91
92 /************************************************************************/
93 /* OGRPoint() */
94 /************************************************************************/
95
96 /**
97 * \brief Create a point.
98 * @param xIn x
99 * @param yIn y
100 */
101
OGRPoint(double xIn,double yIn)102 OGRPoint::OGRPoint( double xIn, double yIn ) :
103 x(xIn),
104 y(yIn),
105 z(0.0),
106 m(0.0)
107 {
108 flags = GetEmptyNonEmptyFlag(xIn, yIn);
109 }
110
111 /************************************************************************/
112 /* OGRPoint() */
113 /************************************************************************/
114
115 /**
116 * \brief Create a point.
117 * @param xIn x
118 * @param yIn y
119 * @param zIn z
120 * @param mIn m
121 */
122
OGRPoint(double xIn,double yIn,double zIn,double mIn)123 OGRPoint::OGRPoint( double xIn, double yIn, double zIn, double mIn ) :
124 x(xIn),
125 y(yIn),
126 z(zIn),
127 m(mIn)
128 {
129 flags = GetEmptyNonEmptyFlag(xIn, yIn) | OGR_G_3D | OGR_G_MEASURED;
130 }
131
132 /************************************************************************/
133 /* createXYM() */
134 /************************************************************************/
135
136 /**
137 * \brief Create a XYM point.
138 * @param x x
139 * @param y y
140 * @param m m
141 * @since GDAL 3.1
142 */
143
createXYM(double x,double y,double m)144 OGRPoint* OGRPoint::createXYM( double x, double y, double m )
145 {
146 auto p = new OGRPoint(x, y, 0, m);
147 p->flags &= ~OGR_G_3D;
148 return p;
149 }
150
151 /************************************************************************/
152 /* OGRPoint( const OGRPoint& ) */
153 /************************************************************************/
154
155 /**
156 * \brief Copy constructor.
157 *
158 * Note: before GDAL 2.1, only the default implementation of the constructor
159 * existed, which could be unsafe to use.
160 *
161 * @since GDAL 2.1
162 */
163
164 OGRPoint::OGRPoint( const OGRPoint& ) = default;
165
166 /************************************************************************/
167 /* ~OGRPoint() */
168 /************************************************************************/
169
170 OGRPoint::~OGRPoint() = default;
171
172 /************************************************************************/
173 /* operator=( const OGRPoint& ) */
174 /************************************************************************/
175
176 /**
177 * \brief Assignment operator.
178 *
179 * Note: before GDAL 2.1, only the default implementation of the operator
180 * existed, which could be unsafe to use.
181 *
182 * @since GDAL 2.1
183 */
184
operator =(const OGRPoint & other)185 OGRPoint& OGRPoint::operator=( const OGRPoint& other )
186 {
187 if( this != &other)
188 {
189 OGRGeometry::operator=( other );
190
191 x = other.x;
192 y = other.y;
193 z = other.z;
194 m = other.m;
195 }
196 return *this;
197 }
198
199 /************************************************************************/
200 /* clone() */
201 /* */
202 /* Make a new object that is a copy of this object. */
203 /************************************************************************/
204
clone() const205 OGRPoint *OGRPoint::clone() const
206
207 {
208 return new (std::nothrow) OGRPoint(*this);
209 }
210
211 /************************************************************************/
212 /* empty() */
213 /************************************************************************/
empty()214 void OGRPoint::empty()
215
216 {
217 x = 0.0;
218 y = 0.0;
219 z = 0.0;
220 m = 0.0;
221 flags &= ~OGR_G_NOT_EMPTY_POINT;
222 }
223
224 /************************************************************************/
225 /* getDimension() */
226 /************************************************************************/
227
getDimension() const228 int OGRPoint::getDimension() const
229
230 {
231 return 0;
232 }
233
234 /************************************************************************/
235 /* getGeometryType() */
236 /************************************************************************/
237
getGeometryType() const238 OGRwkbGeometryType OGRPoint::getGeometryType() const
239
240 {
241 if( (flags & OGR_G_3D) && (flags & OGR_G_MEASURED) )
242 return wkbPointZM;
243 else if( flags & OGR_G_MEASURED )
244 return wkbPointM;
245 else if( flags & OGR_G_3D )
246 return wkbPoint25D;
247 else
248 return wkbPoint;
249 }
250
251 /************************************************************************/
252 /* getGeometryName() */
253 /************************************************************************/
254
getGeometryName() const255 const char * OGRPoint::getGeometryName() const
256
257 {
258 return "POINT";
259 }
260
261 /************************************************************************/
262 /* flattenTo2D() */
263 /************************************************************************/
264
flattenTo2D()265 void OGRPoint::flattenTo2D()
266
267 {
268 z = 0.0;
269 m = 0.0;
270 flags &= ~OGR_G_3D;
271 setMeasured(FALSE);
272 }
273
274 /************************************************************************/
275 /* setCoordinateDimension() */
276 /************************************************************************/
277
setCoordinateDimension(int nNewDimension)278 void OGRPoint::setCoordinateDimension( int nNewDimension )
279
280 {
281 if( nNewDimension == 2 )
282 flattenTo2D();
283 else if( nNewDimension == 3 )
284 flags |= OGR_G_3D;
285
286 setMeasured(FALSE);
287 }
288
289 /************************************************************************/
290 /* WkbSize() */
291 /* */
292 /* Return the size of this object in well known binary */
293 /* representation including the byte order, and type information. */
294 /************************************************************************/
295
WkbSize() const296 size_t OGRPoint::WkbSize() const
297
298 {
299 if( (flags & OGR_G_3D) && (flags & OGR_G_MEASURED) )
300 return 37;
301 else if( (flags & OGR_G_3D) || (flags & OGR_G_MEASURED) )
302 return 29;
303 else
304 return 21;
305 }
306
307 /************************************************************************/
308 /* importFromWkb() */
309 /* */
310 /* Initialize from serialized stream in well known binary */
311 /* format. */
312 /************************************************************************/
313
importFromWkb(const unsigned char * pabyData,size_t nSize,OGRwkbVariant eWkbVariant,size_t & nBytesConsumedOut)314 OGRErr OGRPoint::importFromWkb( const unsigned char *pabyData,
315 size_t nSize,
316 OGRwkbVariant eWkbVariant,
317 size_t& nBytesConsumedOut )
318
319 {
320 nBytesConsumedOut = 0;
321 OGRwkbByteOrder eByteOrder = wkbNDR;
322
323 flags = 0;
324 OGRErr eErr =
325 importPreambleFromWkb( pabyData, nSize, eByteOrder, eWkbVariant );
326 pabyData += 5;
327 if( eErr != OGRERR_NONE )
328 return eErr;
329
330 if( nSize != static_cast<size_t>(-1) )
331 {
332 if( (nSize < 37) && ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED)) )
333 return OGRERR_NOT_ENOUGH_DATA;
334 else if( (nSize < 29) && ((flags & OGR_G_3D) ||
335 (flags & OGR_G_MEASURED)) )
336 return OGRERR_NOT_ENOUGH_DATA;
337 else if( nSize < 21 )
338 return OGRERR_NOT_ENOUGH_DATA;
339 }
340
341 nBytesConsumedOut = 5 + 8 * (2 + ((flags & OGR_G_3D) ? 1 : 0)+
342 ((flags & OGR_G_MEASURED) ? 1 : 0));
343
344 /* -------------------------------------------------------------------- */
345 /* Get the vertex. */
346 /* -------------------------------------------------------------------- */
347 memcpy( &x, pabyData, 8 );
348 pabyData += 8;
349 memcpy( &y, pabyData, 8 );
350 pabyData += 8;
351
352 if( OGR_SWAP( eByteOrder ) )
353 {
354 CPL_SWAPDOUBLE( &x );
355 CPL_SWAPDOUBLE( &y );
356 }
357
358 if( flags & OGR_G_3D )
359 {
360 memcpy( &z, pabyData, 8 );
361 pabyData += 8;
362 if( OGR_SWAP( eByteOrder ) )
363 CPL_SWAPDOUBLE( &z );
364 }
365 else
366 {
367 z = 0;
368 }
369 if( flags & OGR_G_MEASURED )
370 {
371 memcpy( &m, pabyData, 8 );
372 /*pabyData += 8; */
373 if( OGR_SWAP( eByteOrder ) )
374 {
375 CPL_SWAPDOUBLE( &m );
376 }
377 }
378 else
379 {
380 m = 0;
381 }
382
383 // Detect coordinates are not NaN --> NOT EMPTY.
384 if( !(CPLIsNan(x) && CPLIsNan(y)) )
385 flags |= OGR_G_NOT_EMPTY_POINT;
386
387 return OGRERR_NONE;
388 }
389
390 /************************************************************************/
391 /* exportToWkb() */
392 /* */
393 /* Build a well known binary representation of this object. */
394 /************************************************************************/
395
exportToWkb(OGRwkbByteOrder eByteOrder,unsigned char * pabyData,OGRwkbVariant eWkbVariant) const396 OGRErr OGRPoint::exportToWkb( OGRwkbByteOrder eByteOrder,
397 unsigned char * pabyData,
398 OGRwkbVariant eWkbVariant ) const
399
400 {
401 /* -------------------------------------------------------------------- */
402 /* Set the byte order. */
403 /* -------------------------------------------------------------------- */
404 pabyData[0] =
405 DB2_V72_UNFIX_BYTE_ORDER(static_cast<unsigned char>(eByteOrder));
406 pabyData += 1;
407
408 /* -------------------------------------------------------------------- */
409 /* Set the geometry feature type. */
410 /* -------------------------------------------------------------------- */
411
412 GUInt32 nGType = getGeometryType();
413
414 if( eWkbVariant == wkbVariantPostGIS1 )
415 {
416 nGType = wkbFlatten(nGType);
417 if( Is3D() )
418 // Explicitly set wkb25DBit.
419 nGType =
420 static_cast<OGRwkbGeometryType>(nGType | wkb25DBitInternalUse);
421 if( IsMeasured() )
422 nGType = static_cast<OGRwkbGeometryType>(nGType | 0x40000000);
423 }
424 else if( eWkbVariant == wkbVariantIso )
425 {
426 nGType = getIsoGeometryType();
427 }
428
429 if( eByteOrder == wkbNDR )
430 {
431 CPL_LSBPTR32( &nGType );
432 }
433 else
434 {
435 CPL_MSBPTR32( &nGType );
436 }
437
438 memcpy( pabyData, &nGType, 4 );
439 pabyData += 4;
440
441 /* -------------------------------------------------------------------- */
442 /* Copy in the raw data. Swap if needed. */
443 /* -------------------------------------------------------------------- */
444
445 if( IsEmpty() && eWkbVariant == wkbVariantIso )
446 {
447 const double dNan = std::numeric_limits<double>::quiet_NaN();
448 memcpy( pabyData, &dNan, 8 );
449 if( OGR_SWAP( eByteOrder ) )
450 CPL_SWAPDOUBLE( pabyData );
451 pabyData += 8;
452 memcpy( pabyData, &dNan, 8 );
453 if( OGR_SWAP( eByteOrder ) )
454 CPL_SWAPDOUBLE( pabyData );
455 pabyData += 8;
456 if( flags & OGR_G_3D ) {
457 memcpy( pabyData, &dNan, 8 );
458 if( OGR_SWAP( eByteOrder ) )
459 CPL_SWAPDOUBLE( pabyData );
460 pabyData += 8;
461 }
462 if( flags & OGR_G_MEASURED ) {
463 memcpy( pabyData, &dNan, 8 );
464 if( OGR_SWAP( eByteOrder ) )
465 CPL_SWAPDOUBLE( pabyData );
466 }
467 }
468 else
469 {
470 memcpy( pabyData, &x, 8 );
471 if( OGR_SWAP( eByteOrder ) )
472 CPL_SWAPDOUBLE( pabyData );
473 pabyData += 8;
474 memcpy( pabyData, &y, 8 );
475 if( OGR_SWAP( eByteOrder ) )
476 CPL_SWAPDOUBLE( pabyData );
477 pabyData += 8;
478 if( flags & OGR_G_3D ) {
479 memcpy( pabyData, &z, 8 );
480 if( OGR_SWAP( eByteOrder ) )
481 CPL_SWAPDOUBLE( pabyData );
482 pabyData += 8;
483 }
484 if( flags & OGR_G_MEASURED )
485 {
486 memcpy( pabyData, &m, 8 );
487 if( OGR_SWAP( eByteOrder ) )
488 CPL_SWAPDOUBLE( pabyData );
489 }
490 }
491
492 return OGRERR_NONE;
493 }
494
495 /************************************************************************/
496 /* importFromWkt() */
497 /* */
498 /* Instantiate point from well known text format ``POINT */
499 /* (x,y)''. */
500 /************************************************************************/
501
importFromWkt(const char ** ppszInput)502 OGRErr OGRPoint::importFromWkt( const char ** ppszInput )
503
504 {
505 int bHasZ = FALSE;
506 int bHasM = FALSE;
507 bool bIsEmpty = false;
508 OGRErr eErr = importPreambleFromWkt(ppszInput, &bHasZ, &bHasM, &bIsEmpty);
509 flags = 0;
510 if( eErr != OGRERR_NONE )
511 return eErr;
512 if( bHasZ ) flags |= OGR_G_3D;
513 if( bHasM ) flags |= OGR_G_MEASURED;
514 if( bIsEmpty )
515 {
516 return OGRERR_NONE;
517 }
518 else
519 {
520 flags |= OGR_G_NOT_EMPTY_POINT;
521 }
522
523 const char *pszInput = *ppszInput;
524
525 /* -------------------------------------------------------------------- */
526 /* Read the point list which should consist of exactly one point. */
527 /* -------------------------------------------------------------------- */
528 OGRRawPoint *poPoints = nullptr;
529 double *padfZ = nullptr;
530 double *padfM = nullptr;
531 int nMaxPoint = 0;
532 int nPoints = 0;
533 int flagsFromInput = flags;
534
535 pszInput = OGRWktReadPointsM( pszInput, &poPoints, &padfZ, &padfM,
536 &flagsFromInput,
537 &nMaxPoint, &nPoints );
538 if( pszInput == nullptr || nPoints != 1 )
539 {
540 CPLFree( poPoints );
541 CPLFree( padfZ );
542 CPLFree( padfM );
543 return OGRERR_CORRUPT_DATA;
544 }
545 if( (flagsFromInput & OGR_G_3D) && !(flags & OGR_G_3D) )
546 {
547 flags |= OGR_G_3D;
548 bHasZ = TRUE;
549 }
550 if( (flagsFromInput & OGR_G_MEASURED) && !(flags & OGR_G_MEASURED) )
551 {
552 flags |= OGR_G_MEASURED;
553 bHasM = TRUE;
554 }
555
556 x = poPoints[0].x;
557 y = poPoints[0].y;
558
559 CPLFree( poPoints );
560
561 if( bHasZ )
562 {
563 if( padfZ != nullptr )
564 z = padfZ[0];
565 }
566 if( bHasM )
567 {
568 if( padfM != nullptr )
569 m = padfM[0];
570 }
571
572 CPLFree( padfZ );
573 CPLFree( padfM );
574
575 *ppszInput = pszInput;
576
577 return OGRERR_NONE;
578 }
579
580 /************************************************************************/
581 /* exportToWkt() */
582 /* */
583 /* Translate this structure into its well known text format */
584 /* equivalent. */
585 /************************************************************************/
586
exportToWkt(const OGRWktOptions & opts,OGRErr * err) const587 std::string OGRPoint::exportToWkt(const OGRWktOptions& opts, OGRErr *err) const
588 {
589 std::string wkt = getGeometryName() + wktTypeString(opts.variant);
590 if( IsEmpty() )
591 {
592 wkt += "EMPTY";
593 }
594 else
595 {
596 wkt += "(";
597
598 bool measured = ((opts.variant == wkbVariantIso) && IsMeasured());
599 wkt += OGRMakeWktCoordinateM(x, y, z, m, Is3D(), measured, opts);
600
601 wkt += ")";
602 }
603
604 if (err)
605 *err = OGRERR_NONE;
606 return wkt;
607 }
608
609 /************************************************************************/
610 /* getEnvelope() */
611 /************************************************************************/
612
getEnvelope(OGREnvelope * psEnvelope) const613 void OGRPoint::getEnvelope( OGREnvelope * psEnvelope ) const
614
615 {
616 psEnvelope->MinX = getX();
617 psEnvelope->MaxX = getX();
618 psEnvelope->MinY = getY();
619 psEnvelope->MaxY = getY();
620 }
621
622 /************************************************************************/
623 /* getEnvelope() */
624 /************************************************************************/
625
getEnvelope(OGREnvelope3D * psEnvelope) const626 void OGRPoint::getEnvelope( OGREnvelope3D * psEnvelope ) const
627
628 {
629 psEnvelope->MinX = getX();
630 psEnvelope->MaxX = getX();
631 psEnvelope->MinY = getY();
632 psEnvelope->MaxY = getY();
633 psEnvelope->MinZ = getZ();
634 psEnvelope->MaxZ = getZ();
635 }
636
637 /**
638 * \fn double OGRPoint::getX() const;
639 *
640 * \brief Fetch X coordinate.
641 *
642 * Relates to the SFCOM IPoint::get_X() method.
643 *
644 * @return the X coordinate of this point.
645 */
646
647 /**
648 * \fn double OGRPoint::getY() const;
649 *
650 * \brief Fetch Y coordinate.
651 *
652 * Relates to the SFCOM IPoint::get_Y() method.
653 *
654 * @return the Y coordinate of this point.
655 */
656
657 /**
658 * \fn double OGRPoint::getZ() const;
659 *
660 * \brief Fetch Z coordinate.
661 *
662 * Relates to the SFCOM IPoint::get_Z() method.
663 *
664 * @return the Z coordinate of this point, or zero if it is a 2D point.
665 */
666
667 /**
668 * \fn void OGRPoint::setX( double xIn );
669 *
670 * \brief Assign point X coordinate.
671 *
672 * There is no corresponding SFCOM method.
673 */
674
675 /**
676 * \fn void OGRPoint::setY( double yIn );
677 *
678 * \brief Assign point Y coordinate.
679 *
680 * There is no corresponding SFCOM method.
681 */
682
683 /**
684 * \fn void OGRPoint::setZ( double zIn );
685 *
686 * \brief Assign point Z coordinate.
687 * Calling this method will force the geometry
688 * coordinate dimension to 3D (wkbPoint|wkbZ).
689 *
690 * There is no corresponding SFCOM method.
691 */
692
693 /************************************************************************/
694 /* Equal() */
695 /************************************************************************/
696
Equals(const OGRGeometry * poOther) const697 OGRBoolean OGRPoint::Equals( const OGRGeometry * poOther ) const
698
699 {
700 if( poOther== this )
701 return TRUE;
702
703 if( poOther->getGeometryType() != getGeometryType() )
704 return FALSE;
705
706 const auto poOPoint = poOther->toPoint();
707 if( flags != poOPoint->flags )
708 return FALSE;
709
710 if( IsEmpty() )
711 return TRUE;
712
713 // Should eventually test the SRS.
714 if( poOPoint->getX() != getX()
715 || poOPoint->getY() != getY()
716 || poOPoint->getZ() != getZ() )
717 return FALSE;
718
719 return TRUE;
720 }
721
722 /************************************************************************/
723 /* transform() */
724 /************************************************************************/
725
transform(OGRCoordinateTransformation * poCT)726 OGRErr OGRPoint::transform( OGRCoordinateTransformation *poCT )
727
728 {
729 if( poCT->Transform( 1, &x, &y, &z ) )
730 {
731 assignSpatialReference( poCT->GetTargetCS() );
732 return OGRERR_NONE;
733 }
734
735 return OGRERR_FAILURE;
736 }
737
738 /************************************************************************/
739 /* swapXY() */
740 /************************************************************************/
741
swapXY()742 void OGRPoint::swapXY()
743 {
744 std::swap(x, y);
745 }
746
747 /************************************************************************/
748 /* Within() */
749 /************************************************************************/
750
Within(const OGRGeometry * poOtherGeom) const751 OGRBoolean OGRPoint::Within( const OGRGeometry *poOtherGeom ) const
752
753 {
754 if( !IsEmpty() && poOtherGeom != nullptr &&
755 wkbFlatten(poOtherGeom->getGeometryType()) == wkbCurvePolygon )
756 {
757 const auto poCurve = poOtherGeom->toCurvePolygon();
758 return poCurve->Contains(this);
759 }
760
761 return OGRGeometry::Within(poOtherGeom);
762 }
763
764 /************************************************************************/
765 /* Intersects() */
766 /************************************************************************/
767
Intersects(const OGRGeometry * poOtherGeom) const768 OGRBoolean OGRPoint::Intersects( const OGRGeometry *poOtherGeom ) const
769
770 {
771 if( !IsEmpty() && poOtherGeom != nullptr &&
772 wkbFlatten(poOtherGeom->getGeometryType()) == wkbCurvePolygon )
773 {
774 const auto poCurve = poOtherGeom->toCurvePolygon();
775 return poCurve->Intersects(this);
776 }
777
778 return OGRGeometry::Intersects(poOtherGeom);
779 }
780