Approximation of an ellipse on curvature motion

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

Answer 1

[Image Analyst]

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

[Argyrios]

Thank you for noticing that the code had some errors. I fixed them and now it is fully working. However, this paper didn't really help me that much since it mentions techniques for detecting ellipses in an image. Even though I tried to adjust it for the specific case that I have didn't really work as expected..