f

#### Muscle Wrapping over an cylinder - Shortest path

```Dear all,

I have two points, between which I want to calculate the shortest path along the surface of an ellipsoid. In order to obtain the shortest line between P1 and P2 around a cylinder, tangent points R1 and R2 must be calculated. The cylinder is defined by its central axis with position vector O and unit direction vector d, and by its radius r. The shortest line, including P1, P2, R1 and R2 is assumed to be in one plane.

There are several articles published on muscle wrapping, I use the following equations:
* Normal vector n: n = (P1-P2) x d x (P1-P2)   (Eq 1)
* Vector Nr1, perpendicular to the central axis and through R, can be calculated by the following three equations:
(1)    dx*Nr1x +  dy*Nr1y +  dz*Nr1z  =  0             (Eq 2)
(2)    P1x*Nr1x +  P1y*Nr1y +  P1z*Nr1z + 1 =  0   (Eq 3)
(3)    Ox*Nr1x +  Oy*Nr1y +  Oz*Nr1z + r +  1 =  0 (Eq 4)
* Tangent Point R1 = P1 + Lambda1* L1                   (Eq 5)
* Direction vector L1 = n x Nr1                                 (Eq 6)
* R1 is on the surface of the cylinder, which is described by:
[(Rx-Ox)dy  -  (Ry-Oy)dz]^2  +  [(Ry-Oy)dz  -  (Rz-Oz)dy]^2  + [(Rz-Oz)dx  -  (Rx-Ox)dz]^2  - r^2 =0                                                 (Eq 7)

lambda 1 is calculated by substituting Eq5 in Eq 7, and so R1 can be calculated. R2 can be calculated similarly.

-----------------------------------------------------------

Now, my code, when I plot R1 and R2, this is not correct. I hope someone can help me detect my mistake and maybe can tell me how to correctly plot the geodesic between R1 and R2.

----------------------------------------------------------

X1=[1 1 1];         % Startpoint cylinder
X2=[41 5 10];       % End point cylinder
P1=[30,-8,5];P2=[17,10,0];
scatter3(P1(1,1),P1(1,2),P1(1,3))
hold on
scatter3(P2(1,1),P2(1,2),P2(1,3))

d=(X2-X1); %
d=d(:).'/norm(d) % directional vector
O=X1;               % position vector

N1 = [0;0;0];    % Werkt voor aMCL, pop Ten,pmcl
options = optimset('Display','iter','Algorithm',{'levenberg-marquardt',.005},'MaxIter',600,'MaxFunEvals',1400);% Werkt voor aMCL/pmcl,Werkt voor pop Ten
[N]=fsolve(@SolveCylinder, N1,options,P1,d,r,O);

n=cross(cross((P1-P2),d),(P1-P2));   % Normal vector n
n=n(:).'/norm(n);
L1 = cross(n,N);
lambda=[0;0;0];                      % Startvalues lambda
options = optimset('Display','iter','Algorithm',{'levenberg-marquardt',.005},'MaxIter',600,'MaxFunEvals',1400);% Werkt voor aMCL/pmcl,Werkt voor pop Ten
[LL]=fsolve(@SolveCylinderR1, lambda,options,P1,L1,d,O,r);
LL

R1x=P1(1,1)+LL(1)*L1(1,1)
R1y=P1(1,2)+LL(2)*L1(1,2)
R1z=P1(1,3)+LL(3)*L1(1,3)
scatter3(R1x,R1y,R1z,'*','r')
plot3([P1(1,1) R1x],[P1(1,2) R1y],[P1(1,3) R1z],'r' )
plot3([P1(1,1) P2(1,1)],[P1(1,2) P2(1,2)],[P1(1,3) P2(1,3)],'r' )

options = optimset('Display','iter','MaxIter',600);
[N]=fsolve(@SolveCylinder, N1,options,P2,d,r,O);
L1 = cross(n,N);
N
options = optimset('Display','iter','Algorithm',{'levenberg-marquardt',.005},'MaxIter',600,'MaxFunEvals',1400);% Werkt voor aMCL/pmcl,Werkt voor pop Ten
[LL]=fsolve(@SolveCylinderR1, lambda,options,P2,L1,d,O,r);
R2x=P2(1,1)+LL(1)*L1(1,1)
R2y=P2(1,2)+LL(2)*L1(1,2)
R2z=P2(1,3)+LL(3)*L1(1,3)
scatter3(R2x,R2y,R2z,'*','r')
plot3([P2(1,1) R2x],[P2(1,2) R2y],[P2(1,3) R2z],'r' )

R1=[R1x,R1y,R1z];
R2=[R2x,R2y,R2z];
```
 0
