Computing numerical derivatives

Question

I am trying to create a sweeping surface, for which I need the frenet frame of a curve. I am trying to compute this for arbitrary curves but for testing I am just using the parametric unit half circle. i.e. $C(t) = (\cos(t), \sin(t), 0)$

I am computing my derivatives using central differences, i.e.:

pub fn frenet_trihedron(x: f32, f: &dyn Fn(f32) -> Vec3) -> (Vec3, Vec3, Vec3)
{
    let h = 1e-6 as f32;

    let normal = (-2.0 * f(x) + f(x + h) + f(x - h)) / (h * h);
    let tangent = (f(x + h) - f(x - h)) / (2.0 * h);

    let tangent = tangent.normalize();
    let normal = normal.normalize();
    let binormal = tangent.cross(&normal);

    println!("t {}, n {}, b {}", tangent, normal, binormal);

    (tangent, normal, binormal)
}

This however seems to yield very impossible results, for example:

t 
  ┌            ┐
  │ -0.8274758 │
  │ 0.56150144 │
  │          0 │
  └            ┘

, n 
  ┌   ┐
  │ 0 │
  │ 1 │
  │ 0 │
  └   ┘

, b 
  ┌            ┐
  │          0 │
  │          0 │
  │ -0.8274758 │
  └            ┘

That normal is not orthogonal to the tangent not even close What am I doing wrong?

This is the calling code:

let f = |v| Vec3::new(f32::cos(v * PI), f32::sin(v * PI), 0.0);
let (t, n, b) = frenet_trihedron(v, &|x| f(x));

It seems the error is being cause by floating point imprecision. Using 1e-3 asmy epsilon yields remarkably better results. so maybe f32 is just not precise enough.

Answer 1

The unit tangent vector is $\hat t(x)=\frac{C'(x)}{\|C'(x)\|}$. The normal direction is the derivative of this unit tangent, not just the second derivative $$ n(x)=\hat t'(x)=\frac{(\|C'(x)\|^2I-C'(x)C'(x)^T)C''(x)}{\|C'(x)\|^3} =\frac{(I-\hat t(x)\hat t(x)^T)C''(x)}{\|C'(x)\|} $$ then again the unit normal is the normalization after this projection $\hat n(x)=\frac{n(x)}{\|n(x)\|}$. The scalar denominator cancels in the normalization, so can be left out from the start

The winding direction after that is then indeed the cross product.

In Python this can look like

from numpy import array, cross, cos, sin, pi
from numpy.linalg import norm
def C(x): return np.array((cos(pi*x),sin(pi*x),0*x))

def dreibein(C,x):
    h = 0.005*max(1,abs(x))
    Cm,C0,Cp = C(x-h),C(x),C(x+h)
    #tangent
    t = Cp-Cm; t/=norm(t)
    # normal
    n = Cp-2*C0+Cm; n -= n.dot(t)*t; n /= norm(n)
    # winding
    b = cross(t,n)
    return t,n,b

t,n,b =dreibein(C,0.31019)
print(t,n,b)

returns

[-0.82741594  0.56158959  0.        ] [-0.56158959 -0.82741594  0.        ] [0. 0. 1.]

which visibly is an orthonormal basis.

The same basis can be found from the columns of $Q$ in a suitable QR decomposition

Q,R = np.linalg.qr(np.transpose([Cp-Cm,Cp-2*C0+Cm]),"complete")
Q = Q @ np.diag([*np.sign(np.diag(R)),1])
t,n,b = Q.T