orbingol / NURBS-Python

Object-oriented pure Python B-Spline and NURBS library
https://onurraufbingol.com/NURBS-Python/
MIT License
640 stars 156 forks source link

Normal vector NURBS curve wrong #107

Closed andreabonacini closed 4 years ago

andreabonacini commented 4 years ago

Hi, I'm using the curve.normal() method to evaluate the frenet frame on the curve s at a certain parameter u but it seems to be wrong. The normal vector of a curve is always in the direction of the curvature. Plotting the normal vector in output to curve.normal() it is correct only in one verse of the parametrization, not both.

The normal is always defined as the normalized second derivative and then used with the tangent vector to compute the binormal vector in order to obtain a right-hand coordinate frame:

t = (ds/du) / ||ds/du|| n = (d^2 s/du^2) / ||(d^2 s/du^2)|| = (dt/du) / ||(dt/du)|| = (dt/du) / k
b = t x n

where k = ||(dt/du)|| is the curvature of the curve at a certain u value.

The curve.tangent() method works fine and computing the tangent vector manually with those formulas it results the same. The curve.normal() vector instead is wrong. Looking at the code you compute it in this way:

t = ders[1] tn = [0.0, 0.0, 1.0] ---> what this should represent? n = linalg.vector_cross(t, tn) ---> n its not defined in this way

As a quick reference: https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas

I look forward to your feedback Andrea

orbingol commented 4 years ago

Well, I usually prefer not to discuss math on Github and I have been avoiding it for a while. A textual system is not good for such discussions. You We need a blackboard. There are also some discussions on the issue tracker about the curve normal.

Frenet-Serret says that the derivative of the tangent is the normal, which literally corresponds to the curvature vector for the NURBS curves. The definition of being normal to a vector means there would be 90 degrees between the vectors. The curvature vector doesn't have to have 90 degrees with the tangent vector. Another thing is that, when n-dimensional curve is considered having 90 degrees between 2 vectors becomes very indefinite where n > 2. Therefore, I replaced the curve normal with "normal to the tangent vector" idea -- still incorrect when n-dimensions are considered, where n > 3. It is already in the documentation.

I think implementation of a method like Curve.normal or Curve.tangent is incorrect, and they are removed from 6.x. If you need a curvature vector, simply take the second derivative. For the tangent, use first derivative. However, you should know that it won't be normal to the curve, as it might not satisfy the condition of tangent is perpendicular to the normal.

When it comes to surfaces, the normal makes perfect sense. For curves and volumes, it does not.

orbingol commented 4 years ago

I am labeling this issue as "duplicate", as several discussions exist on the issue tracker. I am also labeling as "won't fix" as I have removed these functions from 6.x. This issue is still open to the discussions, please don't worry.

All comments and suggestions are welcome!

andreabonacini commented 4 years ago

I used the math just to support the question, since it's easy math. I read the other questions about normal but it isn't fixed yet, that's why i asked.

Anyway, the curvature k is not a vector, it's a scalar. So i think that it's incorrect to say if curvature has or not 90 degrees with tangent (that is a vector). I suggest this link for reference. https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Vector_Calculus/2%3A_Vector-Valued_Functions_and_Motion_in_Space/2.3%3A_Curvature_and_Normal_Vectors_of_a_Curve

I think implementation of a method like Curve.normal or Curve.tangent is incorrect, and they are removed from 6.x. If you need a curvature vector, simply take the second derivative. For the tangent, use first derivative. However, you should know that it won't be normal to the curve, as it might not satisfy the condition of tangent is perpendicular to the normal.

I don't understand this, beacuse in other languages such as C# i used a library where these functions were implemented. The normal vector is orthogonal to the tangent vector with direction concordant to the curvature of the curve.

Also, as i said before i'm already manually computing tangent and normal vector (binormal is a consequence). Taking the first derivative produce a correct tangent vector, but it's not true for the normal vector. Maybe due to some sort of approximation in derivative computation?

orbingol commented 4 years ago

The normal vector is orthogonal to the tangent vector with direction concordant to the curvature of the curve.

Would you elaborate a little on how a scalar value has a direction?

Edit: I am guessing this question is not that clear. Obviously, the tangent vector is a vector, not a scalar. The curvature is a scalar value. Here, comparing a vector with a scalar value doesn't make sense to me.