Els
6/3/2010 1:06:04 PM
comp.soft-sys.matlab 211266 articles. 19 followers. lunamoonmoon (258) is leader.

10 Replies
848 Views

Similar Articles

[PageSpeed] 7

```I forgot to add the two equation solver functions. Still hope someone can help me.
A lot of thanks in advance.

Best wishes,

Els

_________________________
function [F] = SolveCylinderR1(lambda,P1,L1,d,O,r)

F = (((((P1(1,1)+lambda(1)*L1(1,1))-O(1,1))*d(1,2))-(((P1(1,2)+lambda(2)*L1(1,2))-O(1,2))*d(1,1)))^2   +   ((((P1(1,2)+lambda(2)*L1(1,2))-O(1,2))*d(1,3))-(((P1(1,3)+lambda(3)*L1(1,3))-O(1,3))*d(1,2)))^2   + ((((P1(1,3)+lambda(3)*L1(1,3))-O(1,3))*d(1,1))-(((P1(1,1)+lambda(1)*L1(1,1))-O(1,1))*d(1,3)))^2 -r^2 );
_____________________________
function [F] = SolveCylinder(N1,P1,d,r,O)

F = [ ((P1(1,1)*N1(1))+ (P1(1,2)*N1(2)) +(P1(1,3)*N1(3))+1);
(d(1,1)*N1(1)+ d(1,2)*N1(2) +d(1,3)*N1(3));
((O(1,1)*N1(1))+ (O(1,2)*N1(2)) +(O(1,3)*N1(3))+r+1)];
```
 0
Els
6/4/2010 9:00:30 AM
```Sorry to post here again, but the deadline of my project is approaching, and I still do not see what I am doing wrong. Tried all the - and + thing, but haven't found it.

Sorry, I am a bit desperite.
```
 0
Els
6/4/2010 3:27:04 PM
```"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hu89bs\$ito\$1@fred.mathworks.com>...
> I have two points, between which I want to calculate the shortest path along the surface of an ellipsoid. In order to obtain the shortest line between P1 and P2 around a cylinder, tangent points R1 and R2 must be calculated. The cylinder is defined by its central axis with position vector O and unit direction vector d, and by its radius r. The shortest line, including P1, P2, R1 and R2 is assumed to be in one plane.
> ........
- - - - - - - - - -
I don't understand your reasoning involving what appear to be two discretely separated points P1 and P2 on a cylinder.  Also I don't know what you mean by the two "tangent points" R1 and R2, unless perhaps you mean the points at the ends of unit-length tangent vectors.

However, I do know that the geodesic on a strictly circular cylindrical surface is always a helical path spiraling around the cylinder, and that is not a curve that is contained in any plane.  In particular the statement "the shortest line, including P1, P2, R1 and R2 is assumed to be in one plane" is incorrect.

Demonstrating this fact involves reasoning on the infinitesimal level and is actually a problem in the calculus of variations.  A basic principle for any geodesic on a surface is that the curve's normal - that is the direction of change of its unit-length tangent vector - must always be parallel to the cylinder's surface normal.  In this case that forces the angle of the cylindrical coordinates to be a linear function of the axial coordinate, and that necessarily defines a helix.  It is rather easy to show if the above principle is properly applied.

Roger Stafford
```
 0
Roger
6/5/2010 7:12:06 AM
```Dear Roger,

Thanks for your complete response, very clear.

P1 and P2 are no points on a cylinder, they are random points in space. The R1 and R2 are tangent points on the surface of the cylinder. So there will come a straight line from P1 to R1, then a  geodesic over the surface of the cylinder from R1 to R2 and then again a straight line from R2 to P2. Resulting in the shortest path from P1 to P2.

To get this, you have to make some assumptions, in the article I got from my Professor, was stated that the 4 points should be assumed to be in one plane.

You stated that it was easy to show the above, did you mean my code, or your principle?

Best wishes, Els
```
 0
