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: Bruno Luong on 21 May 2010 13:35 Sorry make, it: % Data, generate ellipse centered at origin, x'*Q*x = c L=randn(2); Q = L'*L; c = rand; % Given Point (to be projected on ellipse) P = randn(2,1); % Engine [U S V] = svd(Q); s = diag(S); A = U*diag(sqrt(c./s)); R = [0 1; -1 0]; B = R*Q*A; alpha = -P.'*R*Q*A; rho = sqrt(s(1)/s(2)); beta = c*(rho-1/rho); Poly = beta*[0 -2 0 2 0] + ... alpha(1)*[-1 0 0 0 1] + ... alpha(2)*[0 2 0 2 0]; tprj = roots(Poly); tprj = tprj(imag(tprj)==0); tprj = reshape(tprj, 1, []); % Candidate of the projections xprj = A*[(1-tprj.^2); 2*tprj]; xprj = bsxfun(@rdivide, xprj, (1+tprj.^2)); % Graphical check % Generate point on ellipse theta = linspace(0,2*pi); x = A*[cos(theta); sin(theta)]; clf; plot(x(1,:),x(2,:)); axis equal; hold on; % Plot P projected to the ellipse plot(P(1),P(2),'or'); for k=1:size(xprj,2) plot([xprj(1,k) P(1)],[xprj(2,k) P(2)],'-r'); end % Bruno
From: Els on 21 May 2010 15:53 Thanks Bruno for all the effort. Though, my ellipsoid is a 3D figure. How difficult is it to adjust your code for 3D?
From: Bruno Luong on 21 May 2010 18:21 "Els " <y.e.t.reeuwijk(a)student.utwente.nl> wrote in message <ht6ob1$20d$1(a)fred.mathworks.com>... > Thanks Bruno for all the effort. Though, my ellipsoid is a 3D figure. How difficult is it to adjust your code for 3D? Oops, my bad I though it's 2D. It's unclear to to me where as the above method can be extended in 3D. Follow up the Lagrangian method proposed by Tortsen, allow me to synthesize the problem: Assuming the ellipsoid is E = { x in R^3 : x'*Q*x = 1 }, Q is (3 x 3) sdp matrix. c is the point to be projected. We want to minimize J(x) = | x - c |^2 such that x'*Q*x = 1 (eqt2) The Lagrangien is L(x) = J(x) + lambda * (x'*Q*x - 1) When you derive wrt lambda, you will find (eqt2). When you derive wrt x, you'll find (x - c) + lambda * Q*x = 0 or c = x + lambda*Q*x (eqt3) Q*x is the normal to the ellipsoid, and the relation (eqt3) is a Lagrange Euler condition. It says nothing more than: c is somewhere on the normal vector placed at x. Now, if we take the dot product of (eqt3) with x, and because x'*Q*x = 1 we find: <c,x> =|x|^2 + lambda, or: lambda = <c,x> - |x|^2. Replace the value of lambda in (eqt3), we get: c = x + (<c,x> - |x|^2) * Q*x This can be solved using FSOLVE. I would probably multiply by inv(Q) to make the problem better conditioned: (<c,x> - |x|^2) * x + P*x + d = 0 where P = inv(Q), and d := -Q\c. You could also try to solve the above with Newton, Gauss-Newton, Levenberg-Marquardt, etc... The function which is needed for FSOLVE is very easy to program (you however should organize the parameters so as to have access to c and Q): function f = myfun(x) f = (dot(x,c) - dot(x,x))*x + Q\(x-c) end Bruno
From: Bruno Luong on 22 May 2010 03:24 Here is the implementation of the idea I put earlier: % Data, generate ellipse centered at origin, { x: x'*Q*x = 1 } n = 3; % dimension L=randn(n); Q = L'*L; % Given point in R^n (to be projected on ellipse) c = randn(n,1); %% % Engine % Solve % (dot(x,c) - dot(x,x))*x + Q\(x-c) = 0 % x'*Q*x = 1 n = size(Q,1); x = c; Qinv = inv(Q); I = eye(n); maxiter = 100; tol = 1e-6; for k=1:maxiter d = c-x; g = dot(x,d); df = g*x - Q\d; M = x*(c-2*x)' + g*I + Qinv; P = null(x'*Q); y = (M*P)\df; xold = x; x = x-P*y; a = x'*Q*x; x = x/sqrt(a); dx = x-xold; if norm(dx)<tol break end end %% % Graphical check % Generate many points on ellipse [U S V] = svd(Q); s = diag(S); radii = sqrt(1./s); A = U*diag(radii); clf; if n==2 theta = linspace(0,2*pi,181); ellipse = A*[cos(theta); sin(theta)]; plot(ellipse(1,:),ellipse(2,:)); axis equal; hold on; % Plot c projected to the ellipse plot(c(1),c(2),'ok'); plot([x(1) c(1)],[x(2) c(2)],'-r.'); elseif n==3 N = 64; [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N); XYZ = U*[X(:) Y(:) Z(:)]'; X = reshape(XYZ(1,:),[N N]+1); Y = reshape(XYZ(2,:),[N N]+1); Z = reshape(XYZ(3,:),[N N]+1); surf(X,Y,Z); axis equal; hold on; % Plot c projected to the ellipse plot3(c(1),c(2),c(3),'ok'); plot3([x(1) c(1)],[x(2) c(2)],[x(3) c(3)],'-r.'); end
From: Roger Stafford on 22 May 2010 23:50
"Els " <y.e.t.reeuwijk(a)student.utwente.nl> wrote in message <ht5ldf$7tf$1(a)fred.mathworks.com>... > I plotted an ellipsoid [x,y,z]=ellipsoid (standard in Matlab) > > Now I have a point C(x,y,z) which lays within this surface. > I am looking for the projection of C on the surface of the ellipsoid,called C', where the normal vector perpendicular to the surface is in line with the vector between C and C'. > > Thanks in advance. The following is an analytic solution of your problem. It is not an iterative method but gets its solutions directly using matlab's 'roots' function. It assumes an ellipsoid with center at the origin and axes aligned with coordinates axes, satifying the equation x^2/a^2 + y^2/b^2 + z^2/c^2 = 1 The point (p,q,r) can be either inside or outside this ellipsoid. Points (x,y,z) on the ellipsoid are found such that (p,q,r) lies on the line through (x,y,z) orthogonal to the surface of the ellipsoid there. In general there will be six solutions to this problem though frequently some of these will be complex-valued. There are always at least two real solutions. It should be noted that if two of the three semi-axes a, b, and c, are equal, giving an oblate or prolate spheroid, then there are in general only five solutions, though again some of these may be complex too. It should be noted that in this case if (p,q,r) lies right on the axis of revolution, infinite rings of solutions about this axis become possible. In the code below in this case the algorithm still produces six solutions, but two of them will be equal - that is, equal to within round off error. The algorithm below computes the appropriate coefficients of a sextic polynomial in the variable t A*t^6 + B*t^5 + C*t^4 + D*t^3 + E*t^2 + F*t + G = 0 and then computes a set of (x,y,z) points on the surface for each root t according to the equations x = a^2*p/(a^2+t) y = b^2*q/(b^2+t) z = c^2*r/(c^2+t) The significance of the t variable is that it is the common ratio (p-x)/(x/p^2) = (q-y)/(y/q^2) = (r-z)/(z/r^2) = t between components of the line between (p,q,r) and (x,y,z) and the components of a line through (x,y,z) orthogonal to the surface there. The code assumes that a, b, c and p, q, r, have been defined. - - - - - - - - - - - - - - - - - - % The program % Compute some convenient intermediate values p2 = p^2; q2 = q^2; r2 = r^2; a2 = a^2; b2 = b^2; c2 = c^2; sbc = b2+c2; sca = c2+a2; sab = a2+b2; mbc = b2*c2; mca = c2*a2; mab = a2*b2; s1 = a2+b2+c2; s2 = mbc+mca+mab; s3 = a2*b2*c2; % Now for the polynomial's actual coefficients A = 1; B = 2*s1; C = s1^2+2*s2-p2*a2-q2*b2-r2*c2; D = 2*(s1*s2+s3-p2*a2*sbc-q2*b2*sca-r2*c2*sab); E = s2^2+2*s3*s1-p^2*a^2*(sbc^2+2*mbc) ... -q2*b2*(sca^2+2*mca)-r2*c2*(sab^2+2*mab); F = 2*s3*(s2-p2*sbc-q2*sca-r2*sab); G = s3*(s3-p2*mbc-q2*mca-r2*mab); % Now find the polynomial's six roots t = roots([A,B,C,D,E,F,G]); % Derive the corresponding sets of coordinates. x = a^2*p./(a^2+t); y = b^2*q./(b^2+t); z = c^2*r./(c^2+t); % The resulting three six-element column vectors give the solutions. - - - - - - - - - - - - - - - - - - It should be pointed out that, depending on the position of the (p,q,r) point relative to the ellipsoid, there can be a comparatively large effect from round off errors. For example, if the ellipsoid is nearly a sphere and the point (p,q,r) is nearly at the origin, all points on the surface come close to being solutions, and this would have the effect of magnifying any round off errors in the solutions actually arrived at. This difficulty is inherent in the nature of the stated problem. Roger Stafford |