Curve curve intersection, take two #199
Conversation
|
One thing to note that might be useful here is that quad beziers can easily be converted into cubic beziers. I haven't started looking through the code, so I'm not sure if a function for this exists already, but if not, take a look at this: |
There was a problem hiding this comment.
Overall, it looks good. I'm not a contributor yet, so I'm not sure that I should change the status of this PR (approval or change-request), but I would recommend considering some of my ideas before merging this PR since I've gone through this before and the edge cases can be killer.
Also, let me know if you need/want me to implement the changes that I've requested, or if you have questions about them.
| let t2 = (domain2.start + domain2.end) * 0.5; | ||
| // There's an unfortunate tendency for curve2 endpoints that end near (but not all | ||
| // that near) to the interior of curve1 to register as intersections, so try to avoid | ||
| // that. (We could be discarding a legitimate intersection here.) |
There was a problem hiding this comment.
I've also come across this issue when I implemented the algorithm in c++. However, when it comes to doing boolean operations of paths (which is what I'll hopefully be contributing to this project), having false intersections is probably better than missing real intersections. Also, I think that this check means we'll miss ALL real end-point intersections, which I believe are valid intersection points.
Perhaps the checks should be for t2 < -end_eps and t2 > 1.0 + end_eps instead? However, one issue I've come across before is that using t parameter comparisons can be an issue if one curve is much longer than the other, because in that scenario, 0..end_eps is much larger in x/y space for the long curve than it is for the short one. It's something to think about.
| Some((min, max)) => (min, max), | ||
| None => return call_count, | ||
| }; | ||
|
|
There was a problem hiding this comment.
One thing that I've found when implementing this algorithm in the past is that you can sometimes get much faster convergence if you also check a perpendicular fat-line. This is the case in scenarios where the curves are someone parallel and when one curve is like a sharp 'u' which is perpendicular to a longer flatter curve (there are probably other cases too).
This article talks a bit about the same things, so I would recommend giving it a quick scan:
https://losingfight.com/blog/2011/07/07/how-to-implement-boolean-operations-on-bezier-paths-part-1/
In case you're interested or curious, here are the other 2 parts to the article:
https://losingfight.com/blog/2011/07/08/how-to-implement-boolean-operations-on-bezier-paths-part-2/
https://losingfight.com/blog/2011/07/09/how-to-implement-boolean-operations-on-bezier-paths-part-3/
| intersections: &mut ArrayVec<[(f64, f64); 9]>, | ||
| flip: bool, | ||
| ) { | ||
| let pt = pt_curve.start(); |
There was a problem hiding this comment.
Since pt_t is the average between the start and end t-value, should this be the average between the start and end points? After all, pt_curve is so small that it should be essentially linear. Alternatively, this could be pt_curve.eval(pt_t) if pt_t is defined sooner.
| pt: Point, | ||
| curve: &T, | ||
| epsilon: f64, | ||
| ) -> ArrayVec<[f64; 9]> { |
There was a problem hiding this comment.
This function is basically finding the t-value(s) that correspond to a given point that is supposedly along a curve. I think this function could be really useful outside of the Bezier clipping, so it might be nice to have it take care of it's duplicates and also to calculate epsilon on it's own since epsilon is simply epsilon_for_point(pt).
There was a problem hiding this comment.
If you don't mind, I think I'll have a crack at implementing (and testing) this. I'll try to make the function as general as possible so that it can be useful to other parts of the code as well. When I'm done I'll pull request it into your branch so that you can review it.
|
|
||
| // Generally speaking |curve| will be quite small at this point, so see if we can get away with | ||
| // just sampling here. | ||
|
|
There was a problem hiding this comment.
point_curve_intersections seems like it would catch scenarios where a pt lies directly on curve at 2 points, but this section doesn't seem like it would catch that scenario. I'm not sure that this case is common enough to worry about, but maybe add a comment here about that in case it ends up being an issue later on.
| t2: f64, | ||
| orig_curve2: &U, | ||
| flip: bool, | ||
| intersections: &mut ArrayVec<[(f64, f64); 9]>, |
There was a problem hiding this comment.
What if the number of intersections is temporarily larger than 9 because of getting repeated intersections of two very close curves?
| for i in 0..intersections.len() { | ||
| let (old_t1, old_t2) = intersections[i]; | ||
| // f64 errors can be particularly bad (over a hundred) if we wind up keeping the "wrong" | ||
| // duplicate intersection, so always keep the one that minimizes sample distance. |
There was a problem hiding this comment.
When it comes to the path boolean operations that I've mentioned before, it might be more useful to keep the start and end intersections when there is a section of many 'duplicate' intersections caused by the curves being very close together. Alternatively, maybe there is a way that the boolean operations could use solve_t_for_x and solve_t_for_y on the start/end points to find the range in which 2 curves are 'very close' before finding the intersections.
|
Thanks for this! Please note that the whole implementation, comments-and-all, is lifted from Lyon: https://github.com/nical/lyon/blob/master/crates/geom/src/cubic_bezier_intersections.rs So if there are improvements, maybe they should head upstream as well? |
|
That's a good question. It looks like lyon is still being actively developed, so I suppose they should head upstream as well. We should test the improvements thoroughly first, but those devs may also have some insight into the changes we're making, which means we'll get even more eyes on it. |
|
Hey @simoncozens. I just wanted to check in and see if you saw the pull requests hat I've made into your repo. I thought it might be best to do it that way so that this pull request remains valid but is updated with the changes that I've made. |
|
I had not! I shall take a look. |
|
You know what, I think I should just be totally honest here: I really want this functionality to exist; but I also really don't want to spend much more of my limited programming time on nasty algorithmic geometry edge-cases. I thought I could get away with grabbing the code from Lyon and making it work, but if people are going to want me to actually understand it and improve it, I don't think I'm ever going to find the time or headspace to do so. So if anyone wants to take over this PR, please do so. |
|
I'm in a similar boat, I want this functionality in kurbo, but also don't have the bandwidth to wrap my head around it. So here's what I think we should do: Merge this code basically as-is, with a disclaimer added to the doc comment that it hasn't been validated thoroughly for edge cases and may not be numerically robust. Then I'll file an issue for the robustification. I have some thoughts on this topic; it's probably a valid research topic in 2D computational geometry, and to my knowledge a numerically robust cubic-cubic intersection problem hasn't been published. Ideally that would attract somebody interested in the topic, or, at some point, I'll probably have a whack at it myself. Sound good? |
|
I'm working on this project because I need this algorithm for my job. I've done this before in c++, so I have quite a bit of experience with the edge cases (and I'm not quite sure that you can ever make this algorithm 100% robust). As I've mentioned to Simon, my main interest is in creating boolean path operations (e.g. union of two closed paths, etc) and Kurbo is a good base for that project. If this doesn't interest others, then I'll just create my own fork/project and work from there. However, I think the others in the community would probably benefit from this project. @simoncozens and @raphlinus If what I've said above sounds good to you, I'm happy to take over this pull request. I've posted 2 pull-requests into this branch just so that you can see what I've been doing, but if I take over this pull request, I'll just merge them in. Let me know what you think! |
|
@tylers-epilog I'm also happy for you to take it over! I do think it's possible to make it completely numerically robust, but I also think that's a hard problem. Basically the main thing I'd ask is that there's enough testing in place to be reasonably confident there aren't major regressions from the Lyon case. It might not be a bad idea to merge this then have Tyler's work as a follow-on PR. |
|
I guess it depends on what you mean by numerically robust. There's always going to be some amount of floating-point error. However, I do think that the improvements I've made make the accuracy pretty good. I think that preventing regression from Lyon might be an issue. Mainly due to the fact that they have a bunch of random epsilon values and they don't capture intersections that happen at the endpoints. In the work that I've done, I've added an input flag that allows developers to choose whether they want to capture endpoint intersections or not, but the definition of an "endpoint" is of course a little fuzzy. I'm trying to make my code as useful to as many people as possible, but I don't know what constraints/edge cases were used in making some of their design choices. As such, I'm not sure that I can make it work for my project AND make it 100% compatible with Lyon unless I have them test my changes. Maybe my code would work for them, maybe it wouldn't. What do you think is the right course of action in this case? I'm totally fine with merging this in and then posting a new pull request with my changes. |
|
The first step in numerical robustness is coming up with a good definition of what that means. That should probably be done in the context of the task, in this case boolean path operations. A reasonable first cut is that the curves must be close (within some epsilon) at the intersection points, and the relative orientations must be consistent (this would need to be precisely defined, I'm just sharing my intuition). It's possible that for "relative orientation" to be well defined, the API should be changed so that it reports at least whether the order of intersection is odd or even, ie whether there is a sign change or not (y = x^2 intersect y = 0 is an example of a second-order intersection, and I believe we can get up to fourth-order by intersecting two cubics). Another approach is to detect the numerically unstable cases and upgrade to higher precision arithmetic; this is what Shewchuk does in his Triangle project. I'd prefer not to do that and think it's not necessary, but that approach shows it is possible. But let's take things in order. Let's land this PR with minimal changes and a disclaimer (it sounds like Tyler's work also makes minor changes to the API), then I'll set aside some time to review Tyler's PR. I also have some ideas on how to generate adversarial examples. |
|
That sounds good. I guess that just means that @simoncozens needs to make a commit adding a disclaimer that this code hasn't been tested against every edge case? Once my PR is up, we can discuss the numerical robustness of my code and how to go about any improvements that you think can be made since you seem to have a pretty good grasp on the subject. We can also then discuss some of the adversarial examples that you have in mind and you can compare them to the test cases that I've added. |
|
I found https://cs.nyu.edu/exact/doc/subdiv1.pdf as an example of the "bigfloat" approach. It might make for good reading :) |
|
That all works for me. I'll add a doc disclaimer in a bit. |
|
I wrote a lot more thoughts on this on the Zulip: https://xi.zulipchat.com/#narrow/stream/260979-kurbo/topic/B.C3.A9zier-B.C3.A9zier.20intersection |
|
Aside from some of it being way over my head :), I like what you've written down and I think you make some good points. Getting a numerically robust algorithm now will save us a lot of time later and will prevent having to code for edge cases. In particular, I'd like to see your work on calculating the Hausdorff distance. In my experience of the path operations, two curves lying very close to each other is one of the biggest issues. In my code that I'll post, I've added a check to see if two curves "lie on top" of each other, but your Hausdorff method is probably a better check. |
|
Hey @simoncozens, I just wanted to check in and see if we could get that commit with the disclaimer. Cheers! |
This borrow's lyon's curve-curve intersection algorithm. Some notes:
curve_intersections.QuadBez, and we would get cubic/quad and quad/quad intersections for free. In practice, I don't understand the algorithm well enough to implement the convex hull thingy on quads.ParamCurveBezierClippingtrait is a bit of a grab-bag of whatever methods were needed to get the algorithm working. Some of these methods (particularlysolve_t_for_xandsolve_t_for_x) may want to go into other traits - but which traits?Again, for discussion.