<div dir="ltr">An interesting idea for sure.<div><br></div><div>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?  </div><div><br></div><div>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)  </div><div><br></div><div>[1] <a href="https://wrfranklin.org/nikola/pubdetails/f-cmopa-90.html">https://wrfranklin.org/nikola/pubdetails/f-cmopa-90.html</a></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Oct 18, 2022 at 9:55 AM Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi,<br>
<br>
While improving the sum resampling method of gdalwarp, to make it really <br>
preserve the sums of pixel values accross the whole raster, I need to do <br>
a lot of intersections of polygons, and get the area of the <br>
intersection. Typically this involves 4 intersections for each target <br>
pixel in a reprojection scenario not modifying the resolution. You <br>
intersect the reprojected shape of each contributing source pixel in <br>
target pixel coordinates with the shape of the target pixel. And they <br>
are just intersections of convex quadrilaterals.<br>
<br>
My first try was with GEOSIntersects_r() and I had the feeling that the <br>
result was quite slow compared to what I felt could be the optimum, <br>
especially when looking at the stack trace which is quite deep, like ~ <br>
15 nested calls inside GEOS, memory allocations, etc.<br>
<br>
I've coded a specialized getConvexPolyIntersection() function in <br>
<a href="https://gist.github.com/rouault/e6a4cfa5952ad1c7601e254b0d25cc5a" rel="noreferrer" target="_blank">https://gist.github.com/rouault/e6a4cfa5952ad1c7601e254b0d25cc5a</a> and it <br>
is ~ 70 times faster than going through GEOS (tested against latest GEOS <br>
master). I may have missed a few edge cases, and the O(N^2) complexity <br>
of the algorithm probably makes it appropriate for very small polygons <br>
(*), but perhaps there's some potential to improve GEOSIntersects for <br>
such simple cases. Or offering a dedicated API.<br>
<br>
Note: for my use case, the convex property is met in nearly all cases, <br>
and in the unlikely event where one quadrilateral isn't convex, I can <br>
just split it into 2 triangles, and fallback doing 2 intersections (one <br>
of the quadrilateral - the target pixel - is always convex, actually a <br>
square).<br>
<br>
Even<br>
<br>
(*) <br>
<a href="https://tildesites.bowdoin.edu/~ltoma/teaching/cs3250-CompGeom/spring17/Lectures/cg-convexintersection.pdf" rel="noreferrer" target="_blank">https://tildesites.bowdoin.edu/~ltoma/teaching/cs3250-CompGeom/spring17/Lectures/cg-convexintersection.pdf</a> <br>
presents an algorithm for convex polygon intersection with O(N1+N2) <br>
complexity.<br>
<br>
-- <br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
My software is free, but my time generally not.<br>
<br>
_______________________________________________<br>
geos-devel mailing list<br>
<a href="mailto:geos-devel@lists.osgeo.org" target="_blank">geos-devel@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/geos-devel" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/geos-devel</a><br>
</blockquote></div>