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

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.

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

• 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. Apr 23 '19 at 15:50
• Thank you. I implemented this for a test case (updated in the question). It works well! Apr 26 '19 at 19:10