Find the 2 coordinates of the 2 control points of a Bezier curve which is an arc of a circle.

Question

I've been searching all day on this topic and I could not figure it out. Hopefully someone is able to dumb it down to my very practical level.

I want to draw an arc. I know the following information:

• Coordinates of Point A: A.x, A.y
• Coordinates of Point B: B.x, B.y
• Coordinate of the center of the circle: 0,0
• The circle Radius: R;

The challenge: I draw my curve with 2 Bezier control points: CA.x, CA.y, CB.x, CB.y. I was able to draw curves for various senarios, but not to find a set of equations that works for any angle in the range of [-360, 360]

Here's where I'm at:

tmp1 = (B.x * A.x) + (B.y * A.y)
tmp2 = ((B.x ^ 2 + B.y ^ 2) * (A.x ^ 2 + A.y ^ 2)) ^ 0.5
RadianAngle = Acos(tmp1 / tmp2)

'Find the directionality of the rotation
D = (A.x * B.y) - (B.x * A.y)
If D < 0 Then
    RadianAngle  = 0 - RadianAngle 
End If

R = 384 ' Calculation code trunkcated for consision

'Up to this point the angle seems to be calculated properly


L = ((4 * R) / 3) * Tan(RadianAngle  / 4)

'Calculate Tengant vectors for points A and B
T1 = Atan2(A.x, A.y)
T2 = Atan2(B.x, B.y)

CA.x = (A.x + L * T1)
CA.y = (A.y + L * T1) 

CB.x = (B.x - L * T2) 
CB.y = (B.y - L * T2) 

Debug values at runtime:

A.x = -384
A.y = 0
B.x = -192
B.y = -332.5537

R = 384

tmp1 = 73727.99
tmp2 = 147456

A = 1.047198 '(60 degrees)
L = 137.19
T1 = 3.141592
T2 = -2.094395

CA.x = 46.99507
CA.y = 430.9951
CB.x = 95.33011
CB.y = -45.22363

I would expect both control points to be negative, most likely somewhere between A and B given there is only 60 degrees between both points.

I've tried to follow the answers:

• from @fang : https://math.stackexchange.com/a/1672700/1146278
• from @robjohn at https://math.stackexchange.com/a/873589/1146278
• from @lemon at https://math.stackexchange.com/a/855541/1146278

... but my success has been very limited and I was not able to implement a working solution

I don't need high precision, this is just to draw a sketch in Excel.

Edit: I found this solution, which seems a lot lenghtier than everything else I see, but it seems to be working ok-ish, until it gets to values near 3/4 of the circle in which case the arc becomes difformed:

https://stackoverflow.com/a/34469229/613365

I'm still be interested in a more robust algorithm.

Answer 1

For this question, I will explain in $\mathbb{R}^3$ at first. To apply to $\mathbb{R}^{2}$ just cut off the third term.

I will divide this answer in two parts:

1. The first one with an exact representation of a circular arc,
2. and the second by approximating the shape from the first part by using Bezier curves

