1 : /******************************************************************************
2 : * $Id: ogrlinearring.cpp 24520 2012-05-30 21:30:11Z 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 24520 2012-05-30 21:30:11Z rouault $");
34 :
35 : /************************************************************************/
36 : /* OGRLinearRing() */
37 : /************************************************************************/
38 :
39 82055 : OGRLinearRing::OGRLinearRing()
40 :
41 : {
42 82055 : }
43 :
44 : /************************************************************************/
45 : /* ~OGRLinearRing() */
46 : /************************************************************************/
47 164721 : OGRLinearRing::~OGRLinearRing()
48 :
49 : {
50 164721 : }
51 :
52 : /************************************************************************/
53 : /* OGRLinearRing() */
54 : /************************************************************************/
55 :
56 82666 : OGRLinearRing::OGRLinearRing( OGRLinearRing * poSrcRing )
57 :
58 : {
59 82666 : if( poSrcRing == NULL )
60 : {
61 0 : CPLDebug( "OGR", "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - passed in ring is NULL!" );
62 0 : return;
63 : }
64 :
65 82666 : setNumPoints( poSrcRing->getNumPoints() );
66 :
67 : memcpy( paoPoints, poSrcRing->paoPoints,
68 82666 : sizeof(OGRRawPoint) * getNumPoints() );
69 :
70 82666 : if( poSrcRing->padfZ )
71 : {
72 194 : Make3D();
73 194 : memcpy( padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints() );
74 : }
75 0 : }
76 :
77 : /************************************************************************/
78 : /* getGeometryName() */
79 : /************************************************************************/
80 :
81 3213 : const char * OGRLinearRing::getGeometryName() const
82 :
83 : {
84 3213 : 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 1681 : OGRErr OGRLinearRing::_importFromWkb( OGRwkbByteOrder eByteOrder, int b3D,
138 : unsigned char * pabyData,
139 : int nBytesAvailable )
140 :
141 : {
142 1681 : 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 1681 : memcpy( &nNewNumPoints, pabyData, 4 );
151 :
152 1681 : if( OGR_SWAP( eByteOrder ) )
153 19 : 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 1681 : int nPointSize = (b3D ? 24 : 16);
160 1681 : if (nNewNumPoints < 0 || nNewNumPoints > INT_MAX / nPointSize)
161 0 : return OGRERR_CORRUPT_DATA;
162 1681 : int nBufferMinSize = nPointSize * nNewNumPoints;
163 :
164 1681 : 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 1681 : setNumPoints( nNewNumPoints );
173 :
174 1681 : if( b3D )
175 297 : Make3D();
176 : else
177 1384 : Make2D();
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Get the vertices */
181 : /* -------------------------------------------------------------------- */
182 1681 : int i = 0;
183 :
184 1681 : if( !b3D )
185 : {
186 1384 : memcpy( paoPoints, pabyData + 4, 16 * nPointCount );
187 : }
188 : else
189 : {
190 5766 : for( int i = 0; i < nPointCount; i++ )
191 : {
192 5469 : memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
193 5469 : memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
194 5469 : memcpy( padfZ + i, pabyData + 4 + 24 * i + 16, 8 );
195 : }
196 : }
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Byte swap if needed. */
200 : /* -------------------------------------------------------------------- */
201 1681 : if( OGR_SWAP( eByteOrder ) )
202 : {
203 98 : for( i = 0; i < nPointCount; i++ )
204 : {
205 79 : CPL_SWAPDOUBLE( &(paoPoints[i].x) );
206 79 : CPL_SWAPDOUBLE( &(paoPoints[i].y) );
207 :
208 79 : if( b3D )
209 : {
210 13 : CPL_SWAPDOUBLE( padfZ + i );
211 : }
212 : }
213 : }
214 :
215 1681 : 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 20985 : 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 20985 : memcpy( pabyData, &nPointCount, 4 );
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Copy in the raw data. */
238 : /* -------------------------------------------------------------------- */
239 20985 : if( b3D )
240 : {
241 73 : nWords = 3 * nPointCount;
242 1547 : for( i = 0; i < nPointCount; i++ )
243 : {
244 1474 : memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
245 1474 : memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
246 1474 : if( padfZ == NULL )
247 0 : memset( pabyData+4+i*24+16, 0, 8 );
248 : else
249 1474 : memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
250 : }
251 : }
252 : else
253 : {
254 20912 : nWords = 2 * nPointCount;
255 20912 : memcpy( pabyData+4, paoPoints, 16 * nPointCount );
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Swap if needed. */
260 : /* -------------------------------------------------------------------- */
261 20985 : if( OGR_SWAP( eByteOrder ) )
262 : {
263 : int nCount;
264 :
265 16 : nCount = CPL_SWAP32( nPointCount );
266 16 : memcpy( pabyData, &nCount, 4 );
267 :
268 161 : for( i = 0; i < nWords; i++ )
269 : {
270 145 : CPL_SWAPDOUBLE( pabyData + 4 + 8 * i );
271 : }
272 : }
273 :
274 20985 : return OGRERR_NONE;
275 : }
276 :
277 : /************************************************************************/
278 : /* _WkbSize() */
279 : /* */
280 : /* Helper method for OGRPolygon. NOT THE NORMAL WkbSize() METHOD! */
281 : /************************************************************************/
282 :
283 45729 : int OGRLinearRing::_WkbSize( int b3D ) const
284 :
285 : {
286 45729 : if( b3D )
287 827 : return 4 + 24 * nPointCount;
288 : else
289 44902 : 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 417 : OGRGeometry *OGRLinearRing::clone() const
300 :
301 : {
302 : OGRLinearRing *poNewLinearRing;
303 :
304 417 : poNewLinearRing = new OGRLinearRing();
305 417 : poNewLinearRing->assignSpatialReference( getSpatialReference() );
306 :
307 417 : poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ );
308 :
309 417 : return poNewLinearRing;
310 : }
311 :
312 :
313 : /************************************************************************/
314 : /* epsilonEqual() */
315 : /************************************************************************/
316 :
317 : static const double EPSILON = 1E-5;
318 :
319 3761 : static inline bool epsilonEqual(double a, double b, double eps)
320 : {
321 3761 : 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 1776 : int OGRLinearRing::isClockwise() const
335 :
336 : {
337 : int i, v, next;
338 : double dx0, dy0, dx1, dy1, crossproduct;
339 1776 : int bUseFallback = FALSE;
340 :
341 1776 : if( nPointCount < 2 )
342 0 : return TRUE;
343 :
344 : /* Find the lowest rightmost vertex */
345 1776 : v = 0;
346 398581 : for ( i = 1; i < nPointCount - 1; i++ )
347 : {
348 : /* => v < end */
349 1029939 : if ( paoPoints[i].y< paoPoints[v].y ||
350 631530 : ( paoPoints[i].y== paoPoints[v].y &&
351 1604 : paoPoints[i].x > paoPoints[v].x ) )
352 : {
353 81579 : v = i;
354 : }
355 : }
356 :
357 : /* previous */
358 1776 : next = v - 1;
359 1776 : if ( next < 0 )
360 : {
361 329 : next = nPointCount - 1 - 1;
362 : }
363 :
364 2054 : if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
365 278 : 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 2 : bUseFallback = TRUE;
370 : }
371 :
372 1776 : dx0 = paoPoints[next].x - paoPoints[v].x;
373 1776 : dy0 = paoPoints[next].y - paoPoints[v].y;
374 :
375 :
376 : /* following */
377 1776 : next = v + 1;
378 1776 : if ( next >= nPointCount - 1 )
379 : {
380 133 : next = 0;
381 : }
382 :
383 1916 : if( epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON) &&
384 140 : 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 3 : bUseFallback = TRUE;
389 : }
390 :
391 1776 : dx1 = paoPoints[next].x - paoPoints[v].x;
392 1776 : dy1 = paoPoints[next].y - paoPoints[v].y;
393 :
394 1776 : crossproduct = dx1 * dy0 - dx0 * dy1;
395 :
396 1776 : if (!bUseFallback)
397 : {
398 1773 : if ( crossproduct > 0 ) /* CCW */
399 1041 : return FALSE;
400 732 : else if ( crossproduct < 0 ) /* CW */
401 732 : 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 3 : double dfSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
410 :
411 1865 : for (i=1; i<nPointCount-1; i++) {
412 1862 : dfSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
413 : }
414 :
415 3 : dfSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
416 :
417 3 : 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 pointA, pointB;
429 :
430 0 : for( int i = 0; i < nPointCount / 2; i++ )
431 : {
432 0 : getPoint( i, &pointA );
433 0 : pos = nPointCount - i - 1;
434 0 : getPoint( pos, &pointB );
435 0 : setPoint( i, &pointB );
436 0 : setPoint( pos, &pointA );
437 0 : }
438 0 : }
439 :
440 : /************************************************************************/
441 : /* closeRing() */
442 : /************************************************************************/
443 :
444 886 : void OGRLinearRing::closeRings()
445 :
446 : {
447 886 : if( nPointCount < 2 )
448 1 : return;
449 :
450 885 : if( getX(0) != getX(nPointCount-1)
451 : || getY(0) != getY(nPointCount-1)
452 : || getZ(0) != getZ(nPointCount-1) )
453 : {
454 261 : OGRPoint oFirstPoint;
455 261 : getPoint( 0, &oFirstPoint );
456 261 : addPoint( &oFirstPoint );
457 : }
458 : }
459 :
460 : /************************************************************************/
461 : /* get_Area() */
462 : /************************************************************************/
463 :
464 : /**
465 : * \brief Compute area of ring.
466 : *
467 : * The area is computed according to Green's Theorem:
468 : *
469 : * Area is "Sum(x(i)*(y(i+1) - y(i-1)))/2" for i = 0 to pointCount-1,
470 : * assuming the last point is a duplicate of the first.
471 : *
472 : * @return computed area.
473 : */
474 :
475 1759 : double OGRLinearRing::get_Area() const
476 :
477 : {
478 1759 : double dfAreaSum = 0.0;
479 : int i;
480 :
481 1759 : if( nPointCount < 2 )
482 1 : return 0;
483 :
484 1758 : dfAreaSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
485 :
486 398664 : for( i = 1; i < nPointCount-1; i++ )
487 : {
488 396906 : dfAreaSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
489 : }
490 :
491 1758 : dfAreaSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
492 :
493 1758 : return 0.5 * fabs(dfAreaSum);
494 : }
495 :
496 : /************************************************************************/
497 : /* isPointInRing() */
498 : /************************************************************************/
499 :
500 280 : OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint, int bTestEnvelope) const
501 : {
502 280 : if ( NULL == poPoint )
503 : {
504 0 : CPLDebug( "OGR", "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - passed point is NULL!" );
505 0 : return 0;
506 : }
507 :
508 280 : const int iNumPoints = getNumPoints();
509 :
510 : // Simple validation
511 280 : if ( iNumPoints < 4 )
512 0 : return 0;
513 :
514 280 : const double dfTestX = poPoint->getX();
515 280 : const double dfTestY = poPoint->getY();
516 :
517 : // Fast test if point is inside extent of the ring
518 280 : if (bTestEnvelope)
519 : {
520 99 : OGREnvelope extent;
521 99 : getEnvelope(&extent);
522 99 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
523 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
524 : {
525 0 : return 0;
526 : }
527 : }
528 :
529 : // For every point p in ring,
530 : // test if ray starting from given point crosses segment (p - 1, p)
531 280 : int iNumCrossings = 0;
532 :
533 78765 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
534 : {
535 78485 : const int iPointPrev = iPoint - 1;
536 :
537 78485 : const double x1 = getX(iPoint) - dfTestX;
538 78485 : const double y1 = getY(iPoint) - dfTestY;
539 :
540 78485 : const double x2 = getX(iPointPrev) - dfTestX;
541 78485 : const double y2 = getY(iPointPrev) - dfTestY;
542 :
543 78485 : if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) )
544 : {
545 : // Check if ray intersects with segment of the ring
546 642 : const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
547 642 : if ( 0.0 < dfIntersection )
548 : {
549 : // Count intersections
550 343 : iNumCrossings++;
551 : }
552 : }
553 : }
554 :
555 : // If iNumCrossings number is even, given point is outside the ring,
556 : // when the crossings number is odd, the point is inside the ring.
557 280 : return ( ( iNumCrossings % 2 ) == 1 ? 1 : 0 );
558 : }
559 :
560 : /************************************************************************/
561 : /* isPointOnRingBoundary() */
562 : /************************************************************************/
563 :
564 7 : OGRBoolean OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint, int bTestEnvelope) const
565 : {
566 7 : if ( NULL == poPoint )
567 : {
568 0 : CPLDebug( "OGR", "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint) - passed point is NULL!" );
569 0 : return 0;
570 : }
571 :
572 7 : const int iNumPoints = getNumPoints();
573 :
574 : // Simple validation
575 7 : if ( iNumPoints < 4 )
576 0 : return 0;
577 :
578 7 : const double dfTestX = poPoint->getX();
579 7 : const double dfTestY = poPoint->getY();
580 :
581 : // Fast test if point is inside extent of the ring
582 7 : OGREnvelope extent;
583 7 : getEnvelope(&extent);
584 7 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
585 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
586 : {
587 0 : return 0;
588 : }
589 :
590 13 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
591 : {
592 13 : const int iPointPrev = iPoint - 1;
593 :
594 13 : const double x1 = getX(iPoint) - dfTestX;
595 13 : const double y1 = getY(iPoint) - dfTestY;
596 :
597 13 : const double x2 = getX(iPointPrev) - dfTestX;
598 13 : const double y2 = getY(iPointPrev) - dfTestY;
599 :
600 : /* If iPoint and iPointPrev are the same, go on */
601 13 : if (x1 == x2 && y1 == y2)
602 : {
603 0 : continue;
604 : }
605 :
606 : /* If the point is on the segment, return immediatly */
607 : /* FIXME? If the test point is not exactly identical to one of */
608 : /* the vertices of the ring, but somewhere on a segment, there's */
609 : /* little chance that we get 0. So that should be tested against some epsilon */
610 13 : if ( x1 * y2 - x2 * y1 == 0 )
611 : {
612 7 : return 1;
613 : }
614 : }
615 :
616 0 : return 0;
617 : }
|
