Calculate a point on an parabola from a nearby point.

Question

For an open source cad software called FreeCAD. In 2d editor I'm making the snapping feature for parabolas curves.

There is an parabola, the user put the mouse pointer close to it. So we have the point A, which is close to the parabola, but not exactly on it. I need to calculate the coordinates of the point B which is on the parabola close to A. (then the point B is used to override the mouse coordinates, such that the parabolais exactly under the mouse pointer).

Do you know how to calculate the point B ?

Thanks !

PS: Here is the same question but for hyperbola: Calculate a point on an hyperbola from a nearby point.

Answer 1

The first step is to generate the implicit equation of the parabola, i.e. find the equation in the form

$ A x^2 + B x y + C y^2 + D x + E y + F = 0 $

To do that, starting from the vertex $V$ and the focus $F$, compute the unit vector

$ u_2 = \dfrac{F - V}{\| F - V \| } $

And calculate the focal distance $p = \| F - V \| $

And finally calculate a unit vector $u_1$ that is perpendicular to $u_2$. Then the parametric equation of the parabola is

$ r = V + t u_1 + \dfrac{1}{4 p} t^2 u_2 $

The vector $u_1$ is computed from $u_2$ as follows:

$u_{1x} = u_{2 y} $

$u_{1y} = - u_{2 x} $

Let

$(x', y') = t u_1 + \dfrac{1}{4p} t^2 u_2 = [u_1, u_2] [ t, \dfrac{1}{4p} t^2]^T $

Let $U = [u_1, u_2]$, note that since $u_1, u_2$ are unit vectors then $U$ is a rotation matrix.

So now we have

$ (x', y') = U (x'', y'') = U [t, \dfrac{1}{4p} t^2 ]^T $

The equation governing $x''$ and $y''$ is $ y'' = \dfrac{1}{4p} x''^2 $

So

$ [x'', y''] Q \begin{bmatrix} x'' \\ y'' \end{bmatrix} + \begin{bmatrix} 0 && - 1 \end{bmatrix} \begin{bmatrix} x'' \\ y'' \end{bmatrix} = 0 $

where $Q = \begin{bmatrix} \dfrac{1}{4p} && 0 \\ 0 && 0 \end{bmatrix} $

But $[x'', y'']^T = U^{-1} [x', y']^T $, so

$ (r - V)^T U Q U^T (r - V) + b^T U^T (r - V) = 0 $

where $b = [0, -1]^T $

This last equation is the equation of the parabola in general format.

Now, we begin finding the closest point $(x_2, y_2)$ on the parabola to a given point $(x_1, y_1)$ that is not on the parabola.

The first equation that $(x_2, y_2)$ has to satisfy is

$ A x_2^2 + B x_2 y_2 + C y_2^2 + D x_2 + E y_2 + F = 0 $

The gradient of this function of two variables is given by

$ \nabla = \begin{bmatrix} 2 A x_2 + B y_2 + D \\ B x_2 + 2 C y_2 + E \end{bmatrix} $

And we want this gradient to have the same direction as $(x_2 - x_1, y_2 - y_1) $

So

$( 2 A x_2 + B y_2 + D , B x_2 + 2 C y_2 + E ) = \alpha (x_2 - x_1, y_2 - y_1 ) $

Eliminating $\alpha$ by dividing the $y$ coordinate by the $x$ coordinate,

$ \dfrac{y_2 - y_1}{x_2 - x_1} = \dfrac{ B x_2 + 2 C y_2 + E }{2 A x_2 + B y_2 + D } $

Cross multiplying,

$ (y_2 - y_1)(2 A x_2 + B y_2 + D) = (x_2 - x_1) (B x_2 + 2 C y_2 + E ) $

Now we have two quadratic equations in the two unknowns $x_2$ and $y_2$. These can be found using the Newton-Raphson root finder.

The following figure shows a sample parabola with a given point off the parabola, and the found point that is closest to it on the parabola.

Below, I've listed the whole program (written in Visual Basic for Applications)

Public Sub solve_2_variable_quadratic_system_2()


Dim xvec(2), yvec(2), jac(2, 2), jacinv(2, 2) As Double
Dim det As Double
Dim dvec(2) As Double
Dim v(3), focus(3) As Double
Dim q(2, 2), u(2, 2), ut(2, 2) As Double
Dim pvec(6) As Double
Dim b0(3) As Double
Dim a, b, c, d, e, f As Double


v(1) = 5
v(2) = 7


focus(1) = 5.5
focus(2) = 6.5


u2(1) = focus(1) - v(1)
u2(2) = focus(2) - v(2)

p = norm(u2, 2)

u2(1) = u2(1) / p
u2(2) = u2(2) / p

u1(1) = u2(2)
u1(2) = -u2(1)

For i = 1 To 2
u(i, 1) = u1(i)
u(i, 2) = u2(i)
Next i

q(1, 1) = 1 / (4 * p)
q(1, 2) = 0
q(2, 1) = 0
q(2, 2) = 0

b0(1) = 0
b0(2) = -1

Call transpose(2, 2, u, ut)

Call multiply(2, 2, 2, u, q, q)
Call multiply(2, 2, 2, q, ut, q)

Call multiplyrv(2, 2, b0, ut, b0)

' now the equation is  (r - V)^T Q (r - V) + b^T (r - V) = 0

a = q(1, 1)

b = 2 * q(1, 2)

c = q(2, 2)

Call multiplyv(2, 2, q, v, u3)

d = b0(1) - 2 * u3(1)

e = b0(2) - 2 * u3(2)

