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

[Martin Davis]

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.rouault at 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 at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/geos-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20221018/59544759/attachment.htm>