Hi,
is there a way to find the normals to a (non-linear) spline at calculated
points?
TIA
Steven

Post a reply to this message

"stevenvh" <nomail@nomail> wrote:
> Hi,
> is there a way to find the normals to a (non-linear) spline at calculated
> points?
> TIA
> Steven

I would recommend finite differences.  The fist derivative is the tangent.  So
if r(t) is the vector defined by the spline, the first derivative is

dr/dt = ( r(t+eps) - r(t-eps) ) / (2*eps)

t = (dr/dt) / |dr/dt|

where epsilon is some small constant.  So normalize the result and you have a
tangent vector.  At this point all you really need to do is cross the tangent
vector with a random vector and you should have a normal vector.  However, if
you want the normal vector in the direction of the curvature, you need to
calculate the second derivative with a finite difference.  As long as the
spline isn't a straight line this should not be collinear with the tangent:

d2r/dt2 = ( r(t+eps) - 2.0*r(t) + r(t-eps) ) / (eps^2)

The problem is that this won't be perpendicular to the tangent vector.  So the
solution is to define the normal vector by subtracting the component of the
second derivative that lies along the tangent vector.  Then all you have is a
vector that's perpendicular to the tangent vector.  Thus

n = d2r/dt2 - t * vdot( d2r/dt2, t )

n = n / |n|

Let me know if that's not clear enough.  I hope that's right, but someone asked
me this before so I even have an OpenGL demo with NURBS if you're interested.

 - Ricky

Post a reply to this message

Hi Ricky,
thanks for the explanation.

I want to place a number of cuboids along a spline path, and rotate them
according to the twisting of the path.
I'll need the tangent as one of the cuboid's axes, but I also need the direction
of curvature to know how the cuboid is to be rotated around the tangent axis.
That's the background.

I've been thinking :-)
Finding the tangent axis is easy: that's r(t+1) - r(t-1), translated to r(t),
right?
Now, as long as r(t-1), r(t) and r(t+1) are not collinear, they define a circle
whose center (c) is the center of curvature in r(t), right? So my normal vector
is r(t) - c(t).
Two points of attention:
1. r(t-1), r(t) and r(t+1) are collinear. In that case I would move c with r,
i.e. c(t) = c(t-1) + r(t) - r(t-1).
2. if the three points are nearly collinear the problem may be ill-conditioned.
I guess I could treat them as collinear if the radius of curvature becomes very
large.

So, what do you think?

As for your OpenGL demo, I have no experience, but I'm interested. Is there a
way on this forum to send you my e-mail address privately?

Again thanks.
Steven

Post a reply to this message

> I've been thinking :-)
> Finding the tangent axis is easy: that's r(t+1) - r(t-1), translated to 
> r(t),
> right?
> Now, as long as r(t-1), r(t) and r(t+1) are not collinear, they define a 
> circle
> whose center (c) is the center of curvature in r(t), right? So my normal 
> vector
> is r(t) - c(t).
> Two points of attention:
> 1. r(t-1), r(t) and r(t+1) are collinear. In that case I would move c with 
> r,
> i.e. c(t) = c(t-1) + r(t) - r(t-1).
> 2. if the three points are nearly collinear the problem may be 
> ill-conditioned.
> I guess I could treat them as collinear if the radius of curvature becomes 
> very
> large.

Given three points you can define a plane, so if you take three
close together points from a path you can find the tangent.
Once you have the tangent, you can take a tangent point, and
two points from the path to define a plane that has a normal
in the direction of curve.

This should work in most situations, but may have odd behavior
along a straight line. You may also need to check the direction
of the normal and reverse it if the curve is flipped, though this
sample doesn't seem to misbehave without that sort of check.

#include "colors.inc"
#include "math.inc"
#include "transforms.inc"

#default {finish{ambient 0.7}}
// --- camera, lights & 
background ------------------------------------------
camera {
   location  <0.0, 1, -4.0>
   direction 1.5*z
   right     x*image_width/image_height
   look_at   <0.0, 0.5,  0.0>
}