Els
6/5/2010 9:34:05 AM
```"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hud5md\$sg3\$1@fred.mathworks.com>...
> Dear Roger,
>
> Thanks for your complete response, very clear.
>
> P1 and P2 are no points on a cylinder, they are random points in space. The R1 and R2 are tangent points on the surface of the cylinder. So there will come a straight line from P1 to R1, then a  geodesic over the surface of the cylinder from R1 to R2 and then again a straight line from R2 to P2. Resulting in the shortest path from P1 to P2.
>
> To get this, you have to make some assumptions, in the article I got from my Professor, was stated that the 4 points should be assumed to be in one plane.
>
> You stated that it was easy to show the above, did you mean my code, or your principle?
>
> Best wishes, Els
- - - - - - - - -
If R1 and R2 are points on the cylinder connected by a geodesic curve on the cylinder, and if P1R1 and P2R2 are lines pointing out along the two geodesic tangent directions at R1 and R2, respectively, then no, the two lines P1R1 and P2R2 cannot be assumed to lie in a common plane.  Go to the extreme of a helix that has wrapped halfway around the cylinder and you can easily see that this would be impossible.

In answer to your question, I was stating that, given the principle that a geodesic curve's normal direction should always be parallel to the surface normal at every point, then the geodesic on a cylinder must be a helix.  That is comparatively easy to show.  I confess I didn't study your code carefully because the initial assumption of a common plane was apparently false.

To show that the geodesic must be a helix, do the following.  Cartesian coordinates, r, t, and z (where t is short for theta) are related to cartesian coordinates by:

x = r*cos(t)
y = r*sin(t)
z = z

Letting s be the arc length along a geodesic, we have

ds = sqrt(dx^2 + dy^2 + dz^2) = sqrt((dx/dz)^2+(dy/dz)^2+1)*dz =
sqrt(r^2*(dt/dz)^2+1)*dz

Therefore the unit tangent is

[dx/ds , dy/ds , dz/ds] = [dx/dz , dy/dz , 1]*dz/ds =
[dx/dz , dy/dz , 1]*/sqrt(r^2*(dt/dz)^2+1) =
[-r*sin(t)*dt/dz , r*cos(t)*dt/dz , 1]/sqrt(r^2*(dt/dz)^2+1)

The geodesic normal is the derivative with respect to s of this unit tangent vector, and it must always be parallel to the surface normal which is [cost) , sin(t) , 0].  The derivative of the third component of the unit tangent will have a factor of d(dt/dz)/dz and since this component is zero in the surface normal, this second derivative must also be zero.  That shows that dt/dz is a constant, and t is therefore a linear function of z.  Hence the geodesic must be a helix.

Roger Stafford
```
 0
Roger
6/5/2010 5:06:04 PM
```Dear Roger,

The assumption of the fact that P1,P2,R1 and R2 are in one plane is correct, because is is about muscle wrapping, and the muscles do not wrap around the whole cylinder, but only have real small part arc length on the cylinder its self. I tried to understand your explanation further, but I must admit, that it is too difficult for me. The dx, dy, dz....are that the distances between R1 and R2? But I do not have R1 and R2, I want to calculate them. The only distance I know on forehand is the distance between P1 and P2, and the position and direction of the axis of the cylinder.

Hope I do not frustrate you, because I do appriciate your help very much. So thanks again a lot.

Best Wishes,

Els
```
 0