• First part: Generalities on $\mathbb{R}^{3}$:

  Let $\mathbf{A}$ and $\mathbf{B}$ two non-coincident points of a circle on $\mathbb{R}^{3}$ which center is $\mathbf{O} \in \mathbb{R}^{3}$ and radius $R$.

  $$\|\mathbf{A}-\mathbf{O}\| = R \ \ \ \ \ \ \ \ \|\mathbf{B}-\mathbf{O}\| = R$$

  There are two unit vectors $\mathbf{u}$ and $\mathbf{v}$ such that it's possible to describe all the points $\mathbf{Q}$ of a circle.

  $$\mathbf{Q}(\theta) = \mathbf{O} + \mathbf{u} \cdot R\cos \theta + \mathbf{v} \cdot R\sin \theta$$

  For convinience, we say that $\mathbf{Q}(0) = \mathbf{A}$, and $\mathbf{Q}(\theta_0) = \mathbf{B}$

  $$\begin{align}\mathbf{A} = \mathbf{Q}(0) & = \mathbf{O} + \mathbf{u} \cdot R \cdot \cos 0 \\ \mathbf{B} = \mathbf{Q}(\theta_0) & = \mathbf{O} + \mathbf{u} \cdot R \cdot \cos \theta_0 + \mathbf{v} \cdot R \cdot \sin \theta_0\end{align}$$

  The angle $\theta_0$ is the angle between $\mathbf{A}-\mathbf{O}$ and $\mathbf{B}-\mathbf{O}$:

  $$\underbrace{\|\mathbf{A}-\mathbf{O}\|}_{R} \cdot \underbrace{\|\mathbf{B}-\mathbf{O}\|}_{R}\cdot \cos(\theta_0) = \left\langle \mathbf{A}-\mathbf{O}, \mathbf{B}-\mathbf{O} \ \right\rangle$$

  As $\cos \theta$ is an even function, there's no sense of direction: If it's clockwise or counter-clockwise. It also causes confusion cause $\arccos$ function maps $\left[-1, \ 1\right] \to \left[0, \ \pi\right]$

  For the next step, we say the arc will always begin from $\mathbf{A}$ and go to $\mathbf{B}$.

  Example 1: Counter-clockwise

  Let $\mathbf{A} = (1, \ 0, \ 0)$, $\mathbf{B}=\left(\dfrac{1}{2}, \ \dfrac{\sqrt{3}}{2}, \ 0\right)$ and $\mathbf{O}=(0, \ 0, \ 0)$ $$\cos \theta_0 = \dfrac{1}{2}$$ Seeing the plane $xy$, if it's clockwise: $$\theta_0 = \dfrac{5\pi}{3}$$ If it's counter-clockwise: $$\theta_0 = \dfrac{\pi}{3}$$

  Example 2: Clockwise

  Let $\mathbf{A} = (1, \ 0, \ 0)$, $\mathbf{B}=\left(\dfrac{1}{2}, \ \dfrac{-\sqrt{3}}{2}, \ 0\right)$ and $\mathbf{O}=(0, \ 0, \ 0)$ $$\cos \theta_0 = \dfrac{\pi}{3}$$ Seeing the plane $xy$, if it's clockwise: $$\theta_0 = \dfrac{\pi}{3}$$ If it's counter-clockwise: $$\theta_0 = \dfrac{5\pi}{3}$$

  Example 3: 180 degrees

  Let $\mathbf{A} = (1, \ 0, \ 0)$, $\mathbf{B}=\left(-1, \ 0, \ 0\right)$ and $\mathbf{O}=(0, \ 0, \ 0)$ $$\cos \theta_0 = -1 $$ Seeing the plane $xy$, it's $\theta_0 = \pi$ clockwise or counter-clockwise

  To specify it, we use the vector $\vec{n}$ which relates to the 'axis' or the circle, and the vectors $\vec{u}$ and $\vec{v}$.

  $$\mathbf{u} = \dfrac{\mathbf{A}-\mathbf{O}}{\|\mathbf{A}-\mathbf{O}\|}$$ $$\mathbf{V} = (\mathbf{B}-\mathbf{O}) - \mathbf{u} \cdot \langle \mathbf{B}-\mathbf{O}, \ \mathbf{u}\rangle$$ $$\mathbf{v} = \dfrac{\mathbf{V}}{\|\mathbf{V}\|}$$ $$\mathbf{n} = \mathbf{u} \times \mathbf{v}$$

  In $\mathbb{R}^2$, if $\mathbf{n} = (0, \ 0, \ 1)$ then it's counter-clockwise, if $\mathbf{n} = (0, \ 0, \ -1)$, then it's clockwise.

  Then, to draw an arc, one can do

  • Define $\mathbf{n}$, unit vector, to say if it's clockwise or counter-clockwise
  • Compute $\mathbf{U} = \mathbf{A}-\mathbf{O}$
  • Compute the radius $R$ by taking the norm $L_2$ of $\mathbf{U}$
  • Compute $\mathbf{u}$ by diving $\mathbf{U}$ by the radius $R$
  • Compute $\mathbf{v} = \mathbf{n} \times \mathbf{u}$
  • Compute $\theta_0 = \arctan_{2}(\langle \mathbf{B}-\mathbf{O}, \ \mathbf{v}\rangle, \ \langle \mathbf{B}-\mathbf{O}, \ \mathbf{u}\rangle)$

  Take care if $\theta_0$ is negative. It's expected to get $\theta_0 \in \left[0, \ 2\pi\right)$

  • For every $\theta \in \left[0, \ \theta_0\right]$ compute the point $\mathbf{Q}$: $$\mathbf{Q}(\theta) = \mathbf{O} + \mathbf{u} \cdot R\cos \theta + \mathbf{v} \cdot R\sin \theta$$

  import numpy as np
  from matplotlib import pyplot as plt
  
  A = np.array([1., 0., 0.])
  B = np.array([0., 1., 0.])
  O = np.array([0., 0., 0.])
  n = np.array([0., 0., -1.])
  
  U = A-O
  R = np.linalg.norm(U)
  u = U/R
  v = np.cross(n, u)
  t0 = np.arctan2(np.inner(B-O, v), np.inner(B-O, u))
  if t0 < 0:
       t0 += 2*np.pi
  print("u = ", u)
  print("v = ", v)
  print("n = ", n)
  print(f"t0 = {t0:.2f} rad = {(180*t0/np.pi):.1f} deg")
  
  theta = np.linspace(0, t0, 1025)
  px = O[0] + u[0]*R*np.cos(theta) + v[0]*R*np.sin(theta)
  py = O[1] + u[1]*R*np.cos(theta) + v[1]*R*np.sin(theta)
  plt.plot(px, py)
  plt.axis("equal")
  plt.show()
  