f = xtqy(2, q, v, v) - dot(b0, v)

pvec(1) = a
pvec(2) = b
pvec(3) = c
pvec(4) = d
pvec(5) = e
pvec(6) = f


x1 = 10
y1 = 4

'------------------- Root finding starts here -------

xvec(1) = x1
xvec(2) = y1

success = 0

For ic = 1 To 20

' display current estimate of root

For i = 1 To 2
ActiveSheet.Cells(ic, i) = xvec(i)
Next i

' Compute the nonlinear functions yvec

Call find_yvec_2(pvec, x1, y1, xvec, yvec)

' norm(yvec, 2) returns the length of yvec, which is SQRT( yvec(1)^2 + yvec(2)^2 )

If norm(yvec, 2) < 0.0000001 Then   ' Tol is 1.E-7

  success = 1
  Exit For

End If

' Compute the jacobian

Call find_jac_2(pvec, x1, y1, xvec, jac)

' find the inverse of the jacobian matrix

det = jac(1, 1) * jac(2, 2) - jac(1, 2) * jac(2, 1)

jacinv(1, 1) = jac(2, 2) / det
jacinv(2, 2) = jac(1, 1) / det

jacinv(1, 2) = -jac(1, 2) / det
jacinv(2, 1) = -jac(2, 1) / det

' multiply jacinv by  yvec

dvec(1) = jacinv(1, 1) * yvec(1) + jacinv(1, 2) * yvec(2)
dvec(2) = jacinv(2, 1) * yvec(1) + jacinv(2, 2) * yvec(2)

' if norm(dvec,2) < 1.e-7 Then exit

If norm(dvec, 2) < 0.0000001 Then

 success = 1
 Exit For

End If

' update xvec

xvec(1) = xvec(1) - dvec(1)
xvec(2) = xvec(2) - dvec(2)


Next ic

If success = 0 Then
   MsgBox ("Did not converge")

Else

   x2 = xvec(1)
   y2 = xvec(2)

   MsgBox ("Converged in " + Str(ic - 1) + " iterations")
   MsgBox ("Snap point is " + Str(x2) + Str(y2))

 End If


 End Sub
 Public Sub find_yvec_2(pvec, x1, y1, xvec, yvec)

Dim x, y As Double
Dim a, b, c, d, e, f As Double

a = pvec(1)
b = pvec(2)
c = pvec(3)
d = pvec(4)
e = pvec(5)
f = pvec(6)

x = xvec(1)
y = xvec(2)

yvec(1) = a * x ^ 2 + b * x * y + c * y ^ 2 + d * x + e * y + f
yvec(2) = (y - y1) * (2 * a * x + b * y + d) - (x - x1) * (b * x + 2 * c * y + e)


End Sub
Public Sub find_jac_2(pvec, x1, y1, xvec, jac)

Dim x, y As Double
Dim a, b, c, d, e, f As Double

a = pvec(1)
b = pvec(2)
c = pvec(3)
d = pvec(4)
e = pvec(5)
f = pvec(6)

x = xvec(1)
y = xvec(2)

jac(1, 1) = 2 * a * x + b * y + d
jac(1, 2) = b * x + 2 * c * y + e

jac(2, 1) = (y - y1) * (2 * a) - (b * x + 2 * c * y + e) - (x - x1) * b

jac(2, 2) = (2 * a * x + b * y + d) + (y - y1) * b - (x - x1) * (2 * c)

End Sub

Answer 2

Let us consider the case of the reference parabola with equations :

$$x=t, \ \ y=at^2 \ \ \iff \ \ y=ax^2\tag{1}$$

Fig. 1 : Different points $A$ (red circles) and their "projections" $H$ (blue circles) onto the parabola.

I propose to replace the point $A$ "clicked" by the user by point $H$ on the parabola as indicated on the figure, that can be described as being "almost" the foot of the perpendicular line drawn from $A$ to line $BC$ in right triangle $ABC$ ; in fact, the real point $H$ is the (slightly different) point obtained as the intersection of this perpendicular line and the parabola.

The different steps are as follows :

Let $A(u,v)$ the point clicked by the user.

• $B$ has coordinates $(u,f(u))$.

• $C$ has coordinates $(\sqrt{v/a},v)$ (check that (1)) is verified).

As the slope of line $BC$ is $s=\frac{v-f(u)}{\sqrt{v/a}-u}$, the slope of the orthogonal line to line $BC$ passing through $A$ is $-\frac{1}{s}$. The equation of this line is therefore :

$$y-v=\underbrace{-\frac{1}{s}}_{slope}(x-u)\tag{2}$$

The determination of the abscissas of the intersection point of this line with the parabola is obtained by plugging (1) into (2).

The result is a quadratic equation whose (well chosen !) solution provides the abscissa of point $H$.

This has to be slightly modified for the case of negative abscissas (see Matlab program below) by detecting the signs of the "protagonists".

The general case is obtained by using a rotation-translation bringing the standard axes we have used above onto the "real axes".

 a=0.5;t=-2:0.01:2;
 f=@(t)(a*t.^2); % function f(t)=at^2
 plot(t,f(t)); % parabola
 for i=1:20
    u=3*rand-1.5;v=f(u)+0.6*rand-0.3; % random points A
    w=sign(u)*sqrt(v/a); % absc. of C
    k=(u-w)/(v-f(u)); % slope of AH
    d=k^2-4*a*(k*u-v); % discriminant of quadratic
    x=(k+sign(u)*sqrt(d))/(2*a); % root of quadratic = absc. of H
    plot([u,x],[v,f(x)]); % line segment AH
 end;