Also, as i said before i'm already manually computing tangent and normal vector (binormal is a consequence). Taking the first derivative produce a correct tangent vector, but it's not true for the normal vector. Maybe due to some sort of approximation in derivative computation?

There is another thing that needs to be considered here. The reference link assumes that the curve is arc-length parameterized, if I understand correctly. Would this change the results?

Edit 2: For me Frenet-Serret is the correct solution (I am pretty sure I answered a question some time ago same way you explained in the original post with the exact link you posted as a reference, but don't remember when or where), but I was unable to numerically validate the tangent vector being perpendicular to the normal vector that I calculate using the Frenet-Serret approach. It is obvious that I am making a mistake in the calculations and I don't know where (here it is easy to blame the derivative function, but I don't think it is the problem). Till I find the correct solution using Frenet-Serret, I have to stick to the current solution -- which I am sure not 100% correct.

orbingol commented 4 years ago

I have a feeling that I might need to answer my questions myself. Let me experiment a bit and try to make Frenet-Serret working. I am closing this issue. Feel free to reopen it, if you like.

andreabonacini commented 4 years ago

Would you elaborate a little on how a scalar value has a direction?

Edit: I am guessing this question is not that clear. Obviously, the tangent vector is a vector, not a scalar. The curvature is a scalar value. Here, comparing a vector with a scalar value doesn't make sense to me.

When i said "with direction concordant to the curvature of the curve" i'm not strictly talking about the curvature parameter k ( a scalar). Think about a circle, the normal vector must point toward the center of the circle, not to outside. And this happen indipendently to how i travel the curve. In your library, this does not happen as a consequence of how you implemented it. Try to create a circle with the well known control points/node vector. If you invert the order of the control point you could see that your normal function will be different ( since its wrong how you computed it), but actually the normal in a point on the curve MUST be always the same, since its depend just on the second derivative.

es) A circle with radius R = 1 has a curvature k = 1, since curvature is the inverse of the radius of the osculating circle (k = 1/R), that for a circle is the circle itself.

There is another thing that needs to be considered here. The reference link assumes that the curve is arc-length parameterized, if I understand correctly. Would this change the results?

As you can see in 2.3.16 and 2.3.19 the tangent vector does not dipend to the parametrization (and that is what we should expect since the geometry of the curve does not change, what is changing is only the parametrization). This is obviusly true also for normal vector.

Anyway, i was trying to find the problem that could be in 2 points of the algorithm:

NURBS library

tangent = np.array(curve.tangent(u)) ---> seems right normal = np.array(curve.normal(u)) ---> wrong binormal = np.array(curve.binormal(u)) ---> wrong

Manually computed

der = np.array(curve.derivatives(u,2)) k= np.linalg.norm(der[2]) #curvature parameter real_tang = der[1]/np.linalg.norm(der[1]) ---> right approach, same result of library real_norm = der[2]/k---> right approach, "inside" the curvature, but not perfectly orthogonal to tangent real_binorm = np.cross(real_tang,real_norm) ---> right approach, right output except for deriving problem of real_norm

Also i made a test creating a NURBS-circle and comparing it to an analytic circle, where in the last one i computed points and their derivatives analytically (without your derivative class method). Both have radius = 1 so curvature must be always equal to 1. Here the results:

#NURBS CIRLCE degree = 2
control_points = [[1, 0, 0], [1, 1, 0], [0, 1, 0],[-1, 1, 0],[-1, 0, 0],[-1, -1, 0],[0, -1, 0],[1, -1, 0],[1, 0, 0]] weights = [1,1/math.sqrt(2),1,1/math.sqrt(2),1,1/math.sqrt(2),1,1/math.sqrt(2),1] knot_vector = [0, 0, 0, 1/4, 1/4, 1/2, 1/2, 3/4, 3/4, 1, 1, 1] curve = Nurbs(degree,control_points,weights,knot_vector) #Nurbs is a child class i created on my own camp = 16 u_step = 1/camp u = 0 list_k = [] for i in range(camp): derivate = np.array(curve.derivatives(u,2)) k = np.linalg.norm(derivate[2]) list_k.append(k) u = u + u_step print("Curvature:") print(list_k)

Output: Curvature: [34.63655040935661, 41.240965995073175, 43.922656064975335, 41.240965995073175, 34.63655040935661, 41.240965995073175, 43.922656064975335, 41.240965995073175, 34.63655040935661, 41.240965995073175, 43.922656064975335, 41.240965995073175, 34.63655040935661, 41.240965995073175, 43.922656064975335, 41.240965995073175] ---> values are wrong

#ANALYTIC CIRCLE camp = 100 uk = 0 uk_step = (2*math.pi)/camp circ_pts = [] list_tk = [] list_Kk = [] list_nk = [] list_bk [] for i in range(camp):

position

xk = math.cos(uk) yk = math.sin(uk) pk = np.array([xk,yk,0]) circ_pts.append(pk)

first derivative

dxk = -math.sin(uk) dyk = math.cos(uk) dpk = np.array([dxk,dyk,0]) tk = dpk/np.linalg.norm(dpk) ---> tangent vector in uk list_tk.append(tk)

second derivative

d2xk = -math.cos(uk) d2yk = -math.sin(uk) d2pk = np.array([d2xk,d2yk,0]) Kk = np.linalg.norm(d2pk) ---> curvature value in uk nk = d2pk/Kk ---> normal vector in uk list_Kk.append(Kk) list_nk.append(nk) bk = np.cross(tk,nk) ---> binormal vector in uk list_bk.append(bk) uk = uk + uk_step

val = ... print("Tangent vector (analytic circle):") print(list_tk[val]) print("Normal vector (analytic circle):") print(list_nk[val]) print("Binormal vector (analytic circle):") print(list_bk[val]) print("Curvature (analytic circle):") print(list_Kk[val])

Output:

val = 10 Tangent vector (analytic circle):
[-0.58778525 0.80901699 0. ] Normal vector (analytic circle):
[-0.80901699 -0.58778525 0. ] Binormal vector (analytic circle):
[0. 0. 1.] Curvature (analytic circle): 0.9999999999999999

val = 20 Tangent vector (analytic circle):
[-0.95105652 0.30901699 0. ] Normal vector (analytic circle):
[-0.30901699 -0.95105652 0. ] Binormal vector (analytic circle):
[0. 0. 1.] Curvature (analytic circle): 1.0

val = 36 Tangent vector (analytic circle):
[-0.77051324 -0.63742399 0. ] Normal vector (analytic circle):
[ 0.63742399 -0.77051324 0. ] Binormal vector (analytic circle):
[0. 0. 1.] Curvature (analytic circle): 1.0

Conclusion: The visualization of the two circle with matplotlib looks like the same ( this should means B-function are computed right?). The curvature of the of the NURBS-circle is wrong since values are != 1 and also they are changing through the curve. The curvature of the ANALYTIC-circle is right, since its always equalt to 1. Also plotting the axis with ax.quiver t,n,b are always ortoghonal and well placed (as i would expect).

I think i gave you a lot of tools and proofs to try to find the issue(s) and maybe it's better to leave the question open. If there are any other doubt let me know it. Andrea

orbingol commented 4 years ago

@andreabonacini thanks! I'll take a look when I find time.

portnov commented 4 years ago

Eeh. Normal vector is always perpendicular to the tangent, that's why it's called normal :) Tangent vector is the derivative of the curve function, normal is the second derivative, and binormal is their cross product. In this formulation, lengths of these vectors depend on parametrization of the curve, but not directions. So usually they are just normalized. One thing about tangent/normal/binormal that confuses people is that direction of normal flips when you pass a point where curve has zero curvature. Binormal flips accordingly.

orbingol commented 4 years ago

Normal vector is always perpendicular to the tangent, that's why it's called normal :)

That is what I am trying to explain. However @andreabonacini says

real_norm = der[2]/k---> right approach, "inside" the curvature, but not perfectly orthogonal to tangent

andreabonacini commented 4 years ago

Normal vector is always perpendicular to the tangent, that's why it's called normal

That is what I am trying to explain.

You both are missing the point that exist infinite vectors that are normal to the "Tangent" versor, but just one that is also the second derivative of the curve at a certain parameter value (lying on the osculating plane) and that is called "Normal" versor.

So i dont care about a method that compute one of the normal vector (orthogonal to the tangent versor) that correspond to the "Normal" versor just if the curve lies on the XY-plane and half of the time depending on how you travel along the curve.

@portnov they are versor, they are always normalized by definition. Also the direction of normal does not change when you pass a curvature with null curvature but when you pass an inflection point.

I'm showing you the visualization of the circle example:

Cattura Cattura