• Second part: Bezier curve:

  A Bezier curve $\mathbf{C}$ is given by a parameter $u \in \left[0, \ 1\right]$ and by $(n+1)$ control points $\mathbf{P}_i$

  $$\mathbf{C}(u) = \sum_{i=0}^{n} B_{i}(u) \cdot \mathbf{P}_i$$ $$B_i(u) = \binom{n}{i} \left(1-u\right)^{n-i} \cdot u^{i}$$

  There is no way to describe a circular path by using Bezier curves. To do it exactly, it's necessary to use Rational polynomial functions (NURBS). For example, a $1/4$ circle is given by

  $$\begin{align}\mathbf{C}(u) & = \left(\dfrac{1-u^2}{1+u^2}, \ \dfrac{2u}{1+u^2}, \ 0\right) = \sum_{i=0}^{2} R_{i2}(u) \cdot \mathbf{P}_i \\ & = \underbrace{\dfrac{(1-u)^2}{1+u^2}}_{R_{02}} \cdot \underbrace{\left(1, \ 0, \ 0\right)}_{\mathbf{P}_0} + \underbrace{\dfrac{2u(1-u)}{1+u^2}}_{R_{12}} \cdot \underbrace{\left(1, \ 1, \ 0\right)}_{\mathbf{P}_1} + \underbrace{\dfrac{2u^2}{1+u^2}}_{R_{22}} \cdot \underbrace{\left(0, \ 1, \ 0\right) }_{\mathbf{P}_0}\end{align}$$

  Although there's no way to do it by using Bezier curves, we can get an approximate shape. The question becomes how to get the $n+1$ control point of Bezier curve.

  The first information is the curve $\mathbf{C}$'s extremities must be $\mathbf{A}$ and $\mathbf{B}$:

  $$\begin{align}\mathbf{A} & = \mathbf{C}(0) \Rightarrow \mathbf{P}_0 = \mathbf{A} \\ \mathbf{B} & = \mathbf{C}(1) \Rightarrow \mathbf{P}_{n} = \mathbf{B}\end{align}$$

  and set

  $$\theta = u \cdot \theta_0$$

  1. One way of doing it is by computing the tangent vector at the extremities when $n=3$:

  $$\begin{align}\left[\dfrac{d\mathbf{C}(u)}{du}\right]_{u=0} & = 3\left(\mathbf{P}_1 - \mathbf{P}_0\right) \\ \left[\dfrac{d\mathbf{C}(u)}{du}\right]_{u=1} & = 3\left(\mathbf{P}_3 - \mathbf{P}_2\right) \end{align}$$ $$\begin{align}\left[\dfrac{d\mathbf{Q}(\theta)}{d\theta}\right]_{\theta=0} & = R \cdot \mathbf{v} \\ \left[\dfrac{d\mathbf{Q}(\theta)}{d\theta}\right]_{\theta=\theta_0} & = R \cdot \mathbf{v} \cos \theta_0 - R \cdot \mathbf{u} \cdot \sin \theta_0\end{align}$$

  Now set $$\begin{align}\mathbf{C}'(0) & = \theta_0 \cdot \mathbf{Q}'(0) \\ \mathbf{C}'(1) & = \theta_0 \cdot \mathbf{Q}'(\theta_0)\end{align}$$

  To get

  $$\begin{align}\mathbf{P}_1 & = \mathbf{A} + \dfrac{\theta_0 R}{3} \cdot \mathbf{v} \\ \mathbf{P}_2 & = \mathbf{B} + \dfrac{\theta_0 R}{3}\left(\mathbf{u} \cdot \sin \theta_0 - \mathbf{v} \cdot \cos \theta_0 \right) \end{align}$$

  2. Using the integral by least square

  Let $\mathbf{D}(u)$ be the distance between the curve $\mathbf{Q}(\theta) = \mathbf{Q}(\theta_0 \cdot u)$ and $\mathbf{C}(u)$. Then, the objective is to reduce the function $J(\mathbf{P}_i)$:

  $$J(\mathbf{P}) = \int_{0}^{1} \|\mathbf{D}\|^2 \ du = \int_{0}^{1} \ \langle \mathbf{D}, \ \mathbf{D}\rangle \ du$$

  Then getting the system of equations $\dfrac{\partial J}{\partial \mathbf{P}_i} = 0 \ \ \forall i$ and solving it.

