22.01.2016, 12:48, Kai Muehlbauer kirjoitti:
Am 22.01.2016 um 10:52 schrieb Ari Jolma:
22.01.2016, 10:06, Kai Muehlbauer kirjoitti:
To quickly check if there is an improvement I used the python shapely
package. Converting the outer loop geometry to an shapely prepared
geometry and calling 'contains()' with the inner loop geometries to
find those fully contained geometries I get an improvement of nearly
40% to the former ogr only code.
As Even pointed out, there is only OGRPreparedGeometryIntersects()
(GEOSPreparedIntersects) exposed within ogrgeometry.cpp. In my use
case I would need GEOSPreparedContains to mimic the shapely behaviour.
I made a small test with the attached diff, which prepares the
geometries of the input layer.
The input layer had two polygons (boundaries of two Finnish
municipalities) and the method layer had 3091 points, all in either one
of the polygons.
The improvement is ~15%.
I think the prepared geometries could be used by default easily.
Prepared contains test would need addition of the respective GEOS
method.
Thanks Ari, I very much appreciate your efforts in testing layer
intersection. If the prepared contains test could be done within the
layer intersection then those features can "just be copied" without
the need of costly intersection calculation. This would give a real
speedup.
It seems so, I added the contains test and in this case the total
speedup is ~80 %.
So, it would need the addition of the OGRPreparedGeometryContains
method (GEOSPreparedContains) and then enhancing your patch.
Would it work the way as in the attached diff? Or did I miss something?
Attached is my diff, which basically copies to add-to-result code, which
might be done better.
What are the next steps to get this feature fully implemented?
The thing is that most of the other layer methods could be enhanced this
way, which means some more work. I believe the C++ API can be enhanced
for version 2.1.
Ari
Cheers,
Kai
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Index: ogr/ogr_geometry.h
===================================================================
--- ogr/ogr_geometry.h (revision 33064)
+++ ogr/ogr_geometry.h (working copy)
@@ -1343,5 +1343,7 @@
void OGRDestroyPreparedGeometry( OGRPreparedGeometry* poPreparedGeom );
int OGRPreparedGeometryIntersects( const OGRPreparedGeometry* poPreparedGeom,
const OGRGeometry* poOtherGeom );
+int OGRPreparedGeometryContains( const OGRPreparedGeometry* poPreparedGeom,
+ const OGRGeometry* poOtherGeom );
#endif /* ndef OGR_GEOMETRY_H_INCLUDED */
Index: ogr/ogrgeometry.cpp
===================================================================
--- ogr/ogrgeometry.cpp (revision 33064)
+++ ogr/ogrgeometry.cpp (working copy)
@@ -4765,6 +4765,28 @@
#endif
}
+int OGRPreparedGeometryContains( UNUSED_IF_NO_GEOS const OGRPreparedGeometry* poPreparedGeom,
+ UNUSED_IF_NO_GEOS const OGRGeometry* poOtherGeom )
+{
+#ifdef HAVE_GEOS_PREPARED_GEOMETRY
+ if( poPreparedGeom == NULL || poOtherGeom == NULL )
+ return FALSE;
+
+ GEOSGeom hGEOSOtherGeom = poOtherGeom->exportToGEOS(poPreparedGeom->hGEOSCtxt);
+ if( hGEOSOtherGeom == NULL )
+ return FALSE;
+
+ int bRet = GEOSPreparedContains_r(poPreparedGeom->hGEOSCtxt,
+ poPreparedGeom->poPreparedGEOSGeom,
+ hGEOSOtherGeom);
+ GEOSGeom_destroy_r( poPreparedGeom->hGEOSCtxt, hGEOSOtherGeom );
+
+ return bRet;
+#else
+ return FALSE;
+#endif
+}
+
/************************************************************************/
/* OGRGeometryFromEWKB() */
/************************************************************************/
Index: ogr/ogrsf_frmts/generic/ogrlayer.cpp
===================================================================
--- ogr/ogrsf_frmts/generic/ogrlayer.cpp (revision 33064)
+++ ogr/ogrsf_frmts/generic/ogrlayer.cpp (working copy)
@@ -2033,6 +2033,8 @@
double progress_ticker = 0;
int bSkipFailures = CPLTestBool(CSLFetchNameValueDef(papszOptions, "SKIP_FAILURES", "NO"));
int bPromoteToMulti = CPLTestBool(CSLFetchNameValueDef(papszOptions, "PROMOTE_TO_MULTI", "NO"));
+ int bUsePreparedGeometries = CPLTestBool(CSLFetchNameValueDef(papszOptions, "USE_PREPARED_GEOMETRIES", "NO"));
+ int bPretestContains = CPLTestBool(CSLFetchNameValueDef(papszOptions, "PRETEST_CONTAINS", "NO"));
// check for GEOS
if (!OGRGeometryFactory::haveGEOS()) {
@@ -2093,10 +2095,49 @@
continue;
}
+ //
+ OGRPreparedGeometry* x_prepared_geom = NULL;
+ if (bUsePreparedGeometries) x_prepared_geom = OGRCreatePreparedGeometry(x_geom);
+
pLayerMethod->ResetReading();
while (OGRFeature *y = pLayerMethod->GetNextFeature()) {
OGRGeometry *y_geom = y->GetGeometryRef();
if (!y_geom) {delete y; continue;}
+
+ if (bUsePreparedGeometries &&
+ bPretestContains &&
+ x_prepared_geom && OGRPreparedGeometryContains(x_prepared_geom, y_geom))
+ {
+ OGRFeature *z = new OGRFeature(poDefnResult);
+ z->SetFieldsFrom(x, mapInput);
+ z->SetFieldsFrom(y, mapMethod);
+ if( bPromoteToMulti )
+ y_geom = promote_to_multi(y_geom);
+ z->SetGeometry(y_geom);
+ delete y;
+ ret = pLayerResult->CreateFeature(z);
+ delete z;
+ if (ret != OGRERR_NONE) {
+ if (!bSkipFailures) {
+ OGRDestroyPreparedGeometry(x_prepared_geom);
+ delete x;
+ goto done;
+ } else {
+ CPLErrorReset();
+ ret = OGRERR_NONE;
+ }
+ }
+ continue;
+ }
+
+ if (bUsePreparedGeometries &&
+ !((x_prepared_geom && OGRPreparedGeometryIntersects(x_prepared_geom, y_geom)) ||
+ x_geom->Intersects(y_geom)))
+ {
+ delete y;
+ continue;
+ }
+
OGRGeometry* poIntersection = x_geom->Intersection(y_geom);
if( poIntersection == NULL || poIntersection->IsEmpty() ||
(x_geom->getDimension() == 2 &&
@@ -2119,6 +2160,7 @@
delete z;
if (ret != OGRERR_NONE) {
if (!bSkipFailures) {
+ OGRDestroyPreparedGeometry(x_prepared_geom);
delete x;
goto done;
} else {
@@ -2129,6 +2171,7 @@
}
}
+ OGRDestroyPreparedGeometry(x_prepared_geom);
delete x;
}
if (pfnProgress && !pfnProgress(1.0, "", pProgressArg)) {
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
