Precise intersection test in ogr provider, made geos export of QgsGeometry 25D aware, fixed a memory leak

git-svn-id: http://svn.osgeo.org/qgis/trunk@6741 c8812cc2-4d05-0410-92ff-de0c093fc19c

src/core/qgsgeometry.cpp

@@ -157,7 +157,6 @@ void QgsGeometry::setWkbAndOwnership(unsigned char * wkb, size_t length)
   mDirtyWkb   = FALSE;
   mDirtyGeos  = TRUE;
   mDirtyWkt   = TRUE;
-
 }
     
 unsigned char * QgsGeometry::wkbBuffer() const
@@ -2190,7 +2189,7 @@ bool QgsGeometry::fast_intersects(const QgsRect& r) const
     rectwkt+="))";
     GEOS_GEOM::GeometryFactory *gf = new GEOS_GEOM::GeometryFactory();
     GEOS_IO::WKTReader *wktReader = new GEOS_IO::WKTReader(gf);
-    GEOS_GEOM::Geometry *geosRect = wktReader->read( qstrdup(rectwkt) );
+    GEOS_GEOM::Geometry *geosRect = wktReader->read(rectwkt.toLocal8Bit().data());
     
     try // geos might throw exception on error
     {
@@ -2209,7 +2208,6 @@ bool QgsGeometry::fast_intersects(const QgsRect& r) const
       QgsLogger::warning("GEOS: " + error);
     }
       
-    delete geosGeom;
     delete geosRect;
     delete gf;
     delete wktReader;
@@ -2501,8 +2499,6 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
     {
       // No need to convert again
       return mGeos;
-      
-      // TODO: make mGeos useful - assign to it and clear mDirty before we return out of this function
     }
 
     if(!mGeometry)
@@ -2522,20 +2518,27 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
     unsigned char *ptr;
     char lsb;
     QgsPoint pt;
-    int wkbtype;
+    QGis::WKBTYPE wkbtype;
+    bool hasZValue = false;
 
-    wkbtype = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4];
+    //wkbtype = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4];//MH: this does not work for 2.5D types
+    wkbtype = wkbType();
     switch(wkbtype)
-    {
+      {
-	case QGis::WKBPoint:
+      case QGis::WKBPoint25D:
+      case QGis::WKBPoint:
 	{
 	   x = (double *) (mGeometry + 5);
 	   y = (double *) (mGeometry + 5 + sizeof(double));
            
-           mDirtyGeos = FALSE;
+	   mDirtyGeos = FALSE;
-	   return geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(*x,*y));
+           mGeos = geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(*x,*y));
+	   mDirtyGeos = FALSE;
+	   return mGeos;
 	}
-	case QGis::WKBMultiPoint:
+      case QGis::WKBMultiPoint25D:
+	hasZValue = true;
+      case QGis::WKBMultiPoint:
 	{
 	    std::vector<GEOS_GEOM::Geometry*>* points=new std::vector<GEOS_GEOM::Geometry*>;
 	    ptr = mGeometry + 5;
@@ -2549,10 +2552,19 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 		y = (double *) ptr;
 		ptr += sizeof(double);
 		points->push_back(geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(*x,*y)));
+		if(hasZValue)
+		  {
+		    ptr += sizeof(double);
+		  }
 	    }
-	    return geosGeometryFactory->createMultiPoint(points);
+	    delete mGeos;
+	    mGeos = geosGeometryFactory->createMultiPoint(points);
+	    mDirtyGeos = FALSE;
+	    return mGeos;
 	}
-	case QGis::WKBLineString:
+      case QGis::WKBLineString25D:
+	hasZValue = true;
+      case QGis::WKBLineString:
 	{
       QgsDebugMsg("QgsGeometry::geosGeometry: Linestring found");
 
@@ -2567,10 +2579,19 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 		y = (double *) ptr;
 		ptr += sizeof(double);
 		sequence->add(GEOS_GEOM::Coordinate(*x,*y));
+		if(hasZValue)
+		  {
+		    ptr += sizeof(double);
+		  }
 	    }
-	    return geosGeometryFactory->createLineString(sequence); 
+	    delete mGeos;
+	    mGeos = geosGeometryFactory->createLineString(sequence);
+	    mDirtyGeos = FALSE;
+	    return mGeos;
 	}