I made a python code to plot the result for $\mathbf{A}=(1, \ 0, \ 0)$, $\mathbf{B}=\left(\frac{-1}{2}, \ \frac{\sqrt{3}}{2}, \ 0\right)$, $\mathbf{C} = \left(0, \ 0, \ 0\right)$ and $\mathbf{n}=\left(0, \ 0, \ -1\right)$

import numpy as np
from matplotlib import pyplot as plt

# Initial configuration
A = np.array([1., 0., 0.])
B = np.array([-0.5, np.sqrt(3)/2, 0.])
O = np.array([0., 0., 0.])
n = np.array([0., 0., -1.])

# Computing basics
U = A-O
R = np.linalg.norm(U)
u = U/R
v = np.cross(n, u)
t0 = np.arctan2(np.inner(B-O, v), np.inner(B-O, u))
if t0 < 0:
    t0 += 2*np.pi

# Computing exact curve
theta = np.linspace(0, t0, 1025)
px = O[0] + u[0]*R*np.cos(theta) + v[0]*R*np.sin(theta)
py = O[1] + u[1]*R*np.cos(theta) + v[1]*R*np.sin(theta)
plt.plot(px, py, label="exact")

# Computing bezier curve, by tangents
P = np.array([A, A + t0*R*v/3, B + (t0*R/3)*(u*np.sin(t0)-v*np.cos(t0)), B])
uplot = np.linspace(0, 1, 1025)
Bez = [(1-uplot)**3, 3*uplot*(1-uplot)**2, 3*(1-uplot)*uplot**2, uplot**3]
Cx = sum([Bez[i]*P[i, 0] for i in range(4)])
Cy = sum([Bez[i]*P[i, 1] for i in range(4)])
plt.plot(Cx, Cy, label="bezier")

# Showing results
plt.legend()
plt.axis("equal")
plt.show()

EDIT:

For the 2D case:

As said before, for the 2D case we can only ignore the third component. Then, the only difference is on computing the $\mathbf{v}$:

