Skip to content

Commit

Permalink
Fix crash with constrained triangulation. Fixes bug #2005, but some n…
Browse files Browse the repository at this point in the history
…umerical problems with constrained delaunay triangulations remain

git-svn-id: http://svn.osgeo.org/qgis/trunk@11941 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mhugent committed Nov 6, 2009
1 parent 4d7a201 commit 05a5712
Show file tree
Hide file tree
Showing 2 changed files with 197 additions and 4 deletions.
10 changes: 10 additions & 0 deletions src/analysis/interpolation/DualEdgeTriangulation.cc
Original file line number Diff line number Diff line change
Expand Up @@ -959,6 +959,11 @@ QList<int>* DualEdgeTriangulation::getSurroundingTriangles( int pointno )

bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, int* n1, Point3D* p2, int* n2, Point3D* p3, int* n3 )
{
if ( mPointVector.size() < 3 )
{
return false;
}

if ( p1 && p2 && p3 )
{
Point3D point( x, y, 0 );
Expand Down Expand Up @@ -1072,6 +1077,11 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, int* n

bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, Point3D* p2, Point3D* p3 )
{
if ( mPointVector.size() < 3 )
{
return false;
}

if ( p1 && p2 && p3 )
{
Point3D point( x, y, 0 );
Expand Down
191 changes: 187 additions & 4 deletions src/analysis/interpolation/qgstininterpolator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -167,10 +167,12 @@ int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputT
}
}

//parse WKB. We cannot use the methods with QgsPoint because they don't contain z-values for 25D types
//parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types
bool hasZValue = false;
double x, y, z;
unsigned char* currentWkbPtr = g->asWkb();
//maybe a structure or break line
Line3D* line = 0;

QGis::WkbType wkbType = g->wkbType();
switch ( wkbType )
Expand Down Expand Up @@ -199,12 +201,36 @@ int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputT
}
break;
}
case QGis::WKBMultiPoint25D:
hasZValue = true;
case QGis::WKBMultiPoint:
{
currentWkbPtr += ( 1 + sizeof( int ) );
int* npoints = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
for ( int index = 0; index < *npoints; ++index )
{
currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and point type
x = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
y = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
if ( hasZValue ) //skip z-coordinate for 25D geometries
{
z = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
}
else
{
z = attributeValue;
}
}
break;
}
case QGis::WKBLineString25D:
hasZValue = true;
case QGis::WKBLineString:
{
//maybe a structure or break line
Line3D* line = 0;
if ( type != POINTS )
{
line = new Line3D();
Expand Down Expand Up @@ -248,8 +274,165 @@ int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputT
}
break;
}
case QGis::WKBMultiLineString25D:
hasZValue = true;
case QGis::WKBMultiLineString:
{
currentWkbPtr += ( 1 + sizeof( int ) );
int* nlines = ( int* )currentWkbPtr;
int* npoints = 0;
currentWkbPtr += sizeof( int );
for ( int index = 0; index < *nlines; ++index )
{
if ( type != POINTS )
{
line = new Line3D();
}
currentWkbPtr += ( sizeof( int ) + 1 );
npoints = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
for ( int index2 = 0; index2 < *npoints; ++index2 )
{
x = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
y = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );

if ( hasZValue ) //skip z-coordinate for 25D geometries
{
z = *(( double* ) currentWkbPtr );
currentWkbPtr += sizeof( double );
}
else
{
z = attributeValue;
}

if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
}
}
if ( type != POINTS )
{
mTriangulation->addLine( line, type == BREAK_LINES );
}
}
break;
}
case QGis::WKBPolygon25D:
hasZValue = true;
case QGis::WKBPolygon:
{
currentWkbPtr += ( 1 + sizeof( int ) );
int* nrings = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
int* npoints;
for ( int index = 0; index < *nrings; ++index )
{
if ( type != POINTS )
{
line = new Line3D();
}

npoints = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
for ( int index2 = 0; index2 < *npoints; ++index2 )
{
x = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
y = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
if ( hasZValue ) //skip z-coordinate for 25D geometries
{
z = *(( double* )currentWkbPtr );;
currentWkbPtr += sizeof( double );
}
else
{
z = attributeValue;
}
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
}
}

if ( type != POINTS )
{
mTriangulation->addLine( line, type == BREAK_LINES );
}
}
break;
}

case QGis::WKBMultiPolygon25D:
hasZValue = true;
case QGis::WKBMultiPolygon:
{
currentWkbPtr += ( 1 + sizeof( int ) );
int* npolys = ( int* )currentWkbPtr;
int* nrings;
int* npoints;
currentWkbPtr += sizeof( int );
for ( int index = 0; index < *npolys; ++index )
{
currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and polygon type
nrings = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
for ( int index2 = 0; index2 < *nrings; ++index2 )
{
if ( type != POINTS )
{
line = new Line3D();
}
npoints = ( int* )currentWkbPtr;
currentWkbPtr += sizeof( int );
for ( int index3 = 0; index3 < *npoints; ++index3 )
{
x = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
y = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
if ( hasZValue ) //skip z-coordinate for 25D geometries
{
z = *(( double* )currentWkbPtr );
currentWkbPtr += sizeof( double );
}
else
{
z = attributeValue;
}
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
}
}
if ( type != POINTS )
{
mTriangulation->addLine( line, type == BREAK_LINES );
}
}
}
break;
}
default:
//todo: add the same for multiline, polygon, multipolygon
//should not happen...
break;
}

Expand Down

0 comments on commit 05a5712

Please sign in to comment.