1 : /******************************************************************************
2 : * $Id: ogrlinearring.cpp 18674 2010-01-27 23:16:09Z 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 18674 2010-01-27 23:16:09Z rouault $");
34 :
35 : /************************************************************************/
36 : /* OGRLinearRing() */
37 : /************************************************************************/
38 :
39 3755 : OGRLinearRing::OGRLinearRing()
40 :
41 : {
42 3755 : }
43 :
44 : /************************************************************************/
45 : /* ~OGRLinearRing() */
46 : /************************************************************************/
47 5771 : OGRLinearRing::~OGRLinearRing()
48 :
49 : {
50 5771 : }
51 :
52 : /************************************************************************/
53 : /* OGRLinearRing() */
54 : /************************************************************************/
55 :
56 2016 : OGRLinearRing::OGRLinearRing( OGRLinearRing * poSrcRing )
57 :
58 : {
59 2016 : if( poSrcRing == NULL )
60 : {
61 0 : CPLDebug( "OGR", "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - passed in ring is NULL!" );
62 0 : return;
63 : }
64 :
65 2016 : setNumPoints( poSrcRing->getNumPoints() );
66 :
67 : memcpy( paoPoints, poSrcRing->paoPoints,
68 2016 : sizeof(OGRRawPoint) * getNumPoints() );
69 :
70 2016 : if( poSrcRing->padfZ )
71 : {
72 120 : Make3D();
73 120 : memcpy( padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints() );
74 : }
75 0 : }
76 :
77 : /************************************************************************/
78 : /* getGeometryName() */
79 : /************************************************************************/
80 :
81 1268 : const char * OGRLinearRing::getGeometryName() const
82 :
83 : {
84 1268 : 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 : OGRErr OGRLinearRing::exportToWkb( OGRwkbByteOrder eByteOrder,
121 0 : 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 : OGRErr OGRLinearRing::_importFromWkb( OGRwkbByteOrder eByteOrder, int b3D,
138 : unsigned char * pabyData,
139 306 : int nBytesAvailable )
140 :
141 : {
142 306 : 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 306 : memcpy( &nNewNumPoints, pabyData, 4 );
151 :
152 306 : if( OGR_SWAP( eByteOrder ) )
153 17 : 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 306 : int nPointSize = (b3D ? 24 : 16);
160 306 : if (nNewNumPoints < 0 || nNewNumPoints > INT_MAX / nPointSize)
161 0 : return OGRERR_CORRUPT_DATA;
162 306 : int nBufferMinSize = nPointSize * nNewNumPoints;
163 :
164 306 : 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 306 : setNumPoints( nNewNumPoints );
173 :
174 306 : if( b3D )
175 19 : Make3D();
176 : else
177 287 : Make2D();
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Get the vertices */
181 : /* -------------------------------------------------------------------- */
182 306 : int i = 0;
183 :
184 306 : if( !b3D )
185 : {
186 287 : memcpy( paoPoints, pabyData + 4, 16 * nPointCount );
187 : }
188 : else
189 : {
190 108 : for( int i = 0; i < nPointCount; i++ )
191 : {
192 89 : memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
193 89 : memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
194 89 : memcpy( padfZ + i, pabyData + 4 + 24 * i + 16, 8 );
195 : }
196 : }
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Byte swap if needed. */
200 : /* -------------------------------------------------------------------- */
201 306 : if( OGR_SWAP( eByteOrder ) )
202 : {
203 88 : for( i = 0; i < nPointCount; i++ )
204 : {
205 71 : CPL_SWAPDOUBLE( &(paoPoints[i].x) );
206 71 : CPL_SWAPDOUBLE( &(paoPoints[i].y) );
207 :
208 71 : if( b3D )
209 : {
210 13 : CPL_SWAPDOUBLE( padfZ + i );
211 : }
212 : }
213 : }
214 :
215 306 : 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 : OGRErr OGRLinearRing::_exportToWkb( OGRwkbByteOrder eByteOrder, int b3D,
226 443 : unsigned char * pabyData ) const
227 :
228 : {
229 : int i, nWords;
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Copy in the raw data. */
233 : /* -------------------------------------------------------------------- */
234 443 : memcpy( pabyData, &nPointCount, 4 );
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Copy in the raw data. */
238 : /* -------------------------------------------------------------------- */
239 443 : if( b3D )
240 : {
241 40 : nWords = 3 * nPointCount;
242 626 : for( i = 0; i < nPointCount; i++ )
243 : {
244 586 : memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
245 586 : memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
246 586 : if( padfZ == NULL )
247 0 : memset( pabyData+4+i*24+16, 0, 8 );
248 : else
249 586 : memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
250 : }
251 : }
252 : else
253 : {
254 403 : nWords = 2 * nPointCount;
255 403 : memcpy( pabyData+4, paoPoints, 16 * nPointCount );
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Swap if needed. */
260 : /* -------------------------------------------------------------------- */
261 443 : if( OGR_SWAP( eByteOrder ) )
262 : {
263 : int nCount;
264 :
265 14 : nCount = CPL_SWAP32( nPointCount );
266 14 : memcpy( pabyData, &nCount, 4 );
267 :
268 143 : for( i = 0; i < nWords; i++ )
269 : {
270 129 : CPL_SWAPDOUBLE( pabyData + 4 + 8 * i );
271 : }
272 : }
273 :
274 443 : return OGRERR_NONE;
275 : }
276 :
277 : /************************************************************************/
278 : /* _WkbSize() */
279 : /* */
280 : /* Helper method for OGRPolygon. NOT THE NORMAL WkbSize() METHOD! */
281 : /************************************************************************/
282 :
283 1668 : int OGRLinearRing::_WkbSize( int b3D ) const
284 :
285 : {
286 1668 : if( b3D )
287 150 : return 4 + 24 * nPointCount;
288 : else
289 1518 : 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 317 : OGRGeometry *OGRLinearRing::clone() const
300 :
301 : {
302 : OGRLinearRing *poNewLinearRing;
303 :
304 317 : poNewLinearRing = new OGRLinearRing();
305 317 : poNewLinearRing->assignSpatialReference( getSpatialReference() );
306 :
307 317 : poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ );
308 :
309 317 : return poNewLinearRing;
310 : }
311 :
312 :
313 : /************************************************************************/
314 : /* epsilonEqual() */
315 : /************************************************************************/
316 :
317 : static const double EPSILON = 1E-5;
318 :
319 695 : static inline bool epsilonEqual(double a, double b, double eps)
320 : {
321 695 : 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 297 : int OGRLinearRing::isClockwise() const
335 :
336 : {
337 : int i, v, next;
338 : double dx0, dy0, dx1, dy1, crossproduct;
339 :
340 297 : if( nPointCount < 2 )
341 0 : return TRUE;
342 :
343 : /* Find the lowest rightmost vertex */
344 297 : v = 0;
345 3753 : for ( i = 1; i < nPointCount - 1; i++ )
346 : {
347 : /* => v < end */
348 3456 : if ( paoPoints[i].y< paoPoints[v].y ||
349 : ( paoPoints[i].y== paoPoints[v].y &&
350 : paoPoints[i].x > paoPoints[v].x ) )
351 : {
352 621 : v = i;
353 : }
354 : }
355 :
356 : /* Vertices may be duplicate, we have to go to nearest different in each direction */
357 : /* preceding */
358 297 : next = v - 1;
359 0 : while ( 1 )
360 : {
361 297 : if ( next < 0 )
362 : {
363 66 : next = nPointCount - 1 - 1;
364 : }
365 :
366 297 : if( !epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON)
367 : || !epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
368 : {
369 297 : break;
370 : }
371 :
372 0 : if ( next == v ) /* So we cannot get into endless loop */
373 : {
374 0 : break;
375 : }
376 :
377 0 : next--;
378 : }
379 :
380 297 : dx0 = paoPoints[next].x - paoPoints[v].x;
381 297 : dy0 = paoPoints[next].y - paoPoints[v].y;
382 :
383 :
384 : /* following */
385 297 : next = v + 1;
386 1 : while ( 1 )
387 : {
388 298 : if ( next >= nPointCount - 1 )
389 : {
390 48 : next = 0;
391 : }
392 :
393 298 : if ( !epsilonEqual(paoPoints[next].x, paoPoints[v].x, EPSILON)
394 : || !epsilonEqual(paoPoints[next].y, paoPoints[v].y, EPSILON) )
395 : {
396 297 : break;
397 : }
398 :
399 1 : if ( next == v ) /* So we cannot get into endless loop */
400 : {
401 0 : break;
402 : }
403 :
404 1 : next++;
405 : }
406 :
407 297 : dx1 = paoPoints[next].x - paoPoints[v].x;
408 297 : dy1 = paoPoints[next].y - paoPoints[v].y;
409 :
410 297 : crossproduct = dx1 * dy0 - dx0 * dy1;
411 :
412 297 : if ( crossproduct > 0 ) /* CCW */
413 73 : return FALSE;
414 224 : else if ( crossproduct < 0 ) /* CW */
415 224 : return TRUE;
416 :
417 : /* ok, this is a degenerate case : the extent of the polygon is less than EPSILON */
418 : /* Try with Green Formula as a fallback, but this is not a guarantee */
419 : /* as we'll probably be affected by numerical instabilities */
420 :
421 0 : double dfSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
422 :
423 0 : for (i=1; i<nPointCount-1; i++) {
424 0 : dfSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
425 : }
426 :
427 0 : dfSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
428 :
429 0 : return dfSum < 0;
430 : }
431 :
432 : /************************************************************************/
433 : /* reverseWindingOrder() */
434 : /************************************************************************/
435 :
436 0 : void OGRLinearRing::reverseWindingOrder()
437 :
438 : {
439 0 : int pos = 0;
440 0 : OGRPoint tempPoint;
441 :
442 0 : for( int i = 0; i < nPointCount / 2; i++ )
443 : {
444 0 : getPoint( i, &tempPoint );
445 0 : pos = nPointCount - i - 1;
446 0 : setPoint( i, getX(pos), getY(pos), getZ(pos) );
447 0 : setPoint( pos, tempPoint.getX(), tempPoint.getY(), tempPoint.getZ() );
448 0 : }
449 0 : }
450 :
451 : /************************************************************************/
452 : /* closeRing() */
453 : /************************************************************************/
454 :
455 115 : void OGRLinearRing::closeRings()
456 :
457 : {
458 115 : if( nPointCount < 2 )
459 0 : return;
460 :
461 115 : if( getX(0) != getX(nPointCount-1)
462 : || getY(0) != getY(nPointCount-1)
463 : || getZ(0) != getZ(nPointCount-1) )
464 : {
465 : /* Avoid implicit change of coordinate dimensionality
466 : * if z=0.0 and dim=2
467 : */
468 21 : if( getCoordinateDimension() == 2 )
469 20 : addPoint( getX(0), getY(0) );
470 : else
471 1 : addPoint( getX(0), getY(0), getZ(0) );
472 :
473 : }
474 : }
475 :
476 : /************************************************************************/
477 : /* get_Area() */
478 : /************************************************************************/
479 :
480 : /**
481 : * \brief Compute area of ring.
482 : *
483 : * The area is computed according to Green's Theorem:
484 : *
485 : * Area is "Sum(x(i)*(y(i+1) - y(i-1)))/2" for i = 0 to pointCount-1,
486 : * assuming the last point is a duplicate of the first.
487 : *
488 : * @return computed area.
489 : */
490 :
491 307 : double OGRLinearRing::get_Area() const
492 :
493 : {
494 307 : double dfAreaSum = 0.0;
495 : int i;
496 :
497 307 : if( nPointCount < 2 )
498 1 : return 0;
499 :
500 306 : dfAreaSum = paoPoints[0].x * (paoPoints[1].y - paoPoints[nPointCount-1].y);
501 :
502 3786 : for( i = 1; i < nPointCount-1; i++ )
503 : {
504 3480 : dfAreaSum += paoPoints[i].x * (paoPoints[i+1].y - paoPoints[i-1].y);
505 : }
506 :
507 306 : dfAreaSum += paoPoints[nPointCount-1].x * (paoPoints[0].y - paoPoints[nPointCount-2].y);
508 :
509 306 : return 0.5 * fabs(dfAreaSum);
510 : }
511 :
512 : /************************************************************************/
513 : /* isPointInRing() */
514 : /************************************************************************/
515 :
516 259 : OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint, int bTestEnvelope) const
517 : {
518 259 : if ( NULL == poPoint )
519 : {
520 0 : CPLDebug( "OGR", "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - passed point is NULL!" );
521 0 : return 0;
522 : }
523 :
524 259 : const int iNumPoints = getNumPoints();
525 :
526 : // Simple validation
527 259 : if ( iNumPoints < 4 )
528 0 : return 0;
529 :
530 259 : const double dfTestX = poPoint->getX();
531 259 : const double dfTestY = poPoint->getY();
532 :
533 : // Fast test if point is inside extent of the ring
534 259 : if (bTestEnvelope)
535 : {
536 99 : OGREnvelope extent;
537 99 : getEnvelope(&extent);
538 99 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
539 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
540 : {
541 0 : return 0;
542 : }
543 : }
544 :
545 : // For every point p in ring,
546 : // test if ray starting from given point crosses segment (p - 1, p)
547 259 : int iNumCrossings = 0;
548 :
549 78660 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
550 : {
551 78401 : const int iPointPrev = iPoint - 1;
552 :
553 78401 : const double x1 = getX(iPoint) - dfTestX;
554 78401 : const double y1 = getY(iPoint) - dfTestY;
555 :
556 78401 : const double x2 = getX(iPointPrev) - dfTestX;
557 78401 : const double y2 = getY(iPointPrev) - dfTestY;
558 :
559 78401 : if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) )
560 : {
561 : // Check if ray intersects with segment of the ring
562 600 : const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
563 600 : if ( 0.0 < dfIntersection )
564 : {
565 : // Count intersections
566 322 : iNumCrossings++;
567 : }
568 : }
569 : }
570 :
571 : // If iNumCrossings number is even, given point is outside the ring,
572 : // when the crossings number is odd, the point is inside the ring.
573 259 : return ( ( iNumCrossings % 2 ) == 1 ? 1 : 0 );
574 : }
575 :
576 : /************************************************************************/
577 : /* isPointOnRingBoundary() */
578 : /************************************************************************/
579 :
580 7 : OGRBoolean OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint, int bTestEnvelope) const
581 : {
582 7 : if ( NULL == poPoint )
583 : {
584 0 : CPLDebug( "OGR", "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint) - passed point is NULL!" );
585 0 : return 0;
586 : }
587 :
588 7 : const int iNumPoints = getNumPoints();
589 :
590 : // Simple validation
591 7 : if ( iNumPoints < 4 )
592 0 : return 0;
593 :
594 7 : const double dfTestX = poPoint->getX();
595 7 : const double dfTestY = poPoint->getY();
596 :
597 : // Fast test if point is inside extent of the ring
598 7 : OGREnvelope extent;
599 7 : getEnvelope(&extent);
600 7 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
601 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
602 : {
603 0 : return 0;
604 : }
605 :
606 13 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
607 : {
608 13 : const int iPointPrev = iPoint - 1;
609 :
610 13 : const double x1 = getX(iPoint) - dfTestX;
611 13 : const double y1 = getY(iPoint) - dfTestY;
612 :
613 13 : const double x2 = getX(iPointPrev) - dfTestX;
614 13 : const double y2 = getY(iPointPrev) - dfTestY;
615 :
616 : /* If iPoint and iPointPrev are the same, go on */
617 13 : if (x1 == x2 && y1 == y2)
618 : {
619 0 : continue;
620 : }
621 :
622 : /* If the point is on the segment, return immediatly */
623 : /* FIXME? If the test point is not exactly identical to one of */
624 : /* the vertices of the ring, but somewhere on a segment, there's */
625 : /* little chance that we get 0. So that should be tested against some epsilon */
626 13 : if ( x1 * y2 - x2 * y1 == 0 )
627 : {
628 7 : return 1;
629 : }
630 : }
631 :
632 0 : return 0;
633 : }
|