Els
6/5/2010 7:34:03 PM
```"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hue8rb\$80i\$1@fred.mathworks.com>...
> Dear Roger,
>
> The assumption of the fact that P1,P2,R1 and R2 are in one plane is correct, because is is about muscle wrapping, and the muscles do not wrap around the whole cylinder, but only have real small part arc length on the cylinder its self. I tried to understand your explanation further, but I must admit, that it is too difficult for me. The dx, dy, dz....are that the distances between R1 and R2? But I do not have R1 and R2, I want to calculate them. The only distance I know on forehand is the distance between P1 and P2, and the position and direction of the axis of the cylinder.
>
> Hope I do not frustrate you, because I do appriciate your help very much. So thanks again a lot.
>
> Best Wishes,
>
> Els

No, I must firmly disagree with you.  Even with a very small travel along a helix, the two lines R1P1 and R2P2 will not be coplanar.  I only asked you to envision things halfway around because the statement should be glaringly obvious in that case.  If you doubt this, try plugging some numbers into those formulae I gave you for a much smaller separation and check if the four points are coplanar.  Four points, [x1,y1,z1], [x2,y2,z2], [x3,y3,z3], [x4,y4,z4], are coplanar if and only if the determinant

det([x1,x2,x3,x4;y1,y2,y3,y4;z1,z2,z3,z4;1,1,1,1])

is equal to zero.

However, your intuition ought to tell you that if the points are not colinear when the rotation is half way around the cylinder, that situation is not likely to suddenly clear up even if they are only part way around the cylinder.  Actually they can't be coplanar unless P1 and P2 are the same points, or else the angle around the helix is some multiple of 2*pi (360 degrees).

Roger Stafford
```
 0
Roger
6/5/2010 10:32:04 PM
```Ok, and now I do have the two points on the surface. How do I calculate the distance between them? Because the problem is that my cylinder is rotated along a line, and that I can't use the calculation below.

------------------------

If P1 and P2 have rectangular coordinates (x1,y1,z1) and (x2,y2,z2), their cylindrical
coordinates are (r*cos(t1),r*sin(t1),z1) and (r*cos(t2),r*sin(t2),z2), where t1 and t2 are in
radians.  If we now unroll the cylinder onto the plane x=r, P1 goes to Q1 and P2 goes to
Q2.  The rectangular coordinates of Q1 and Q2 are
Q1=(r,r*t1,z1)
Q2=(r,r*t2,z2)

The distance d from P1 to P2, on the cylinder, is equal to the distance from Q1 to Q2,
namely,

d = sqrt[ (r*t1 - r*t2)^2 + (z1-z2)^2 ]
-------------------------
```
 0
Els
6/6/2010 2:23:04 PM
```"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hugb08\$ns6\$1@fred.mathworks.com>...
> Ok, and now I do have the two points on the surface. How do I calculate the distance between them? Because the problem is that my cylinder is rotated along a line, and that I can't use the calculation below.
>
> ------------------------
>
> If P1 and P2 have rectangular coordinates (x1,y1,z1) and (x2,y2,z2), their cylindrical
> coordinates are (r*cos(t1),r*sin(t1),z1) and (r*cos(t2),r*sin(t2),z2), where t1 and t2 are in
> radians.  If we now unroll the cylinder onto the plane x=r, P1 goes to Q1 and P2 goes to
> Q2.  The rectangular coordinates of Q1 and Q2 are
> Q1=(r,r*t1,z1)
> Q2=(r,r*t2,z2)
>
> The distance d from P1 to P2, on the cylinder, is equal to the distance from Q1 to Q2,
> namely,
>
> d = sqrt[ (r*t1 - r*t2)^2 + (z1-z2)^2 ]
> -------------------------

Yes, that last formula is indeed the best way to calculate the kind of cylindrical "distance" you are seeking.  It should however be modified a bit.  Something like this:

d = sqrt((r*(mod(t2-t1+pi,2*pi)-pi))^2+(z1-z2)^2)

The unrolled cylinder has to be conceived as having at least two other unrolled copies rather than just the one extending from 0 to 2*pi*r, and you need to select the shortest straight line path between the copies of P1 and P2.  The mod operation above does that.

As for making use of this method, you first have to do a translation of some point on your axis line over to the origin and bring P1 and P2 along with it.  Then you have to rotate this translated axis so as to lie along the z-axis, again rotating P1 and P2 along with it.  I think you've had advice in the past about how to do this rotation.  Next you convert the new P1 and P2 cartesian coordinates to cylindrical coordinates.  Then finally you can use the above method.

By the way, this method of unrolling a cylinder should be a powerful indicator that indeed the geodesics on circular cylinders are the helices (except of course when the angles are equal and they are straight lines, or when the z coordinates are equal and they are circles.)

Roger Stafford
```
 0
