Vectors on ellipsoid surface
in [Matlab]
|
Prev: calling value of variable from main program into myfun while using fsolve
Next: Matlab GUI: jump from one callback to another
From: Els on 24 May 2010 08:40 Dear Bruno, Really did my best to understand your code and your explenations, but I am really loosing it now. I just don't get what you mean and how I should do this. My input values are P1 P2 [xr,yr,zr] % which are the surface points of the ellipsoid, which is rotated along the z axis. The z axis is along vector t. [cx,cy,cz] % the center of the ellipsoid [ax,ay,az] % length of the axis in the unrotated x,y,z direction C % point on line P1-P2 Output should be c' on the surface of my ellipsoid. I have the feeling that we are so close to a solution, and really hope you could help me with a code. That would be great. Or just a starter, because you allready have done so much. Best wishes, Els
From: Bruno Luong on 24 May 2010 11:08 You want C' on (E) and CC' orthogonal to P1P2, right? Let us denote by (P) the plane passing through C and perpendicular to vector P1P2. This plane must contain C' (because you want CC' to be orthogonal to P1P2) and C as well. This plane (P) is 2D affine space. It intersects the 3D ellipsoid (E) in a ellipse contour. Let us call the contour (F). So let's summarize: There is a 2D-plane (P) that contains a point C and a ellipse contour (F). Everything happens in this plane. If you take the projection of C on (F) and call it C', it verifies everything you want. No? If yes, the problem is reduces to 2D problem. How can you project on (P) and reduce mathematically to two coordinates? By selecting an orthonormal basis of (P), like this: B = null((P1(:)-P2(:))'). B is 3 x 2, each columns is a vector // to (P). B is orthonormal. The projection any point in R^3 in (P) can be carried out by the conversion formulas: x -> xi = B'*(x - C) (eqt1); and the reverse R^2 -> R^3 xi -> x = C + B*xi (eqt2). Assuming the ellipsoid equation (E) is { x such that x'*A*x + b'*x + c = 0 } (eqt3) Injecting (eqt2) in (eqt3), we get: xi'*B'*A*B*xi + 2*C'*A*B*xi + C'*A*C + b'*B*xi + b'*C + c = 0. This is nothing than the equation of the ellipse (F), the intersection of (E) and (P) in 2D xi-coordinates: { xi : xi'*A2*xi + b2'*xi + c2 = 0 } where A2 = B'*A*B, b2 = 2*B'*A*C + B'*b, c2 = C'*A*C + b'*C + c. The xi-coordinates of C is (0,0)' (from eqt1). Now we need to project C = (0,0)' on (F) using the method we have showed you. That gives the xi-coordinates of C'. Convert in 3D coordinates using (eqt2) and you are done. If something is not clear you should ask a specific question. I can't guess what is clear what is not or where you get stuck. Bruno
From: Bruno Luong on 24 May 2010 11:51 I clean up the code and better structure it. The convolution part is now carried out in an optimal manner (not at all important for 2D or 3D, but it's much nicer that way). The submission on FEX is here: http://www.mathworks.com/matlabcentral/fileexchange/27711-euclidian-projection-on-ellipsoid Bruno
From: Els on 24 May 2010 13:45 Dear Bruno, This is amazing, i must agree with Mark. Such a 'simple' question evaluates to a great contribution. I do not know how to thank you. My thesis on kinematics of ligaments in the knee will improve really much with this. I added this code to EllGeo2Alg radii = reshape(radii, 1, []); W = bsxfun(@rdivide, U, radii); A = W*W'; b = -2*A*x0; c = x0'*A*x0-1; B=null((P1(:)-P2(:))') ; C=[0,0]'; A2=B'*A*B; b2=2*B'*A*C + B'*b; c2=C'*A*C + b'*C + C; So that the output of the function is now A2,b2 and C2. So by doing the next code, I thought I would get C'. ------------------------------------------------------------------- [A2 b2 c] = EllGeo2Alg(radii, U, x0,P1,P2) x = EllPrj(c, radii, U); % Compute l = |x-c| l = sqrt(sum(bsxfun(@minus, x, c).^2, 1)); % Where the distance is minimal [lmin imin] = min(l); ---------------------------------------------------------------------- Only it allready goes wrong in b2=2*B'*A*C + B'*b; Should I transform C(0,0) to 3D before I do this? Or is my whole order of doing things wrong? By transforming C(0,0) to 3D I should use eq 2 (C+ B*xi), where I do know B and xi, but what parameter should I choose for C? And how can your file in the file exchanger help me with this? Hope these questions are specific enough. Again, thanks for all the effort and help. It is getting much clearer now. Best wishes, Els
From: Bruno Luong on 24 May 2010 14:09
Here are few mistakes I catch: > > I added this code to EllGeo2Alg > > radii = reshape(radii, 1, []); > W = bsxfun(@rdivide, U, radii); > A = W*W'; > b = -2*A*x0; > c = x0'*A*x0-1; > > B=null((P1(:)-P2(:))') ; > C=[0,0]'; % <- THIS IS WRONG, you should use 3D coordinates > A2=B'*A*B; > b2=2*B'*A*C + B'*b; > c2=C'*A*C + b'*C + C; % <- This is wrong you should use the lower-case c (of A, b c) It should goes like this (not check): radii = reshape(radii, 1, []); W = bsxfun(@rdivide, U, radii); A = W*W'; b = -2*A*x0; c = x0'*A*x0-1; B=null((P1(:)-P2(:))') ; A2 = B'*A*B; b2 = 2*B'*A*C+B'*b; c2 = C'*A*C + b'*C + c; [radii U x0] = EllAlg2Geo(A2, b2, c2); C2 = [0; 0]; % xi-coordinates of C Cp2 = EllPrj(C2, radii, U, x0); % C' l = sqrt(sum(Cp2.^2, 1)); % Where the distance is minimal [lmin imin] = min(l); Cp2 = Cp2(:,imin); CP = C + B*Cp2; % <- This is the 3D coordinates of the projection C' % Bruno |