1 : /******************************************************************************
2 : * $Id: ogrlinearring.cpp 16853 2009-04-26 12:22:10Z 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 16853 2009-04-26 12:22:10Z rouault $");
34 :
35 : /************************************************************************/
36 : /* OGRLinearRing() */
37 : /************************************************************************/
38 :
39 2786 : OGRLinearRing::OGRLinearRing()
40 :
41 : {
42 2786 : }
43 :
44 : /************************************************************************/
45 : /* ~OGRLinearRing() */
46 : /************************************************************************/
47 8402 : OGRLinearRing::~OGRLinearRing()
48 :
49 : {
50 8402 : }
51 :
52 : /************************************************************************/
53 : /* OGRLinearRing() */
54 : /************************************************************************/
55 :
56 1434 : OGRLinearRing::OGRLinearRing( OGRLinearRing * poSrcRing )
57 :
58 : {
59 1434 : if( poSrcRing == NULL )
60 : {
61 0 : CPLDebug( "OGR", "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - passed in ring is NULL!" );
62 0 : return;
63 : }
64 :
65 1434 : setNumPoints( poSrcRing->getNumPoints() );
66 :
67 : memcpy( paoPoints, poSrcRing->paoPoints,
68 1434 : sizeof(OGRRawPoint) * getNumPoints() );
69 :
70 1434 : if( poSrcRing->padfZ )
71 : {
72 73 : Make3D();
73 73 : memcpy( padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints() );
74 : }
75 0 : }
76 :
77 : /************************************************************************/
78 : /* getGeometryName() */
79 : /************************************************************************/
80 :
81 832 : const char * OGRLinearRing::getGeometryName() const
82 :
83 : {
84 832 : return "LINEARRING";
85 : }
86 :
87 : /************************************************************************/
88 : /* WkbSize() */
89 : /* */
90 : /* Disable this method. */
91 : /************************************************************************/
92 :
93 194 : int OGRLinearRing::WkbSize() const
94 :
95 : {
96 194 : 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 223 : OGRErr OGRLinearRing::_importFromWkb( OGRwkbByteOrder eByteOrder, int b3D,
138 : unsigned char * pabyData,
139 : int nBytesAvailable )
140 :
141 : {
142 223 : 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 223 : memcpy( &nNewNumPoints, pabyData, 4 );
151 :
152 223 : 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 223 : int nPointSize = (b3D ? 24 : 16);
160 223 : if (nNewNumPoints < 0 || nNewNumPoints > INT_MAX / nPointSize)
161 0 : return OGRERR_CORRUPT_DATA;
162 223 : int nBufferMinSize = nPointSize * nNewNumPoints;
163 :
164 223 : 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 223 : setNumPoints( nNewNumPoints );
173 :
174 223 : if( b3D )
175 16 : Make3D();
176 : else
177 207 : Make2D();
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Get the vertices */
181 : /* -------------------------------------------------------------------- */
182 223 : int i = 0;
183 :
184 223 : if( !b3D )
185 : {
186 207 : memcpy( paoPoints, pabyData + 4, 16 * nPointCount );
187 : }
188 : else
189 : {
190 90 : for( int i = 0; i < nPointCount; i++ )
191 : {
192 74 : memcpy( &(paoPoints[i].x), pabyData + 4 + 24 * i, 8 );
193 74 : memcpy( &(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8 );
194 74 : memcpy( padfZ + i, pabyData + 4 + 24 * i + 16, 8 );
195 : }
196 : }
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Byte swap if needed. */
200 : /* -------------------------------------------------------------------- */
201 223 : 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 223 : 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 179 : 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 179 : memcpy( pabyData, &nPointCount, 4 );
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Copy in the raw data. */
238 : /* -------------------------------------------------------------------- */
239 179 : if( b3D )
240 : {
241 20 : nWords = 3 * nPointCount;
242 311 : for( i = 0; i < nPointCount; i++ )
243 : {
244 291 : memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
245 291 : memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
246 291 : if( padfZ == NULL )
247 0 : memset( pabyData+4+i*24+16, 0, 8 );
248 : else
249 291 : memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
250 : }
251 : }
252 : else
253 : {
254 159 : nWords = 2 * nPointCount;
255 159 : memcpy( pabyData+4, paoPoints, 16 * nPointCount );
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Swap if needed. */
260 : /* -------------------------------------------------------------------- */
261 179 : 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 179 : return OGRERR_NONE;
275 : }
276 :
277 : /************************************************************************/
278 : /* _WkbSize() */
279 : /* */
280 : /* Helper method for OGRPolygon. NOT THE NORMAL WkbSize() METHOD! */
281 : /************************************************************************/
282 :
283 907 : int OGRLinearRing::_WkbSize( int b3D ) const
284 :
285 : {
286 907 : if( b3D )
287 96 : return 4 + 24 * nPointCount;
288 : else
289 811 : 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 279 : OGRGeometry *OGRLinearRing::clone() const
300 :
301 : {
302 : OGRLinearRing *poNewLinearRing;
303 :
304 279 : poNewLinearRing = new OGRLinearRing();
305 279 : poNewLinearRing->assignSpatialReference( getSpatialReference() );
306 :
307 279 : poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ );
308 :
309 279 : return poNewLinearRing;
310 : }
311 :
312 : /************************************************************************/
313 : /* isClockwise() */
314 : /************************************************************************/
315 :
316 : /**
317 : * \brief Returns TRUE if the ring has clockwise winding (or less than 2 points)
318 : *
319 : * @return TRUE if clockwise otherwise FALSE.
320 : */
321 :
322 283 : int OGRLinearRing::isClockwise() const
323 :
324 : {
325 283 : double dfSum = 0.0;
326 :
327 283 : if( nPointCount < 2 )
328 0 : return TRUE;
329 :
330 3961 : for( int iVert = 0; iVert < nPointCount-1; iVert++ )
331 : {
332 7356 : dfSum += paoPoints[iVert].x * paoPoints[iVert+1].y
333 7356 : - paoPoints[iVert].y * paoPoints[iVert+1].x;
334 : }
335 :
336 566 : dfSum += paoPoints[nPointCount-1].x * paoPoints[0].y
337 566 : - paoPoints[nPointCount-1].y * paoPoints[0].x;
338 :
339 283 : return dfSum < 0.0;
340 : }
341 :
342 : /************************************************************************/
343 : /* reverseWindingOrder() */
344 : /************************************************************************/
345 :
346 0 : void OGRLinearRing::reverseWindingOrder()
347 :
348 : {
349 0 : int pos = 0;
350 0 : OGRPoint tempPoint;
351 :
352 0 : for( int i = 0; i < nPointCount / 2; i++ )
353 : {
354 0 : getPoint( i, &tempPoint );
355 0 : pos = nPointCount - i - 1;
356 0 : setPoint( i, getX(pos), getY(pos), getZ(pos) );
357 0 : setPoint( pos, tempPoint.getX(), tempPoint.getY(), tempPoint.getZ() );
358 0 : }
359 0 : }
360 :
361 : /************************************************************************/
362 : /* closeRing() */
363 : /************************************************************************/
364 :
365 95 : void OGRLinearRing::closeRings()
366 :
367 : {
368 95 : if( nPointCount < 2 )
369 0 : return;
370 :
371 95 : if( getX(0) != getX(nPointCount-1)
372 : || getY(0) != getY(nPointCount-1)
373 : || getZ(0) != getZ(nPointCount-1) )
374 : {
375 : /* Avoid implicit change of coordinate dimensionality
376 : * if z=0.0 and dim=2
377 : */
378 21 : if( getCoordinateDimension() == 2 )
379 20 : addPoint( getX(0), getY(0) );
380 : else
381 1 : addPoint( getX(0), getY(0), getZ(0) );
382 :
383 : }
384 : }
385 :
386 : /************************************************************************/
387 : /* get_Area() */
388 : /************************************************************************/
389 :
390 : /**
391 : * \brief Compute area of ring.
392 : *
393 : * The area is computed according to Green's Theorem:
394 : *
395 : * Area is "Sum(x(i)*y(i+1) - x(i+1)*y(i))/2" for i = 0 to pointCount-1,
396 : * assuming the last point is a duplicate of the first.
397 : *
398 : * @return computed area.
399 : */
400 :
401 292 : double OGRLinearRing::get_Area() const
402 :
403 : {
404 292 : double dfAreaSum = 0.0;
405 : int i;
406 :
407 292 : if( nPointCount < 2 )
408 1 : return 0;
409 :
410 3998 : for( i = 0; i < nPointCount-1; i++ )
411 : {
412 7414 : dfAreaSum += 0.5 * ( paoPoints[i].x * paoPoints[i+1].y
413 7414 : - paoPoints[i+1].x * paoPoints[i].y );
414 : }
415 :
416 582 : dfAreaSum += 0.5 * ( paoPoints[nPointCount-1].x * paoPoints[0].y
417 582 : - paoPoints[0].x * paoPoints[nPointCount-1].y );
418 :
419 291 : return fabs(dfAreaSum);
420 : }
421 :
422 : /************************************************************************/
423 : /* isPointInRing() */
424 : /************************************************************************/
425 :
426 259 : OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint, int bTestEnvelope) const
427 : {
428 259 : if ( NULL == poPoint )
429 : {
430 0 : CPLDebug( "OGR", "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - passed point is NULL!" );
431 0 : return 0;
432 : }
433 :
434 259 : const int iNumPoints = getNumPoints();
435 :
436 : // Simple validation
437 259 : if ( iNumPoints < 4 )
438 0 : return 0;
439 :
440 259 : const double dfTestX = poPoint->getX();
441 259 : const double dfTestY = poPoint->getY();
442 :
443 : // Fast test if point is inside extent of the ring
444 259 : if (bTestEnvelope)
445 : {
446 99 : OGREnvelope extent;
447 99 : getEnvelope(&extent);
448 99 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
449 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
450 : {
451 0 : return 0;
452 : }
453 : }
454 :
455 : // For every point p in ring,
456 : // test if ray starting from given point crosses segment (p - 1, p)
457 259 : int iNumCrossings = 0;
458 :
459 78660 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
460 : {
461 78401 : const int iPointPrev = iPoint - 1;
462 :
463 78401 : const double x1 = getX(iPoint) - dfTestX;
464 78401 : const double y1 = getY(iPoint) - dfTestY;
465 :
466 78401 : const double x2 = getX(iPointPrev) - dfTestX;
467 78401 : const double y2 = getY(iPointPrev) - dfTestY;
468 :
469 78401 : if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) )
470 : {
471 : // Check if ray intersects with segment of the ring
472 600 : const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
473 600 : if ( 0.0 < dfIntersection )
474 : {
475 : // Count intersections
476 322 : iNumCrossings++;
477 : }
478 : }
479 : }
480 :
481 : // If iNumCrossings number is even, given point is outside the ring,
482 : // when the crossings number is odd, the point is inside the ring.
483 259 : return ( ( iNumCrossings % 2 ) == 1 ? 1 : 0 );
484 : }
485 :
486 : /************************************************************************/
487 : /* isPointOnRingBoundary() */
488 : /************************************************************************/
489 :
490 7 : OGRBoolean OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint, int bTestEnvelope) const
491 : {
492 7 : if ( NULL == poPoint )
493 : {
494 0 : CPLDebug( "OGR", "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* poPoint) - passed point is NULL!" );
495 0 : return 0;
496 : }
497 :
498 7 : const int iNumPoints = getNumPoints();
499 :
500 : // Simple validation
501 7 : if ( iNumPoints < 4 )
502 0 : return 0;
503 :
504 7 : const double dfTestX = poPoint->getX();
505 7 : const double dfTestY = poPoint->getY();
506 :
507 : // Fast test if point is inside extent of the ring
508 7 : OGREnvelope extent;
509 7 : getEnvelope(&extent);
510 7 : if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
511 : && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
512 : {
513 0 : return 0;
514 : }
515 :
516 13 : for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ )
517 : {
518 13 : const int iPointPrev = iPoint - 1;
519 :
520 13 : const double x1 = getX(iPoint) - dfTestX;
521 13 : const double y1 = getY(iPoint) - dfTestY;
522 :
523 13 : const double x2 = getX(iPointPrev) - dfTestX;
524 13 : const double y2 = getY(iPointPrev) - dfTestY;
525 :
526 : /* If iPoint and iPointPrev are the same, go on */
527 13 : if (x1 == x2 && y1 == y2)
528 : {
529 0 : continue;
530 : }
531 :
532 : /* If the point is on the segment, return immediatly */
533 : /* FIXME? If the test point is not exactly identical to one of */
534 : /* the vertices of the ring, but somewhere on a segment, there's */
535 : /* little chance that we get 0. So that should be tested against some epsilon */
536 13 : if ( x1 * y2 - x2 * y1 == 0 )
537 : {
538 7 : return 1;
539 : }
540 : }
541 :
542 0 : return 0;
543 : }
|