1 : /******************************************************************************
2 : * $Id: ogrlinearring.cpp 24284 2012-04-21 17:31:57Z rouault $
3 : *
4 : * Project: OpenGIS Simple Features Reference Implementation
5 : * Purpose: The OGRLinearRing geometry class.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
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 "ogr_geometry.h"
31 : #include "ogr_p.h"
32 :
33 : CPL_CVSID("$Id: ogrlinearring.cpp 24284 2012-04-21 17:31:57Z rouault $");
34 :
35 : /************************************************************************/
36 : /* OGRLinearRing() */
37 : /************************************************************************/
38 :
39 127881 : OGRLinearRing::OGRLinearRing()
40 :
41 : {
42 127881 : }
43 :
44 : /************************************************************************/
45 : /* ~OGRLinearRing() */
46 : /************************************************************************/
47 258051 : OGRLinearRing::~OGRLinearRing()
48 :
49 : {
50 258051 : }
51 :
52 : /************************************************************************/
53 : /* OGRLinearRing() */
54 : /************************************************************************/
55 :
56 130170 : OGRLinearRing::OGRLinearRing( OGRLinearRing * poSrcRing )
57 :
58 : {
59 130170 : if( poSrcRing == NULL )
60 : {
61 0 : CPLDebug( "OGR", "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - passed in ring is NULL!" );
62 0 : return;
63 : }
64 :
65 130170 : setNumPoints( poSrcRing->getNumPoints() );
66 :
67 : memcpy( paoPoints, poSrcRing->paoPoints,
68 130170 : sizeof(OGRRawPoint) * getNumPoints() );
69 :
70 130170 : if( poSrcRing->padfZ )
71 : {
72 360 : Make3D();
73 360 : memcpy( padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints() );
74 : }
75 0 : }
76 :
77 : /************************************************************************/
78 : /* getGeometryName() */
79 : /************************************************************************/
80 :
81 4090 : const char * OGRLinearRing::getGeometryName() const
82 :
83 : {
84 4090 : return "LINEARRING";
85 : }
86 :
87 : /************************************************************************/
88 : /* WkbSize() */
89 : /* */
90 : /* Disable this method. */
91 : /************************************************************************/
92 :
93 0 : int OGRLinearRing::WkbSize() const
94 :
95 : {
96 0 : return 0;
97 : }
98 :
99 : /************************************************************************/
100 : /* importFromWkb() */
101 : /* */
102 : /* Disable method for this class. */
103 : /************************************************************************/
104 :
105 0 : OGRErr OGRLinearRing::importFromWkb( unsigned char *pabyData, int nSize )
106 :
107 : {
108 : (void) pabyData;
109 : (void) nSize;
110 :
111 0 : return OGRERR_UNSUPPORTED_OPERATION;
112 : }
113 :
114 : /************************************************************************/
115 : /* exportToWkb() */
116 : /* */
117 : /* Disable method for this class. */
118 : /************************************************************************/
119 :
120 0 : OGRErr OGRLinearRing::exportToWkb( OGRwkbByteOrder eByteOrder,
121 : unsigned char * pabyData ) const
122 :
123 : {
124 : (void) eByteOrder;
125 : (void) pabyData;
126 :
127 0 : return OGRERR_UNSUPPORTED_OPERATION;
128 : }
129 :
130 : /************************************************************************/
131 : /* _importFromWkb() */
132 : /* */
133 : /* Helper method for OGRPolygon. NOT A NORMAL importFromWkb() */
134 : /* method! */
135 : /************************************************************************/
136 :
137 2208 : OGRErr OGRLinearRing::_importFromWkb( OGRwkbByteOrder eByteOrder, int b3D,
138 : unsigned char * pabyData,
139 : int nBytesAvailable )
140 :
141 : {
142 2208 : if( nBytesAvailable < 4 && nBytesAvailable != -1 )
143 0 : return OGRERR_NOT_ENOUGH_DATA;
144 :
145 : /* -------------------------------------------------------------------- */
146 : /* Get the vertex count. */
147 : /* -------------------------------------------------------------------- */
148 : int nNewNumPoints;
149 :
150 2208 : memcpy( &nNewNumPoints, pabyData, 4 );
151 :
152 2208 : if( OGR_SWAP( eByteOrder ) )
153 38 : nNewNumPoints = CPL_SWAP32(nNewNumPoints);
154 :
155 : /* Check if the wkb stream buffer is big enough to store
156 : * fetched number of points.
157 : * 16 or 24 - size of point structure
158 : */
159 2208 : int nPointSize = (b3D ? 24 : 16);
160 2208 : if (nNewNumPoints < 0 || nNewNumPoints > INT_MAX / nPointSize)
161 0 : return OGRERR_CORRUPT_DATA;
162 2208 : int nBufferMinSize = nPointSize * nNewNumPoints;
163 :
164 2208 : if( nBytesAvailable != -1 && nBufferMinSize > nBytesAvailable - 4 )
165 : {
166 : CPLError( CE_Failure, CPLE_AppDefined,
167 0 : "Length of input WKB is too small" );
168 0 : return OGRERR_NOT_ENOUGH_DATA;
169 : }
170 :
171 : /* (Re)Allocation of paoPoints buffer. */
172 2208 : setNumPoints( nNewNumPoints );
173 :
174 2208 : if( b3D )
175 546 : Make3D();
176 : else
177 1662 : Make2D();
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Get the vertices */
181 : /* -------------------------------------------------------------------- */
182 2208 : int i = 0;
183 :
184 2208 : if( !b3D )
185 : {
186 1662 : memcpy( paoPoints, pabyData + 4, 16 * nPointCount );
187 : }
188 : else
189 : {
190 10850 : for( int i = 0; i < nPointCount; i++ )
191 : {
192 10304 : memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
193 10304 : memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
194 10304 : memcpy( padfZ + i, pabyData + 4 + 24 * i + 16, 8 );
195 : }
196 : }
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Byte swap if needed. */
200 : /* -------------------------------------------------------------------- */
201 2208 : if( OGR_SWAP( eByteOrder ) )
202 : {
203 196 : for( i = 0; i < nPointCount; i++ )
204 : {
205 158 : CPL_SWAPDOUBLE( &(paoPoints[i].x) );
206 158 : CPL_SWAPDOUBLE( &(paoPoints[i].y) );
207 :
208 158 : if( b3D )
209 : {
210 26 : CPL_SWAPDOUBLE( padfZ + i );
211 : }
212 : }
213 : }
214 :
215 2208 : return OGRERR_NONE;
216 : }
217 :
218 : /************************************************************************/
219 : /* _exportToWkb() */
220 : /* */
221 : /* Helper method for OGRPolygon. THIS IS NOT THE NORMAL */
222 : /* exportToWkb() METHOD! */
223 : /************************************************************************/
224 :
225 4808 : OGRErr OGRLinearRing::_exportToWkb( OGRwkbByteOrder eByteOrder, int b3D,
226 : unsigned char * pabyData ) const
227 :
228 : {
229 : int i, nWords;
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Copy in the raw data. */
233 : /* -------------------------------------------------------------------- */
234 4808 : memcpy( pabyData, &nPointCount, 4 );
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Copy in the raw data. */
238 : /* -------------------------------------------------------------------- */
239 4808 : if( b3D )
240 : {
241 148 : nWords = 3 * nPointCount;
242 3106 : for( i = 0; i < nPointCount; i++ )
243 : {
244 2958 : memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
245 2958 : memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
246 2958 : if( padfZ == NULL )
247 0 : memset( pabyData+4+i*24+16, 0, 8 );
248 : else
249 2958 : memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
250 : }
251 : }
252 : else
253 : {
254 4660 : nWords = 2 * nPointCount;
255 4660 : memcpy( pabyData+4, paoPoints, 16 * nPointCount );
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Swap if needed. */
260 : /* -------------------------------------------------------------------- */
261 4808 : if( OGR_SWAP( eByteOrder ) )
262 : {
263 : int nCount;
264 :
265 32 : nCount = CPL_SWAP32( nPointCount );
266 32 : memcpy( pabyData, &nCount, 4 );
267 :
268 322 : for( i = 0; i < nWords; i++ )
269 : {
270 290 : CPL_SWAPDOUBLE( pabyData + 4 + 8 * i );
271 : }
272 : }
273 :
274 4808 : return OGRERR_NONE;
275 : }
276 :
277 : /************************************************************************/
278 : /* _WkbSize() */
279 : /* */
280 : /* Helper method for OGRPolygon. NOT THE NORMAL WkbSize() METHOD! */
281 : /************************************************************************/
282 :
283 14760 : int OGRLinearRing::_WkbSize( int b3D ) const
284 :
285 : {
286 14760 : if( b3D )
287 1532 : return 4 + 24 * nPointCount;
288 : else
289 13228 : return 4 + 16 * nPointCount;
290 : }
291 :
292 : /************************************************************************/
293 : /* clone() */
294 : /* */
295 : /* We override the OGRCurve clone() to ensure that we get the */
296 : /* correct virtual table. */
297 : /************************************************************************/
298 :
299 746 : OGRGeometry *OGRLinearRing::clone() const
300 :
301 : {
302 : OGRLinearRing *poNewLinearRing;
303 :
304 746 : poNewLinearRing = new OGRLinearRing();
305 746 : poNewLinearRing->assignSpatialReference( getSpatialReference() );
306 :
307 746 : poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ );
308 :
309 746 : return poNewLinearRing;
310 : }
311 :
312 :
313 : /************************************************************************/
314 : /* epsilonEqual() */
315 : /************************************************************************/
316 :
317 : static const double EPSILON = 1E-5;
318 :
319 1788 : static inline bool epsilonEqual(double a, double b, double eps)
320 : {
321 1788 : return (::fabs(a - b) < eps);
322 : }
323 :
324 : /************************************************************************/
325 : /* isClockwise() */
326 : /************************************************************************/
327 :
328 : /**
329 : * \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
330 : *
331 : * @return TRUE if clockwise otherwise FALSE.
332 : */
333 :
334 740 : int OGRLinearRing::isClockwise() const
335 :
336 : {
337 : int i, v, next;
338 : double dx0, dy0, dx1, dy1, crossproduct;
339 740 : int bUseFallback = FALSE;
340 :
341 740 : if( nPointCount < 2 )
342 0 : return TRUE;
343 :
344 : /* Find the lowest rightmost vertex */
345 740 : v = 0;
346 8924 : for ( i = 1; i < nPointCount - 1; i++ )
347 : {
348 : /* => v < end */
349 22020 : if ( paoPoints[i].y< paoPoints[v].y ||
350 13228 : ( paoPoints[i].y== paoPoints[v].y &&
351 608 : paoPoints[i].x > paoPoints[v].x ) )
352 : {
353 1798 : v = i;
354 : }
355 : }
356 :
357 : /* previous */
358 740 : next = v - 1;
359 740 : if ( next < 0 )
360 : {
361 136 : next = nPointCount - 1 - 1;
362 : }
363 :
364 1252 : if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
365 512 : epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
366 : {
367 : /* Don't try to be too clever by retrying with a next point */
368 : /* This can lead to false results as in the case of #3356 */
369 4 : bUseFallback = TRUE;
370 : }
371 :
372 740 : dx0 = paoPoints[next].x - paoPoints[v].x;
373 740 : dy0 = paoPoints[next].y - paoPoints[v].y;
374 :
375 :
376 : /* following */
377 740 : next = v + 1;
378 740 : if ( next >= nPointCount - 1 )
379 : {
380 202 : next = 0;
381 : }
382 :
383 844 : if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
384 104 : epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
385 : {
386 : /* Don't try to be too clever by retrying with a next point */
387 : /* This can lead to false results as in the case of #3356 */
388 6 : bUseFallback = TRUE;
389 : }
390 :
391 740 : dx1 = paoPoints[next].x - paoPoints[v].x;
392 740 : dy1 = paoPoints[next].y - paoPoints[v].y;
393 :
394 740 : crossproduct = dx1 * dy0 - dx0 * dy1;
395 :
396 740 : if (!bUseFallback)
397 : {
398 734 : if ( crossproduct > 0 ) /* CCW */
399 152 : return FALSE;
400 582 : else if ( crossproduct < 0 ) /* CW */
401 582 : return TRUE;
402 : }
403 :
404 : /* ok, this is a degenerate case : the extent of the polygon is less than EPSILON */
405 : /* or 2 nearly identical points were found */
406 : /* Try with Green Formula as a fallback, but this is not a guarantee */
407 : /* as we'll probably be affected by numerical instabilities */
408 :
409 6 : double dfSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
410 :
411 3730 : for (i=1; i<nPointCount-1; i++) {
412 3724 : dfSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
413 : }
414 :
415 6 : dfSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
416 :
417 6 : return dfSum < 0;
418 : }
419 :
420 : /************************************************************************/
421 : /* reverseWindingOrder() */
422 : /************************************************************************/
423 :
424 0 : void OGRLinearRing::reverseWindingOrder()
425 :
426 : {
427 0 : int pos = 0;
428 0 : OGRPoint tempPoint;
429 :
430 0 : for( int i = 0; i < nPointCount / 2; i++ )
431 : {
432 0 : getPoint( i, &tempPoint );
433 0 : pos = nPointCount - i - 1;
434 0 : if( getCoordinateDimension() == 2 )
435 : {
436 0 : setPoint( i, getX(pos), getY(pos) );
437 0 : setPoint( pos, tempPoint.getX(), tempPoint.getY() );
438 : }
439 : else
440 : {
441 0 : setPoint( i, getX(pos), getY(pos), getZ(pos) );
442 0 : setPoint( pos, tempPoint.getX(), tempPoint.getY(), tempPoint.getZ() );
443 : }
444 0 : }
445 0 : }
446 :
447 : /************************************************************************/
448 : /* closeRing() */
449 : /************************************************************************/
450 :
451 1684 : void OGRLinearRing::closeRings()
452 :
453 : {
454 1684 : if( nPointCount < 2 )
455 0 : return;
456 :
457 1684 : if( getX(0) != getX(nPointCount-1)
458 : || getY(0) != getY(nPointCount-1)
459 : || getZ(0) != getZ(nPointCount-1) )
460 : {
461 : /* Avoid implicit change of coordinate dimensionality
462 : * if z=0.0 and dim=2
463 : */
464 476 : if( getCoordinateDimension() == 2 )
465 474 : addPoint( getX(0), getY(0) );
466 : else
467 2 : addPoint( getX(0), getY(0), getZ(0) );
468 :
469 : }
470 : }
471 :
472 : /************************************************************************/
473 : /* get_Area() */
474 : /************************************************************************/
475 :
476 : /**
477 : * \brief Compute area of ring.
478 : *
479 : * The area is computed according to Green's Theorem:
480 : *
481 : * Area is "Sum(x(i)*(y(i+1) - y(i-1)))/2" for i = 0 to pointCount-1,
482 : * assuming the last point is a duplicate of the first.
483 : *
484 : * @return computed area.
485 : */
486 :
487 706 : double OGRLinearRing::get_Area() const
488 :
489 : {
490 706 : double dfAreaSum = 0.0;
491 : int i;
492 :
493 706 : if( nPointCount < 2 )
494 2 : return 0;
495 :
496 704 : dfAreaSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
497 :
498 9090 : for( i = 1; i < nPointCount-1; i++ )
499 : {
500 8386 : dfAreaSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
501 : }
502 :
503 704 : dfAreaSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
504 :
505 704 : return 0.5 * fabs(dfAreaSum);
506 : }
507 :
508 : /************************************************************************/
509 : /* isPointInRing() */
510 : /************************************************************************/
511 :
512 520 : OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint, int bTestEnvelope) const
513 : {
514 520 : if ( NULL == poPoint )
515 : {
516 0 : CPLDebug( "OGR", "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - passed point is NULL!" );
517 0 : return 0;
518 : }
519 :
520 520 : const int iNumPoints = getNumPoints();
521 :
522 : // Simple validation
523 520 : if ( iNumPoints < 4 )
524 0 : return 0;
525 :
526 520 : const double dfTestX = poPoint->getX();
527 520 : const double dfTestY = poPoint->getY();
528 :
529 : // Fast test if point is inside extent of the ring
530 520 : if (bTestEnvelope)
531 : {
532 198 : OGREnvelope extent;
533 198 : getEnvelope(&extent);
534 198 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
535 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
536 : {
537 0 : return 0;
538 : }
539 : }
540 :
541 : // For every point p in ring,
542 : // test if ray starting from given point crosses segment (p - 1, p)
543 520 : int iNumCrossings = 0;
544 :
545 157330 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
546 : {
547 156810 : const int iPointPrev = iPoint - 1;
548 :
549 156810 : const double x1 = getX(iPoint) - dfTestX;
550 156810 : const double y1 = getY(iPoint) - dfTestY;
551 :
552 156810 : const double x2 = getX(iPointPrev) - dfTestX;
553 156810 : const double y2 = getY(iPointPrev) - dfTestY;
554 :
555 156810 : if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) )
556 : {
557 : // Check if ray intersects with segment of the ring
558 1204 : const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
559 1204 : if ( 0.0 < dfIntersection )
560 : {
561 : // Count intersections
562 646 : iNumCrossings++;
563 : }
564 : }
565 : }
566 :
567 : // If iNumCrossings number is even, given point is outside the ring,
568 : // when the crossings number is odd, the point is inside the ring.
569 520 : return ( ( iNumCrossings % 2 ) == 1 ? 1 : 0 );
570 : }
571 :
572 : /************************************************************************/
573 : /* isPointOnRingBoundary() */
574 : /************************************************************************/
575 :
576 14 : OGRBoolean OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint, int bTestEnvelope) const
577 : {
578 14 : if ( NULL == poPoint )
579 : {
580 0 : CPLDebug( "OGR", "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint) - passed point is NULL!" );
581 0 : return 0;
582 : }
583 :
584 14 : const int iNumPoints = getNumPoints();
585 :
586 : // Simple validation
587 14 : if ( iNumPoints < 4 )
588 0 : return 0;
589 :
590 14 : const double dfTestX = poPoint->getX();
591 14 : const double dfTestY = poPoint->getY();
592 :
593 : // Fast test if point is inside extent of the ring
594 14 : OGREnvelope extent;
595 14 : getEnvelope(&extent);
596 14 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
597 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
598 : {
599 0 : return 0;
600 : }
601 :
602 26 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
603 : {
604 26 : const int iPointPrev = iPoint - 1;
605 :
606 26 : const double x1 = getX(iPoint) - dfTestX;
607 26 : const double y1 = getY(iPoint) - dfTestY;
608 :
609 26 : const double x2 = getX(iPointPrev) - dfTestX;
610 26 : const double y2 = getY(iPointPrev) - dfTestY;
611 :
612 : /* If iPoint and iPointPrev are the same, go on */
613 26 : if (x1 == x2 && y1 == y2)
614 : {
615 0 : continue;
616 : }
617 :
618 : /* If the point is on the segment, return immediatly */
619 : /* FIXME? If the test point is not exactly identical to one of */
620 : /* the vertices of the ring, but somewhere on a segment, there's */
621 : /* little chance that we get 0. So that should be tested against some epsilon */
622 26 : if ( x1 * y2 - x2 * y1 == 0 )
623 : {
624 14 : return 1;
625 : }
626 : }
627 :
628 0 : return 0;
629 : }
|
