1

I have several coordinates for curves in 3D space and can calculate the curvature and torsion via splipy:

from splipy import curve_factory
points = np.load('my_coordinates.npy')
curve = curve_factory.cubic_curve(points)
sample_points = 1000 
t = numpy.linspace(curve.start()[0], curve.end()[0],sample_points)
curvature = curve.curvature(t)
torsion = curve.torsion(t)

I tried to find functions that now reconstruct my original curve via the curvature and torsion. The only thing I found was this: reconstruct curve mathematically where it is explained that to reconstruct the curve you have to solve a system of ordinary differential equations.

I tried to solve the ode system that is given in the above link with the method that is shown here: how to solve ode system However, my problem is that my curvature and torsion are in discrete form (as in contrast to the given example). Here my attempt (Note, that I save the x[0] variable for gamma):

At first I defined the right-hand-side with the equations from the link:

def rhs(s, x):
    k = lambda s : kappa(s)  #curvature
    t = lambda s : tao(s)    #torsion
    return [x[1], k(s)*x[2], -k(s)*x[1]+t(s)*x[3], -t(s)*x[2]]

I then interpolated the curvature and torsion with scipy to make them work as functions so that they can deal with continuous input. Furthermore, I set the initial conditions. Because are values are vectors, I set T,N and B as the unit vectors and gamma as the origin of the coordinate system:

kappa = np.load("curvature.npy")
tao = np.load("torsion.npy")
length = len(kappa)
x = np.linspace(0,1,length)
kappa = interpolate.interp1d(x, kappa)
tao = interpolate.interp1d(x, tao)
initial_cond = np.array([[0,0,0],[0,1,0],[1,0,0],[0,0,1]])

and then tried to solve the whole thing with scipy.integrate:

res = solve_ivp(rhs, (0,1), initial_cond)

And this is the error I get: Error

Would be even possible to solve vector-odes with solve_ivp or do I have to deconstruct it into 4*3=12 different equations?

Thanks for the help!

MangoMat
  • 21
  • 4
  • *"... therefore my code is not running somehow."* What does that mean? Do you get an error? If so, include the complete traceback (i.e. the complete error message) in the question. Does the program given unexpected output? If so, show the output, and explain what you expected. – Warren Weckesser Oct 24 '22 at 13:25
  • The initial vectors T and N should not be all 0. The vectors T, N, and B form an orthogonal set, the *Frenet frame*, that evolves along the curve. So the initial values T and N must be unit vectors that are perpendicular; T is the initial tangent vector, and N is the initial normal vector. See [Frenet–Serret formulas](https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas) for more details. – Warren Weckesser Oct 24 '22 at 13:31
  • thanks for the comment! I edited the question accordingly – MangoMat Oct 25 '22 at 02:39
  • I have an example of this in my "experiments" repository on github: https://github.com/WarrenWeckesser/experiments/tree/master/python/scipy/frenet. In that example, the state of the system consists of the the point `x` and the tangent and normal vectors, flattened into a single 9-d state vector. The binormal is just the cross product T ⨯ N, so it is not include in the state vector. In my code, it is assumed that the curvature and torsion are given functions. In your case, those would be the interpolators of your data. – Warren Weckesser Oct 25 '22 at 04:16
  • Thanks! I tried your git-code and it worked. However, the generated curve is not the same as the original curve. I think its because my curvature and torsion are very "spiky" and non-regular, like a signal, maybe this causes some problems or information loss. – MangoMat Oct 25 '22 at 08:54

0 Answers0