-	case QGis::WKBMultiLineString:
+      case QGis::WKBMultiLineString25D:
+	hasZValue = true;
+      case QGis::WKBMultiLineString:
 	{
 	    std::vector<GEOS_GEOM::Geometry*>* lines=new std::vector<GEOS_GEOM::Geometry*>;
 	    numLineStrings = (int) (mGeometry[5]);
@@ -2590,14 +2611,23 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 		    y = (double *) ptr;
 		    ptr += sizeof(double);
 		    sequence->add(GEOS_GEOM::Coordinate(*x,*y));
+		    if(hasZValue)
+		      {
+			ptr += sizeof(double);
+		      }
 		}
 		lines->push_back(geosGeometryFactory->createLineString(sequence));
 	    }
-	    return geosGeometryFactory->createMultiLineString(lines);
+	    delete mGeos;
+	    mGeos = geosGeometryFactory->createMultiLineString(lines);
+	    mDirtyGeos = FALSE;
+	    return mGeos;
 	}
-	case QGis::WKBPolygon: 
+      case QGis::WKBPolygon25D:
+	hasZValue = true;
+      case QGis::WKBPolygon: 
 	{
-      QgsDebugMsg("Polygon found");
+	  QgsDebugMsg("Polygon found");
 
 	    // get number of rings in the polygon
 	    numRings = (int *) (mGeometry + 1 + sizeof(int));
@@ -2623,6 +2653,10 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 		    y = (double *) ptr;
 		    ptr += sizeof(double);
 		    sequence->add(GEOS_GEOM::Coordinate(*x,*y));
+		    if(hasZValue)
+		      {
+			ptr += sizeof(double);
+		      }
 		}
 		GEOS_GEOM::LinearRing* ring=geosGeometryFactory->createLinearRing(sequence);
 		if(idx==0)
@@ -2634,10 +2668,14 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 		    inner->push_back(ring);
 		}
 	    }
-	    return geosGeometryFactory->createPolygon(outer,inner);
+	    delete mGeos;
+	    mGeos = geosGeometryFactory->createPolygon(outer,inner);
+	    mDirtyGeos = FALSE;
+	    return mGeos;
 	}
-	
+      case QGis::WKBMultiPolygon25D:
-	case QGis::WKBMultiPolygon:
+	hasZValue = true;
+      case QGis::WKBMultiPolygon:
 	{
 	    QgsDebugMsg("Multipolygon found");
 
@@ -2675,6 +2713,10 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 			y = (double *) ptr;
 			ptr += sizeof(double);
 			sequence->add(GEOS_GEOM::Coordinate(*x,*y));
+			if(hasZValue)
+			  {
+			    ptr += sizeof(double);
+			  }
 		    }
 		    GEOS_GEOM::LinearRing* ring=geosGeometryFactory->createLinearRing(sequence);
 		    if(idx==0)
@@ -2689,10 +2731,13 @@ GEOS_GEOM::Geometry* QgsGeometry::geosGeometry() const
 	    
 		polygons->push_back(geosGeometryFactory->createPolygon(outer,inner));
 	    }
-	    return (geosGeometryFactory->createMultiPolygon(polygons));
+	    delete mGeos;
+	    mGeos = (geosGeometryFactory->createMultiPolygon(polygons));
+	    mDirtyGeos = false;
+	    return mGeos;
 	}
-	default:
+    default:
-	    return 0;
+      return 0;
     }
     
 }

src/core/qgsgeometry.h

@@ -237,7 +237,8 @@ class CORE_EXPORT QgsGeometry {
     /** Test for containment of a point (uses GEOS) */
     bool contains(QgsPoint* p) const;
 
-    /**Creates a geos geometry from this features geometry. Note, that the returned object needs to be deleted*/
+    /**Returns a geos geometry from the features geometry. Note that object pointer is cached (mGoes) and
+     should _not_ be deleted by the calling method*/
     GEOS_GEOM::Geometry* geosGeometry() const;
 
 

src/providers/ogr/qgsogrprovider.cpp

@@ -175,6 +175,7 @@ QgsOgrProvider::~QgsOgrProvider()
   extent_ = 0;
   delete geometryFactory;
   delete wktReader;
+  delete mSelectionRectangle;
 }
 
 void QgsOgrProvider::setEncoding(const QString& e)
