bj,
- Code: Tout sélectionner
tu détermines facilement les éléments de la conique
c'est pour moi le plus difficile...je développe, pe y a til plus simple?
à partir de le rep matricielle
X'AX -2B'X + c=0 (A dim 2 def positive et sym)
donc avec A et B à trouver
si tout se passait bien on aurait
X'DX = 1 (avec D(1,1)=1/a^2, D(2,2)=1/b^2)
donc j'impose c = 1, et A et B sont scalés de ce facteur c
translation:
soit Y + U= X
developpant
(1) Y'AY
(2) 2Y'AU - 2B'Y
(3) U'AU -2B'U + c
posant U = inv(A)B, (2) fait 0 et donc
avec
U est le vecteur translation (centre de lellipse)
en appliquant l'homothétie? de rapport -c2, on a
d'où
Puis on diagonalise pour trouver les axes:
et pour l'angle, en constatant que P (orthonormée..) est une matrice de rotation:
theta = atan2(D(1,2), D(1,1))
la méthode semble ok m'enfin..
- Code: Tout sélectionner
1;
clear all;
% generate an ellipse, with some rotation
% ---------------------------------------
a = 5;
b = 3;
xc = 22;%translation center
yc = 8;
theta = pi/3;%towards a axis
function [ xs,ys ] = genEllipse(a,b,xc,yc, theta)
ef = sqrt(1-b^2/a^2)
c = a*ef
npoints = 200
angles = (1:npoints)/npoints*2*pi;
r = a*(1-ef^2)./(1+ef*cos(angles));
xRef = c+r.*cos(angles);
yRef = r.*sin(angles);
xRefRot = cos(theta)*xRef - sin(theta)*yRef;
yRefRot = sin(theta)*xRef + cos(theta)*yRef;
xs = xc + xRefRot;
ys = yc + yRefRot;
%plot(xRef, yRef,'o', xs, ys, 'x')
endfunction
[ xs, ys ] = genEllipse(a,b,xc,yc,theta);
%% add some noise onto some points
% ---------------------------------------
nRegPoints = 25
indices = (1:nRegPoints)*floor(length(xs)/nRegPoints);
percentError=0.1%
points = [xs(indices); ys(indices)];
noiseScal=1+(rand(1, nRegPoints)*2-1)*percentError;
noisePoints=[
(points(1,:)-xc).*noiseScal+xc;
(points(2,:)-yc).*noiseScal+yc
]';
%plot(points(1,:), points(2,:),'o',noisePoints(:,1), noisePoints(:,2) '.','x')
% regress the parameters from those
% ---------------------------------------
x0 = [0,0,0,0,0,1];
function obj=phi(x, points)
A = [
x(1),x(2);
x(2),x(3)
];
B = [x(4),x(5)];
c = x(6);
r=[];
for i = 1:length(points)
p = points(i,:);
r(i) = p*A*p' - 2*B*p' -1;%assumes, c!=0
end
obj = sumsq(r);
endfunction
[x, obj, info, iter, nf, lambda] = sqp(x0, @(x) phi(x, noisePoints))
% get back the original ellipse
% ---------------------------------------
A = [
x(1),x(2);
x(2),x(3)
];
B = [x(4),x(5)]';
c = x(6);
U = inv(A)*B;
c2 = -1-B'*U;
A = A/(-c2);
[P,D] = eig(A) % ordered decreasing value. must be def positive
aBack = 1/sqrt(D(1,1))
bBack = 1/sqrt(D(2,2))
thetaBack = atan2(P(1,2),P(1,1))
xcBack = U(1)
ycBack = U(2)
% compare
% ---------------------------------------
[xs2, ys2] = genEllipse(aBack, bBack, xcBack, ycBack, thetaBack);
plot(xs, ys, 'x', xs2, ys2, 'o', noisePoints(:,1), noisePoints(:,2), '.')
ex:
bleu, référence,
rouge: les points à interpoler
vert: ellipse interpolante