1
$\begingroup$

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).

enter image description here

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 !

$\endgroup$
2
  • $\begingroup$ As written this isn't quite enough information: Is $B$ supposed to be the point on the curve which is closest to $A$? (You could also imagine criteria like "snap to the point on the curve which is directly above/below A".) $\endgroup$ Commented Feb 8, 2023 at 16:58
  • $\begingroup$ Thanks for your feedback. Yes, ideally the point on the curve the closest to A. But less ideal solutions are welcome too. $\endgroup$ Commented Feb 8, 2023 at 17:04

2 Answers 2

1
$\begingroup$

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.

enter image description here

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
$\endgroup$
8
  • $\begingroup$ Thanks very much for your reply. However I have some issues interpreting your reply to get to the solution I seek. How to calculate 'a' and 'b' ? I have the coordinates of : centerPoint, majorAxis and minorAxis. Also could you perhaps write me the algorithm in a for loop? I added to the question a template of such algorithm but it isn't correct. Thanks! $\endgroup$ Commented Feb 9, 2023 at 9:26
  • $\begingroup$ $a$ and $b$ should be known to you, they are the axes of the hyperbola that I assume is given. $\endgroup$ Commented Feb 9, 2023 at 9:34
  • $\begingroup$ I will write the algorithm to find $x_2$ and $y_2$ shortly. $\endgroup$ Commented Feb 9, 2023 at 9:37
  • $\begingroup$ I've included a listing of the program to find the snap point $(x_2,y_2)$ given the axes $a$ and $b$ and the initial point $(x_1, y_1)$. $\endgroup$ Commented Feb 9, 2023 at 10:17
  • $\begingroup$ Thanks very much ! It is now working as expected. For reference x1 and y1 needs to be in the hyperbola coordinate system. I will post the c++ code here too for reference for future readers. $\endgroup$ Commented Feb 9, 2023 at 13:16
1
$\begingroup$

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;
}
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .