Hi,
While improving the sum resampling method of gdalwarp, to make it really
preserve the sums of pixel values accross the whole raster, I need to do
a lot of intersections of polygons, and get the area of the
intersection. Typically this involves 4 intersections for each target
pixel in a reprojection scenario not modifying the resolution. You
intersect the reprojected shape of each contributing source pixel in
target pixel coordinates with the shape of the target pixel. And they
are just intersections of convex quadrilaterals.
My first try was with GEOSIntersects_r() and I had the feeling that the
result was quite slow compared to what I felt could be the optimum,
especially when looking at the stack trace which is quite deep, like ~
15 nested calls inside GEOS, memory allocations, etc.
I've coded a specialized getConvexPolyIntersection() function in
https://gist.github.com/rouault/e6a4cfa5952ad1c7601e254b0d25cc5a and it
is ~ 70 times faster than going through GEOS (tested against latest GEOS
master). I may have missed a few edge cases, and the O(N^2) complexity
of the algorithm probably makes it appropriate for very small polygons
(*), but perhaps there's some potential to improve GEOSIntersects for
such simple cases. Or offering a dedicated API.
Note: for my use case, the convex property is met in nearly all cases,
and in the unlikely event where one quadrilateral isn't convex, I can
just split it into 2 triangles, and fallback doing 2 intersections (one
of the quadrilateral - the target pixel - is always convex, actually a
square).
Even
(*)
https://tildesites.bowdoin.edu/~ltoma/teaching/cs3250-CompGeom/spring17/Lectures/cg-convexintersection.pdf
presents an algorithm for convex polygon intersection with O(N1+N2)
complexity.
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
geos-devel mailing list
geos-devel@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/geos-devel
