Halfplane Intersection

[Chillee]

Currently, this code is somewhat a mix of Zimpha's code, MIT's code, and SJTU's code.

One thing in particular I thought was interesting was this comparator from MIT. Considering they have the more "obvious" comparator commented out, I presume this one has superior numerical precision.

This code has been tested on an easy half plane problem.

[ecnerwala]

Quick question: do you support infinite intersections? That may or may not be good to include.

[Chillee]

I don't think we do. I think supporting infinite intersections can be done fairly easily manually, no? Just add a really big bounding box.

Ie:

lines.push_back({(INF, INF), (-INF, INF)});
lines.push_back({(-INF, INF), (-INF, -INF)});
...

[ecnerwala]

Sure, good solution. I remember Gennady had some code that supported it, but bounding box works fine.

[Chillee]

I believe this is now ready for review.

Previously I was using MIT's code to benchmark, but MIT's code seems to fail on this case: https://ideone.com/gPbt7X, which corresponds to these half planes (with a bounding box of length 20 centered at 0). I suspect the issue is that MIT's code doesn't handle lines with the same slope correctly.

MIT's code outputs -14.2857.

I've switched to using TankEngineer's half plane code for fuzz testing.

Re: Precision issues, there are 3 types of precision issues which may occur.

1. lineInter calls: These are all in the form lineInter(input line, input line). I don't see any way to get around these calls. At the very least, you need them for the final answer.
2. sideOf calls: These are all in the form sideOf(input line, computed point from lineInter). Once again, I don't see any way to get around these calls. These are how you determine whether the newly added half plane eliminates a previous point or not. We already have a version of sideOf that supports some epsilon handling, so I don't think there's much to do here.
3. angDiff sorting/equality: This is where we might be able to get some gains. As @simonlindholm pointed out, using atan2 is only accurate for points up to 2^26, as angles of the form x/y and z/w can differ by as little as 1/yw. We may want to provide instructions for how to do a more accurate sort.

[Chillee]

Thanks to @xyzhan, I found a bug in the current half-plane code. Need to think about how to address it.

I would also have found the bug if I ran the fuzz tests for a higher number of iterations.

[Chillee]

Precision starts to become a large issue even with very small inputs. See this example

    vector<Line> t({{P(8,9),P(8,2)},{P(3,9),P(5,2)},{P(8,2),P(8,3)},{P(7,2),P(1,7)},{P(1,0),P(7,1)},{P(9,2),P(5,6)},{P(10,10),P(-10,10)},{P(-10,10),P(-10,-10)},{P(-10,-10),P(10,-10)},{P(10,-10),P(10,10)}});
    auto res = halfPlaneIntersection(t);
    cout<<polygonArea2(res)<<endl;
    cout << sjtu::halfPlaneIntersection(t)<<endl;

Our code currently outputs -561.8, when the correct answer is 0. This can be fixed in a number of ways:

1. Add epsilon checks to sideOf.
2. Use the lineDistance function from SJTU (currently in lineIntersection under the name lineInter2. I don't know if this is because of better precision or it just happened to work out.

I think we should probably add epsilon checks to sideOf. It appears that half plane intersection can fail disastrously even with small inputs.

With either one of those fixes, the function passes as much fuzz testing as I have been able to run.

[Chillee]

Somewhat interestingly, it seems like lineInter2 barely makes a difference in accuracy, but it seems to be enough to matter (ie: lineinter2 passes all the fuzz-tests, lineinter doesn't). I'll note that SJTU's and the (updated) MIT Halfplane code seem to always output the same results.

[Chillee]

Ok, I've updated with a stable version. It currently passes some fairly exhaustive fuzz testing.

[simonlindholm]

Is the only fix to the precision issue the new lineIntersection? Naively, that would seem to fix only the case with integer coordinates <= ~2^16, since lineIntersection uses products of three numbers, and above that the divisions are no longer precise. (If it needs mathematically unequal divisions to generate unequal results it might be lower, and some property of the code might well make it higher). That's pretty low...

Do we have a nice characterization of the failure case? Is it only for area 0, for one?

[Chillee]

Ok, so I modified the halfplane.h to use fuzzy results for line intersection, and I played around with various values of EPS for sideOf. The code looks like this:

const double EPS = 1e-4;
const double mult = 2;
std::uniform_real_distribution<> dist(-EPS, EPS);
while (ah<at && sideOf(sp(vs[i]),ans[at-1], mult*EPS) < 0) at--;
while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah], mult*EPS)<0) ah++;
auto res = lineInter(sp(vs[i]), sp(deq[at]));
res.second = res.second + P(dist(e2), dist(e2));

With mult=0, literally any EPS > 0 caused halfplane to fail catastrophically. I did a wide range of testing, and I strongly suspect that for any EPS, we need a mult >= sqrt(2)/2 ~= 1.414 in order to prevent catastrophic failure. With an mult=1.4, I was able to induce catastrophic failures for EPS=1e-4, 1e-8, 1e-20. With a mult=1.414, I was not able to induce any catastrophic failures.

[simonlindholm]

Apparently the implementation is a bit slow; it got within 66% of the TL on https://ncpc20.kattis.com/problems/bigbrother. Probably due to calling atan2 in the comparator.