00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <limits>
00018 #include <cstdarg>
00019 #include <cstdio>
00020 #include <cmath>
00021
00022 #include "qgis.h"
00023 #include "qgsgeometry.h"
00024 #include "qgsapplication.h"
00025 #include "qgslogger.h"
00026 #include "qgspoint.h"
00027 #include "qgsrectangle.h"
00028 #include "qgslogger.h"
00029
00030 #include "qgsmaplayerregistry.h"
00031 #include "qgsvectorlayer.h"
00032 #include "qgsproject.h"
00033
00034 #define DEFAULT_QUADRANT_SEGMENTS 8
00035
00036 #define CATCH_GEOS(r) \
00037 catch (GEOSException &e) \
00038 { \
00039 Q_UNUSED(e); \
00040 QgsDebugMsg("GEOS: " + QString( e.what() ) ); \
00041 return r; \
00042 }
00043
00044 class GEOSException
00045 {
00046 public:
00047 GEOSException( QString theMsg )
00048 {
00049 if ( theMsg == "Unknown exception thrown" && lastMsg.isNull() )
00050 {
00051 msg = theMsg;
00052 }
00053 else
00054 {
00055 msg = theMsg;
00056 lastMsg = msg;
00057 }
00058 }
00059
00060
00061 GEOSException( const GEOSException &rhs )
00062 {
00063 *this = rhs;
00064 }
00065
00066 ~GEOSException()
00067 {
00068 if ( lastMsg == msg )
00069 lastMsg = QString::null;
00070 }
00071
00072 QString what()
00073 {
00074 return msg;
00075 }
00076
00077 private:
00078 QString msg;
00079 static QString lastMsg;
00080 };
00081
00082 QString GEOSException::lastMsg;
00083
00084 static void throwGEOSException( const char *fmt, ... )
00085 {
00086 va_list ap;
00087 char buffer[1024];
00088
00089 va_start( ap, fmt );
00090 vsnprintf( buffer, sizeof buffer, fmt, ap );
00091 va_end( ap );
00092
00093 QgsDebugMsg( QString( "GEOS exception encountered: %1" ).arg( buffer ) );
00094
00095 throw GEOSException( QString::fromUtf8( buffer ) );
00096 }
00097
00098 static void printGEOSNotice( const char *fmt, ... )
00099 {
00100 #if defined(QGISDEBUG)
00101 va_list ap;
00102 char buffer[1024];
00103
00104 va_start( ap, fmt );
00105 vsnprintf( buffer, sizeof buffer, fmt, ap );
00106 va_end( ap );
00107
00108 QgsDebugMsg( QString( "GEOS notice: %1" ).arg( QString::fromUtf8( buffer ) ) );
00109 #endif
00110 }
00111
00112 class GEOSInit
00113 {
00114 public:
00115 GEOSInit()
00116 {
00117 initGEOS( printGEOSNotice, throwGEOSException );
00118 }
00119
00120 ~GEOSInit()
00121 {
00122 finishGEOS();
00123 }
00124 };
00125
00126 static GEOSInit geosinit;
00127
00128
00129 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR<3)
00130 #define GEOSGeom_getCoordSeq(g) GEOSGeom_getCoordSeq( (GEOSGeometry *) g )
00131 #define GEOSGetExteriorRing(g) GEOSGetExteriorRing( (GEOSGeometry *)g )
00132 #define GEOSGetNumInteriorRings(g) GEOSGetNumInteriorRings( (GEOSGeometry *)g )
00133 #define GEOSGetInteriorRingN(g,i) GEOSGetInteriorRingN( (GEOSGeometry *)g, i )
00134 #define GEOSDisjoint(g0,g1) GEOSDisjoint( (GEOSGeometry *)g0, (GEOSGeometry*)g1 )
00135 #define GEOSIntersection(g0,g1) GEOSIntersection( (GEOSGeometry*) g0, (GEOSGeometry*)g1 )
00136 #define GEOSBuffer(g, d, s) GEOSBuffer( (GEOSGeometry*) g, d, s )
00137 #define GEOSArea(g, a) GEOSArea( (GEOSGeometry*) g, a )
00138 #define GEOSSimplify(g, t) GEOSSimplify( (GEOSGeometry*) g, t )
00139 #define GEOSGetCentroid(g) GEOSGetCentroid( (GEOSGeometry*) g )
00140
00141 #define GEOSCoordSeq_getSize(cs,n) GEOSCoordSeq_getSize( (GEOSCoordSequence *) cs, n )
00142 #define GEOSCoordSeq_getX(cs,i,x) GEOSCoordSeq_getX( (GEOSCoordSequence *)cs, i, x )
00143 #define GEOSCoordSeq_getY(cs,i,y) GEOSCoordSeq_getY( (GEOSCoordSequence *)cs, i, y )
00144
00145 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms );
00146
00147 static GEOSGeometry *cloneGeosGeom( const GEOSGeometry *geom )
00148 {
00149
00150
00151 int type = GEOSGeomTypeId(( GEOSGeometry * ) geom );
00152
00153 if ( type == GEOS_MULTIPOINT || type == GEOS_MULTILINESTRING || type == GEOS_MULTIPOLYGON )
00154 {
00155 QVector<GEOSGeometry *> geoms;
00156
00157 try
00158 {
00159
00160 for ( int i = 0; i < GEOSGetNumGeometries(( GEOSGeometry * )geom ); ++i )
00161 geoms << GEOSGeom_clone(( GEOSGeometry * ) GEOSGetGeometryN(( GEOSGeometry * ) geom, i ) );
00162
00163 return createGeosCollection( type, geoms );
00164 }
00165 catch ( GEOSException &e )
00166 {
00167 Q_UNUSED( e );
00168 for ( int i = 0; i < geoms.count(); i++ )
00169 GEOSGeom_destroy( geoms[i] );
00170
00171 return 0;
00172 }
00173 }
00174 else
00175 {
00176 return GEOSGeom_clone(( GEOSGeometry * ) geom );
00177 }
00178 }
00179
00180 #define GEOSGeom_clone(g) cloneGeosGeom(g)
00181 #endif
00182
00183 QgsGeometry::QgsGeometry()
00184 : mGeometry( 0 ),
00185 mGeometrySize( 0 ),
00186 mGeos( 0 ),
00187 mDirtyWkb( false ),
00188 mDirtyGeos( false )
00189 {
00190 }
00191
00192 QgsGeometry::QgsGeometry( QgsGeometry const & rhs )
00193 : mGeometry( 0 ),
00194 mGeometrySize( rhs.mGeometrySize ),
00195 mDirtyWkb( rhs.mDirtyWkb ),
00196 mDirtyGeos( rhs.mDirtyGeos )
00197 {
00198 if ( mGeometrySize && rhs.mGeometry )
00199 {
00200 mGeometry = new unsigned char[mGeometrySize];
00201 memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
00202 }
00203
00204
00205 if ( rhs.mGeos )
00206 {
00207 mGeos = GEOSGeom_clone( rhs.mGeos );
00208 }
00209 else
00210 {
00211 mGeos = 0;
00212 }
00213 }
00214
00216 QgsGeometry::~QgsGeometry()
00217 {
00218 if ( mGeometry )
00219 {
00220 delete [] mGeometry;
00221 }
00222
00223 if ( mGeos )
00224 {
00225 GEOSGeom_destroy( mGeos );
00226 }
00227 }
00228
00229 static unsigned int getNumGeosPoints( const GEOSGeometry *geom )
00230 {
00231 unsigned int n;
00232 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( geom );
00233 GEOSCoordSeq_getSize( cs, &n );
00234 return n;
00235 }
00236
00237 static GEOSGeometry *createGeosPoint( const QgsPoint &point )
00238 {
00239 GEOSCoordSequence *coord = GEOSCoordSeq_create( 1, 2 );
00240 GEOSCoordSeq_setX( coord, 0, point.x() );
00241 GEOSCoordSeq_setY( coord, 0, point.y() );
00242 return GEOSGeom_createPoint( coord );
00243 }
00244
00245 static GEOSCoordSequence *createGeosCoordSequence( const QgsPolyline& points )
00246 {
00247 GEOSCoordSequence *coord = 0;
00248
00249 try
00250 {
00251 coord = GEOSCoordSeq_create( points.count(), 2 );
00252 int i;
00253 for ( i = 0; i < points.count(); i++ )
00254 {
00255 GEOSCoordSeq_setX( coord, i, points[i].x() );
00256 GEOSCoordSeq_setY( coord, i, points[i].y() );
00257 }
00258 return coord;
00259 }
00260 catch ( GEOSException &e )
00261 {
00262 Q_UNUSED( e );
00263
00264
00265 throw;
00266 }
00267 }
00268
00269 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms )
00270 {
00271 GEOSGeometry **geomarr = new GEOSGeometry*[ geoms.size()];
00272 if ( !geomarr )
00273 return 0;
00274
00275 for ( int i = 0; i < geoms.size(); i++ )
00276 geomarr[i] = geoms[i];
00277
00278 GEOSGeometry *geom = 0;
00279
00280 try
00281 {
00282 geom = GEOSGeom_createCollection( typeId, geomarr, geoms.size() );
00283 }
00284 catch ( GEOSException &e )
00285 {
00286 Q_UNUSED( e );
00287 }
00288
00289 delete [] geomarr;
00290
00291 return geom;
00292 }
00293
00294 static GEOSGeometry *createGeosLineString( const QgsPolyline& polyline )
00295 {
00296 GEOSCoordSequence *coord = 0;
00297
00298 try
00299 {
00300 coord = createGeosCoordSequence( polyline );
00301 return GEOSGeom_createLineString( coord );
00302 }
00303 catch ( GEOSException &e )
00304 {
00305 Q_UNUSED( e );
00306
00307
00308
00309 return 0;
00310 }
00311 }
00312
00313 static GEOSGeometry *createGeosLinearRing( const QgsPolyline& polyline )
00314 {
00315 GEOSCoordSequence *coord = 0;
00316
00317 if ( polyline.count() == 0 )
00318 return 0;
00319
00320 try
00321 {
00322 if ( polyline[0] != polyline[polyline.size()-1] )
00323 {
00324
00325 QgsPolyline closed( polyline );
00326 closed << closed[0];
00327 coord = createGeosCoordSequence( closed );
00328 }
00329 else
00330 {
00331 coord = createGeosCoordSequence( polyline );
00332 }
00333
00334 return GEOSGeom_createLinearRing( coord );
00335 }
00336 catch ( GEOSException &e )
00337 {
00338 Q_UNUSED( e );
00339
00340
00341
00342 return 0;
00343 }
00344 }
00345
00346 static GEOSGeometry *createGeosPolygon( const QVector<GEOSGeometry*> &rings )
00347 {
00348 GEOSGeometry *shell = rings[0];
00349 GEOSGeometry **holes = NULL;
00350
00351 if ( rings.size() > 1 )
00352 {
00353 holes = new GEOSGeometry*[ rings.size()-1 ];
00354 if ( !holes )
00355 return 0;
00356
00357 for ( int i = 0; i < rings.size() - 1; i++ )
00358 holes[i] = rings[i+1];
00359 }
00360
00361 GEOSGeometry *geom = GEOSGeom_createPolygon( shell, holes, rings.size() - 1 );
00362
00363 if ( holes )
00364 delete [] holes;
00365
00366 return geom;
00367 }
00368
00369 static GEOSGeometry *createGeosPolygon( GEOSGeometry *shell )
00370 {
00371 return createGeosPolygon( QVector<GEOSGeometry*>() << shell );
00372 }
00373
00374 static GEOSGeometry *createGeosPolygon( const QgsPolygon& polygon )
00375 {
00376 if ( polygon.count() == 0 )
00377 return 0;
00378
00379 QVector<GEOSGeometry *> geoms;
00380
00381 try
00382 {
00383 for ( int i = 0; i < polygon.count(); i++ )
00384 geoms << createGeosLinearRing( polygon[i] );
00385
00386 return createGeosPolygon( geoms );
00387 }
00388 catch ( GEOSException &e )
00389 {
00390 Q_UNUSED( e );
00391 for ( int i = 0; i < geoms.count(); i++ )
00392 GEOSGeom_destroy( geoms[i] );
00393 return 0;
00394 }
00395 }
00396
00397 static QgsGeometry *fromGeosGeom( GEOSGeometry *geom )
00398 {
00399 if ( !geom )
00400 return 0;
00401
00402 QgsGeometry* g = new QgsGeometry;
00403 g->fromGeos( geom );
00404 return g;
00405 }
00406
00407 QgsGeometry* QgsGeometry::fromWkt( QString wkt )
00408 {
00409 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR>=3)
00410 GEOSWKTReader *reader = GEOSWKTReader_create();;
00411 QgsGeometry *g = fromGeosGeom( GEOSWKTReader_read( reader, wkt.toLocal8Bit().data() ) );
00412 GEOSWKTReader_destroy( reader );
00413 return g;
00414 #else
00415 return fromGeosGeom( GEOSGeomFromWKT( wkt.toLocal8Bit().data() ) );
00416 #endif
00417 }
00418
00419 QgsGeometry* QgsGeometry::fromPoint( const QgsPoint& point )
00420 {
00421 return fromGeosGeom( createGeosPoint( point ) );
00422 }
00423
00424 QgsGeometry* QgsGeometry::fromPolyline( const QgsPolyline& polyline )
00425 {
00426 return fromGeosGeom( createGeosLineString( polyline ) );
00427 }
00428
00429 QgsGeometry* QgsGeometry::fromPolygon( const QgsPolygon& polygon )
00430 {
00431 return fromGeosGeom( createGeosPolygon( polygon ) );
00432 }
00433
00434 QgsGeometry* QgsGeometry::fromMultiPoint( const QgsMultiPoint& multipoint )
00435 {
00436 QVector<GEOSGeometry *> geoms;
00437
00438 try
00439 {
00440 for ( int i = 0; i < multipoint.size(); ++i )
00441 {
00442 geoms << createGeosPoint( multipoint[i] );
00443 }
00444
00445 return fromGeosGeom( createGeosCollection( GEOS_MULTIPOINT, geoms ) );
00446 }
00447 catch ( GEOSException &e )
00448 {
00449 Q_UNUSED( e );
00450 for ( int i = 0; i < geoms.size(); ++i )
00451 GEOSGeom_destroy( geoms[i] );
00452
00453 return 0;
00454 }
00455 }
00456
00457 QgsGeometry* QgsGeometry::fromMultiPolyline( const QgsMultiPolyline& multiline )
00458 {
00459 QVector<GEOSGeometry *> geoms;
00460
00461 try
00462 {
00463 for ( int i = 0; i < multiline.count(); i++ )
00464 geoms << createGeosLineString( multiline[i] );
00465
00466 return fromGeosGeom( createGeosCollection( GEOS_MULTILINESTRING, geoms ) );
00467 }
00468 catch ( GEOSException &e )
00469 {
00470 Q_UNUSED( e );
00471 for ( int i = 0; i < geoms.count(); i++ )
00472 GEOSGeom_destroy( geoms[i] );
00473
00474 return 0;
00475 }
00476 }
00477
00478 QgsGeometry* QgsGeometry::fromMultiPolygon( const QgsMultiPolygon& multipoly )
00479 {
00480 if ( multipoly.count() == 0 )
00481 return 0;
00482
00483 QVector<GEOSGeometry *> geoms;
00484
00485 try
00486 {
00487 for ( int i = 0; i < multipoly.count(); i++ )
00488 geoms << createGeosPolygon( multipoly[i] );
00489
00490 return fromGeosGeom( createGeosCollection( GEOS_MULTIPOLYGON, geoms ) );
00491 }
00492 catch ( GEOSException &e )
00493 {
00494 Q_UNUSED( e );
00495 for ( int i = 0; i < geoms.count(); i++ )
00496 GEOSGeom_destroy( geoms[i] );
00497
00498 return 0;
00499 }
00500 }
00501
00502 QgsGeometry* QgsGeometry::fromRect( const QgsRectangle& rect )
00503 {
00504 QgsPolyline ring;
00505 ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
00506 ring.append( QgsPoint( rect.xMaximum(), rect.yMinimum() ) );
00507 ring.append( QgsPoint( rect.xMaximum(), rect.yMaximum() ) );
00508 ring.append( QgsPoint( rect.xMinimum(), rect.yMaximum() ) );
00509 ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
00510
00511 QgsPolygon polygon;
00512 polygon.append( ring );
00513
00514 return fromPolygon( polygon );
00515 }
00516
00517 QgsGeometry & QgsGeometry::operator=( QgsGeometry const & rhs )
00518 {
00519 if ( &rhs == this )
00520 {
00521 return *this;
00522 }
00523
00524
00525 if ( mGeometry )
00526 {
00527 delete [] mGeometry;
00528 mGeometry = 0;
00529 }
00530
00531 mGeometrySize = rhs.mGeometrySize;
00532
00533
00534 GEOSGeom_destroy( mGeos );
00535 mGeos = rhs.mGeos ? GEOSGeom_clone( rhs.mGeos ) : 0;
00536
00537 mDirtyGeos = rhs.mDirtyGeos;
00538 mDirtyWkb = rhs.mDirtyWkb;
00539
00540 if ( mGeometrySize && rhs.mGeometry )
00541 {
00542 mGeometry = new unsigned char[mGeometrySize];
00543 memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
00544 }
00545
00546 return *this;
00547 }
00548
00549
00550 void QgsGeometry::fromWkb( unsigned char * wkb, size_t length )
00551 {
00552
00553 if ( mGeometry )
00554 {
00555 delete [] mGeometry;
00556 mGeometry = 0;
00557 }
00558 if ( mGeos )
00559 {
00560 GEOSGeom_destroy( mGeos );
00561 mGeos = 0;
00562 }
00563
00564 mGeometry = wkb;
00565 mGeometrySize = length;
00566
00567 mDirtyWkb = false;
00568 mDirtyGeos = true;
00569 }
00570
00571 unsigned char * QgsGeometry::asWkb()
00572 {
00573 if ( mDirtyWkb )
00574 {
00575
00576 exportGeosToWkb();
00577 }
00578
00579 return mGeometry;
00580 }
00581
00582 size_t QgsGeometry::wkbSize()
00583 {
00584 if ( mDirtyWkb )
00585 {
00586 exportGeosToWkb();
00587 }
00588
00589 return mGeometrySize;
00590 }
00591
00592 GEOSGeometry* QgsGeometry::asGeos()
00593 {
00594 if ( mDirtyGeos )
00595 {
00596 if ( !exportWkbToGeos() )
00597 {
00598 return 0;
00599 }
00600 }
00601 return mGeos;
00602 }
00603
00604
00605 QGis::WkbType QgsGeometry::wkbType()
00606 {
00607 unsigned char *geom = asWkb();
00608 if ( geom )
00609 {
00610 unsigned int wkbType;
00611 memcpy( &wkbType, ( geom + 1 ), sizeof( wkbType ) );
00612 return ( QGis::WkbType ) wkbType;
00613 }
00614 else
00615 {
00616 return QGis::WKBUnknown;
00617 }
00618 }
00619
00620
00621 QGis::GeometryType QgsGeometry::type()
00622 {
00623 if ( mDirtyWkb )
00624 {
00625
00626 exportGeosToWkb();
00627 }
00628
00629 QGis::WkbType type = wkbType();
00630 if ( type == QGis::WKBPoint || type == QGis::WKBPoint25D ||
00631 type == QGis::WKBMultiPoint || type == QGis::WKBMultiPoint25D )
00632 return QGis::Point;
00633 if ( type == QGis::WKBLineString || type == QGis::WKBLineString25D ||
00634 type == QGis::WKBMultiLineString || type == QGis::WKBMultiLineString25D )
00635 return QGis::Line;
00636 if ( type == QGis::WKBPolygon || type == QGis::WKBPolygon25D ||
00637 type == QGis::WKBMultiPolygon || type == QGis::WKBMultiPolygon25D )
00638 return QGis::Polygon;
00639
00640 return QGis::UnknownGeometry;
00641 }
00642
00643 bool QgsGeometry::isMultipart()
00644 {
00645 if ( mDirtyWkb )
00646 {
00647
00648 exportGeosToWkb();
00649 }
00650
00651 QGis::WkbType type = wkbType();
00652 if ( type == QGis::WKBMultiPoint ||
00653 type == QGis::WKBMultiPoint25D ||
00654 type == QGis::WKBMultiLineString ||
00655 type == QGis::WKBMultiLineString25D ||
00656 type == QGis::WKBMultiPolygon ||
00657 type == QGis::WKBMultiPolygon25D )
00658 return true;
00659
00660 return false;
00661 }
00662
00663
00664 void QgsGeometry::fromGeos( GEOSGeometry* geos )
00665 {
00666
00667
00668 if ( mGeos )
00669 {
00670 GEOSGeom_destroy( mGeos );
00671 mGeos = 0;
00672 }
00673 if ( mGeometry )
00674 {
00675 delete [] mGeometry;
00676 mGeometry = 0;
00677 }
00678
00679 mGeos = geos;
00680
00681 mDirtyWkb = true;
00682 mDirtyGeos = false;
00683 }
00684
00685 QgsPoint QgsGeometry::closestVertex( const QgsPoint& point, int& atVertex, int& beforeVertex, int& afterVertex, double& sqrDist )
00686 {
00687
00688 if ( mDirtyWkb )
00689 {
00690 exportGeosToWkb();
00691 }
00692
00693 if ( !mGeometry )
00694 {
00695 QgsDebugMsg( "WKB geometry not available!" );
00696 return QgsPoint( 0, 0 );
00697 }
00698
00699 int vertexnr = -1;
00700 int vertexcounter = 0;
00701 QGis::WkbType wkbType;
00702 double actdist = std::numeric_limits<double>::max();
00703 double x = 0;
00704 double y = 0;
00705 double *tempx, *tempy;
00706 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
00707 beforeVertex = -1;
00708 afterVertex = -1;
00709 bool hasZValue = false;
00710
00711 switch ( wkbType )
00712 {
00713 case QGis::WKBPoint25D:
00714 case QGis::WKBPoint:
00715 {
00716 x = *(( double * )( mGeometry + 5 ) );
00717 y = *(( double * )( mGeometry + 5 + sizeof( double ) ) );
00718 actdist = point.sqrDist( x, y );
00719 vertexnr = 0;
00720 break;
00721 }
00722 case QGis::WKBLineString25D:
00723 hasZValue = true;
00724
00725
00726
00727 case QGis::WKBLineString:
00728 {
00729 unsigned char* ptr = mGeometry + 5;
00730 int* npoints = ( int* )ptr;
00731 ptr += sizeof( int );
00732 for ( int index = 0; index < *npoints; ++index )
00733 {
00734 tempx = ( double* )ptr;
00735 ptr += sizeof( double );
00736 tempy = ( double* )ptr;
00737 if ( point.sqrDist( *tempx, *tempy ) < actdist )
00738 {
00739 x = *tempx;
00740 y = *tempy;
00741 actdist = point.sqrDist( *tempx, *tempy );
00742 vertexnr = index;
00743 if ( index == 0 )
00744 {
00745 beforeVertex = -1;
00746 }
00747 else
00748 {
00749 beforeVertex = index - 1;
00750 }
00751 if ( index == ( *npoints - 1 ) )
00752 {
00753 afterVertex = -1;
00754 }
00755 else
00756 {
00757 afterVertex = index + 1;
00758 }
00759 }
00760 ptr += sizeof( double );
00761 if ( hasZValue )
00762 {
00763 ptr += sizeof( double );
00764 }
00765 }
00766 break;
00767 }
00768 case QGis::WKBPolygon25D:
00769 hasZValue = true;
00770 case QGis::WKBPolygon:
00771 {
00772 int* nrings = ( int* )( mGeometry + 5 );
00773 int* npoints;
00774 unsigned char* ptr = mGeometry + 9;
00775 for ( int index = 0; index < *nrings; ++index )
00776 {
00777 npoints = ( int* )ptr;
00778 ptr += sizeof( int );
00779 for ( int index2 = 0; index2 < *npoints; ++index2 )
00780 {
00781 tempx = ( double* )ptr;
00782 ptr += sizeof( double );
00783 tempy = ( double* )ptr;
00784 if ( point.sqrDist( *tempx, *tempy ) < actdist )
00785 {
00786 x = *tempx;
00787 y = *tempy;
00788 actdist = point.sqrDist( *tempx, *tempy );
00789 vertexnr = vertexcounter;
00790
00791 if ( index2 == 0 )
00792 {
00793 beforeVertex = vertexcounter + ( *npoints - 2 );
00794 afterVertex = vertexcounter + 1;
00795 }
00796 else if ( index2 == ( *npoints - 1 ) )
00797 {
00798 beforeVertex = vertexcounter - 1;
00799 afterVertex = vertexcounter - ( *npoints - 2 );
00800 }
00801 else
00802 {
00803 beforeVertex = vertexcounter - 1;
00804 afterVertex = vertexcounter + 1;
00805 }
00806 }
00807 ptr += sizeof( double );
00808 if ( hasZValue )
00809 {
00810 ptr += sizeof( double );
00811 }
00812 ++vertexcounter;
00813 }
00814 }
00815 break;
00816 }
00817 case QGis::WKBMultiPoint25D:
00818 hasZValue = true;
00819 case QGis::WKBMultiPoint:
00820 {
00821 unsigned char* ptr = mGeometry + 5;
00822 int* npoints = ( int* )ptr;
00823 ptr += sizeof( int );
00824 for ( int index = 0; index < *npoints; ++index )
00825 {
00826 ptr += ( 1 + sizeof( int ) );
00827 tempx = ( double* )ptr;
00828 tempy = ( double* )( ptr + sizeof( double ) );
00829 if ( point.sqrDist( *tempx, *tempy ) < actdist )
00830 {
00831 x = *tempx;
00832 y = *tempy;
00833 actdist = point.sqrDist( *tempx, *tempy );
00834 vertexnr = index;
00835 }
00836 ptr += ( 2 * sizeof( double ) );
00837 if ( hasZValue )
00838 {
00839 ptr += sizeof( double );
00840 }
00841 }
00842 break;
00843 }
00844 case QGis::WKBMultiLineString25D:
00845 hasZValue = true;
00846 case QGis::WKBMultiLineString:
00847 {
00848 unsigned char* ptr = mGeometry + 5;
00849 int* nlines = ( int* )ptr;
00850 int* npoints = 0;
00851 ptr += sizeof( int );
00852 for ( int index = 0; index < *nlines; ++index )
00853 {
00854 ptr += ( sizeof( int ) + 1 );
00855 npoints = ( int* )ptr;
00856 ptr += sizeof( int );
00857 for ( int index2 = 0; index2 < *npoints; ++index2 )
00858 {
00859 tempx = ( double* )ptr;
00860 ptr += sizeof( double );
00861 tempy = ( double* )ptr;
00862 ptr += sizeof( double );
00863 if ( point.sqrDist( *tempx, *tempy ) < actdist )
00864 {
00865 x = *tempx;
00866 y = *tempy;
00867 actdist = point.sqrDist( *tempx, *tempy );
00868 vertexnr = vertexcounter;
00869
00870 if ( index2 == 0 )
00871 {
00872 beforeVertex = -1;
00873 }
00874 else
00875 {
00876 beforeVertex = vertexnr - 1;
00877 }
00878 if ( index2 == ( *npoints ) - 1 )
00879 {
00880 afterVertex = -1;
00881 }
00882 else
00883 {
00884 afterVertex = vertexnr + 1;
00885 }
00886 }
00887 if ( hasZValue )
00888 {
00889 ptr += sizeof( double );
00890 }
00891 ++vertexcounter;
00892 }
00893 }
00894 break;
00895 }
00896 case QGis::WKBMultiPolygon25D:
00897 hasZValue = true;
00898 case QGis::WKBMultiPolygon:
00899 {
00900 unsigned char* ptr = mGeometry + 5;
00901 int* npolys = ( int* )ptr;
00902 int* nrings;
00903 int* npoints;
00904 ptr += sizeof( int );
00905 for ( int index = 0; index < *npolys; ++index )
00906 {
00907 ptr += ( 1 + sizeof( int ) );
00908 nrings = ( int* )ptr;
00909 ptr += sizeof( int );
00910 for ( int index2 = 0; index2 < *nrings; ++index2 )
00911 {
00912 npoints = ( int* )ptr;
00913 ptr += sizeof( int );
00914 for ( int index3 = 0; index3 < *npoints; ++index3 )
00915 {
00916 tempx = ( double* )ptr;
00917 ptr += sizeof( double );
00918 tempy = ( double* )ptr;
00919 if ( point.sqrDist( *tempx, *tempy ) < actdist )
00920 {
00921 x = *tempx;
00922 y = *tempy;
00923 actdist = point.sqrDist( *tempx, *tempy );
00924 vertexnr = vertexcounter;
00925
00926
00927 if ( index3 == 0 )
00928 {
00929 beforeVertex = vertexcounter + ( *npoints - 2 );
00930 afterVertex = vertexcounter + 1;
00931 }
00932 else if ( index3 == ( *npoints - 1 ) )
00933 {
00934 beforeVertex = vertexcounter - 1;
00935 afterVertex = vertexcounter - ( *npoints - 2 );
00936 }
00937 else
00938 {
00939 beforeVertex = vertexcounter - 1;
00940 afterVertex = vertexcounter + 1;
00941 }
00942 }
00943 ptr += sizeof( double );
00944 if ( hasZValue )
00945 {
00946 ptr += sizeof( double );
00947 }
00948 ++vertexcounter;
00949 }
00950 }
00951 }
00952 break;
00953 }
00954 default:
00955 break;
00956 }
00957 sqrDist = actdist;
00958 atVertex = vertexnr;
00959 return QgsPoint( x, y );
00960 }
00961
00962
00963 void QgsGeometry::adjacentVertices( int atVertex, int& beforeVertex, int& afterVertex )
00964 {
00965
00966 if ( mDirtyWkb )
00967 {
00968 exportGeosToWkb();
00969 }
00970
00971 beforeVertex = -1;
00972 afterVertex = -1;
00973
00974 if ( !mGeometry )
00975 {
00976 QgsDebugMsg( "WKB geometry not available!" );
00977 return;
00978 }
00979
00980 int vertexcounter = 0;
00981
00982 QGis::WkbType wkbType;
00983 bool hasZValue = false;
00984
00985 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
00986
00987 switch ( wkbType )
00988 {
00989 case QGis::WKBPoint:
00990 {
00991
00992 break;
00993 }
00994 case QGis::WKBLineString25D:
00995 case QGis::WKBLineString:
00996 {
00997 unsigned char* ptr = mGeometry + 5;
00998 int* npoints = ( int* ) ptr;
00999
01000 const int index = atVertex;
01001
01002
01003
01004 if ( index == 0 )
01005 {
01006 beforeVertex = -1;
01007 }
01008 else
01009 {
01010 beforeVertex = index - 1;
01011 }
01012
01013 if ( index == ( *npoints - 1 ) )
01014 {
01015 afterVertex = -1;
01016 }
01017 else
01018 {
01019 afterVertex = index + 1;
01020 }
01021
01022 break;
01023 }
01024 case QGis::WKBPolygon25D:
01025 hasZValue = true;
01026 case QGis::WKBPolygon:
01027 {
01028 int* nrings = ( int* )( mGeometry + 5 );
01029 int* npoints;
01030 unsigned char* ptr = mGeometry + 9;
01031
01032
01033
01034 for ( int index0 = 0; index0 < *nrings; ++index0 )
01035 {
01036 npoints = ( int* )ptr;
01037 ptr += sizeof( int );
01038
01039 for ( int index1 = 0; index1 < *npoints; ++index1 )
01040 {
01041 ptr += sizeof( double );
01042 ptr += sizeof( double );
01043 if ( hasZValue )
01044 {
01045 ptr += sizeof( double );
01046 }
01047 if ( vertexcounter == atVertex )
01048 {
01049 if ( index1 == 0 )
01050 {
01051 beforeVertex = vertexcounter + ( *npoints - 2 );
01052 afterVertex = vertexcounter + 1;
01053 }
01054 else if ( index1 == ( *npoints - 1 ) )
01055 {
01056 beforeVertex = vertexcounter - 1;
01057 afterVertex = vertexcounter - ( *npoints - 2 );
01058 }
01059 else
01060 {
01061 beforeVertex = vertexcounter - 1;
01062 afterVertex = vertexcounter + 1;
01063 }
01064 }
01065
01066 ++vertexcounter;
01067 }
01068 }
01069 break;
01070 }
01071 case QGis::WKBMultiPoint25D:
01072 case QGis::WKBMultiPoint:
01073 {
01074
01075 break;
01076 }
01077
01078 case QGis::WKBMultiLineString25D:
01079 hasZValue = true;
01080 case QGis::WKBMultiLineString:
01081 {
01082 unsigned char* ptr = mGeometry + 5;
01083 int* nlines = ( int* )ptr;
01084 int* npoints = 0;
01085 ptr += sizeof( int );
01086
01087 for ( int index0 = 0; index0 < *nlines; ++index0 )
01088 {
01089 ptr += ( sizeof( int ) + 1 );
01090 npoints = ( int* )ptr;
01091 ptr += sizeof( int );
01092
01093 for ( int index1 = 0; index1 < *npoints; ++index1 )
01094 {
01095 ptr += sizeof( double );
01096 ptr += sizeof( double );
01097 if ( hasZValue )
01098 {
01099 ptr += sizeof( double );
01100 }
01101
01102 if ( vertexcounter == atVertex )
01103 {
01104
01105 if ( index1 == 0 )
01106 {
01107 beforeVertex = -1;
01108 }
01109 else
01110 {
01111 beforeVertex = vertexcounter - 1;
01112 }
01113 if ( index1 == ( *npoints ) - 1 )
01114 {
01115 afterVertex = -1;
01116 }
01117 else
01118 {
01119 afterVertex = vertexcounter + 1;
01120 }
01121 }
01122 ++vertexcounter;
01123 }
01124 }
01125 break;
01126 }
01127 case QGis::WKBMultiPolygon25D:
01128 hasZValue = true;
01129 case QGis::WKBMultiPolygon:
01130 {
01131 unsigned char* ptr = mGeometry + 5;
01132 int* npolys = ( int* )ptr;
01133 int* nrings;
01134 int* npoints;
01135 ptr += sizeof( int );
01136
01137 for ( int index0 = 0; index0 < *npolys; ++index0 )
01138 {
01139 ptr += ( 1 + sizeof( int ) );
01140 nrings = ( int* )ptr;
01141 ptr += sizeof( int );
01142
01143 for ( int index1 = 0; index1 < *nrings; ++index1 )
01144 {
01145 npoints = ( int* )ptr;
01146 ptr += sizeof( int );
01147
01148 for ( int index2 = 0; index2 < *npoints; ++index2 )
01149 {
01150 ptr += sizeof( double );
01151 ptr += sizeof( double );
01152 if ( hasZValue )
01153 {
01154 ptr += sizeof( double );
01155 }
01156 if ( vertexcounter == atVertex )
01157 {
01158
01159
01160
01161 if ( index2 == 0 )
01162 {
01163 beforeVertex = vertexcounter + ( *npoints - 2 );
01164 afterVertex = vertexcounter + 1;
01165 }
01166 else if ( index2 == ( *npoints - 1 ) )
01167 {
01168 beforeVertex = vertexcounter - 1;
01169 afterVertex = vertexcounter - ( *npoints - 2 );
01170 }
01171 else
01172 {
01173 beforeVertex = vertexcounter - 1;
01174 afterVertex = vertexcounter + 1;
01175 }
01176 }
01177 ++vertexcounter;
01178 }
01179 }
01180 }
01181
01182 break;
01183 }
01184
01185 default:
01186 break;
01187 }
01188 }
01189
01190
01191
01192 bool QgsGeometry::insertVertex( double x, double y,
01193 int beforeVertex,
01194 const GEOSCoordSequence* old_sequence,
01195 GEOSCoordSequence** new_sequence )
01196
01197 {
01198
01199 if ( beforeVertex < 0 )
01200 {
01201 *new_sequence = 0;
01202 return false;
01203 }
01204
01205 unsigned int numPoints;
01206 GEOSCoordSeq_getSize( old_sequence, &numPoints );
01207
01208 *new_sequence = GEOSCoordSeq_create( numPoints + 1, 2 );
01209 if ( !*new_sequence )
01210 return false;
01211
01212 bool inserted = false;
01213 for ( unsigned int i = 0, j = 0; i < numPoints; i++, j++ )
01214 {
01215
01216 if ( beforeVertex == static_cast<int>( i ) )
01217 {
01218 GEOSCoordSeq_setX( *new_sequence, j, x );
01219 GEOSCoordSeq_setY( *new_sequence, j, y );
01220 j++;
01221 inserted = true;
01222 }
01223
01224 double aX, aY;
01225 GEOSCoordSeq_getX( old_sequence, i, &aX );
01226 GEOSCoordSeq_getY( old_sequence, i, &aY );
01227
01228 GEOSCoordSeq_setX( *new_sequence, j, aX );
01229 GEOSCoordSeq_setY( *new_sequence, j, aY );
01230 }
01231
01232 if ( !inserted )
01233 {
01234
01235
01236 GEOSCoordSeq_setX( *new_sequence, numPoints, x );
01237 GEOSCoordSeq_setY( *new_sequence, numPoints, y );
01238 }
01239
01240
01241 return inserted;
01242 }
01243
01244 bool QgsGeometry::moveVertex( double x, double y, int atVertex )
01245 {
01246 int vertexnr = atVertex;
01247
01248
01249 if ( mDirtyWkb )
01250 {
01251 exportGeosToWkb();
01252 }
01253
01254 if ( !mGeometry )
01255 {
01256 QgsDebugMsg( "WKB geometry not available!" );
01257 return false;
01258 }
01259
01260 QGis::WkbType wkbType;
01261 bool hasZValue = false;
01262 unsigned char* ptr = mGeometry + 1;
01263 memcpy( &wkbType, ptr, sizeof( wkbType ) );
01264 ptr += sizeof( wkbType );
01265
01266 switch ( wkbType )
01267 {
01268 case QGis::WKBPoint25D:
01269 case QGis::WKBPoint:
01270 {
01271 if ( vertexnr == 0 )
01272 {
01273 memcpy( ptr, &x, sizeof( double ) );
01274 ptr += sizeof( double );
01275 memcpy( ptr, &y, sizeof( double ) );
01276 mDirtyGeos = true;
01277 return true;
01278 }
01279 else
01280 {
01281 return false;
01282 }
01283 }
01284 case QGis::WKBMultiPoint25D:
01285 hasZValue = true;
01286 case QGis::WKBMultiPoint:
01287 {
01288 int* nrPoints = ( int* )ptr;
01289 if ( vertexnr > *nrPoints || vertexnr < 0 )
01290 {
01291 return false;
01292 }
01293 ptr += sizeof( int );
01294 if ( hasZValue )
01295 {
01296 ptr += ( 3 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr;
01297 }
01298 else
01299 {
01300 ptr += ( 2 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr;
01301 }
01302 ptr += ( 1 + sizeof( int ) );
01303 memcpy( ptr, &x, sizeof( double ) );
01304 ptr += sizeof( double );
01305 memcpy( ptr, &y, sizeof( double ) );
01306 mDirtyGeos = true;
01307 return true;
01308 }
01309 case QGis::WKBLineString25D:
01310 hasZValue = true;
01311 case QGis::WKBLineString:
01312 {
01313 int* nrPoints = ( int* )ptr;
01314 if ( vertexnr > *nrPoints || vertexnr < 0 )
01315 {
01316 return false;
01317 }
01318 ptr += sizeof( int );
01319 ptr += 2 * sizeof( double ) * vertexnr;
01320 if ( hasZValue )
01321 {
01322 ptr += sizeof( double ) * vertexnr;
01323 }
01324 memcpy( ptr, &x, sizeof( double ) );
01325 ptr += sizeof( double );
01326 memcpy( ptr, &y, sizeof( double ) );
01327 mDirtyGeos = true;
01328 return true;
01329 }
01330 case QGis::WKBMultiLineString25D:
01331 hasZValue = true;
01332 case QGis::WKBMultiLineString:
01333 {
01334 int* nrLines = ( int* )ptr;
01335 ptr += sizeof( int );
01336 int* nrPoints = 0;
01337 int pointindex = 0;
01338 for ( int linenr = 0; linenr < *nrLines; ++linenr )
01339 {
01340 ptr += sizeof( int ) + 1;
01341 nrPoints = ( int* )ptr;
01342 ptr += sizeof( int );
01343 if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nrPoints ) )
01344 {
01345 ptr += ( vertexnr - pointindex ) * 2 * sizeof( double );
01346 if ( hasZValue )
01347 {
01348 ptr += ( vertexnr - pointindex ) * sizeof( double );
01349 }
01350 memcpy( ptr, &x, sizeof( double ) );
01351 memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
01352 mDirtyGeos = true;
01353 return true;
01354 }
01355 pointindex += ( *nrPoints );
01356 ptr += 2 * sizeof( double ) * ( *nrPoints );
01357 if ( hasZValue )
01358 {
01359 ptr += sizeof( double ) * ( *nrPoints );
01360 }
01361 }
01362 return false;
01363 }
01364 case QGis::WKBPolygon25D:
01365 hasZValue = true;
01366 case QGis::WKBPolygon:
01367 {
01368 int* nrRings = ( int* )ptr;
01369 ptr += sizeof( int );
01370 int* nrPoints = 0;
01371 int pointindex = 0;
01372
01373 for ( int ringnr = 0; ringnr < *nrRings; ++ringnr )
01374 {
01375 nrPoints = ( int* )ptr;
01376 ptr += sizeof( int );
01377 if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )
01378 {
01379 memcpy( ptr, &x, sizeof( double ) );
01380 memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
01381 if ( hasZValue )
01382 {
01383 memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
01384 }
01385 else
01386 {
01387 memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
01388 }
01389 if ( hasZValue )
01390 {
01391 memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
01392 }
01393 else
01394 {
01395 memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
01396 }
01397 mDirtyGeos = true;
01398 return true;
01399 }
01400 else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )
01401 {
01402 ptr += 2 * sizeof( double ) * ( vertexnr - pointindex );
01403 if ( hasZValue )
01404 {
01405 ptr += sizeof( double ) * ( vertexnr - pointindex );
01406 }
01407 memcpy( ptr, &x, sizeof( double ) );
01408 ptr += sizeof( double );
01409 memcpy( ptr, &y, sizeof( double ) );
01410 mDirtyGeos = true;
01411 return true;
01412 }
01413 ptr += 2 * sizeof( double ) * ( *nrPoints );
01414 if ( hasZValue )
01415 {
01416 ptr += sizeof( double ) * ( *nrPoints );
01417 }
01418 pointindex += *nrPoints;
01419 }
01420 return false;
01421 }
01422 case QGis::WKBMultiPolygon25D:
01423 hasZValue = true;
01424 case QGis::WKBMultiPolygon:
01425 {
01426 int* nrPolygons = ( int* )ptr;
01427 ptr += sizeof( int );
01428 int* nrRings = 0;
01429 int* nrPoints = 0;
01430 int pointindex = 0;
01431
01432 for ( int polynr = 0; polynr < *nrPolygons; ++polynr )
01433 {
01434 ptr += ( 1 + sizeof( int ) );
01435 nrRings = ( int* )ptr;
01436 ptr += sizeof( int );
01437 for ( int ringnr = 0; ringnr < *nrRings; ++ringnr )
01438 {
01439 nrPoints = ( int* )ptr;
01440 ptr += sizeof( int );
01441 if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )
01442 {
01443 memcpy( ptr, &x, sizeof( double ) );
01444 memcpy( ptr + sizeof( double ), &y, sizeof( double ) );
01445 if ( hasZValue )
01446 {
01447 memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
01448 }
01449 else
01450 {
01451 memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) );
01452 }
01453 if ( hasZValue )
01454 {
01455 memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
01456 }
01457 else
01458 {
01459 memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) );
01460 }
01461 mDirtyGeos = true;
01462 return true;
01463 }
01464 else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )
01465 {
01466 ptr += 2 * sizeof( double ) * ( vertexnr - pointindex );
01467 if ( hasZValue )
01468 {
01469 ptr += 3 * sizeof( double ) * ( vertexnr - pointindex );
01470 }
01471 memcpy( ptr, &x, sizeof( double ) );
01472 ptr += sizeof( double );
01473 memcpy( ptr, &y, sizeof( double ) );
01474 mDirtyGeos = true;
01475 return true;
01476 }
01477 ptr += 2 * sizeof( double ) * ( *nrPoints );
01478 if ( hasZValue )
01479 {
01480 ptr += sizeof( double ) * ( *nrPoints );
01481 }
01482 pointindex += *nrPoints;
01483 }
01484 }
01485 return false;
01486 }
01487
01488 default:
01489 return false;
01490 }
01491 }
01492
01493 bool QgsGeometry::deleteVertex( int atVertex )
01494 {
01495 int vertexnr = atVertex;
01496 bool success = false;
01497
01498
01499 if ( mDirtyWkb )
01500 {
01501 exportGeosToWkb();
01502 }
01503
01504 if ( !mGeometry )
01505 {
01506 QgsDebugMsg( "WKB geometry not available!" );
01507 return false;
01508 }
01509
01510
01511 unsigned char* newbuffer;
01512 int pointindex = 0;
01513 QGis::WkbType wkbType;
01514 bool hasZValue = false;
01515 unsigned char* ptr = mGeometry + 1;
01516 memcpy( &wkbType, ptr, sizeof( wkbType ) );
01517
01518 switch ( wkbType )
01519 {
01520 case QGis::WKBPoint25D:
01521 case QGis::WKBLineString25D:
01522 case QGis::WKBPolygon25D:
01523 case QGis::WKBMultiPoint25D:
01524 case QGis::WKBMultiLineString25D:
01525 case QGis::WKBMultiPolygon25D:
01526 newbuffer = new unsigned char[mGeometrySize-3*sizeof( double )];
01527 break;
01528 default:
01529 newbuffer = new unsigned char[mGeometrySize-2*sizeof( double )];
01530 }
01531
01532 memcpy( newbuffer, mGeometry, 1 + sizeof( int ) );
01533
01534 ptr += sizeof( wkbType );
01535 unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int );
01536
01537 switch ( wkbType )
01538 {
01539 case QGis::WKBPoint25D:
01540 case QGis::WKBPoint:
01541 {
01542 break;
01543 }
01544 case QGis::WKBMultiPoint25D:
01545 hasZValue = true;
01546 case QGis::WKBMultiPoint:
01547 {
01548
01549 break;
01550 }
01551 case QGis::WKBLineString25D:
01552 hasZValue = true;
01553 case QGis::WKBLineString:
01554 {
01555 int* nPoints = ( int* )ptr;
01556 if (( *nPoints ) < 3 || vertexnr > ( *nPoints ) - 1 || vertexnr < 0 )
01557 {
01558 delete newbuffer;
01559 return false;
01560 }
01561 int newNPoints = ( *nPoints ) - 1;
01562 memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
01563 ptr += sizeof( int );
01564 newBufferPtr += sizeof( int );
01565 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01566 {
01567 if ( vertexnr != pointindex )
01568 {
01569 memcpy( newBufferPtr, ptr, sizeof( double ) );
01570 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01571 newBufferPtr += 2 * sizeof( double );
01572 if ( hasZValue )
01573 {
01574 newBufferPtr += sizeof( double );
01575 }
01576 }
01577 else
01578 {
01579 success = true;
01580 }
01581 ptr += 2 * sizeof( double );
01582 if ( hasZValue )
01583 {
01584 ptr += sizeof( double );
01585 }
01586 ++pointindex;
01587 }
01588 break;
01589 }
01590 case QGis::WKBMultiLineString25D:
01591 hasZValue = true;
01592 case QGis::WKBMultiLineString:
01593 {
01594 int* nLines = ( int* )ptr;
01595 memcpy( newBufferPtr, nLines, sizeof( int ) );
01596 newBufferPtr += sizeof( int );
01597 ptr += sizeof( int );
01598 int* nPoints = 0;
01599 int pointindex = 0;
01600 for ( int linenr = 0; linenr < *nLines; ++linenr )
01601 {
01602 memcpy( newBufferPtr, ptr, sizeof( int ) + 1 );
01603 ptr += ( sizeof( int ) + 1 );
01604 newBufferPtr += ( sizeof( int ) + 1 );
01605 nPoints = ( int* )ptr;
01606 ptr += sizeof( int );
01607 int newNPoint;
01608
01609
01610 if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nPoints ) )
01611 {
01612 if ( *nPoints < 3 )
01613 {
01614 delete newbuffer;
01615 return false;
01616 }
01617 newNPoint = ( *nPoints ) - 1;
01618 }
01619 else
01620 {
01621 newNPoint = *nPoints;
01622 }
01623 memcpy( newBufferPtr, &newNPoint, sizeof( int ) );
01624 newBufferPtr += sizeof( int );
01625
01626 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01627 {
01628 if ( vertexnr != pointindex )
01629 {
01630 memcpy( newBufferPtr, ptr, sizeof( double ) );
01631 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01632 newBufferPtr += 2 * sizeof( double );
01633 if ( hasZValue )
01634 {
01635 newBufferPtr += sizeof( double );
01636 }
01637 }
01638 else
01639 {
01640 success = true;
01641 }
01642 ptr += 2 * sizeof( double );
01643 if ( hasZValue )
01644 {
01645 ptr += sizeof( double );
01646 }
01647 ++pointindex;
01648 }
01649 }
01650 break;
01651 }
01652 case QGis::WKBPolygon25D:
01653 hasZValue = true;
01654 case QGis::WKBPolygon:
01655 {
01656 int* nRings = ( int* )ptr;
01657 memcpy( newBufferPtr, nRings, sizeof( int ) );
01658 ptr += sizeof( int );
01659 newBufferPtr += sizeof( int );
01660 int* nPoints = 0;
01661 int pointindex = 0;
01662
01663 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
01664 {
01665 nPoints = ( int* )ptr;
01666 ptr += sizeof( int );
01667 int newNPoints;
01668 if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )
01669 {
01670 if ( *nPoints < 5 )
01671 {
01672 delete newbuffer;
01673 return false;
01674 }
01675 newNPoints = *nPoints - 1;
01676 }
01677 else
01678 {
01679 newNPoints = *nPoints;
01680 }
01681 memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
01682 newBufferPtr += sizeof( int );
01683
01684 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01685 {
01686 if ( vertexnr != pointindex )
01687 {
01688 memcpy( newBufferPtr, ptr, sizeof( double ) );
01689 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01690 newBufferPtr += 2 * sizeof( double );
01691 if ( hasZValue )
01692 {
01693 newBufferPtr += sizeof( double );
01694 }
01695 }
01696 else
01697 {
01698 if ( pointnr == 0 )
01699 {
01700 memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) );
01701 memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) );
01702 }
01703 if ( pointnr == ( *nPoints ) - 1 )
01704 {
01705 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) );
01706 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) );
01707 }
01708 success = true;
01709 }
01710 ptr += 2 * sizeof( double );
01711 if ( hasZValue )
01712 {
01713 ptr += sizeof( double );
01714 }
01715 ++pointindex;
01716 }
01717 }
01718 break;
01719 }
01720 case QGis::WKBMultiPolygon25D:
01721 hasZValue = true;
01722 case QGis::WKBMultiPolygon:
01723 {
01724 int* nPolys = ( int* )ptr;
01725 memcpy( newBufferPtr, nPolys, sizeof( int ) );
01726 newBufferPtr += sizeof( int );
01727 ptr += sizeof( int );
01728 int* nRings = 0;
01729 int* nPoints = 0;
01730 int pointindex = 0;
01731
01732 for ( int polynr = 0; polynr < *nPolys; ++polynr )
01733 {
01734 memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) );
01735 ptr += ( 1 + sizeof( int ) );
01736 newBufferPtr += ( 1 + sizeof( int ) );
01737 nRings = ( int* )ptr;
01738 memcpy( newBufferPtr, nRings, sizeof( int ) );
01739 newBufferPtr += sizeof( int );
01740 ptr += sizeof( int );
01741 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
01742 {
01743 nPoints = ( int* )ptr;
01744 ptr += sizeof( int );
01745 int newNPoints;
01746 if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )
01747 {
01748 if ( *nPoints < 5 )
01749 {
01750 delete newbuffer;
01751 return false;
01752 }
01753 newNPoints = *nPoints - 1;
01754 }
01755 else
01756 {
01757 newNPoints = *nPoints;
01758 }
01759 memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
01760 newBufferPtr += sizeof( int );
01761
01762 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01763 {
01764 if ( vertexnr != pointindex )
01765 {
01766 memcpy( newBufferPtr, ptr, sizeof( double ) );
01767 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01768 newBufferPtr += 2 * sizeof( double );
01769 if ( hasZValue )
01770 {
01771 newBufferPtr += sizeof( double );
01772 }
01773 }
01774 else
01775 {
01776 if ( pointnr == 0 )
01777 {
01778 memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) );
01779 memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) );
01780 }
01781 if ( pointnr == ( *nPoints ) - 1 )
01782 {
01783 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) );
01784 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) );
01785 }
01786 success = true;
01787 }
01788 ptr += 2 * sizeof( double );
01789 if ( hasZValue )
01790 {
01791 ptr += sizeof( double );
01792 }
01793 ++pointindex;
01794 }
01795 }
01796 }
01797 break;
01798 }
01799 case QGis::WKBNoGeometry:
01800 case QGis::WKBUnknown:
01801 break;
01802 }
01803 if ( success )
01804 {
01805 delete [] mGeometry;
01806 mGeometry = newbuffer;
01807 mGeometrySize -= ( 2 * sizeof( double ) );
01808 if ( hasZValue )
01809 {
01810 mGeometrySize -= sizeof( double );
01811 }
01812 mDirtyGeos = true;
01813 return true;
01814 }
01815 else
01816 {
01817 delete[] newbuffer;
01818 return false;
01819 }
01820 }
01821
01822 bool QgsGeometry::insertVertex( double x, double y, int beforeVertex )
01823 {
01824 int vertexnr = beforeVertex;
01825 bool success = false;
01826
01827
01828 if ( mDirtyWkb )
01829 {
01830 exportGeosToWkb();
01831 }
01832
01833 if ( !mGeometry )
01834 {
01835 QgsDebugMsg( "WKB geometry not available!" );
01836 return false;
01837 }
01838
01839
01840 unsigned char* newbuffer;
01841
01842 int pointindex = 0;
01843 QGis::WkbType wkbType;
01844 bool hasZValue = false;
01845
01846 unsigned char* ptr = mGeometry + 1;
01847 memcpy( &wkbType, ptr, sizeof( wkbType ) );
01848
01849 switch ( wkbType )
01850 {
01851 case QGis::WKBPoint25D:
01852 case QGis::WKBLineString25D:
01853 case QGis::WKBPolygon25D:
01854 case QGis::WKBMultiPoint25D:
01855 case QGis::WKBMultiLineString25D:
01856 case QGis::WKBMultiPolygon25D:
01857 newbuffer = new unsigned char[mGeometrySize+3*sizeof( double )];
01858 break;
01859 default:
01860 newbuffer = new unsigned char[mGeometrySize+2*sizeof( double )];
01861 }
01862 memcpy( newbuffer, mGeometry, 1 + sizeof( int ) );
01863
01864 ptr += sizeof( wkbType );
01865 unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int );
01866
01867 switch ( wkbType )
01868 {
01869 case QGis::WKBPoint25D:
01870 case QGis::WKBPoint:
01871 {
01872 delete newbuffer;
01873 return false;
01874 }
01875 case QGis::WKBMultiPoint25D:
01876 hasZValue = true;
01877 case QGis::WKBMultiPoint:
01878 {
01879
01880 break;
01881 }
01882 case QGis::WKBLineString25D:
01883 hasZValue = true;
01884 case QGis::WKBLineString:
01885 {
01886 int* nPoints = ( int* )ptr;
01887 if ( vertexnr > *nPoints || vertexnr < 0 )
01888 {
01889 break;
01890 }
01891 int newNPoints = ( *nPoints ) + 1;
01892 memcpy( newBufferPtr, &newNPoints, sizeof( int ) );
01893 newBufferPtr += sizeof( int );
01894 ptr += sizeof( int );
01895
01896 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01897 {
01898 memcpy( newBufferPtr, ptr, sizeof( double ) );
01899 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01900 ptr += 2 * sizeof( double );
01901 if ( hasZValue )
01902 {
01903 ptr += sizeof( double );
01904 }
01905 newBufferPtr += 2 * sizeof( double );
01906 if ( hasZValue )
01907 {
01908 newBufferPtr += sizeof( double );
01909 }
01910 ++pointindex;
01911 if ( pointindex == vertexnr )
01912 {
01913 memcpy( newBufferPtr, &x, sizeof( double ) );
01914 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
01915 newBufferPtr += 2 * sizeof( double );
01916 if ( hasZValue )
01917 {
01918 newBufferPtr += sizeof( double );
01919 }
01920 success = true;
01921 }
01922 }
01923 break;
01924 }
01925 case QGis::WKBMultiLineString25D:
01926 hasZValue = true;
01927 case QGis::WKBMultiLineString:
01928 {
01929 int* nLines = ( int* )ptr;
01930 int* nPoints = 0;
01931 ptr += sizeof( int );
01932 memcpy( newBufferPtr, nLines, sizeof( int ) );
01933 newBufferPtr += sizeof( int );
01934 int pointindex = 0;
01935
01936 for ( int linenr = 0; linenr < *nLines; ++linenr )
01937 {
01938 memcpy( newBufferPtr, ptr, sizeof( int ) + 1 );
01939 ptr += ( sizeof( int ) + 1 );
01940 newBufferPtr += ( sizeof( int ) + 1 );
01941 nPoints = ( int* )ptr;
01942 int newNPoints;
01943 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )
01944 {
01945 newNPoints = ( *nPoints ) + 1;
01946 }
01947 else
01948 {
01949 newNPoints = *nPoints;
01950 }
01951 memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
01952 newBufferPtr += sizeof( int );
01953 ptr += sizeof( int );
01954
01955 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
01956 {
01957 memcpy( newBufferPtr, ptr, sizeof( double ) );
01958 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
01959 ptr += 2 * sizeof( double );
01960 newBufferPtr += 2 * sizeof( double );
01961 if ( hasZValue )
01962 {
01963 ptr += sizeof( double );
01964 newBufferPtr += sizeof( double );
01965 }
01966 ++pointindex;
01967 if ( pointindex == vertexnr )
01968 {
01969 memcpy( newBufferPtr, &x, sizeof( double ) );
01970 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
01971 newBufferPtr += 2 * sizeof( double );
01972 if ( hasZValue )
01973 {
01974 newBufferPtr += sizeof( double );
01975 }
01976 success = true;
01977 }
01978 }
01979 }
01980 break;
01981 }
01982 case QGis::WKBPolygon25D:
01983 hasZValue = true;
01984 case QGis::WKBPolygon:
01985 {
01986 int* nRings = ( int* )ptr;
01987 int* nPoints = 0;
01988 ptr += sizeof( int );
01989 memcpy( newBufferPtr, nRings, sizeof( int ) );
01990 newBufferPtr += sizeof( int );
01991 int pointindex = 0;
01992
01993 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
01994 {
01995 nPoints = ( int* )ptr;
01996 int newNPoints;
01997 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )
01998 {
01999 newNPoints = ( *nPoints ) + 1;
02000 }
02001 else
02002 {
02003 newNPoints = *nPoints;
02004 }
02005 memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
02006 newBufferPtr += sizeof( int );
02007 ptr += sizeof( int );
02008
02009 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02010 {
02011 memcpy( newBufferPtr, ptr, sizeof( double ) );
02012 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
02013 ptr += 2 * sizeof( double );
02014 newBufferPtr += 2 * sizeof( double );
02015 if ( hasZValue )
02016 {
02017 ptr += sizeof( double );
02018 newBufferPtr += sizeof( double );
02019 }
02020 ++pointindex;
02021 if ( pointindex == vertexnr )
02022 {
02023 memcpy( newBufferPtr, &x, sizeof( double ) );
02024 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
02025 newBufferPtr += 2 * sizeof( double );
02026 if ( hasZValue )
02027 {
02028 newBufferPtr += sizeof( double );
02029 }
02030 success = true;
02031 }
02032 }
02033 }
02034 break;
02035 }
02036 case QGis::WKBMultiPolygon25D:
02037 hasZValue = true;
02038 case QGis::WKBMultiPolygon:
02039 {
02040 int* nPolys = ( int* )ptr;
02041 int* nRings = 0;
02042 int* nPoints = 0;
02043 memcpy( newBufferPtr, nPolys, sizeof( int ) );
02044 ptr += sizeof( int );
02045 newBufferPtr += sizeof( int );
02046 int pointindex = 0;
02047
02048 for ( int polynr = 0; polynr < *nPolys; ++polynr )
02049 {
02050 memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) );
02051 ptr += ( 1 + sizeof( int ) );
02052 newBufferPtr += ( 1 + sizeof( int ) );
02053 nRings = ( int* )ptr;
02054 ptr += sizeof( int );
02055 memcpy( newBufferPtr, nRings, sizeof( int ) );
02056 newBufferPtr += sizeof( int );
02057
02058 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
02059 {
02060 nPoints = ( int* )ptr;
02061 int newNPoints;
02062 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )
02063 {
02064 newNPoints = ( *nPoints ) + 1;
02065 }
02066 else
02067 {
02068 newNPoints = *nPoints;
02069 }
02070 memcpy( newBufferPtr, &newNPoints, sizeof( double ) );
02071 newBufferPtr += sizeof( int );
02072 ptr += sizeof( int );
02073
02074 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02075 {
02076 memcpy( newBufferPtr, ptr, sizeof( double ) );
02077 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );
02078 ptr += 2 * sizeof( double );
02079 newBufferPtr += 2 * sizeof( double );
02080 if ( hasZValue )
02081 {
02082 ptr += sizeof( double );
02083 newBufferPtr += sizeof( double );
02084 }
02085 ++pointindex;
02086 if ( pointindex == vertexnr )
02087 {
02088 memcpy( newBufferPtr, &x, sizeof( double ) );
02089 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) );
02090 newBufferPtr += 2 * sizeof( double );
02091 if ( hasZValue )
02092 {
02093 newBufferPtr += sizeof( double );
02094 }
02095 success = true;
02096 }
02097 }
02098 }
02099
02100 }
02101 break;
02102 }
02103 case QGis::WKBNoGeometry:
02104 case QGis::WKBUnknown:
02105 break;
02106 }
02107
02108 if ( success )
02109 {
02110 delete [] mGeometry;
02111 mGeometry = newbuffer;
02112 mGeometrySize += 2 * sizeof( double );
02113 if ( hasZValue )
02114 {
02115 mGeometrySize += sizeof( double );
02116 }
02117 mDirtyGeos = true;
02118 return true;
02119 }
02120 else
02121 {
02122 delete newbuffer;
02123 return false;
02124 }
02125 }
02126
02127 QgsPoint QgsGeometry::vertexAt( int atVertex )
02128 {
02129 double x, y;
02130
02131 if ( mDirtyWkb )
02132 {
02133 exportGeosToWkb();
02134 }
02135
02136 if ( !mGeometry )
02137 {
02138 QgsDebugMsg( "WKB geometry not available!" );
02139 return QgsPoint( 0, 0 );
02140 }
02141
02142 QGis::WkbType wkbType;
02143 bool hasZValue = false;
02144 unsigned char* ptr;
02145
02146 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
02147 switch ( wkbType )
02148 {
02149 case QGis::WKBPoint25D:
02150 case QGis::WKBPoint:
02151 {
02152 if ( atVertex == 0 )
02153 {
02154 ptr = mGeometry + 1 + sizeof( int );
02155 memcpy( &x, ptr, sizeof( double ) );
02156 ptr += sizeof( double );
02157 memcpy( &y, ptr, sizeof( double ) );
02158 return QgsPoint( x, y );
02159 }
02160 else
02161 {
02162 return QgsPoint( 0, 0 );
02163 }
02164 }
02165 case QGis::WKBLineString25D:
02166 hasZValue = true;
02167 case QGis::WKBLineString:
02168 {
02169 int *nPoints;
02170
02171 ptr = mGeometry + 1 + sizeof( int );
02172 nPoints = ( int * ) ptr;
02173
02174
02175 if ( 0 > atVertex || *nPoints <= atVertex )
02176 {
02177 return QgsPoint( 0, 0 );
02178 }
02179
02180
02181 if ( hasZValue )
02182 {
02183 ptr = mGeometry + 9 + ( atVertex * 3 * sizeof( double ) );
02184 }
02185 else
02186 {
02187 ptr = mGeometry + 9 + ( atVertex * 2 * sizeof( double ) );
02188 }
02189 memcpy( &x, ptr, sizeof( double ) );
02190 ptr += sizeof( double );
02191 memcpy( &y, ptr, sizeof( double ) );
02192 return QgsPoint( x, y );
02193 }
02194 case QGis::WKBPolygon25D:
02195 hasZValue = true;
02196 case QGis::WKBPolygon:
02197 {
02198 int *nRings;
02199 int *nPoints = 0;
02200 ptr = mGeometry + 1 + sizeof( int );
02201 nRings = ( int* )ptr;
02202 ptr += sizeof( int );
02203 int pointindex = 0;
02204 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
02205 {
02206 nPoints = ( int* )ptr;
02207 ptr += sizeof( int );
02208 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02209 {
02210 if ( pointindex == atVertex )
02211 {
02212 memcpy( &x, ptr, sizeof( double ) );
02213 ptr += sizeof( double );
02214 memcpy( &y, ptr, sizeof( double ) );
02215 return QgsPoint( x, y );
02216 }
02217 ptr += 2 * sizeof( double );
02218 if ( hasZValue )
02219 {
02220 ptr += sizeof( double );
02221 }
02222 ++pointindex;
02223 }
02224 }
02225 return QgsPoint( 0, 0 );
02226 }
02227 case QGis::WKBMultiPoint25D:
02228 hasZValue = true;
02229 case QGis::WKBMultiPoint:
02230 {
02231 ptr = mGeometry + 1 + sizeof( int );
02232 int* nPoints = ( int* )ptr;
02233 if ( atVertex < 0 || atVertex >= *nPoints )
02234 {
02235 return QgsPoint( 0, 0 );
02236 }
02237 if ( hasZValue )
02238 {
02239 ptr += atVertex * ( 3 * sizeof( double ) + 1 + sizeof( int ) );
02240 }
02241 else
02242 {
02243 ptr += atVertex * ( 2 * sizeof( double ) + 1 + sizeof( int ) );
02244 }
02245 ptr += 1 + sizeof( int );
02246 memcpy( &x, ptr, sizeof( double ) );
02247 ptr += sizeof( double );
02248 memcpy( &y, ptr, sizeof( double ) );
02249 return QgsPoint( x, y );
02250 }
02251 case QGis::WKBMultiLineString25D:
02252 hasZValue = true;
02253 case QGis::WKBMultiLineString:
02254 {
02255 ptr = mGeometry + 1 + sizeof( int );
02256 int* nLines = ( int* )ptr;
02257 int* nPoints = 0;
02258 int pointindex = 0;
02259 ptr += sizeof( int );
02260 for ( int linenr = 0; linenr < *nLines; ++linenr )
02261 {
02262 ptr += sizeof( int ) + 1;
02263 nPoints = ( int* )ptr;
02264 ptr += sizeof( int );
02265 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02266 {
02267 if ( pointindex == atVertex )
02268 {
02269 memcpy( &x, ptr, sizeof( double ) );
02270 ptr += sizeof( double );
02271 memcpy( &y, ptr, sizeof( double ) );
02272 return QgsPoint( x, y );
02273 }
02274 ptr += 2 * sizeof( double );
02275 if ( hasZValue )
02276 {
02277 ptr += sizeof( double );
02278 }
02279 ++pointindex;
02280 }
02281 }
02282 return QgsPoint( 0, 0 );
02283 }
02284 case QGis::WKBMultiPolygon25D:
02285 hasZValue = true;
02286 case QGis::WKBMultiPolygon:
02287 {
02288 ptr = mGeometry + 1 + sizeof( int );
02289 int* nRings = 0;
02290 int* nPoints = 0;
02291 int pointindex = 0;
02292 int* nPolygons = ( int* )ptr;
02293 ptr += sizeof( int );
02294 for ( int polynr = 0; polynr < *nPolygons; ++polynr )
02295 {
02296 ptr += ( 1 + sizeof( int ) );
02297 nRings = ( int* )ptr;
02298 ptr += sizeof( int );
02299 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
02300 {
02301 nPoints = ( int* )ptr;
02302 ptr += sizeof( int );
02303 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02304 {
02305 if ( pointindex == atVertex )
02306 {
02307 memcpy( &x, ptr, sizeof( double ) );
02308 ptr += sizeof( double );
02309 memcpy( &y, ptr, sizeof( double ) );
02310 return QgsPoint( x, y );
02311 }
02312 ++pointindex;
02313 ptr += 2 * sizeof( double );
02314 if ( hasZValue )
02315 {
02316 ptr += sizeof( double );
02317 }
02318 }
02319 }
02320 }
02321 return QgsPoint( 0, 0 );
02322 }
02323 default:
02324 QgsDebugMsg( "error: mGeometry type not recognized" );
02325 return QgsPoint( 0, 0 );
02326 }
02327 }
02328
02329
02330 double QgsGeometry::sqrDistToVertexAt( QgsPoint& point, int atVertex )
02331 {
02332 QgsPoint pnt = vertexAt( atVertex );
02333 if ( pnt != QgsPoint( 0, 0 ) )
02334 {
02335 QgsDebugMsg( "Exiting with distance to " + pnt.toString() );
02336 return point.sqrDist( pnt );
02337 }
02338 else
02339 {
02340 QgsDebugMsg( "Exiting with std::numeric_limits<double>::max()." );
02341
02342 return std::numeric_limits<double>::max();
02343 }
02344 }
02345
02346
02347 double QgsGeometry::closestVertexWithContext( const QgsPoint& point, int& atVertex )
02348 {
02349 double sqrDist = std::numeric_limits<double>::max();
02350
02351 try
02352 {
02353
02354 int closestVertexIndex = 0;
02355
02356
02357 if ( !exportWkbToGeos() )
02358 return -1;
02359
02360 const GEOSGeometry *g = GEOSGetExteriorRing( mGeos );
02361 if ( !g )
02362 return -1;
02363
02364 const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( g );
02365
02366 unsigned int n;
02367 GEOSCoordSeq_getSize( sequence, &n );
02368
02369 for ( unsigned int i = 0; i < n; i++ )
02370 {
02371 double x, y;
02372 GEOSCoordSeq_getX( sequence, i, &x );
02373 GEOSCoordSeq_getY( sequence, i, &y );
02374
02375 double testDist = point.sqrDist( x, y );
02376 if ( testDist < sqrDist )
02377 {
02378 closestVertexIndex = i;
02379 sqrDist = testDist;
02380 }
02381 }
02382
02383 atVertex = closestVertexIndex;
02384 }
02385 catch ( GEOSException &e )
02386 {
02387 Q_UNUSED( e );
02388 return -1;
02389 }
02390
02391 return sqrDist;
02392 }
02393
02394
02395 double QgsGeometry::closestSegmentWithContext(
02396 const QgsPoint& point,
02397 QgsPoint& minDistPoint,
02398 int& beforeVertex )
02399 {
02400 QgsPoint distPoint;
02401
02402 QGis::WkbType wkbType;
02403 bool hasZValue = false;
02404 double *thisx = NULL;
02405 double *thisy = NULL;
02406 double *prevx = NULL;
02407 double *prevy = NULL;
02408 double testdist;
02409 int closestSegmentIndex = 0;
02410
02411
02412 double sqrDist = std::numeric_limits<double>::max();
02413
02414
02415 if ( mDirtyWkb )
02416 {
02417 exportGeosToWkb();
02418 }
02419
02420 if ( !mGeometry )
02421 {
02422 QgsDebugMsg( "WKB geometry not available!" );
02423 return -1;
02424 }
02425
02426 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) );
02427
02428 switch ( wkbType )
02429 {
02430 case QGis::WKBPoint25D:
02431 case QGis::WKBPoint:
02432 case QGis::WKBMultiPoint25D:
02433 case QGis::WKBMultiPoint:
02434 {
02435
02436 return -1;
02437 }
02438 case QGis::WKBLineString25D:
02439 hasZValue = true;
02440 case QGis::WKBLineString:
02441 {
02442 unsigned char* ptr = mGeometry + 1 + sizeof( int );
02443 int* npoints = ( int* ) ptr;
02444 ptr += sizeof( int );
02445 for ( int index = 0; index < *npoints; ++index )
02446 {
02447 if ( index > 0 )
02448 {
02449 prevx = thisx;
02450 prevy = thisy;
02451 }
02452 thisx = ( double* ) ptr;
02453 ptr += sizeof( double );
02454 thisy = ( double* ) ptr;
02455
02456 if ( index > 0 )
02457 {
02458 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
02459 {
02460 closestSegmentIndex = index;
02461 sqrDist = testdist;
02462 minDistPoint = distPoint;
02463 }
02464 }
02465 ptr += sizeof( double );
02466 if ( hasZValue )
02467 {
02468 ptr += sizeof( double );
02469 }
02470 }
02471 beforeVertex = closestSegmentIndex;
02472 break;
02473 }
02474 case QGis::WKBMultiLineString25D:
02475 hasZValue = true;
02476 case QGis::WKBMultiLineString:
02477 {
02478 unsigned char* ptr = mGeometry + 1 + sizeof( int );
02479 int* nLines = ( int* )ptr;
02480 ptr += sizeof( int );
02481 int* nPoints = 0;
02482 int pointindex = 0;
02483 for ( int linenr = 0; linenr < *nLines; ++linenr )
02484 {
02485 ptr += sizeof( int ) + 1;
02486 nPoints = ( int* )ptr;
02487 ptr += sizeof( int );
02488 prevx = 0;
02489 prevy = 0;
02490 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02491 {
02492 thisx = ( double* ) ptr;
02493 ptr += sizeof( double );
02494 thisy = ( double* ) ptr;
02495 ptr += sizeof( double );
02496 if ( hasZValue )
02497 {
02498 ptr += sizeof( double );
02499 }
02500 if ( prevx && prevy )
02501 {
02502 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
02503 {
02504 closestSegmentIndex = pointindex;
02505 sqrDist = testdist;
02506 minDistPoint = distPoint;
02507 }
02508 }
02509 prevx = thisx;
02510 prevy = thisy;
02511 ++pointindex;
02512 }
02513 }
02514 beforeVertex = closestSegmentIndex;
02515 break;
02516 }
02517 case QGis::WKBPolygon25D:
02518 hasZValue = true;
02519 case QGis::WKBPolygon:
02520 {
02521 int index = 0;
02522 unsigned char* ptr = mGeometry + 1 + sizeof( int );
02523 int* nrings = ( int* )ptr;
02524 int* npoints = 0;
02525 ptr += sizeof( int );
02526 for ( int ringnr = 0; ringnr < *nrings; ++ringnr )
02527 {
02528 npoints = ( int* )ptr;
02529 ptr += sizeof( int );
02530 prevx = 0;
02531 prevy = 0;
02532 for ( int pointnr = 0; pointnr < *npoints; ++pointnr )
02533 {
02534 thisx = ( double* )ptr;
02535 ptr += sizeof( double );
02536 thisy = ( double* )ptr;
02537 ptr += sizeof( double );
02538 if ( hasZValue )
02539 {
02540 ptr += sizeof( double );
02541 }
02542 if ( prevx && prevy )
02543 {
02544 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
02545 {
02546 closestSegmentIndex = index;
02547 sqrDist = testdist;
02548 minDistPoint = distPoint;
02549 }
02550 }
02551 prevx = thisx;
02552 prevy = thisy;
02553 ++index;
02554 }
02555 }
02556 beforeVertex = closestSegmentIndex;
02557 break;
02558 }
02559 case QGis::WKBMultiPolygon25D:
02560 hasZValue = true;
02561 case QGis::WKBMultiPolygon:
02562 {
02563 unsigned char* ptr = mGeometry + 1 + sizeof( int );
02564 int* nRings = 0;
02565 int* nPoints = 0;
02566 int pointindex = 0;
02567 int* nPolygons = ( int* )ptr;
02568 ptr += sizeof( int );
02569 for ( int polynr = 0; polynr < *nPolygons; ++polynr )
02570 {
02571 ptr += ( 1 + sizeof( int ) );
02572 nRings = ( int* )ptr;
02573 ptr += sizeof( int );
02574 for ( int ringnr = 0; ringnr < *nRings; ++ringnr )
02575 {
02576 nPoints = ( int* )ptr;
02577 ptr += sizeof( int );
02578 prevx = 0;
02579 prevy = 0;
02580 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr )
02581 {
02582 thisx = ( double* )ptr;
02583 ptr += sizeof( double );
02584 thisy = ( double* )ptr;
02585 ptr += sizeof( double );
02586 if ( hasZValue )
02587 {
02588 ptr += sizeof( double );
02589 }
02590 if ( prevx && prevy )
02591 {
02592 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint ) ) < sqrDist )
02593 {
02594 closestSegmentIndex = pointindex;
02595 sqrDist = testdist;
02596 minDistPoint = distPoint;
02597 }
02598 }
02599 prevx = thisx;
02600 prevy = thisy;
02601 ++pointindex;
02602 }
02603 }
02604 }
02605 beforeVertex = closestSegmentIndex;
02606 break;
02607 }
02608 case QGis::WKBUnknown:
02609 default:
02610 return -1;
02611 break;
02612 }
02613
02614
02615 QgsDebugMsg( "Exiting with nearest point " + point.toString() +
02616 ", dist: " + QString::number( sqrDist ) + "." );
02617
02618 return sqrDist;
02619 }
02620
02621 int QgsGeometry::addRing( const QList<QgsPoint>& ring )
02622 {
02623
02624 if ( type() != QGis::Polygon )
02625 return 1;
02626
02627
02628 if ( ring.size() < 4 )
02629 return 3;
02630
02631
02632 if ( ring.first() != ring.last() )
02633 return 2;
02634
02635
02636 if ( !mGeos || mDirtyGeos )
02637 if ( !exportWkbToGeos() )
02638 return 6;
02639
02640 int type = GEOSGeomTypeId( mGeos );
02641
02642
02643 QVector<const GEOSGeometry*> polygonList;
02644
02645 if ( wkbType() == QGis::WKBPolygon )
02646 {
02647 if ( type != GEOS_POLYGON )
02648 return 1;
02649
02650 polygonList << mGeos;
02651 }
02652 else if ( wkbType() == QGis::WKBMultiPolygon )
02653 {
02654 if ( type != GEOS_MULTIPOLYGON )
02655 return 1;
02656
02657 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); ++i )
02658 polygonList << GEOSGetGeometryN( mGeos, i );
02659 }
02660
02661
02662 GEOSGeometry *newRing = 0;
02663 GEOSGeometry *newRingPolygon = 0;
02664
02665 try
02666 {
02667 newRing = createGeosLinearRing( ring.toVector() );
02668 if ( !GEOSisValid( newRing ) )
02669 {
02670 throwGEOSException( "ring is invalid" );
02671 }
02672
02673 newRingPolygon = createGeosPolygon( newRing );
02674 if ( !GEOSisValid( newRingPolygon ) )
02675 {
02676 throwGEOSException( "ring is invalid" );
02677 }
02678 }
02679 catch ( GEOSException &e )
02680 {
02681 Q_UNUSED( e );
02682 if ( newRingPolygon )
02683 GEOSGeom_destroy( newRingPolygon );
02684 else if ( newRing )
02685 GEOSGeom_destroy( newRing );
02686
02687 return 3;
02688 }
02689
02690 QVector<GEOSGeometry*> rings;
02691
02692 int i;
02693 for ( i = 0; i < polygonList.size(); i++ )
02694 {
02695 for ( int j = 0; j < rings.size(); j++ )
02696 GEOSGeom_destroy( rings[j] );
02697 rings.clear();
02698
02699 GEOSGeometry *shellRing = 0;
02700 GEOSGeometry *shell = 0;
02701 try
02702 {
02703 shellRing = GEOSGeom_clone( GEOSGetExteriorRing( polygonList[i] ) );
02704 shell = createGeosPolygon( shellRing );
02705
02706 if ( !GEOSWithin( newRingPolygon, shell ) )
02707 {
02708 GEOSGeom_destroy( shell );
02709 continue;
02710 }
02711 }
02712 catch ( GEOSException &e )
02713 {
02714 Q_UNUSED( e );
02715 if ( shell )
02716 GEOSGeom_destroy( shell );
02717 else if ( shellRing )
02718 GEOSGeom_destroy( shellRing );
02719
02720 GEOSGeom_destroy( newRingPolygon );
02721
02722 return 4;
02723 }
02724
02725
02726 rings << GEOSGeom_clone( shellRing );
02727
02728 GEOSGeom_destroy( shell );
02729
02730
02731 int n = GEOSGetNumInteriorRings( polygonList[i] );
02732
02733 int j;
02734 for ( j = 0; j < n; j++ )
02735 {
02736 GEOSGeometry *holeRing = 0;
02737 GEOSGeometry *hole = 0;
02738 try
02739 {
02740 holeRing = GEOSGeom_clone( GEOSGetInteriorRingN( polygonList[i], j ) );
02741 hole = createGeosPolygon( holeRing );
02742
02743 if ( !GEOSDisjoint( hole, newRingPolygon ) )
02744 {
02745 GEOSGeom_destroy( hole );
02746 break;
02747 }
02748 }
02749 catch ( GEOSException &e )
02750 {
02751 Q_UNUSED( e );
02752 if ( hole )
02753 GEOSGeom_destroy( hole );
02754 else if ( holeRing )
02755 GEOSGeom_destroy( holeRing );
02756
02757 break;
02758 }
02759
02760 rings << GEOSGeom_clone( holeRing );
02761 GEOSGeom_destroy( hole );
02762 }
02763
02764 if ( j == n )
02765
02766 break;
02767 }
02768
02769 if ( i == polygonList.size() )
02770 {
02771
02772 for ( int j = 0; j < rings.size(); j++ )
02773 GEOSGeom_destroy( rings[j] );
02774 rings.clear();
02775
02776 GEOSGeom_destroy( newRingPolygon );
02777
02778
02779 return 5;
02780 }
02781
02782 rings << GEOSGeom_clone( newRing );
02783 GEOSGeom_destroy( newRingPolygon );
02784
02785 GEOSGeometry *newPolygon = createGeosPolygon( rings );
02786
02787 if ( wkbType() == QGis::WKBPolygon )
02788 {
02789 GEOSGeom_destroy( mGeos );
02790 mGeos = newPolygon;
02791 }
02792 else if ( wkbType() == QGis::WKBMultiPolygon )
02793 {
02794 QVector<GEOSGeometry*> newPolygons;
02795
02796 for ( int j = 0; j < polygonList.size(); j++ )
02797 {
02798 newPolygons << ( i == j ? newPolygon : GEOSGeom_clone( polygonList[j] ) );
02799 }
02800
02801 GEOSGeom_destroy( mGeos );
02802 mGeos = createGeosCollection( GEOS_MULTIPOLYGON, newPolygons );
02803 }
02804
02805 mDirtyWkb = true;
02806 mDirtyGeos = false;
02807 return 0;
02808 }
02809
02810 int QgsGeometry::addIsland( const QList<QgsPoint>& ring )
02811 {
02812
02813 if ( ring.size() < 4 )
02814 {
02815 return 2;
02816 }
02817
02818
02819 if ( ring.first() != ring.last() )
02820 {
02821 return 2;
02822 }
02823
02824 if ( wkbType() == QGis::WKBPolygon || wkbType() == QGis::WKBPolygon25D )
02825 {
02826 if ( !convertToMultiType() )
02827 {
02828 return 1;
02829 }
02830 }
02831
02832
02833 if ( wkbType() != QGis::WKBMultiPolygon && wkbType() != QGis::WKBMultiPolygon25D )
02834 {
02835 return 1;
02836 }
02837
02838
02839 if ( !mGeos || mDirtyGeos )
02840 {
02841 if ( !exportWkbToGeos() )
02842 return 4;
02843 }
02844
02845 if ( GEOSGeomTypeId( mGeos ) != GEOS_MULTIPOLYGON )
02846 return 1;
02847
02848
02849
02850
02851 GEOSGeometry *newRing = 0;
02852 GEOSGeometry *newPolygon = 0;
02853
02854 try
02855 {
02856 newRing = createGeosLinearRing( ring.toVector() );
02857 if ( !GEOSisValid( newRing ) )
02858 throw GEOSException( "ring invalid" );
02859 newPolygon = createGeosPolygon( newRing );
02860 if ( !GEOSisValid( newPolygon ) )
02861 throw GEOSException( "polygon invalid" );
02862 }
02863 catch ( GEOSException &e )
02864 {
02865 Q_UNUSED( e );
02866 if ( newPolygon )
02867 GEOSGeom_destroy( newPolygon );
02868 else if ( newRing )
02869 GEOSGeom_destroy( newRing );
02870 return 2;
02871 }
02872
02873 QVector<GEOSGeometry*> polygons;
02874
02875
02876 int n = GEOSGetNumGeometries( mGeos );
02877 int i;
02878 for ( i = 0; i < n; ++i )
02879 {
02880 const GEOSGeometry *polygonN = GEOSGetGeometryN( mGeos, i );
02881
02882 if ( !GEOSDisjoint( polygonN, newPolygon ) )
02883
02884 break;
02885
02886 polygons << GEOSGeom_clone( polygonN );
02887 }
02888
02889 if ( i < n )
02890 {
02891
02892 for ( int i = 0; i < polygons.size(); i++ )
02893 GEOSGeom_destroy( polygons[i] );
02894 return 3;
02895 }
02896
02897 polygons << newPolygon;
02898
02899 GEOSGeom_destroy( mGeos );
02900 mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
02901 mDirtyWkb = true;
02902 mDirtyGeos = false;
02903 return 0;
02904 }
02905
02906 int QgsGeometry::translate( double dx, double dy )
02907 {
02908 if ( mDirtyWkb )
02909 {
02910 exportGeosToWkb();
02911 }
02912
02913 if ( !mGeometry )
02914 {
02915 QgsDebugMsg( "WKB geometry not available!" );
02916 return 1;
02917 }
02918
02919 QGis::WkbType wkbType;
02920 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
02921 bool hasZValue = false;
02922 int wkbPosition = 5;
02923
02924 switch ( wkbType )
02925 {
02926 case QGis::WKBPoint25D:
02927 case QGis::WKBPoint:
02928 {
02929 translateVertex( wkbPosition, dx, dy, hasZValue );
02930 }
02931 break;
02932
02933 case QGis::WKBLineString25D:
02934 hasZValue = true;
02935 case QGis::WKBLineString:
02936 {
02937 int* npoints = ( int* )( &mGeometry[wkbPosition] );
02938 wkbPosition += sizeof( int );
02939 for ( int index = 0; index < *npoints; ++index )
02940 {
02941 translateVertex( wkbPosition, dx, dy, hasZValue );
02942 }
02943 break;
02944 }
02945
02946 case QGis::WKBPolygon25D:
02947 hasZValue = true;
02948 case QGis::WKBPolygon:
02949 {
02950 int* nrings = ( int* )( &( mGeometry[wkbPosition] ) );
02951 wkbPosition += sizeof( int );
02952 int* npoints;
02953
02954 for ( int index = 0; index < *nrings; ++index )
02955 {
02956 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
02957 wkbPosition += sizeof( int );
02958 for ( int index2 = 0; index2 < *npoints; ++index2 )
02959 {
02960 translateVertex( wkbPosition, dx, dy, hasZValue );
02961 }
02962 }
02963 break;
02964 }
02965
02966 case QGis::WKBMultiPoint25D:
02967 hasZValue = true;
02968 case QGis::WKBMultiPoint:
02969 {
02970 int* npoints = ( int* )( &( mGeometry[wkbPosition] ) );
02971 wkbPosition += sizeof( int );
02972 for ( int index = 0; index < *npoints; ++index )
02973 {
02974 wkbPosition += ( sizeof( int ) + 1 );
02975 translateVertex( wkbPosition, dx, dy, hasZValue );
02976 }
02977 break;
02978 }
02979
02980 case QGis::WKBMultiLineString25D:
02981 hasZValue = true;
02982 case QGis::WKBMultiLineString:
02983 {
02984 int* nlines = ( int* )( &( mGeometry[wkbPosition] ) );
02985 int* npoints = 0;
02986 wkbPosition += sizeof( int );
02987 for ( int index = 0; index < *nlines; ++index )
02988 {
02989 wkbPosition += ( sizeof( int ) + 1 );
02990 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
02991 wkbPosition += sizeof( int );
02992 for ( int index2 = 0; index2 < *npoints; ++index2 )
02993 {
02994 translateVertex( wkbPosition, dx, dy, hasZValue );
02995 }
02996 }
02997 break;
02998 }
02999
03000 case QGis::WKBMultiPolygon25D:
03001 hasZValue = true;
03002 case QGis::WKBMultiPolygon:
03003 {
03004 int* npolys = ( int* )( &( mGeometry[wkbPosition] ) );
03005 int* nrings;
03006 int* npoints;
03007 wkbPosition += sizeof( int );
03008 for ( int index = 0; index < *npolys; ++index )
03009 {
03010 wkbPosition += ( 1 + sizeof( int ) );
03011 nrings = ( int* )( &( mGeometry[wkbPosition] ) );
03012 wkbPosition += sizeof( int );
03013 for ( int index2 = 0; index2 < *nrings; ++index2 )
03014 {
03015 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
03016 wkbPosition += sizeof( int );
03017 for ( int index3 = 0; index3 < *npoints; ++index3 )
03018 {
03019 translateVertex( wkbPosition, dx, dy, hasZValue );
03020 }
03021 }
03022 }
03023 }
03024
03025 default:
03026 break;
03027 }
03028 mDirtyGeos = true;
03029 return 0;
03030 }
03031
03032 int QgsGeometry::transform( const QgsCoordinateTransform& ct )
03033 {
03034 if ( mDirtyWkb )
03035 {
03036 exportGeosToWkb();
03037 }
03038
03039 if ( !mGeometry )
03040 {
03041 QgsDebugMsg( "WKB geometry not available!" );
03042 return 1;
03043 }
03044
03045 QGis::WkbType wkbType;
03046 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
03047 bool hasZValue = false;
03048 int wkbPosition = 5;
03049
03050 switch ( wkbType )
03051 {
03052 case QGis::WKBPoint25D:
03053 case QGis::WKBPoint:
03054 {
03055 transformVertex( wkbPosition, ct, hasZValue );
03056 }
03057 break;
03058
03059 case QGis::WKBLineString25D:
03060 hasZValue = true;
03061 case QGis::WKBLineString:
03062 {
03063 int* npoints = ( int* )( &mGeometry[wkbPosition] );
03064 wkbPosition += sizeof( int );
03065 for ( int index = 0; index < *npoints; ++index )
03066 {
03067 transformVertex( wkbPosition, ct, hasZValue );
03068 }
03069 break;
03070 }
03071
03072 case QGis::WKBPolygon25D:
03073 hasZValue = true;
03074 case QGis::WKBPolygon:
03075 {
03076 int* nrings = ( int* )( &( mGeometry[wkbPosition] ) );
03077 wkbPosition += sizeof( int );
03078 int* npoints;
03079
03080 for ( int index = 0; index < *nrings; ++index )
03081 {
03082 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
03083 wkbPosition += sizeof( int );
03084 for ( int index2 = 0; index2 < *npoints; ++index2 )
03085 {
03086 transformVertex( wkbPosition, ct, hasZValue );
03087 }
03088 }
03089 break;
03090 }
03091
03092 case QGis::WKBMultiPoint25D:
03093 hasZValue = true;
03094 case QGis::WKBMultiPoint:
03095 {
03096 int* npoints = ( int* )( &( mGeometry[wkbPosition] ) );
03097 wkbPosition += sizeof( int );
03098 for ( int index = 0; index < *npoints; ++index )
03099 {
03100 wkbPosition += ( sizeof( int ) + 1 );
03101 transformVertex( wkbPosition, ct, hasZValue );
03102 }
03103 break;
03104 }
03105
03106 case QGis::WKBMultiLineString25D:
03107 hasZValue = true;
03108 case QGis::WKBMultiLineString:
03109 {
03110 int* nlines = ( int* )( &( mGeometry[wkbPosition] ) );
03111 int* npoints = 0;
03112 wkbPosition += sizeof( int );
03113 for ( int index = 0; index < *nlines; ++index )
03114 {
03115 wkbPosition += ( sizeof( int ) + 1 );
03116 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
03117 wkbPosition += sizeof( int );
03118 for ( int index2 = 0; index2 < *npoints; ++index2 )
03119 {
03120 transformVertex( wkbPosition, ct, hasZValue );
03121 }
03122 }
03123 break;
03124 }
03125
03126 case QGis::WKBMultiPolygon25D:
03127 hasZValue = true;
03128 case QGis::WKBMultiPolygon:
03129 {
03130 int* npolys = ( int* )( &( mGeometry[wkbPosition] ) );
03131 int* nrings;
03132 int* npoints;
03133 wkbPosition += sizeof( int );
03134 for ( int index = 0; index < *npolys; ++index )
03135 {
03136 wkbPosition += ( 1 + sizeof( int ) );
03137 nrings = ( int* )( &( mGeometry[wkbPosition] ) );
03138 wkbPosition += sizeof( int );
03139 for ( int index2 = 0; index2 < *nrings; ++index2 )
03140 {
03141 npoints = ( int* )( &( mGeometry[wkbPosition] ) );
03142 wkbPosition += sizeof( int );
03143 for ( int index3 = 0; index3 < *npoints; ++index3 )
03144 {
03145 transformVertex( wkbPosition, ct, hasZValue );
03146 }
03147 }
03148 }
03149 }
03150
03151 default:
03152 break;
03153 }
03154 mDirtyGeos = true;
03155 return 0;
03156 }
03157
03158 int QgsGeometry::splitGeometry( const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint> &topologyTestPoints )
03159 {
03160 int returnCode = 0;
03161
03162
03163 if ( type() == QGis::Point )
03164 {
03165 return 1;
03166 }
03167
03168
03169 if ( mDirtyWkb )
03170 {
03171 exportGeosToWkb();
03172 }
03173 if ( !mGeos || mDirtyGeos )
03174 {
03175 if ( !exportWkbToGeos() )
03176 return 1;
03177 }
03178
03179
03180 if ( splitLine.size() < 2 )
03181 {
03182 return 1;
03183 }
03184
03185 newGeometries.clear();
03186
03187 try
03188 {
03189 GEOSGeometry *splitLineGeos = createGeosLineString( splitLine.toVector() );
03190 if ( !GEOSisValid( splitLineGeos ) || !GEOSisSimple( splitLineGeos ) )
03191 {
03192 GEOSGeom_destroy( splitLineGeos );
03193 return 1;
03194 }
03195
03196 if ( topological )
03197 {
03198
03199 if ( topologicalTestPointsSplit( splitLineGeos, topologyTestPoints ) != 0 )
03200 {
03201 return 1;
03202 }
03203 }
03204
03205
03206 if ( type() == QGis::Line )
03207 {
03208 returnCode = splitLinearGeometry( splitLineGeos, newGeometries );
03209 GEOSGeom_destroy( splitLineGeos );
03210 }
03211 else if ( type() == QGis::Polygon )
03212 {
03213 returnCode = splitPolygonGeometry( splitLineGeos, newGeometries );
03214 GEOSGeom_destroy( splitLineGeos );
03215 }
03216 else
03217 {
03218 return 1;
03219 }
03220 }
03221 CATCH_GEOS( 2 )
03222
03223 return returnCode;
03224 }
03225
03227 int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine )
03228 {
03229 if ( reshapeWithLine.size() < 2 )
03230 {
03231 return 1;
03232 }
03233 if ( type() == QGis::Point )
03234 {
03235 return 1;
03236 }
03237
03238 GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() );
03239
03240
03241 if ( !mGeos || mDirtyGeos )
03242 {
03243 exportWkbToGeos();
03244 }
03245
03246
03247 int numGeoms = GEOSGetNumGeometries( mGeos );
03248 if ( numGeoms == -1 )
03249 {
03250 return 1;
03251 }
03252
03253 bool isMultiGeom = false;
03254 int geosTypeId = GEOSGeomTypeId( mGeos );
03255 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
03256 {
03257 isMultiGeom = true;
03258 }
03259
03260 bool isLine = ( type() == QGis::Line );
03261
03262
03263 if ( !isMultiGeom )
03264 {
03265 GEOSGeometry* reshapedGeometry;
03266 if ( isLine )
03267 {
03268 reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos );
03269 }
03270 else
03271 {
03272 reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos );
03273 }
03274
03275 GEOSGeom_destroy( reshapeLineGeos );
03276 if ( reshapedGeometry )
03277 {
03278 GEOSGeom_destroy( mGeos );
03279 mGeos = reshapedGeometry;
03280 mDirtyWkb = true;
03281 return 0;
03282 }
03283 else
03284 {
03285 return 1;
03286 }
03287 }
03288 else
03289 {
03290
03291 bool reshapeTookPlace = false;
03292
03293 GEOSGeometry* currentReshapeGeometry = 0;
03294 GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms];
03295
03296 for ( int i = 0; i < numGeoms; ++i )
03297 {
03298 if ( isLine )
03299 {
03300 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos );
03301 }
03302 else
03303 {
03304 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos );
03305 }
03306
03307 if ( currentReshapeGeometry )
03308 {
03309 newGeoms[i] = currentReshapeGeometry;
03310 reshapeTookPlace = true;
03311 }
03312 else
03313 {
03314 newGeoms[i] = GEOSGeom_clone( GEOSGetGeometryN( mGeos, i ) );
03315 }
03316 }
03317 GEOSGeom_destroy( reshapeLineGeos );
03318
03319 GEOSGeometry* newMultiGeom = 0;
03320 if ( isLine )
03321 {
03322 newMultiGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, newGeoms, numGeoms );
03323 }
03324 else
03325 {
03326 newMultiGeom = GEOSGeom_createCollection( GEOS_MULTIPOLYGON, newGeoms, numGeoms );
03327 }
03328
03329 delete[] newGeoms;
03330 if ( ! newMultiGeom )
03331 {
03332 return 3;
03333 }
03334
03335 if ( reshapeTookPlace )
03336 {
03337 GEOSGeom_destroy( mGeos );
03338 mGeos = newMultiGeom;
03339 mDirtyWkb = true;
03340 return 0;
03341 }
03342 else
03343 {
03344 GEOSGeom_destroy( newMultiGeom );
03345 return 1;
03346 }
03347 }
03348 }
03349
03350 int QgsGeometry::makeDifference( QgsGeometry* other )
03351 {
03352
03353 if ( !mGeos || mDirtyGeos )
03354 {
03355 exportWkbToGeos();
03356 }
03357
03358 if ( !mGeos )
03359 {
03360 return 1;
03361 }
03362
03363 if ( !GEOSisValid( mGeos ) )
03364 {
03365 return 2;
03366 }
03367
03368 if ( !GEOSisSimple( mGeos ) )
03369 {
03370 return 3;
03371 }
03372
03373
03374 if ( !other->mGeos || other->mDirtyGeos )
03375 {
03376 other->exportWkbToGeos();
03377 }
03378
03379 if ( !other->mGeos )
03380 {
03381 return 4;
03382 }
03383
03384
03385 try
03386 {
03387 if ( GEOSIntersects( mGeos, other->mGeos ) )
03388 {
03389
03390 bool multiType = isMultipart();
03391
03392 mGeos = GEOSDifference( mGeos, other->mGeos );
03393 mDirtyWkb = true;
03394
03395 if ( multiType && !isMultipart() )
03396 {
03397 convertToMultiType();
03398 exportWkbToGeos();
03399 }
03400 }
03401 else
03402 {
03403 return 0;
03404 }
03405 }
03406 CATCH_GEOS( 5 )
03407
03408 if ( !mGeos )
03409 {
03410 mDirtyGeos = true;
03411 return 6;
03412 }
03413
03414 return 0;
03415 }
03416
03417 QgsRectangle QgsGeometry::boundingBox()
03418 {
03419 double xmin = std::numeric_limits<double>::max();
03420 double ymin = std::numeric_limits<double>::max();
03421 double xmax = -std::numeric_limits<double>::max();
03422 double ymax = -std::numeric_limits<double>::max();
03423
03424 double *x;
03425 double *y;
03426 int *nPoints;
03427 int *numRings;
03428 int *numPolygons;
03429 int numLineStrings;
03430 int idx, jdx, kdx;
03431 unsigned char *ptr;
03432 char lsb;
03433 QgsPoint pt;
03434 QGis::WkbType wkbType;
03435 bool hasZValue = false;
03436
03437
03438 if ( mDirtyWkb )
03439 {
03440 exportGeosToWkb();
03441 }
03442
03443 if ( !mGeometry )
03444 {
03445 QgsDebugMsg( "WKB geometry not available!" );
03446 return QgsRectangle( 0, 0, 0, 0 );
03447 }
03448
03449
03450 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
03451 switch ( wkbType )
03452 {
03453 case QGis::WKBPoint25D:
03454 case QGis::WKBPoint:
03455 x = ( double * )( mGeometry + 5 );
03456 y = ( double * )( mGeometry + 5 + sizeof( double ) );
03457 if ( *x < xmin )
03458 {
03459 xmin = *x;
03460 }
03461 if ( *x > xmax )
03462 {
03463 xmax = *x;
03464 }
03465 if ( *y < ymin )
03466 {
03467 ymin = *y;
03468 }
03469 if ( *y > ymax )
03470 {
03471 ymax = *y;
03472 }
03473 break;
03474 case QGis::WKBMultiPoint25D:
03475 hasZValue = true;
03476 case QGis::WKBMultiPoint:
03477 {
03478 ptr = mGeometry + 1 + sizeof( int );
03479 nPoints = ( int * ) ptr;
03480 ptr += sizeof( int );
03481 for ( idx = 0; idx < *nPoints; idx++ )
03482 {
03483 ptr += ( 1 + sizeof( int ) );
03484 x = ( double * ) ptr;
03485 ptr += sizeof( double );
03486 y = ( double * ) ptr;
03487 ptr += sizeof( double );
03488 if ( hasZValue )
03489 {
03490 ptr += sizeof( double );
03491 }
03492 if ( *x < xmin )
03493 {
03494 xmin = *x;
03495 }
03496 if ( *x > xmax )
03497 {
03498 xmax = *x;
03499 }
03500 if ( *y < ymin )
03501 {
03502 ymin = *y;
03503 }
03504 if ( *y > ymax )
03505 {
03506 ymax = *y;
03507 }
03508 }
03509 break;
03510 }
03511 case QGis::WKBLineString25D:
03512 hasZValue = true;
03513 case QGis::WKBLineString:
03514 {
03515
03516 ptr = mGeometry + 5;
03517 nPoints = ( int * ) ptr;
03518 ptr = mGeometry + 1 + 2 * sizeof( int );
03519 for ( idx = 0; idx < *nPoints; idx++ )
03520 {
03521 x = ( double * ) ptr;
03522 ptr += sizeof( double );
03523 y = ( double * ) ptr;
03524 ptr += sizeof( double );
03525 if ( hasZValue )
03526 {
03527 ptr += sizeof( double );
03528 }
03529 if ( *x < xmin )
03530 {
03531 xmin = *x;
03532 }
03533 if ( *x > xmax )
03534 {
03535 xmax = *x;
03536 }
03537 if ( *y < ymin )
03538 {
03539 ymin = *y;
03540 }
03541 if ( *y > ymax )
03542 {
03543 ymax = *y;
03544 }
03545 }
03546 break;
03547 }
03548 case QGis::WKBMultiLineString25D:
03549 hasZValue = true;
03550 case QGis::WKBMultiLineString:
03551 {
03552 numLineStrings = ( int )( mGeometry[5] );
03553 ptr = mGeometry + 9;
03554 for ( jdx = 0; jdx < numLineStrings; jdx++ )
03555 {
03556
03557 lsb = *ptr;
03558 ptr += 5;
03559 nPoints = ( int * ) ptr;
03560 ptr += sizeof( int );
03561 for ( idx = 0; idx < *nPoints; idx++ )
03562 {
03563 x = ( double * ) ptr;
03564 ptr += sizeof( double );
03565 y = ( double * ) ptr;
03566 ptr += sizeof( double );
03567 if ( hasZValue )
03568 {
03569 ptr += sizeof( double );
03570 }
03571 if ( *x < xmin )
03572 {
03573 xmin = *x;
03574 }
03575 if ( *x > xmax )
03576 {
03577 xmax = *x;
03578 }
03579 if ( *y < ymin )
03580 {
03581 ymin = *y;
03582 }
03583 if ( *y > ymax )
03584 {
03585 ymax = *y;
03586 }
03587 }
03588 }
03589 break;
03590 }
03591 case QGis::WKBPolygon25D:
03592 hasZValue = true;
03593 case QGis::WKBPolygon:
03594 {
03595
03596 numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
03597 ptr = mGeometry + 1 + 2 * sizeof( int );
03598 for ( idx = 0; idx < *numRings; idx++ )
03599 {
03600
03601 nPoints = ( int * ) ptr;
03602 ptr += 4;
03603 for ( jdx = 0; jdx < *nPoints; jdx++ )
03604 {
03605
03606 x = ( double * ) ptr;
03607 ptr += sizeof( double );
03608 y = ( double * ) ptr;
03609 ptr += sizeof( double );
03610 if ( hasZValue )
03611 {
03612 ptr += sizeof( double );
03613 }
03614 if ( *x < xmin )
03615 {
03616 xmin = *x;
03617 }
03618 if ( *x > xmax )
03619 {
03620 xmax = *x;
03621 }
03622 if ( *y < ymin )
03623 {
03624 ymin = *y;
03625 }
03626 if ( *y > ymax )
03627 {
03628 ymax = *y;
03629 }
03630 }
03631 }
03632 break;
03633 }
03634 case QGis::WKBMultiPolygon25D:
03635 hasZValue = true;
03636 case QGis::WKBMultiPolygon:
03637 {
03638
03639 ptr = mGeometry + 5;
03640 numPolygons = ( int * ) ptr;
03641 ptr += 4;
03642
03643 for ( kdx = 0; kdx < *numPolygons; kdx++ )
03644 {
03645
03646
03647 ptr += 5;
03648 numRings = ( int * ) ptr;
03649 ptr += 4;
03650 for ( idx = 0; idx < *numRings; idx++ )
03651 {
03652
03653 nPoints = ( int * ) ptr;
03654 ptr += 4;
03655 for ( jdx = 0; jdx < *nPoints; jdx++ )
03656 {
03657
03658 x = ( double * ) ptr;
03659 ptr += sizeof( double );
03660 y = ( double * ) ptr;
03661 ptr += sizeof( double );
03662 if ( hasZValue )
03663 {
03664 ptr += sizeof( double );
03665 }
03666 if ( *x < xmin )
03667 {
03668 xmin = *x;
03669 }
03670 if ( *x > xmax )
03671 {
03672 xmax = *x;
03673 }
03674 if ( *y < ymin )
03675 {
03676 ymin = *y;
03677 }
03678 if ( *y > ymax )
03679 {
03680 ymax = *y;
03681 }
03682 }
03683 }
03684 }
03685 break;
03686 }
03687
03688 default:
03689 QgsDebugMsg( "Unknown WkbType ENCOUNTERED" );
03690 return QgsRectangle( 0, 0, 0, 0 );
03691 break;
03692
03693 }
03694 return QgsRectangle( xmin, ymin, xmax, ymax );
03695 }
03696
03697 bool QgsGeometry::intersects( const QgsRectangle& r )
03698 {
03699 QgsGeometry* g = fromRect( r );
03700 bool res = intersects( g );
03701 delete g;
03702 return res;
03703 }
03704
03705 bool QgsGeometry::intersects( QgsGeometry* geometry )
03706 {
03707 try
03708 {
03709
03710 exportWkbToGeos();
03711 geometry->exportWkbToGeos();
03712
03713 if ( !mGeos || !geometry->mGeos )
03714 {
03715 QgsDebugMsg( "GEOS geometry not available!" );
03716 return false;
03717 }
03718
03719 return GEOSIntersects( mGeos, geometry->mGeos );
03720 }
03721 CATCH_GEOS( false )
03722 }
03723
03724
03725 bool QgsGeometry::contains( QgsPoint* p )
03726 {
03727 exportWkbToGeos();
03728
03729 if ( !mGeos )
03730 {
03731 QgsDebugMsg( "GEOS geometry not available!" );
03732 return false;
03733 }
03734
03735 GEOSGeometry *geosPoint = 0;
03736
03737 bool returnval = false;
03738
03739 try
03740 {
03741 geosPoint = createGeosPoint( *p );
03742 returnval = GEOSContains( mGeos, geosPoint );
03743 }
03744 catch ( GEOSException &e )
03745 {
03746 Q_UNUSED( e );
03747 returnval = false;
03748 }
03749
03750 if ( geosPoint )
03751 GEOSGeom_destroy( geosPoint );
03752
03753 return returnval;
03754 }
03755
03756 bool QgsGeometry::geosRelOp(
03757 char( *op )( const GEOSGeometry*, const GEOSGeometry * ),
03758 QgsGeometry *a,
03759 QgsGeometry *b )
03760 {
03761 try
03762 {
03763
03764 a->exportWkbToGeos();
03765 b->exportWkbToGeos();
03766
03767 if ( !a->mGeos || !b->mGeos )
03768 {
03769 QgsDebugMsg( "GEOS geometry not available!" );
03770 return false;
03771 }
03772 return op( a->mGeos, b->mGeos );
03773 }
03774 CATCH_GEOS( false )
03775 }
03776
03777 bool QgsGeometry::contains( QgsGeometry* geometry )
03778 {
03779 return geosRelOp( GEOSContains, this, geometry );
03780 }
03781
03782 bool QgsGeometry::disjoint( QgsGeometry* geometry )
03783 {
03784 return geosRelOp( GEOSDisjoint, this, geometry );
03785 }
03786
03787 bool QgsGeometry::equals( QgsGeometry* geometry )
03788 {
03789 return geosRelOp( GEOSEquals, this, geometry );
03790 }
03791
03792 bool QgsGeometry::touches( QgsGeometry* geometry )
03793 {
03794 return geosRelOp( GEOSTouches, this, geometry );
03795 }
03796
03797 bool QgsGeometry::overlaps( QgsGeometry* geometry )
03798 {
03799 return geosRelOp( GEOSOverlaps, this, geometry );
03800 }
03801
03802 bool QgsGeometry::within( QgsGeometry* geometry )
03803 {
03804 return geosRelOp( GEOSWithin, this, geometry );
03805 }
03806
03807 bool QgsGeometry::crosses( QgsGeometry* geometry )
03808 {
03809 return geosRelOp( GEOSCrosses, this, geometry );
03810 }
03811
03812 QString QgsGeometry::exportToWkt()
03813 {
03814 QgsDebugMsg( "entered." );
03815
03816
03817 if ( mDirtyWkb )
03818 {
03819 exportGeosToWkb();
03820 }
03821
03822 if ( !mGeometry )
03823 {
03824 QgsDebugMsg( "WKB geometry not available!" );
03825 return false;
03826 }
03827
03828 QGis::WkbType wkbType;
03829 bool hasZValue = false;
03830 double *x, *y;
03831
03832 QString mWkt;
03833
03834
03835
03836 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) );
03837
03838 switch ( wkbType )
03839 {
03840 case QGis::WKBPoint25D:
03841 case QGis::WKBPoint:
03842 {
03843 mWkt += "POINT(";
03844 x = ( double * )( mGeometry + 5 );
03845 mWkt += QString::number( *x, 'f', 6 );
03846 mWkt += " ";
03847 y = ( double * )( mGeometry + 5 + sizeof( double ) );
03848 mWkt += QString::number( *y, 'f', 6 );
03849 mWkt += ")";
03850 return mWkt;
03851 }
03852
03853 case QGis::WKBLineString25D:
03854 hasZValue = true;
03855 case QGis::WKBLineString:
03856 {
03857 QgsDebugMsg( "LINESTRING found" );
03858 unsigned char *ptr;
03859 int *nPoints;
03860 int idx;
03861
03862 mWkt += "LINESTRING(";
03863
03864 ptr = mGeometry + 5;
03865 nPoints = ( int * ) ptr;
03866 ptr = mGeometry + 1 + 2 * sizeof( int );
03867 for ( idx = 0; idx < *nPoints; ++idx )
03868 {
03869 if ( idx != 0 )
03870 {
03871 mWkt += ", ";
03872 }
03873 x = ( double * ) ptr;
03874 mWkt += QString::number( *x, 'f', 6 );
03875 mWkt += " ";
03876 ptr += sizeof( double );
03877 y = ( double * ) ptr;
03878 mWkt += QString::number( *y, 'f', 6 );
03879 ptr += sizeof( double );
03880 if ( hasZValue )
03881 {
03882 ptr += sizeof( double );
03883 }
03884 }
03885 mWkt += ")";
03886 return mWkt;
03887 }
03888
03889 case QGis::WKBPolygon25D:
03890 hasZValue = true;
03891 case QGis::WKBPolygon:
03892 {
03893 QgsDebugMsg( "POLYGON found" );
03894 unsigned char *ptr;
03895 int idx, jdx;
03896 int *numRings, *nPoints;
03897
03898 mWkt += "POLYGON(";
03899
03900 numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
03901 if ( !( *numRings ) )
03902 {
03903 return QString();
03904 }
03905 int *ringStart;
03906 int *ringNumPoints;
03907 ringStart = new int[*numRings];
03908 ringNumPoints = new int[*numRings];
03909 ptr = mGeometry + 1 + 2 * sizeof( int );
03910 for ( idx = 0; idx < *numRings; idx++ )
03911 {
03912 if ( idx != 0 )
03913 {
03914 mWkt += ",";
03915 }
03916 mWkt += "(";
03917
03918 nPoints = ( int * ) ptr;
03919 ringNumPoints[idx] = *nPoints;
03920 ptr += 4;
03921
03922 for ( jdx = 0; jdx < *nPoints; jdx++ )
03923 {
03924 if ( jdx != 0 )
03925 {
03926 mWkt += ",";
03927 }
03928 x = ( double * ) ptr;
03929 mWkt += QString::number( *x, 'f', 6 );
03930 mWkt += " ";
03931 ptr += sizeof( double );
03932 y = ( double * ) ptr;
03933 mWkt += QString::number( *y, 'f', 6 );
03934 ptr += sizeof( double );
03935 if ( hasZValue )
03936 {
03937 ptr += sizeof( double );
03938 }
03939 }
03940 mWkt += ")";
03941 }
03942 mWkt += ")";
03943 delete [] ringStart;
03944 delete [] ringNumPoints;
03945 return mWkt;
03946 }
03947
03948 case QGis::WKBMultiPoint25D:
03949 hasZValue = true;
03950 case QGis::WKBMultiPoint:
03951 {
03952 unsigned char *ptr;
03953 int idx;
03954 int *nPoints;
03955
03956 mWkt += "MULTIPOINT(";
03957 nPoints = ( int* )( mGeometry + 5 );
03958 ptr = mGeometry + 5 + sizeof( int );
03959 for ( idx = 0; idx < *nPoints; ++idx )
03960 {
03961 ptr += ( 1 + sizeof( int ) );
03962 if ( idx != 0 )
03963 {
03964 mWkt += ", ";
03965 }
03966 x = ( double * )( ptr );
03967 mWkt += QString::number( *x, 'f', 6 );
03968 mWkt += " ";
03969 ptr += sizeof( double );
03970 y = ( double * )( ptr );
03971 mWkt += QString::number( *y, 'f', 6 );
03972 ptr += sizeof( double );
03973 if ( hasZValue )
03974 {
03975 ptr += sizeof( double );
03976 }
03977 }
03978 mWkt += ")";
03979 return mWkt;
03980 }
03981
03982 case QGis::WKBMultiLineString25D:
03983 hasZValue = true;
03984 case QGis::WKBMultiLineString:
03985 {
03986 QgsDebugMsg( "MULTILINESTRING found" );
03987 unsigned char *ptr;
03988 int idx, jdx, numLineStrings;
03989 int *nPoints;
03990
03991 mWkt += "MULTILINESTRING(";
03992 numLineStrings = ( int )( mGeometry[5] );
03993 ptr = mGeometry + 9;
03994 for ( jdx = 0; jdx < numLineStrings; jdx++ )
03995 {
03996 if ( jdx != 0 )
03997 {
03998 mWkt += ", ";
03999 }
04000 mWkt += "(";
04001 ptr += 5;
04002 nPoints = ( int * ) ptr;
04003 ptr += sizeof( int );
04004 for ( idx = 0; idx < *nPoints; idx++ )
04005 {
04006 if ( idx != 0 )
04007 {
04008 mWkt += ", ";
04009 }
04010 x = ( double * ) ptr;
04011 mWkt += QString::number( *x, 'f', 6 );
04012 ptr += sizeof( double );
04013 mWkt += " ";
04014 y = ( double * ) ptr;
04015 mWkt += QString::number( *y, 'f', 6 );
04016 ptr += sizeof( double );
04017 if ( hasZValue )
04018 {
04019 ptr += sizeof( double );
04020 }
04021 }
04022 mWkt += ")";
04023 }
04024 mWkt += ")";
04025 return mWkt;
04026 }
04027
04028 case QGis::WKBMultiPolygon25D:
04029 hasZValue = true;
04030 case QGis::WKBMultiPolygon:
04031 {
04032 QgsDebugMsg( "MULTIPOLYGON found" );
04033 unsigned char *ptr;
04034 int idx, jdx, kdx;
04035 int *numPolygons, *numRings, *nPoints;
04036
04037 mWkt += "MULTIPOLYGON(";
04038 ptr = mGeometry + 5;
04039 numPolygons = ( int * ) ptr;
04040 ptr = mGeometry + 9;
04041 for ( kdx = 0; kdx < *numPolygons; kdx++ )
04042 {
04043 if ( kdx != 0 )
04044 {
04045 mWkt += ",";
04046 }
04047 mWkt += "(";
04048 ptr += 5;
04049 numRings = ( int * ) ptr;
04050 ptr += 4;
04051 for ( idx = 0; idx < *numRings; idx++ )
04052 {
04053 if ( idx != 0 )
04054 {
04055 mWkt += ",";
04056 }
04057 mWkt += "(";
04058 nPoints = ( int * ) ptr;
04059 ptr += 4;
04060 for ( jdx = 0; jdx < *nPoints; jdx++ )
04061 {
04062 if ( jdx != 0 )
04063 {
04064 mWkt += ",";
04065 }
04066 x = ( double * ) ptr;
04067 mWkt += QString::number( *x, 'f', 6 );
04068 ptr += sizeof( double );
04069 mWkt += " ";
04070 y = ( double * ) ptr;
04071 mWkt += QString::number( *y, 'f', 6 );
04072 ptr += sizeof( double );
04073 if ( hasZValue )
04074 {
04075 ptr += sizeof( double );
04076 }
04077 }
04078 mWkt += ")";
04079 }
04080 mWkt += ")";
04081 }
04082 mWkt += ")";
04083 return mWkt;
04084 }
04085
04086 default:
04087 QgsDebugMsg( "error: mGeometry type not recognized" );
04088 return QString();
04089 }
04090 }
04091
04092 bool QgsGeometry::exportWkbToGeos()
04093 {
04094 QgsDebugMsgLevel( "entered.", 3 );
04095
04096 if ( !mDirtyGeos )
04097 {
04098
04099 return true;
04100 }
04101
04102 if ( mGeos )
04103 {
04104 GEOSGeom_destroy( mGeos );
04105 mGeos = 0;
04106 }
04107
04108
04109 if ( !mGeometry )
04110 {
04111
04112 mDirtyGeos = false;
04113 return true;
04114 }
04115
04116 double *x;
04117 double *y;
04118 int *nPoints;
04119 int *numRings;
04120 int *numPolygons;
04121 int numLineStrings;
04122 int idx, jdx, kdx;
04123 unsigned char *ptr;
04124 char lsb;
04125 QgsPoint pt;
04126 QGis::WkbType wkbtype;
04127 bool hasZValue = false;
04128
04129
04130 memcpy( &wkbtype, &( mGeometry[1] ), sizeof( int ) );
04131
04132 try
04133 {
04134 switch ( wkbtype )
04135 {
04136 case QGis::WKBPoint25D:
04137 case QGis::WKBPoint:
04138 {
04139 x = ( double * )( mGeometry + 5 );
04140 y = ( double * )( mGeometry + 5 + sizeof( double ) );
04141
04142 mGeos = createGeosPoint( QgsPoint( *x, *y ) );
04143 mDirtyGeos = false;
04144 break;
04145 }
04146
04147 case QGis::WKBMultiPoint25D:
04148 hasZValue = true;
04149 case QGis::WKBMultiPoint:
04150 {
04151 QVector<GEOSGeometry *> points;
04152
04153 ptr = mGeometry + 5;
04154 nPoints = ( int * ) ptr;
04155 ptr = mGeometry + 1 + 2 * sizeof( int );
04156 for ( idx = 0; idx < *nPoints; idx++ )
04157 {
04158 ptr += ( 1 + sizeof( int ) );
04159 x = ( double * ) ptr;
04160 ptr += sizeof( double );
04161 y = ( double * ) ptr;
04162 ptr += sizeof( double );
04163 if ( hasZValue )
04164 {
04165 ptr += sizeof( double );
04166 }
04167 points << createGeosPoint( QgsPoint( *x, *y ) );
04168 }
04169 mGeos = createGeosCollection( GEOS_MULTIPOINT, points );
04170 mDirtyGeos = false;
04171 break;
04172 }
04173
04174 case QGis::WKBLineString25D:
04175 hasZValue = true;
04176 case QGis::WKBLineString:
04177 {
04178 QgsDebugMsgLevel( "Linestring found", 3 );
04179
04180 QgsPolyline sequence;
04181
04182 ptr = mGeometry + 5;
04183 nPoints = ( int * ) ptr;
04184 ptr = mGeometry + 1 + 2 * sizeof( int );
04185 for ( idx = 0; idx < *nPoints; idx++ )
04186 {
04187 x = ( double * ) ptr;
04188 ptr += sizeof( double );
04189 y = ( double * ) ptr;
04190 ptr += sizeof( double );
04191 if ( hasZValue )
04192 {
04193 ptr += sizeof( double );
04194 }
04195
04196 sequence << QgsPoint( *x, *y );
04197 }
04198 mDirtyGeos = false;
04199 mGeos = createGeosLineString( sequence );
04200 break;
04201 }
04202
04203 case QGis::WKBMultiLineString25D:
04204 hasZValue = true;
04205 case QGis::WKBMultiLineString:
04206 {
04207 QVector<GEOSGeometry*> lines;
04208 numLineStrings = ( int )( mGeometry[5] );
04209 ptr = ( mGeometry + 9 );
04210 for ( jdx = 0; jdx < numLineStrings; jdx++ )
04211 {
04212 QgsPolyline sequence;
04213
04214
04215 lsb = *ptr;
04216 ptr += 5;
04217 nPoints = ( int * ) ptr;
04218 ptr += sizeof( int );
04219 for ( idx = 0; idx < *nPoints; idx++ )
04220 {
04221 x = ( double * ) ptr;
04222 ptr += sizeof( double );
04223 y = ( double * ) ptr;
04224 ptr += sizeof( double );
04225 if ( hasZValue )
04226 {
04227 ptr += sizeof( double );
04228 }
04229 sequence << QgsPoint( *x, *y );
04230 }
04231 lines << createGeosLineString( sequence );
04232 }
04233 mGeos = createGeosCollection( GEOS_MULTILINESTRING, lines );
04234 mDirtyGeos = false;
04235 break;
04236 }
04237
04238 case QGis::WKBPolygon25D:
04239 hasZValue = true;
04240 case QGis::WKBPolygon:
04241 {
04242 QgsDebugMsgLevel( "Polygon found", 3 );
04243
04244
04245 numRings = ( int * )( mGeometry + 1 + sizeof( int ) );
04246 ptr = mGeometry + 1 + 2 * sizeof( int );
04247
04248 QVector<GEOSGeometry*> rings;
04249
04250 for ( idx = 0; idx < *numRings; idx++ )
04251 {
04252
04253
04254 QgsPolyline sequence;
04255
04256
04257 nPoints = ( int * ) ptr;
04258 ptr += 4;
04259 for ( jdx = 0; jdx < *nPoints; jdx++ )
04260 {
04261
04262 x = ( double * ) ptr;
04263 ptr += sizeof( double );
04264 y = ( double * ) ptr;
04265 ptr += sizeof( double );
04266 if ( hasZValue )
04267 {
04268 ptr += sizeof( double );
04269 }
04270 sequence << QgsPoint( *x, *y );
04271 }
04272
04273 rings << createGeosLinearRing( sequence );
04274 }
04275 mGeos = createGeosPolygon( rings );
04276 mDirtyGeos = false;
04277 break;
04278 }
04279
04280 case QGis::WKBMultiPolygon25D:
04281 hasZValue = true;
04282 case QGis::WKBMultiPolygon:
04283 {
04284 QgsDebugMsgLevel( "Multipolygon found", 3 );
04285
04286 QVector<GEOSGeometry*> polygons;
04287
04288
04289 ptr = mGeometry + 5;
04290 numPolygons = ( int * ) ptr;
04291 ptr = mGeometry + 9;
04292 for ( kdx = 0; kdx < *numPolygons; kdx++ )
04293 {
04294
04295 QVector<GEOSGeometry*> rings;
04296
04297
04298
04299 ptr += 5;
04300 numRings = ( int * ) ptr;
04301 ptr += 4;
04302 for ( idx = 0; idx < *numRings; idx++ )
04303 {
04304
04305
04306 QgsPolyline sequence;
04307
04308
04309 nPoints = ( int * ) ptr;
04310 ptr += 4;
04311 for ( jdx = 0; jdx < *nPoints; jdx++ )
04312 {
04313
04314 x = ( double * ) ptr;
04315 ptr += sizeof( double );
04316 y = ( double * ) ptr;
04317 ptr += sizeof( double );
04318 if ( hasZValue )
04319 {
04320 ptr += sizeof( double );
04321 }
04322 sequence << QgsPoint( *x, *y );
04323 }
04324
04325 rings << createGeosLinearRing( sequence );
04326 }
04327
04328 polygons << createGeosPolygon( rings );
04329 }
04330 mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
04331 mDirtyGeos = false;
04332 break;
04333 }
04334
04335 default:
04336 return false;
04337 }
04338 }
04339 CATCH_GEOS( false )
04340
04341 return true;
04342 }
04343
04344 bool QgsGeometry::exportGeosToWkb()
04345 {
04346
04347 if ( !mDirtyWkb )
04348 {
04349
04350 return true;
04351 }
04352
04353
04354 if ( mGeometry )
04355 {
04356 delete [] mGeometry;
04357 mGeometry = 0;
04358 }
04359
04360 if ( !mGeos )
04361 {
04362
04363 mDirtyWkb = false;
04364 return true;
04365 }
04366
04367
04368 char byteOrder = QgsApplication::endian();
04369
04370 switch ( GEOSGeomTypeId( mGeos ) )
04371 {
04372 case GEOS_POINT:
04373 {
04374 mGeometrySize = 1 +
04375 4 +
04376 2 * sizeof( double );
04377 mGeometry = new unsigned char[mGeometrySize];
04378
04379
04380 memcpy( mGeometry, &byteOrder, 1 );
04381
04382
04383 int wkbType = QGis::WKBPoint;
04384 memcpy( mGeometry + 1, &wkbType, 4 );
04385
04386 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos );
04387
04388 double x, y;
04389 GEOSCoordSeq_getX( cs, 0, &x );
04390 GEOSCoordSeq_getY( cs, 0, &y );
04391
04392 memcpy( mGeometry + 5, &x, sizeof( double ) );
04393 memcpy( mGeometry + 13, &y, sizeof( double ) );
04394
04395 mDirtyWkb = false;
04396 return true;
04397 }
04398
04399 case GEOS_LINESTRING:
04400 {
04401
04402
04403 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos );
04404 unsigned int numPoints;
04405 GEOSCoordSeq_getSize( cs, &numPoints );
04406
04407
04408 mGeometrySize = 1 +
04409 4 +
04410 4 +
04411 (( sizeof( double ) +
04412 sizeof( double ) ) * numPoints );
04413
04414 mGeometry = new unsigned char[mGeometrySize];
04415
04416 unsigned char* ptr = mGeometry;
04417
04418
04419 memcpy( ptr, &byteOrder, 1 );
04420 ptr += 1;
04421
04422
04423 int wkbType = QGis::WKBLineString;
04424 memcpy( ptr, &wkbType, 4 );
04425 ptr += 4;
04426
04427
04428 memcpy( ptr, &numPoints, 4 );
04429 ptr += 4;
04430
04431 const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( mGeos );
04432
04433
04434 for ( unsigned int n = 0; n < numPoints; n++ )
04435 {
04436
04437 GEOSCoordSeq_getX( sequence, n, ( double * )ptr );
04438 ptr += sizeof( double );
04439
04440
04441 GEOSCoordSeq_getY( sequence, n, ( double * )ptr );
04442 ptr += sizeof( double );
04443 }
04444
04445 mDirtyWkb = false;
04446 return true;
04447
04448
04449 }
04450
04451 case GEOS_LINEARRING:
04452 {
04453
04454 break;
04455 }
04456
04457 case GEOS_POLYGON:
04458 {
04459 int geometrySize;
04460 unsigned int nPointsInRing = 0;
04461
04462
04463 geometrySize = 1 + 2 * sizeof( int );
04464 const GEOSGeometry *theRing = GEOSGetExteriorRing( mGeos );
04465 if ( theRing )
04466 {
04467 geometrySize += sizeof( int );
04468 geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
04469 }
04470 for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); ++i )
04471 {
04472 geometrySize += sizeof( int );
04473 theRing = GEOSGetInteriorRingN( mGeos, i );
04474 if ( theRing )
04475 {
04476 geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
04477 }
04478 }
04479
04480 mGeometry = new unsigned char[geometrySize];
04481 mGeometrySize = geometrySize;
04482
04483
04484 int position = 0;
04485
04486 memcpy( mGeometry, &byteOrder, 1 );
04487 position += 1;
04488 int wkbtype = QGis::WKBPolygon;
04489 memcpy( &mGeometry[position], &wkbtype, sizeof( int ) );
04490 position += sizeof( int );
04491 int nRings = GEOSGetNumInteriorRings( mGeos ) + 1;
04492 memcpy( &mGeometry[position], &nRings, sizeof( int ) );
04493 position += sizeof( int );
04494
04495
04496 theRing = GEOSGetExteriorRing( mGeos );
04497 if ( theRing )
04498 {
04499 nPointsInRing = getNumGeosPoints( theRing );
04500 memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) );
04501 position += sizeof( int );
04502
04503 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
04504 unsigned int n;
04505 GEOSCoordSeq_getSize( cs, &n );
04506
04507 for ( unsigned int j = 0; j < n; ++j )
04508 {
04509 GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] );
04510 position += sizeof( double );
04511 GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] );
04512 position += sizeof( double );
04513 }
04514 }
04515
04516
04517 for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); i++ )
04518 {
04519 theRing = GEOSGetInteriorRingN( mGeos, i );
04520
04521 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
04522 GEOSCoordSeq_getSize( cs, &nPointsInRing );
04523
04524 memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) );
04525 position += sizeof( int );
04526
04527 for ( unsigned int j = 0; j < nPointsInRing; j++ )
04528 {
04529 GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] );
04530 position += sizeof( double );
04531 GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] );
04532 position += sizeof( double );
04533 }
04534 }
04535 mDirtyWkb = false;
04536 return true;
04537 }
04538 break;
04539
04540 case GEOS_MULTIPOINT:
04541 {
04542
04543 int geometrySize = 1 + 2 * sizeof( int );
04544 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04545 {
04546 geometrySize += 1 + sizeof( int ) + 2 * sizeof( double );
04547 }
04548
04549 mGeometry = new unsigned char[geometrySize];
04550 mGeometrySize = geometrySize;
04551 int wkbPosition = 0;
04552
04553 memcpy( mGeometry, &byteOrder, 1 );
04554 wkbPosition += 1;
04555 int wkbtype = QGis::WKBMultiPoint;
04556 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
04557 wkbPosition += sizeof( int );
04558 int numPoints = GEOSGetNumGeometries( mGeos );
04559 memcpy( &mGeometry[wkbPosition], &numPoints, sizeof( int ) );
04560 wkbPosition += sizeof( int );
04561
04562 int pointType = QGis::WKBPoint;
04563 const GEOSGeometry *currentPoint = 0;
04564
04565 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04566 {
04567
04568 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
04569 wkbPosition += 1;
04570 memcpy( &mGeometry[wkbPosition], &pointType, sizeof( int ) );
04571 wkbPosition += sizeof( int );
04572
04573 currentPoint = GEOSGetGeometryN( mGeos, i );
04574
04575 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( currentPoint );
04576
04577 GEOSCoordSeq_getX( cs, 0, ( double* )&mGeometry[wkbPosition] );
04578 wkbPosition += sizeof( double );
04579 GEOSCoordSeq_getY( cs, 0, ( double* )&mGeometry[wkbPosition] );
04580 wkbPosition += sizeof( double );
04581 }
04582 mDirtyWkb = false;
04583 return true;
04584 }
04585
04586 case GEOS_MULTILINESTRING:
04587 {
04588
04589 int geometrySize = 1 + 2 * sizeof( int );
04590 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04591 {
04592 geometrySize += 1 + 2 * sizeof( int );
04593 geometrySize += getNumGeosPoints( GEOSGetGeometryN( mGeos, i ) ) * 2 * sizeof( double );
04594 }
04595
04596 mGeometry = new unsigned char[geometrySize];
04597 mGeometrySize = geometrySize;
04598 int wkbPosition = 0;
04599
04600 memcpy( mGeometry, &byteOrder, 1 );
04601 wkbPosition += 1;
04602 int wkbtype = QGis::WKBMultiLineString;
04603 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
04604 wkbPosition += sizeof( int );
04605 int numLines = GEOSGetNumGeometries( mGeos );
04606 memcpy( &mGeometry[wkbPosition], &numLines, sizeof( int ) );
04607 wkbPosition += sizeof( int );
04608
04609
04610 int lineType = QGis::WKBLineString;
04611 const GEOSCoordSequence *cs = 0;
04612 unsigned int lineSize;
04613
04614 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04615 {
04616
04617 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
04618 wkbPosition += 1;
04619 memcpy( &mGeometry[wkbPosition], &lineType, sizeof( int ) );
04620 wkbPosition += sizeof( int );
04621
04622 cs = GEOSGeom_getCoordSeq( GEOSGetGeometryN( mGeos, i ) );
04623
04624
04625 GEOSCoordSeq_getSize( cs, &lineSize );
04626 memcpy( &mGeometry[wkbPosition], &lineSize, sizeof( int ) );
04627 wkbPosition += sizeof( int );
04628
04629
04630 for ( unsigned int j = 0; j < lineSize; ++j )
04631 {
04632 GEOSCoordSeq_getX( cs, j, ( double* )&mGeometry[wkbPosition] );
04633 wkbPosition += sizeof( double );
04634 GEOSCoordSeq_getY( cs, j, ( double* )&mGeometry[wkbPosition] );
04635 wkbPosition += sizeof( double );
04636 }
04637 }
04638 mDirtyWkb = false;
04639 return true;
04640 }
04641
04642 case GEOS_MULTIPOLYGON:
04643 {
04644
04645 int geometrySize = 1 + ( 2 * sizeof( int ) );
04646 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04647 {
04648 const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i );
04649 geometrySize += 1 + 2 * sizeof( int );
04650
04651 geometrySize += sizeof( int );
04652 const GEOSGeometry *exRing = GEOSGetExteriorRing( thePoly );
04653 geometrySize += 2 * sizeof( double ) * getNumGeosPoints( exRing );
04654
04655 const GEOSGeometry *intRing = 0;
04656 for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ )
04657 {
04658 geometrySize += sizeof( int );
04659 intRing = GEOSGetInteriorRingN( thePoly, j );
04660 geometrySize += 2 * sizeof( double ) * getNumGeosPoints( intRing );
04661 }
04662 }
04663
04664 mGeometry = new unsigned char[geometrySize];
04665 mGeometrySize = geometrySize;
04666 int wkbPosition = 0;
04667
04668 memcpy( mGeometry, &byteOrder, 1 );
04669 wkbPosition += 1;
04670 int wkbtype = QGis::WKBMultiPolygon;
04671 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) );
04672 wkbPosition += sizeof( int );
04673 int numPolygons = GEOSGetNumGeometries( mGeos );
04674 memcpy( &mGeometry[wkbPosition], &numPolygons, sizeof( int ) );
04675 wkbPosition += sizeof( int );
04676
04677
04678 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ )
04679 {
04680 const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i );
04681 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 );
04682 wkbPosition += 1;
04683 int polygonType = QGis::WKBPolygon;
04684 memcpy( &mGeometry[wkbPosition], &polygonType, sizeof( int ) );
04685 wkbPosition += sizeof( int );
04686 int numRings = GEOSGetNumInteriorRings( thePoly ) + 1;
04687 memcpy( &mGeometry[wkbPosition], &numRings, sizeof( int ) );
04688 wkbPosition += sizeof( int );
04689
04690
04691 const GEOSGeometry *theRing = GEOSGetExteriorRing( thePoly );
04692 int nPointsInRing = getNumGeosPoints( theRing );
04693 memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) );
04694 wkbPosition += sizeof( int );
04695 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
04696
04697 for ( int k = 0; k < nPointsInRing; ++k )
04698 {
04699 GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] );
04700 wkbPosition += sizeof( double );
04701 GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] );
04702 wkbPosition += sizeof( double );
04703 }
04704
04705
04706 for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ )
04707 {
04708 theRing = GEOSGetInteriorRingN( thePoly, j );
04709 nPointsInRing = getNumGeosPoints( theRing );
04710 memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) );
04711 wkbPosition += sizeof( int );
04712 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing );
04713
04714 for ( int k = 0; k < nPointsInRing; ++k )
04715 {
04716 GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] );
04717 wkbPosition += sizeof( double );
04718 GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] );
04719 wkbPosition += sizeof( double );
04720 }
04721 }
04722 }
04723 mDirtyWkb = false;
04724 return true;
04725 }
04726
04727 case GEOS_GEOMETRYCOLLECTION:
04728 {
04729
04730 QgsDebugMsg( "geometry collection - not supported" );
04731 break;
04732 }
04733
04734 }
04735
04736 return false;
04737 }
04738
04739 bool QgsGeometry::convertToMultiType()
04740 {
04741 if ( !mGeometry )
04742 {
04743 return false;
04744 }
04745
04746 QGis::WkbType geomType = wkbType();
04747
04748 if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D ||
04749 geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D ||
04750 geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown )
04751 {
04752 return false;
04753 }
04754
04755 int newGeomSize = mGeometrySize + 1 + 2 * sizeof( int );
04756 unsigned char* newGeometry = new unsigned char[newGeomSize];
04757
04758 int currentWkbPosition = 0;
04759
04760 char byteOrder = QgsApplication::endian();
04761 memcpy( &newGeometry[currentWkbPosition], &byteOrder, 1 );
04762 currentWkbPosition += 1;
04763
04764
04765
04766 QGis::WkbType newMultiType;
04767 switch ( geomType )
04768 {
04769 case QGis::WKBPoint:
04770 newMultiType = QGis::WKBMultiPoint;
04771 break;
04772 case QGis::WKBPoint25D:
04773 newMultiType = QGis::WKBMultiPoint25D;
04774 break;
04775 case QGis::WKBLineString:
04776 newMultiType = QGis::WKBMultiLineString;
04777 break;
04778 case QGis::WKBLineString25D:
04779 newMultiType = QGis::WKBMultiLineString25D;
04780 break;
04781 case QGis::WKBPolygon:
04782 newMultiType = QGis::WKBMultiPolygon;
04783 break;
04784 case QGis::WKBPolygon25D:
04785 newMultiType = QGis::WKBMultiPolygon25D;
04786 break;
04787 default:
04788 delete newGeometry;
04789 return false;
04790 }
04791 memcpy( &newGeometry[currentWkbPosition], &newMultiType, sizeof( int ) );
04792 currentWkbPosition += sizeof( int );
04793
04794
04795 int nGeometries = 1;
04796 memcpy( &newGeometry[currentWkbPosition], &nGeometries, sizeof( int ) );
04797 currentWkbPosition += sizeof( int );
04798
04799
04800 memcpy( &newGeometry[currentWkbPosition], mGeometry, mGeometrySize );
04801
04802 delete [] mGeometry;
04803 mGeometry = newGeometry;
04804 mGeometrySize = newGeomSize;
04805 mDirtyGeos = true;
04806 return true;
04807 }
04808
04809 void QgsGeometry::translateVertex( int& wkbPosition, double dx, double dy, bool hasZValue )
04810 {
04811 double x, y, translated_x, translated_y;
04812
04813
04814 x = *(( double * )( &( mGeometry[wkbPosition] ) ) );
04815 translated_x = x + dx;
04816 memcpy( &( mGeometry[wkbPosition] ), &translated_x, sizeof( double ) );
04817 wkbPosition += sizeof( double );
04818
04819
04820 y = *(( double * )( &( mGeometry[wkbPosition] ) ) );
04821 translated_y = y + dy;
04822 memcpy( &( mGeometry[wkbPosition] ), &translated_y, sizeof( double ) );
04823 wkbPosition += sizeof( double );
04824
04825 if ( hasZValue )
04826 {
04827 wkbPosition += sizeof( double );
04828 }
04829 }
04830
04831 void QgsGeometry::transformVertex( int& wkbPosition, const QgsCoordinateTransform& ct, bool hasZValue )
04832 {
04833 double x, y, z;
04834
04835
04836 x = *(( double * )( &( mGeometry[wkbPosition] ) ) );
04837 y = *(( double * )( &( mGeometry[wkbPosition + sizeof( double )] ) ) );
04838 z = 0.0;
04839
04840 ct.transformInPlace( x, y, z );
04841
04842
04843 memcpy( &( mGeometry[wkbPosition] ), &x, sizeof( double ) );
04844 wkbPosition += sizeof( double );
04845
04846
04847 memcpy( &( mGeometry[wkbPosition] ), &y, sizeof( double ) );
04848 wkbPosition += sizeof( double );
04849
04850 if ( hasZValue )
04851 {
04852 wkbPosition += sizeof( double );
04853 }
04854 }
04855
04856 int QgsGeometry::splitLinearGeometry( GEOSGeometry *splitLine, QList<QgsGeometry*>& newGeometries )
04857 {
04858 if ( !splitLine )
04859 {
04860 return 2;
04861 }
04862
04863 if ( !mGeos || mDirtyGeos )
04864 {
04865 if ( !exportWkbToGeos() )
04866 return 5;
04867 }
04868
04869
04870 if ( !GEOSIntersects( splitLine, mGeos ) )
04871 {
04872 return 1;
04873 }
04874
04875
04876 GEOSGeometry* nodedGeometry = nodeGeometries( splitLine, mGeos );
04877 if ( !nodedGeometry )
04878 {
04879 return 3;
04880 }
04881
04882 GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
04883 if ( !mergedLines )
04884 {
04885 GEOSGeom_destroy( nodedGeometry );
04886 return 4;
04887 }
04888
04889 QVector<GEOSGeometry*> testedGeometries;
04890
04891 for ( int i = 0; i < GEOSGetNumGeometries( mergedLines ); i++ )
04892 {
04893 const GEOSGeometry *testing = GEOSGetGeometryN( mergedLines, i );
04894 if ( lineContainedInLine( testing, mGeos ) == 1 )
04895 {
04896 testedGeometries << GEOSGeom_clone( testing );
04897 }
04898 }
04899
04900 mergeGeometriesMultiTypeSplit( testedGeometries );
04901
04902 if ( testedGeometries.size() > 0 )
04903 {
04904 GEOSGeom_destroy( mGeos );
04905 mGeos = testedGeometries[0];
04906 mDirtyWkb = true;
04907 }
04908
04909 for ( int i = 1; i < testedGeometries.size(); ++i )
04910 newGeometries << fromGeosGeom( testedGeometries[i] );
04911
04912 GEOSGeom_destroy( nodedGeometry );
04913 GEOSGeom_destroy( mergedLines );
04914 return 0;
04915 }
04916
04917 int QgsGeometry::splitPolygonGeometry( GEOSGeometry* splitLine, QList<QgsGeometry*>& newGeometries )
04918 {
04919 if ( !splitLine )
04920 {
04921 return 2;
04922 }
04923
04924 if ( !mGeos || mDirtyGeos )
04925 {
04926 if ( !exportWkbToGeos() )
04927 return 5;
04928 }
04929
04930
04931 if ( !GEOSIntersects( splitLine, mGeos ) )
04932 {
04933 return 1;
04934 }
04935
04936
04937 GEOSGeometry *nodedGeometry = nodeGeometries( splitLine, mGeos );
04938 if ( !nodedGeometry )
04939 {
04940 return 2;
04941 }
04942
04943 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
04944 ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=1)))
04945 GEOSGeometry *cutEdges = GEOSPolygonizer_getCutEdges( &nodedGeometry, 1 );
04946 if ( cutEdges )
04947 {
04948 if ( numberOfGeometries( cutEdges ) > 0 )
04949 {
04950 GEOSGeom_destroy( cutEdges );
04951 GEOSGeom_destroy( nodedGeometry );
04952 return 3;
04953 }
04954
04955 GEOSGeom_destroy( cutEdges );
04956 }
04957 #endif
04958
04959 GEOSGeometry *polygons = GEOSPolygonize( &nodedGeometry, 1 );
04960 if ( !polygons || numberOfGeometries( polygons ) == 0 )
04961 {
04962 if ( polygons )
04963 GEOSGeom_destroy( polygons );
04964
04965 GEOSGeom_destroy( nodedGeometry );
04966
04967 return 4;
04968 }
04969
04970 GEOSGeom_destroy( nodedGeometry );
04971
04972
04973
04974 QVector<GEOSGeometry*> testedGeometries;
04975 GEOSGeometry *intersectGeometry = 0;
04976
04977
04978
04979
04980 for ( int i = 0; i < numberOfGeometries( polygons ); i++ )
04981 {
04982 const GEOSGeometry *polygon = GEOSGetGeometryN( polygons, i );
04983 intersectGeometry = GEOSIntersection( mGeos, polygon );
04984 if ( !intersectGeometry )
04985 {
04986 QgsDebugMsg( "intersectGeometry is NULL" );
04987 continue;
04988 }
04989
04990 double intersectionArea;
04991 GEOSArea( intersectGeometry, &intersectionArea );
04992
04993 double polygonArea;
04994 GEOSArea( polygon, &polygonArea );
04995
04996 const double areaRatio = intersectionArea / polygonArea;
04997 if ( areaRatio > 0.99 && areaRatio < 1.01 )
04998 testedGeometries << GEOSGeom_clone( polygon );
04999
05000 GEOSGeom_destroy( intersectGeometry );
05001 }
05002
05003 bool splitDone = true;
05004 int nGeometriesThis = numberOfGeometries( mGeos );
05005 if ( testedGeometries.size() == nGeometriesThis )
05006 {
05007 splitDone = false;
05008 }
05009
05010 mergeGeometriesMultiTypeSplit( testedGeometries );
05011
05012
05013 if ( !splitDone )
05014 {
05015 for ( int i = 0; i < testedGeometries.size(); ++i )
05016 {
05017 GEOSGeom_destroy( testedGeometries[i] );
05018 }
05019 return 1;
05020 }
05021 else if ( testedGeometries.size() > 0 )
05022 {
05023 GEOSGeom_destroy( mGeos );
05024 mGeos = testedGeometries[0];
05025 mDirtyWkb = true;
05026 }
05027
05028 for ( int i = 1; i < testedGeometries.size(); ++i )
05029 {
05030 newGeometries << fromGeosGeom( testedGeometries[i] );
05031 }
05032
05033 GEOSGeom_destroy( polygons );
05034 return 0;
05035 }
05036
05037 GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos )
05038 {
05039
05040 int nIntersections = 0;
05041 int lastIntersectingRing = -2;
05042 const GEOSGeometry* lastIntersectingGeom = 0;
05043
05044 int nRings = GEOSGetNumInteriorRings( polygon );
05045 if ( nRings < 0 )
05046 {
05047 return 0;
05048 }
05049
05050
05051 const GEOSGeometry* outerRing = GEOSGetExteriorRing( polygon );
05052 if ( GEOSIntersects( outerRing, reshapeLineGeos ) == 1 )
05053 {
05054 ++nIntersections;
05055 lastIntersectingRing = -1;
05056 lastIntersectingGeom = outerRing;
05057 }
05058
05059
05060 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
05061 for ( int i = 0; i < nRings; ++i )
05062 {
05063 innerRings[i] = GEOSGetInteriorRingN( polygon, i );
05064 if ( GEOSIntersects( innerRings[i], reshapeLineGeos ) == 1 )
05065 {
05066 ++nIntersections;
05067 lastIntersectingRing = i;
05068 lastIntersectingGeom = innerRings[i];
05069 }
05070 }
05071
05072 if ( nIntersections != 1 )
05073 {
05074 delete [] innerRings;
05075 return 0;
05076 }
05077
05078
05079 GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos );
05080 if ( !reshapeResult )
05081 {
05082 delete [] innerRings;
05083 return 0;
05084 }
05085
05086
05087 GEOSGeometry* newRing = 0;
05088 const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq( reshapeResult );
05089 GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone( reshapeSequence );
05090
05091 GEOSGeom_destroy( reshapeResult );
05092
05093 newRing = GEOSGeom_createLinearRing( newCoordSequence );
05094 if ( !newRing )
05095 {
05096 delete [] innerRings;
05097 return 0;
05098 }
05099
05100
05101 GEOSGeometry* newOuterRing = 0;
05102 if ( lastIntersectingRing == -1 )
05103 {
05104 newOuterRing = newRing;
05105 }
05106 else
05107 {
05108 newOuterRing = GEOSGeom_clone( outerRing );
05109 }
05110
05111
05112 QList<GEOSGeometry*> ringList;
05113 if ( nRings > 0 )
05114 {
05115 GEOSGeometry* outerRingPoly = GEOSGeom_createPolygon( GEOSGeom_clone( newOuterRing ), 0, 0 );
05116 if ( outerRingPoly )
05117 {
05118 GEOSGeometry* currentRing = 0;
05119 for ( int i = 0; i < nRings; ++i )
05120 {
05121 if ( lastIntersectingRing == i )
05122 {
05123 currentRing = newRing;
05124 }
05125 else
05126 {
05127 currentRing = GEOSGeom_clone( innerRings[i] );
05128 }
05129
05130
05131 if ( GEOSContains( outerRingPoly, currentRing ) == 1 )
05132 {
05133 ringList.push_back( currentRing );
05134 }
05135 else
05136 {
05137 GEOSGeom_destroy( currentRing );
05138 }
05139 }
05140 }
05141 GEOSGeom_destroy( outerRingPoly );
05142 }
05143
05144 GEOSGeometry** newInnerRings = new GEOSGeometry*[ringList.size()];
05145 QList<GEOSGeometry*>::const_iterator it = ringList.constBegin();
05146 for ( int i = 0; i < ringList.size(); ++i )
05147 {
05148 newInnerRings[i] = ringList.at( i );
05149 }
05150
05151 delete [] innerRings;
05152
05153 GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon( newOuterRing, newInnerRings, ringList.size() );
05154 delete[] newInnerRings;
05155 if ( !reshapedPolygon )
05156 {
05157 return 0;
05158 }
05159 return reshapedPolygon;
05160 }
05161
05162 GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos )
05163 {
05164 if ( !line || !reshapeLineGeos )
05165 {
05166 return 0;
05167 }
05168
05169
05170 GEOSGeometry* intersectGeom = GEOSIntersection( line, reshapeLineGeos );
05171 bool atLeastTwoIntersections = ( GEOSGeomTypeId( intersectGeom ) == GEOS_MULTIPOINT && GEOSGetNumGeometries( intersectGeom ) > 1 );
05172 GEOSGeom_destroy( intersectGeom );
05173 if ( !atLeastTwoIntersections )
05174 {
05175 return 0;
05176 }
05177
05178
05179 const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq( line );
05180 if ( !lineCoordSeq )
05181 {
05182 return 0;
05183 }
05184 unsigned int lineCoordSeqSize;
05185 if ( GEOSCoordSeq_getSize( lineCoordSeq, &lineCoordSeqSize ) == 0 )
05186 {
05187 return 0;
05188 }
05189 if ( lineCoordSeqSize < 2 )
05190 {
05191 return 0;
05192 }
05193
05194 double x1, y1, x2, y2;
05195 GEOSCoordSeq_getX( lineCoordSeq, 0, &x1 );
05196 GEOSCoordSeq_getY( lineCoordSeq, 0, &y1 );
05197 GEOSCoordSeq_getX( lineCoordSeq, lineCoordSeqSize - 1, &x2 );
05198 GEOSCoordSeq_getY( lineCoordSeq, lineCoordSeqSize - 1, &y2 );
05199 GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) );
05200 GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) );
05201
05202 bool isRing = false;
05203 if ( GEOSGeomTypeId( line ) == GEOS_LINEARRING || GEOSEquals( beginLineVertex, endLineVertex ) == 1 )
05204 {
05205 isRing = true;
05206 }
05207
05208
05209 GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line );
05210 if ( !nodedGeometry )
05211 {
05212 GEOSGeom_destroy( beginLineVertex );
05213 GEOSGeom_destroy( endLineVertex );
05214 return 0;
05215 }
05216
05217
05218 GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
05219 GEOSGeom_destroy( nodedGeometry );
05220 if ( !mergedLines )
05221 {
05222 GEOSGeom_destroy( beginLineVertex );
05223 GEOSGeom_destroy( endLineVertex );
05224 return 0;
05225 }
05226
05227 int numMergedLines = GEOSGetNumGeometries( mergedLines );
05228 if ( numMergedLines < 2 )
05229 {
05230 GEOSGeom_destroy( beginLineVertex );
05231 GEOSGeom_destroy( endLineVertex );
05232 if ( numMergedLines == 1 )
05233 {
05234 return GEOSGeom_clone( reshapeLineGeos );
05235 }
05236 else
05237 {
05238 return 0;
05239 }
05240 }
05241
05242 QList<GEOSGeometry*> resultLineParts;
05243 QList<GEOSGeometry*> probableParts;
05244
05245 for ( int i = 0; i < numMergedLines; ++i )
05246 {
05247 const GEOSGeometry* currentGeom;
05248
05249 currentGeom = GEOSGetGeometryN( mergedLines, i );
05250 const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq( currentGeom );
05251 unsigned int currentCoordSeqSize;
05252 GEOSCoordSeq_getSize( currentCoordSeq, ¤tCoordSeqSize );
05253 if ( currentCoordSeqSize < 2 )
05254 {
05255 continue;
05256 }
05257
05258
05259 double xBegin, xEnd, yBegin, yEnd;
05260 GEOSCoordSeq_getX( currentCoordSeq, 0, &xBegin );
05261 GEOSCoordSeq_getY( currentCoordSeq, 0, &yBegin );
05262 GEOSCoordSeq_getX( currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
05263 GEOSCoordSeq_getY( currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
05264 GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) );
05265 GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) );
05266
05267
05268 int nEndpointsOnOriginalLine = 0;
05269 if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 )
05270 {
05271 nEndpointsOnOriginalLine += 1;
05272 }
05273
05274 if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 )
05275 {
05276 nEndpointsOnOriginalLine += 1;
05277 }
05278
05279
05280 int nEndpointsSameAsOriginalLine = 0;
05281 if ( GEOSEquals( beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( beginCurrentGeomVertex, endLineVertex ) == 1 )
05282 {
05283 nEndpointsSameAsOriginalLine += 1;
05284 }
05285 if ( GEOSEquals( endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( endCurrentGeomVertex, endLineVertex ) == 1 )
05286 {
05287 nEndpointsSameAsOriginalLine += 1;
05288 }
05289
05290
05291 bool currentGeomOverlapsOriginalGeom = false;
05292 bool currentGeomOverlapsReshapeLine = false;
05293 if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 )
05294 {
05295 currentGeomOverlapsOriginalGeom = true;
05296 }
05297 if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
05298 {
05299 currentGeomOverlapsReshapeLine = true;
05300 }
05301
05302
05303
05304 if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
05305 {
05306 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
05307 }
05308
05309 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
05310 {
05311 probableParts.push_back( GEOSGeom_clone( currentGeom ) );
05312 }
05313 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
05314 {
05315 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
05316 }
05317 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
05318 {
05319 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
05320 }
05321 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
05322 {
05323 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
05324 }
05325
05326 GEOSGeom_destroy( beginCurrentGeomVertex );
05327 GEOSGeom_destroy( endCurrentGeomVertex );
05328 }
05329
05330
05331 if ( isRing && probableParts.size() > 0 )
05332 {
05333 GEOSGeometry* maxGeom = 0;
05334 GEOSGeometry* currentGeom = 0;
05335 double maxLength = -DBL_MAX;
05336 double currentLength = 0;
05337 for ( int i = 0; i < probableParts.size(); ++i )
05338 {
05339 currentGeom = probableParts.at( i );
05340 GEOSLength( currentGeom, ¤tLength );
05341 if ( currentLength > maxLength )
05342 {
05343 maxLength = currentLength;
05344 GEOSGeom_destroy( maxGeom );
05345 maxGeom = currentGeom;
05346 }
05347 else
05348 {
05349 GEOSGeom_destroy( currentGeom );
05350 }
05351 }
05352 resultLineParts.push_back( maxGeom );
05353 }
05354
05355 GEOSGeom_destroy( beginLineVertex );
05356 GEOSGeom_destroy( endLineVertex );
05357 GEOSGeom_destroy( mergedLines );
05358
05359 GEOSGeometry* result = 0;
05360 if ( resultLineParts.size() < 1 )
05361 {
05362 return 0;
05363 }
05364 if ( resultLineParts.size() == 1 )
05365 {
05366 result = resultLineParts[0];
05367 }
05368 else
05369 {
05370 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
05371 for ( int i = 0; i < resultLineParts.size(); ++i )
05372 {
05373 lineArray[i] = resultLineParts[i];
05374 }
05375
05376
05377 GEOSGeometry* multiLineGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, lineArray, resultLineParts.size() );
05378 delete [] lineArray;
05379
05380
05381 result = GEOSLineMerge( multiLineGeom );
05382 GEOSGeom_destroy( multiLineGeom );
05383 }
05384
05385
05386 if ( GEOSGeomTypeId( result ) != GEOS_LINESTRING )
05387 {
05388 GEOSGeom_destroy( result );
05389 return 0;
05390 }
05391 return result;
05392 }
05393
05394 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const
05395 {
05396
05397
05398
05399
05400 testPoints.clear();
05401 GEOSGeometry* intersectionGeom = GEOSIntersection( mGeos, splitLine );
05402 if ( !intersectionGeom )
05403 {
05404 return 1;
05405 }
05406
05407 bool simple = false;
05408 int nIntersectGeoms = 1;
05409 if ( GEOSGeomTypeId( intersectionGeom ) == GEOS_LINESTRING || GEOSGeomTypeId( intersectionGeom ) == GEOS_POINT )
05410 {
05411 simple = true;
05412 }
05413
05414 if ( !simple )
05415 {
05416 nIntersectGeoms = GEOSGetNumGeometries( intersectionGeom );
05417 }
05418
05419 for ( int i = 0; i < nIntersectGeoms; ++i )
05420 {
05421 const GEOSGeometry* currentIntersectGeom;
05422 if ( simple )
05423 {
05424 currentIntersectGeom = intersectionGeom;
05425 }
05426 else
05427 {
05428 currentIntersectGeom = GEOSGetGeometryN( intersectionGeom, i );
05429 }
05430
05431 const GEOSCoordSequence* lineSequence = GEOSGeom_getCoordSeq( currentIntersectGeom );
05432 unsigned int sequenceSize = 0;
05433 double x, y;
05434 if ( GEOSCoordSeq_getSize( lineSequence, &sequenceSize ) != 0 )
05435 {
05436 for ( unsigned int i = 0; i < sequenceSize; ++i )
05437 {
05438 if ( GEOSCoordSeq_getX( lineSequence, i, &x ) != 0 )
05439 {
05440 if ( GEOSCoordSeq_getY( lineSequence, i, &y ) != 0 )
05441 {
05442 testPoints.push_back( QgsPoint( x, y ) );
05443 }
05444 }
05445 }
05446 }
05447 }
05448 GEOSGeom_destroy( intersectionGeom );
05449 return 0;
05450 }
05451
05452 GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
05453 {
05454 if ( !splitLine || !geom )
05455 {
05456 return 0;
05457 }
05458
05459 GEOSGeometry *geometryBoundary = 0;
05460 if ( GEOSGeomTypeId( geom ) == GEOS_POLYGON || GEOSGeomTypeId( geom ) == GEOS_MULTIPOLYGON )
05461 {
05462 geometryBoundary = GEOSBoundary( geom );
05463 }
05464 else
05465 {
05466 geometryBoundary = GEOSGeom_clone( geom );
05467 }
05468
05469 GEOSGeometry *splitLineClone = GEOSGeom_clone( splitLine );
05470 GEOSGeometry *unionGeometry = GEOSUnion( splitLineClone, geometryBoundary );
05471 GEOSGeom_destroy( splitLineClone );
05472
05473 GEOSGeom_destroy( geometryBoundary );
05474 return unionGeometry;
05475 }
05476
05477 int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 )
05478 {
05479 if ( !line1 || !line2 )
05480 {
05481 return -1;
05482 }
05483
05484
05485 double bufferDistance = 0.00001;
05486 if ( geomInDegrees( line2 ) )
05487 {
05488 bufferDistance = 0.00000001;
05489 }
05490 GEOSGeometry* bufferGeom = GEOSBuffer( line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS );
05491 if ( !bufferGeom )
05492 {
05493 return -2;
05494 }
05495
05496 GEOSGeometry* intersectionGeom = GEOSIntersection( bufferGeom, line1 );
05497
05498
05499 double intersectGeomLength;
05500 double line1Length;
05501
05502 GEOSLength( intersectionGeom, &intersectGeomLength );
05503 GEOSLength( line1, &line1Length );
05504
05505 GEOSGeom_destroy( bufferGeom );
05506 GEOSGeom_destroy( intersectionGeom );
05507
05508 double intersectRatio = line1Length / intersectGeomLength;
05509 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
05510 {
05511 return 1;
05512 }
05513 return 0;
05514 }
05515
05516 int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line )
05517 {
05518 if ( !point || !line )
05519 {
05520 return -1;
05521 }
05522
05523 double bufferDistance = 0.000001;
05524 if ( geomInDegrees( line ) )
05525 {
05526 bufferDistance = 0.00000001;
05527 }
05528 GEOSGeometry* lineBuffer = GEOSBuffer( line, bufferDistance, 8 );
05529 if ( !lineBuffer )
05530 {
05531 return -2;
05532 }
05533
05534 bool contained = false;
05535 if ( GEOSContains( lineBuffer, point ) == 1 )
05536 {
05537 contained = true;
05538 }
05539
05540 GEOSGeom_destroy( lineBuffer );
05541 return contained;
05542 }
05543
05544 bool QgsGeometry::geomInDegrees( const GEOSGeometry* geom )
05545 {
05546 GEOSGeometry* bbox = GEOSEnvelope( geom );
05547 if ( !bbox )
05548 {
05549 return false;
05550 }
05551
05552 const GEOSGeometry* bBoxRing = GEOSGetExteriorRing( bbox );
05553 if ( !bBoxRing )
05554 {
05555 return false;
05556 }
05557 const GEOSCoordSequence* bBoxCoordSeq = GEOSGeom_getCoordSeq( bBoxRing );
05558
05559 if ( !bBoxCoordSeq )
05560 {
05561 return false;
05562 }
05563
05564 unsigned int nCoords = 0;
05565 if ( !GEOSCoordSeq_getSize( bBoxCoordSeq, &nCoords ) )
05566 {
05567 return false;
05568 }
05569
05570 double x, y;
05571 for ( unsigned int i = 0; i < ( nCoords - 1 ); ++i )
05572 {
05573 GEOSCoordSeq_getX( bBoxCoordSeq, i, &x );
05574 if ( x > 180 || x < -180 )
05575 {
05576 return false;
05577 }
05578 GEOSCoordSeq_getY( bBoxCoordSeq, i, &y );
05579 if ( y > 90 || y < -90 )
05580 {
05581 return false;
05582 }
05583 }
05584
05585 return true;
05586 }
05587
05588 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const
05589 {
05590 if ( !g )
05591 {
05592 return 0;
05593 }
05594 int geometryType = GEOSGeomTypeId( g );
05595 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
05596 || geometryType == GEOS_POLYGON )
05597 {
05598 return 1;
05599 }
05600
05601
05602 return GEOSGetNumGeometries( g );
05603 }
05604
05605 int QgsGeometry::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry*>& splitResult )
05606 {
05607 if ( !mGeos || mDirtyGeos )
05608 {
05609 if ( !exportWkbToGeos() )
05610 return 1;
05611 }
05612
05613
05614 int type = GEOSGeomTypeId( mGeos );
05615 if ( type != GEOS_GEOMETRYCOLLECTION &&
05616 type != GEOS_MULTILINESTRING &&
05617 type != GEOS_MULTIPOLYGON &&
05618 type != GEOS_MULTIPOINT )
05619 {
05620 return 0;
05621 }
05622
05623 QVector<GEOSGeometry*> copyList = splitResult;
05624 splitResult.clear();
05625
05626
05627 QVector<GEOSGeometry*> unionGeom;
05628
05629 for ( int i = 0; i < copyList.size(); ++i )
05630 {
05631
05632 bool isPart = false;
05633 for ( int j = 0; j < GEOSGetNumGeometries( mGeos ); j++ )
05634 {
05635 if ( GEOSEquals( copyList[i], GEOSGetGeometryN( mGeos, j ) ) )
05636 {
05637 isPart = true;
05638 break;
05639 }
05640 }
05641
05642 if ( isPart )
05643 {
05644 unionGeom << copyList[i];
05645 }
05646 else
05647 {
05648 QVector<GEOSGeometry*> geomVector;
05649 geomVector << copyList[i];
05650
05651 if ( type == GEOS_MULTILINESTRING )
05652 {
05653 splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector );
05654 }
05655 else if ( type == GEOS_MULTIPOLYGON )
05656 {
05657 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
05658 }
05659 else
05660 {
05661 GEOSGeom_destroy( copyList[i] );
05662 }
05663 }
05664 }
05665
05666
05667 if ( unionGeom.size() > 0 )
05668 {
05669 if ( type == GEOS_MULTILINESTRING )
05670 {
05671 splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom );
05672 }
05673 else if ( type == GEOS_MULTIPOLYGON )
05674 {
05675 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom );
05676 }
05677 }
05678 else
05679 {
05680 unionGeom.clear();
05681 }
05682
05683 return 0;
05684 }
05685
05686 QgsPoint QgsGeometry::asPoint( unsigned char*& ptr, bool hasZValue )
05687 {
05688 ptr += 5;
05689 double* x = ( double * )( ptr );
05690 double* y = ( double * )( ptr + sizeof( double ) );
05691 ptr += 2 * sizeof( double );
05692
05693 if ( hasZValue )
05694 ptr += sizeof( double );
05695
05696 return QgsPoint( *x, *y );
05697 }
05698
05699
05700 QgsPolyline QgsGeometry::asPolyline( unsigned char*& ptr, bool hasZValue )
05701 {
05702 double x, y;
05703 ptr += 5;
05704 unsigned int nPoints = *(( int* )ptr );
05705 ptr += 4;
05706
05707 QgsPolyline line( nPoints );
05708
05709
05710 for ( uint i = 0; i < nPoints; ++i )
05711 {
05712 x = *(( double * ) ptr );
05713 y = *(( double * )( ptr + sizeof( double ) ) );
05714
05715 ptr += 2 * sizeof( double );
05716
05717 line[i] = QgsPoint( x, y );
05718
05719 if ( hasZValue )
05720 ptr += sizeof( double );
05721 }
05722
05723 return line;
05724 }
05725
05726
05727 QgsPolygon QgsGeometry::asPolygon( unsigned char*& ptr, bool hasZValue )
05728 {
05729 double x, y;
05730
05731 ptr += 5;
05732
05733
05734 unsigned int numRings = *(( int* )ptr );
05735 ptr += 4;
05736
05737 if ( numRings == 0 )
05738 return QgsPolygon();
05739
05740 QgsPolygon rings( numRings );
05741
05742 for ( uint idx = 0; idx < numRings; idx++ )
05743 {
05744 uint nPoints = *(( int* )ptr );
05745 ptr += 4;
05746
05747 QgsPolyline ring( nPoints );
05748
05749 for ( uint jdx = 0; jdx < nPoints; jdx++ )
05750 {
05751 x = *(( double * ) ptr );
05752 y = *(( double * )( ptr + sizeof( double ) ) );
05753
05754 ptr += 2 * sizeof( double );
05755
05756 if ( hasZValue )
05757 ptr += sizeof( double );
05758
05759 ring[jdx] = QgsPoint( x, y );
05760 }
05761
05762 rings[idx] = ring;
05763 }
05764
05765 return rings;
05766 }
05767
05768
05769 QgsPoint QgsGeometry::asPoint()
05770 {
05771 QGis::WkbType type = wkbType();
05772 if ( type != QGis::WKBPoint && type != QGis::WKBPoint25D )
05773 return QgsPoint( 0, 0 );
05774
05775 unsigned char* ptr = mGeometry;
05776 return asPoint( ptr, type == QGis::WKBPoint25D );
05777 }
05778
05779 QgsPolyline QgsGeometry::asPolyline()
05780 {
05781 QGis::WkbType type = wkbType();
05782 if ( type != QGis::WKBLineString && type != QGis::WKBLineString25D )
05783 return QgsPolyline();
05784
05785 unsigned char *ptr = mGeometry;
05786 return asPolyline( ptr, type == QGis::WKBLineString25D );
05787 }
05788
05789 QgsPolygon QgsGeometry::asPolygon()
05790 {
05791 QGis::WkbType type = wkbType();
05792 if ( type != QGis::WKBPolygon && type != QGis::WKBPolygon25D )
05793 return QgsPolygon();
05794
05795 unsigned char *ptr = mGeometry;
05796 return asPolygon( ptr, type == QGis::WKBPolygon25D );
05797 }
05798
05799 QgsMultiPoint QgsGeometry::asMultiPoint()
05800 {
05801 QGis::WkbType type = wkbType();
05802 if ( type != QGis::WKBMultiPoint && type != QGis::WKBMultiPoint25D )
05803 return QgsMultiPoint();
05804
05805 bool hasZValue = ( type == QGis::WKBMultiPoint25D );
05806
05807 unsigned char* ptr = mGeometry + 5;
05808 unsigned int nPoints = *(( int* )ptr );
05809 ptr += 4;
05810
05811 QgsMultiPoint points( nPoints );
05812 for ( uint i = 0; i < nPoints; i++ )
05813 {
05814 points[i] = asPoint( ptr, hasZValue );
05815 }
05816
05817 return points;
05818 }
05819
05820 QgsMultiPolyline QgsGeometry::asMultiPolyline()
05821 {
05822 QGis::WkbType type = wkbType();
05823 if ( type != QGis::WKBMultiLineString && type != QGis::WKBMultiLineString25D )
05824 return QgsMultiPolyline();
05825
05826 bool hasZValue = ( type == QGis::WKBMultiLineString25D );
05827
05828 unsigned char* ptr = mGeometry + 5;
05829 unsigned int numLineStrings = *(( int* )ptr );
05830 ptr += 4;
05831
05832 QgsMultiPolyline lines( numLineStrings );
05833
05834 for ( uint i = 0; i < numLineStrings; i++ )
05835 {
05836 lines[i] = asPolyline( ptr, hasZValue );
05837 }
05838
05839 return lines;
05840 }
05841
05842 QgsMultiPolygon QgsGeometry::asMultiPolygon()
05843 {
05844 QGis::WkbType type = wkbType();
05845 if ( type != QGis::WKBMultiPolygon && type != QGis::WKBMultiPolygon25D )
05846 return QgsMultiPolygon();
05847
05848 bool hasZValue = ( type == QGis::WKBMultiPolygon25D );
05849
05850 unsigned char* ptr = mGeometry + 5;
05851 unsigned int numPolygons = *(( int* )ptr );
05852 ptr += 4;
05853
05854 QgsMultiPolygon polygons( numPolygons );
05855
05856 for ( uint i = 0; i < numPolygons; i++ )
05857 {
05858 polygons[i] = asPolygon( ptr, hasZValue );
05859 }
05860
05861 return polygons;
05862 }
05863
05864 double QgsGeometry::area()
05865 {
05866 if ( !mGeos )
05867 {
05868 exportWkbToGeos();
05869 }
05870
05871 double area;
05872
05873 try
05874 {
05875 if ( GEOSArea( mGeos, &area ) == 0 )
05876 return -1.0;
05877 }
05878 CATCH_GEOS( -1.0 )
05879
05880 return area;
05881 }
05882
05883 double QgsGeometry::length()
05884 {
05885 if ( !mGeos )
05886 {
05887 exportWkbToGeos();
05888 }
05889
05890 double length;
05891
05892 try
05893 {
05894 if ( GEOSLength( mGeos, &length ) == 0 )
05895 return -1.0;
05896 }
05897 CATCH_GEOS( -1.0 )
05898
05899 return length;
05900 }
05901 double QgsGeometry::distance( QgsGeometry& geom )
05902 {
05903 if ( !mGeos )
05904 {
05905 exportWkbToGeos();
05906 }
05907
05908 if ( !geom.mGeos )
05909 {
05910 geom.exportWkbToGeos();
05911 }
05912
05913 if ( !mGeos || !geom.mGeos )
05914 return -1.0;
05915
05916 double dist = -1.0;
05917
05918 try
05919 {
05920 GEOSDistance( mGeos, geom.mGeos, &dist );
05921 }
05922 CATCH_GEOS( -1.0 )
05923
05924 return dist;
05925 }
05926
05927
05928 QgsGeometry* QgsGeometry::buffer( double distance, int segments )
05929 {
05930 if ( !mGeos )
05931 {
05932 exportWkbToGeos();
05933 }
05934 if ( !mGeos )
05935 {
05936 return 0;
05937 }
05938
05939 try
05940 {
05941 return fromGeosGeom( GEOSBuffer( mGeos, distance, segments ) );
05942 }
05943 CATCH_GEOS( 0 )
05944 }
05945
05946 QgsGeometry* QgsGeometry::simplify( double tolerance )
05947 {
05948 if ( !mGeos )
05949 {
05950 exportWkbToGeos();
05951 }
05952 if ( !mGeos )
05953 {
05954 return 0;
05955 }
05956 try
05957 {
05958 return fromGeosGeom( GEOSSimplify( mGeos, tolerance ) );
05959 }
05960 CATCH_GEOS( 0 )
05961 }
05962
05963 QgsGeometry* QgsGeometry::centroid()
05964 {
05965 if ( !mGeos )
05966 {
05967 exportWkbToGeos();
05968 }
05969 if ( !mGeos )
05970 {
05971 return 0;
05972 }
05973 try
05974 {
05975 return fromGeosGeom( GEOSGetCentroid( mGeos ) );
05976 }
05977 CATCH_GEOS( 0 )
05978 }
05979
05980 QgsGeometry* QgsGeometry::convexHull()
05981 {
05982 if ( !mGeos )
05983 {
05984 exportWkbToGeos();
05985 }
05986 if ( !mGeos )
05987 {
05988 return 0;
05989 }
05990
05991 try
05992 {
05993 return fromGeosGeom( GEOSConvexHull( mGeos ) );
05994 }
05995 CATCH_GEOS( 0 )
05996 }
05997
05998 QgsGeometry* QgsGeometry::intersection( QgsGeometry* geometry )
05999 {
06000 if ( !geometry )
06001 {
06002 return NULL;
06003 }
06004 if ( !mGeos )
06005 {
06006 exportWkbToGeos();
06007 }
06008 if ( !geometry->mGeos )
06009 {
06010 geometry->exportWkbToGeos();
06011 }
06012 if ( !mGeos || !geometry->mGeos )
06013 {
06014 return 0;
06015 }
06016
06017 try
06018 {
06019 return fromGeosGeom( GEOSIntersection( mGeos, geometry->mGeos ) );
06020 }
06021 CATCH_GEOS( 0 )
06022 }
06023
06024 QgsGeometry* QgsGeometry::combine( QgsGeometry* geometry )
06025 {
06026 if ( !geometry )
06027 {
06028 return NULL;
06029 }
06030 if ( !mGeos )
06031 {
06032 exportWkbToGeos();
06033 }
06034 if ( !geometry->mGeos )
06035 {
06036 geometry->exportWkbToGeos();
06037 }
06038 if ( !mGeos || !geometry->mGeos )
06039 {
06040 return 0;
06041 }
06042
06043 try
06044 {
06045 GEOSGeometry* unionGeom = GEOSUnion( mGeos, geometry->mGeos );
06046 QGis::WkbType thisGeomType = wkbType();
06047 QGis::WkbType otherGeomType = geometry->wkbType();
06048 if (( thisGeomType == QGis::WKBLineString || thisGeomType == QGis::WKBLineString25D ) \
06049 && ( otherGeomType == QGis::WKBLineString || otherGeomType == QGis::WKBLineString25D ) )
06050 {
06051 GEOSGeometry* mergedGeom = GEOSLineMerge( unionGeom );
06052 if ( mergedGeom )
06053 {
06054 GEOSGeom_destroy( unionGeom );
06055 unionGeom = mergedGeom;
06056 }
06057 }
06058 return fromGeosGeom( unionGeom );
06059 }
06060 CATCH_GEOS( new QgsGeometry( *this ) )
06061 }
06062
06063 QgsGeometry* QgsGeometry::difference( QgsGeometry* geometry )
06064 {
06065 if ( !geometry )
06066 {
06067 return NULL;
06068 }
06069 if ( !mGeos )
06070 {
06071 exportWkbToGeos();
06072 }
06073 if ( !geometry->mGeos )
06074 {
06075 geometry->exportWkbToGeos();
06076 }
06077 if ( !mGeos || !geometry->mGeos )
06078 {
06079 return 0;
06080 }
06081
06082 try
06083 {
06084 return fromGeosGeom( GEOSDifference( mGeos, geometry->mGeos ) );
06085 }
06086 CATCH_GEOS( 0 )
06087 }
06088
06089 QgsGeometry* QgsGeometry::symDifference( QgsGeometry* geometry )
06090 {
06091 if ( !geometry )
06092 {
06093 return NULL;
06094 }
06095 if ( !mGeos )
06096 {
06097 exportWkbToGeos();
06098 }
06099 if ( !geometry->mGeos )
06100 {
06101 geometry->exportWkbToGeos();
06102 }
06103 if ( !mGeos || !geometry->mGeos )
06104 {
06105 return 0;
06106 }
06107
06108 try
06109 {
06110 return fromGeosGeom( GEOSSymDifference( mGeos, geometry->mGeos ) );
06111 }
06112 CATCH_GEOS( 0 )
06113 }
06114
06115 QList<QgsGeometry*> QgsGeometry::asGeometryCollection()
06116 {
06117 if ( !mGeos )
06118 {
06119 exportWkbToGeos();
06120 if ( !mGeos )
06121 return QList<QgsGeometry*>();
06122 }
06123
06124 int type = GEOSGeomTypeId( mGeos );
06125 QgsDebugMsg( "geom type: " + QString::number( type ) );
06126
06127 QList<QgsGeometry*> geomCollection;
06128
06129 if ( type != GEOS_MULTIPOINT &&
06130 type != GEOS_MULTILINESTRING &&
06131 type != GEOS_MULTIPOLYGON &&
06132 type != GEOS_GEOMETRYCOLLECTION )
06133 {
06134
06135 geomCollection.append( new QgsGeometry( *this ) );
06136 return geomCollection;
06137 }
06138
06139 int count = GEOSGetNumGeometries( mGeos );
06140 QgsDebugMsg( "geom count: " + QString::number( count ) );
06141
06142 for ( int i = 0; i < count; ++i )
06143 {
06144 const GEOSGeometry * geometry = GEOSGetGeometryN( mGeos, i );
06145 geomCollection.append( fromGeosGeom( GEOSGeom_clone( geometry ) ) );
06146 }
06147
06148 return geomCollection;
06149 }
06150
06151
06152 bool QgsGeometry::deleteRing( int ringNum, int partNum )
06153 {
06154 if ( ringNum <= 0 || partNum < 0 )
06155 return false;
06156
06157 switch ( wkbType() )
06158 {
06159 case QGis::WKBPolygon25D:
06160 case QGis::WKBPolygon:
06161 {
06162 if ( partNum != 0 )
06163 return false;
06164
06165 QgsPolygon polygon = asPolygon();
06166 if ( ringNum >= polygon.count() )
06167 return false;
06168
06169 polygon.remove( ringNum );
06170
06171 QgsGeometry* g2 = QgsGeometry::fromPolygon( polygon );
06172 *this = *g2;
06173 delete g2;
06174 return true;
06175 }
06176
06177 case QGis::WKBMultiPolygon25D:
06178 case QGis::WKBMultiPolygon:
06179 {
06180 QgsMultiPolygon mpolygon = asMultiPolygon();
06181
06182 if ( partNum >= mpolygon.count() )
06183 return false;
06184
06185 if ( ringNum >= mpolygon[partNum].count() )
06186 return false;
06187
06188 mpolygon[partNum].remove( ringNum );
06189
06190 QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
06191 *this = *g2;
06192 delete g2;
06193 return true;
06194 }
06195
06196 default:
06197 return false;
06198 }
06199 }
06200
06201
06202 bool QgsGeometry::deletePart( int partNum )
06203 {
06204 if ( partNum < 0 )
06205 return false;
06206
06207 switch ( wkbType() )
06208 {
06209 case QGis::WKBMultiPoint25D:
06210 case QGis::WKBMultiPoint:
06211 {
06212 QgsMultiPoint mpoint = asMultiPoint();
06213
06214 if ( partNum >= mpoint.size() || mpoint.size() == 1 )
06215 return false;
06216
06217 mpoint.remove( partNum );
06218
06219 QgsGeometry* g2 = QgsGeometry::fromMultiPoint( mpoint );
06220 *this = *g2;
06221 delete g2;
06222 break;
06223 }
06224
06225 case QGis::WKBMultiLineString25D:
06226 case QGis::WKBMultiLineString:
06227 {
06228 QgsMultiPolyline mline = asMultiPolyline();
06229
06230 if ( partNum >= mline.size() || mline.size() == 1 )
06231 return false;
06232
06233 mline.remove( partNum );
06234
06235 QgsGeometry* g2 = QgsGeometry::fromMultiPolyline( mline );
06236 *this = *g2;
06237 delete g2;
06238 break;
06239 }
06240
06241 case QGis::WKBMultiPolygon25D:
06242 case QGis::WKBMultiPolygon:
06243 {
06244 QgsMultiPolygon mpolygon = asMultiPolygon();
06245
06246 if ( partNum >= mpolygon.size() || mpolygon.size() == 1 )
06247 return false;
06248
06249 mpolygon.remove( partNum );
06250
06251 QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
06252 *this = *g2;
06253 delete g2;
06254 break;
06255 }
06256
06257 default:
06258
06259 return false;
06260 }
06261
06262 return true;
06263 }
06264
06265 int QgsGeometry::avoidIntersections()
06266 {
06267 int returnValue = 0;
06268
06269
06270 if ( type() != QGis::Polygon )
06271 {
06272 return 1;
06273 }
06274
06275 QGis::WkbType geomTypeBeforeModification = wkbType();
06276
06277
06278 bool listReadOk;
06279 QStringList avoidIntersectionsList = QgsProject::instance()->readListEntry( "Digitizing", "/AvoidIntersectionsList", &listReadOk );
06280 if ( !listReadOk )
06281 {
06282 return true;
06283 }
06284
06285
06286 QgsVectorLayer* currentLayer = 0;
06287 QStringList::const_iterator aIt = avoidIntersectionsList.constBegin();
06288 for ( ; aIt != avoidIntersectionsList.constEnd(); ++aIt )
06289 {
06290 currentLayer = dynamic_cast<QgsVectorLayer*>( QgsMapLayerRegistry::instance()->mapLayer( *aIt ) );
06291 if ( currentLayer )
06292 {
06293 if ( currentLayer->removePolygonIntersections( this ) != 0 )
06294 {
06295 returnValue = 3;
06296 }
06297 }
06298 }
06299
06300
06301 if ( wkbType() != geomTypeBeforeModification )
06302 {
06303 return 2;
06304 }
06305
06306 return returnValue;
06307 }
06308
06309
06310
06311
06312
06313
06314 static double distLine2Point( QgsPoint p, QgsVector v, QgsPoint q )
06315 {
06316 if ( v.length() == 0 )
06317 {
06318 throw QgsException( QObject::tr( "invalid line" ) );
06319 }
06320
06321 return ( v.x()*( q.y() - p.y() ) - v.y()*( q.x() - p.x() ) ) / v.length();
06322 }
06323
06324 static bool intersectLines( QgsPoint p, QgsVector v, QgsPoint q, QgsVector w, QgsPoint &s )
06325 {
06326 double d = v.y() * w.x() - v.x() * w.y();
06327
06328 if ( d == 0 )
06329 return false;
06330
06331 double dx = q.x() - p.x();
06332 double dy = q.y() - p.y();
06333 double k = ( dy * w.x() - dx * w.y() ) / d;
06334
06335 s = p + v * k;
06336
06337 return true;
06338 }
06339
06340 bool pointInRing( const QgsPolyline &ring, const QgsPoint &p )
06341 {
06342 bool inside = false;
06343 int j = ring.size() - 1;
06344
06345 for ( int i = 0; i < ring.size(); i++ )
06346 {
06347 if ( ring[i].x() == p.x() && ring[i].y() == p.y() )
06348 return true;
06349
06350 if (( ring[i].y() < p.y() && ring[j].y() >= p.y() ) ||
06351 ( ring[j].y() < p.y() && ring[i].y() >= p.y() ) )
06352 {
06353 if ( ring[i].x() + ( p.y() - ring[i].y() ) / ( ring[j].y() - ring[i].y() )*( ring[j].x() - ring[i].x() ) <= p.x() )
06354 inside = !inside;
06355 }
06356
06357 j = i;
06358 }
06359
06360 return inside;
06361 }
06362
06363 static bool ringInRing( const QgsPolyline &inside, const QgsPolyline &outside )
06364 {
06365 for ( int i = 0; i < inside.size(); i++ )
06366 {
06367 if ( !pointInRing( outside, inside[i] ) )
06368 return false;
06369 }
06370
06371 return true;
06372 }
06373
06374 void QgsGeometry::checkRingIntersections( QList<Error> &errors,
06375 int p0, int i0, const QgsPolyline &ring0,
06376 int p1, int i1, const QgsPolyline &ring1 )
06377 {
06378 for ( int i = 0; i < ring0.size() - 1; i++ )
06379 {
06380 QgsVector v = ring0[i+1] - ring0[i];
06381
06382 for ( int j = 0; j < ring1.size() - 1; j++ )
06383 {
06384 QgsVector w = ring1[j+1] - ring1[j];
06385
06386 QgsPoint s;
06387 if ( intersectLines( ring0[i], v, ring1[j], w, s ) )
06388 {
06389 double d = -distLine2Point( ring0[i], v.perpVector(), s );
06390
06391 if ( d >= 0 && d <= v.length() )
06392 {
06393 d = -distLine2Point( ring1[j], w.perpVector(), s );
06394 if ( d >= 0 && d <= w.length() )
06395 {
06396 QString msg = QObject::tr( "segment %1 of ring %2 of polygon %3 intersects segment %4 of ring %5 of polygon %6 at %7" )
06397 .arg( i0 ).arg( i ).arg( p0 )
06398 .arg( i1 ).arg( j ).arg( p1 )
06399 .arg( s.toString() );
06400 QgsDebugMsg( msg );
06401 errors << Error( msg, s );
06402 if ( errors.size() > 100 )
06403 {
06404 QString msg = QObject::tr( "stopping validation after more than 100 errors" );
06405 QgsDebugMsg( msg );
06406 errors << Error( msg );
06407 return;
06408 }
06409 }
06410 }
06411 }
06412 }
06413 }
06414 }
06415
06416 void QgsGeometry::validatePolyline( QList<Error> &errors, int i, QgsPolyline line, bool ring )
06417 {
06418 if ( ring )
06419 {
06420 if ( line.size() < 3 )
06421 {
06422 QString msg = QObject::tr( "ring %1 with less than three points" ).arg( i );
06423 QgsDebugMsg( msg );
06424 errors << Error( msg );
06425 return;
06426 }
06427
06428 if ( line[0] != line[ line.size()-1 ] )
06429 {
06430 QString msg = QObject::tr( "ring %1 not closed" ).arg( i );
06431 QgsDebugMsg( msg );
06432 errors << Error( msg );
06433 return;
06434 }
06435 }
06436 else if ( line.size() < 2 )
06437 {
06438 QString msg = QObject::tr( "line %1 with less than two points" ).arg( i );
06439 QgsDebugMsg( msg );
06440 errors << Error( msg );
06441 return;
06442 }
06443
06444 int j = 0;
06445 while ( j < line.size() - 1 )
06446 {
06447 int n = 0;
06448 while ( j < line.size() - 1 && line[j] == line[j+1] )
06449 {
06450 line.remove( j );
06451 n++;
06452 }
06453
06454 if ( n > 0 )
06455 {
06456 QString msg = QObject::tr( "line %1 contains %n duplicate node(s) at %2", "number of duplicate nodes", n ).arg( i ).arg( j );
06457 QgsDebugMsg( msg );
06458 errors << Error( msg, line[j] );
06459 }
06460
06461 j++;
06462 }
06463
06464 for ( j = 0; j < line.size() - 3; j++ )
06465 {
06466 QgsVector v = line[j+1] - line[j];
06467 double vl = v.length();
06468
06469 int n = ( j == 0 && ring ) ? line.size() - 2 : line.size() - 1;
06470
06471 for ( int k = j + 2; k < n; k++ )
06472 {
06473 QgsVector w = line[k+1] - line[k];
06474
06475 QgsPoint s;
06476 if ( !intersectLines( line[j], v, line[k], w, s ) )
06477 continue;
06478
06479 double d = -distLine2Point( line[j], v.perpVector(), s );
06480 if ( d < 0 || d > vl )
06481 continue;
06482
06483 d = -distLine2Point( line[k], w.perpVector(), s );
06484 if ( d < 0 || d > w.length() )
06485 continue;
06486
06487 QString msg = QObject::tr( "segments %1 and %2 of line %3 intersect at %4" ).arg( j ).arg( k ).arg( i ).arg( s.toString() );
06488 QgsDebugMsg( msg );
06489 errors << Error( msg, s );
06490 if ( errors.size() > 100 )
06491 {
06492 QString msg = QObject::tr( "stopping validation after more than 100 errors" );
06493 QgsDebugMsg( msg );
06494 errors << Error( msg );
06495 return;
06496 }
06497 }
06498 }
06499 }
06500
06501 void QgsGeometry::validatePolygon( QList<Error> &errors, int idx, const QgsPolygon &polygon )
06502 {
06503
06504 for ( int i = 1; i < polygon.size(); i++ )
06505 {
06506 if ( !ringInRing( polygon[i], polygon[0] ) )
06507 {
06508 QString msg = QObject::tr( "ring %1 of polygon %2 not in exterior ring" ).arg( i ).arg( idx );
06509 QgsDebugMsg( msg );
06510 errors << Error( msg );
06511 }
06512 }
06513
06514
06515 for ( int i = 1; i < polygon.size(); i++ )
06516 {
06517 for ( int j = i + 1; j < polygon.size(); j++ )
06518 {
06519 checkRingIntersections( errors, idx, i, polygon[i], idx, j, polygon[j] );
06520 }
06521 }
06522
06523
06524 for ( int i = 0; i < polygon.size(); i++ )
06525 {
06526 validatePolyline( errors, i, polygon[i], true );
06527 }
06528 }
06529
06530 void QgsGeometry::validateGeometry( QList<Error> &errors )
06531 {
06532 errors.clear();
06533
06534 switch ( wkbType() )
06535 {
06536 case QGis::WKBPoint:
06537 case QGis::WKBPoint25D:
06538 case QGis::WKBMultiPoint:
06539 case QGis::WKBMultiPoint25D:
06540 break;
06541
06542 case QGis::WKBLineString:
06543 case QGis::WKBLineString25D:
06544 validatePolyline( errors, 0, asPolyline() );
06545 break;
06546
06547 case QGis::WKBMultiLineString:
06548 case QGis::WKBMultiLineString25D:
06549 {
06550 QgsMultiPolyline mp = asMultiPolyline();
06551 for ( int i = 0; i < mp.size(); i++ )
06552 validatePolyline( errors, i, mp[i] );
06553 }
06554 break;
06555
06556 case QGis::WKBPolygon:
06557 case QGis::WKBPolygon25D:
06558 {
06559 validatePolygon( errors, 0, asPolygon() );
06560 }
06561 break;
06562
06563 case QGis::WKBMultiPolygon:
06564 case QGis::WKBMultiPolygon25D:
06565 {
06566 QgsMultiPolygon mp = asMultiPolygon();
06567 for ( int i = 0; i < mp.size(); i++ )
06568 {
06569 validatePolygon( errors, i, mp[i] );
06570 }
06571
06572 for ( int i = 0; i < mp.size(); i++ )
06573 {
06574 for ( int j = i + 1; j < mp.size(); j++ )
06575 {
06576 if ( ringInRing( mp[i][0], mp[j][0] ) )
06577 {
06578 errors << Error( QObject::tr( "polygon %1 inside polygon %2" ).arg( i ).arg( j ) );
06579 }
06580 else if ( ringInRing( mp[j][0], mp[i][0] ) )
06581 {
06582 errors << Error( QObject::tr( "polygon %1 inside polygon %2" ).arg( j ).arg( i ) );
06583 }
06584 else
06585 {
06586 checkRingIntersections( errors, i, 0, mp[i][0], j, 0, mp[j][0] );
06587 }
06588 }
06589 }
06590 }
06591 break;
06592
06593 case QGis::WKBNoGeometry:
06594 case QGis::WKBUnknown:
06595 QgsDebugMsg( QObject::tr( "Unknown geometry type" ) );
06596 errors << Error( QObject::tr( "Unknown geometry type" ) );
06597 break;
06598 }
06599 }
06600
06601 bool QgsGeometry::isGeosValid()
06602 {
06603 try
06604 {
06605 GEOSGeometry *g = asGeos();
06606
06607 if ( !g )
06608 return false;
06609
06610 return GEOSisValid( g );
06611 }
06612 catch ( GEOSException &e )
06613 {
06614
06615 Q_UNUSED( e );
06616 QgsDebugMsg( QString( "GEOS exception caught: %1" ).arg( e.what() ) );
06617 return false;
06618 }
06619 }
06620
06621 bool QgsGeometry::isGeosEqual( QgsGeometry &g )
06622 {
06623 return geosRelOp( GEOSEquals, this, &g );
06624 }
06625
06626 bool QgsGeometry::isGeosEmpty()
06627 {
06628 try
06629 {
06630 GEOSGeometry *g = asGeos();
06631
06632 if ( !g )
06633 return false;
06634
06635 return GEOSisEmpty( g );
06636 }
06637 catch ( GEOSException &e )
06638 {
06639
06640 Q_UNUSED( e );
06641 QgsDebugMsg( QString( "GEOS exception caught: %1" ).arg( e.what() ) );
06642 return false;
06643 }
06644 }