Prev: i need matlab code for finding qrs interval in ecg signal
Next: need to find clustering efficiency
From: Gus Lott on 1 May 2010 10:17
I'm working on an image processing problem where I'm analyzing the boundary of an object that I've segmented.

I have the object's boundary sorted starting from the point closest to an angle of zero in the coordinate system and advancing to 2pi around the boundary. I then treat the boundary as a periodic signal (i.e. fourier descriptor technique) and I take the first component (describes the ellipse fit to the perimeter).

Then what I'm trying to do is determine the angle of orientation of that ellipse (phi). Of course, this ellipse will have a 180 degree symmetry on phi. But I can't seem to get it.

I've got a & b terms for the fourier series representation of the border in both the x and y component and it clearly fits the data when I re-build the boundary just from that component, but I think that the angle of the ellipse is somehow confounded with the angle at which I start the path around the perimeter in my mathematics. Does anyone know of a good walk through relating the first fourier descriptor to the ellipse that it encodes and how to calculate the phase angle of the ellipse in 2D space from the 1D fourier descriptors for x and y of the perimeter?
From: Gus Lott on 1 May 2010 10:56
I'll try to clarify:

I have a perimeter of an object in an image. p(x,y)

p(x,y) is periodic in both x and y as functions of the index of points in the perimeter as you travel around the perimeter.

I can subtract the average x and y position to remove the DC component of the perimeter. I then order the perimeter points by sorting them based on their angle relative to the object centroid (so they're nominally sorted in order around the perimeter from 0 to 2pi). This wouldn't work for a complex shape, but my objects are fairly close to ellipse to start with. This sort step is part of the edge detection since my perimeter is not necessarily in order when I obtain points from the edge detector.

The first fourier descriptor for x and y is defined by:

n=1;

Lx = numel(x)/2;
xp = 1:numel(x);
ax = sum(x.*cos(n*pi*xp/Lx))/Lx;
bx = sum(x.*sin(n*pi*xp/Lx))/Lx;

% and

Ly = numel(y)/2;
yp = 1:numel(y);
ay = sum(y.*cos(n*pi*yp/Ly))/Ly;
by = sum(y.*sin(n*pi*yp/Ly))/Ly;

Now, an ellipse can be regenerated from the data in the following fashion:

theta = linspace(-pi,pi,100);
fit_x =ax*cos(theta) + bx*sin(theta);
fit_y =ay*cos(theta) + by*sin(theta);

plot(fit_x,fit_y);

The trigonometry is frustrating me here. It's not clear to me how to back out the orientation of the ellipse based on these components, but it must be in there somewhere.

Any thoughts?
From: Bruno Luong on 1 May 2010 13:27
"Gus Lott" <lottg.nospam(a)janelia.hhmi.org> wrote in message <hrhfe6$e0h$1(a)fred.mathworks.com>...
> I'll try to clarify:
>
> I have a perimeter of an object in an image. p(x,y)
>
> p(x,y) is periodic in both x and y as functions of the index of points in the perimeter as you travel around the perimeter.
>
> I can subtract the average x and y position to remove the DC component of the perimeter. I then order the perimeter points by sorting them based on their angle relative to the object centroid (so they're nominally sorted in order around the perimeter from 0 to 2pi). This wouldn't work for a complex shape, but my objects are fairly close to ellipse to start with. This sort step is part of the edge detection since my perimeter is not necessarily in order when I obtain points from the edge detector.
>
> The first fourier descriptor for x and y is defined by:
>
> n=1;
>
> Lx = numel(x)/2;
> xp = 1:numel(x);
> ax = sum(x.*cos(n*pi*xp/Lx))/Lx;
> bx = sum(x.*sin(n*pi*xp/Lx))/Lx;
>
> % and
>
> Ly = numel(y)/2;
> yp = 1:numel(y);
> ay = sum(y.*cos(n*pi*yp/Ly))/Ly;
> by = sum(y.*sin(n*pi*yp/Ly))/Ly;
>
> Now, an ellipse can be regenerated from the data in the following fashion:
>
> theta = linspace(-pi,pi,100);
> fit_x =ax*cos(theta) + bx*sin(theta);
> fit_y =ay*cos(theta) + by*sin(theta);
>

% You can insert next in your program like this:

d = (ay*bx-ax*by);
H = [ay^2+by^2 -(ax*ay+bx*by); -(ax*ay+bx*by) ax^2+bx^2]/(d^2);
[V D]=eig(H);
d = diag(D);
% Two axis of the ellipse
W = V.*repmat(1./sqrt(d'),2,1)

plot(fit_x, fit_y);
hold on
plot([0 W(1,1)],[0 W(2,1)],'r')
plot([0 W(1,2)],[0 W(2,2)],'r')

% Bruno
From: Bruno Luong on 1 May 2010 13:57
Another leaner way is

A = [bx ax; by ay];
[U S V] = svd(A);
d = diag(S);
W = U.*repmat(d.',2,1)

plot(fit_x, fit_y);
hold on
plot([0 W(1,1)],[0 W(2,1)],'r')
plot([0 W(1,2)],[0 W(2,2)],'r')
axis equal

% Bruno
From: Gus Lott on 1 May 2010 15:19
Thanks a ton Bruno. I see THAT that works, now I gotta go figure out how :)

Cheers
 |  Next  |  Last
Pages: 1 2
Prev: i need matlab code for finding qrs interval in ecg signal
Next: need to find clustering efficiency