Geometry Boost List: [ggl] intersection of two vectors of polygons

Subject: [ggl] intersection of two vectors of polygons
From: Barend Gehrels (barend)
Date: 2011-04-29 17:52:57

Hi Vishnu,

OK, now I've studied your problem closely.

Indeed, the solution by first creating two multipolygons is wrong in 
this case, because you will lose your ID's of the individual polygons.

> The reason it's important to get 9 discrete polygons in the intersection is
> as follows. Say that 4 adjacent squares in the first collection are 4 farms
> labeled by ownership and the 4 adjacent squares in the 2nd set are labeled
> according to mineral rights. Let's say one wants the polygons corresponding
> to farm with owner "1" and mineral rights "3".

Sure, I understand it now.

> The problem, as I see it is that the multi_polygons are too restrictive in
> that they don't allow a collection of 4 adjacent squares such as a 2x2 cell
> on graph paper to be treated as discrete polygons.

It is not only this. A multi_polygon must be seen as a single geometry. 
It is one geometry, with one attribute (or one set of attributes). All 
the individual polygons making up that geometry share the same attribute(s).

If you're into (spatial) databases, see it as one row, with a 
geometry-column (the multi-polygon) and one or more other columns (the 
attributes).

Therefore, in your case, it is important to keep them in single 
polygons, each with its own ID.

> The alternative I am thinking of is to intersect each polygon in the first
> vector with each polygon in the second vector and then put the resulting
> polygons into a 3rd vector. This would be an O(n^2) operation, hence bad for
> large vectors of inputs.
Sure, this is a solution and for a small amount of input probably the 
best solution. I worked out your problem and will create either a sample 
or a blog of it. What makes it interesting is the additional attributes, 
we don't have many samples of that.

So here is the n^2 solution http://codepad.org/Zn35UPEx
and it gives this output:

n^2 solution
Owner 1 and mineral rights 1 intersect with area 0.0625
Owner 1 and mineral rights 3 intersect with area 0.0625
Owner 2 and mineral rights 1 intersect with area 0.0625
Owner 2 and mineral rights 2 intersect with area 0.0625
Owner 2 and mineral rights 3 intersect with area 0.0625
Owner 2 and mineral rights 4 intersect with area 0.0625
Owner 3 and mineral rights 3 intersect with area 0.0625
Owner 4 and mineral rights 3 intersect with area 0.0625
Owner 4 and mineral rights 4 intersect with area 0.0625

So it is probably comparable to your solution and output.

> Does anyone see another way? I think that having a way to do intersections
> on vectors of polygons and not just multi_polygons would be a valuable
> feature in boost-geometry.

All Boost.Geometry functions are based on geometry / geometry, and 
nothing on a collection. This is how most geometry libraries work. But 
that does not mean that there is a better solution.

Boost.Geometry has also the partition algorithm, which can visit a 
collection of geometries or (in this case) two collections of 
algorithms, relating them together spatially.

I worked out that as well and that solution is 
here:http://codepad.org/MBRPDhwj

As you see the loop is completely disappeared, the partition takes three 
policies, one to determine the boxes, one check if geometries are 
related (! disjoint), and one to do the actual thing (the visitor). It 
gives this output:

n-log-n solution
Owner 2 and mineral rights 1 intersect with area 0.0625
Owner 2 and mineral rights 3 intersect with area 0.0625
Owner 4 and mineral rights 3 intersect with area 0.0625
Owner 2 and mineral rights 2 intersect with area 0.0625
Owner 2 and mineral rights 4 intersect with area 0.0625
Owner 4 and mineral rights 4 intersect with area 0.0625
Owner 1 and mineral rights 1 intersect with area 0.0625
Owner 1 and mineral rights 3 intersect with area 0.0625
Owner 3 and mineral rights 3 intersect with area 0.0625

As you see the order is different but the contents is exactly the same.

I also made an image (originally SVG) of the use case:

showing the nine intersections together with their ideas. The complete 
program (with both approaches and SVG output) is here:

http://codepad.org/P2L7eyWD

(I realize there is code duplication of the registration of the polygon, 
etc, but with codepad it is convenient to have all just in one paste).

Note that if you are using Boost.Geometry 0.9, replace "return_envelope" 
by "make_envelope", same for centroid. If you are using Boost.Geometry 
0.8, I'm not sure but I think the partition didn't exist yet as a 
separate algorithm.

Thanks for this use case.

Finally, shortly on Luke's answer
> I support M different property types on M different multi polygons in a single sweep and avoid quadratic runtime

Similar story here, no property concept necessary...

> I have implemented it in Boost.Polygon (...) Unfortunately, the implementation details of boost.geometry polygon intersection make extending it to N multi polygons with N properties non-trivial.

... and I have to ask you kindly to please stop competing us on our own 
mailing list, and to verify what you are posting.

Regards, Barend

-------------- next part --------------
Skipped content of type multipart/related