Circles In Angles

I did find an approach that feels more elegant. To start with, let’s assume that we have normalized unit-length vectors for the direction of each ray:

\[ \hat{U_1} = ( P_1 - V ) / \lVert P_1 - V \rVert \] \[ \hat{U_2} = ( P_2 - V ) / \lVert P_2 - V \rVert \]

For at least one of my applications, I already had these on hand from other calculations so using these didn’t cost anything extra. Because they have equal length, their vector sum, \(\hat{U_1}+\hat{U_2}\), bisects \(\angle P_1VP_2\) and by symmetry the center of the circle must be somewhere along that line. Now consider the line segment from \(V+\hat{U_1}\) to \(V+\hat{U_1}+\hat{U_2}\). Because we got it by adding a unit-length vector, its length must be one. And by parallel lines, the angle it forms with line \(\overline{VP_1}\) is congruent to \(\angle P_1VP_2\). Taking this line segment as a hypoteneuse, the distance of \(V+\hat{U_1}+\hat{U_2}\) from the line \(\overline{VP_1}\) (and by symmetry from \(\overline{VP_2}\)) must therefore be the sine, \(s\), of \(\angle P_1VP_2\):

It’s well known that the cosine of the angle between two vectors is related to their dot product. What’s less well known (though fairly easy to see) is that the sine of the angle between two plane vectors is related to their “perp dot product”, where the first vector is rotated 90 degrees counterclockwise. Because both of these are unit vectors, their perp dot product simply is the sine. We want \(s\) positive and don’t care about whether \(\hat{U_2}\) lies to the left or right of \(\hat{U_1}\) so we take the absolute value:

\[ s = | \hat{U}_1^\perp \cdot \hat{U}_2 | \]

Finally, the ratio of the radius \(r\) to \(s\) tells us how far along the vector sum we must go to find the circle’s center:

\[ C = V + \frac{ r }{ s } ( \hat{U}_1 + \hat{U}_2 ) \] \[ \mathopen{\hphantom{C}} = V + r \frac{ \hat{U}_1 + \hat{U}_2 }{ | \hat{U}_1^\perp \cdot \hat{U}_2 | } \]

And that’s it! Fairly elegant in my opinion. Here’s how that looks in code:

void circle_in_angle(
    float   r,
    float P1x, float P1y,
    float  Vx, float  Vy,
    float P2x, float P2y,
    float &Cx, float &Cy )
{
    float VP1 = hypotf( P1x - Vx, P1y - Vy );
    float VP2 = hypotf( P2x - Vx, P2y - Vy );
    float U1x = ( P1x - Vx ) / VP1;
    float U1y = ( P1y - Vy ) / VP1;
    float U2x = ( P2x - Vx ) / VP2;
    float U2y = ( P2y - Vy ) / VP2;
    float s = fabsf( -U1y * U2x + U1x * U2y );
    Cx = Vx + r / s * ( U1x + U2x );
    Cy = Vy + r / s * ( U1y + U2y );
}

Why is this useful? One of the things I’ve been doing for fun over the past year is writing a tiny 2D vector graphics renderer. I was originally trying to implement something like the arcTo method from the HTML5 2D canvas API. If the last point on a sub-path was at \(P_1\) and you call ctx.arcTo( Vx, Vy, P2x, P2y, r ); then it will extend the path with a line and an arc like so:

You can then use a ctx.lineTo( P2x, P2y ) to extend it the rest of the way to \(P_2\) and the result when drawn looks like just drawing a line from \(P_1\) to \(V\) to \(P_2\), but with sharp corner smooothed off into an arc of the given radius. The nice thing is that the canvas method derives the location of the circle center automatically from the other factors.

The other place that I found it useful is in path expansion for stroking lines. If we have line segments on the path from \(P_1\) to \(V\) to \(P_2\), and the line half-width is \(r\), then point \(C\) is the location of the “knee” on the inside of the line join:

If we instead subtract the main term then we get the location of the end of the miter join,

\[ M = V - \frac{ r }{ s } ( \hat{U}_1 + \hat{U}_2 ). \]

Closely related, \(\lVert \hat{U}_1 + \hat{U}_2 \rVert / s\) gives us the miter join ratio which we can compare to the miter limit to determine whether to draw a proper miter join here or truncate it to a bevel join.