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