Roger
6/6/2010 4:46:03 PM
```"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hugjcb\$23v\$1@fred.mathworks.com>...
> "Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <hugb08\$ns6\$1@fred.mathworks.com>...
> > Ok, and now I do have the two points on the surface. How do I calculate the distance between them? Because the problem is that my cylinder is rotated along a line, and that I can't use the calculation below.
> >
> > ------------------------
> >
> > If P1 and P2 have rectangular coordinates (x1,y1,z1) and (x2,y2,z2), their cylindrical
> > coordinates are (r*cos(t1),r*sin(t1),z1) and (r*cos(t2),r*sin(t2),z2), where t1 and t2 are in
> > radians.  If we now unroll the cylinder onto the plane x=r, P1 goes to Q1 and P2 goes to
> > Q2.  The rectangular coordinates of Q1 and Q2 are
> > Q1=(r,r*t1,z1)
> > Q2=(r,r*t2,z2)
> >
> > The distance d from P1 to P2, on the cylinder, is equal to the distance from Q1 to Q2,
> > namely,
> >
> > d = sqrt[ (r*t1 - r*t2)^2 + (z1-z2)^2 ]
> > -------------------------
>
>   Yes, that last formula is indeed the best way to calculate the kind of cylindrical "distance" you are seeking.  It should however be modified a bit.  Something like this:
>
>  d = sqrt((r*(mod(t2-t1+pi,2*pi)-pi))^2+(z1-z2)^2)
>
> The unrolled cylinder has to be conceived as having at least two other unrolled copies rather than just the one extending from 0 to 2*pi*r, and you need to select the shortest straight line path between the copies of P1 and P2.  The mod operation above does that.
>
>   As for making use of this method, you first have to do a translation of some point on your axis line over to the origin and bring P1 and P2 along with it.  Then you have to rotate this translated axis so as to lie along the z-axis, again rotating P1 and P2 along with it.  I think you've had advice in the past about how to do this rotation.  Next you convert the new P1 and P2 cartesian coordinates to cylindrical coordinates.  Then finally you can use the above method.
>
>   By the way, this method of unrolling a cylinder should be a powerful indicator that indeed the geodesics on circular cylinders are the helices (except of course when the angles are equal and they are straight lines, or when the z coordinates are equal and they are circles.)
>
> Roger Stafford

Dear Roger,

For my understanding, you are implying that I have to rotate and translate the axis of the cylinder. And do the same with points P1 and P2. Then I can calculate the cylindrical coordinates and use your formula.

Concerning the ellipsoids, you gave me the following code to do the rotation:
-------------------------------------
% Calculate necessary rotation data
t = t/norm(t); % Make t a unit vector
u = [1,0,0]; % Unit vector along the positive x-axis
w = cross(u,t); % The axis of rotation
cosa = dot(u,t); % Cosine of the angle of rotation
sina = norm(w); % Sine of the angle of rotation
w = w/sina; % Make w a unit vector
v = cross(w,u); % Unit vector orthogonal to u and w

% Proceed with the rotation transformation
n = numel(x); % Number of points in surface
p = [x(:),y(:),z(:)]; % Create n x 3 array of ellipsoid surface points
u = repmat(u,n,1); % Match sizes of u, v, and w to that of p
v = repmat(v,n,1);
w = repmat(w,n,1);
wp = cross(w,p,2);
q = dot(p,w,2)*w + cosa*cross(wp,w,2) + sina*wp; % Rotate p
xr = q(:,1); yr = q(:,2); zr = q(:,3); % Extract components
xr = reshape(xr,size(x)); % These are the rotated coordinates
yr = reshape(yr,size(x)); % in mesh format again
zr = reshape(zr,size(x));

% Now plot the rotated surface
surf(xr,yr,zr)
---------------------------

But if I am correct, I have to rotate now over 3 axis, to get the axis of my cylinder on the origin. The translation is simple, but the rotation afterwards I still don't get in Matlab. As for the above code, t should be the vector from the origin to the point on the z-axis which corresponds with the length of the cylinder (t=[0 0 LengthCylinder]) and u should be [0,0,1], still my cylinder is pointed along the original (but translated) axis.

Best Wishes,

Els
```
 0
6/6/2010 5:29:04 PM