background {Gray85}

light_source {
  <-30, 30, -30>
  color rgb <1, 1, 1>
}

// ---  
macros ---------------------------------------------------------------
#declare Path_Spline = spline {
   natural_spline
   0, <-2,0,0>
   1, <-1,0,0>
   2, <-0.5,0.5,0>
   3, <0,1,1>
   4, <0.5,0.5,0>
   5, <1,0,0>
   6, <2,0,0>
};

#macro Path_Normal(N, rez)
   #local P1 = Path_Spline(N-rez);
   #local P2 = Path_Spline(N);
   #local P3 = Path_Spline(N+rez);
   #local HPerp = VPerp_To_Plane(vnormalize(P2-P1),vnormalize(P3-P2));
   #local Perp = VPerp_To_Plane(vnormalize(HPerp-P2),vnormalize(P3-P2));
   (Perp)
#end

// show spline
#local C = 0;
#while (C <= 6)
   sphere{Path_Spline(C),0.01 pigment {Yellow}}
   #local C=C+0.005;
#end

// inside of curve
#local C = 0.5;
#while (C <= 5.5)
   #local Perp = Path_Normal(C,0.1);
   sphere {Path_Spline(C)+Perp/10,0.01 pigment{Orange}}
   #local C=C+0.005;
#end

Post a reply to this message

"Tim Attwood" <tim### [at] comcastnet> wrote:
> Given three points you can define a plane, so if you take three
> close together points from a path you can find the tangent.
> Once you have the tangent, you can take a tangent point, and
> two points from the path to define a plane that has a normal
> in the direction of curve.
>
> This should work in most situations, but may have odd behavior
> along a straight line. You may also need to check the direction
> of the normal and reverse it if the curve is flipped, though this
> sample doesn't seem to misbehave without that sort of check.
>
> #include "colors.inc"
> #include "math.inc"
> #include "transforms.inc"
>
> #default {finish{ambient 0.7}}
> // --- camera, lights &
> background ------------------------------------------
> camera {
>    location  <0.0, 1, -4.0>
>    direction 1.5*z
>    right     x*image_width/image_height
>    look_at   <0.0, 0.5,  0.0>
> }
>
> background {Gray85}
>
> light_source {
>   <-30, 30, -30>
>   color rgb <1, 1, 1>
> }
>
> // ---
> macros ---------------------------------------------------------------
> #declare Path_Spline = spline {
>    natural_spline
>    0, <-2,0,0>
>    1, <-1,0,0>
>    2, <-0.5,0.5,0>
>    3, <0,1,1>
>    4, <0.5,0.5,0>
>    5, <1,0,0>
>    6, <2,0,0>
> };
>
> #macro Path_Normal(N, rez)
>    #local P1 = Path_Spline(N-rez);
>    #local P2 = Path_Spline(N);
>    #local P3 = Path_Spline(N+rez);
>    #local HPerp = VPerp_To_Plane(vnormalize(P2-P1),vnormalize(P3-P2));
>    #local Perp = VPerp_To_Plane(vnormalize(HPerp-P2),vnormalize(P3-P2));
>    (Perp)
> #end
>
> // show spline
> #local C = 0;
> #while (C <= 6)
>    sphere{Path_Spline(C),0.01 pigment {Yellow}}
>    #local C=C+0.005;
> #end
>
> // inside of curve
> #local C = 0.5;
> #while (C <= 5.5)
>    #local Perp = Path_Normal(C,0.1);
>    sphere {Path_Spline(C)+Perp/10,0.01 pigment{Orange}}
>    #local C=C+0.005;
> #end

