An interesting idea for sure. Do you need only the area of the intersection? If so, perhaps the "area-only overlay" idea of William Franklin [1] would provide an even faster/simpler implementation?
I have a prototype implementation of this in JTS. One advantage might be that it will work for any polygons, although will only be really efficient for smaller/simpler polygons. It would be worth investigating how fast an implementation for low-vertex-count polygons would be (i.e that would using simple loops rather than indexed structures) [1] https://wrfranklin.org/nikola/pubdetails/f-cmopa-90.html On Tue, Oct 18, 2022 at 9:55 AM Even Rouault <even.roua...@spatialys.com> wrote: > 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 >
_______________________________________________ geos-devel mailing list geos-devel@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/geos-devel