• If clockwise: $$\mathbf{v}=(u_y, \ -u_x)$$

• If counter-clockwise: $$\mathbf{v}=(-u_y, \ u_x)$$

Answer 2

Alright, this answer is clearly more about VBA programming than Math per say, but given Carlos provided a really detailed answer, I felt it would be a good contribution regardless for future users coming from Google.

While the implementation of Carlos's methodology worked great, I had trouble with some edge cases where, I think (still investigating), negative rotation above 90 degrees cause issues. I'm sure it is my implementation which lack finesse, but I have also been working on implementing another approach, which at this point in time seems more reliable.

So I figure I'd share it, even if it's significantly messier.

Note that I'm still working towards drawing rotation from -359 to +359. The application is basically calculating an object rotation on a pivot hinge, so arcs greater than 180 will need to be drawn. Which is why the code has some ToDo's in it and use the concept of Segments. Currently the code will draw the smallest of the 2 possible arcs, which 180 is the maximum you can do with: a starting point, an end point and 2 control points. Future implementations aims to allow for multiple segments.

Both functions make use the following custom types and constants:

Const PI As Double = 3.14159265358979
Const DEG_TO_RAD As Double = PI / 180#
Const RAD_TO_DEG As Double = 180# / PI

Private Type Coordinates
    X As Single
    Y As Single
End Type

Private Type BezierSegment
    Vertex As Coordinates
    ControlPoint1 As Coordinates
    ControlPoint2 As Coordinates
End Type

So here's one implementation which seems to give the best results so far:

Private Function BezierPointsArrayBuilder4( _
        StartPoint As Coordinates, _
        EndPoint As Coordinates, _
        Axis As Coordinates, _
        ByRef ArrBezierSegments As Variant, _
        Optional n As Integer = 1 _
) As Coordinates 'Return the coordinate of the startpoint (which will vary if the rotation angle is positive or negative

    Dim t1x As Double
    Dim t1y As Double
    Dim t2x As Double
    Dim t2y As Double
    Dim dx As Double
    Dim dy As Double
    Dim k As Double
    Dim Tx As Double
    Dim Ty As Double
    Dim D As Double
    Dim a As Double
    Dim B As Double
    Dim C As Double
    
    Dim Segment As BezierSegment
    Set Segment = New BezierSegment
    
    Dim Radius As Single
    Radius = Axis.DistanceFromCoordinate(EndPoint)

    Dim Angle As Single
    'Currently this is calculated, however the angle is a user-fed data, which is why I'm looking for the range of -359 to +359
    Angle = Axis.AngleBetweenCoordinates(StartPoint, EndPoint)

    t1x = Axis.Y - StartPoint.Y
    t1y = StartPoint.X - Axis.X
    t2x = EndPoint.Y - Axis.Y
    t2y = Axis.X - EndPoint.X
    dx = (StartPoint.X + EndPoint.X) / 2 - Axis.X
    dy = (StartPoint.Y + EndPoint.Y) / 2 - Axis.Y
    Tx = 3 / 8 * (t1x + t2x)
    Ty = 3 / 8 * (t1y + t2y)
    a = Tx * Tx + Ty * Ty
    B = dx * Tx + dy * Ty
    C = dx * dx + dy * dy - Radius * Radius
    D = B * B - a * C

    If D < 0 Then
        Err.Raise vbObjectError + 1000, , "Unable to calculate the curve as D < 0"
        Exit Function
    End If
    
    k = Round((D ^ 0.5 - B) / a, 4)
    If a = 0 Then 
        
        Err.Raise vbObjectError + 1000, , "This arc can't be drawn."
        Exit Function
        
    ElseIf Angle > 0 And Angle <= 180 Then
        
        Segment.ControlPoint1.X = StartPoint.X + Round(k * t1x)
        Segment.ControlPoint1.Y = StartPoint.Y + Round(k * t1y)
        
        Segment.ControlPoint2.X = EndPoint.X + Round(k * t2x)
        Segment.ControlPoint2.Y = EndPoint.Y + Round(k * t2y)
        
        Segment.Vertex.X = EndPoint.X
        Segment.Vertex.Y = EndPoint.Y
        
        ReDim Preserve ArrBezierSegments(n)
        Set ArrBezierSegments(n) = Segment
        Set BezierPointsArrayBuilder4 = EndPoint
    
    ElseIf Angle > 180 And Angle < 360 Then
        'TODO : Redim the array and recursively call this function with n+1
    ElseIf Angle < 0 And Angle >= -180 Then
    
        Segment.ControlPoint1.X = EndPoint.X + Round(k * t2x)
        Segment.ControlPoint1.Y = EndPoint.Y + Round(k * t2y)
        
        Segment.ControlPoint2.X = StartPoint.X + Round(k * t1x)
        Segment.ControlPoint2.Y = StartPoint.Y + Round(k * t1y)
        
        Segment.Vertex.X = StartPoint.X
        Segment.Vertex.Y = StartPoint.Y
        
        ReDim Preserve ArrBezierSegments(n)
        Set ArrBezierSegments(n) = Segment
        Set BezierPointsArrayBuilder4 = StartPoint
        
    ElseIf Angle < -180 And Angle > -360 Then
        'TODO : Redim the array and recursively call this function with n+1
    End If
    
    

