Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
qgsgeometry.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometry.cpp - Geometry (stored as Open Geospatial Consortium WKB)
3  -------------------------------------------------------------------
4 Date : 02 May 2005
5 Copyright : (C) 2005 by Brendan Morley
6 email : morb at ozemail dot com dot au
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 /* $Id$ */
16 
17 #include <limits>
18 #include <cstdarg>
19 #include <cstdio>
20 #include <cmath>
21 
22 #include "qgis.h"
23 #include "qgsgeometry.h"
24 #include "qgsapplication.h"
25 #include "qgslogger.h"
26 #include "qgspoint.h"
27 #include "qgsrectangle.h"
28 
29 #include "qgsmaplayerregistry.h"
30 #include "qgsvectorlayer.h"
31 #include "qgsproject.h"
32 
33 #define DEFAULT_QUADRANT_SEGMENTS 8
34 
35 #define CATCH_GEOS(r) \
36  catch (GEOSException &e) \
37  { \
38  Q_UNUSED(e); \
39  QgsDebugMsg("GEOS: " + QString( e.what() ) ); \
40  return r; \
41  }
42 
44 {
45  public:
46  GEOSException( QString theMsg )
47  {
48  if ( theMsg == "Unknown exception thrown" && lastMsg.isNull() )
49  {
50  msg = theMsg;
51  }
52  else
53  {
54  msg = theMsg;
55  lastMsg = msg;
56  }
57  }
58 
59  // copy constructor
61  {
62  *this = rhs;
63  }
64 
66  {
67  if ( lastMsg == msg )
68  lastMsg = QString::null;
69  }
70 
71  QString what()
72  {
73  return msg;
74  }
75 
76  private:
77  QString msg;
78  static QString lastMsg;
79 };
80 
82 
83 static void throwGEOSException( const char *fmt, ... )
84 {
85  va_list ap;
86  char buffer[1024];
87 
88  va_start( ap, fmt );
89  vsnprintf( buffer, sizeof buffer, fmt, ap );
90  va_end( ap );
91 
92  QgsDebugMsg( QString( "GEOS exception encountered: %1" ).arg( buffer ) );
93 
94  throw GEOSException( QString::fromUtf8( buffer ) );
95 }
96 
97 static void printGEOSNotice( const char *fmt, ... )
98 {
99 #if defined(QGISDEBUG)
100  va_list ap;
101  char buffer[1024];
102 
103  va_start( ap, fmt );
104  vsnprintf( buffer, sizeof buffer, fmt, ap );
105  va_end( ap );
106 
107  QgsDebugMsg( QString( "GEOS notice: %1" ).arg( QString::fromUtf8( buffer ) ) );
108 #endif
109 }
110 
111 class GEOSInit
112 {
113  public:
115  {
116  initGEOS( printGEOSNotice, throwGEOSException );
117  }
118 
120  {
121  finishGEOS();
122  }
123 };
124 
126 
127 
128 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR<3)
129 #define GEOSGeom_getCoordSeq(g) GEOSGeom_getCoordSeq( (GEOSGeometry *) g )
130 #define GEOSGetExteriorRing(g) GEOSGetExteriorRing( (GEOSGeometry *)g )
131 #define GEOSGetNumInteriorRings(g) GEOSGetNumInteriorRings( (GEOSGeometry *)g )
132 #define GEOSGetInteriorRingN(g,i) GEOSGetInteriorRingN( (GEOSGeometry *)g, i )
133 #define GEOSDisjoint(g0,g1) GEOSDisjoint( (GEOSGeometry *)g0, (GEOSGeometry*)g1 )
134 #define GEOSIntersection(g0,g1) GEOSIntersection( (GEOSGeometry*) g0, (GEOSGeometry*)g1 )
135 #define GEOSBuffer(g, d, s) GEOSBuffer( (GEOSGeometry*) g, d, s )
136 #define GEOSArea(g, a) GEOSArea( (GEOSGeometry*) g, a )
137 #define GEOSSimplify(g, t) GEOSSimplify( (GEOSGeometry*) g, t )
138 #define GEOSGetCentroid(g) GEOSGetCentroid( (GEOSGeometry*) g )
139 
140 #define GEOSCoordSeq_getSize(cs,n) GEOSCoordSeq_getSize( (GEOSCoordSequence *) cs, n )
141 #define GEOSCoordSeq_getX(cs,i,x) GEOSCoordSeq_getX( (GEOSCoordSequence *)cs, i, x )
142 #define GEOSCoordSeq_getY(cs,i,y) GEOSCoordSeq_getY( (GEOSCoordSequence *)cs, i, y )
143 
144 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms );
145 
146 static GEOSGeometry *cloneGeosGeom( const GEOSGeometry *geom )
147 {
148  // for GEOS < 3.0 we have own cloning function
149  // because when cloning multipart geometries they're copied into more general geometry collection instance
150  int type = GEOSGeomTypeId(( GEOSGeometry * ) geom );
151 
152  if ( type == GEOS_MULTIPOINT || type == GEOS_MULTILINESTRING || type == GEOS_MULTIPOLYGON )
153  {
154  QVector<GEOSGeometry *> geoms;
155 
156  try
157  {
158 
159  for ( int i = 0; i < GEOSGetNumGeometries(( GEOSGeometry * )geom ); ++i )
160  geoms << GEOSGeom_clone(( GEOSGeometry * ) GEOSGetGeometryN(( GEOSGeometry * ) geom, i ) );
161 
162  return createGeosCollection( type, geoms );
163  }
164  catch ( GEOSException &e )
165  {
166  Q_UNUSED( e );
167  for ( int i = 0; i < geoms.count(); i++ )
168  GEOSGeom_destroy( geoms[i] );
169 
170  return 0;
171  }
172  }
173  else
174  {
175  return GEOSGeom_clone(( GEOSGeometry * ) geom );
176  }
177 }
178 
179 #define GEOSGeom_clone(g) cloneGeosGeom(g)
180 #endif
181 
183  : mGeometry( 0 ),
184  mGeometrySize( 0 ),
185  mGeos( 0 ),
186  mDirtyWkb( false ),
187  mDirtyGeos( false )
188 {
189 }
190 
192  : mGeometry( 0 ),
193  mGeometrySize( rhs.mGeometrySize ),
194  mDirtyWkb( rhs.mDirtyWkb ),
195  mDirtyGeos( rhs.mDirtyGeos )
196 {
197  if ( mGeometrySize && rhs.mGeometry )
198  {
199  mGeometry = new unsigned char[mGeometrySize];
200  memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
201  }
202 
203  // deep-copy the GEOS Geometry if appropriate
204  if ( rhs.mGeos )
205  {
206  mGeos = GEOSGeom_clone( rhs.mGeos );
207  }
208  else
209  {
210  mGeos = 0;
211  }
212 }
213 
216 {
217  if ( mGeometry )
218  {
219  delete [] mGeometry;
220  }
221 
222  if ( mGeos )
223  {
224  GEOSGeom_destroy( mGeos );
225  }
226 }
227 
228 static unsigned int getNumGeosPoints( const GEOSGeometry *geom )
229 {
230  unsigned int n;
231  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( geom );
232  GEOSCoordSeq_getSize( cs, &n );
233  return n;
234 }
235 
236 static GEOSGeometry *createGeosPoint( const QgsPoint &point )
237 {
238  GEOSCoordSequence *coord = GEOSCoordSeq_create( 1, 2 );
239  GEOSCoordSeq_setX( coord, 0, point.x() );
240  GEOSCoordSeq_setY( coord, 0, point.y() );
241  return GEOSGeom_createPoint( coord );
242 }
243 
244 static GEOSCoordSequence *createGeosCoordSequence( const QgsPolyline& points )
245 {
246  GEOSCoordSequence *coord = 0;
247 
248  try
249  {
250  coord = GEOSCoordSeq_create( points.count(), 2 );
251  int i;
252  for ( i = 0; i < points.count(); i++ )
253  {
254  GEOSCoordSeq_setX( coord, i, points[i].x() );
255  GEOSCoordSeq_setY( coord, i, points[i].y() );
256  }
257  return coord;
258  }
259  catch ( GEOSException &e )
260  {
261  Q_UNUSED( e );
262  /*if ( coord )
263  GEOSCoordSeq_destroy( coord );*/
264  throw;
265  }
266 }
267 
268 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms )
269 {
270  GEOSGeometry **geomarr = new GEOSGeometry*[ geoms.size()];
271  if ( !geomarr )
272  return 0;
273 
274  for ( int i = 0; i < geoms.size(); i++ )
275  geomarr[i] = geoms[i];
276 
277  GEOSGeometry *geom = 0;
278 
279  try
280  {
281  geom = GEOSGeom_createCollection( typeId, geomarr, geoms.size() );
282  }
283  catch ( GEOSException &e )
284  {
285  Q_UNUSED( e );
286  }
287 
288  delete [] geomarr;
289 
290  return geom;
291 }
292 
293 static GEOSGeometry *createGeosLineString( const QgsPolyline& polyline )
294 {
295  GEOSCoordSequence *coord = 0;
296 
297  try
298  {
299  coord = createGeosCoordSequence( polyline );
300  return GEOSGeom_createLineString( coord );
301  }
302  catch ( GEOSException &e )
303  {
304  Q_UNUSED( e );
305  //MH: for strange reasons, geos3 crashes when removing the coordinate sequence
306  //if ( coord )
307  //GEOSCoordSeq_destroy( coord );
308  return 0;
309  }
310 }
311 
312 static GEOSGeometry *createGeosLinearRing( const QgsPolyline& polyline )
313 {
314  GEOSCoordSequence *coord = 0;
315 
316  if ( polyline.count() == 0 )
317  return 0;
318 
319  try
320  {
321  if ( polyline[0] != polyline[polyline.size()-1] )
322  {
323  // Ring not closed
324  QgsPolyline closed( polyline );
325  closed << closed[0];
326  coord = createGeosCoordSequence( closed );
327  }
328  else
329  {
330  coord = createGeosCoordSequence( polyline );
331  }
332 
333  return GEOSGeom_createLinearRing( coord );
334  }
335  catch ( GEOSException &e )
336  {
337  Q_UNUSED( e );
338  /* as MH has noticed ^, this crashes geos
339  if ( coord )
340  GEOSCoordSeq_destroy( coord );*/
341  return 0;
342  }
343 }
344 
345 static GEOSGeometry *createGeosPolygon( const QVector<GEOSGeometry*> &rings )
346 {
347  GEOSGeometry *shell = rings[0];
348  GEOSGeometry **holes = NULL;
349 
350  if ( rings.size() > 1 )
351  {
352  holes = new GEOSGeometry*[ rings.size()-1 ];
353  if ( !holes )
354  return 0;
355 
356  for ( int i = 0; i < rings.size() - 1; i++ )
357  holes[i] = rings[i+1];
358  }
359 
360  GEOSGeometry *geom = GEOSGeom_createPolygon( shell, holes, rings.size() - 1 );
361 
362  if ( holes )
363  delete [] holes;
364 
365  return geom;
366 }
367 
368 static GEOSGeometry *createGeosPolygon( GEOSGeometry *shell )
369 {
370  return createGeosPolygon( QVector<GEOSGeometry*>() << shell );
371 }
372 
373 static GEOSGeometry *createGeosPolygon( const QgsPolygon& polygon )
374 {
375  if ( polygon.count() == 0 )
376  return 0;
377 
378  QVector<GEOSGeometry *> geoms;
379 
380  try
381  {
382  for ( int i = 0; i < polygon.count(); i++ )
383  geoms << createGeosLinearRing( polygon[i] );
384 
385  return createGeosPolygon( geoms );
386  }
387  catch ( GEOSException &e )
388  {
389  Q_UNUSED( e );
390  for ( int i = 0; i < geoms.count(); i++ )
391  GEOSGeom_destroy( geoms[i] );
392  return 0;
393  }
394 }
395 
396 static QgsGeometry *fromGeosGeom( GEOSGeometry *geom )
397 {
398  if ( !geom )
399  return 0;
400 
401  QgsGeometry* g = new QgsGeometry;
402  g->fromGeos( geom );
403  return g;
404 }
405 
407 {
408  try
409  {
410 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR>=3)
411  GEOSWKTReader *reader = GEOSWKTReader_create();
412  QgsGeometry *g = fromGeosGeom( GEOSWKTReader_read( reader, wkt.toLocal8Bit().data() ) );
413  GEOSWKTReader_destroy( reader );
414  return g;
415 #else
416  return fromGeosGeom( GEOSGeomFromWKT( wkt.toLocal8Bit().data() ) );
417 #endif
418  }
419  catch ( GEOSException &e )
420  {
421  Q_UNUSED( e );
422  return 0;
423  }
424 }
425 
427 {
428  return fromGeosGeom( createGeosPoint( point ) );
429 }
430 
432 {
433  return fromGeosGeom( createGeosLineString( polyline ) );
434 }
435 
437 {
438  return fromGeosGeom( createGeosPolygon( polygon ) );
439 }
440 
442 {
443  QVector<GEOSGeometry *> geoms;
444 
445  try
446  {
447  for ( int i = 0; i < multipoint.size(); ++i )
448  {
449  geoms << createGeosPoint( multipoint[i] );
450  }
451 
452  return fromGeosGeom( createGeosCollection( GEOS_MULTIPOINT, geoms ) );
453  }
454  catch ( GEOSException &e )
455  {
456  Q_UNUSED( e );
457  for ( int i = 0; i < geoms.size(); ++i )
458  GEOSGeom_destroy( geoms[i] );
459 
460  return 0;
461  }
462 }
463 
465 {
466  QVector<GEOSGeometry *> geoms;
467 
468  try
469  {
470  for ( int i = 0; i < multiline.count(); i++ )
471  geoms << createGeosLineString( multiline[i] );
472 
473  return fromGeosGeom( createGeosCollection( GEOS_MULTILINESTRING, geoms ) );
474  }
475  catch ( GEOSException &e )
476  {
477  Q_UNUSED( e );
478  for ( int i = 0; i < geoms.count(); i++ )
479  GEOSGeom_destroy( geoms[i] );
480 
481  return 0;
482  }
483 }
484 
486 {
487  if ( multipoly.count() == 0 )
488  return 0;
489 
490  QVector<GEOSGeometry *> geoms;
491 
492  try
493  {
494  for ( int i = 0; i < multipoly.count(); i++ )
495  geoms << createGeosPolygon( multipoly[i] );
496 
497  return fromGeosGeom( createGeosCollection( GEOS_MULTIPOLYGON, geoms ) );
498  }
499  catch ( GEOSException &e )
500  {
501  Q_UNUSED( e );
502  for ( int i = 0; i < geoms.count(); i++ )
503  GEOSGeom_destroy( geoms[i] );
504 
505  return 0;
506  }
507 }
508 
510 {
511  QgsPolyline ring;
512  ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
513  ring.append( QgsPoint( rect.xMaximum(), rect.yMinimum() ) );
514  ring.append( QgsPoint( rect.xMaximum(), rect.yMaximum() ) );
515  ring.append( QgsPoint( rect.xMinimum(), rect.yMaximum() ) );
516  ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
517 
518  QgsPolygon polygon;
519  polygon.append( ring );
520 
521  return fromPolygon( polygon );
522 }
523 
525 {
526  if ( &rhs == this )
527  {
528  return *this;
529  }
530 
531  // remove old geometry if it exists
532  if ( mGeometry )
533  {
534  delete [] mGeometry;
535  mGeometry = 0;
536  }
537 
539 
540  // deep-copy the GEOS Geometry if appropriate
541  GEOSGeom_destroy( mGeos );
542  mGeos = rhs.mGeos ? GEOSGeom_clone( rhs.mGeos ) : 0;
543 
544  mDirtyGeos = rhs.mDirtyGeos;
545  mDirtyWkb = rhs.mDirtyWkb;
546 
547  if ( mGeometrySize && rhs.mGeometry )
548  {
549  mGeometry = new unsigned char[mGeometrySize];
550  memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
551  }
552 
553  return *this;
554 } // QgsGeometry::operator=( QgsGeometry const & rhs )
555 
556 
557 void QgsGeometry::fromWkb( unsigned char * wkb, size_t length )
558 {
559  // delete any existing WKB geometry before assigning new one
560  if ( mGeometry )
561  {
562  delete [] mGeometry;
563  mGeometry = 0;
564  }
565  if ( mGeos )
566  {
567  GEOSGeom_destroy( mGeos );
568  mGeos = 0;
569  }
570 
571  mGeometry = wkb;
573 
574  mDirtyWkb = false;
575  mDirtyGeos = true;
576 }
577 
578 unsigned char * QgsGeometry::asWkb()
579 {
580  if ( mDirtyWkb )
581  {
582  // convert from GEOS
583  exportGeosToWkb();
584  }
585 
586  return mGeometry;
587 }
588 
590 {
591  if ( mDirtyWkb )
592  {
593  exportGeosToWkb();
594  }
595 
596  return mGeometrySize;
597 }
598 
599 GEOSGeometry* QgsGeometry::asGeos()
600 {
601  if ( mDirtyGeos )
602  {
603  if ( !exportWkbToGeos() )
604  {
605  return 0;
606  }
607  }
608  return mGeos;
609 }
610 
611 
613 {
614  unsigned char *geom = asWkb(); // ensure that wkb representation exists
615  if ( geom )
616  {
617  unsigned int wkbType;
618  memcpy( &wkbType, ( geom + 1 ), sizeof( wkbType ) );
619  return ( QGis::WkbType ) wkbType;
620  }
621  else
622  {
623  return QGis::WKBUnknown;
624  }
625 }
626 
627 
629 {
630  if ( mDirtyWkb )
631  {
632  // convert from GEOS
633  exportGeosToWkb();
634  }
635 
636  QGis::WkbType type = wkbType();
637  if ( type == QGis::WKBPoint || type == QGis::WKBPoint25D ||
638  type == QGis::WKBMultiPoint || type == QGis::WKBMultiPoint25D )
639  return QGis::Point;
640  if ( type == QGis::WKBLineString || type == QGis::WKBLineString25D ||
642  return QGis::Line;
643  if ( type == QGis::WKBPolygon || type == QGis::WKBPolygon25D ||
645  return QGis::Polygon;
646 
647  return QGis::UnknownGeometry;
648 }
649 
651 {
652  if ( mDirtyWkb )
653  {
654  // convert from GEOS
655  exportGeosToWkb();
656  }
657 
658  QGis::WkbType type = wkbType();
659  if ( type == QGis::WKBMultiPoint ||
660  type == QGis::WKBMultiPoint25D ||
661  type == QGis::WKBMultiLineString ||
662  type == QGis::WKBMultiLineString25D ||
663  type == QGis::WKBMultiPolygon ||
664  type == QGis::WKBMultiPolygon25D )
665  return true;
666 
667  return false;
668 }
669 
670 
671 void QgsGeometry::fromGeos( GEOSGeometry* geos )
672 {
673  // TODO - make this more heap-friendly
674 
675  if ( mGeos )
676  {
677  GEOSGeom_destroy( mGeos );
678  mGeos = 0;
679  }
680  if ( mGeometry )
681  {
682  delete [] mGeometry;
683  mGeometry = 0;
684  }
685 
686  mGeos = geos;
687 
688  mDirtyWkb = true;
689  mDirtyGeos = false;
690 }
691 
692 QgsPoint QgsGeometry::closestVertex( const QgsPoint& point, int& atVertex, int& beforeVertex, int& afterVertex, double& sqrDist )
693 {
694  // TODO: implement with GEOS
695  if ( mDirtyWkb )
696  {
697  exportGeosToWkb();
698  }
699 
700  if ( !mGeometry )
701  {
702  QgsDebugMsg( "WKB geometry not available!" );
703  return QgsPoint( 0, 0 );
704  }
705 
706  int vertexnr = -1;
707  int vertexcounter = 0;
709  double actdist = std::numeric_limits<double>::max();
710  double x = 0;
711  double y = 0;
712  double *tempx, *tempy;
713  memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
714  beforeVertex = -1;
715  afterVertex = -1;
716  bool hasZValue = false;
717 
718  switch ( wkbType )
719  {
720  case QGis::WKBPoint25D:
721  case QGis::WKBPoint:
722  {
723  x = *(( double * )( mGeometry + 5 ) );
724  y = *(( double * )( mGeometry + 5 + sizeof( double ) ) );
725  actdist = point.sqrDist( x, y );
726  vertexnr = 0;
727  break;
728  }
730  hasZValue = true;
731 
732  // fall-through
733 
734  case QGis::WKBLineString:
735  {
736  unsigned char* ptr = mGeometry + 5;
737  int* npoints = ( int* )ptr;
738  ptr += sizeof( int );
739  for ( int index = 0; index < *npoints; ++index )
740  {
741  tempx = ( double* )ptr;
742  ptr += sizeof( double );
743  tempy = ( double* )ptr;
744  if ( point.sqrDist( *tempx, *tempy ) < actdist )
745  {
746  x = *tempx;
747  y = *tempy;
748  actdist = point.sqrDist( *tempx, *tempy );
749  vertexnr = index;
750  if ( index == 0 )//assign the rubber band indices
751  {
752  beforeVertex = -1;
753  }
754  else
755  {
756  beforeVertex = index - 1;
757  }
758  if ( index == ( *npoints - 1 ) )
759  {
760  afterVertex = -1;
761  }
762  else
763  {
764  afterVertex = index + 1;
765  }
766  }
767  ptr += sizeof( double );
768  if ( hasZValue ) //skip z-coordinate for 25D geometries
769  {
770  ptr += sizeof( double );
771  }
772  }
773  break;
774  }
775  case QGis::WKBPolygon25D:
776  hasZValue = true;
777  case QGis::WKBPolygon:
778  {
779  int* nrings = ( int* )( mGeometry + 5 );
780  int* npoints;
781  unsigned char* ptr = mGeometry + 9;
782  for ( int index = 0; index < *nrings; ++index )
783  {
784  npoints = ( int* )ptr;
785  ptr += sizeof( int );
786  for ( int index2 = 0; index2 < *npoints; ++index2 )
787  {
788  tempx = ( double* )ptr;
789  ptr += sizeof( double );
790  tempy = ( double* )ptr;
791  if ( point.sqrDist( *tempx, *tempy ) < actdist )
792  {
793  x = *tempx;
794  y = *tempy;
795  actdist = point.sqrDist( *tempx, *tempy );
796  vertexnr = vertexcounter;
797  //assign the rubber band indices
798  if ( index2 == 0 )
799  {
800  beforeVertex = vertexcounter + ( *npoints - 2 );
801  afterVertex = vertexcounter + 1;
802  }
803  else if ( index2 == ( *npoints - 1 ) )
804  {
805  beforeVertex = vertexcounter - 1;
806  afterVertex = vertexcounter - ( *npoints - 2 );
807  }
808  else
809  {
810  beforeVertex = vertexcounter - 1;
811  afterVertex = vertexcounter + 1;
812  }
813  }
814  ptr += sizeof( double );
815  if ( hasZValue ) //skip z-coordinate for 25D geometries
816  {
817  ptr += sizeof( double );
818  }
819  ++vertexcounter;
820  }
821  }
822  break;
823  }
825  hasZValue = true;
826  case QGis::WKBMultiPoint:
827  {
828  unsigned char* ptr = mGeometry + 5;
829  int* npoints = ( int* )ptr;
830  ptr += sizeof( int );
831  for ( int index = 0; index < *npoints; ++index )
832  {
833  ptr += ( 1 + sizeof( int ) ); //skip endian and point type
834  tempx = ( double* )ptr;
835  tempy = ( double* )( ptr + sizeof( double ) );
836  if ( point.sqrDist( *tempx, *tempy ) < actdist )
837  {
838  x = *tempx;
839  y = *tempy;
840  actdist = point.sqrDist( *tempx, *tempy );
841  vertexnr = index;
842  }
843  ptr += ( 2 * sizeof( double ) );
844  if ( hasZValue ) //skip z-coordinate for 25D geometries
845  {
846  ptr += sizeof( double );
847  }
848  }
849  break;
850  }
852  hasZValue = true;
854  {
855  unsigned char* ptr = mGeometry + 5;
856  int* nlines = ( int* )ptr;
857  int* npoints = 0;
858  ptr += sizeof( int );
859  for ( int index = 0; index < *nlines; ++index )
860  {
861  ptr += ( sizeof( int ) + 1 );
862  npoints = ( int* )ptr;
863  ptr += sizeof( int );
864  for ( int index2 = 0; index2 < *npoints; ++index2 )
865  {
866  tempx = ( double* )ptr;
867  ptr += sizeof( double );
868  tempy = ( double* )ptr;
869  ptr += sizeof( double );
870  if ( point.sqrDist( *tempx, *tempy ) < actdist )
871  {
872  x = *tempx;
873  y = *tempy;
874  actdist = point.sqrDist( *tempx, *tempy );
875  vertexnr = vertexcounter;
876 
877  if ( index2 == 0 )//assign the rubber band indices
878  {
879  beforeVertex = -1;
880  }
881  else
882  {
883  beforeVertex = vertexnr - 1;
884  }
885  if ( index2 == ( *npoints ) - 1 )
886  {
887  afterVertex = -1;
888  }
889  else
890  {
891  afterVertex = vertexnr + 1;
892  }
893  }
894  if ( hasZValue ) //skip z-coordinate for 25D geometries
895  {
896  ptr += sizeof( double );
897  }
898  ++vertexcounter;
899  }
900  }
901  break;
902  }
904  hasZValue = true;
906  {
907  unsigned char* ptr = mGeometry + 5;
908  int* npolys = ( int* )ptr;
909  int* nrings;
910  int* npoints;
911  ptr += sizeof( int );
912  for ( int index = 0; index < *npolys; ++index )
913  {
914  ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
915  nrings = ( int* )ptr;
916  ptr += sizeof( int );
917  for ( int index2 = 0; index2 < *nrings; ++index2 )
918  {
919  npoints = ( int* )ptr;
920  ptr += sizeof( int );
921  for ( int index3 = 0; index3 < *npoints; ++index3 )
922  {
923  tempx = ( double* )ptr;
924  ptr += sizeof( double );
925  tempy = ( double* )ptr;
926  if ( point.sqrDist( *tempx, *tempy ) < actdist )
927  {
928  x = *tempx;
929  y = *tempy;
930  actdist = point.sqrDist( *tempx, *tempy );
931  vertexnr = vertexcounter;
932 
933  //assign the rubber band indices
934  if ( index3 == 0 )
935  {
936  beforeVertex = vertexcounter + ( *npoints - 2 );
937  afterVertex = vertexcounter + 1;
938  }
939  else if ( index3 == ( *npoints - 1 ) )
940  {
941  beforeVertex = vertexcounter - 1;
942  afterVertex = vertexcounter - ( *npoints - 2 );
943  }
944  else
945  {
946  beforeVertex = vertexcounter - 1;
947  afterVertex = vertexcounter + 1;
948  }
949  }
950  ptr += sizeof( double );
951  if ( hasZValue ) //skip z-coordinate for 25D geometries
952  {
953  ptr += sizeof( double );
954  }
955  ++vertexcounter;
956  }
957  }
958  }
959  break;
960  }
961  default:
962  break;
963  }
964  sqrDist = actdist;
965  atVertex = vertexnr;
966  return QgsPoint( x, y );
967 }
968 
969 
970 void QgsGeometry::adjacentVertices( int atVertex, int& beforeVertex, int& afterVertex )
971 {
972  // TODO: implement with GEOS
973  if ( mDirtyWkb )
974  {
975  exportGeosToWkb();
976  }
977 
978  beforeVertex = -1;
979  afterVertex = -1;
980 
981  if ( !mGeometry )
982  {
983  QgsDebugMsg( "WKB geometry not available!" );
984  return;
985  }
986 
987  int vertexcounter = 0;
988 
990  bool hasZValue = false;
991 
992  memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
993 
994  switch ( wkbType )
995  {
996  case QGis::WKBPoint:
997  {
998  // NOOP - Points do not have adjacent verticies
999  break;
1000  }
1002  case QGis::WKBLineString:
1003  {
1004  unsigned char* ptr = mGeometry + 5;
1005  int* npoints = ( int* ) ptr;
1006 
1007  const int index = atVertex;
1008 
1009  // assign the rubber band indices
1010 
1011  if ( index == 0 )
1012  {
1013  beforeVertex = -1;
1014  }
1015  else
1016  {
1017  beforeVertex = index - 1;
1018  }
1019 
1020  if ( index == ( *npoints - 1 ) )
1021  {
1022  afterVertex = -1;
1023  }
1024  else
1025  {
1026  afterVertex = index + 1;
1027  }
1028 
1029  break;
1030  }
1031  case QGis::WKBPolygon25D:
1032  hasZValue = true;
1033  case QGis::WKBPolygon:
1034  {
1035  int* nrings = ( int* )( mGeometry + 5 );
1036  int* npoints;
1037  unsigned char* ptr = mGeometry + 9;
1038 
1039  // Walk through the POLYGON WKB
1040 
1041  for ( int index0 = 0; index0 < *nrings; ++index0 )
1042  {
1043  npoints = ( int* )ptr;
1044  ptr += sizeof( int );
1045 
1046  for ( int index1 = 0; index1 < *npoints; ++index1 )
1047  {
1048  ptr += sizeof( double );
1049  ptr += sizeof( double );
1050  if ( hasZValue ) //skip z-coordinate for 25D geometries
1051  {
1052  ptr += sizeof( double );
1053  }
1054  if ( vertexcounter == atVertex )
1055  {
1056  if ( index1 == 0 )
1057  {
1058  beforeVertex = vertexcounter + ( *npoints - 2 );
1059  afterVertex = vertexcounter + 1;
1060  }
1061  else if ( index1 == ( *npoints - 1 ) )
1062  {
1063  beforeVertex = vertexcounter - 1;
1064  afterVertex = vertexcounter - ( *npoints - 2 );
1065  }
1066  else
1067  {
1068  beforeVertex = vertexcounter - 1;
1069  afterVertex = vertexcounter + 1;
1070  }
1071  }
1072 
1073  ++vertexcounter;
1074  }
1075  }
1076  break;
1077  }
1079  case QGis::WKBMultiPoint:
1080  {
1081  // NOOP - Points do not have adjacent verticies
1082  break;
1083  }
1084 
1086  hasZValue = true;
1088  {
1089  unsigned char* ptr = mGeometry + 5;
1090  int* nlines = ( int* )ptr;
1091  int* npoints = 0;
1092  ptr += sizeof( int );
1093 
1094  for ( int index0 = 0; index0 < *nlines; ++index0 )
1095  {
1096  ptr += ( sizeof( int ) + 1 );
1097  npoints = ( int* )ptr;
1098  ptr += sizeof( int );
1099 
1100  for ( int index1 = 0; index1 < *npoints; ++index1 )
1101  {
1102  ptr += sizeof( double );
1103  ptr += sizeof( double );
1104  if ( hasZValue ) //skip z-coordinate for 25D geometries
1105  {
1106  ptr += sizeof( double );
1107  }
1108 
1109  if ( vertexcounter == atVertex )
1110  {
1111  // Found the vertex of the linestring we were looking for.
1112  if ( index1 == 0 )
1113  {
1114  beforeVertex = -1;
1115  }
1116  else
1117  {
1118  beforeVertex = vertexcounter - 1;
1119  }
1120  if ( index1 == ( *npoints ) - 1 )
1121  {
1122  afterVertex = -1;
1123  }
1124  else
1125  {
1126  afterVertex = vertexcounter + 1;
1127  }
1128  }
1129  ++vertexcounter;
1130  }
1131  }
1132  break;
1133  }
1135  hasZValue = true;
1136  case QGis::WKBMultiPolygon:
1137  {
1138  unsigned char* ptr = mGeometry + 5;
1139  int* npolys = ( int* )ptr;
1140  int* nrings;
1141  int* npoints;
1142  ptr += sizeof( int );
1143 
1144  for ( int index0 = 0; index0 < *npolys; ++index0 )
1145  {
1146  ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
1147  nrings = ( int* )ptr;
1148  ptr += sizeof( int );
1149 
1150  for ( int index1 = 0; index1 < *nrings; ++index1 )
1151  {
1152  npoints = ( int* )ptr;
1153  ptr += sizeof( int );
1154 
1155  for ( int index2 = 0; index2 < *npoints; ++index2 )
1156  {
1157  ptr += sizeof( double );
1158  ptr += sizeof( double );
1159  if ( hasZValue ) //skip z-coordinate for 25D geometries
1160  {
1161  ptr += sizeof( double );
1162  }
1163  if ( vertexcounter == atVertex )
1164  {
1165  // Found the vertex of the linear-ring of the polygon we were looking for.
1166  // assign the rubber band indices
1167 
1168  if ( index2 == 0 )
1169  {
1170  beforeVertex = vertexcounter + ( *npoints - 2 );
1171  afterVertex = vertexcounter + 1;
1172  }
1173  else if ( index2 == ( *npoints - 1 ) )
1174  {
1175  beforeVertex = vertexcounter - 1;
1176  afterVertex = vertexcounter - ( *npoints - 2 );
1177  }
1178  else
1179  {
1180  beforeVertex = vertexcounter - 1;
1181  afterVertex = vertexcounter + 1;
1182  }
1183  }
1184  ++vertexcounter;
1185  }
1186  }
1187  }
1188 
1189  break;
1190  }
1191 
1192  default:
1193  break;
1194  } // switch (wkbType)
1195 }
1196 
1197 
1198 
1199 bool QgsGeometry::insertVertex( double x, double y,
1200  int beforeVertex,
1201  const GEOSCoordSequence* old_sequence,
1202  GEOSCoordSequence** new_sequence )
1203 
1204 {
1205  // Bounds checking
1206  if ( beforeVertex < 0 )
1207  {
1208  *new_sequence = 0;
1209  return false;
1210  }
1211 
1212  unsigned int numPoints;
1213  GEOSCoordSeq_getSize( old_sequence, &numPoints );
1214 
1215  *new_sequence = GEOSCoordSeq_create( numPoints + 1, 2 );
1216  if ( !*new_sequence )
1217  return false;
1218 
1219  bool inserted = false;
1220  for ( unsigned int i = 0, j = 0; i < numPoints; i++, j++ )
1221  {
1222  // Do we insert the new vertex here?
1223  if ( beforeVertex == static_cast<int>( i ) )
1224  {
1225  GEOSCoordSeq_setX( *new_sequence, j, x );
1226  GEOSCoordSeq_setY( *new_sequence, j, y );
1227  j++;
1228  inserted = true;
1229  }
1230 
1231  double aX, aY;
1232  GEOSCoordSeq_getX( old_sequence, i, &aX );
1233  GEOSCoordSeq_getY( old_sequence, i, &aY );
1234 
1235  GEOSCoordSeq_setX( *new_sequence, j, aX );
1236  GEOSCoordSeq_setY( *new_sequence, j, aY );
1237  }
1238 
1239  if ( !inserted )
1240  {
1241  // The beforeVertex is greater than the actual number of vertices
1242  // in the geometry - append it.
1243  GEOSCoordSeq_setX( *new_sequence, numPoints, x );
1244  GEOSCoordSeq_setY( *new_sequence, numPoints, y );
1245  }
1246  // TODO: Check that the sequence is still simple, e.g. with GEOS_GEOM::Geometry->isSimple()
1247 
1248  return inserted;
1249 }
1250 
1251 bool QgsGeometry::moveVertex( double x, double y, int atVertex )
1252 {
1253  int vertexnr = atVertex;
1254 
1255  // TODO: implement with GEOS
1256  if ( mDirtyWkb )
1257  {
1258  exportGeosToWkb();
1259  }
1260 
1261  if ( !mGeometry )
1262  {
1263  QgsDebugMsg( "WKB geometry not available!" );
1264  return false;
1265  }
1266 
1268  bool hasZValue = false;
1269  unsigned char* ptr = mGeometry + 1;
1270  memcpy( &wkbType, ptr, sizeof( wkbType ) );
1271  ptr += sizeof( wkbType );
1272 
1273  switch ( wkbType )
1274  {
1275  case QGis::WKBPoint25D:
1276  case QGis::WKBPoint:
1277  {
1278  if ( vertexnr == 0 )
1279  {
1280  memcpy( ptr, &x, sizeof( double ) );
1281  ptr += sizeof( double );
1282  memcpy( ptr, &y, sizeof( double ) );
1283  mDirtyGeos = true;
1284  return true;
1285  }
1286  else
1287  {
1288  return false;
1289  }
1290  }
1292  hasZValue = true;
1293  case QGis::WKBMultiPoint:
1294  {
1295  int* nrPoints = ( int* )ptr;
1296  if ( vertexnr > *nrPoints || vertexnr < 0 )
1297  {
1298  return false;
1299  }
1300  ptr += sizeof( int );
1301  if ( hasZValue )
1302  {
1303  ptr += ( 3 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr;
1304  }
1305  else
1306  {
1307  ptr += ( 2 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr;
1308  }
1309  ptr += ( 1 + sizeof( int ) );
1310  memcpy( ptr, &x, sizeof( double ) );
1311  ptr += sizeof( double );
1312  memcpy( ptr, &y, sizeof( double ) );
1313  mDirtyGeos = true;
1314  return true;
1315  }
1317  hasZValue = true;
1318  case QGis::WKBLineString:
1319  {
1320  int* nrPoints = ( int* )ptr;
1321  if ( vertexnr > *nrPoints || vertexnr < 0 )
1322  {
1323  return false;
1324  }
1325  ptr += sizeof( int );
1326  ptr += 2 * sizeof( double ) * vertexnr;
1327  if ( hasZValue )
1328  {
1329  ptr += sizeof( double ) * vertexnr;
1330  }
1331  memcpy( ptr, &x, sizeof( double ) );
1332  ptr += sizeof( double );
1333  memcpy( ptr, &y, sizeof( double ) );
1334  mDirtyGeos = true;
1335  return true;
1336  }
1338  hasZValue = true;
1340  {
1341  int* nrLines = ( int* )ptr;
1342  ptr += sizeof( int );
1343  int* nrPoints = 0; //numer of points in a line
1344  int pointindex = 0;
1345  for ( int linenr = 0; linenr < *nrLines; ++linenr )
1346  {
1347  ptr += sizeof( int ) + 1;
1348  nrPoints = ( int* )ptr;
1349  ptr += sizeof( int );
1350  if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nrPoints ) )
1351  {
1352  ptr += ( vertexnr - pointindex ) * 2 * sizeof( double );
1353  if ( hasZValue )
1354  {
1355  ptr += ( vertexnr - pointindex ) * sizeof( double );
1356  }
1357  memcpy( ptr, &x, sizeof( double ) );
1358  memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
1359  mDirtyGeos = true;
1360  return true;
1361  }
1362  pointindex += ( *nrPoints );
1363  ptr += 2 * sizeof( double ) * ( *nrPoints );
1364  if ( hasZValue )
1365  {
1366  ptr += sizeof( double ) * ( *nrPoints );
1367  }
1368  }
1369  return false;
1370  }
1371  case QGis::WKBPolygon25D:
1372  hasZValue = true;
1373  case QGis::WKBPolygon:
1374  {
1375  int* nrRings = ( int* )ptr;
1376  ptr += sizeof( int );
1377  int* nrPoints = 0; //numer of points in a ring
1378  int pointindex = 0;
1379 
1380  for ( int ringnr = 0; ringnr < *nrRings; ++ringnr )
1381  {
1382  nrPoints = ( int* )ptr;
1383  ptr += sizeof( int );
1384  if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )//move two points
1385  {
1386  memcpy( ptr, &x, sizeof( double ) );
1387  memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
1388  if ( hasZValue )
1389  {
1390  memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
1391  }
1392  else
1393  {
1394  memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
1395  }
1396  if ( hasZValue )
1397  {
1398  memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
1399  }
1400  else
1401  {
1402  memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
1403  }
1404  mDirtyGeos = true;
1405  return true;
1406  }
1407  else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )//move only one point
1408  {
1409  ptr += 2 * sizeof( double ) * ( vertexnr - pointindex );
1410  if ( hasZValue )
1411  {
1412  ptr += sizeof( double ) * ( vertexnr - pointindex );
1413  }
1414  memcpy( ptr, &x, sizeof( double ) );
1415  ptr += sizeof( double );
1416  memcpy( ptr, &y, sizeof( double ) );
1417  mDirtyGeos = true;
1418  return true;
1419  }
1420  ptr += 2 * sizeof( double ) * ( *nrPoints );
1421  if ( hasZValue )
1422  {
1423  ptr += sizeof( double ) * ( *nrPoints );
1424  }
1425  pointindex += *nrPoints;
1426  }
1427  return false;
1428  }
1430  hasZValue = true;
1431  case QGis::WKBMultiPolygon:
1432  {
1433  int* nrPolygons = ( int* )ptr;
1434  ptr += sizeof( int );
1435  int* nrRings = 0; //number of rings in a polygon
1436  int* nrPoints = 0; //number of points in a ring
1437  int pointindex = 0;
1438 
1439  for ( int polynr = 0; polynr < *nrPolygons; ++polynr )
1440  {
1441  ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
1442  nrRings = ( int* )ptr;
1443  ptr += sizeof( int );
1444  for ( int ringnr = 0; ringnr < *nrRings; ++ringnr )
1445  {
1446  nrPoints = ( int* )ptr;
1447  ptr += sizeof( int );
1448  if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )//move two points
1449  {
1450  memcpy( ptr, &x, sizeof( double ) );
1451  memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
1452  if ( hasZValue )
1453  {
1454  memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
1455  }
1456  else
1457  {
1458  memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
1459  }
1460  if ( hasZValue )
1461  {
1462  memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
1463  }
1464  else
1465  {
1466  memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
1467  }
1468  mDirtyGeos = true;
1469  return true;
1470  }
1471  else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )//move only one point
1472  {
1473  ptr += 2 * sizeof( double ) * ( vertexnr - pointindex );
1474  if ( hasZValue )
1475  {
1476  ptr += sizeof( double ) * ( vertexnr - pointindex );
1477  }
1478  memcpy( ptr, &x, sizeof( double ) );
1479  ptr += sizeof( double );
1480  memcpy( ptr, &y, sizeof( double ) );
1481  mDirtyGeos = true;
1482  return true;
1483  }
1484  ptr += 2 * sizeof( double ) * ( *nrPoints );
1485  if ( hasZValue )
1486  {
1487  ptr += sizeof( double ) * ( *nrPoints );
1488  }
1489  pointindex += *nrPoints;
1490  }
1491  }
1492  return false;
1493  }
1494 
1495  default:
1496  return false;
1497  }
1498 }
1499 
1500 bool QgsGeometry::deleteVertex( int atVertex )
1501 {
1502  int vertexnr = atVertex;
1503  bool success = false;
1504 
1505  // TODO: implement with GEOS
1506  if ( mDirtyWkb )
1507  {
1508  exportGeosToWkb();
1509  }
1510 
1511  if ( !mGeometry )
1512  {
1513  QgsDebugMsg( "WKB geometry not available!" );
1514  return false;
1515  }
1516 
1517  //create a new geometry buffer for the modified geometry
1518  unsigned char* newbuffer;
1519  int pointindex = 0;
1521  bool hasZValue = false;
1522  unsigned char* ptr = mGeometry + 1;
1523  memcpy( &wkbType, ptr, sizeof( wkbType ) );
1524 
1525  switch ( wkbType )
1526  {
1527  case QGis::WKBPoint25D:
1529  case QGis::WKBPolygon25D:
1533  newbuffer = new unsigned char[mGeometrySize-3*sizeof( double )];
1534  break;
1535  default:
1536  newbuffer = new unsigned char[mGeometrySize-2*sizeof( double )];
1537  }
1538 
1539  memcpy( newbuffer, mGeometry, 1 + sizeof( int ) ); //endian and type are the same
1540 
1541  ptr += sizeof( wkbType );
1542  unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int );
1543 
1544  switch ( wkbType )
1545  {
1546  case QGis::WKBPoint25D:
1547  case QGis::WKBPoint:
1548  {
1549  break; //cannot remove the only point vertex
1550  }
1552  hasZValue = true;
1553  case QGis::WKBMultiPoint:
1554  {
1555  //todo
1556  break;
1557  }
1559  hasZValue = true;
1560  case QGis::WKBLineString:
1561  {
1562  int* nPoints = ( int* )ptr;
1563  if (( *nPoints ) < 3 || vertexnr > ( *nPoints ) - 1 || vertexnr < 0 ) //line needs at least 2 vertices
1564  {
1565  delete newbuffer;
1566  return false;
1567  }
1568  int newNPoints = ( *nPoints ) - 1; //new number of points
1569  memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
1570  ptr += sizeof( int );
1571  newBufferPtr += sizeof( int );
1572  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1573  {
1574  if ( vertexnr != pointindex )
1575  {
1576  memcpy( newBufferPtr, ptr, sizeof( double ) );
1577  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
1578  newBufferPtr += 2 * sizeof( double );
1579  if ( hasZValue )
1580  {
1581  newBufferPtr += sizeof( double );
1582  }
1583  }
1584  else
1585  {
1586  success = true;
1587  }
1588  ptr += 2 * sizeof( double );
1589  if ( hasZValue )
1590  {
1591  ptr += sizeof( double );
1592  }
1593  ++pointindex;
1594  }
1595  break;
1596  }
1598  hasZValue = true;
1600  {
1601  int* nLines = ( int* )ptr;
1602  memcpy( newBufferPtr, nLines, sizeof( int ) );
1603  newBufferPtr += sizeof( int );
1604  ptr += sizeof( int );
1605  int* nPoints = 0; //number of points in a line
1606  int pointindex = 0;
1607  for ( int linenr = 0; linenr < *nLines; ++linenr )
1608  {
1609  memcpy( newBufferPtr, ptr, sizeof( int ) + 1 );
1610  ptr += ( sizeof( int ) + 1 );
1611  newBufferPtr += ( sizeof( int ) + 1 );
1612  nPoints = ( int* )ptr;
1613  ptr += sizeof( int );
1614  int newNPoint;
1615 
1616  //find out if the vertex is in this line
1617  if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nPoints ) )
1618  {
1619  if ( *nPoints < 3 ) //line needs at least 2 vertices
1620  {
1621  delete newbuffer;
1622  return false;
1623  }
1624  newNPoint = ( *nPoints ) - 1;
1625  }
1626  else
1627  {
1628  newNPoint = *nPoints;
1629  }
1630  memcpy( newBufferPtr, &newNPoint, sizeof( int ) );
1631  newBufferPtr += sizeof( int );
1632 
1633  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1634  {
1635  if ( vertexnr != pointindex )
1636  {
1637  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
1638  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
1639  newBufferPtr += 2 * sizeof( double );
1640  if ( hasZValue )
1641  {
1642  newBufferPtr += sizeof( double );
1643  }
1644  }
1645  else
1646  {
1647  success = true;
1648  }
1649  ptr += 2 * sizeof( double );
1650  if ( hasZValue )
1651  {
1652  ptr += sizeof( double );
1653  }
1654  ++pointindex;
1655  }
1656  }
1657  break;
1658  }
1659  case QGis::WKBPolygon25D:
1660  hasZValue = true;
1661  case QGis::WKBPolygon:
1662  {
1663  int* nRings = ( int* )ptr;
1664  memcpy( newBufferPtr, nRings, sizeof( int ) );
1665  ptr += sizeof( int );
1666  newBufferPtr += sizeof( int );
1667  int* nPoints = 0; //number of points in a ring
1668  int pointindex = 0;
1669 
1670  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
1671  {
1672  nPoints = ( int* )ptr;
1673  ptr += sizeof( int );
1674  int newNPoints;
1675  if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )//vertex to delete is in this ring
1676  {
1677  if ( *nPoints < 5 ) //a ring has at least 3 points
1678  {
1679  delete newbuffer;
1680  return false;
1681  }
1682  newNPoints = *nPoints - 1;
1683  }
1684  else
1685  {
1686  newNPoints = *nPoints;
1687  }
1688  memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
1689  newBufferPtr += sizeof( int );
1690 
1691  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1692  {
1693  if ( vertexnr != pointindex )
1694  {
1695  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
1696  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
1697  newBufferPtr += 2 * sizeof( double );
1698  if ( hasZValue )
1699  {
1700  newBufferPtr += sizeof( double );
1701  }
1702  }
1703  else
1704  {
1705  if ( pointnr == 0 )//adjust the last point of the ring
1706  {
1707  memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) );
1708  memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) );
1709  }
1710  if ( pointnr == ( *nPoints ) - 1 )
1711  {
1712  memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) );
1713  memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) );
1714  }
1715  success = true;
1716  }
1717  ptr += 2 * sizeof( double );
1718  if ( hasZValue )
1719  {
1720  ptr += sizeof( double );
1721  }
1722  ++pointindex;
1723  }
1724  }
1725  break;
1726  }
1728  hasZValue = true;
1729  case QGis::WKBMultiPolygon:
1730  {
1731  int* nPolys = ( int* )ptr;
1732  memcpy( newBufferPtr, nPolys, sizeof( int ) );
1733  newBufferPtr += sizeof( int );
1734  ptr += sizeof( int );
1735  int* nRings = 0;
1736  int* nPoints = 0;
1737  int pointindex = 0;
1738 
1739  for ( int polynr = 0; polynr < *nPolys; ++polynr )
1740  {
1741  memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) );
1742  ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
1743  newBufferPtr += ( 1 + sizeof( int ) );
1744  nRings = ( int* )ptr;
1745  memcpy( newBufferPtr, nRings, sizeof( int ) );
1746  newBufferPtr += sizeof( int );
1747  ptr += sizeof( int );
1748  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
1749  {
1750  nPoints = ( int* )ptr;
1751  ptr += sizeof( int );
1752  int newNPoints;
1753  if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )//vertex to delete is in this ring
1754  {
1755  if ( *nPoints < 5 ) //a ring has at least 3 points
1756  {
1757  delete newbuffer;
1758  return false;
1759  }
1760  newNPoints = *nPoints - 1;
1761  }
1762  else
1763  {
1764  newNPoints = *nPoints;
1765  }
1766  memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
1767  newBufferPtr += sizeof( int );
1768 
1769  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1770  {
1771  if ( vertexnr != pointindex )
1772  {
1773  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
1774  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
1775  newBufferPtr += 2 * sizeof( double );
1776  if ( hasZValue )
1777  {
1778  newBufferPtr += sizeof( double );
1779  }
1780  }
1781  else
1782  {
1783  if ( pointnr == 0 )//adjust the last point of the ring
1784  {
1785  memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) );
1786  memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) );
1787  }
1788  if ( pointnr == ( *nPoints ) - 1 )
1789  {
1790  memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) );
1791  memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) );
1792  }
1793  success = true;
1794  }
1795  ptr += 2 * sizeof( double );
1796  if ( hasZValue )
1797  {
1798  ptr += sizeof( double );
1799  }
1800  ++pointindex;
1801  }
1802  }
1803  }
1804  break;
1805  }
1806  case QGis::WKBNoGeometry:
1807  case QGis::WKBUnknown:
1808  break;
1809  }
1810  if ( success )
1811  {
1812  delete [] mGeometry;
1813  mGeometry = newbuffer;
1814  mGeometrySize -= ( 2 * sizeof( double ) );
1815  if ( hasZValue )
1816  {
1817  mGeometrySize -= sizeof( double );
1818  }
1819  mDirtyGeos = true;
1820  return true;
1821  }
1822  else
1823  {
1824  delete[] newbuffer;
1825  return false;
1826  }
1827 }
1828 
1829 bool QgsGeometry::insertVertex( double x, double y, int beforeVertex )
1830 {
1831  int vertexnr = beforeVertex;
1832  bool success = false;
1833 
1834  // TODO: implement with GEOS
1835  if ( mDirtyWkb )
1836  {
1837  exportGeosToWkb();
1838  }
1839 
1840  if ( !mGeometry )
1841  {
1842  QgsDebugMsg( "WKB geometry not available!" );
1843  return false;
1844  }
1845 
1846  //create a new geometry buffer for the modified geometry
1847  unsigned char* newbuffer;
1848 
1849  int pointindex = 0;
1851  bool hasZValue = false;
1852 
1853  unsigned char* ptr = mGeometry + 1;
1854  memcpy( &wkbType, ptr, sizeof( wkbType ) );
1855 
1856  switch ( wkbType )
1857  {
1858  case QGis::WKBPoint25D:
1860  case QGis::WKBPolygon25D:
1864  newbuffer = new unsigned char[mGeometrySize+3*sizeof( double )];
1865  break;
1866  default:
1867  newbuffer = new unsigned char[mGeometrySize+2*sizeof( double )];
1868  }
1869  memcpy( newbuffer, mGeometry, 1 + sizeof( int ) ); //endian and type are the same
1870 
1871  ptr += sizeof( wkbType );
1872  unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int );
1873 
1874  switch ( wkbType )
1875  {
1876  case QGis::WKBPoint25D:
1877  case QGis::WKBPoint://cannot insert a vertex before another one on point types
1878  {
1879  delete newbuffer;
1880  return false;
1881  }
1883  hasZValue = true;
1884  case QGis::WKBMultiPoint:
1885  {
1886  //todo
1887  break;
1888  }
1890  hasZValue = true;
1891  case QGis::WKBLineString:
1892  {
1893  int* nPoints = ( int* )ptr;
1894  if ( vertexnr > *nPoints || vertexnr < 0 )
1895  {
1896  break;
1897  }
1898  int newNPoints = ( *nPoints ) + 1;
1899  memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
1900  newBufferPtr += sizeof( int );
1901  ptr += sizeof( int );
1902 
1903  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1904  {
1905  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
1906  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//x
1907  ptr += 2 * sizeof( double );
1908  if ( hasZValue )
1909  {
1910  ptr += sizeof( double );
1911  }
1912  newBufferPtr += 2 * sizeof( double );
1913  if ( hasZValue )
1914  {
1915  newBufferPtr += sizeof( double );
1916  }
1917  ++pointindex;
1918  if ( pointindex == vertexnr )
1919  {
1920  memcpy( newBufferPtr, &x, sizeof( double ) );
1921  memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
1922  newBufferPtr += 2 * sizeof( double );
1923  if ( hasZValue )
1924  {
1925  newBufferPtr += sizeof( double );
1926  }
1927  success = true;
1928  }
1929  }
1930  break;
1931  }
1933  hasZValue = true;
1935  {
1936  int* nLines = ( int* )ptr;
1937  int* nPoints = 0; //number of points in a line
1938  ptr += sizeof( int );
1939  memcpy( newBufferPtr, nLines, sizeof( int ) );
1940  newBufferPtr += sizeof( int );
1941  int pointindex = 0;
1942 
1943  for ( int linenr = 0; linenr < *nLines; ++linenr )
1944  {
1945  memcpy( newBufferPtr, ptr, sizeof( int ) + 1 );
1946  ptr += ( sizeof( int ) + 1 );
1947  newBufferPtr += ( sizeof( int ) + 1 );
1948  nPoints = ( int* )ptr;
1949  int newNPoints;
1950  if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring
1951  {
1952  newNPoints = ( *nPoints ) + 1;
1953  }
1954  else
1955  {
1956  newNPoints = *nPoints;
1957  }
1958  memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
1959  newBufferPtr += sizeof( int );
1960  ptr += sizeof( int );
1961 
1962  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
1963  {
1964  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
1965  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
1966  ptr += 2 * sizeof( double );
1967  newBufferPtr += 2 * sizeof( double );
1968  if ( hasZValue )
1969  {
1970  ptr += sizeof( double );
1971  newBufferPtr += sizeof( double );
1972  }
1973  ++pointindex;
1974  if ( pointindex == vertexnr )
1975  {
1976  memcpy( newBufferPtr, &x, sizeof( double ) );
1977  memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
1978  newBufferPtr += 2 * sizeof( double );
1979  if ( hasZValue )
1980  {
1981  newBufferPtr += sizeof( double );
1982  }
1983  success = true;
1984  }
1985  }
1986  }
1987  break;
1988  }
1989  case QGis::WKBPolygon25D:
1990  hasZValue = true;
1991  case QGis::WKBPolygon:
1992  {
1993  int* nRings = ( int* )ptr;
1994  int* nPoints = 0; //number of points in a ring
1995  ptr += sizeof( int );
1996  memcpy( newBufferPtr, nRings, sizeof( int ) );
1997  newBufferPtr += sizeof( int );
1998  int pointindex = 0;
1999 
2000  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
2001  {
2002  nPoints = ( int* )ptr;
2003  int newNPoints;
2004  if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring
2005  {
2006  newNPoints = ( *nPoints ) + 1;
2007  }
2008  else
2009  {
2010  newNPoints = *nPoints;
2011  }
2012  memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
2013  newBufferPtr += sizeof( int );
2014  ptr += sizeof( int );
2015 
2016  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2017  {
2018  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
2019  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
2020  ptr += 2 * sizeof( double );
2021  newBufferPtr += 2 * sizeof( double );
2022  if ( hasZValue )
2023  {
2024  ptr += sizeof( double );
2025  newBufferPtr += sizeof( double );
2026  }
2027  ++pointindex;
2028  if ( pointindex == vertexnr )
2029  {
2030  memcpy( newBufferPtr, &x, sizeof( double ) );
2031  memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
2032  newBufferPtr += 2 * sizeof( double );
2033  if ( hasZValue )
2034  {
2035  newBufferPtr += sizeof( double );
2036  }
2037  success = true;
2038  }
2039  }
2040  }
2041  break;
2042  }
2044  hasZValue = true;
2045  case QGis::WKBMultiPolygon:
2046  {
2047  int* nPolys = ( int* )ptr;
2048  int* nRings = 0; //number of rings in a polygon
2049  int* nPoints = 0; //number of points in a ring
2050  memcpy( newBufferPtr, nPolys, sizeof( int ) );
2051  ptr += sizeof( int );
2052  newBufferPtr += sizeof( int );
2053  int pointindex = 0;
2054 
2055  for ( int polynr = 0; polynr < *nPolys; ++polynr )
2056  {
2057  memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) );
2058  ptr += ( 1 + sizeof( int ) );//skip endian and polygon type
2059  newBufferPtr += ( 1 + sizeof( int ) );
2060  nRings = ( int* )ptr;
2061  ptr += sizeof( int );
2062  memcpy( newBufferPtr, nRings, sizeof( int ) );
2063  newBufferPtr += sizeof( int );
2064 
2065  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
2066  {
2067  nPoints = ( int* )ptr;
2068  int newNPoints;
2069  if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring
2070  {
2071  newNPoints = ( *nPoints ) + 1;
2072  }
2073  else
2074  {
2075  newNPoints = *nPoints;
2076  }
2077  memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
2078  newBufferPtr += sizeof( int );
2079  ptr += sizeof( int );
2080 
2081  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2082  {
2083  memcpy( newBufferPtr, ptr, sizeof( double ) );//x
2084  memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y
2085  ptr += 2 * sizeof( double );
2086  newBufferPtr += 2 * sizeof( double );
2087  if ( hasZValue )
2088  {
2089  ptr += sizeof( double );
2090  newBufferPtr += sizeof( double );
2091  }
2092  ++pointindex;
2093  if ( pointindex == vertexnr )
2094  {
2095  memcpy( newBufferPtr, &x, sizeof( double ) );
2096  memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
2097  newBufferPtr += 2 * sizeof( double );
2098  if ( hasZValue )
2099  {
2100  newBufferPtr += sizeof( double );
2101  }
2102  success = true;
2103  }
2104  }
2105  }
2106 
2107  }
2108  break;
2109  }
2110  case QGis::WKBNoGeometry:
2111  case QGis::WKBUnknown:
2112  break;
2113  }
2114 
2115  if ( success )
2116  {
2117  delete [] mGeometry;
2118  mGeometry = newbuffer;
2119  mGeometrySize += 2 * sizeof( double );
2120  if ( hasZValue )
2121  {
2122  mGeometrySize += sizeof( double );
2123  }
2124  mDirtyGeos = true;
2125  return true;
2126  }
2127  else
2128  {
2129  delete newbuffer;
2130  return false;
2131  }
2132 }
2133 
2135 {
2136  double x, y;
2137 
2138  if ( mDirtyWkb )
2139  {
2140  exportGeosToWkb();
2141  }
2142 
2143  if ( !mGeometry )
2144  {
2145  QgsDebugMsg( "WKB geometry not available!" );
2146  return QgsPoint( 0, 0 );
2147  }
2148 
2150  bool hasZValue = false;
2151  unsigned char* ptr;
2152 
2153  memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
2154  switch ( wkbType )
2155  {
2156  case QGis::WKBPoint25D:
2157  case QGis::WKBPoint:
2158  {
2159  if ( atVertex == 0 )
2160  {
2161  ptr = mGeometry + 1 + sizeof( int );
2162  memcpy( &x, ptr, sizeof( double ) );
2163  ptr += sizeof( double );
2164  memcpy( &y, ptr, sizeof( double ) );
2165  return QgsPoint( x, y );
2166  }
2167  else
2168  {
2169  return QgsPoint( 0, 0 );
2170  }
2171  }
2173  hasZValue = true;
2174  case QGis::WKBLineString:
2175  {
2176  int *nPoints;
2177  // get number of points in the line
2178  ptr = mGeometry + 1 + sizeof( int ); // now at mGeometry.numPoints
2179  nPoints = ( int * ) ptr;
2180 
2181  // return error if underflow
2182  if ( 0 > atVertex || *nPoints <= atVertex )
2183  {
2184  return QgsPoint( 0, 0 );
2185  }
2186 
2187  // copy the vertex coordinates
2188  if ( hasZValue )
2189  {
2190  ptr = mGeometry + 9 + ( atVertex * 3 * sizeof( double ) );
2191  }
2192  else
2193  {
2194  ptr = mGeometry + 9 + ( atVertex * 2 * sizeof( double ) );
2195  }
2196  memcpy( &x, ptr, sizeof( double ) );
2197  ptr += sizeof( double );
2198  memcpy( &y, ptr, sizeof( double ) );
2199  return QgsPoint( x, y );
2200  }
2201  case QGis::WKBPolygon25D:
2202  hasZValue = true;
2203  case QGis::WKBPolygon:
2204  {
2205  int *nRings;
2206  int *nPoints = 0;
2207  ptr = mGeometry + 1 + sizeof( int );
2208  nRings = ( int* )ptr;
2209  ptr += sizeof( int );
2210  int pointindex = 0;
2211  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
2212  {
2213  nPoints = ( int* )ptr;
2214  ptr += sizeof( int );
2215  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2216  {
2217  if ( pointindex == atVertex )
2218  {
2219  memcpy( &x, ptr, sizeof( double ) );
2220  ptr += sizeof( double );
2221  memcpy( &y, ptr, sizeof( double ) );
2222  return QgsPoint( x, y );
2223  }
2224  ptr += 2 * sizeof( double );
2225  if ( hasZValue )
2226  {
2227  ptr += sizeof( double );
2228  }
2229  ++pointindex;
2230  }
2231  }
2232  return QgsPoint( 0, 0 );
2233  }
2235  hasZValue = true;
2236  case QGis::WKBMultiPoint:
2237  {
2238  ptr = mGeometry + 1 + sizeof( int );
2239  int* nPoints = ( int* )ptr;
2240  if ( atVertex < 0 || atVertex >= *nPoints )
2241  {
2242  return QgsPoint( 0, 0 );
2243  }
2244  if ( hasZValue )
2245  {
2246  ptr += atVertex * ( 3 * sizeof( double ) + 1 + sizeof( int ) );
2247  }
2248  else
2249  {
2250  ptr += atVertex * ( 2 * sizeof( double ) + 1 + sizeof( int ) );
2251  }
2252  ptr += 1 + sizeof( int );
2253  memcpy( &x, ptr, sizeof( double ) );
2254  ptr += sizeof( double );
2255  memcpy( &y, ptr, sizeof( double ) );
2256  return QgsPoint( x, y );
2257  }
2259  hasZValue = true;
2261  {
2262  ptr = mGeometry + 1 + sizeof( int );
2263  int* nLines = ( int* )ptr;
2264  int* nPoints = 0; //number of points in a line
2265  int pointindex = 0; //global point counter
2266  ptr += sizeof( int );
2267  for ( int linenr = 0; linenr < *nLines; ++linenr )
2268  {
2269  ptr += sizeof( int ) + 1;
2270  nPoints = ( int* )ptr;
2271  ptr += sizeof( int );
2272  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2273  {
2274  if ( pointindex == atVertex )
2275  {
2276  memcpy( &x, ptr, sizeof( double ) );
2277  ptr += sizeof( double );
2278  memcpy( &y, ptr, sizeof( double ) );
2279  return QgsPoint( x, y );
2280  }
2281  ptr += 2 * sizeof( double );
2282  if ( hasZValue )
2283  {
2284  ptr += sizeof( double );
2285  }
2286  ++pointindex;
2287  }
2288  }
2289  return QgsPoint( 0, 0 );
2290  }
2292  hasZValue = true;
2293  case QGis::WKBMultiPolygon:
2294  {
2295  ptr = mGeometry + 1 + sizeof( int );
2296  int* nRings = 0;//number of rings in a polygon
2297  int* nPoints = 0;//number of points in a ring
2298  int pointindex = 0; //global point counter
2299  int* nPolygons = ( int* )ptr;
2300  ptr += sizeof( int );
2301  for ( int polynr = 0; polynr < *nPolygons; ++polynr )
2302  {
2303  ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
2304  nRings = ( int* )ptr;
2305  ptr += sizeof( int );
2306  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
2307  {
2308  nPoints = ( int* )ptr;
2309  ptr += sizeof( int );
2310  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2311  {
2312  if ( pointindex == atVertex )
2313  {
2314  memcpy( &x, ptr, sizeof( double ) );
2315  ptr += sizeof( double );
2316  memcpy( &y, ptr, sizeof( double ) );
2317  return QgsPoint( x, y );
2318  }
2319  ++pointindex;
2320  ptr += 2 * sizeof( double );
2321  if ( hasZValue )
2322  {
2323  ptr += sizeof( double );
2324  }
2325  }
2326  }
2327  }
2328  return QgsPoint( 0, 0 );
2329  }
2330  default:
2331  QgsDebugMsg( "error: mGeometry type not recognized" );
2332  return QgsPoint( 0, 0 );
2333  }
2334 }
2335 
2336 
2337 double QgsGeometry::sqrDistToVertexAt( QgsPoint& point, int atVertex )
2338 {
2339  QgsPoint pnt = vertexAt( atVertex );
2340  if ( pnt != QgsPoint( 0, 0 ) )
2341  {
2342  QgsDebugMsg( "Exiting with distance to " + pnt.toString() );
2343  return point.sqrDist( pnt );
2344  }
2345  else
2346  {
2347  QgsDebugMsg( "Exiting with std::numeric_limits<double>::max()." );
2348  // probably safest to bail out with a very large number
2350  }
2351 }
2352 
2353 
2354 double QgsGeometry::closestVertexWithContext( const QgsPoint& point, int& atVertex )
2355 {
2356  double sqrDist = std::numeric_limits<double>::max();
2357 
2358  try
2359  {
2360  // Initialise some stuff
2361  int closestVertexIndex = 0;
2362 
2363  // set up the GEOS geometry
2364  if ( !exportWkbToGeos() )
2365  return -1;
2366 
2367  const GEOSGeometry *g = GEOSGetExteriorRing( mGeos );
2368  if ( !g )
2369  return -1;
2370 
2371  const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( g );
2372 
2373  unsigned int n;
2374  GEOSCoordSeq_getSize( sequence, &n );
2375 
2376  for ( unsigned int i = 0; i < n; i++ )
2377  {
2378  double x, y;
2379  GEOSCoordSeq_getX( sequence, i, &x );
2380  GEOSCoordSeq_getY( sequence, i, &y );
2381 
2382  double testDist = point.sqrDist( x, y );
2383  if ( testDist < sqrDist )
2384  {
2385  closestVertexIndex = i;
2386  sqrDist = testDist;
2387  }
2388  }
2389 
2390  atVertex = closestVertexIndex;
2391  }
2392  catch ( GEOSException &e )
2393  {
2394  Q_UNUSED( e );
2395  return -1;
2396  }
2397 
2398  return sqrDist;
2399 }
2400 
2401 
2403  const QgsPoint& point,
2404  QgsPoint& minDistPoint,
2405  int& beforeVertex )
2406 {
2407  QgsPoint distPoint;
2408 
2410  bool hasZValue = false;
2411  double *thisx = NULL;
2412  double *thisy = NULL;
2413  double *prevx = NULL;
2414  double *prevy = NULL;
2415  double testdist;
2416  int closestSegmentIndex = 0;
2417 
2418  // Initialise some stuff
2419  double sqrDist = std::numeric_limits<double>::max();
2420 
2421  // TODO: implement with GEOS
2422  if ( mDirtyWkb ) //convert latest geos to mGeometry
2423  {
2424  exportGeosToWkb();
2425  }
2426 
2427  if ( !mGeometry )
2428  {
2429  QgsDebugMsg( "WKB geometry not available!" );
2430  return -1;
2431  }
2432 
2433  memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
2434 
2435  switch ( wkbType )
2436  {
2437  case QGis::WKBPoint25D:
2438  case QGis::WKBPoint:
2440  case QGis::WKBMultiPoint:
2441  {
2442  // Points have no lines
2443  return -1;
2444  }
2446  hasZValue = true;
2447  case QGis::WKBLineString:
2448  {
2449  unsigned char* ptr = mGeometry + 1 + sizeof( int );
2450  int* npoints = ( int* ) ptr;
2451  ptr += sizeof( int );
2452  for ( int index = 0; index < *npoints; ++index )
2453  {
2454  if ( index > 0 )
2455  {
2456  prevx = thisx;
2457  prevy = thisy;
2458  }
2459  thisx = ( double* ) ptr;
2460  ptr += sizeof( double );
2461  thisy = ( double* ) ptr;
2462 
2463  if ( index > 0 )
2464  {
2465  if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
2466  {
2467  closestSegmentIndex = index;
2468  sqrDist = testdist;
2469  minDistPoint = distPoint;
2470  }
2471  }
2472  ptr += sizeof( double );
2473  if ( hasZValue )
2474  {
2475  ptr += sizeof( double );
2476  }
2477  }
2478  beforeVertex = closestSegmentIndex;
2479  break;
2480  }
2482  hasZValue = true;
2484  {
2485  unsigned char* ptr = mGeometry + 1 + sizeof( int );
2486  int* nLines = ( int* )ptr;
2487  ptr += sizeof( int );
2488  int* nPoints = 0; //number of points in a line
2489  int pointindex = 0;//global pointindex
2490  for ( int linenr = 0; linenr < *nLines; ++linenr )
2491  {
2492  ptr += sizeof( int ) + 1;
2493  nPoints = ( int* )ptr;
2494  ptr += sizeof( int );
2495  prevx = 0;
2496  prevy = 0;
2497  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2498  {
2499  thisx = ( double* ) ptr;
2500  ptr += sizeof( double );
2501  thisy = ( double* ) ptr;
2502  ptr += sizeof( double );
2503  if ( hasZValue )
2504  {
2505  ptr += sizeof( double );
2506  }
2507  if ( prevx && prevy )
2508  {
2509  if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
2510  {
2511  closestSegmentIndex = pointindex;
2512  sqrDist = testdist;
2513  minDistPoint = distPoint;
2514  }
2515  }
2516  prevx = thisx;
2517  prevy = thisy;
2518  ++pointindex;
2519  }
2520  }
2521  beforeVertex = closestSegmentIndex;
2522  break;
2523  }
2524  case QGis::WKBPolygon25D:
2525  hasZValue = true;
2526  case QGis::WKBPolygon:
2527  {
2528  int index = 0;
2529  unsigned char* ptr = mGeometry + 1 + sizeof( int );
2530  int* nrings = ( int* )ptr;
2531  int* npoints = 0; //number of points in a ring
2532  ptr += sizeof( int );
2533  for ( int ringnr = 0; ringnr < *nrings; ++ringnr )//loop over rings
2534  {
2535  npoints = ( int* )ptr;
2536  ptr += sizeof( int );
2537  prevx = 0;
2538  prevy = 0;
2539  for ( int pointnr = 0; pointnr < *npoints; ++pointnr )//loop over points in a ring
2540  {
2541  thisx = ( double* )ptr;
2542  ptr += sizeof( double );
2543  thisy = ( double* )ptr;
2544  ptr += sizeof( double );
2545  if ( hasZValue )
2546  {
2547  ptr += sizeof( double );
2548  }
2549  if ( prevx && prevy )
2550  {
2551  if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
2552  {
2553  closestSegmentIndex = index;
2554  sqrDist = testdist;
2555  minDistPoint = distPoint;
2556  }
2557  }
2558  prevx = thisx;
2559  prevy = thisy;
2560  ++index;
2561  }
2562  }
2563  beforeVertex = closestSegmentIndex;
2564  break;
2565  }
2567  hasZValue = true;
2568  case QGis::WKBMultiPolygon:
2569  {
2570  unsigned char* ptr = mGeometry + 1 + sizeof( int );
2571  int* nRings = 0;
2572  int* nPoints = 0;
2573  int pointindex = 0;
2574  int* nPolygons = ( int* )ptr;
2575  ptr += sizeof( int );
2576  for ( int polynr = 0; polynr < *nPolygons; ++polynr )
2577  {
2578  ptr += ( 1 + sizeof( int ) );
2579  nRings = ( int* )ptr;
2580  ptr += sizeof( int );
2581  for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
2582  {
2583  nPoints = ( int* )ptr;
2584  ptr += sizeof( int );
2585  prevx = 0;
2586  prevy = 0;
2587  for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
2588  {
2589  thisx = ( double* )ptr;
2590  ptr += sizeof( double );
2591  thisy = ( double* )ptr;
2592  ptr += sizeof( double );
2593  if ( hasZValue )
2594  {
2595  ptr += sizeof( double );
2596  }
2597  if ( prevx && prevy )
2598  {
2599  if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
2600  {
2601  closestSegmentIndex = pointindex;
2602  sqrDist = testdist;
2603  minDistPoint = distPoint;
2604  }
2605  }
2606  prevx = thisx;
2607  prevy = thisy;
2608  ++pointindex;
2609  }
2610  }
2611  }
2612  beforeVertex = closestSegmentIndex;
2613  break;
2614  }
2615  case QGis::WKBUnknown:
2616  default:
2617  return -1;
2618  break;
2619  } // switch (wkbType)
2620 
2621 
2622  QgsDebugMsg( "Exiting with nearest point " + point.toString() +
2623  ", dist: " + QString::number( sqrDist ) + "." );
2624 
2625  return sqrDist;
2626 }
2627 
2628 int QgsGeometry::addRing( const QList<QgsPoint>& ring )
2629 {
2630  //bail out if this geometry is not polygon/multipolygon
2631  if ( type() != QGis::Polygon )
2632  return 1;
2633 
2634  //test for invalid geometries
2635  if ( ring.size() < 4 )
2636  return 3;
2637 
2638  //ring must be closed
2639  if ( ring.first() != ring.last() )
2640  return 2;
2641 
2642  //create geos geometry from wkb if not already there
2643  if ( !mGeos || mDirtyGeos )
2644  if ( !exportWkbToGeos() )
2645  return 6;
2646 
2647  int type = GEOSGeomTypeId( mGeos );
2648 
2649  //Fill GEOS Polygons of the feature into list
2650  QVector<const GEOSGeometry*> polygonList;
2651 
2652  if ( wkbType() == QGis::WKBPolygon )
2653  {
2654  if ( type != GEOS_POLYGON )
2655  return 1;
2656 
2657  polygonList << mGeos;
2658  }
2659  else if ( wkbType() == QGis::WKBMultiPolygon )
2660  {
2661  if ( type != GEOS_MULTIPOLYGON )
2662  return 1;
2663 
2664  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); ++i )
2665  polygonList << GEOSGetGeometryN( mGeos, i );
2666  }
2667 
2668  //create new ring
2669  GEOSGeometry *newRing = 0;
2670  GEOSGeometry *newRingPolygon = 0;
2671 
2672  try
2673  {
2674  newRing = createGeosLinearRing( ring.toVector() );
2675  if ( !GEOSisValid( newRing ) )
2676  {
2677  throwGEOSException( "ring is invalid" );
2678  }
2679 
2680  newRingPolygon = createGeosPolygon( newRing );
2681  if ( !GEOSisValid( newRingPolygon ) )
2682  {
2683  throwGEOSException( "ring is invalid" );
2684  }
2685  }
2686  catch ( GEOSException &e )
2687  {
2688  Q_UNUSED( e );
2689  if ( newRingPolygon )
2690  GEOSGeom_destroy( newRingPolygon );
2691  else if ( newRing )
2692  GEOSGeom_destroy( newRing );
2693 
2694  return 3;
2695  }
2696 
2697  QVector<GEOSGeometry*> rings;
2698 
2699  int i;
2700  for ( i = 0; i < polygonList.size(); i++ )
2701  {
2702  for ( int j = 0; j < rings.size(); j++ )
2703  GEOSGeom_destroy( rings[j] );
2704  rings.clear();
2705 
2706  GEOSGeometry *shellRing = 0;
2707  GEOSGeometry *shell = 0;
2708  try
2709  {
2710  shellRing = GEOSGeom_clone( GEOSGetExteriorRing( polygonList[i] ) );
2711  shell = createGeosPolygon( shellRing );
2712 
2713  if ( !GEOSWithin( newRingPolygon, shell ) )
2714  {
2715  GEOSGeom_destroy( shell );
2716  continue;
2717  }
2718  }
2719  catch ( GEOSException &e )
2720  {
2721  Q_UNUSED( e );
2722  if ( shell )
2723  GEOSGeom_destroy( shell );
2724  else if ( shellRing )
2725  GEOSGeom_destroy( shellRing );
2726 
2727  GEOSGeom_destroy( newRingPolygon );
2728 
2729  return 4;
2730  }
2731 
2732  // add outer ring
2733  rings << GEOSGeom_clone( shellRing );
2734 
2735  GEOSGeom_destroy( shell );
2736 
2737  // check inner rings
2738  int n = GEOSGetNumInteriorRings( polygonList[i] );
2739 
2740  int j;
2741  for ( j = 0; j < n; j++ )
2742  {
2743  GEOSGeometry *holeRing = 0;
2744  GEOSGeometry *hole = 0;
2745  try
2746  {
2747  holeRing = GEOSGeom_clone( GEOSGetInteriorRingN( polygonList[i], j ) );
2748  hole = createGeosPolygon( holeRing );
2749 
2750  if ( !GEOSDisjoint( hole, newRingPolygon ) )
2751  {
2752  GEOSGeom_destroy( hole );
2753  break;
2754  }
2755  }
2756  catch ( GEOSException &e )
2757  {
2758  Q_UNUSED( e );
2759  if ( hole )
2760  GEOSGeom_destroy( hole );
2761  else if ( holeRing )
2762  GEOSGeom_destroy( holeRing );
2763 
2764  break;
2765  }
2766 
2767  rings << GEOSGeom_clone( holeRing );
2768  GEOSGeom_destroy( hole );
2769  }
2770 
2771  if ( j == n )
2772  // this is it...
2773  break;
2774  }
2775 
2776  if ( i == polygonList.size() )
2777  {
2778  // clear rings
2779  for ( int j = 0; j < rings.size(); j++ )
2780  GEOSGeom_destroy( rings[j] );
2781  rings.clear();
2782 
2783  GEOSGeom_destroy( newRingPolygon );
2784 
2785  // no containing polygon found
2786  return 5;
2787  }
2788 
2789  rings << GEOSGeom_clone( newRing );
2790  GEOSGeom_destroy( newRingPolygon );
2791 
2792  GEOSGeometry *newPolygon = createGeosPolygon( rings );
2793 
2794  if ( wkbType() == QGis::WKBPolygon )
2795  {
2796  GEOSGeom_destroy( mGeos );
2797  mGeos = newPolygon;
2798  }
2799  else if ( wkbType() == QGis::WKBMultiPolygon )
2800  {
2801  QVector<GEOSGeometry*> newPolygons;
2802 
2803  for ( int j = 0; j < polygonList.size(); j++ )
2804  {
2805  newPolygons << ( i == j ? newPolygon : GEOSGeom_clone( polygonList[j] ) );
2806  }
2807 
2808  GEOSGeom_destroy( mGeos );
2809  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, newPolygons );
2810  }
2811 
2812  mDirtyWkb = true;
2813  mDirtyGeos = false;
2814  return 0;
2815 }
2816 
2817 int QgsGeometry::addIsland( const QList<QgsPoint>& ring )
2818 {
2819  //Ring needs to have at least three points and must be closed
2820  if ( ring.size() < 4 )
2821  {
2822  return 2;
2823  }
2824 
2825  //ring must be closed
2826  if ( ring.first() != ring.last() )
2827  {
2828  return 2;
2829  }
2830 
2832  {
2833  if ( !convertToMultiType() )
2834  {
2835  return 1;
2836  }
2837  }
2838 
2839  //bail out if wkbtype is not multipolygon
2841  {
2842  return 1;
2843  }
2844 
2845  //create geos geometry from wkb if not already there
2846  if ( !mGeos || mDirtyGeos )
2847  {
2848  if ( !exportWkbToGeos() )
2849  return 4;
2850  }
2851 
2852  if ( GEOSGeomTypeId( mGeos ) != GEOS_MULTIPOLYGON )
2853  return 1;
2854 
2855  //create new polygon from ring
2856 
2857  //coordinate sequence first
2858  GEOSGeometry *newRing = 0;
2859  GEOSGeometry *newPolygon = 0;
2860 
2861  try
2862  {
2863  newRing = createGeosLinearRing( ring.toVector() );
2864  if ( !GEOSisValid( newRing ) )
2865  throw GEOSException( "ring invalid" );
2866  newPolygon = createGeosPolygon( newRing );
2867  if ( !GEOSisValid( newPolygon ) )
2868  throw GEOSException( "polygon invalid" );
2869  }
2870  catch ( GEOSException &e )
2871  {
2872  Q_UNUSED( e );
2873  if ( newPolygon )
2874  GEOSGeom_destroy( newPolygon );
2875  else if ( newRing )
2876  GEOSGeom_destroy( newRing );
2877  return 2;
2878  }
2879 
2880  QVector<GEOSGeometry*> polygons;
2881 
2882  //create new multipolygon
2883  int n = GEOSGetNumGeometries( mGeos );
2884  int i;
2885  for ( i = 0; i < n; ++i )
2886  {
2887  const GEOSGeometry *polygonN = GEOSGetGeometryN( mGeos, i );
2888 
2889  if ( !GEOSDisjoint( polygonN, newPolygon ) )
2890  //bail out if new polygon is not disjoint with existing ones
2891  break;
2892 
2893  polygons << GEOSGeom_clone( polygonN );
2894  }
2895 
2896  if ( i < n )
2897  {
2898  // bailed out
2899  for ( int i = 0; i < polygons.size(); i++ )
2900  GEOSGeom_destroy( polygons[i] );
2901  return 3;
2902  }
2903 
2904  polygons << newPolygon;
2905 
2906  GEOSGeom_destroy( mGeos );
2907  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
2908  mDirtyWkb = true;
2909  mDirtyGeos = false;
2910  return 0;
2911 }
2912 
2913 int QgsGeometry::translate( double dx, double dy )
2914 {
2915  if ( mDirtyWkb )
2916  {
2917  exportGeosToWkb();
2918  }
2919 
2920  if ( !mGeometry )
2921  {
2922  QgsDebugMsg( "WKB geometry not available!" );
2923  return 1;
2924  }
2925 
2927  memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
2928  bool hasZValue = false;
2929  int wkbPosition = 5;
2930 
2931  switch ( wkbType )
2932  {
2933  case QGis::WKBPoint25D:
2934  case QGis::WKBPoint:
2935  {
2936  translateVertex( wkbPosition, dx, dy, hasZValue );
2937  }
2938  break;
2939 
2941  hasZValue = true;
2942  case QGis::WKBLineString:
2943  {
2944  int* npoints = ( int* )( &mGeometry[wkbPosition] );
2945  wkbPosition += sizeof( int );
2946  for ( int index = 0; index < *npoints; ++index )
2947  {
2948  translateVertex( wkbPosition, dx, dy, hasZValue );
2949  }
2950  break;
2951  }
2952 
2953  case QGis::WKBPolygon25D:
2954  hasZValue = true;
2955  case QGis::WKBPolygon:
2956  {
2957  int* nrings = ( int* )( &( mGeometry[wkbPosition] ) );
2958  wkbPosition += sizeof( int );
2959  int* npoints;
2960 
2961  for ( int index = 0; index < *nrings; ++index )
2962  {
2963  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
2964  wkbPosition += sizeof( int );
2965  for ( int index2 = 0; index2 < *npoints; ++index2 )
2966  {
2967  translateVertex( wkbPosition, dx, dy, hasZValue );
2968  }
2969  }
2970  break;
2971  }
2972 
2974  hasZValue = true;
2975  case QGis::WKBMultiPoint:
2976  {
2977  int* npoints = ( int* )( &( mGeometry[wkbPosition] ) );
2978  wkbPosition += sizeof( int );
2979  for ( int index = 0; index < *npoints; ++index )
2980  {
2981  wkbPosition += ( sizeof( int ) + 1 );
2982  translateVertex( wkbPosition, dx, dy, hasZValue );
2983  }
2984  break;
2985  }
2986 
2988  hasZValue = true;
2990  {
2991  int* nlines = ( int* )( &( mGeometry[wkbPosition] ) );
2992  int* npoints = 0;
2993  wkbPosition += sizeof( int );
2994  for ( int index = 0; index < *nlines; ++index )
2995  {
2996  wkbPosition += ( sizeof( int ) + 1 );
2997  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
2998  wkbPosition += sizeof( int );
2999  for ( int index2 = 0; index2 < *npoints; ++index2 )
3000  {
3001  translateVertex( wkbPosition, dx, dy, hasZValue );
3002  }
3003  }
3004  break;
3005  }
3006 
3008  hasZValue = true;
3009  case QGis::WKBMultiPolygon:
3010  {
3011  int* npolys = ( int* )( &( mGeometry[wkbPosition] ) );
3012  int* nrings;
3013  int* npoints;
3014  wkbPosition += sizeof( int );
3015  for ( int index = 0; index < *npolys; ++index )
3016  {
3017  wkbPosition += ( 1 + sizeof( int ) ); //skip endian and polygon type
3018  nrings = ( int* )( &( mGeometry[wkbPosition] ) );
3019  wkbPosition += sizeof( int );
3020  for ( int index2 = 0; index2 < *nrings; ++index2 )
3021  {
3022  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
3023  wkbPosition += sizeof( int );
3024  for ( int index3 = 0; index3 < *npoints; ++index3 )
3025  {
3026  translateVertex( wkbPosition, dx, dy, hasZValue );
3027  }
3028  }
3029  }
3030  }
3031 
3032  default:
3033  break;
3034  }
3035  mDirtyGeos = true;
3036  return 0;
3037 }
3038 
3040 {
3041  if ( mDirtyWkb )
3042  {
3043  exportGeosToWkb();
3044  }
3045 
3046  if ( !mGeometry )
3047  {
3048  QgsDebugMsg( "WKB geometry not available!" );
3049  return 1;
3050  }
3051 
3053  memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
3054  bool hasZValue = false;
3055  int wkbPosition = 5;
3056 
3057  switch ( wkbType )
3058  {
3059  case QGis::WKBPoint25D:
3060  case QGis::WKBPoint:
3061  {
3062  transformVertex( wkbPosition, ct, hasZValue );
3063  }
3064  break;
3065 
3067  hasZValue = true;
3068  case QGis::WKBLineString:
3069  {
3070  int* npoints = ( int* )( &mGeometry[wkbPosition] );
3071  wkbPosition += sizeof( int );
3072  for ( int index = 0; index < *npoints; ++index )
3073  {
3074  transformVertex( wkbPosition, ct, hasZValue );
3075  }
3076  break;
3077  }
3078 
3079  case QGis::WKBPolygon25D:
3080  hasZValue = true;
3081  case QGis::WKBPolygon:
3082  {
3083  int* nrings = ( int* )( &( mGeometry[wkbPosition] ) );
3084  wkbPosition += sizeof( int );
3085  int* npoints;
3086 
3087  for ( int index = 0; index < *nrings; ++index )
3088  {
3089  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
3090  wkbPosition += sizeof( int );
3091  for ( int index2 = 0; index2 < *npoints; ++index2 )
3092  {
3093  transformVertex( wkbPosition, ct, hasZValue );
3094  }
3095  }
3096  break;
3097  }
3098 
3100  hasZValue = true;
3101  case QGis::WKBMultiPoint:
3102  {
3103  int* npoints = ( int* )( &( mGeometry[wkbPosition] ) );
3104  wkbPosition += sizeof( int );
3105  for ( int index = 0; index < *npoints; ++index )
3106  {
3107  wkbPosition += ( sizeof( int ) + 1 );
3108  transformVertex( wkbPosition, ct, hasZValue );
3109  }
3110  break;
3111  }
3112 
3114  hasZValue = true;
3116  {
3117  int* nlines = ( int* )( &( mGeometry[wkbPosition] ) );
3118  int* npoints = 0;
3119  wkbPosition += sizeof( int );
3120  for ( int index = 0; index < *nlines; ++index )
3121  {
3122  wkbPosition += ( sizeof( int ) + 1 );
3123  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
3124  wkbPosition += sizeof( int );
3125  for ( int index2 = 0; index2 < *npoints; ++index2 )
3126  {
3127  transformVertex( wkbPosition, ct, hasZValue );
3128  }
3129  }
3130  break;
3131  }
3132 
3134  hasZValue = true;
3135  case QGis::WKBMultiPolygon:
3136  {
3137  int* npolys = ( int* )( &( mGeometry[wkbPosition] ) );
3138  int* nrings;
3139  int* npoints;
3140  wkbPosition += sizeof( int );
3141  for ( int index = 0; index < *npolys; ++index )
3142  {
3143  wkbPosition += ( 1 + sizeof( int ) ); //skip endian and polygon type
3144  nrings = ( int* )( &( mGeometry[wkbPosition] ) );
3145  wkbPosition += sizeof( int );
3146  for ( int index2 = 0; index2 < *nrings; ++index2 )
3147  {
3148  npoints = ( int* )( &( mGeometry[wkbPosition] ) );
3149  wkbPosition += sizeof( int );
3150  for ( int index3 = 0; index3 < *npoints; ++index3 )
3151  {
3152  transformVertex( wkbPosition, ct, hasZValue );
3153  }
3154  }
3155  }
3156  }
3157 
3158  default:
3159  break;
3160  }
3161  mDirtyGeos = true;
3162  return 0;
3163 }
3164 
3165 int QgsGeometry::splitGeometry( const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint> &topologyTestPoints )
3166 {
3167  int returnCode = 0;
3168 
3169  //return if this type is point/multipoint
3170  if ( type() == QGis::Point )
3171  {
3172  return 1; //cannot split points
3173  }
3174 
3175  //make sure, mGeos and mWkb are there and up-to-date
3176  if ( mDirtyWkb )
3177  {
3178  exportGeosToWkb();
3179  }
3180 
3181  if ( !mGeos || mDirtyGeos )
3182  {
3183  if ( !exportWkbToGeos() )
3184  return 1;
3185  }
3186 
3187  if ( !GEOSisValid( mGeos ) )
3188  {
3189  return 7;
3190  }
3191 
3192  //make sure splitLine is valid
3193  if ( splitLine.size() < 2 )
3194  {
3195  return 1;
3196  }
3197 
3198  newGeometries.clear();
3199 
3200  try
3201  {
3202  GEOSGeometry *splitLineGeos = createGeosLineString( splitLine.toVector() );
3203  if ( !GEOSisValid( splitLineGeos ) || !GEOSisSimple( splitLineGeos ) )
3204  {
3205  GEOSGeom_destroy( splitLineGeos );
3206  return 1;
3207  }
3208 
3209  if ( topological )
3210  {
3211  //find out candidate points for topological corrections
3212  if ( topologicalTestPointsSplit( splitLineGeos, topologyTestPoints ) != 0 )
3213  {
3214  return 1;
3215  }
3216  }
3217 
3218  //call split function depending on geometry type
3219  if ( type() == QGis::Line )
3220  {
3221  returnCode = splitLinearGeometry( splitLineGeos, newGeometries );
3222  GEOSGeom_destroy( splitLineGeos );
3223  }
3224  else if ( type() == QGis::Polygon )
3225  {
3226  returnCode = splitPolygonGeometry( splitLineGeos, newGeometries );
3227  GEOSGeom_destroy( splitLineGeos );
3228  }
3229  else
3230  {
3231  return 1;
3232  }
3233  }
3234  CATCH_GEOS( 2 )
3235 
3236  return returnCode;
3237 }
3238 
3240 int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine )
3241 {
3242  if ( reshapeWithLine.size() < 2 )
3243  {
3244  return 1;
3245  }
3246 
3247  if ( type() == QGis::Point )
3248  {
3249  return 1; //cannot reshape points
3250  }
3251 
3252  GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() );
3253 
3254  //make sure this geos geometry is up-to-date
3255  if ( !mGeos || mDirtyGeos )
3256  {
3257  exportWkbToGeos();
3258  }
3259 
3260  //single or multi?
3261  int numGeoms = GEOSGetNumGeometries( mGeos );
3262  if ( numGeoms == -1 )
3263  {
3264  return 1;
3265  }
3266 
3267  bool isMultiGeom = false;
3268  int geosTypeId = GEOSGeomTypeId( mGeos );
3269  if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
3270  {
3271  isMultiGeom = true;
3272  }
3273 
3274  bool isLine = ( type() == QGis::Line );
3275 
3276  //polygon or multipolygon?
3277  if ( !isMultiGeom )
3278  {
3279  GEOSGeometry* reshapedGeometry;
3280  if ( isLine )
3281  {
3282  reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos );
3283  }
3284  else
3285  {
3286  reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos );
3287  }
3288 
3289  GEOSGeom_destroy( reshapeLineGeos );
3290  if ( reshapedGeometry )
3291  {
3292  GEOSGeom_destroy( mGeos );
3293  mGeos = reshapedGeometry;
3294  mDirtyWkb = true;
3295  return 0;
3296  }
3297  else
3298  {
3299  return 1;
3300  }
3301  }
3302  else
3303  {
3304  //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
3305  bool reshapeTookPlace = false;
3306 
3307  GEOSGeometry* currentReshapeGeometry = 0;
3308  GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms];
3309 
3310  for ( int i = 0; i < numGeoms; ++i )
3311  {
3312  if ( isLine )
3313  {
3314  currentReshapeGeometry = reshapeLine( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos );
3315  }
3316  else
3317  {
3318  currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos );
3319  }
3320 
3321  if ( currentReshapeGeometry )
3322  {
3323  newGeoms[i] = currentReshapeGeometry;
3324  reshapeTookPlace = true;
3325  }
3326  else
3327  {
3328  newGeoms[i] = GEOSGeom_clone( GEOSGetGeometryN( mGeos, i ) );
3329  }
3330  }
3331  GEOSGeom_destroy( reshapeLineGeos );
3332 
3333  GEOSGeometry* newMultiGeom = 0;
3334  if ( isLine )
3335  {
3336  newMultiGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, newGeoms, numGeoms );
3337  }
3338  else //multipolygon
3339  {
3340  newMultiGeom = GEOSGeom_createCollection( GEOS_MULTIPOLYGON, newGeoms, numGeoms );
3341  }
3342 
3343  delete[] newGeoms;
3344  if ( ! newMultiGeom )
3345  {
3346  return 3;
3347  }
3348 
3349  if ( reshapeTookPlace )
3350  {
3351  GEOSGeom_destroy( mGeos );
3352  mGeos = newMultiGeom;
3353  mDirtyWkb = true;
3354  return 0;
3355  }
3356  else
3357  {
3358  GEOSGeom_destroy( newMultiGeom );
3359  return 1;
3360  }
3361  }
3362 }
3363 
3365 {
3366  //make sure geos geometry is up to date
3367  if ( !mGeos || mDirtyGeos )
3368  {
3369  exportWkbToGeos();
3370  }
3371 
3372  if ( !mGeos )
3373  {
3374  return 1;
3375  }
3376 
3377  if ( !GEOSisValid( mGeos ) )
3378  {
3379  return 2;
3380  }
3381 
3382  if ( !GEOSisSimple( mGeos ) )
3383  {
3384  return 3;
3385  }
3386 
3387  //convert other geometry to geos
3388  if ( !other->mGeos || other->mDirtyGeos )
3389  {
3390  other->exportWkbToGeos();
3391  }
3392 
3393  if ( !other->mGeos )
3394  {
3395  return 4;
3396  }
3397 
3398  //make geometry::difference
3399  try
3400  {
3401  if ( GEOSIntersects( mGeos, other->mGeos ) )
3402  {
3403  //check if multitype before and after
3404  bool multiType = isMultipart();
3405 
3406  mGeos = GEOSDifference( mGeos, other->mGeos );
3407  mDirtyWkb = true;
3408 
3409  if ( multiType && !isMultipart() )
3410  {
3412  exportWkbToGeos();
3413  }
3414  }
3415  else
3416  {
3417  return 0; //nothing to do
3418  }
3419  }
3420  CATCH_GEOS( 5 )
3421 
3422  if ( !mGeos )
3423  {
3424  mDirtyGeos = true;
3425  return 6;
3426  }
3427 
3428  return 0;
3429 }
3430 
3432 {
3433  double xmin = std::numeric_limits<double>::max();
3434  double ymin = std::numeric_limits<double>::max();
3435  double xmax = -std::numeric_limits<double>::max();
3436  double ymax = -std::numeric_limits<double>::max();
3437 
3438  double *x;
3439  double *y;
3440  int *nPoints;
3441  int *numRings;
3442  int *numPolygons;
3443  int numLineStrings;
3444  int idx, jdx, kdx;
3445  unsigned char *ptr;
3446  QgsPoint pt;
3448  bool hasZValue = false;
3449 
3450  // TODO: implement with GEOS
3451  if ( mDirtyWkb )
3452  {
3453  exportGeosToWkb();
3454  }
3455 
3456  if ( !mGeometry )
3457  {
3458  QgsDebugMsg( "WKB geometry not available!" );
3459  return QgsRectangle( 0, 0, 0, 0 );
3460  }
3461  // consider endian when fetching feature type
3462  //wkbType = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4]; //MH: Does not work for 25D geometries
3463  memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
3464  switch ( wkbType )
3465  {
3466  case QGis::WKBPoint25D:
3467  case QGis::WKBPoint:
3468  x = ( double * )( mGeometry + 5 );
3469  y = ( double * )( mGeometry + 5 + sizeof( double ) );
3470  if ( *x < xmin )
3471  {
3472  xmin = *x;
3473  }
3474  if ( *x > xmax )
3475  {
3476  xmax = *x;
3477  }
3478  if ( *y < ymin )
3479  {
3480  ymin = *y;
3481  }
3482  if ( *y > ymax )
3483  {
3484  ymax = *y;
3485  }
3486  break;
3488  hasZValue = true;
3489  case QGis::WKBMultiPoint:
3490  {
3491  ptr = mGeometry + 1 + sizeof( int );
3492  nPoints = ( int * ) ptr;
3493  ptr += sizeof( int );
3494  for ( idx = 0; idx < *nPoints; idx++ )
3495  {
3496  ptr += ( 1 + sizeof( int ) );
3497  x = ( double * ) ptr;
3498  ptr += sizeof( double );
3499  y = ( double * ) ptr;
3500  ptr += sizeof( double );
3501  if ( hasZValue )
3502  {
3503  ptr += sizeof( double );
3504  }
3505  if ( *x < xmin )
3506  {
3507  xmin = *x;
3508  }
3509  if ( *x > xmax )
3510  {
3511  xmax = *x;
3512  }
3513  if ( *y < ymin )
3514  {
3515  ymin = *y;
3516  }
3517  if ( *y > ymax )
3518  {
3519  ymax = *y;
3520  }
3521  }
3522  break;
3523  }
3525  hasZValue = true;
3526  case QGis::WKBLineString:
3527  {
3528  // get number of points in the line
3529  ptr = mGeometry + 5;
3530  nPoints = ( int * ) ptr;
3531  ptr = mGeometry + 1 + 2 * sizeof( int );
3532  for ( idx = 0; idx < *nPoints; idx++ )
3533  {
3534  x = ( double * ) ptr;
3535  ptr += sizeof( double );
3536  y = ( double * ) ptr;
3537  ptr += sizeof( double );
3538  if ( hasZValue )
3539  {
3540  ptr += sizeof( double );
3541  }
3542  if ( *x < xmin )
3543  {
3544  xmin = *x;
3545  }
3546  if ( *x > xmax )
3547  {
3548  xmax = *x;
3549  }
3550  if ( *y < ymin )
3551  {
3552  ymin = *y;
3553  }
3554  if ( *y > ymax )
3555  {
3556  ymax = *y;
3557  }
3558  }
3559  break;
3560  }
3562  hasZValue = true;
3564  {
3565  numLineStrings = ( int )( mGeometry[5] );
3566  ptr = mGeometry + 9;
3567  for ( jdx = 0; jdx < numLineStrings; jdx++ )
3568  {
3569  // each of these is a wbklinestring so must handle as such
3570  ptr += 5; // skip type since we know its 2
3571  nPoints = ( int * ) ptr;
3572  ptr += sizeof( int );
3573  for ( idx = 0; idx < *nPoints; idx++ )
3574  {
3575  x = ( double * ) ptr;
3576  ptr += sizeof( double );
3577  y = ( double * ) ptr;
3578  ptr += sizeof( double );
3579  if ( hasZValue )
3580  {
3581  ptr += sizeof( double );
3582  }
3583  if ( *x < xmin )
3584  {
3585  xmin = *x;
3586  }
3587  if ( *x > xmax )
3588  {
3589  xmax = *x;
3590  }
3591  if ( *y < ymin )
3592  {
3593  ymin = *y;
3594  }
3595  if ( *y > ymax )
3596  {
3597  ymax = *y;
3598  }
3599  }
3600  }
3601  break;
3602  }
3603  case QGis::WKBPolygon25D:
3604  hasZValue = true;
3605  case QGis::WKBPolygon:
3606  {
3607  // get number of rings in the polygon
3608  numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
3609  ptr = mGeometry + 1 + 2 * sizeof( int );
3610  for ( idx = 0; idx < *numRings; idx++ )
3611  {
3612  // get number of points in the ring
3613  nPoints = ( int * ) ptr;
3614  ptr += 4;
3615  for ( jdx = 0; jdx < *nPoints; jdx++ )
3616  {
3617  // add points to a point array for drawing the polygon
3618  x = ( double * ) ptr;
3619  ptr += sizeof( double );
3620  y = ( double * ) ptr;
3621  ptr += sizeof( double );
3622  if ( hasZValue )
3623  {
3624  ptr += sizeof( double );
3625  }
3626  if ( *x < xmin )
3627  {
3628  xmin = *x;
3629  }
3630  if ( *x > xmax )
3631  {
3632  xmax = *x;
3633  }
3634  if ( *y < ymin )
3635  {
3636  ymin = *y;
3637  }
3638  if ( *y > ymax )
3639  {
3640  ymax = *y;
3641  }
3642  }
3643  }
3644  break;
3645  }
3647  hasZValue = true;
3648  case QGis::WKBMultiPolygon:
3649  {
3650  // get the number of polygons
3651  ptr = mGeometry + 5;
3652  numPolygons = ( int * ) ptr;
3653  ptr += 4;
3654 
3655  for ( kdx = 0; kdx < *numPolygons; kdx++ )
3656  {
3657  //skip the endian and mGeometry type info and
3658  // get number of rings in the polygon
3659  ptr += 5;
3660  numRings = ( int * ) ptr;
3661  ptr += 4;
3662  for ( idx = 0; idx < *numRings; idx++ )
3663  {
3664  // get number of points in the ring
3665  nPoints = ( int * ) ptr;
3666  ptr += 4;
3667  for ( jdx = 0; jdx < *nPoints; jdx++ )
3668  {
3669  // add points to a point array for drawing the polygon
3670  x = ( double * ) ptr;
3671  ptr += sizeof( double );
3672  y = ( double * ) ptr;
3673  ptr += sizeof( double );
3674  if ( hasZValue )
3675  {
3676  ptr += sizeof( double );
3677  }
3678  if ( *x < xmin )
3679  {
3680  xmin = *x;
3681  }
3682  if ( *x > xmax )
3683  {
3684  xmax = *x;
3685  }
3686  if ( *y < ymin )
3687  {
3688  ymin = *y;
3689  }
3690  if ( *y > ymax )
3691  {
3692  ymax = *y;
3693  }
3694  }
3695  }
3696  }
3697  break;
3698  }
3699 
3700  default:
3701  QgsDebugMsg( "Unknown WkbType ENCOUNTERED" );
3702  return QgsRectangle( 0, 0, 0, 0 );
3703  break;
3704 
3705  }
3706  return QgsRectangle( xmin, ymin, xmax, ymax );
3707 }
3708 
3710 {
3711  QgsGeometry* g = fromRect( r );
3712  bool res = intersects( g );
3713  delete g;
3714  return res;
3715 }
3716 
3718 {
3719  try // geos might throw exception on error
3720  {
3721  // ensure that both geometries have geos geometry
3722  exportWkbToGeos();
3723  geometry->exportWkbToGeos();
3724 
3725  if ( !mGeos || !geometry->mGeos )
3726  {
3727  QgsDebugMsg( "GEOS geometry not available!" );
3728  return false;
3729  }
3730 
3731  return GEOSIntersects( mGeos, geometry->mGeos );
3732  }
3733  CATCH_GEOS( false )
3734 }
3735 
3736 
3738 {
3739  exportWkbToGeos();
3740 
3741  if ( !mGeos )
3742  {
3743  QgsDebugMsg( "GEOS geometry not available!" );
3744  return false;
3745  }
3746 
3747  GEOSGeometry *geosPoint = 0;
3748 
3749  bool returnval = false;
3750 
3751  try
3752  {
3753  geosPoint = createGeosPoint( *p );
3754  returnval = GEOSContains( mGeos, geosPoint );
3755  }
3756  catch ( GEOSException &e )
3757  {
3758  Q_UNUSED( e );
3759  returnval = false;
3760  }
3761 
3762  if ( geosPoint )
3763  GEOSGeom_destroy( geosPoint );
3764 
3765  return returnval;
3766 }
3767 
3769  char( *op )( const GEOSGeometry*, const GEOSGeometry * ),
3770  QgsGeometry *a,
3771  QgsGeometry *b )
3772 {
3773  try // geos might throw exception on error
3774  {
3775  // ensure that both geometries have geos geometry
3776  a->exportWkbToGeos();
3777  b->exportWkbToGeos();
3778 
3779  if ( !a->mGeos || !b->mGeos )
3780  {
3781  QgsDebugMsg( "GEOS geometry not available!" );
3782  return false;
3783  }
3784  return op( a->mGeos, b->mGeos );
3785  }
3786  CATCH_GEOS( false )
3787 }
3788 
3790 {
3791  return geosRelOp( GEOSContains, this, geometry );
3792 }
3793 
3795 {
3796  return geosRelOp( GEOSDisjoint, this, geometry );
3797 }
3798 
3800 {
3801  return geosRelOp( GEOSEquals, this, geometry );
3802 }
3803 
3805 {
3806  return geosRelOp( GEOSTouches, this, geometry );
3807 }
3808 
3810 {
3811  return geosRelOp( GEOSOverlaps, this, geometry );
3812 }
3813 
3815 {
3816  return geosRelOp( GEOSWithin, this, geometry );
3817 }
3818 
3820 {
3821  return geosRelOp( GEOSCrosses, this, geometry );
3822 }
3823 
3825 {
3826  QgsDebugMsg( "entered." );
3827 
3828  // TODO: implement with GEOS
3829  if ( mDirtyWkb )
3830  {
3831  exportGeosToWkb();
3832  }
3833 
3834  if ( !mGeometry )
3835  {
3836  QgsDebugMsg( "WKB geometry not available!" );
3837  return QString::null;
3838  }
3839 
3841  bool hasZValue = false;
3842  double *x, *y;
3843 
3844  QString mWkt; // TODO: rename
3845 
3846  // Will this really work when mGeometry[0] == 0 ???? I (gavin) think not.
3847  //wkbType = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4];
3848  memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
3849 
3850  switch ( wkbType )
3851  {
3852  case QGis::WKBPoint25D:
3853  case QGis::WKBPoint:
3854  {
3855  mWkt += "POINT(";
3856  x = ( double * )( mGeometry + 5 );
3857  mWkt += QString::number( *x, 'f', 6 );
3858  mWkt += " ";
3859  y = ( double * )( mGeometry + 5 + sizeof( double ) );
3860  mWkt += QString::number( *y, 'f', 6 );
3861  mWkt += ")";
3862  return mWkt;
3863  }
3864 
3866  hasZValue = true;
3867  case QGis::WKBLineString:
3868  {
3869  QgsDebugMsg( "LINESTRING found" );
3870  unsigned char *ptr;
3871  int *nPoints;
3872  int idx;
3873 
3874  mWkt += "LINESTRING(";
3875  // get number of points in the line
3876  ptr = mGeometry + 5;
3877  nPoints = ( int * ) ptr;
3878  ptr = mGeometry + 1 + 2 * sizeof( int );
3879  for ( idx = 0; idx < *nPoints; ++idx )
3880  {
3881  if ( idx != 0 )
3882  {
3883  mWkt += ", ";
3884  }
3885  x = ( double * ) ptr;
3886  mWkt += QString::number( *x, 'f', 6 );
3887  mWkt += " ";
3888  ptr += sizeof( double );
3889  y = ( double * ) ptr;
3890  mWkt += QString::number( *y, 'f', 6 );
3891  ptr += sizeof( double );
3892  if ( hasZValue )
3893  {
3894  ptr += sizeof( double );
3895  }
3896  }
3897  mWkt += ")";
3898  return mWkt;
3899  }
3900 
3901  case QGis::WKBPolygon25D:
3902  hasZValue = true;
3903  case QGis::WKBPolygon:
3904  {
3905  QgsDebugMsg( "POLYGON found" );
3906  unsigned char *ptr;
3907  int idx, jdx;
3908  int *numRings, *nPoints;
3909 
3910  mWkt += "POLYGON(";
3911  // get number of rings in the polygon
3912  numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
3913  if ( !( *numRings ) ) // sanity check for zero rings in polygon
3914  {
3915  return QString();
3916  }
3917  int *ringStart; // index of first point for each ring
3918  int *ringNumPoints; // number of points in each ring
3919  ringStart = new int[*numRings];
3920  ringNumPoints = new int[*numRings];
3921  ptr = mGeometry + 1 + 2 * sizeof( int ); // set pointer to the first ring
3922  for ( idx = 0; idx < *numRings; idx++ )
3923  {
3924  if ( idx != 0 )
3925  {
3926  mWkt += ",";
3927  }
3928  mWkt += "(";
3929  // get number of points in the ring
3930  nPoints = ( int * ) ptr;
3931  ringNumPoints[idx] = *nPoints;
3932  ptr += 4;
3933 
3934  for ( jdx = 0; jdx < *nPoints; jdx++ )
3935  {
3936  if ( jdx != 0 )
3937  {
3938  mWkt += ",";
3939  }
3940  x = ( double * ) ptr;
3941  mWkt += QString::number( *x, 'f', 6 );
3942  mWkt += " ";
3943  ptr += sizeof( double );
3944  y = ( double * ) ptr;
3945  mWkt += QString::number( *y, 'f', 6 );
3946  ptr += sizeof( double );
3947  if ( hasZValue )
3948  {
3949  ptr += sizeof( double );
3950  }
3951  }
3952  mWkt += ")";
3953  }
3954  mWkt += ")";
3955  delete [] ringStart;
3956  delete [] ringNumPoints;
3957  return mWkt;
3958  }
3959 
3961  hasZValue = true;
3962  case QGis::WKBMultiPoint:
3963  {
3964  unsigned char *ptr;
3965  int idx;
3966  int *nPoints;
3967 
3968  mWkt += "MULTIPOINT(";
3969  nPoints = ( int* )( mGeometry + 5 );
3970  ptr = mGeometry + 5 + sizeof( int );
3971  for ( idx = 0; idx < *nPoints; ++idx )
3972  {
3973  ptr += ( 1 + sizeof( int ) );
3974  if ( idx != 0 )
3975  {
3976  mWkt += ", ";
3977  }
3978  x = ( double * )( ptr );
3979  mWkt += QString::number( *x, 'f', 6 );
3980  mWkt += " ";
3981  ptr += sizeof( double );
3982  y = ( double * )( ptr );
3983  mWkt += QString::number( *y, 'f', 6 );
3984  ptr += sizeof( double );
3985  if ( hasZValue )
3986  {
3987  ptr += sizeof( double );
3988  }
3989  }
3990  mWkt += ")";
3991  return mWkt;
3992  }
3993 
3995  hasZValue = true;
3997  {
3998  QgsDebugMsg( "MULTILINESTRING found" );
3999  unsigned char *ptr;
4000  int idx, jdx, numLineStrings;
4001  int *nPoints;
4002 
4003  mWkt += "MULTILINESTRING(";
4004  numLineStrings = ( int )( mGeometry[5] );
4005  ptr = mGeometry + 9;
4006  for ( jdx = 0; jdx < numLineStrings; jdx++ )
4007  {
4008  if ( jdx != 0 )
4009  {
4010  mWkt += ", ";
4011  }
4012  mWkt += "(";
4013  ptr += 5; // skip type since we know its 2
4014  nPoints = ( int * ) ptr;
4015  ptr += sizeof( int );
4016  for ( idx = 0; idx < *nPoints; idx++ )
4017  {
4018  if ( idx != 0 )
4019  {
4020  mWkt += ", ";
4021  }
4022  x = ( double * ) ptr;
4023  mWkt += QString::number( *x, 'f', 6 );
4024  ptr += sizeof( double );
4025  mWkt += " ";
4026  y = ( double * ) ptr;
4027  mWkt += QString::number( *y, 'f', 6 );
4028  ptr += sizeof( double );
4029  if ( hasZValue )
4030  {
4031  ptr += sizeof( double );
4032  }
4033  }
4034  mWkt += ")";
4035  }
4036  mWkt += ")";
4037  return mWkt;
4038  }
4039 
4041  hasZValue = true;
4042  case QGis::WKBMultiPolygon:
4043  {
4044  QgsDebugMsg( "MULTIPOLYGON found" );
4045  unsigned char *ptr;
4046  int idx, jdx, kdx;
4047  int *numPolygons, *numRings, *nPoints;
4048 
4049  mWkt += "MULTIPOLYGON(";
4050  ptr = mGeometry + 5;
4051  numPolygons = ( int * ) ptr;
4052  ptr = mGeometry + 9;
4053  for ( kdx = 0; kdx < *numPolygons; kdx++ )
4054  {
4055  if ( kdx != 0 )
4056  {
4057  mWkt += ",";
4058  }
4059  mWkt += "(";
4060  ptr += 5;
4061  numRings = ( int * ) ptr;
4062  ptr += 4;
4063  for ( idx = 0; idx < *numRings; idx++ )
4064  {
4065  if ( idx != 0 )
4066  {
4067  mWkt += ",";
4068  }
4069  mWkt += "(";
4070  nPoints = ( int * ) ptr;
4071  ptr += 4;
4072  for ( jdx = 0; jdx < *nPoints; jdx++ )
4073  {
4074  if ( jdx != 0 )
4075  {
4076  mWkt += ",";
4077  }
4078  x = ( double * ) ptr;
4079  mWkt += QString::number( *x, 'f', 6 );
4080  ptr += sizeof( double );
4081  mWkt += " ";
4082  y = ( double * ) ptr;
4083  mWkt += QString::number( *y, 'f', 6 );
4084  ptr += sizeof( double );
4085  if ( hasZValue )
4086  {
4087  ptr += sizeof( double );
4088  }
4089  }
4090  mWkt += ")";
4091  }
4092  mWkt += ")";
4093  }
4094  mWkt += ")";
4095  return mWkt;
4096  }
4097 
4098  default:
4099  QgsDebugMsg( "error: mGeometry type not recognized" );
4100  return QString::null;
4101  }
4102 }
4103 
4105 {
4106  QgsDebugMsgLevel( "entered.", 3 );
4107 
4108  if ( !mDirtyGeos )
4109  {
4110  // No need to convert again
4111  return true;
4112  }
4113 
4114  if ( mGeos )
4115  {
4116  GEOSGeom_destroy( mGeos );
4117  mGeos = 0;
4118  }
4119 
4120  // this probably shouldn't return true
4121  if ( !mGeometry )
4122  {
4123  // no WKB => no GEOS
4124  mDirtyGeos = false;
4125  return true;
4126  }
4127 
4128  double *x;
4129  double *y;
4130  int *nPoints;
4131  int *numRings;
4132  int *numPolygons;
4133  int numLineStrings;
4134  int idx, jdx, kdx;
4135  unsigned char *ptr;
4136  QgsPoint pt;
4137  QGis::WkbType wkbtype;
4138  bool hasZValue = false;
4139 
4140  //wkbtype = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4];
4141  memcpy( &wkbtype, &( mGeometry[1] ), sizeof( int ) );
4142 
4143  try
4144  {
4145  switch ( wkbtype )
4146  {
4147  case QGis::WKBPoint25D:
4148  case QGis::WKBPoint:
4149  {
4150  x = ( double * )( mGeometry + 5 );
4151  y = ( double * )( mGeometry + 5 + sizeof( double ) );
4152 
4153  mGeos = createGeosPoint( QgsPoint( *x, *y ) );
4154  mDirtyGeos = false;
4155  break;
4156  }
4157 
4159  hasZValue = true;
4160  case QGis::WKBMultiPoint:
4161  {
4162  QVector<GEOSGeometry *> points;
4163 
4164  ptr = mGeometry + 5;
4165  nPoints = ( int * ) ptr;
4166  ptr = mGeometry + 1 + 2 * sizeof( int );
4167  for ( idx = 0; idx < *nPoints; idx++ )
4168  {
4169  ptr += ( 1 + sizeof( int ) );
4170  x = ( double * ) ptr;
4171  ptr += sizeof( double );
4172  y = ( double * ) ptr;
4173  ptr += sizeof( double );
4174  if ( hasZValue )
4175  {
4176  ptr += sizeof( double );
4177  }
4178  points << createGeosPoint( QgsPoint( *x, *y ) );
4179  }
4180  mGeos = createGeosCollection( GEOS_MULTIPOINT, points );
4181  mDirtyGeos = false;
4182  break;
4183  }
4184 
4186  hasZValue = true;
4187  case QGis::WKBLineString:
4188  {
4189  QgsDebugMsgLevel( "Linestring found", 3 );
4190 
4191  QgsPolyline sequence;
4192 
4193  ptr = mGeometry + 5;
4194  nPoints = ( int * ) ptr;
4195  ptr = mGeometry + 1 + 2 * sizeof( int );
4196  for ( idx = 0; idx < *nPoints; idx++ )
4197  {
4198  x = ( double * ) ptr;
4199  ptr += sizeof( double );
4200  y = ( double * ) ptr;
4201  ptr += sizeof( double );
4202  if ( hasZValue )
4203  {
4204  ptr += sizeof( double );
4205  }
4206 
4207  sequence << QgsPoint( *x, *y );
4208  }
4209  mDirtyGeos = false;
4210  mGeos = createGeosLineString( sequence );
4211  break;
4212  }
4213 
4215  hasZValue = true;
4217  {
4218  QVector<GEOSGeometry*> lines;
4219  numLineStrings = ( int )( mGeometry[5] );
4220  ptr = ( mGeometry + 9 );
4221  for ( jdx = 0; jdx < numLineStrings; jdx++ )
4222  {
4223  QgsPolyline sequence;
4224 
4225  // each of these is a wbklinestring so must handle as such
4226  ptr += 5; // skip type since we know its 2
4227  nPoints = ( int * ) ptr;
4228  ptr += sizeof( int );
4229  for ( idx = 0; idx < *nPoints; idx++ )
4230  {
4231  x = ( double * ) ptr;
4232  ptr += sizeof( double );
4233  y = ( double * ) ptr;
4234  ptr += sizeof( double );
4235  if ( hasZValue )
4236  {
4237  ptr += sizeof( double );
4238  }
4239  sequence << QgsPoint( *x, *y );
4240  }
4241  lines << createGeosLineString( sequence );
4242  }
4243  mGeos = createGeosCollection( GEOS_MULTILINESTRING, lines );
4244  mDirtyGeos = false;
4245  break;
4246  }
4247 
4248  case QGis::WKBPolygon25D:
4249  hasZValue = true;
4250  case QGis::WKBPolygon:
4251  {
4252  QgsDebugMsgLevel( "Polygon found", 3 );
4253 
4254  // get number of rings in the polygon
4255  numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
4256  ptr = mGeometry + 1 + 2 * sizeof( int );
4257 
4258  QVector<GEOSGeometry*> rings;
4259 
4260  for ( idx = 0; idx < *numRings; idx++ )
4261  {
4262  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4263 
4264  QgsPolyline sequence;
4265 
4266  // get number of points in the ring
4267  nPoints = ( int * ) ptr;
4268  ptr += 4;
4269  for ( jdx = 0; jdx < *nPoints; jdx++ )
4270  {
4271  // add points to a point array for drawing the polygon
4272  x = ( double * ) ptr;
4273  ptr += sizeof( double );
4274  y = ( double * ) ptr;
4275  ptr += sizeof( double );
4276  if ( hasZValue )
4277  {
4278  ptr += sizeof( double );
4279  }
4280  sequence << QgsPoint( *x, *y );
4281  }
4282 
4283  rings << createGeosLinearRing( sequence );
4284  }
4285  mGeos = createGeosPolygon( rings );
4286  mDirtyGeos = false;
4287  break;
4288  }
4289 
4291  hasZValue = true;
4292  case QGis::WKBMultiPolygon:
4293  {
4294  QgsDebugMsgLevel( "Multipolygon found", 3 );
4295 
4296  QVector<GEOSGeometry*> polygons;
4297 
4298  // get the number of polygons
4299  ptr = mGeometry + 5;
4300  numPolygons = ( int * ) ptr;
4301  ptr = mGeometry + 9;
4302  for ( kdx = 0; kdx < *numPolygons; kdx++ )
4303  {
4304  //QgsDebugMsg("Polygon nr: "+QString::number(kdx));
4305  QVector<GEOSGeometry*> rings;
4306 
4307  //skip the endian and mGeometry type info and
4308  // get number of rings in the polygon
4309  ptr += 5;
4310  numRings = ( int * ) ptr;
4311  ptr += 4;
4312  for ( idx = 0; idx < *numRings; idx++ )
4313  {
4314  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4315 
4316  QgsPolyline sequence;
4317 
4318  // get number of points in the ring
4319  nPoints = ( int * ) ptr;
4320  ptr += 4;
4321  for ( jdx = 0; jdx < *nPoints; jdx++ )
4322  {
4323  // add points to a point array for drawing the polygon
4324  x = ( double * ) ptr;
4325  ptr += sizeof( double );
4326  y = ( double * ) ptr;
4327  ptr += sizeof( double );
4328  if ( hasZValue )
4329  {
4330  ptr += sizeof( double );
4331  }
4332  sequence << QgsPoint( *x, *y );
4333  }
4334 
4335  rings << createGeosLinearRing( sequence );
4336  }
4337 
4338  polygons << createGeosPolygon( rings );
4339  }
4340  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
4341  mDirtyGeos = false;
4342  break;
4343  }
4344 
4345  default:
4346  return false;
4347  }
4348  }
4349  CATCH_GEOS( false )
4350 
4351  return true;
4352 }
4353 
4355 {
4356  //QgsDebugMsg("entered.");
4357  if ( !mDirtyWkb )
4358  {
4359  // No need to convert again
4360  return true;
4361  }
4362 
4363  // clear the WKB, ready to replace with the new one
4364  if ( mGeometry )
4365  {
4366  delete [] mGeometry;
4367  mGeometry = 0;
4368  }
4369 
4370  if ( !mGeos )
4371  {
4372  // GEOS is null, therefore WKB is null.
4373  mDirtyWkb = false;
4374  return true;
4375  }
4376 
4377  // set up byteOrder
4378  char byteOrder = QgsApplication::endian();
4379 
4380  switch ( GEOSGeomTypeId( mGeos ) )
4381  {
4382  case GEOS_POINT: // a point
4383  {
4384  mGeometrySize = 1 + // sizeof(byte)
4385  4 + // sizeof(uint32)
4386  2 * sizeof( double );
4387  mGeometry = new unsigned char[mGeometrySize];
4388 
4389  // assign byteOrder
4390  memcpy( mGeometry, &byteOrder, 1 );
4391 
4392  // assign wkbType
4393  int wkbType = QGis::WKBPoint;
4394  memcpy( mGeometry + 1, &wkbType, 4 );
4395 
4396  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos );
4397 
4398  double x, y;
4399  GEOSCoordSeq_getX( cs, 0, &x );
4400  GEOSCoordSeq_getY( cs, 0, &y );
4401 
4402  memcpy( mGeometry + 5, &x, sizeof( double ) );
4403  memcpy( mGeometry + 13, &y, sizeof( double ) );
4404 
4405  mDirtyWkb = false;
4406  return true;
4407  } // case GEOS_GEOM::GEOS_POINT
4408 
4409  case GEOS_LINESTRING: // a linestring
4410  {
4411  //QgsDebugMsg("Got a geos::GEOS_LINESTRING.");
4412 
4413  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos );
4414  unsigned int numPoints;
4415  GEOSCoordSeq_getSize( cs, &numPoints );
4416 
4417  // allocate some space for the WKB
4418  mGeometrySize = 1 + // sizeof(byte)
4419  4 + // sizeof(uint32)
4420  4 + // sizeof(uint32)
4421  (( sizeof( double ) +
4422  sizeof( double ) ) * numPoints );
4423 
4424  mGeometry = new unsigned char[mGeometrySize];
4425 
4426  unsigned char* ptr = mGeometry;
4427 
4428  // assign byteOrder
4429  memcpy( ptr, &byteOrder, 1 );
4430  ptr += 1;
4431 
4432  // assign wkbType
4434  memcpy( ptr, &wkbType, 4 );
4435  ptr += 4;
4436 
4437  // assign numPoints
4438  memcpy( ptr, &numPoints, 4 );
4439  ptr += 4;
4440 
4441  const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( mGeos );
4442 
4443  // assign points
4444  for ( unsigned int n = 0; n < numPoints; n++ )
4445  {
4446  // assign x
4447  GEOSCoordSeq_getX( sequence, n, ( double * )ptr );
4448  ptr += sizeof( double );
4449 
4450  // assign y
4451  GEOSCoordSeq_getY( sequence, n, ( double * )ptr );
4452  ptr += sizeof( double );
4453  }
4454 
4455  mDirtyWkb = false;
4456  return true;
4457 
4458  // TODO: Deal with endian-ness
4459  } // case GEOS_GEOM::GEOS_LINESTRING
4460 
4461  case GEOS_LINEARRING: // a linear ring (linestring with 1st point == last point)
4462  {
4463  // TODO
4464  break;
4465  } // case GEOS_GEOM::GEOS_LINEARRING
4466 
4467  case GEOS_POLYGON: // a polygon
4468  {
4469  int geometrySize;
4470  unsigned int nPointsInRing = 0;
4471 
4472  //first calculate the geometry size
4473  geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of rings
4474  const GEOSGeometry *theRing = GEOSGetExteriorRing( mGeos );
4475  if ( theRing )
4476  {
4477  geometrySize += sizeof( int );
4478  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4479  }
4480  for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); ++i )
4481  {
4482  geometrySize += sizeof( int ); //number of points in ring
4483  theRing = GEOSGetInteriorRingN( mGeos, i );
4484  if ( theRing )
4485  {
4486  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4487  }
4488  }
4489 
4490  mGeometry = new unsigned char[geometrySize];
4491  mGeometrySize = geometrySize;
4492 
4493  //then fill the geometry itself into the wkb
4494  int position = 0;
4495  // assign byteOrder
4496  memcpy( mGeometry, &byteOrder, 1 );
4497  position += 1;
4498  int wkbtype = QGis::WKBPolygon;
4499  memcpy( &mGeometry[position], &wkbtype, sizeof( int ) );
4500  position += sizeof( int );
4501  int nRings = GEOSGetNumInteriorRings( mGeos ) + 1;
4502  memcpy( &mGeometry[position], &nRings, sizeof( int ) );
4503  position += sizeof( int );
4504 
4505  //exterior ring first
4506  theRing = GEOSGetExteriorRing( mGeos );
4507  if ( theRing )
4508  {
4509  nPointsInRing = getNumGeosPoints( theRing );
4510  memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) );
4511  position += sizeof( int );
4512 
4513  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
4514  unsigned int n;
4515  GEOSCoordSeq_getSize( cs, &n );
4516 
4517  for ( unsigned int j = 0; j < n; ++j )
4518  {
4519  GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] );
4520  position += sizeof( double );
4521  GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] );
4522  position += sizeof( double );
4523  }
4524  }
4525 
4526  //interior rings after
4527  for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); i++ )
4528  {
4529  theRing = GEOSGetInteriorRingN( mGeos, i );
4530 
4531  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
4532  GEOSCoordSeq_getSize( cs, &nPointsInRing );
4533 
4534  memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) );
4535  position += sizeof( int );
4536 
4537  for ( unsigned int j = 0; j < nPointsInRing; j++ )
4538  {
4539  GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] );
4540  position += sizeof( double );
4541  GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] );
4542  position += sizeof( double );
4543  }
4544  }
4545  mDirtyWkb = false;
4546  return true;
4547  } // case GEOS_GEOM::GEOS_POLYGON
4548  break;
4549 
4550  case GEOS_MULTIPOINT: // a collection of points
4551  {
4552  // determine size of geometry
4553  int geometrySize = 1 + 2 * sizeof( int );
4554  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4555  {
4556  geometrySize += 1 + sizeof( int ) + 2 * sizeof( double );
4557  }
4558 
4559  mGeometry = new unsigned char[geometrySize];
4560  mGeometrySize = geometrySize;
4561  int wkbPosition = 0; //current position in the byte array
4562 
4563  memcpy( mGeometry, &byteOrder, 1 );
4564  wkbPosition += 1;
4565  int wkbtype = QGis::WKBMultiPoint;
4566  memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
4567  wkbPosition += sizeof( int );
4568  int numPoints = GEOSGetNumGeometries( mGeos );
4569  memcpy( &mGeometry[wkbPosition], &numPoints, sizeof( int ) );
4570  wkbPosition += sizeof( int );
4571 
4572  int pointType = QGis::WKBPoint;
4573  const GEOSGeometry *currentPoint = 0;
4574 
4575  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4576  {
4577  //copy endian and point type
4578  memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
4579  wkbPosition += 1;
4580  memcpy( &mGeometry[wkbPosition], &pointType, sizeof( int ) );
4581  wkbPosition += sizeof( int );
4582 
4583  currentPoint = GEOSGetGeometryN( mGeos, i );
4584 
4585  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( currentPoint );
4586 
4587  GEOSCoordSeq_getX( cs, 0, ( double* )&mGeometry[wkbPosition] );
4588  wkbPosition += sizeof( double );
4589  GEOSCoordSeq_getY( cs, 0, ( double* )&mGeometry[wkbPosition] );
4590  wkbPosition += sizeof( double );
4591  }
4592  mDirtyWkb = false;
4593  return true;
4594  } // case GEOS_GEOM::GEOS_MULTIPOINT
4595 
4596  case GEOS_MULTILINESTRING: // a collection of linestrings
4597  {
4598  // determine size of geometry
4599  int geometrySize = 1 + 2 * sizeof( int );
4600  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4601  {
4602  geometrySize += 1 + 2 * sizeof( int );
4603  geometrySize += getNumGeosPoints( GEOSGetGeometryN( mGeos, i ) ) * 2 * sizeof( double );
4604  }
4605 
4606  mGeometry = new unsigned char[geometrySize];
4607  mGeometrySize = geometrySize;
4608  int wkbPosition = 0; //current position in the byte array
4609 
4610  memcpy( mGeometry, &byteOrder, 1 );
4611  wkbPosition += 1;
4612  int wkbtype = QGis::WKBMultiLineString;
4613  memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
4614  wkbPosition += sizeof( int );
4615  int numLines = GEOSGetNumGeometries( mGeos );
4616  memcpy( &mGeometry[wkbPosition], &numLines, sizeof( int ) );
4617  wkbPosition += sizeof( int );
4618 
4619  //loop over lines
4620  int lineType = QGis::WKBLineString;
4621  const GEOSCoordSequence *cs = 0;
4622  unsigned int lineSize;
4623 
4624  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4625  {
4626  //endian and type WKBLineString
4627  memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
4628  wkbPosition += 1;
4629  memcpy( &mGeometry[wkbPosition], &lineType, sizeof( int ) );
4630  wkbPosition += sizeof( int );
4631 
4632  cs = GEOSGeom_getCoordSeq( GEOSGetGeometryN( mGeos, i ) );
4633 
4634  //line size
4635  GEOSCoordSeq_getSize( cs, &lineSize );
4636  memcpy( &mGeometry[wkbPosition], &lineSize, sizeof( int ) );
4637  wkbPosition += sizeof( int );
4638 
4639  //vertex coordinates
4640  for ( unsigned int j = 0; j < lineSize; ++j )
4641  {
4642  GEOSCoordSeq_getX( cs, j, ( double* )&mGeometry[wkbPosition] );
4643  wkbPosition += sizeof( double );
4644  GEOSCoordSeq_getY( cs, j, ( double* )&mGeometry[wkbPosition] );
4645  wkbPosition += sizeof( double );
4646  }
4647  }
4648  mDirtyWkb = false;
4649  return true;
4650  } // case GEOS_GEOM::GEOS_MULTILINESTRING
4651 
4652  case GEOS_MULTIPOLYGON: // a collection of polygons
4653  {
4654  //first determine size of geometry
4655  int geometrySize = 1 + ( 2 * sizeof( int ) ); //endian, type, number of polygons
4656  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4657  {
4658  const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i );
4659  geometrySize += 1 + 2 * sizeof( int ); //endian, type, number of rings
4660  //exterior ring
4661  geometrySize += sizeof( int ); //number of points in exterior ring
4662  const GEOSGeometry *exRing = GEOSGetExteriorRing( thePoly );
4663  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( exRing );
4664 
4665  const GEOSGeometry *intRing = 0;
4666  for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ )
4667  {
4668  geometrySize += sizeof( int ); //number of points in ring
4669  intRing = GEOSGetInteriorRingN( thePoly, j );
4670  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( intRing );
4671  }
4672  }
4673 
4674  mGeometry = new unsigned char[geometrySize];
4675  mGeometrySize = geometrySize;
4676  int wkbPosition = 0; //current position in the byte array
4677 
4678  memcpy( mGeometry, &byteOrder, 1 );
4679  wkbPosition += 1;
4680  int wkbtype = QGis::WKBMultiPolygon;
4681  memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
4682  wkbPosition += sizeof( int );
4683  int numPolygons = GEOSGetNumGeometries( mGeos );
4684  memcpy( &mGeometry[wkbPosition], &numPolygons, sizeof( int ) );
4685  wkbPosition += sizeof( int );
4686 
4687  //loop over polygons
4688  for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
4689  {
4690  const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i );
4691  memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
4692  wkbPosition += 1;
4693  int polygonType = QGis::WKBPolygon;
4694  memcpy( &mGeometry[wkbPosition], &polygonType, sizeof( int ) );
4695  wkbPosition += sizeof( int );
4696  int numRings = GEOSGetNumInteriorRings( thePoly ) + 1;
4697  memcpy( &mGeometry[wkbPosition], &numRings, sizeof( int ) );
4698  wkbPosition += sizeof( int );
4699 
4700  //exterior ring
4701  const GEOSGeometry *theRing = GEOSGetExteriorRing( thePoly );
4702  int nPointsInRing = getNumGeosPoints( theRing );
4703  memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) );
4704  wkbPosition += sizeof( int );
4705  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
4706 
4707  for ( int k = 0; k < nPointsInRing; ++k )
4708  {
4709  GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] );
4710  wkbPosition += sizeof( double );
4711  GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] );
4712  wkbPosition += sizeof( double );
4713  }
4714 
4715  //interior rings
4716  for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ )
4717  {
4718  theRing = GEOSGetInteriorRingN( thePoly, j );
4719  nPointsInRing = getNumGeosPoints( theRing );
4720  memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) );
4721  wkbPosition += sizeof( int );
4722  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
4723 
4724  for ( int k = 0; k < nPointsInRing; ++k )
4725  {
4726  GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] );
4727  wkbPosition += sizeof( double );
4728  GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] );
4729  wkbPosition += sizeof( double );
4730  }
4731  }
4732  }
4733  mDirtyWkb = false;
4734  return true;
4735  } // case GEOS_GEOM::GEOS_MULTIPOLYGON
4736 
4737  case GEOS_GEOMETRYCOLLECTION: // a collection of heterogeneus geometries
4738  {
4739  // TODO
4740  QgsDebugMsg( "geometry collection - not supported" );
4741  break;
4742  } // case GEOS_GEOM::GEOS_GEOMETRYCOLLECTION
4743 
4744  } // switch (mGeos->getGeometryTypeId())
4745 
4746  return false;
4747 }
4748 
4750 {
4751  if ( !mGeometry )
4752  {
4753  return false;
4754  }
4755 
4756  QGis::WkbType geomType = wkbType();
4757 
4758  if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D ||
4759  geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D ||
4760  geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown )
4761  {
4762  return false; //no need to convert
4763  }
4764 
4765  int newGeomSize = mGeometrySize + 1 + 2 * sizeof( int ); //endian: 1, multitype: sizeof(int), number of geometries: sizeof(int)
4766  unsigned char* newGeometry = new unsigned char[newGeomSize];
4767 
4768  int currentWkbPosition = 0;
4769  //copy endian
4770  char byteOrder = QgsApplication::endian();
4771  memcpy( &newGeometry[currentWkbPosition], &byteOrder, 1 );
4772  currentWkbPosition += 1;
4773 
4774  //copy wkbtype
4775  //todo
4776  QGis::WkbType newMultiType;
4777  switch ( geomType )
4778  {
4779  case QGis::WKBPoint:
4780  newMultiType = QGis::WKBMultiPoint;
4781  break;
4782  case QGis::WKBPoint25D:
4783  newMultiType = QGis::WKBMultiPoint25D;
4784  break;
4785  case QGis::WKBLineString:
4786  newMultiType = QGis::WKBMultiLineString;
4787  break;
4789  newMultiType = QGis::WKBMultiLineString25D;
4790  break;
4791  case QGis::WKBPolygon:
4792  newMultiType = QGis::WKBMultiPolygon;
4793  break;
4794  case QGis::WKBPolygon25D:
4795  newMultiType = QGis::WKBMultiPolygon25D;
4796  break;
4797  default:
4798  delete newGeometry;
4799  return false;
4800  }
4801  memcpy( &newGeometry[currentWkbPosition], &newMultiType, sizeof( int ) );
4802  currentWkbPosition += sizeof( int );
4803 
4804  //copy number of geometries
4805  int nGeometries = 1;
4806  memcpy( &newGeometry[currentWkbPosition], &nGeometries, sizeof( int ) );
4807  currentWkbPosition += sizeof( int );
4808 
4809  //copy the existing single geometry
4810  memcpy( &newGeometry[currentWkbPosition], mGeometry, mGeometrySize );
4811 
4812  delete [] mGeometry;
4813  mGeometry = newGeometry;
4814  mGeometrySize = newGeomSize;
4815  mDirtyGeos = true;
4816  return true;
4817 }
4818 
4819 void QgsGeometry::translateVertex( int& wkbPosition, double dx, double dy, bool hasZValue )
4820 {
4821  double x, y, translated_x, translated_y;
4822 
4823  //x-coordinate
4824  x = *(( double * )( &( mGeometry[wkbPosition] ) ) );
4825  translated_x = x + dx;
4826  memcpy( &( mGeometry[wkbPosition] ), &translated_x, sizeof( double ) );
4827  wkbPosition += sizeof( double );
4828 
4829  //y-coordinate
4830  y = *(( double * )( &( mGeometry[wkbPosition] ) ) );
4831  translated_y = y + dy;
4832  memcpy( &( mGeometry[wkbPosition] ), &translated_y, sizeof( double ) );
4833  wkbPosition += sizeof( double );
4834 
4835  if ( hasZValue )
4836  {
4837  wkbPosition += sizeof( double );
4838  }
4839 }
4840 
4841 void QgsGeometry::transformVertex( int& wkbPosition, const QgsCoordinateTransform& ct, bool hasZValue )
4842 {
4843  double x, y, z;
4844 
4845 
4846  x = *(( double * )( &( mGeometry[wkbPosition] ) ) );
4847  y = *(( double * )( &( mGeometry[wkbPosition + sizeof( double )] ) ) );
4848  z = 0.0; // Ignore Z for now.
4849 
4850  ct.transformInPlace( x, y, z );
4851 
4852  // new x-coordinate
4853  memcpy( &( mGeometry[wkbPosition] ), &x, sizeof( double ) );
4854  wkbPosition += sizeof( double );
4855 
4856  // new y-coordinate
4857  memcpy( &( mGeometry[wkbPosition] ), &y, sizeof( double ) );
4858  wkbPosition += sizeof( double );
4859 
4860  if ( hasZValue )
4861  {
4862  wkbPosition += sizeof( double );
4863  }
4864 }
4865 
4866 int QgsGeometry::splitLinearGeometry( GEOSGeometry *splitLine, QList<QgsGeometry*>& newGeometries )
4867 {
4868  if ( !splitLine )
4869  {
4870  return 2;
4871  }
4872 
4873  if ( !mGeos || mDirtyGeos )
4874  {
4875  if ( !exportWkbToGeos() )
4876  return 5;
4877  }
4878 
4879  //first test if linestring intersects geometry. If not, return straight away
4880  if ( !GEOSIntersects( splitLine, mGeos ) )
4881  {
4882  return 1;
4883  }
4884 
4885  //first union all the polygon rings together (to get them noded, see JTS developer guide)
4886  GEOSGeometry* nodedGeometry = nodeGeometries( splitLine, mGeos );
4887  if ( !nodedGeometry )
4888  {
4889  return 3; //an error occured during noding
4890  }
4891 
4892  GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
4893  if ( !mergedLines )
4894  {
4895  GEOSGeom_destroy( nodedGeometry );
4896  return 4;
4897  }
4898 
4899  QVector<GEOSGeometry*> testedGeometries;
4900 
4901  for ( int i = 0; i < GEOSGetNumGeometries( mergedLines ); i++ )
4902  {
4903  const GEOSGeometry *testing = GEOSGetGeometryN( mergedLines, i );
4904  if ( lineContainedInLine( testing, mGeos ) == 1 )
4905  {
4906  testedGeometries << GEOSGeom_clone( testing );
4907  }
4908  }
4909 
4910  mergeGeometriesMultiTypeSplit( testedGeometries );
4911 
4912  if ( testedGeometries.size() > 0 )
4913  {
4914  GEOSGeom_destroy( mGeos );
4915  mGeos = testedGeometries[0];
4916  mDirtyWkb = true;
4917  }
4918 
4919  for ( int i = 1; i < testedGeometries.size(); ++i )
4920  newGeometries << fromGeosGeom( testedGeometries[i] );
4921 
4922  GEOSGeom_destroy( nodedGeometry );
4923  GEOSGeom_destroy( mergedLines );
4924  return 0;
4925 }
4926 
4927 int QgsGeometry::splitPolygonGeometry( GEOSGeometry* splitLine, QList<QgsGeometry*>& newGeometries )
4928 {
4929  if ( !splitLine )
4930  {
4931  return 2;
4932  }
4933 
4934  if ( !mGeos || mDirtyGeos )
4935  {
4936  if ( !exportWkbToGeos() )
4937  return 5;
4938  }
4939 
4940  //first test if linestring intersects geometry. If not, return straight away
4941  if ( !GEOSIntersects( splitLine, mGeos ) )
4942  {
4943  return 1;
4944  }
4945 
4946  //first union all the polygon rings together (to get them noded, see JTS developer guide)
4947  GEOSGeometry *nodedGeometry = nodeGeometries( splitLine, mGeos );
4948  if ( !nodedGeometry )
4949  {
4950  return 2; //an error occured during noding
4951  }
4952 
4953 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
4954  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=1)))
4955  GEOSGeometry *cutEdges = GEOSPolygonizer_getCutEdges( &nodedGeometry, 1 );
4956  if ( cutEdges )
4957  {
4958  if ( numberOfGeometries( cutEdges ) > 0 )
4959  {
4960  GEOSGeom_destroy( cutEdges );
4961  GEOSGeom_destroy( nodedGeometry );
4962  return 3;
4963  }
4964 
4965  GEOSGeom_destroy( cutEdges );
4966  }
4967 #endif
4968 
4969  GEOSGeometry *polygons = GEOSPolygonize( &nodedGeometry, 1 );
4970  if ( !polygons || numberOfGeometries( polygons ) == 0 )
4971  {
4972  if ( polygons )
4973  GEOSGeom_destroy( polygons );
4974 
4975  GEOSGeom_destroy( nodedGeometry );
4976 
4977  return 4;
4978  }
4979 
4980  GEOSGeom_destroy( nodedGeometry );
4981 
4982  //test every polygon if contained in original geometry
4983  //include in result if yes
4984  QVector<GEOSGeometry*> testedGeometries;
4985  GEOSGeometry *intersectGeometry = 0;
4986 
4987  //ratio intersect geometry / geometry. This should be close to 1
4988  //if the polygon belongs to the input geometry
4989 
4990  for ( int i = 0; i < numberOfGeometries( polygons ); i++ )
4991  {
4992  const GEOSGeometry *polygon = GEOSGetGeometryN( polygons, i );
4993  intersectGeometry = GEOSIntersection( mGeos, polygon );
4994  if ( !intersectGeometry )
4995  {
4996  QgsDebugMsg( "intersectGeometry is NULL" );
4997  continue;
4998  }
4999 
5000  double intersectionArea;
5001  GEOSArea( intersectGeometry, &intersectionArea );
5002 
5003  double polygonArea;
5004  GEOSArea( polygon, &polygonArea );
5005 
5006  const double areaRatio = intersectionArea / polygonArea;
5007  if ( areaRatio > 0.99 && areaRatio < 1.01 )
5008  testedGeometries << GEOSGeom_clone( polygon );
5009 
5010  GEOSGeom_destroy( intersectGeometry );
5011  }
5012 
5013  bool splitDone = true;
5014  int nGeometriesThis = numberOfGeometries( mGeos ); //original number of geometries
5015  if ( testedGeometries.size() == nGeometriesThis )
5016  {
5017  splitDone = false;
5018  }
5019 
5020  mergeGeometriesMultiTypeSplit( testedGeometries );
5021 
5022  //no split done, preserve original geometry
5023  if ( !splitDone )
5024  {
5025  for ( int i = 0; i < testedGeometries.size(); ++i )
5026  {
5027  GEOSGeom_destroy( testedGeometries[i] );
5028  }
5029  return 1;
5030  }
5031  else if ( testedGeometries.size() > 0 ) //split successfull
5032  {
5033  GEOSGeom_destroy( mGeos );
5034  mGeos = testedGeometries[0];
5035  mDirtyWkb = true;
5036  }
5037 
5038  for ( int i = 1; i < testedGeometries.size(); ++i )
5039  {
5040  newGeometries << fromGeosGeom( testedGeometries[i] );
5041  }
5042 
5043  GEOSGeom_destroy( polygons );
5044  return 0;
5045 }
5046 
5047 GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos )
5048 {
5049  //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
5050  int nIntersections = 0;
5051  int lastIntersectingRing = -2;
5052  const GEOSGeometry* lastIntersectingGeom = 0;
5053 
5054  int nRings = GEOSGetNumInteriorRings( polygon );
5055  if ( nRings < 0 )
5056  {
5057  return 0;
5058  }
5059 
5060  //does outer ring intersect?
5061  const GEOSGeometry* outerRing = GEOSGetExteriorRing( polygon );
5062  if ( GEOSIntersects( outerRing, reshapeLineGeos ) == 1 )
5063  {
5064  ++nIntersections;
5065  lastIntersectingRing = -1;
5066  lastIntersectingGeom = outerRing;
5067  }
5068 
5069  //do inner rings intersect?
5070  const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
5071 
5072  try
5073  {
5074  for ( int i = 0; i < nRings; ++i )
5075  {
5076  innerRings[i] = GEOSGetInteriorRingN( polygon, i );
5077  if ( GEOSIntersects( innerRings[i], reshapeLineGeos ) == 1 )
5078  {
5079  ++nIntersections;
5080  lastIntersectingRing = i;
5081  lastIntersectingGeom = innerRings[i];
5082  }
5083  }
5084  }
5085  catch ( GEOSException &e )
5086  {
5087  Q_UNUSED( e );
5088  nIntersections = 0;
5089  }
5090 
5091  if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
5092  {
5093  delete [] innerRings;
5094  return 0;
5095  }
5096 
5097  //we have one intersecting ring, let's try to reshape it
5098  GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos );
5099  if ( !reshapeResult )
5100  {
5101  delete [] innerRings;
5102  return 0;
5103  }
5104 
5105  //if reshaping took place, we need to reassemble the polygon and its rings
5106  GEOSGeometry* newRing = 0;
5107  const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq( reshapeResult );
5108  GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone( reshapeSequence );
5109 
5110  GEOSGeom_destroy( reshapeResult );
5111 
5112  newRing = GEOSGeom_createLinearRing( newCoordSequence );
5113  if ( !newRing )
5114  {
5115  delete [] innerRings;
5116  return 0;
5117  }
5118 
5119 
5120  GEOSGeometry* newOuterRing = 0;
5121  if ( lastIntersectingRing == -1 )
5122  {
5123  newOuterRing = newRing;
5124  }
5125  else
5126  {
5127  newOuterRing = GEOSGeom_clone( outerRing );
5128  }
5129 
5130  //check if all the rings are still inside the outer boundary
5131  QList<GEOSGeometry*> ringList;
5132  if ( nRings > 0 )
5133  {
5134  GEOSGeometry* outerRingPoly = GEOSGeom_createPolygon( GEOSGeom_clone( newOuterRing ), 0, 0 );
5135  if ( outerRingPoly )
5136  {
5137  GEOSGeometry* currentRing = 0;
5138  for ( int i = 0; i < nRings; ++i )
5139  {
5140  if ( lastIntersectingRing == i )
5141  {
5142  currentRing = newRing;
5143  }
5144  else
5145  {
5146  currentRing = GEOSGeom_clone( innerRings[i] );
5147  }
5148 
5149  //possibly a ring is no longer contained in the result polygon after reshape
5150  if ( GEOSContains( outerRingPoly, currentRing ) == 1 )
5151  {
5152  ringList.push_back( currentRing );
5153  }
5154  else
5155  {
5156  GEOSGeom_destroy( currentRing );
5157  }
5158  }
5159  }
5160  GEOSGeom_destroy( outerRingPoly );
5161  }
5162 
5163  GEOSGeometry** newInnerRings = new GEOSGeometry*[ringList.size()];
5164  for ( int i = 0; i < ringList.size(); ++i )
5165  {
5166  newInnerRings[i] = ringList.at( i );
5167  }
5168 
5169  delete [] innerRings;
5170 
5171  GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon( newOuterRing, newInnerRings, ringList.size() );
5172  delete[] newInnerRings;
5173  if ( !reshapedPolygon )
5174  {
5175  return 0;
5176  }
5177  return reshapedPolygon;
5178 }
5179 
5180 GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos )
5181 {
5182  if ( !line || !reshapeLineGeos )
5183  {
5184  return 0;
5185  }
5186 
5187  bool atLeastTwoIntersections;
5188 
5189  try
5190  {
5191  //make sure there are at least two intersection between line and reshape geometry
5192  GEOSGeometry* intersectGeom = GEOSIntersection( line, reshapeLineGeos );
5193  atLeastTwoIntersections = ( GEOSGeomTypeId( intersectGeom ) == GEOS_MULTIPOINT && GEOSGetNumGeometries( intersectGeom ) > 1 );
5194  GEOSGeom_destroy( intersectGeom );
5195  }
5196  catch ( GEOSException &e )
5197  {
5198  Q_UNUSED( e );
5199  atLeastTwoIntersections = false;
5200  }
5201 
5202  if ( !atLeastTwoIntersections )
5203  {
5204  return 0;
5205  }
5206 
5207  //begin and end point of original line
5208  const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq( line );
5209  if ( !lineCoordSeq )
5210  {
5211  return 0;
5212  }
5213  unsigned int lineCoordSeqSize;
5214  if ( GEOSCoordSeq_getSize( lineCoordSeq, &lineCoordSeqSize ) == 0 )
5215  {
5216  return 0;
5217  }
5218  if ( lineCoordSeqSize < 2 )
5219  {
5220  return 0;
5221  }
5222  //first and last vertex of line
5223  double x1, y1, x2, y2;
5224  GEOSCoordSeq_getX( lineCoordSeq, 0, &x1 );
5225  GEOSCoordSeq_getY( lineCoordSeq, 0, &y1 );
5226  GEOSCoordSeq_getX( lineCoordSeq, lineCoordSeqSize - 1, &x2 );
5227  GEOSCoordSeq_getY( lineCoordSeq, lineCoordSeqSize - 1, &y2 );
5228  GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) );
5229  GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) );
5230 
5231  bool isRing = false;
5232  if ( GEOSGeomTypeId( line ) == GEOS_LINEARRING || GEOSEquals( beginLineVertex, endLineVertex ) == 1 )
5233  {
5234  isRing = true;
5235  }
5236 
5237 //node line and reshape line
5238  GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line );
5239  if ( !nodedGeometry )
5240  {
5241  GEOSGeom_destroy( beginLineVertex );
5242  GEOSGeom_destroy( endLineVertex );
5243  return 0;
5244  }
5245 
5246  //and merge them together
5247  GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
5248  GEOSGeom_destroy( nodedGeometry );
5249  if ( !mergedLines )
5250  {
5251  GEOSGeom_destroy( beginLineVertex );
5252  GEOSGeom_destroy( endLineVertex );
5253  return 0;
5254  }
5255 
5256  int numMergedLines = GEOSGetNumGeometries( mergedLines );
5257  if ( numMergedLines < 2 ) //some special cases. Normally it is >2
5258  {
5259  GEOSGeom_destroy( beginLineVertex );
5260  GEOSGeom_destroy( endLineVertex );
5261  if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
5262  {
5263  return GEOSGeom_clone( reshapeLineGeos );
5264  }
5265  else
5266  {
5267  return 0;
5268  }
5269  }
5270 
5271  QList<GEOSGeometry*> resultLineParts; //collection with the line segments that will be contained in result
5272  QList<GEOSGeometry*> probableParts; //parts where we can decide on inclusion only after going through all the candidates
5273 
5274  for ( int i = 0; i < numMergedLines; ++i )
5275  {
5276  const GEOSGeometry* currentGeom;
5277 
5278  currentGeom = GEOSGetGeometryN( mergedLines, i );
5279  const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq( currentGeom );
5280  unsigned int currentCoordSeqSize;
5281  GEOSCoordSeq_getSize( currentCoordSeq, &currentCoordSeqSize );
5282  if ( currentCoordSeqSize < 2 )
5283  {
5284  continue;
5285  }
5286 
5287  //get the two endpoints of the current line merge result
5288  double xBegin, xEnd, yBegin, yEnd;
5289  GEOSCoordSeq_getX( currentCoordSeq, 0, &xBegin );
5290  GEOSCoordSeq_getY( currentCoordSeq, 0, &yBegin );
5291  GEOSCoordSeq_getX( currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
5292  GEOSCoordSeq_getY( currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
5293  GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) );
5294  GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) );
5295 
5296  //check how many endpoints of the line merge result are on the (original) line
5297  int nEndpointsOnOriginalLine = 0;
5298  if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 )
5299  {
5300  nEndpointsOnOriginalLine += 1;
5301  }
5302 
5303  if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 )
5304  {
5305  nEndpointsOnOriginalLine += 1;
5306  }
5307 
5308  //check how many endpoints equal the endpoints of the original line
5309  int nEndpointsSameAsOriginalLine = 0;
5310  if ( GEOSEquals( beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( beginCurrentGeomVertex, endLineVertex ) == 1 )
5311  {
5312  nEndpointsSameAsOriginalLine += 1;
5313  }
5314  if ( GEOSEquals( endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( endCurrentGeomVertex, endLineVertex ) == 1 )
5315  {
5316  nEndpointsSameAsOriginalLine += 1;
5317  }
5318 
5319  //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
5320  bool currentGeomOverlapsOriginalGeom = false;
5321  bool currentGeomOverlapsReshapeLine = false;
5322  if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 )
5323  {
5324  currentGeomOverlapsOriginalGeom = true;
5325  }
5326  if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
5327  {
5328  currentGeomOverlapsReshapeLine = true;
5329  }
5330 
5331 
5332  //logic to decide if this part belongs to the result
5333  if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
5334  {
5335  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
5336  }
5337  //for closed rings, we take one segment from the candidate list
5338  else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
5339  {
5340  probableParts.push_back( GEOSGeom_clone( currentGeom ) );
5341  }
5342  else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
5343  {
5344  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
5345  }
5346  else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
5347  {
5348  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
5349  }
5350  else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
5351  {
5352  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
5353  }
5354 
5355  GEOSGeom_destroy( beginCurrentGeomVertex );
5356  GEOSGeom_destroy( endCurrentGeomVertex );
5357  }
5358 
5359  //add the longest segment from the probable list for rings (only used for polygon rings)
5360  if ( isRing && probableParts.size() > 0 )
5361  {
5362  GEOSGeometry* maxGeom = 0; //the longest geometry in the probabla list
5363  GEOSGeometry* currentGeom = 0;
5364  double maxLength = -DBL_MAX;
5365  double currentLength = 0;
5366  for ( int i = 0; i < probableParts.size(); ++i )
5367  {
5368  currentGeom = probableParts.at( i );
5369  GEOSLength( currentGeom, &currentLength );
5370  if ( currentLength > maxLength )
5371  {
5372  maxLength = currentLength;
5373  GEOSGeom_destroy( maxGeom );
5374  maxGeom = currentGeom;
5375  }
5376  else
5377  {
5378  GEOSGeom_destroy( currentGeom );
5379  }
5380  }
5381  resultLineParts.push_back( maxGeom );
5382  }
5383 
5384  GEOSGeom_destroy( beginLineVertex );
5385  GEOSGeom_destroy( endLineVertex );
5386  GEOSGeom_destroy( mergedLines );
5387 
5388  GEOSGeometry* result = 0;
5389  if ( resultLineParts.size() < 1 )
5390  {
5391  return 0;
5392  }
5393  if ( resultLineParts.size() == 1 ) //the whole result was reshaped
5394  {
5395  result = resultLineParts[0];
5396  }
5397  else //>1
5398  {
5399  GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
5400  for ( int i = 0; i < resultLineParts.size(); ++i )
5401  {
5402  lineArray[i] = resultLineParts[i];
5403  }
5404 
5405  //create multiline from resultLineParts
5406  GEOSGeometry* multiLineGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, lineArray, resultLineParts.size() );
5407  delete [] lineArray;
5408 
5409  //then do a linemerge with the newly combined partstrings
5410  result = GEOSLineMerge( multiLineGeom );
5411  GEOSGeom_destroy( multiLineGeom );
5412  }
5413 
5414  //now test if the result is a linestring. Otherwise something went wrong
5415  if ( GEOSGeomTypeId( result ) != GEOS_LINESTRING )
5416  {
5417  GEOSGeom_destroy( result );
5418  return 0;
5419  }
5420  return result;
5421 }
5422 
5423 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const
5424 {
5425  //Find out the intersection points between splitLineGeos and this geometry.
5426  //These points need to be tested for topological correctness by the calling function
5427  //if topological editing is enabled
5428 
5429  testPoints.clear();
5430  GEOSGeometry* intersectionGeom = GEOSIntersection( mGeos, splitLine );
5431  if ( !intersectionGeom )
5432  {
5433  return 1;
5434  }
5435 
5436  bool simple = false;
5437  int nIntersectGeoms = 1;
5438  if ( GEOSGeomTypeId( intersectionGeom ) == GEOS_LINESTRING || GEOSGeomTypeId( intersectionGeom ) == GEOS_POINT )
5439  {
5440  simple = true;
5441  }
5442 
5443  if ( !simple )
5444  {
5445  nIntersectGeoms = GEOSGetNumGeometries( intersectionGeom );
5446  }
5447 
5448  for ( int i = 0; i < nIntersectGeoms; ++i )
5449  {
5450  const GEOSGeometry* currentIntersectGeom;
5451  if ( simple )
5452  {
5453  currentIntersectGeom = intersectionGeom;
5454  }
5455  else
5456  {
5457  currentIntersectGeom = GEOSGetGeometryN( intersectionGeom, i );
5458  }
5459 
5460  const GEOSCoordSequence* lineSequence = GEOSGeom_getCoordSeq( currentIntersectGeom );
5461  unsigned int sequenceSize = 0;
5462  double x, y;
5463  if ( GEOSCoordSeq_getSize( lineSequence, &sequenceSize ) != 0 )
5464  {
5465  for ( unsigned int i = 0; i < sequenceSize; ++i )
5466  {
5467  if ( GEOSCoordSeq_getX( lineSequence, i, &x ) != 0 )
5468  {
5469  if ( GEOSCoordSeq_getY( lineSequence, i, &y ) != 0 )
5470  {
5471  testPoints.push_back( QgsPoint( x, y ) );
5472  }
5473  }
5474  }
5475  }
5476  }
5477  GEOSGeom_destroy( intersectionGeom );
5478  return 0;
5479 }
5480 
5481 GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
5482 {
5483  if ( !splitLine || !geom )
5484  {
5485  return 0;
5486  }
5487 
5488  GEOSGeometry *geometryBoundary = 0;
5489  if ( GEOSGeomTypeId( geom ) == GEOS_POLYGON || GEOSGeomTypeId( geom ) == GEOS_MULTIPOLYGON )
5490  {
5491  geometryBoundary = GEOSBoundary( geom );
5492  }
5493  else
5494  {
5495  geometryBoundary = GEOSGeom_clone( geom );
5496  }
5497 
5498  GEOSGeometry *splitLineClone = GEOSGeom_clone( splitLine );
5499  GEOSGeometry *unionGeometry = GEOSUnion( splitLineClone, geometryBoundary );
5500  GEOSGeom_destroy( splitLineClone );
5501 
5502  GEOSGeom_destroy( geometryBoundary );
5503  return unionGeometry;
5504 }
5505 
5506 int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 )
5507 {
5508  if ( !line1 || !line2 )
5509  {
5510  return -1;
5511  }
5512 
5513 
5514  double bufferDistance = 0.00001;
5515  if ( geomInDegrees( line2 ) ) //use more accurate tolerance for degrees
5516  {
5517  bufferDistance = 0.00000001;
5518  }
5519  GEOSGeometry* bufferGeom = GEOSBuffer( line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS );
5520  if ( !bufferGeom )
5521  {
5522  return -2;
5523  }
5524 
5525  GEOSGeometry* intersectionGeom = GEOSIntersection( bufferGeom, line1 );
5526 
5527  //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
5528  double intersectGeomLength;
5529  double line1Length;
5530 
5531  GEOSLength( intersectionGeom, &intersectGeomLength );
5532  GEOSLength( line1, &line1Length );
5533 
5534  GEOSGeom_destroy( bufferGeom );
5535  GEOSGeom_destroy( intersectionGeom );
5536 
5537  double intersectRatio = line1Length / intersectGeomLength;
5538  if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
5539  {
5540  return 1;
5541  }
5542  return 0;
5543 }
5544 
5545 int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line )
5546 {
5547  if ( !point || !line )
5548  {
5549  return -1;
5550  }
5551 
5552  double bufferDistance = 0.000001;
5553  if ( geomInDegrees( line ) )
5554  {
5555  bufferDistance = 0.00000001;
5556  }
5557  GEOSGeometry* lineBuffer = GEOSBuffer( line, bufferDistance, 8 );
5558  if ( !lineBuffer )
5559  {
5560  return -2;
5561  }
5562 
5563  bool contained = false;
5564  if ( GEOSContains( lineBuffer, point ) == 1 )
5565  {
5566  contained = true;
5567  }
5568 
5569  GEOSGeom_destroy( lineBuffer );
5570  return contained;
5571 }
5572 
5573 bool QgsGeometry::geomInDegrees( const GEOSGeometry* geom )
5574 {
5575  GEOSGeometry* bbox = GEOSEnvelope( geom );
5576  if ( !bbox )
5577  {
5578  return false;
5579  }
5580 
5581  const GEOSGeometry* bBoxRing = GEOSGetExteriorRing( bbox );
5582  if ( !bBoxRing )
5583  {
5584  return false;
5585  }
5586  const GEOSCoordSequence* bBoxCoordSeq = GEOSGeom_getCoordSeq( bBoxRing );
5587 
5588  if ( !bBoxCoordSeq )
5589  {
5590  return false;
5591  }
5592 
5593  unsigned int nCoords = 0;
5594  if ( !GEOSCoordSeq_getSize( bBoxCoordSeq, &nCoords ) )
5595  {
5596  return false;
5597  }
5598 
5599  double x, y;
5600  for ( unsigned int i = 0; i < ( nCoords - 1 ); ++i )
5601  {
5602  GEOSCoordSeq_getX( bBoxCoordSeq, i, &x );
5603  if ( x > 180 || x < -180 )
5604  {
5605  return false;
5606  }
5607  GEOSCoordSeq_getY( bBoxCoordSeq, i, &y );
5608  if ( y > 90 || y < -90 )
5609  {
5610  return false;
5611  }
5612  }
5613 
5614  return true;
5615 }
5616 
5617 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const
5618 {
5619  if ( !g )
5620  {
5621  return 0;
5622  }
5623  int geometryType = GEOSGeomTypeId( g );
5624  if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
5625  || geometryType == GEOS_POLYGON )
5626  {
5627  return 1;
5628  }
5629 
5630  //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
5631  return GEOSGetNumGeometries( g );
5632 }
5633 
5634 int QgsGeometry::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry*>& splitResult )
5635 {
5636  if ( !mGeos || mDirtyGeos )
5637  {
5638  if ( !exportWkbToGeos() )
5639  return 1;
5640  }
5641 
5642  //convert mGeos to geometry collection
5643  int type = GEOSGeomTypeId( mGeos );
5644  if ( type != GEOS_GEOMETRYCOLLECTION &&
5645  type != GEOS_MULTILINESTRING &&
5646  type != GEOS_MULTIPOLYGON &&
5647  type != GEOS_MULTIPOINT )
5648  {
5649  return 0;
5650  }
5651 
5652  QVector<GEOSGeometry*> copyList = splitResult;
5653  splitResult.clear();
5654 
5655  //collect all the geometries that belong to the initial multifeature
5656  QVector<GEOSGeometry*> unionGeom;
5657 
5658  for ( int i = 0; i < copyList.size(); ++i )
5659  {
5660  //is this geometry a part of the original multitype?
5661  bool isPart = false;
5662  for ( int j = 0; j < GEOSGetNumGeometries( mGeos ); j++ )
5663  {
5664  if ( GEOSEquals( copyList[i], GEOSGetGeometryN( mGeos, j ) ) )
5665  {
5666  isPart = true;
5667  break;
5668  }
5669  }
5670 
5671  if ( isPart )
5672  {
5673  unionGeom << copyList[i];
5674  }
5675  else
5676  {
5677  QVector<GEOSGeometry*> geomVector;
5678  geomVector << copyList[i];
5679 
5680  if ( type == GEOS_MULTILINESTRING )
5681  {
5682  splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector );
5683  }
5684  else if ( type == GEOS_MULTIPOLYGON )
5685  {
5686  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
5687  }
5688  else
5689  {
5690  GEOSGeom_destroy( copyList[i] );
5691  }
5692  }
5693  }
5694 
5695  //make multifeature out of unionGeom
5696  if ( unionGeom.size() > 0 )
5697  {
5698  if ( type == GEOS_MULTILINESTRING )
5699  {
5700  splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom );
5701  }
5702  else if ( type == GEOS_MULTIPOLYGON )
5703  {
5704  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom );
5705  }
5706  }
5707  else
5708  {
5709  unionGeom.clear();
5710  }
5711 
5712  return 0;
5713 }
5714 
5715 QgsPoint QgsGeometry::asPoint( unsigned char*& ptr, bool hasZValue )
5716 {
5717  ptr += 5;
5718  double* x = ( double * )( ptr );
5719  double* y = ( double * )( ptr + sizeof( double ) );
5720  ptr += 2 * sizeof( double );
5721 
5722  if ( hasZValue )
5723  ptr += sizeof( double );
5724 
5725  return QgsPoint( *x, *y );
5726 }
5727 
5728 
5729 QgsPolyline QgsGeometry::asPolyline( unsigned char*& ptr, bool hasZValue )
5730 {
5731  double x, y;
5732  ptr += 5;
5733  unsigned int nPoints = *(( int* )ptr );
5734  ptr += 4;
5735 
5736  QgsPolyline line( nPoints );
5737 
5738  // Extract the points from the WKB format into the x and y vectors.
5739  for ( uint i = 0; i < nPoints; ++i )
5740  {
5741  x = *(( double * ) ptr );
5742  y = *(( double * )( ptr + sizeof( double ) ) );
5743 
5744  ptr += 2 * sizeof( double );
5745 
5746  line[i] = QgsPoint( x, y );
5747 
5748  if ( hasZValue ) // ignore Z value
5749  ptr += sizeof( double );
5750  }
5751 
5752  return line;
5753 }
5754 
5755 
5756 QgsPolygon QgsGeometry::asPolygon( unsigned char*& ptr, bool hasZValue )
5757 {
5758  double x, y;
5759 
5760  ptr += 5;
5761 
5762  // get number of rings in the polygon
5763  unsigned int numRings = *(( int* )ptr );
5764  ptr += 4;
5765 
5766  if ( numRings == 0 ) // sanity check for zero rings in polygon
5767  return QgsPolygon();
5768 
5769  QgsPolygon rings( numRings );
5770 
5771  for ( uint idx = 0; idx < numRings; idx++ )
5772  {
5773  uint nPoints = *(( int* )ptr );
5774  ptr += 4;
5775 
5776  QgsPolyline ring( nPoints );
5777 
5778  for ( uint jdx = 0; jdx < nPoints; jdx++ )
5779  {
5780  x = *(( double * ) ptr );
5781  y = *(( double * )( ptr + sizeof( double ) ) );
5782 
5783  ptr += 2 * sizeof( double );
5784 
5785  if ( hasZValue )
5786  ptr += sizeof( double );
5787 
5788  ring[jdx] = QgsPoint( x, y );
5789  }
5790 
5791  rings[idx] = ring;
5792  }
5793 
5794  return rings;
5795 }
5796 
5797 
5799 {
5800  QGis::WkbType type = wkbType();
5801  if ( type != QGis::WKBPoint && type != QGis::WKBPoint25D )
5802  return QgsPoint( 0, 0 );
5803 
5804  unsigned char* ptr = mGeometry;
5805  return asPoint( ptr, type == QGis::WKBPoint25D );
5806 }
5807 
5809 {
5810  QGis::WkbType type = wkbType();
5811  if ( type != QGis::WKBLineString && type != QGis::WKBLineString25D )
5812  return QgsPolyline();
5813 
5814  unsigned char *ptr = mGeometry;
5815  return asPolyline( ptr, type == QGis::WKBLineString25D );
5816 }
5817 
5819 {
5820  QGis::WkbType type = wkbType();
5821  if ( type != QGis::WKBPolygon && type != QGis::WKBPolygon25D )
5822  return QgsPolygon();
5823 
5824  unsigned char *ptr = mGeometry;
5825  return asPolygon( ptr, type == QGis::WKBPolygon25D );
5826 }
5827 
5829 {
5830  QGis::WkbType type = wkbType();
5831  if ( type != QGis::WKBMultiPoint && type != QGis::WKBMultiPoint25D )
5832  return QgsMultiPoint();
5833 
5834  bool hasZValue = ( type == QGis::WKBMultiPoint25D );
5835 
5836  unsigned char* ptr = mGeometry + 5;
5837  unsigned int nPoints = *(( int* )ptr );
5838  ptr += 4;
5839 
5840  QgsMultiPoint points( nPoints );
5841  for ( uint i = 0; i < nPoints; i++ )
5842  {
5843  points[i] = asPoint( ptr, hasZValue );
5844  }
5845 
5846  return points;
5847 }
5848 
5850 {
5851  QGis::WkbType type = wkbType();
5852  if ( type != QGis::WKBMultiLineString && type != QGis::WKBMultiLineString25D )
5853  return QgsMultiPolyline();
5854 
5855  bool hasZValue = ( type == QGis::WKBMultiLineString25D );
5856 
5857  unsigned char* ptr = mGeometry + 5;
5858  unsigned int numLineStrings = *(( int* )ptr );
5859  ptr += 4;
5860 
5861  QgsMultiPolyline lines( numLineStrings );
5862 
5863  for ( uint i = 0; i < numLineStrings; i++ )
5864  {
5865  lines[i] = asPolyline( ptr, hasZValue );
5866  }
5867 
5868  return lines;
5869 }
5870 
5872 {
5873  QGis::WkbType type = wkbType();
5874  if ( type != QGis::WKBMultiPolygon && type != QGis::WKBMultiPolygon25D )
5875  return QgsMultiPolygon();
5876 
5877  bool hasZValue = ( type == QGis::WKBMultiPolygon25D );
5878 
5879  unsigned char* ptr = mGeometry + 5;
5880  unsigned int numPolygons = *(( int* )ptr );
5881  ptr += 4;
5882 
5883  QgsMultiPolygon polygons( numPolygons );
5884 
5885  for ( uint i = 0; i < numPolygons; i++ )
5886  {
5887  polygons[i] = asPolygon( ptr, hasZValue );
5888  }
5889 
5890  return polygons;
5891 }
5892 
5894 {
5895  if ( !mGeos )
5896  {
5897  exportWkbToGeos();
5898  }
5899 
5900  double area;
5901 
5902  try
5903  {
5904  if ( GEOSArea( mGeos, &area ) == 0 )
5905  return -1.0;
5906  }
5907  CATCH_GEOS( -1.0 )
5908 
5909  return area;
5910 }
5911 
5913 {
5914  if ( !mGeos )
5915  {
5916  exportWkbToGeos();
5917  }
5918 
5919  double length;
5920 
5921  try
5922  {
5923  if ( GEOSLength( mGeos, &length ) == 0 )
5924  return -1.0;
5925  }
5926  CATCH_GEOS( -1.0 )
5927 
5928  return length;
5929 }
5931 {
5932  if ( !mGeos )
5933  {
5934  exportWkbToGeos();
5935  }
5936 
5937  if ( !geom.mGeos )
5938  {
5939  geom.exportWkbToGeos();
5940  }
5941 
5942  if ( !mGeos || !geom.mGeos )
5943  return -1.0;
5944 
5945  double dist = -1.0;
5946 
5947  try
5948  {
5949  GEOSDistance( mGeos, geom.mGeos, &dist );
5950  }
5951  CATCH_GEOS( -1.0 )
5952 
5953  return dist;
5954 }
5955 
5956 
5957 QgsGeometry* QgsGeometry::buffer( double distance, int segments )
5958 {
5959  if ( !mGeos )
5960  {
5961  exportWkbToGeos();
5962  }
5963  if ( !mGeos )
5964  {
5965  return 0;
5966  }
5967 
5968  try
5969  {
5970  return fromGeosGeom( GEOSBuffer( mGeos, distance, segments ) );
5971  }
5972  CATCH_GEOS( 0 )
5973 }
5974 
5976 {
5977  if ( !mGeos )
5978  {
5979  exportWkbToGeos();
5980  }
5981  if ( !mGeos )
5982  {
5983  return 0;
5984  }
5985  try
5986  {
5987  return fromGeosGeom( GEOSSimplify( mGeos, tolerance ) );
5988  }
5989  CATCH_GEOS( 0 )
5990 }
5991 
5993 {
5994  if ( !mGeos )
5995  {
5996  exportWkbToGeos();
5997  }
5998  if ( !mGeos )
5999  {
6000  return 0;
6001  }
6002  try
6003  {
6004  return fromGeosGeom( GEOSGetCentroid( mGeos ) );
6005  }
6006  CATCH_GEOS( 0 )
6007 }
6008 
6010 {
6011  if ( !mGeos )
6012  {
6013  exportWkbToGeos();
6014  }
6015  if ( !mGeos )
6016  {
6017  return 0;
6018  }
6019 
6020  try
6021  {
6022  return fromGeosGeom( GEOSConvexHull( mGeos ) );
6023  }
6024  CATCH_GEOS( 0 )
6025 }
6026 
6028 {
6029  if ( !geometry )
6030  {
6031  return NULL;
6032  }
6033  if ( !mGeos )
6034  {
6035  exportWkbToGeos();
6036  }
6037  if ( !geometry->mGeos )
6038  {
6039  geometry->exportWkbToGeos();
6040  }
6041  if ( !mGeos || !geometry->mGeos )
6042  {
6043  return 0;
6044  }
6045 
6046  try
6047  {
6048  return fromGeosGeom( GEOSIntersection( mGeos, geometry->mGeos ) );
6049  }
6050  CATCH_GEOS( 0 )
6051 }
6052 
6054 {
6055  if ( !geometry )
6056  {
6057  return NULL;
6058  }
6059  if ( !mGeos )
6060  {
6061  exportWkbToGeos();
6062  }
6063  if ( !geometry->mGeos )
6064  {
6065  geometry->exportWkbToGeos();
6066  }
6067  if ( !mGeos || !geometry->mGeos )
6068  {
6069  return 0;
6070  }
6071 
6072  try
6073  {
6074  GEOSGeometry* unionGeom = GEOSUnion( mGeos, geometry->mGeos );
6075  QGis::WkbType thisGeomType = wkbType();
6076  QGis::WkbType otherGeomType = geometry->wkbType();
6077  if (( thisGeomType == QGis::WKBLineString || thisGeomType == QGis::WKBLineString25D ) \
6078  && ( otherGeomType == QGis::WKBLineString || otherGeomType == QGis::WKBLineString25D ) )
6079  {
6080  GEOSGeometry* mergedGeom = GEOSLineMerge( unionGeom );
6081  if ( mergedGeom )
6082  {
6083  GEOSGeom_destroy( unionGeom );
6084  unionGeom = mergedGeom;
6085  }
6086  }
6087  return fromGeosGeom( unionGeom );
6088  }
6089  CATCH_GEOS( new QgsGeometry( *this ) ) //return this geometry if union not possible
6090 }
6091 
6093 {
6094  if ( !geometry )
6095  {
6096  return NULL;
6097  }
6098  if ( !mGeos )
6099  {
6100  exportWkbToGeos();
6101  }
6102  if ( !geometry->mGeos )
6103  {
6104  geometry->exportWkbToGeos();
6105  }
6106  if ( !mGeos || !geometry->mGeos )
6107  {
6108  return 0;
6109  }
6110 
6111  try
6112  {
6113  return fromGeosGeom( GEOSDifference( mGeos, geometry->mGeos ) );
6114  }
6115  CATCH_GEOS( 0 )
6116 }
6117 
6119 {
6120  if ( !geometry )
6121  {
6122  return NULL;
6123  }
6124  if ( !mGeos )
6125  {
6126  exportWkbToGeos();
6127  }
6128  if ( !geometry->mGeos )
6129  {
6130  geometry->exportWkbToGeos();
6131  }
6132  if ( !mGeos || !geometry->mGeos )
6133  {
6134  return 0;
6135  }
6136 
6137  try
6138  {
6139  return fromGeosGeom( GEOSSymDifference( mGeos, geometry->mGeos ) );
6140  }
6141  CATCH_GEOS( 0 )
6142 }
6143 
6145 {
6146  if ( !mGeos )
6147  {
6148  exportWkbToGeos();
6149  if ( !mGeos )
6150  return QList<QgsGeometry*>();
6151  }
6152 
6153  int type = GEOSGeomTypeId( mGeos );
6154  QgsDebugMsg( "geom type: " + QString::number( type ) );
6155 
6156  QList<QgsGeometry*> geomCollection;
6157 
6158  if ( type != GEOS_MULTIPOINT &&
6159  type != GEOS_MULTILINESTRING &&
6160  type != GEOS_MULTIPOLYGON &&
6161  type != GEOS_GEOMETRYCOLLECTION )
6162  {
6163  // we have a single-part geometry - put there a copy of this one
6164  geomCollection.append( new QgsGeometry( *this ) );
6165  return geomCollection;
6166  }
6167 
6168  int count = GEOSGetNumGeometries( mGeos );
6169  QgsDebugMsg( "geom count: " + QString::number( count ) );
6170 
6171  for ( int i = 0; i < count; ++i )
6172  {
6173  const GEOSGeometry * geometry = GEOSGetGeometryN( mGeos, i );
6174  geomCollection.append( fromGeosGeom( GEOSGeom_clone( geometry ) ) );
6175  }
6176 
6177  return geomCollection;
6178 }
6179 
6180 
6181 bool QgsGeometry::deleteRing( int ringNum, int partNum )
6182 {
6183  if ( ringNum <= 0 || partNum < 0 )
6184  return false;
6185 
6186  switch ( wkbType() )
6187  {
6188  case QGis::WKBPolygon25D:
6189  case QGis::WKBPolygon:
6190  {
6191  if ( partNum != 0 )
6192  return false;
6193 
6194  QgsPolygon polygon = asPolygon();
6195  if ( ringNum >= polygon.count() )
6196  return false;
6197 
6198  polygon.remove( ringNum );
6199 
6200  QgsGeometry* g2 = QgsGeometry::fromPolygon( polygon );
6201  *this = *g2;
6202  delete g2;
6203  return true;
6204  }
6205 
6207  case QGis::WKBMultiPolygon:
6208  {
6209  QgsMultiPolygon mpolygon = asMultiPolygon();
6210 
6211  if ( partNum >= mpolygon.count() )
6212  return false;
6213 
6214  if ( ringNum >= mpolygon[partNum].count() )
6215  return false;
6216 
6217  mpolygon[partNum].remove( ringNum );
6218 
6219  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
6220  *this = *g2;
6221  delete g2;
6222  return true;
6223  }
6224 
6225  default:
6226  return false; // only makes sense with polygons and multipolygons
6227  }
6228 }
6229 
6230 
6231 bool QgsGeometry::deletePart( int partNum )
6232 {
6233  if ( partNum < 0 )
6234  return false;
6235 
6236  switch ( wkbType() )
6237  {
6239  case QGis::WKBMultiPoint:
6240  {
6241  QgsMultiPoint mpoint = asMultiPoint();
6242 
6243  if ( partNum >= mpoint.size() || mpoint.size() == 1 )
6244  return false;
6245 
6246  mpoint.remove( partNum );
6247 
6248  QgsGeometry* g2 = QgsGeometry::fromMultiPoint( mpoint );
6249  *this = *g2;
6250  delete g2;
6251  break;
6252  }
6253 
6256  {
6258 
6259  if ( partNum >= mline.size() || mline.size() == 1 )
6260  return false;
6261 
6262  mline.remove( partNum );
6263 
6265  *this = *g2;
6266  delete g2;
6267  break;
6268  }
6269 
6271  case QGis::WKBMultiPolygon:
6272  {
6273  QgsMultiPolygon mpolygon = asMultiPolygon();
6274 
6275  if ( partNum >= mpolygon.size() || mpolygon.size() == 1 )
6276  return false;
6277 
6278  mpolygon.remove( partNum );
6279 
6280  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
6281  *this = *g2;
6282  delete g2;
6283  break;
6284  }
6285 
6286  default:
6287  // single part geometries are ignored
6288  return false;
6289  }
6290 
6291  return true;
6292 }
6293 
6295 {
6296  int returnValue = 0;
6297 
6298  //check if g has polygon type
6299  if ( type() != QGis::Polygon )
6300  {
6301  return 1;
6302  }
6303 
6304  QGis::WkbType geomTypeBeforeModification = wkbType();
6305 
6306  //read avoid intersections list from project properties
6307  bool listReadOk;
6308  QStringList avoidIntersectionsList = QgsProject::instance()->readListEntry( "Digitizing", "/AvoidIntersectionsList", &listReadOk );
6309  if ( !listReadOk )
6310  {
6311  return true; //no intersections stored in project does not mean error
6312  }
6313 
6314  //go through list, convert each layer to vector layer and call QgsVectorLayer::removePolygonIntersections for each
6315  QgsVectorLayer* currentLayer = 0;
6316  QStringList::const_iterator aIt = avoidIntersectionsList.constBegin();
6317  for ( ; aIt != avoidIntersectionsList.constEnd(); ++aIt )
6318  {
6319  currentLayer = dynamic_cast<QgsVectorLayer*>( QgsMapLayerRegistry::instance()->mapLayer( *aIt ) );
6320  if ( currentLayer )
6321  {
6322  if ( currentLayer->removePolygonIntersections( this ) != 0 )
6323  {
6324  returnValue = 3;
6325  }
6326  }
6327  }
6328 
6329  //make sure the geometry still has the same type (e.g. no change from polygon to multipolygon)
6330  if ( wkbType() != geomTypeBeforeModification )
6331  {
6332  return 2;
6333  }
6334 
6335  return returnValue;
6336 }
6337 
6338 //
6339 // distance of point q from line through p in direction v
6340 // return >0 => q lies left of the line
6341 // <0 => q lies right of the line
6342 //
6343 static double distLine2Point( QgsPoint p, QgsVector v, QgsPoint q )
6344 {
6345  if ( v.length() == 0 )
6346  {
6347  throw QgsException( QObject::tr( "invalid line" ) );
6348  }
6349 
6350  return ( v.x()*( q.y() - p.y() ) - v.y()*( q.x() - p.x() ) ) / v.length();
6351 }
6352 
6354 {
6355  double d = v.y() * w.x() - v.x() * w.y();
6356 
6357  if ( d == 0 )
6358  return false;
6359 
6360  double dx = q.x() - p.x();
6361  double dy = q.y() - p.y();
6362  double k = ( dy * w.x() - dx * w.y() ) / d;
6363 
6364  s = p + v * k;
6365 
6366  return true;
6367 }
6368 
6369 bool pointInRing( const QgsPolyline &ring, const QgsPoint &p )
6370 {
6371  bool inside = false;
6372  int j = ring.size() - 1;
6373 
6374  for ( int i = 0; i < ring.size(); i++ )
6375  {
6376  if ( ring[i].x() == p.x() && ring[i].y() == p.y() )
6377  return true;
6378 
6379  if (( ring[i].y() < p.y() && ring[j].y() >= p.y() ) ||
6380  ( ring[j].y() < p.y() && ring[i].y() >= p.y() ) )
6381  {
6382  if ( ring[i].x() + ( p.y() - ring[i].y() ) / ( ring[j].y() - ring[i].y() )*( ring[j].x() - ring[i].x() ) <= p.x() )
6383  inside = !inside;
6384  }
6385 
6386  j = i;
6387  }
6388 
6389  return inside;
6390 }
6391 
6392 static bool ringInRing( const QgsPolyline &inside, const QgsPolyline &outside )
6393 {
6394  for ( int i = 0; i < inside.size(); i++ )
6395  {
6396  if ( !pointInRing( outside, inside[i] ) )
6397  return false;
6398  }
6399 
6400  return true;
6401 }
6402 
6403 void QgsGeometry::checkRingIntersections( QList<Error> &errors,
6404  int p0, int i0, const QgsPolyline &ring0,
6405  int p1, int i1, const QgsPolyline &ring1 )
6406 {
6407  for ( int i = 0; i < ring0.size() - 1; i++ )
6408  {
6409  QgsVector v = ring0[i+1] - ring0[i];
6410 
6411  for ( int j = 0; j < ring1.size() - 1; j++ )
6412  {
6413  QgsVector w = ring1[j+1] - ring1[j];
6414 
6415  QgsPoint s;
6416  if ( intersectLines( ring0[i], v, ring1[j], w, s ) )
6417  {
6418  double d = -distLine2Point( ring0[i], v.perpVector(), s );
6419 
6420  if ( d >= 0 && d <= v.length() )
6421  {
6422  d = -distLine2Point( ring1[j], w.perpVector(), s );
6423  if ( d >= 0 && d <= w.length() )
6424  {
6425  QString msg = QObject::tr( "segment %1 of ring %2 of polygon %3 intersects segment %4 of ring %5 of polygon %6 at %7" )
6426  .arg( i0 ).arg( i ).arg( p0 )
6427  .arg( i1 ).arg( j ).arg( p1 )
6428  .arg( s.toString() );
6429  QgsDebugMsg( msg );
6430  errors << Error( msg, s );
6431  if ( errors.size() > 100 )
6432  {
6433  QString msg = QObject::tr( "stopping validation after more than 100 errors" );
6434  QgsDebugMsg( msg );
6435  errors << Error( msg );
6436  return;
6437  }
6438  }
6439  }
6440  }
6441  }
6442  }
6443 }
6444 
6445 void QgsGeometry::validatePolyline( QList<Error> &errors, int i, QgsPolyline line, bool ring )
6446 {
6447  if ( ring )
6448  {
6449  if ( line.size() < 3 )
6450  {
6451  QString msg = QObject::tr( "ring %1 with less than three points" ).arg( i );
6452  QgsDebugMsg( msg );
6453  errors << Error( msg );
6454  return;
6455  }
6456 
6457  if ( line[0] != line[ line.size()-1 ] )
6458  {
6459  QString msg = QObject::tr( "ring %1 not closed" ).arg( i );
6460  QgsDebugMsg( msg );
6461  errors << Error( msg );
6462  return;
6463  }
6464  }
6465  else if ( line.size() < 2 )
6466  {
6467  QString msg = QObject::tr( "line %1 with less than two points" ).arg( i );
6468  QgsDebugMsg( msg );
6469  errors << Error( msg );
6470  return;
6471  }
6472 
6473  int j = 0;
6474  while ( j < line.size() - 1 )
6475  {
6476  int n = 0;
6477  while ( j < line.size() - 1 && line[j] == line[j+1] )
6478  {
6479  line.remove( j );
6480  n++;
6481  }
6482 
6483  if ( n > 0 )
6484  {
6485  QString msg = QObject::tr( "line %1 contains %n duplicate node(s) at %2", "number of duplicate nodes", n ).arg( i ).arg( j );
6486  QgsDebugMsg( msg );
6487  errors << Error( msg, line[j] );
6488  }
6489 
6490  j++;
6491  }
6492 
6493  for ( j = 0; j < line.size() - 3; j++ )
6494  {
6495  QgsVector v = line[j+1] - line[j];
6496  double vl = v.length();
6497 
6498  int n = ( j == 0 && ring ) ? line.size() - 2 : line.size() - 1;
6499 
6500  for ( int k = j + 2; k < n; k++ )
6501  {
6502  QgsVector w = line[k+1] - line[k];
6503 
6504  QgsPoint s;
6505  if ( !intersectLines( line[j], v, line[k], w, s ) )
6506  continue;
6507 
6508  double d = -distLine2Point( line[j], v.perpVector(), s );
6509  if ( d < 0 || d > vl )
6510  continue;
6511 
6512  d = -distLine2Point( line[k], w.perpVector(), s );
6513  if ( d < 0 || d > w.length() )
6514  continue;
6515 
6516  QString msg = QObject::tr( "segments %1 and %2 of line %3 intersect at %4" ).arg( j ).arg( k ).arg( i ).arg( s.toString() );
6517  QgsDebugMsg( msg );
6518  errors << Error( msg, s );
6519  if ( errors.size() > 100 )
6520  {
6521  QString msg = QObject::tr( "stopping validation after more than 100 errors" );
6522  QgsDebugMsg( msg );
6523  errors << Error( msg );
6524  return;
6525  }
6526  }
6527  }
6528 }
6529 
6530 void QgsGeometry::validatePolygon( QList<Error> &errors, int idx, const QgsPolygon &polygon )
6531 {
6532  // check if holes are inside polygon
6533  for ( int i = 1; i < polygon.size(); i++ )
6534  {
6535  if ( !ringInRing( polygon[i], polygon[0] ) )
6536  {
6537  QString msg = QObject::tr( "ring %1 of polygon %2 not in exterior ring" ).arg( i ).arg( idx );
6538  QgsDebugMsg( msg );
6539  errors << Error( msg );
6540  }
6541  }
6542 
6543  // check holes for intersections
6544  for ( int i = 1; i < polygon.size(); i++ )
6545  {
6546  for ( int j = i + 1; j < polygon.size(); j++ )
6547  {
6548  checkRingIntersections( errors, idx, i, polygon[i], idx, j, polygon[j] );
6549  }
6550  }
6551 
6552  // check if rings are self-intersecting
6553  for ( int i = 0; i < polygon.size(); i++ )
6554  {
6555  validatePolyline( errors, i, polygon[i], true );
6556  }
6557 }
6558 
6559 void QgsGeometry::validateGeometry( QList<Error> &errors )
6560 {
6561  errors.clear();
6562 
6563  switch ( wkbType() )
6564  {
6565  case QGis::WKBPoint:
6566  case QGis::WKBPoint25D:
6567  case QGis::WKBMultiPoint:
6569  break;
6570 
6571  case QGis::WKBLineString:
6573  validatePolyline( errors, 0, asPolyline() );
6574  break;
6575 
6578  {
6580  for ( int i = 0; i < mp.size(); i++ )
6581  validatePolyline( errors, i, mp[i] );
6582  }
6583  break;
6584 
6585  case QGis::WKBPolygon:
6586  case QGis::WKBPolygon25D:
6587  {
6588  validatePolygon( errors, 0, asPolygon() );
6589  }
6590  break;
6591 
6592  case QGis::WKBMultiPolygon:
6594  {
6596  for ( int i = 0; i < mp.size(); i++ )
6597  {
6598  validatePolygon( errors, i, mp[i] );
6599  }
6600 
6601  for ( int i = 0; i < mp.size(); i++ )
6602  {
6603  for ( int j = i + 1; j < mp.size(); j++ )
6604  {
6605  if ( ringInRing( mp[i][0], mp[j][0] ) )
6606  {
6607  errors << Error( QObject::tr( "polygon %1 inside polygon %2" ).arg( i ).arg( j ) );
6608  }
6609  else if ( ringInRing( mp[j][0], mp[i][0] ) )
6610  {
6611  errors << Error( QObject::tr( "polygon %1 inside polygon %2" ).arg( j ).arg( i ) );
6612  }
6613  else
6614  {
6615  checkRingIntersections( errors, i, 0, mp[i][0], j, 0, mp[j][0] );
6616  }
6617  }
6618  }
6619  }
6620  break;
6621 
6622  case QGis::WKBNoGeometry:
6623  case QGis::WKBUnknown:
6624  QgsDebugMsg( QObject::tr( "Unknown geometry type" ) );
6625  errors << Error( QObject::tr( "Unknown geometry type" ) );
6626  break;
6627  }
6628 }
6629 
6631 {
6632  try
6633  {
6634  GEOSGeometry *g = asGeos();
6635 
6636  if ( !g )
6637  return false;
6638 
6639  return GEOSisValid( g );
6640  }
6641  catch ( GEOSException &e )
6642  {
6643  // looks like geometry is fubar
6644  Q_UNUSED( e );
6645  QgsDebugMsg( QString( "GEOS exception caught: %1" ).arg( e.what() ) );
6646  return false;
6647  }
6648 }
6649 
6651 {
6652  return geosRelOp( GEOSEquals, this, &g );
6653 }
6654 
6656 {
6657  try
6658  {
6659  GEOSGeometry *g = asGeos();
6660 
6661  if ( !g )
6662  return false;
6663 
6664  return GEOSisEmpty( g );
6665  }
6666  catch ( GEOSException &e )
6667  {
6668  // looks like geometry is fubar
6669  Q_UNUSED( e );
6670  QgsDebugMsg( QString( "GEOS exception caught: %1" ).arg( e.what() ) );
6671  return false;
6672  }
6673 }