Intersecting a patch with a path3 can be extremely slow

[charlesstaats]

Here's one simple example:

Welcome to Asymptote version 2.36-19 (to view the manual, type help)
> import three;
> patch plane = patch(O -- (0,1,1) -- (1,1,2) -- (1,0,1) -- cycle);
> path3 line = (0,0,-1e-5) -- (1,1,2);
> cputime();
0.08u 0.01s 0.14U 0.03S 
> cputime();
0.00u 0.00s 0.14U 0.03S 
> real[][] isections = intersections(line, plane);
> cputime()
46.06u 0.04s 46.20U 0.07S 
> isections.length
1

Intersecting a line with a plane should not take 46 seconds, even if the line is nearly parallel to the plane. Here's a somewhat more complex example:

import three;
triple[][] P = {{(-0.866025403784439,0.5,0), (-0.866025403784439,0.5,-0.552284749830793), (-1.25375818409268,0.723857625084604,-1), (-1.73205080756888,1,-1)},
        {(-0.953793735509369,0.347980790156861,0), (-0.968163779838173,0.352863275764704,-0.552284749830793), (-1.39026076109554,0.506984178095615,-1), (-1.90758747101874,0.695961580313722,-1)},
        {(-1,0.175536663449861,0), (-1.01488606624617,0.17849332904426,-0.552284749830793), (-1.45749322604123,0.256069203000192,-1), (-2,0.351073326899723,-1)},
        {(-1,2.77555756156289e-16,0), (-1,2.77555756156289e-16,-0.552284749830793), (-1.44771525016921,4.01821700959705e-16,-1), (-2,5.55111512312578e-16,-1)}};
path3 line = (-1.61075487116513,4.47074286248665e-16,-0.921415042944955)--(500,400,200);
write(cputime());
intersections(line, P);
write(cputime());

The resulting output:

0.08u 0.02s 0.15U 0.03S 
32.66u 0.03s 32.81U 0.06S

In other words, this intersection took over 32 seconds to compute.

Possible strategy for fixing: in the final intersections function in path3.cc, before computing bounding boxes, apply a rotation (to both the path and the patch) so that three of the patches four corners lie in a plane parallel to the xy-plane. When the patch subdivision gets small enough, the patches will be nearly planar, and applying the suggested rotation can drastically reduce the volume of the bounding box (and hopefully also the extraneous path segments it catches).

Another suggestion, specifically targeted at the second example, is to keep halving the path3 until both halves intersect the bounding box of the patch. (Then put those two halves back together?) The idea is to make sure that extremely long path segments get cut down to size in the very first step; I don't know how well it would really work.

[johncbowman]

These are good ideas to think about for the future. Fortunately, commits f82e767 and
6e2f7c5
have reduced these execution times substantially: the second example now runs about 15 times faster than it did.

[johncbowman]

Unfortunately, in ddf49a3 I had to revert
some of the changes in f82e767 to avoid duplicate interesections. So the second example is slow again.

[charlesstaats]

A thought I had on this problem this morning: Ideally, you would test whether the convex hulls of the control points intersect; using the bounding boxes instead of the control points is a quick-and-dirty way to get a valid result, but it can also give rise to lots of false positives that make for many unnecessary iterations. The "rotation" suggestion is meant to reduce these false positives, but it's not particularly elegant.

But consider this: To tell whether two convex hulls intersect, you don't actually have to compute the convex hulls -- you only have to determine whether the two sets of points are linearly separable. And that can be done using linear programming. This should provide a fairly straightforward solution that is both more powerful and more appealing than the "rotation trick" suggestion.