function ShowAveraging()
% Illustrates functions use to examine repeated polygon smoothing.
close all
n = 20; % The number of polygon edges
speedy = 1; % Setting this to zero slows the simulation
its = 0; % Keeps track of the number of averagings
[x,y,u,v,d,itMax] = GenEG(n); % Generate an example
Show(x,y,its,d,u,v) % Display the polygon.
pause(1)
close
[x,y] = AverageOnce(x,y); % Average once
its = its+1;
Show(x,y,its,d,u,v)
pause(1)
close
[x,y] = AverageOnce(x,y); % Average again
its = its+1;
Show(x,y,its,d,u,v)
pause(1)
[x,y] = Simulate(x,y,d,its,itMax,speedy,u,v); % Repeatedly average
% To suppress the display of the limiting ellipse, change the each
% Show(x,y,its,d,u,v) to Show(x,y,its,d) and change the last line
% to [x,y] = Simulate(x,y,d,its,itMax,speedy);
% For your convenience, the following functions are also available in
% UntangleGUI.m.
% It is recommended that you study how AverageOnce, Simulate, and Show
% work. GenEg involves more advanced math--do not spend time on how
% it works.
function [x,y] = Simulate(x0,y0,d,k,itMax,speedy,u,v)
% x0 and y0 are column n-vectors that define the vertices of a normalized
% polygon P.
% k is a positive integer that represents the number of averagings that
% have been performed on some original polygon to obtain P.
% itMax a positive integer >=k.
% x and y are column n-vectors that define the vertices of a normalized
% polygon that is the result of itMax averagings of the same original
% polygon.
% u and v are vectors with the property that plot(u,v) displays the
% limiting ellipse.
% speedy has the value of 1 if a fast simulation is desired, 0 otherwise.
% d is a real number with the property that the simulation takes place in
% the square with vertices (d,d), (-d,d), (-d,-d), (d,-d).
% Set up the time delay based on the value of speedy...
if speedy
tDelay = 0.01;
else
tDelay = .15;
end
% The simulation...
x = x0;
y = y0;
for it=k+1:itMax
% Average and display..
[x,y] = AverageOnce(x,y);
if nargin==6
% Do not include the limiting ellipse...
Show(x,y,it,d)
else
% Include the limiting ellipse...
Show(x,y,it,d,u,v)
end
pause(tDelay)
end
function Show(x,y,k,d,u,v)
% Displays the polygon whose vertices are defined by the column n-vectors
% x and y.
% k is an integer, typically the number of completed averaging steps.
% d is a positive real. The x and y ranges will be -d to +d.
% If nargin==6, then the u and v are vector of the limiting ellipse and it
% will be displayed.
% A call of the form Show(x,y,k,d) is also valid but does not display the
% limiting ellipse.
plot([x;x(1)],[y;y(1)],'r',x,y,'.r','Markersize',20)
axis([-d d -d d])
axis square off
title(sprintf('n = %2d Number of Averagings = %1d',length(x),k),'FontSize',18)
hold on
plot([d d -d -d d],[-d d d -d -d],'linewidth',10)
if nargin==6
plot(u,v,'k','linewidth',2')
plot([-d d],[-d d],'--k','linewidth',2)
plot([-d d],[d -d],'--k','linewidth',2)
end
hold off
function [xNew,yNew] = AverageOnce(x,y)
% x and y are column n-vectors that define the vertices of a normalized
% polygon P.
% xNew and yNew are column n-vectors that define the vertices of a normalized
% polygon Q obtained by connecting the midpoints of P's edges.
n = length(x);
xNew = (x + [x(2:n);x(1)])/2; xNew = xNew - mean(xNew); xNew = xNew/norm(xNew);
yNew = (y + [y(2:n);y(1)])/2; yNew = yNew - mean(yNew); yNew = yNew/norm(yNew);
function [x,y,u,v,d,itMax] = GenEG(n)
% n is an integer >= 4.
% x and y are column n-vectors that define the vertices of a normalized polygon.
% u and v are vectors with the property that plot(u,v) displays the
% limiting ellipse.
% d is a positive real with the property that the simulation takes
% place in the square with vertices (d,d), (-d,d), (-d,-d), (d,-d).
% itMax is a positive integer with the property that with that many
% iterations, the polygon has visiably converged to the limiting ellipse.
% x and y vectors define the normalized random polygon P...
x = rand(n,1); x = x-mean(x); x = x/norm(x);
y = rand(n,1); y = y-mean(y); y = y/norm(y);
% u and v vectors...
tau = (2*pi/n)*(0:n-1)';
c = sqrt(2/n)*cos(tau);
s = sqrt(2/n)*sin(tau);
cos_theta_u = (c'*x)/sqrt((c'*x)^2+(s'*x)^2);
sin_theta_u = (s'*x)/sqrt((c'*x)^2+(s'*x)^2);
cos_theta_v = (c'*y)/sqrt((c'*y)^2+(s'*y)^2);
sin_theta_v = (s'*y)/sqrt((c'*y)^2+(s'*y)^2);
R = sqrt(2/n)*[cos_theta_u sin_theta_u ; cos_theta_v sin_theta_v];
t = linspace(0,2*pi);
M = R*[cos(t);sin(t)];
u = M(1,:)';
v = M(2,:)';
% Simulation parameters...
d = 1.6*max(max(u),max(v));
itMax = ceil(3./(-log10( 1-(3/2)*(pi/n)^2 )));