11
$\begingroup$

I need to compute the intersection(s) between two circular arcs. Each arc is specified by its endpoints and their height. The height is the perpendicular distance from the chord connecting the endpoints to the middle of the arc. I use this representation because it is numerically robust for very slightly bent arcs, as well as straight line segments, for which the height is zero. In these cases, representing an arc using the center of its circle could lead to the center being far far away from the endpoints of the arc, and hence, numerically unstable.

My question at the highest level is how I would go about computing the intersection points, given that the centers of the circles of the arcs cannot necessarily be computed robustly. At a lower level, I am wondering if there is a parameterization of an arc using only the information I have stated above (which does not include the circle center). Of course, keep in mind numerical robustness is my principal concern here; otherwise I would just do the naive thing and compute the circle center for all non-linear arcs and hope for the best.

Edit: Formula for computing center of circle of arcs:

Suppose the chord length is $2t$, and the height is $h$. The distance from the chord to the circle center is $c$, so that $r=h+c$. Then it follows that $c=(t^2-h^2)/2h$, which breaks down when $h$ is very small. Computing the location of the circle center is some simple vector arithmetic using the chord vector and its perpendicular.

  • 0
    Cool, so how are you drawing your arcs (i.e., what parametric equations are you using to draw arcs)?2010-12-30
  • 2
    Haven't thought about it much, but if you have line-arc intersection available, I believe you can (robustly) reduce arc-arc to line-arc using inversion. btw, two endpoints + height determines two arcs.2010-12-30
  • 1
    @Moron: the height is a signed number, and taken in a definite orientation depending on the orientation of the endpoints. In short, the height uniquely defines a single arc.2010-12-30
  • 1
    Yeah, just nitpicking (distance usually is unsigned). I was pretty sure you had something to uniquely identify. Just pointed that to reduce possible ambiguity people might have come across...2010-12-30

2 Answers 2

4

This is an interesting conundrum, and what I will suggest is just one possible approach, by no means a definitive answer. Perhaps what you could do is compute Bézier curve segment approximating your arcs, and then compute the intersection between Bézier segments. To compute the Bézier segment for an arc, you need the tangent vector at one endpoint (the tangent at the other endpoint is obtained by symmetry). So this approach reduces to finding a tangent angle at one endpoint.

Let the angle subtended by an arc at the circle center be $\theta$. If the arc's chord length is $c$ and its height $h$, then, if I've calculated correctly, $$\frac{\theta}{2} = \cos^{-1} \left( \frac{c^2 - 4h^2}{c^2+4h^2} \right) \;.$$ This avoids calculating the circle radius $r = (c^2 + 4h^2)/(8h)$. Now from $\theta$ you could compute the needed tangent angle, and from there a Bézier segment.

I have not analyzed this to see if it is indeed robust.

  • 0
    If Victor doesn't need an exact solution anyway, I suppose crossing Béziers ought to do (though if you look at it, it's a bit funny that you're solving a problem of intersecting quadratics by intersecting cubics... :) ). +1 anyway!2010-12-30
  • 0
    @J.M.: Yes, you are right, this approach is a bit counterintuitive, and I am uncertain if it worthwhile. And you are right that if he wants an exact solution, it cannot suffice. However, nearly every circle we see (e.g., in Adobe Illustrator/Photoshop/Flash/etc.) is a Bézier approximation to a true circle.2010-12-30
  • 0
    This is an interesting solution, and I particularly like the formula for the tangent angle. I think that a rational quadratic B-spline can represent circular arcs exactly, but I haven't looked into how hard it is to compute intersections with them.2010-12-30
  • 0
    @Victor: They do, but I'm getting the feeling that the numerical instability you were trying so hard to avoid for arcs of tiny curvature would find a way to manifest itself when you use rational B-splines... anyway, you have to experiment.2010-12-30
  • 0
    @Joseph: I'm not sure about the last statement in your comment. In most of the programs you mention, ellipses and Bézier splines are distinct primitives, and the program writers would want to exploit [the fast known methods for drawing circles and ellipses](http://en.wikipedia.org/wiki/Midpoint_circle_algorithm). I expect a Bézier approximation only takes place when the user explicitly asks to convert an ellipse to an editable Bézier spline.2010-12-30
  • 1
    @Rahul: Create a circle in Adobe Illustrator, and note that it is composed of four arcs with tangents at each endpoint. These tangents can be stretched. Circles in Illustrator are four Bézier segments.2010-12-30
  • 0
    @Joseph: I stand corrected. I don't have access to Illustrator, so I had assumed it functioned similarly to Inkscape, in which ellipses and Bézier splines are distinct. For what it's worth, the same goes for Powerpoint and a number of 2D graphics APIs.2010-12-30
5

Have you considered finding the intersections using an implicit form for the circles, $$\frac{x^2}{r^2} + \frac{y^2}{r^2} + ax + by + c = 0?$$ This representation doesn't have any coefficients that diverge as the circle approaches a straight line. To find intersections, you'll have to solve a quadratic equation whose leading coefficient could be zero or arbitrarily close to it, but the alternative form of the quadratic formula should be able to deal with that robustly.