Thanks, Tim.
While I agree with your code, I'm a bit puzzled by the results.
I changed point 3 of the spline to <0,1,0>, so that the spline lies in the Z=0
plane.
First observation: if the spline lies in the Z=0 plane, so should the orange
locus, but it doesn't. I don't know why.
OTOH the orange curve is discontinuous where you expect it, i.e. at the
inflection points of the spline.

I also would write (P3-P1) instead of (P3-P2) in the calculation of Perp, but
since P1, P2 and P3 are nearly collinear this shouldn't make much difference
(and it doesn't).

Thanks again.
Steven

(Learned this lesson: "VPerp_To_Plane" :-))

Post a reply to this message

> Thanks, Tim.
> While I agree with your code, I'm a bit puzzled by the results.
> I changed point 3 of the spline to <0,1,0>, so that the spline lies in the 
> Z=0
> plane.
> First observation: if the spline lies in the Z=0 plane, so should the 
> orange
> locus, but it doesn't. I don't know why.

Oops, the new point is P2-HPerp, not HPerp...

> OTOH the orange curve is discontinuous where you expect it, i.e. at the
> inflection points of the spline.
>
> I also would write (P3-P1) instead of (P3-P2) in the calculation of Perp, 
> but
> since P1, P2 and P3 are nearly collinear this shouldn't make much 
> difference
> (and it doesn't).

Depending on which order you define the points for HPerp, the
Path_Normal will be on the inside of the curve, or on the outside.

#macro Path_Normal(N, rez)
   #local P1 = Path_Spline(N-rez);
   #local P2 = Path_Spline(N);
   #local P3 = Path_Spline(N+rez);
   #local HPerp = VPerp_To_Plane(P1-P2,P3-P2);
   //#local HPerp = VPerp_To_Plane(P2-P1,P3-P1);
   #local Perp = VPerp_To_Plane((P2-HPerp)-P3,P2-P3);
   (Perp)
#end

Post a reply to this message

stevenvh nous apporta ses lumieres en ce 2008/04/09 16:20:
> Hi Ricky,
> thanks for the explanation.
> 
> I want to place a number of cuboids along a spline path, and rotate them
> according to the twisting of the path.
> I'll need the tangent as one of the cuboid's axes, but I also need the direction
> of curvature to know how the cuboid is to be rotated around the tangent axis.
> That's the background.
> 
> I've been thinking :-)
> Finding the tangent axis is easy: that's r(t+1) - r(t-1), translated to r(t),
> right?
> Now, as long as r(t-1), r(t) and r(t+1) are not collinear, they define a circle
> whose center (c) is the center of curvature in r(t), right? So my normal vector
> is r(t) - c(t).
> Two points of attention:
> 1. r(t-1), r(t) and r(t+1) are collinear. In that case I would move c with r,
> i.e. c(t) = c(t-1) + r(t) - r(t-1).
> 2. if the three points are nearly collinear the problem may be ill-conditioned.
> I guess I could treat them as collinear if the radius of curvature becomes very
> large.
> 
> So, what do you think?
> 
> As for your OpenGL demo, I have no experience, but I'm interested. Is there a
> way on this forum to send you my e-mail address privately?
> 
> Again thanks.
> Steven
> 
> 
> 
> 
An easy way: use the spline follow macro.
You define a spline, then call the macro with the spline name and a position as 
parameters. It returns the coordinate at that point along with the orientation. 
There is a sample scene showing it's working and application.
You'll find it in /scene/animation/splinefollow/splinefollow.pov

-- 
Alain
-------------------------------------------------
You know you've been raytracing too long when you prefer bald romatic partners, 
because they're easier to model.
John VanSickle

Post a reply to this message

stevenvh nous apporta ses lumieres en ce 2008/04/09 16:20:


> As for your OpenGL demo, I have no experience, but I'm interested. Is there a
> way on this forum to send you my e-mail address privately?
> 
> Again thanks.
> Steven
> 
> 
The common way to send an e-mail address on a news group is to obfuscate it. It 
often look like:

some.name AT domain DISCARD DOT com

You can also add some trhowaway character in the address. For example, I've 
added a single DOT in my user name. Try to send me an e-mail with my address as 
it apears, and it will bounce, remove the dot after "electro" and I will receive it.

-- 
Alain
-------------------------------------------------
Please hassle me, I thrive on stress.

Post a reply to this message

"stevenvh" <nomail@nomail> wrote:

> I've been thinking :-)
> Finding the tangent axis is easy: that's r(t+1) - r(t-1), translated to r(t),
> right?
....
> So, what do you think?

I think you pretty much have it, and by now everyone else has pretty much hit
the nail on the head.  Had a small attack of real life these last couple days,
but here's the code anyway, and here's the OpenGL demo:

http://rsreusser.googlepages.com/nurbs3d.tar.gz

Fun to play around with, but not a great example of good coding.  I'm not sure
how to compile that on Windows, although it should be perfectly compatible.  On
Linux or Mac OS X it should be trivial.

 - Ricky

And just for kicks, this one calculates and plots the radius of curvature too:

camera{ location <3.5,2,-3.5> look_at <2.7,0,0> }
plane{y,0 pigment{rgb 1}}
background{rgb 0.5}
light_source{<25,25,-25>, rgb 1
    area_light 10*x,10*y,5,5 circular orient jitter
}

#declare sp=
    spline {
 natural_spline
 0.0, <0.0, 0,  0>,
 0.2, <1.0, 0.2, 0.5>,
 0.4, <2.0, 0.4,-0.9>,
 0.6, <3.5, 0.2, 1.0>,
 0.8, <4.0, 0.1,-0.3>,
 1.0, <5.0, 0.0, 0.0>
    }

cylinder{0,x,0.01 pigment{rgb x}}
cylinder{0,y,0.01 pigment{rgb y}}
cylinder{0,z,0.01 pigment{rgb z}}
union{
    #declare npt = 100;
    #declare c=0;
    #declare p0 = sp(0.0);
    #declare rad=0.02;
    #while(c<npt)
 #declare tv = c/(npt-1);
 #declare p1 = sp(tv);
 sphere{p0,rad pigment{rgb <tv,0,1-tv>}}
 #if(c>0) cylinder{p0,p1,rad pigment{rgb <tv,0,1-tv>}} #end
 #declare p0=p1;
 #declare c=c+1;
    #end
}
union{
    #declare points = array[4]{0.17, 0.37, 0.61, 0.85}
    #declare eps = 0.0001;
    #declare ncirc = dimension_size(points,1);
    #declare c=0;
    #while(c<ncirc)
 #declare tv=points[c];
 #declare pm=sp(tv-eps);
 #declare p0=sp(tv);
 #declare pp=sp(tv+eps);
 #declare drdt = (pp-pm)/(2.0*eps);
 #declare velocity = vlength(drdt);
 #declare tangent = vnormalize(drdt);
 #declare d2rdt2 = (pp-2.0*p0+pm)/(eps*eps);
 #declare accel = vlength(d2rdt2);
 #declare perp = vnormalize(d2rdt2-vdot(tangent,d2rdt2)*tangent);
 #declare crad = velocity*velocity/accel;
 #declare perp2 = vcross(tangent,perp);
 #declare e1=perp;
 #declare e2=perp2;
 #declare e3=tangent;
 #declare c0=p0+perp*crad;
 torus{crad,rad/2
     matrix < e1.x, e1.y, e1.z,
       e2.x, e2.y, e2.z,
       e3.x, e3.y, e3.z,
       c0.x, c0.y, c0.z >
 }
 cylinder{p0,p0+e1/5,rad/4 pigment{rgb x}}
 cylinder{p0,p0+e2/5,rad/4 pigment{rgb y}}
 cylinder{p0,p0+e3/5,rad/4 pigment{rgb z}}
 #declare c=c+1;
    #end
    pigment{rgb 2}
}

Post a reply to this message