Prev: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)
Next: Using interp1 for timeseries data
From: Walter Roberson on 29 Apr 2010 21:53
GAURAV wrote:
> Hey Walter,
>
> Can you tell me how you solved these four equations?

Using a symbolic package not directly supported by Matlab:


F := (X,Y) -> (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) - 1;

#that's the basic equation, but it is more efficient to rationalize it
#and then get rid of the denominator

F1 := (X,Y) -> numer(normal(F(X,Y));

solve([F1(a1,b1),F1(a2,b2),F1(a3,b3),F1(a4,b4)],[a,b,Xc,Yc]);

Based upon the length of time it is taking, I suspect that when I did
this earlier, I must have mis-written the equations. I'll post later
when it comes up with an answer.
From: Bruno Luong on 30 Apr 2010 02:35
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hrd86k$fje$1(a)fred.mathworks.com>...

>
> Even if it runs near the origin, forcing F to be 1 will encounter awkward accuracy problems.

Roger, I don't think Roger it should not be any accuracy problem by fixing F=1. It is like divide the original implicit ellipsoid equation by F, and the linear equation to be solved is:

M = [x.^2,y.^2,x,y];
X = - M(:,1:4) \ ones(4,1)

A = X(1)
C = X(2)
D = X(3)
E = X(4)
F = 1

In contrary, it should be more accurate (if not equal), because the new matrix above is submatrix of the 5-column matrix, thus the ratio of singular values becomes smaller because it is the projection.

I have my full code to fit ellipsoid in ndimensional that I have posted long ago in the newsgroup ( http://www.mathworks.com/matlabcentral/newsreader/view_thread/168603 ) where I just revised for axis-alignement:

function [H X0 W] = solveellipse(X, axisaligned)
% function [H X0 W] = solveellipse(X)
%
% function [H X0 W] = solveellipse(X, axisaligned)
%
% INPUT:
% X : (m x n) : m data points in R^n, they are supposed to belong
% to an ellipsoide (or more generally a second-order
% implicite hyper-surface)
% if axisaligned is TRUE, elleipsoid is forced to be aligned with axis
%
% OUTPUT:
% H : Matrix of the Bilinear form associated with the ellipsoide
% Elippsoide = { X in R^n : (X-X0)' * H * (X-X0) = 1 }
% X0 : (n x 1), center of the ellipsoide
% W : (n x n) square matrix where each column is the axis of ellipsoide
%

[ndata ndim]=size(X);

if ndim<2
error('solveellipse: dimension number must be greater than 1');
end

% This vector will be used in few places for reshapping purpose
uno = ones(1,ndim);

%
% Generate all combinations of polynomial-order <=2 for ndim variables
%
order=cell(1,ndim);
order(:)={(0:2)};
ORDER=cell(1,ndim); % 1 x ndim cell of vectors (0:2)
[ORDER{:}]=ndgrid(order{:}); % Set {(0:2)}^ndim
ORDER=reshape(cat(ndim+1,ORDER{:}),[],ndim);
ORDER=ORDER(sum(ORDER,2)<=2,:); % second order only

%
% Remove constant term
%
ORDER=ORDER(2:end,:); % similar to ORDER(~any(ORDER,2),:)=[];

% Remove cross-term
if nargin>=2 && axisaligned
ORDER(sum(ORDER==1,2)==2,:) = [];
end

npol=size(ORDER,1); % number of basis 2nd order-polynomials


if ndata<npol
error('solveellipse: Not enough data');
end

function loc=findloc(order) % nested function
%
% Look for the location of the 'order' in the set of {basis
% 2nd-order polynomials}
%
[~, loc]=ismember(reshape(order,1,[]),ORDER,'rows');
end

%
% Both of these are 3d array of (ndata x ndim x npol)
%
BIGX=repmat(X,[1 1 npol]);
BIGORDER=repmat(reshape(ORDER',[1 ndim npol]),[ndata 1 1]);
M=squeeze(prod(BIGX.^BIGORDER,2)); % product of power in every dimensions
clear BIGX BIGORDER

%
% Solve for polynomial coefficiens that lead to 1 on input points
%
rhs=ones(ndata,1);
P=M\rhs;
clear M rhs;

%
% Extract the Hessian from the solution
%
I=repmat(eye(ndim),ndim,1);
I=reshape(I,ndim*[1 1 1]);
J=permute(I,[2 1 3]); % swap the first two dimensions
K=squeeze(mat2cell(I+J,uno,uno,ndim));
P2_loc=cellfun(@findloc,K);
P2 = zeros(size(P2_loc));
P2(P2_loc>0) = P(P2_loc>0); % second order term
A=(tril(P2)+triu(P2))/2; % devide by 2, accept diagonal terms

%
% Extract the gradient
%
I=mat2cell(eye(ndim),uno,ndim);
P1_loc=cellfun(@findloc,I);
P1=-0.5*P(P1_loc); % -0.5 * first order term

X0=A\P1;
lambda=1/(1+X0'*P1);
H=lambda*A;

if nargout>=3 % Compute main-axis of ellipsoide
[V D]=eig(H);
d=diag(D);
if any(d<0)
warning('solveellipse:npdH', ...
'solveellipse: non positive Hessian matrix');
end
W = V.*repmat(1./sqrt(d'),ndim,1);
end

end

% Bruno
From: Torsten Hennig on 29 Apr 2010 22:47
> Hello Eveyone,
>
> I am trying to get an ellipse passing through the
> four points I have which are not symmetric.
> So, I am using the equation (X-Xc)^2/(a^2) +
> (Y-Yc)^2/(b^2) =1
>
> I have four points which wil go in X and Y of the
> equation above.
> Let the points be (a1,b1), (a2,b2), (a3,b3) and
> (a4,b4)..
>
> How do I get a,b,Xc and Yc from this??
>
> Please help me out. Really appreciate..
>
> Thanks,
> Gaurav

The equation of the ellipse is given under

http://mathworld.wolfram.com/Ellipse.html

Best wishes
Torsten.
From: Torsten Hennig on 29 Apr 2010 22:46
> Hello Eveyone,
>
> I am trying to get an ellipse passing through the
> four points I have which are not symmetric.
> So, I am using the equation (X-Xc)^2/(a^2) +
> (Y-Yc)^2/(b^2) =1
>
> I have four points which wil go in X and Y of the
> equation above.
> Let the points be (a1,b1), (a2,b2), (a3,b3) and
> (a4,b4)..
>
> How do I get a,b,Xc and Yc from this??
>
> Please help me out. Really appreciate..
>
> Thanks,
> Gaurav

The eqation of the ellipse is given under

http://mathworld.wolfram.com/Ellipse.html

Best wishes
Torsten.
From: Bob Taylor on 30 Apr 2010 17:10
"GAURAV " <gsharda(a)engineering.uiowa.edu> wrote in message <hrcrds$7e0$1(a)fred.mathworks.com>...
> Hello Eveyone,
>
> I am trying to get an ellipse passing through the four points I have which are not symmetric.
> So, I am using the equation (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) =1
>
> I have four points which wil go in X and Y of the equation above.
> Let the points be (a1,b1), (a2,b2), (a3,b3) and (a4,b4)..
>
> How do I get a,b,Xc and Yc from this??
>
> Please help me out. Really appreciate..
>
> Thanks,
> Gaurav

You can set this up as a nonlinear optimization problem with fmincon. The variables to solve for are (a, b, Xc, Yc). You minimize the sum of the left side of your ellipse equation over each of the four points, set lower bounds for a and b to 0, and form 4 nonlinear equality constraints for each point (a1, b1), etc. to require the point to fall on the ellipse. If you use fmincon (interior-point method seems to work), an initial point is just (sqrt(2)*a1, sqrt(2)*b1, 0, 0) and fmincon will get a solution (if one exists since 4 arbitrary points may not fall on a single ellipse - consider three points that form a triangle with the fourth point inside the triangle).

Minimize
(a1 - Xc)^2/a^2 + (b1 - Yc)^2/b^2 + ... + (a4 - Xc)^2/a^2 + (b4 - Yc)^2/b^2
subject to
a >= 0
b >= 0
(a1 - Xc)^2/a^2 + (b1 - Yc)^2/b^2 = 1
....
(a4 - Xc)^2/a^2 + (b4 - Yc)^2/b^2 = 1.
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)
Next: Using interp1 for timeseries data