How to calculate the geodesic curvature of a discrete 3D curve?

Question

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

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

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.

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