Today I'm wrapping up my analysis of the results of my Cody problem on eliminating unnecessary polygon vertices.

Contents

A solution by Cris

Here is Cris' first correct solution.

function P = simplify_polygon_108897_luengo(P)
  N = size(P,1);
  if N<=2
     return
  end
  % assume last point is same as first one -- remove
  %    in production code I'd run through the whole list and remove
  %    duplicate points.
  P(end,:) = [];
  N = N-1
  % I use 0-based indexing to use MOD more easily, always index with ii+1!
  ii = 0; % The point under consideration
  first = NaN; % The first point we did not remove
  while N>=3 && ii~=first
    p1 = P(mod(ii-1,N)+1,:);
    p2 = P(ii+1,:);
    p3 = P(mod(ii+1,N)+1,:);
    v1 = p2-p1;
    v2 = p3-p2;
    disp([v1;v2])
    if ( v1(1)*v2(2) == v1(2)*v2(1) ) && any(v1.*v2>0)
       % the two vectors are in the same direction: delete the middle point
       P(ii+1,:) = [];
       N = N-1;
       ii = mod(ii,N); % in case this was the last point
    else
       % keep the point & move on
       if isnan(first)
           first = ii;
       end
       ii = mod(ii+1,N);
    end
  end
  % add first point to end again
  P(end+1,:) = P(1,:);
end

This code walks forward (in the while loop) through the polygon vertices. Each iteration of the while loop checks one vertex to see if it lies on the line segment between the previous vertex and the next one. If so, the code deletes that vertex.

When I first formulated this problem, I wasn't certain whether a fully vectorized solution was possible. I thought it might be necessary to remove vertices one at a time, as Cris does here.

Cris uses a few different techniques that I have often used in my own code.

Handling Edge Cases

I prefer to write code that handles edge cases naturally, without conditional logic. However, when edge cases can't be handled smoothly by the main algorithm without adding conditional logic, I often handle them with special-case code and an early return at the beginning.

N = size(P,1);
if N<=2
   return
end

Preprocess and Postprocess the Vertex List

For his algorithm, Cris didn't want the first polygon vertex to appear again at the end of the list. So he removed it at the beginning and added it back at the end.

P(end,:) = [];

     -SNIP-

% add first point to end again
P(end+1,:) = P(1,:);

Use the mod Function for Circular Indexing

Suppose you have an index into a vector, and you want the "next" element in the vector, where "next" means go back to the beginning if you are at the end. One could use conditional logic to achieve this:

if i < length(v)
    i = i + 1;
else
    i = 1;
end

Or you could use the mod (modulo) function. You just have to account for 1-based indexing in MATLAB.

i = mod(i - 1,length(v)) + 1;

The MATLAB function circshift is implemented using exactly this idea. Here's the key part of circshift:

% Loop through each dimension of the input matrix to
% calculate shifted indices
for k = 1:numDimsA
  m      = sizeA(k);
  idx{k} = mod((0:m-1)-p(k), m)+1;
end

% Perform the actual conversion by indexing
% into the input matrix
b = a(idx{:});

Using the MATLAB File Compare Tool

Have you ever used the MATLAB File Comparison Tool? It's can be very useful.

I was interested to see how the code changed when Cris submitted his second correct solution. Here's what the File Comparison Tool (available from the Current Folder Browser) showed:

I don't see any algorithm changes. It looks like Cris has started to tweak his solution to make it smaller (as measured by the Cody solution scoring system).

A Solution by Sven

The first correct solution was this one by Sven:

function Ps = simplify_polygon(P)
  % Handle degenerate cases!
  if size(P,1)<3
      Ps = P;
      return;
  end
  % Normalise the difference between successive points
  pdiff = diff(P([end-1 1:end],:)); % Append the 2nd-to-last
  npdiff = bsxfun(@rdivide, pdiff, sqrt(sum(pdiff.^2,2)));
  % Any successive normalised differences that don't change can be discarded
  keepIdxs = find(any(diff(npdiff),2));
  % Append the first/last points (the question states that P(1,:)==P(end,:))
  Ps = P(keepIdxs([1:end 1]),:);
end

Sven uses diff and bsxfun to get a vectorized solution. He also handles edge cases up front and postprocesses the vertex list.

A Solution by Richard

Richard got into the game early with Cris and Sven. Here is his first correct solution.

function Ps = simplify_polygon(P)
  if ismember(size(P, 1), [0 1])
    Ps = P;
    return
  end

   dP = diff(P);
   dP = bsxfun(@times, 1./hypot(dP(:, 1), dP(:, 2)), dP);
   idx = find(abs(1 - dot(dP, circshift(dP, 1), 2)) > 1e-10);
   Ps = P([idx; idx(1)], :);
 end

This solution is somewhat similar to Sven's. Note the use of hypot and circshift.

Tag Team Effort

Sven and Richard engaged in a tag team effort to make their solutions smaller and smaller. Their creative burst certainly contributed to one of the highest "average solutions per solver" seen on Cody:

Here's the one of the smallest solutions (at least, it's one of the smallest that doesn't use a Neural Networks toolbox function):

function P = simplify_polygon(P)
  try
  diag(sum(abs(diff(P)),2)) \ diff(P);
  P(any(ans - circshift(ans,1),2),:);
  P = vertcat(ans, ans(1,:));
  end

No, I wouldn't write production code like this. But it's clever, and understanding how it works can teach us useful things about MATLAB.

Output Argument = Input Argument

Did you know that you can use the same variable name in both the input argument list and the output argument list of a function? In other words, you can do this:

  function P = simplify_polygon(P)

The effect is to initialize the output argument value to the value of the input argument. It's equivalent to this:

  function P_out = simplify_polygon(P_in)
  P_out = P_in;

My favorite example of using this technique shows up as a solution to another of my Cody problems, swapping two values.

  function [b,a] = swap(a,b)

That's the entire function!

Richard and Sven use this technique in combination with a try block to avoid having to explicity write code to handle the special cases.

Automatic Output Variable, ans

When a function returns an output argument and you don't assign it to a variable, MATLAB automatically creates a variable called ans to hold the result. You have probably seen this behavior in the Command Window:

>> magic(3)

ans =

     8     1     6
     3     5     7
     4     9     2

>> sum(ans)

ans =

    15    15    15

But this behavior also happens when you are running functions. The behavior is exploited in these two lines:

   diag(sum(abs(diff(P)),2)) \ diff(P);
   P(any(ans - circshift(ans,1),2),:);

The call to diag produces an output argument that isn't stored explicitly in a variable. The next line can get that result by using the automatic variable ans.

The Importance of Formulating Your Problem Carefully

This Cody problem was originally inspired by the problem of postprocessing the output of the Image Processing Toolbox function bwboundaries in order to remove unneeded vertices. This function produces polygons containing horizontal and vertical line segments vertices whose coordinates are integers. My selection of test cases was heavily influenced by this use case. However, I posed the problem on Cody as a more general problem: "Eliminate unnecessary polygon vertices." As I pointed out last week, solving this problem robustly in general can be quite complicated because of floating-point arithmetic issues.

Richard and Sven were very aware of this:

Thanks to All the Solvers!

Although I have focused on a few particular solutions from three different people today, this problem sparked a very creative collection of solutions from many people. I wish could have spent more time analyzing and discussing all of the solutions.

Thanks for giving it a try, everyone, and I hope you enjoyed it!

Cody solvers, if you'd like to add a comment to this post, please tell us about a particularly interesting Cody problem that you think we should go look at.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14