I have coordinates of a set of points that form a closed loop that lies in a 3D surface. I know the equation of the surface and I can calculate it's surface normal at any point. I found that for a unit speed curve $\alpha(t)$ I can calculate the geodesic curvature as $$ \kappa_g(t) = \alpha''(t)\cdot (n(t)\times\alpha'(t)) $$

where $n(t)$ is the surface normal. I am confused about how to calculate the unit-speed derivatives $\alpha'(t)$ and $\alpha''(t)$. I think I should use some finite difference scheme like $$ \alpha'(t) \approx \frac{\alpha(t + h) - \alpha(t - h)}{2h} $$

So my algorithm for calculating the geodesic curvature at any point $P$ of this curve should be as follows:

  1. Identify a reference point on the curve to measure arc-length from.
  2. For the point $P$ calculate the arc-length by adding up the Euclidean distances between successive points from the reference point to the point P. This arc length gives me $t$.
  3. Find the position vectors of points at $t + h$ and $t - h$ by interpolation and then use the above finite difference equation to calculate $\alpha'(t)$.

Is my algorithm correct? I am not clear about how to get arc-length parameterization from the points on the curve.

Are there any standard algorithms (or even Python libraries) for such a calculation?

UPDATE 1: I implemented @Futurologist 's wonderful answer for the test case of a circle on a cone as shown below

Cone with circle

I get the Frenet-Serret frame along the circle as shown below

TNB Frame

The total curvature at each point on the circle is $k_i = 1.0$ which is correct because I designed the circle to have a radius of $1.0$. I get the geodesic curvature at all points as $-1/\sqrt{2}$. So everything looks good. Can someone confirm if the sign of the geodesic curvature or the $N_i$ vectors in the TNB frames are correct? I know that I can flip it by changing the order in the cross-product for $b_i$ but I am asking if there are is standard convention.

UPDATE 2: I had to use $b_i = (\alpha_{i} - \alpha_{i-1})\times(\alpha_{i+1} - \alpha_i)$ to get the correct sign for $N_i$ and a right-handed TNB frame. Also at some of the points the tangent $T_i$ was pointing in the opposite direction of the increasing arc-length parameter due to angle between $b_i$ and $n_i$ increasing above $\pi$. In those cases, I had to flip the $T_i$. Now I get the geodesic curvature at all points as $1/\sqrt{2}$ and the vectors $N_i$ point inwards as desired.


1 Answer 1


Why don't you try something geometric rather than numerical. I propose the following approach.

Let the points from the loop form the sequence $\alpha_i \,\, : \,\, i = 1, 2, 3 ... I$ and as you said, all of them lie on a given smooth surface. Furthermore, you know the unit normal ${n}$ everywhere on the surface. Then you know the unit normal to the surface at every point $\alpha_i$ from the loop. Denote that normal by $n_i$.

For each $i = 2,3,... I-1$:

1. First calculate the vector $$b_i = (\alpha_{i+1} - \alpha_i)\times(\alpha_i - \alpha_{i-1})$$.

2. Then calculate the unit tangent to the loop at the point $\alpha_i$ as follows: $$T_i = \frac{n_i \times b_i}{{\|}n_i \times b_i\|}$$ The idea is that the three points $\alpha_{i+1}, \, \alpha_i, \, \alpha_{i-1}$ form an approximation of the osculating plane of the loop at $\alpha_i$ and the tangent $T_i$ must lie simultaneously in the osculating plane and in the plane tangent to the surface at $\alpha_i$. The vectors $\alpha_{i+1} - \alpha_i$ and $\alpha_{i} - \alpha_{i-1}$ span the approximate osculating plane so their cross product $b_i$ is normal to the osculating plane. Furthermore, $n_i$ is normal to the tangent plane to the surface at $\alpha_i$. Therefore, the tangent must be normal to both vectors $b_i$ and $n_i$ in order to be aligned with the intersection line of the two planes. Hence the proposed double cross product between $b_i$ and $n_i$. The denominator normalizes the vector to make it of unit length.

3. Next calculate the unit binormal vector to the loop at $\alpha_i$ as follows $$B_i = \frac{b_i}{\|b_i\|}$$

4. The unit normal then can be computed as $$N_i = B_i \times T_i$$

5. To calculate the (total) curvature of the loop, compute the scalar $$\kappa_i = \frac{1}{2}\left( \frac{\|T_{i+1} - T_i\|}{\|\alpha_{i+1} - \alpha_i\|} + \frac{\|T_{i} - T_{i-1}\|}{\|\alpha_{i} - \alpha_{i-1}\|}\right)$$
The idea is that we approximate $\Big\|\frac{dT}{ds}\Big\|$, where $s$ is the arc-length parametrization of the loop $\alpha$, by calculating the magnitute of the change of the unit tangent from $T_i$ to $T_{i+1}$ and then dividing by the distance traveled $\|\alpha_{i+1} - \alpha_i\|$ from point $\alpha_i$ to point $\alpha_{i+1}$. The norm $\|\alpha_{i+1} - \alpha_i\|$ serves as an approximation of arc-length. We do the same for the loop points with indices $i-1$ to $i$ and finally we average the result.

6. Now, you can interpret the vector $\kappa_i N_i$ as an approximation of the derivative $\frac{dT}{ds} = \frac{d^2\alpha}{ds^2}$ where $s$ is the arc-length parametrization of the loop $\alpha$

7. The geodesic curvature of the loop at point $\alpha_i$ can be calculated as $$\kappa_g(i) = \kappa_i \Big(N_i \cdot (n_i \times T_i)\Big)$$

  • $\begingroup$ Thank you very much for your elaborate response. I will implement this solution and get back to you on the results tomorrow. I don't have enough reputation to upvote but I really appreciate your time and effort. $\endgroup$
    – Amit Singh
    Commented Apr 23, 2019 at 15:50
  • $\begingroup$ Thank you. I implemented this for a test case (updated in the question). It works well! $\endgroup$
    – Amit Singh
    Commented Apr 26, 2019 at 19:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.