Guest User

Staking Out the Integers Validator

a guest
Dec 10th, 2014
193
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function [As,PsPolys] = stakeareas( Ps )
  2.     n = size( Ps, 1 );
  3.     PIs = allperms( 1:n );
  4.     nP = size( PIs, 1 );
  5.     if anycolinear( Ps )
  6.         As = [];
  7.         return;
  8.     end
  9.  
  10.     As = [];
  11.     PsPolys_ = [];
  12.     for i = 1:nP
  13.         PsPoly = Ps(PIs(i,:),:);
  14.         if ~isjordan( PsPoly )
  15.             continue;
  16.         end
  17.         A = abs( area( PsPoly ) );
  18.         if ~ismember( A, As )
  19.             PsPolys_ = [PsPolys_; PsPoly]; %#ok<AGROW>
  20.             As = [As, A]; %#ok<AGROW>
  21.         end
  22.     end
  23.     As = sort( As );
  24.     if nargout >= 2
  25.         PsPolys = PsPolys_;
  26.     end
  27. end
  28.  
  29. function PIs = allperms( Is )
  30.     n = length( Is );
  31.     if n == 1
  32.         PIs = Is;
  33.     elseif n == 2
  34.         PIs = [Is; fliplr( Is )];
  35.     else
  36.         PIs = [];
  37.         for i = 1:n
  38.             i1 = Is(i);
  39.             IsRem = allperms( Is(setdiff( 1:n, i )) );
  40.             PIs = [PIs; i1*ones( size( IsRem, 1 ), 1 ), IsRem]; %#ok<AGROW>
  41.         end
  42.     end
  43. end
  44.  
  45. function b = isjordan( Ps )
  46.     n = size( Ps, 1 );
  47.     b = true;
  48.     for i = 1:n-2
  49.         for j = i+2:n
  50.             if intersect( Ps(i,:), Ps(i+1,:), Ps(j,:), Ps(mod( j, n )+1,:) )
  51.                 b = false;
  52.                 return;
  53.             end
  54.         end
  55.     end
  56.  
  57.     % Detects edges that "sneak through" a vertex. A complete hack, but guaranteed to work for integer coordinates.
  58.     for i = 1:n
  59.         P_nr = (Ps(i,:) + Ps(mod( i, n )+1,:)*63)/64;
  60.         P_fr = (Ps(mod( i+1, n )+1,:) + Ps(mod( i, n )+1,:)*63)/64;
  61.         for j = setdiff( 1:n, [i, mod( i, n )+1] )
  62.             if intersect( P_nr, P_fr, Ps(j,:), Ps(mod( j, n )+1,:) )
  63.                 b = false;
  64.                 return;
  65.             end
  66.         end
  67.     end
  68. end
  69.  
  70. function b = intersect( P11, P12, P21, P22 )
  71.     x1 = P12(1) - P11(1);
  72.     y1 = P12(2) - P11(2);
  73.     x2 = P22(1) - P21(1);
  74.     y2 = P22(2) - P21(2);
  75.  
  76.     b = ((y1*(P21(1) - P11(1))-(P21(2) - P11(2))*x1)* ...
  77.         (y1*(P22(1) - P11(1))-(P22(2) - P11(2))*x1) < 0) && ...
  78.             ((y2*(P11(1) - P21(1))-(P11(2) - P21(2))*x2)* ...
  79.         (y2*(P12(1) - P21(1))-(P12(2) - P21(2))*x2) < 0);
  80. end
  81.  
  82. function A = area( Ps )
  83.     n = size( Ps, 1 );
  84.     A = Ps(n,1)*Ps(1,2) - Ps(n,2)*Ps(1,1);
  85.     for i = 2:n
  86.         A = A + Ps(i-1,1)*Ps(i,2) - Ps(i-1,2)*Ps(i,1);
  87.     end
  88.     A = A/2;
  89. end
  90.  
  91. function b = anycolinear( Ps )
  92.     b = false;
  93.     n = size( Ps, 1 );
  94.     for i = 1:n-2
  95.         for j = i+1:n-1
  96.             for k = j+1:n
  97.                 if colinear( Ps(i,:), Ps(j,:), Ps(k,:) )
  98.                     b = true;
  99.                     return;
  100.                 end
  101.             end
  102.         end
  103.     end
  104. end
  105.  
  106. function b = colinear( P_m1, P, P_p1 )
  107.     y21 = P(2) - P_m1(2);
  108.     y13 = P_m1(2) - P_p1(2);
  109.     y32 = P_p1(2) - P(2);
  110.     L = y21*y21 + (P(1) - P_m1(1))*(P(1) - P_m1(1));
  111.     L2 = y13*y13 + (P_m1(1) - P_p1(1))*(P_m1(1) - P_p1(1));
  112.  
  113.     if L2 > L
  114.         L = L2;
  115.     end
  116.     L2 = y32*y32 + (P_p1(1) - P(1))*(P_p1(1) - P(1));
  117.     if L2 > L
  118.         L = L2;
  119.     end
  120.  
  121.     b = L < 1e-9 || (abs( P(1)*y13 + P_p1(1)*y21 + P_m1(1)*y32 ) < 1e-9*sqrt( L ));
  122. end
RAW Paste Data