Calculate a point on an hyperbola from a nearby point.

Question

I am working on an open source cad software called FreeCAD. In the sketcher I'm making the snapping feature for hyperbolas curves.

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

Base::Vector2d A, centerPoint,  majorAxis, minorAxis are known.

//This is appantly wrong : 
Base::Vector2d vec = A - centerPoint;
double U = atanh(vec.y / minorAxis.Length()) + acosh(vec.x / majorAxis.Length());

//This seems to be OK as it create points on the hyperbola :
B = centerPoint + majorAxis * cosh(U) + minorAxis * sinh(U);

The last line calculates points that are on the hyperbola. So it seems OK, if so the U that is fed into would be wrong.

Do you know how to calculate the point B ? Thanks !

Answer 1

Let the given hyperbola be

$ \dfrac{x^2}{a^2} - \dfrac{y^2}{b^2} = 1 $

And assume that $A = (x_1, y_1) $. We want to find the closest point $B$ on the hyperbola to point $A$. Let $B = (x_2, y_2) $, then we can write two equations for $x_2$ and $y_2$

$\dfrac{x_2^2}{a^2} - \dfrac{y_2^2}{b^2} = 1 \hspace{25pt}(1) $

and

$ ( x_2 - x_1 , y_2 - y_1 ) = \alpha \vec{n} = \alpha ( \dfrac{2 x_2}{a^2} , \dfrac{ -2 y_2}{b^2} ) \hspace{25pt}(2) $

The second equation has an extra variable $\alpha$ that we can eliminate by dividing the $y$ coordinate by the $x$ coordinate. And this gives us

$ \dfrac{ y_2 - y_1 }{ x_2 - x_1 } = - \dfrac{ a^2 y_2 }{ b^2 x_2 } \hspace{25pt} (3)$

Cross multiplying, yields,

$ b^2 (y_2 - y_1) x_2 = - a^2 (x_2 - x_1) y_2 \hspace{25pt}(4)$

Equations $(1)$ and $(4)$ are two quadratic equations in the two unknowns $x_2$ and $y_2$, and can be solved with the help of a math app like Mathematica for example. In the absence of a math app, you can solve numerically, using the multivariate Newton-Raphson root finder algorithm, which is very fast and converges to a solution (from the initial guess of $x_2 = x_1$ and $y_2 = y_1$) with very few iterations (about $4$ or $5$ iterations).

The following figure shows the hyperbola with the initial point off the parabola and the closest point to it on the hyperbola.

In the following, I list the Visual Basic Code to implement the Newton-Raphson root finder for this specific problem.

Public Sub solve_2_variable_quadratic_system()

' system of nonlinear (quadratic) equations given by

' x^2/a^2 - y^2/b^2 = 1

' a^2 (y - y1) y = - b^2 (x - x1) x

Dim xvec(2), yvec(2), jac(2, 2), jacinv(2, 2) As Double
Dim det As Double
Dim dvec(2) As Double

a = 5
b = 2

x1 = 9
y1 = 2.5

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

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

success = 0

For ic = 1 To 100

' Compute the nonlinear functions yvec

Call find_yvec(a, b, 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(a, b, 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(a, b, x1, y1, xvec, yvec)

Dim x, y As Double

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

yvec(1) = x ^ 2 / a ^ 2 - y ^ 2 / b ^ 2 - 1
yvec(2) = b ^ 2 * (y - y1) * x + a ^ 2 * (x - x1) * y


End Sub
Public Sub find_jac(a, b, x1, y1, xvec, jac)

Dim x, y As Double

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

jac(1, 1) = 2 * x / a ^ 2
jac(1, 2) = -2 * y / b ^ 2

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

End Sub

Answer 2

Big thanks to Glorious Erin to solve this. Here is the code in C++ for reference for future readers :

const Part::GeomArcOfHyperbola* hyperbola = static_cast<const Part::GeomArcOfHyperbola*>(geo);
Base::Vector2d centerPoint = Base::Vector2d(hyperbola->getCenter().x, hyperbola->getCenter().y);
Base::Vector2d majorAxis = Base::Vector2d(hyperbola->getMajorAxisDir().x, hyperbola->getMajorAxisDir().y);
Base::Vector2d minorAxis = Base::Vector2d(hyperbola->getMajorAxisDir().y, - hyperbola->getMajorAxisDir().x);

Base::Vector2d xvec, yvec, dvec;
double jac[2][2], jacinv[2][2];
double det;
double a = hyperbola->getMajorRadius();
double b = hyperbola->getMinorRadius();

//First we need to know the coordinates of pointToOverride in the coordinate system of the hyperbola.
Base::Vector2d p = pointToOverride - centerPoint;
Base::Vector2d px, py;
px.ProjectToLine(p, majorAxis);
py.ProjectToLine(p, minorAxis);
bool reverseY = py * minorAxis < 0;
p = Base::Vector2d(px.Length(), py.Length());

//Initialization xvec
xvec = p;

bool success = false;
for (int i = 0; i < 100; ++i) {  // Maximum 100 iterations
    //Compute the nonlinear functions yvec
    yvec.x = xvec.x * xvec.x / (a * a) - xvec.y * xvec.y / (b * b) - 1;
    yvec.y = a * a * (xvec.y - p.y) * xvec.y + b * b * (xvec.x - p.x) * xvec.x;

    if (yvec.Length() < 1e-7) {
        success = true;
        break;
    }

    //Compute the jacobian
    jac[0][0] = 2 * xvec.x / (a * a);
    jac[0][1] = - 2 * xvec.y / (b * b);
    jac[1][0] = a * a * (2 * xvec.y - p.y);
    jac[1][1] = b * b * (2 * xvec.x - p.x);

    //find the inverse of the jacobian matrix
    det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
    jacinv[0][0] = jac[1][1] / det;
    jacinv[1][1] = jac[0][0] / det;
    jacinv[0][1] = -jac[0][1] / det;
    jacinv[1][0] = -jac[1][0] / det;

    //multiply jacinv by  yvec
    dvec.x = jacinv[0][0] * yvec.x + jacinv[0][1] * yvec.y;
    dvec.y = jacinv[1][0] * yvec.x + jacinv[1][1] * yvec.y;

    if (dvec.Length() < 1e-7) {
        success = true;
        break;
    }

    //update xvec
    xvec = xvec - dvec;
}

if (success) {
    pointToOverride = centerPoint + xvec.x * majorAxis + (reverseY ? -1 : 1) * xvec.y * minorAxis;
}