You'll then have to do some jiggery-pokery to figure out whether the intersection points lie within the arcs. If the arc's bending angle is smaller than $\pi$, a projection onto the line joining the endpoints will suffice.

(Disclaimer: While all of this feels like it should work, I haven't analyzed it in any detail. Also, there could still be a problem when the circle is close to a line and you want the longer arc. But I can't imagine that's a case that would turn up in any practical application.)

Update: For a concrete example, here is the equation for a circular arc passing through the three points $(0,0)$, $(0.5, h)$, and $(1,0)$: $$\kappa^2 x^2 + \kappa^2 y^2 - \kappa^2 x - 2\eta y = 0,$$ where $$\begin{align}\kappa &= \frac{8h}{4h^2 + 1}, \\ \eta &= \frac{8h(4h^2-1)}{(4h^2+1)^2}.\end{align}$$ As you can see, the coefficients remain bounded as $h \to 0$.

Update 2: Wait, that equation becomes trivial if $h = 0$, which is bad. We really want something like $x^2/r + y^2/r + ax + by + c,$ i.e. multiply the previous expression through by $r$. Then for the same example, our equation becomes $$\kappa x^2 + \kappa y^2 - \kappa x - 2\eta' y = 0,$$ where $\eta' = (4h^2-1)/(4h^2+1)$. Here are some explicit values.

$h = 1/2$: $$2 x^2 + 2 y^2 - 2 x = 0,$$ $h = 0.01$: $$0.07997 x^2 + 0.07997 y^2 - 0.07997 x + 1.998 y = 0,$$ $h = 0$: $$2 y = 0.$$

By the way, in this format, the linear terms will always be simply $-2(x_0/r)x$ and $-2(y_0/r)y$, where the center of the circle is at $(x_0,y_0)$. As the center goes to infinity but the endpoints remain fixed, these coefficients remain bounded and nonzero (i.e. not both zero).

  • 0
    I gathered from his problem that he's only considering small arcs and not big arcs... anyway, it theoretically works, but I'm not that confident that the ill condition he was trying so hard to avoid would not manifest in the coefficients of the quadratic.2010-12-30
  • 0
    @J.M.: I'm pretty sure it wouldn't. I'll see if I can work out the coefficients explicitly and add them to my answer when I have time.2010-12-30
  • 0
    @Rahul: I'll have to work out the coefficients to see if this is a stable route. Taking O'Rourke's advice, I am looking into Bezier intersection (http://tom.cs.byu.edu/~tom/papers/C3CIC.pdf) and using implicitization, I feel like it might reduce to your solution. Also, you can disregard the case where the arc is almost a complete circle since I would simply perform a circle-arc or circle-circle test in that case. See also my comment above about the uniqueness of the arc specification.2010-12-30
  • 0
    @Victor: I added an equation for the coefficients in a concrete example.2010-12-30
  • 0
    @J.M.: Please see my edit. Unfortunately [only one person can be notified per comment](http://meta.stackexchange.com/questions/43019/how-do-comment-replies-work/43020#43020)...2010-12-30
  • 0
    I'm trying out values of $\kappa$ near the square root of machine epsilon, and it does look like a good approach... :)2010-12-30
  • 0
    @Rahul: Your coefficient for the $x$ term doesn't look right, by dimensional analysis. I worked on calculating the coefficients briefly last night and I obtained expressions for the linear terms which depend on the radius or position of the origin of the circle (the coefficients of the linear terms must have units of 1/distance).2010-12-30
  • 0
    @Victor: That's because I fixed the $x$-coordinate of the center at 0.5 units, so that doesn't "show up" in the final expression. For what it's worth, before I simplified it, the linear term in $x$ was $2\kappa^2 x_0 x$, which does have the right units. Also, I was a little careless before; please see my edit.2010-12-30
  • 0
    @Victor: If you wanted to use this approach, I'd suggest using the formula above and simply transforming it to fit each arc by appropriate translation, rotation, and scaling. Quadratic forms are well-behaved under such similarity transformations. Also, if your arcs can be far from the origin, it's probably best to work in the coordinate frame of one of the arcs to minimize numerical error in the constant term.2010-12-30
  • 0
    @Rahul: Under your new scaling, the constant term $c$ must have units of distance, and scales with the radius of the circle (albeit there's a cancellation which might make it well behaved). In your examples, you always choose to place one endpoint at the origin, but for intersection tests, this is not possible. My expression for $c$ is $(x_0^2+y_0^2)/r-r$. I'm not sure how I would compute that robustly.2010-12-30
  • 0
    @Victor: Our comments must have gone past each other. I'd suggest forming the quadratic in a coordinate frame attached to the arc in some standard way (so the points are $(0,0), (\ell/2, h), (\ell,0)$ or maybe $(-\ell/2,0), (0, h), (\ell/2,0)$), in which the coefficients can be calculated robustly, and then transforming it by translation and rotation to the actual placement of the arcs in world space.2010-12-30