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 23 May 2010 02:33
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hta8lc$e4f$1(a)fred.mathworks.com>...
>
>
> 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
>

That's a way to go Roger. Very good that you though about this. Bravo.

Bruno
From: Roger Stafford on 23 May 2010 02:37
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hta8lc$e4f$1(a)fred.mathworks.com>...
> ........
> 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
> ........

I erred in writing the equation explaining the ratio, t. It should read:

"The significance of the t variable is that it is the common ratio

(p-x)/(x/a^2) = (q-y)/(y/b^2) = (r-z)/(z/c^2) = t "

This error does not affect the code itself, just the explanation.

Roger Stafford
From: Bruno Luong on 23 May 2010 05:01
I have extended Roger's idea of parametrization with respect to Lagrange parameter (he call it "t", I call it "lambda") to n-dimensional ellipsoid. The polynomial to be solved is 2*n in order.

Here is the code:

%%%

% 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 = 1*randn(n,1);

%%
[U S V] = svd(Q);
s = diag(S);
radii = sqrt(1./s);
A = U*diag(radii);

%%
% Engine

% Change variable
% c = U*d
% x = U*y
d = U'*c;
% solve for lambda
% d-y = lambda*S*y
% y'*S*y = 1

% This are polynomial for (s*lambda+1)^2
pd = [s.^2 2*s ones(size(s))];
a = s.*d.^2;
P = 0;
for i=1:n
% Pi := Prod pd(j) for j#i
Pi = 1;
for j=1:n
if j~=i
Pi = conv(Pi, pd(j,:));
end
end
% P = sum(Pi)
P = P + a(i)*Pi;
end
% PC = Prod pd(j) for j#i
Pc = 1;
for j=1:n
Pc = conv(Pc, pd(j,:));
end
P = [0 0 P]-Pc;
lambda = roots(P);
lambda = reshape(lambda, 1, []);
lambda = lambda(imag(lambda)==0);

% Compute: y = d / (lambda*s + 1)
ls = bsxfun(@times, lambda, s);
y = bsxfun(@rdivide, d, ls+1);
x = U*y;

% Compute l = |x-c|
l = sqrt(sum(bsxfun(@minus, x, c).^2, 1));
% Where the distance is minimal
[lmin imin] = min(l);
lmin

%%
% Graphical check
% Generate many points on ellipse
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');
for k=1:size(x,2)
if k==imin
linespec = '-r.';
else
linespec = '-c.';
end
plot([x(1,k) c(1)],[x(2,k) c(2)],linespec);
end
elseif n==3
N = 64;
[X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N);
XYZ = U*[X(:) Y(:) Z(:)]';
min(sqrt(sum(bsxfun(@minus,XYZ,c).^2,1)))
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,'Edgecolor','none','FaceColor','flat');
colormap gray;
axis equal;
hold on;
plot3(c(1),c(2),c(3),'ok');
% Plot c projected to the ellipse
for k=1:size(x,2)
if k==imin
linespec = '-r.';
else
linespec = '-c.';
end
plot3([x(1,k) c(1)],[x(2,k) c(2)],[x(3,k) c(3)],linespec);
end
else
fprintf('Sorry I can''t plot\n');
end

% Bruno

There is a detail about the convolution part is somewhat redundant coded, since there is a lot of overlapping of polynomial product. I have not though how to do it in an optimal manner.

Bruno
From: Els on 23 May 2010 09:21
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <htaqsk$5h1$1(a)fred.mathworks.com>...
> I have extended Roger's idea of parametrization with respect to Lagrange parameter (he call it "t", I call it "lambda") to n-dimensional ellipsoid. The polynomial to be solved is 2*n in order.
>
> Here is the code:
>
> %%%
>
> % 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 = 1*randn(n,1);
>
> %%
> [U S V] = svd(Q);
> s = diag(S);
> radii = sqrt(1./s);
> A = U*diag(radii);
>
> %%
> % Engine
>
> % Change variable
> % c = U*d
> % x = U*y
> d = U'*c;
> % solve for lambda
> % d-y = lambda*S*y
> % y'*S*y = 1
>
> % This are polynomial for (s*lambda+1)^2
> pd = [s.^2 2*s ones(size(s))];
> a = s.*d.^2;
> P = 0;
> for i=1:n
> % Pi := Prod pd(j) for j#i
> Pi = 1;
> for j=1:n
> if j~=i
> Pi = conv(Pi, pd(j,:));
> end
> end
> % P = sum(Pi)
> P = P + a(i)*Pi;
> end
> % PC = Prod pd(j) for j#i
> Pc = 1;
> for j=1:n
> Pc = conv(Pc, pd(j,:));
> end
> P = [0 0 P]-Pc;
> lambda = roots(P);
> lambda = reshape(lambda, 1, []);
> lambda = lambda(imag(lambda)==0);
>
> % Compute: y = d / (lambda*s + 1)
> ls = bsxfun(@times, lambda, s);
> y = bsxfun(@rdivide, d, ls+1);
> x = U*y;
>
> % Compute l = |x-c|
> l = sqrt(sum(bsxfun(@minus, x, c).^2, 1));
> % Where the distance is minimal
> [lmin imin] = min(l);
> lmin
>
> %%
> % Graphical check
> % Generate many points on ellipse
> 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');
> for k=1:size(x,2)
> if k==imin
> linespec = '-r.';
> else
> linespec = '-c.';
> end
> plot([x(1,k) c(1)],[x(2,k) c(2)],linespec);
> end
> elseif n==3
> N = 64;
> [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N);
> XYZ = U*[X(:) Y(:) Z(:)]';
> min(sqrt(sum(bsxfun(@minus,XYZ,c).^2,1)))
> 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,'Edgecolor','none','FaceColor','flat');
> colormap gray;
> axis equal;
> hold on;
> plot3(c(1),c(2),c(3),'ok');
> % Plot c projected to the ellipse
> for k=1:size(x,2)
> if k==imin
> linespec = '-r.';
> else
> linespec = '-c.';
> end
> plot3([x(1,k) c(1)],[x(2,k) c(2)],[x(3,k) c(3)],linespec);
> end
> else
> fprintf('Sorry I can''t plot\n');
> end
>
> % Bruno
>
> There is a detail about the convolution part is somewhat redundant coded, since there is a lot of overlapping of polynomial product. I have not though how to do it in an optimal manner.
>
> Bruno

Dear Bruno,

Again, thanks for all the effort, this is great. Now I am trying to incorporate it into my own program. Therefore I am trying to understand your code, which is really advanced for me. But what I am not getting is how to give in P1 and P2 (the two points between which C lies), and is the projected C in 90 degrees of this line?

What should I do if my origin is not on [0,0,0], just change [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N); ?

And what should I do if my axis are not straigt but rotated?

Really hope you can help me with this final step.

Best wishes,

Els
From: Mark Shore on 23 May 2010 09:51
I have to say results can be very impressive when a challenging problem catches the interest of CSSM contributors.
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5 6
Prev: calling value of variable from main program into myfun while using fsolve
Next: Matlab GUI: jump from one callback to another