Re: Equivalence between (multi-)polygons with spikes (4)
vs> The resulting multi-polygon polygonC is the empty multi-polygon. This is consistent with the above: The part of the polygon, that is left after the difference operation, is mathematically non-empty, but in its int representation it is a spike, and thus is removed entirely.
bg> OK, yes we now remove spikes in spatial intersection operations. Because they make the output invalid. I think that should be done.
vs> On the other hand, the result seems to be at odds with the fact that boost::geometry::covered_by(polygonA, rectB) returns false. (Note: As of 1.57.0, covered_by is not implemented for polygon and box, so for this test I converted the box to a polygon. I assume that this doesn't affect the argument of this discussion.) As far as I understand, that's "by design" in boost::geometry.
bg> We should come back to that later. Indeed - if the output would have been a spike which is then removed, the results might be mutually inconsistent... I see what you mean.
Exactly. I'm not saying that this is a huge problem, but it certainly is a problem that warrants some thought and needs to be made explicit. You will always have to make compromises simply because we do not have the means to accurately represent sets of points in a continuous space. Calculation on doubles always involves rounding error and inaccuracies. We've come to live with this, but we have to make explicit the compromises that we make.
Which brings up a very interesting idea that could solve these kinds of issues for int-based geometrics: Apparently, right now, in boost::geometry int-based geometrics are no different from floating-point geometrics, except that coordinates always happen to be at integral numbers. This is meant to say that you are actually representing non-rectilinear shapes with int coordinates. Let's consider a simple triangle: (1 1)(1 2)(3 1). It is currently my understanding, that in boost::geometry this data represents a set of points shaped like a triangle. The point (1 1.6) is included in this set, and the point (1 1.7) is not. The point (2 1.3) is included, and the point (2 1.4) is not.
What if -- for the purpose of int-based geometrics -- you were regarding the world as composed from discreet square bricks of 1x1? The above outline could still be a valid representation, but it would be equivalent to (1 1)(1 2)(2 2)(2 1)(3 1), which contains a spike and is equivalent to (1 1)(1 2)(2 2)(2 1). In this very simple example, you end up with a square. Larger structures would map to shapes with a "saw tooth" profile. The nice part about this: There are actual, discreet, denumerable n-dimensional points, and for every point you can tell with 100% certainty whether or not it is part of a given geometry. That's fundamentally different from what you try to achieve in floating-point geometrics, and a lot of the problems we are currently discussing would just "go away".
My apologies for the rant. Maybe it's worth exploring, maybe not.