Prev: Growing matrix
Next: Loading a large dataset from a database-- java.lang.OutOfMemoryError: Java heap space
From: Els on 3 Jun 2010 09:06
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
r=5; % Radius
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];
From: Els on 4 Jun 2010 05:00
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)];
From: Els on 4 Jun 2010 11:27
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.
From: Roger Stafford on 5 Jun 2010 03:12
"Els " <y.e.t.reeuwijk(a)student.utwente.nl> wrote in message <hu89bs$ito$1(a)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
From: Els on 5 Jun 2010 05:34
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
 |  Next  |  Last
Pages: 1 2 3
Prev: Growing matrix
Next: Loading a large dataset from a database-- java.lang.OutOfMemoryError: Java heap space