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:
(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