Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Halfplane Intersection #90

Open
wants to merge 48 commits into
base: main
Choose a base branch
from

Conversation

Chillee
Copy link
Collaborator

@Chillee Chillee commented May 3, 2019

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
Copy link

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

@Chillee
Copy link
Collaborator Author

Chillee commented May 3, 2019

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
Copy link

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

@Chillee Chillee force-pushed the halfplane branch 2 times, most recently from b246350 to 95b770b Compare May 4, 2019 21:02
@Chillee
Copy link
Collaborator Author

Chillee commented May 6, 2019

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
Copy link
Collaborator Author

Chillee commented May 6, 2019

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
Copy link
Collaborator Author

Chillee commented May 6, 2019

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

@simonlindholm
Copy link
Member

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
Copy link
Collaborator Author

Chillee commented Oct 31, 2019

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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants