  function [xf,yf] = ShowPBetter(x,y,nSmoothings,display)
% Polygon Smoothing
% x and y are column n-vectors with centroid zero
% nSmoothings is the number of smoothing steps taken.
% If display == 0, only the final polygon is displayed
% If display == 1, entire polygon sequence displayed.
% Uses axis equal off instead of axis tight off to ensure that circles
% are displayed as circles
% xf and yf are column n-vectors that define the vertics of the last
% polygon.

  close all
  figure
  set(gcf,'position',[100 100 1000 500])
  subplot(1,2,1)
  a = .6;
  plot([x;x(1)],[y;y(1)],'k',x,y,'or',[-a a],[a -a],'m',[-a a],[-a a],'m')
  axis equal off
  title('Original Polygon','FontWeight','bold','FontSize',14)
 
% Repeatedly smooth and display the polygon... 
  for k=1:nSmoothings
     subplot(1,2,2)
     [x,y] = Smooth(x,y);
     if display==1 || k == nSmoothings
        plot([x ; x(1)],[y ; y(1)],'k',x,y,'or',[-a a],[a -a],'m',[-a a],[-a a],'m')
        axis equal off
        title(sprintf('After %1d Smoothings',k),'FontWeight','bold','FontSize',14)
        shg
     end
  end
  xf = x;
  yf = y;

  function [xNew,yNew] = Smooth(x,y)
% x and y are column n-vectors that define the vertices of a polygon P. 
% xNew and yNew are column n-vectors that define a polygon Q
% whose edges are the midpoints of P's edges.

  n = length(x); xNew = zeros(n,1); yNew = zeros(n,1);
  for i=1:n-1
     xNew(i) = (x(i) + x(i+1))/2;
     yNew(i) = (y(i) + y(i+1))/2;
  end
  xNew(n) = (x(n)+x(1))/2;  xNew = xNew/norm(xNew);
  yNew(n) = (y(n)+y(1))/2;  yNew = yNew/norm(yNew);