End Function

For the record, my current implementation of Carlos methodogy is as follow:

Private Function BezierPointsArrayBuilder3( _
        StartPoint As Coordinates, _
        EndPoint As Coordinates, _
        Axis As Coordinates, _
        ByRef ArrBezierSegments As Variant, _
        Optional n As Integer = 1 _
) As Coordinates 'Return the coord
    Dim F As Double 
    Dim R As Single 'Radius
    Dim A_deg As Single
    Dim A_rad As Single
    
    'Currently this is calculated, however the angle is a user-fed data, which is why I'm looking for the range of -359 to +359
    A_deg = Axis.AngleBetweenCoordinates(StartPoint, EndPoint, False)
    A_rad = A_deg * DEG_TO_RAD
    F = 4 / 3 * Tan(A_rad / 4)
    R = Axis.DistanceFromCoordinate(StartPoint) ' Find the distance between the axis and the a point on the arc
    
    Dim U, V, Q As Coordinates
    Set U = New Coordinates
    Set V = New Coordinates
    Set Q = New Coordinates
    
    
    U.X = StartPoint.X - Axis.X
    U.Y = StartPoint.Y - Axis.Y
    
    V.X = U.Y
    V.Y = 0 - U.X
    
    
    Q.X = Axis.X + (U.X / R) * R * Cos(A_rad) + V.X * R * Sin(A_rad)
    Q.Y = Axis.Y + (U.Y / R) * R * Cos(A_rad) + V.Y * R * Sin(A_rad)
    
    Dim i As Integer
    i = UBound(ArrBezierSegments)
    If i < 0 Then 'Non-intanciated array, position the first index at 1
        i = 0
    End If
    
    ReDim Preserve ArrBezierSegments(i + 1)
    
    Dim Segment As BezierSegment
    Dim t1, t2 As Single
    Set Segment = New BezierSegment
    
    Segment.ControlPoint1.X = Axis.X - R
    Segment.ControlPoint1.Y = Axis.Y - R * F
    
    t1 = R * (Cos(A_rad) + F * Sin(A_rad))
    t2 = R * (Sin(A_rad) - F * Cos(A_rad))
    Segment.ControlPoint2.X = Axis.X - t1
    Segment.ControlPoint2.Y = Axis.Y - t2
    
    'The following code does not work.
    t1 = Cos(A_rad)
    t2 = Sin(A_rad)
    Segment.Vertex.X = Axis.X - R * t2
    Segment.Vertex.Y = Axis.Y - R * t1
    
    'This lazy alternative works, but is problematic to link multiple segments.
    Segment.Vertex.X = EndPoint.X 
    Segment.Vertex.Y = EndPoint.Y
        
    Set ArrBezierSegments(i + 1) = Segment
    Set BezierPointsArrayBuilder3 = Segment.Vertex
End Function