po_simple(polygon)

function simple = po_simple(polygon)
%po_simple: test of closed polygon (non) self-intersection
%   po_simple(p) is equal to true if the closed polygon p is not
%   self-intersecting.
%
%See also polygon.
%
%Polygon Toolbox by Eric Debreuve
%Last update: June 13, 2006

number_of_edges = size(polygon,2) - 1;

edges = diff(polygon, 1, 2);
lengths = sum(edges.^2);

for current_edge = 1:(number_of_edges - 2)
   current_vertices = polygon(:,[current_edge, current_edge+1]);
   bounding_box = sort(current_vertices, 2);

   if current_edge == 1
      against = 3:(number_of_edges - 1);
   else
      against = current_edge + 2:number_of_edges;
   end
   before = (polygon(1,[against, against(end)+1]) < bounding_box(1,1));
   after  = (polygon(1,[against, against(end)+1]) > bounding_box(1,2));
   out_of_bb = ((before(1:(end-1)) & before(2:end)) | (after(1:(end-1)) & after(2:end)));
   against(out_of_bb) = [];
   if isempty(against)
      continue
   end
   remaining_edges = unique([against, against+1]);
   before = (polygon(2,remaining_edges) < bounding_box(2,1));
   after  = (polygon(2,remaining_edges) > bounding_box(2,2));
   indices = cumsum([1, (diff(against)>1)+1]);
   out_of_bb = ((before(indices) & before(indices+1)) | (after(indices) & after(indices+1)));
   against(out_of_bb) = [];

   for candidate_edge = against
      v1 = current_vertices(:,1) - polygon(:,candidate_edge);
      v2 = current_vertices(:,2) - polygon(:,candidate_edge);

      d1 = v1(1) * edges(2,candidate_edge) - v1(2) * edges(1,candidate_edge);
      d2 = v2(1) * edges(2,candidate_edge) - v2(2) * edges(1,candidate_edge);

      if (d1 ~= 0) && (d1 * d2 <= 0)
         a = abs(d1) / (abs(d1) + abs(d2));
         intersection = current_vertices(:,1) + (a * edges(:,current_edge)) - polygon(:,candidate_edge);

         if (dot(edges(:,candidate_edge), intersection) > 0) && (lengths(candidate_edge) >= sum(intersection.^2))
            simple = false;
            return
         end
      end
   end
end

simple = true;