@@ -245,90 +246,58 @@ bool QgsOgrProvider::getNextFeature(QgsFeature& feature,
   }
 
   OGRFeature *fet;
+  QgsRect selectionRect;
 
-  // skip features without geometry
   while ((fet = ogrLayer->GetNextFeature()) != NULL)
   {
+    // skip features without geometry
     if (fet->GetGeometryRef() != NULL || mFetchFeaturesWithoutGeom)
     {
+      OGRFeatureDefn * featureDefinition = fet->GetDefnRef();
+      QString featureTypeName = featureDefinition ? QString(featureDefinition->GetName()) : QString("");
+      feature.setFeatureId(fet->GetFID());
+      feature.setTypeName(featureTypeName);
+
+      if (fetchGeometry)
+	{
+	  OGRGeometry *geom = fet->GetGeometryRef();
+    
+	  // get the wkb representation
+	  unsigned char *wkb = new unsigned char[geom->WkbSize()];
+	  geom->exportToWkb((OGRwkbByteOrder) QgsApplication::endian(), wkb);
+	  
+	  feature.setGeometryAndOwnership(wkb, geom->WkbSize());
+
+	  if(mUseIntersect)
+	    {
+	      //precise test for intersection with search rectangle
+	      //first make QgsRect from OGRPolygon
+	      OGREnvelope env;
+	      mSelectionRectangle->getEnvelope(&env);
+	      if(env.IsInit()) //if envelope is invalid, skip the precise intersection test
+		{
+		  selectionRect.set(env.MinX, env.MinY, env.MaxX, env.MaxY);
+		  if(!feature.geometry()->fast_intersects(selectionRect))
+		    {
+		      delete fet;
+		      continue;
+		    }
+		}
+	    }
+	}
+
+      for(QgsAttributeList::iterator it = fetchAttributes.begin(); it != fetchAttributes.end(); ++it)
+	{
+	  getFeatureAttribute(fet,feature,*it);
+	}
+    
       break;
     }
-      
-    delete fet;
-  
-    /* TODO: [MD]
-    if (mUseIntersect)
-    {
-      GEOS_GEOM::Geometry *geosGeom = 0;
-      geosGeom = f->geometry()->geosGeometry();
-      assert(geosGeom != 0);
-         
-      char *sWkt = new char[2 * mSelectionRectangle->WkbSize()];
-      mSelectionRectangle->exportToWkt(&sWkt);  
-      GEOS_GEOM::Geometry *geosRect = wktReader->read(sWkt);
-      assert(geosRect != 0);
-      try // geos might throw exception on error 
-      { 
-        if(geosGeom->intersects(geosRect)) 
-        { 
-          returnval=true; 
-        } 
-      } 
-      catch (GEOS_UTIL::TopologyException* e) 
-      { 
-#if GEOS_VERSION_MAJOR < 3
-        QString error = e->toString().c_str();
-#else
-        QString error = e->what();
-#endif
-        QString error = e->toString().c_str(); 
-        QgsLogger::warning("GEOS: " + error);
-      }
-      
-      if (returnval)
-      {
-        QgsDebugMsg("intersection found");
-        delete[] sWkt;
-        delete geosGeom;
-        break;
-      }
-      else
-      {
-        QgsDebugMsg("no intersection found");
-        delete[] sWkt;
-        delete geosGeom;
-        delete f;
-        f=0;
-      }
-    }*/
   }
-  
+   
   if (fet)
   {
-    // get type name
-    OGRFeatureDefn * featureDefinition = fet->GetDefnRef();
-    QString featureTypeName =   
-        featureDefinition ? QString(featureDefinition->GetName()) : QString("");
-    
-    feature = QgsFeature(fet->GetFID(), featureTypeName);
-
-    if (fetchGeometry)
-    {
-      OGRGeometry *geom = fet->GetGeometryRef();
-    
-      // get the wkb representation
-      unsigned char *wkb = new unsigned char[geom->WkbSize()];
-      geom->exportToWkb((OGRwkbByteOrder) QgsApplication::endian(), wkb);
-    
-      feature.setGeometryAndOwnership(wkb, geom->WkbSize());
-    }
-
-    for(QgsAttributeList::iterator it = fetchAttributes.begin(); it != fetchAttributes.end(); ++it)
-    {
-      getFeatureAttribute(fet,feature,*it);
-    }
     delete fet;
-
     return true;
   }
   else