[QGIS Commit] r11941 - trunk/qgis/src/analysis/interpolation

[svn_qgis at osgeo.org]

Author: mhugent
Date: 2009-11-06 11:14:42 -0500 (Fri, 06 Nov 2009)
New Revision: 11941

Modified:
   trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc
   trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp
Log:
Fix crash with constrained triangulation. Fixes bug #2005, but some numerical problems with constrained delaunay triangulations remain

Modified: trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc
===================================================================
--- trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc	2009-11-06 15:29:53 UTC (rev 11940)
+++ trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc	2009-11-06 16:14:42 UTC (rev 11941)
@@ -959,6 +959,11 @@
 
 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 );
@@ -1072,6 +1077,11 @@
 
 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 );

Modified: trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp	2009-11-06 15:29:53 UTC (rev 11940)
+++ trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp	2009-11-06 16:14:42 UTC (rev 11941)
@@ -167,10 +167,12 @@
     }
   }
 
-  //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 )
@@ -199,12 +201,36 @@
       }
       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();
@@ -248,8 +274,165 @@
       }
       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;
   }