Convex polygon smoothing via DFT

indicative_output.md

Examples

Square

Freeform

main.m

% Data Acquisition
% Uncomment the following block to allow the user 
% to enter a polygon as a series of points on a 
% cartesian plane.
% [x, y] = ginput(8);

% Uncomment the following block to centre the 
% entered points around 0,0
% x = x - mean(x);
% y = y - mean(y);

% A shape we prepared earlier (It's a square).
% Comment the following block if you are trying to 
% process a user defined shape (from above).
x = [-1; 1; 1; -1];
y = [1; 1; -1; -1];


% Main Processing
K = 3;    % Recursive subdivision levels
fc = 6;   % "Filter Cut Off Frequency". That's the textbook name, in this context
          % it determines how many elementary ellipsoids to retain in reconstructing
          % the polygon and it thereore determines the "smoothing" factor.
n = 1;    % This is the "order" of the filter. In this context it determines how 
          % "strict" the filtering is.

% Due to the use of the Discrete Fourier Transform, this 
% periodicity is implied, but the first point is added here
% to the vectors to make it clearer visibly.
x = [x; x(1)]
y = [y; y(1)]

% Main Smoothing Step
% Simply produce a sequence of complex numbers.
% From this point onwards, the curve can be treated as a 
% vector of numbers.
s = x + i*y;

% Subdivide each piecewise segment to increase the spatial resolution 
% of the curve.
v = [s(1), subdiv(s(1), s(2), K)];
for u = 2:(length(s)-1)
    v = [v,[s(u), subdiv(s(u), s(u+1), K)]];
end;

% Construct the filter
% The simplest low-pass filter is really a Sinc filter that only retains 
% a smaller subset of the ellipsoids. In this context, using a smoother 
% transition band such as Butterworth's gives a better result.
h = (1./(1+(cumsum(ones(1, length(v)./2))./fc).^(2*n)));

% Apply the filter 
% (The filter has to be symmetric, in the same way the DFT is,
%  hence the mirroring of the `h` vector here)
q = ifft(fft(v).*[h, h(end:-1:1)]);

% Display everything
%
plot(s);        % The original polygon
hold on;        
plot(q,'.-');   % THe smoothed polygon

subdiv.m

function m = subdiv(m1,m2, N)
    % Plain simple recursive mid-point subdivision of a line segment. 
    % The line segment is defined by points m1,m2
    % 
    % m1: Starting point of a line segment as [x0,y0]
    % m2: Ending point of a line segment as [x1,y1]
    % N : Subdivision level
    %
    % Calculate the mid-point
    m = m1 + (m2-m1)/2.0;
    % Return the mid-point itself or the next level subdivision.
    if N>0
        p1 = subdiv(m1,m,N-1);
        p2 = subdiv(m,m2,N-1);
        m = [p1, m, p2];
end;