Bruce LandGhosh and Jain (1) showed how to charactize a simple, closed, piece-wise linear curve as a Fourier series in a single parametric variable. Such a formulation allows (see Tsow and Land below):
- interpolation of curves by sampling the parameter more densely.
- morphing of one curve into another by locating matching parametric points on the two curves and smoothly interpolating the Fourier components.
- simplification of curves by removing Fourier components. This maintains the overall shape, but eliminates sharp corners.
- location of the curve in the plane as the DC components of the series.
- ability to "add" curves.
- ability to add band-limited noise to a curve.
- ability to connect serveral curves (possibly representing serial sections) to make a 3D object.
We intend to use the Fourier fit scheme as a way to simplify and characterize the outlines of moving animals. The animals will be videotaped. The videotape will then be analysed to extract the outline of the animal, Fourier fit the outline, then use the low bandwidth components to characterize the motion. For example, the DC components give the position of the animal. The fundamental frequency coefficients of the fit give the ellipse (2) which best fits the animal outline.
Code
A matlab program was written to animate a test object, extract its outline, then do the Fourier fit.
This program Fourier fits a square:
%Fourier fit of motion
clear all
figure(1)
clf; clc
set(gcf,'position',[50,600,400,200])
axis([0 400 0 200])
% axis tight
set(gcf,'units','pixels');
set(gca,'units','pixels');
w_pos = get(gcf, 'position');
set(gca, 'position', [0 0 w_pos(3) w_pos(4)]);
%number of Fourier components to use
%vary this to make a better fit
N = 11;
%make an rectangle
figure(1)
cla;axis off
rectangle('Position', [30 30 100 100]);
im = getframe(gcf);
ed = im.cdata;
% figure(2)
% imshow(ed)
%order the edge pixels
%find them in the array
[xpos,ypos] = find(ed(:,:,1)==0);
npts = length(xpos);
%keep track of visit order
ordering=[1:npts] * 0;
visited=[1:npts] * 0;
%pick the first point on the edge and calculate the normal
%move the next pixel which is:
%--on the edge
%--closeest
%now just move to the nearest pixel which has not been visited
nextindex=1; ordering(1)=1; visited(1)=1;
count = 1;
while (count<npts)
currentindex=nextindex;
currentmin=inf;
for i=1:npts
d = (xpos(currentindex)-xpos(i)).^2 + (ypos(currentindex)-ypos(i)).^2 ;
if (d<currentmin) & (visited(i)==0)
currentmin = d;
nextindex = i;
end
end
%at this point, we have the index to the nearest point
visited(nextindex)=1;
count = count + 1;
ordering(count)=nextindex;
end
%define a piece-wise linear, closed curve
clear p
p(:,1) = ypos(ordering);
p(:,2) = xpos(ordering);
K = length(p);
%close the loop
p(K+1,:) = p(1,:);
%get fourier transform of the closed curve
%Ghosh and Jain, IEEE conputer graphics and appl, 1994
dx = diff(p(:,1));
dy = diff(p(:,2));
dL = sqrt( dx.^2 + dy.^2 );
L = sum(dL);
Lcum = cumsum(dL);
%get a0, the x-offset
temp = p(1,1)*Lcum(1)+ dx(1)*Lcum(1)^2/dL(1)*0.5; %for i=1
for i=2:K
temp = temp + ...
(p(i,1)-dx(i)*Lcum(i-1)/dL(i))*(Lcum(i)-Lcum(i-1)) + ...
0.5*dx(i)/dL(i)*(Lcum(i)^2 - Lcum(i-1)^2);
end
a0=temp/L;
%get c0, the y-offset
temp = p(1,2)*Lcum(1)+ dy(1)*Lcum(1)^2/dL(1)*0.5; %for i=1
for i=2:K
temp = temp + ...
(p(i,2)-dy(i)*Lcum(i-1)/dL(i))*(Lcum(i)-Lcum(i-1)) + ...
0.5*dy(i)/dL(i)*(Lcum(i)^2 - Lcum(i-1)^2);
end
c0=temp/L;
for n=1:N
tempa = dx(1)/dL(1)*(cos(2*pi*n*Lcum(1)/L)-1);
tempb = dx(1)/dL(1)*sin(2*pi*n*Lcum(1)/L);
tempc = dy(1)/dL(1)*(cos(2*pi*n*Lcum(1)/L)-1);
tempd = dy(1)/dL(1)*sin(2*pi*n*Lcum(1)/L);
for i=2:K
tempa = tempa + ...
dx(i)/dL(i)*(cos(2*pi*n*Lcum(i)/L)-cos(2*pi*n*Lcum(i-1)/L));
tempb = tempb + ...
dx(i)/dL(i)*(sin(2*pi*n*Lcum(i)/L)-sin(2*pi*n*Lcum(i-1)/L));
tempc = tempc + ...
dy(i)/dL(i)*(cos(2*pi*n*Lcum(i)/L)-cos(2*pi*n*Lcum(i-1)/L));
tempd = tempd + ...
dy(i)/dL(i)*(sin(2*pi*n*Lcum(i)/L)-sin(2*pi*n*Lcum(i-1)/L));
end
a(n) = tempa * L /(2*pi^2*n^2);
b(n) = tempb * L /(2*pi^2*n^2);
c(n) = tempc * L /(2*pi^2*n^2);
d(n) = tempd * L /(2*pi^2*n^2);
end
%draw the reconstruction
figure(1);
axis ij
tt=0:.02:1;
for t=tt
x = a0 + sum(a.*cos(2*pi*(1:N)*t) + b.*sin(2*pi*(1:N)*t));
y = c0 - 40 + sum(c.*cos(2*pi*(1:N)*t) + d.*sin(2*pi*(1:N)*t));
line(x,y,'markersize',5,'marker','x')
end
%57-y is needed because of inversion in image axis
% line(a0,c0,'markersize',5,'marker','x','markeredgecolor','red')
% line(sum(ypos)/K, sum(xpos)/K,...
% 'markersize',5,'marker','x','markeredgecolor','k')
%
Results
The shape that was animated was chosen to be a very simple approximation of a insect larvae turning. There are 4 reconstructions linked below:
- When the reconstruction is limited to the fundamental (n=1), we get an ellipse. The x in the center of the elipse is the program's estimate of the DC component (position).
- Adding in n=2 shows some bending.
- Adding n=3 shows some of the rectangular shape, since right angle information is carried in the 3rd harmonic.
- Adding n=4 and 5 shows fairly good reconstruction. The Fourier series terms diminish in amplitude as 1/n2, so 5 harmonics should get us within 4% or so of the actual curve.
References
- Ghosh, PK and Jain, PK (1993) An Algebra of Geometric Shapes, IEEE Computer Graphics and Applications, vol 13 pp 50-59, issue 5
- Rouben Rostamian, Equation of an Ellipse, http://mathforum.org/epigone/geom.puzzles/27/ce2iei180me9@forum.mathforum.com
- Tsow R, Land B (1994), Applications of transformed curves, Cornell Theory Center Technical Report 94-182, Cornell University,