Approximation of an ellipse on curvature motion
13 views (last 30 days)
Show older comments
Hello,
I am trying to approximate an ellipse (centered at (0.5,0.5) with some major and minor axis Rx and Ry) that is shrinking by curvature motion and it seems that I have some trouble in approximating the proper curvature and/or possibly the minor and major axis. (I am really confident for the normal unit vector estimation).
The method that I am using to do it so far is the following
For the minor and major axis I use
[~,K] = min(abs(y-0.5));
[~,L] = min(abs(x-0.5));
Rx = sqrt((cpxvec(K)-0.5).^2+(cpyvec(K)-0.5).^2);
Ry = sqrt((cpxvec(L)-0.5).^2+(cpyvec(L)-0.5).^2);
in order to find the axis, where Rx is the axis in the x direction and Ry the one in the y direction.
For the curvature I use the formula
curv = -Rx*Ry./(((Ry*((x-0.5)/Rx)).^2+(Rx*((y-0.5)/Ry)).^2).^(3/2));
which is supposed to be the exact curvature.
The quantities x and y are found in the form (x,y)=(Rx*cos(t)+0.5,Ry*sin(t)+0.5) and the parameter t is computed with really high accuracy.
However, when I run the code, an ellipse with initial axis Rx=0.2 and Ry=0.1 is shrinking a little bit faster than the expected final time (which should be Rx*Ry/2). There seems to be some rounding error or another error that I could not find so far and so my question would be if there is any better way to approximate these values mentioned before.
Thank you.
(Below is the full coding)
Rx = 0.2;
Ry = 0.1;
Nx = 5*2^3;
Ny = Nx;
dx = 1/Nx;
dy = 1/Ny;
dt = 0.4*dx^2;
tend = Rx*Ry/2;
nt = round(tend/dt);
dt = tend/nt;
x = (0:dx:1)';
y = x;
[xx,yy] = meshgrid(x,y);
p = 3;
order = 2;
dim = 2;
bw = 1.0001*sqrt((dim-1)*((p+1)/2)^2 + ((order/2+(p+1)/2)^2));
xvec = xx(:);
yvec = yy(:);
cpxvec = zeros(length(xvec),1);
cpyvec = cpxvec;etaX = cpxvec;etaY = cpxvec;curv = cpxvec;dist = cpxvec;
parfor j=1:length(xvec)
t = fminsearch(@(t) 1/2*sqrt((xvec(j)-0.5-Rx*cos(t)).^2+(yvec(j)-0.5-Ry*sin(t)).^2),0,...
optimset('TolX',1e-11));
cpxvec(j) = Rx*cos(t)+0.5;
cpyvec(j) = Ry*sin(t)+0.5;
etaX(j) = -Ry*cos(t)/sqrt(Ry^2*(cos(t))^2+Rx^2*(sin(t))^2);
etaY(j) = -Rx*sin(t)/sqrt(Ry^2*(cos(t))^2+Rx^2*(sin(t))^2);
curv(j) = -Rx*Ry/((Ry^2*(cos(t))^2+Rx^2*(sin(t))^2)^(3/2));
dist(j) = sqrt((xvec(j)-cpxvec(j))^2+(yvec(j)-cpyvec(j))^2);
end
temp = sqrt(etaX.^2+etaY.^2);
etaX = etaX./temp;
etaY = etaY./temp;
band = find( abs(dist)<=bw*dx);
cpxvec = cpxvec(band);
cpyvec = cpyvec(band);
etaX = etaX(band);
etaY = etaY(band);
curv = curv(band);
xb = xvec(band);
yb = yvec(band);
figure(1)
plot(xx(band),yy(band),'bo',cpxvec,cpyvec,'r+')
axis([0.2 0.8 0.2 0.8])
% hold on
% quiver(cpxvec,cpyvec,-curv.*etaX,-curv.*etaY)
% hold off
errorEx = zeros(nt,1);
ind = 1;
R = [Rx Ry];
dA = [];
array2 = [];
for T = 1:nt
cpxvec = cpxvec-dt*curv.*etaX;
cpyvec = cpyvec-dt*curv.*etaY;
[~,K] = min(abs(cpyvec-0.5));
[~,L] = min(abs(cpxvec-0.5));
Rx = sqrt((cpxvec(K)-0.5).^2+(cpyvec(K)-0.5).^2);
Ry = sqrt((cpxvec(L)-0.5).^2+(cpyvec(L)-0.5).^2);
R = [R;Rx Ry];
errorEx(T) = norm(sqrt(((cpxvec-0.5)/Rx).^2+((cpyvec-0.5)/Ry).^2)-1,inf);
A = pi*Rx*Ry;
dA = [dA;((pi*0.2*0.1)-A)/(T*dt)];
cpxvec = zeros(length(xvec),1);
cpyvec = cpxvec;etaX = cpxvec;etaY = cpxvec;curv = cpxvec;dist = cpxvec;
parfor j=1:length(xvec)
t = fminsearch(@(t) 1/2*sqrt((xvec(j)-0.5-Rx*cos(t)).^2+(yvec(j)-0.5-Ry*sin(t)).^2),0,...
optimset('TolX',1e-14));
cpxvec(j) = Rx*cos(t)+0.5;
cpyvec(j) = Ry*sin(t)+0.5;
etaX(j) = -Ry*cos(t)/sqrt(Ry^2*(cos(t))^2+Rx^2*(sin(t))^2);
etaY(j) = -Rx*sin(t)/sqrt(Ry^2*(cos(t))^2+Rx^2*(sin(t))^2);
curv(j) = -Rx*Ry/((Ry^2*(cos(t))^2+Rx^2*(sin(t))^2)^(3/2));
dist(j) = sqrt((xvec(j)-cpxvec(j))^2+(yvec(j)-cpyvec(j))^2);
end
temp = sqrt(etaX.^2+etaY.^2);
etaX = etaX./temp;
etaY = etaY./temp;
band = find( abs(dist)<=bw*dx);
cpxvec = cpxvec(band);
cpyvec = cpyvec(band);
etaX = etaX(band);
etaY = etaY(band);
curv = curv(band);
xb = xvec(band);
yb = yvec(band);
figure(1)
subplot(1,3,1)
plot(xx(band),yy(band),'bo',cpxvec,cpyvec,'r+')
axis([0.1 0.9 0.1 0.9])
title(['The ellipse at time t=' num2str(T*dt)])
% hold on
% quiver(cpxvec,cpyvec,etaX,etaY)
% hold off
subplot(1,3,2)
plot((0:T)*dt,R(:,1),'b',(0:T)*dt,R(:,2),'r')
title('The major and minor axis');
xlabel('t');ylabel('Rx , Ry');
figure(1)
subplot(1,3,3)
plot((1:T)*dt,dA,'b',(1:T)*dt,2*pi*ones(T,1),'r')
title('The derivative of the area over time');
xlabel('t');ylabel('The derivative of the area');
end
0 Comments
Answers (1)
Image Analyst
on 2 Mar 2014
I have no idea what that code does (won't even run on my computer). But, if you, or anyone else, wants to find ellipses in an image, see this paper http://www.ecse.rpi.edu/homepages/qji/Papers/ellipse_det_icpr02.pdf
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!