[geos-devel] (Much) faster intersection of convex quadrilaterals

[Even Rouault]

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.