function pts = pt_in_poly_plot(N,P);

% plots results for new pt_in_poly procedure. The idea here is to 
% count the number of crossings of first coords for points in P with
% that of the reference point, x.

%PT_IN_POLY_PLOT constructs N random points in poly and plots them.

%INPUTS:    N = number of points to generate
%           P = polygon in the form 
%                 [x1,y1
%                  x2,y2
%                  .
%                  .
%                  x1,y1]; last point same as first.

% GOOD EXAMPLE:    P = [-1 -1
%                        1  0
%                       -1  1
%                        0  2
%                        2  1
%                        2 -1
%                        0 -2
%                       -1 -1];

if size(P,2) ~= 2
    
    error('poly must be nx2 matrix');
    
end

%******************

%rand('seed',12345);

%*******************


%Form bounding box

v = min(P);
xmin = v(1);
ymin = v(2);

v = max(P);
xmax = v(1);
ymax = v(2);

pts = zeros(N,2);

poly = P;

P = P(1:end-1,:); %Remove repeat element

n = length(P);

u = ones(n,1);

X = P(:,1); %first coords of P
Y = P(:,2); %second coords of P

for p = 1:N    
   
    flag = 0;
    
    ctr = 1;
    
    while flag == 0
        
        x = xmin + rand(1,1)*(xmax - xmin);
    
        y = ymin + rand(1,1)*(ymax - ymin);        
             
        X0 = (X > x); 
        
        Y0 = (Y > y);
        
        Y1 = [Y0(end);Y0(1:end-1)];  %Shift Y0
        
        I = find(Y0 + Y1 == 1);  %Set of Y-crossing points for y
        
        k = length(I);
        
        if k == 0
            
            hit = 0;            
                                               
        else %must have some crossings
            
            rc = zeros(k,1); %indicator for right crossings
            
            for i = 1:k
                
                if I(i) == 1
                    
                    i0 = n;  %Previous index in P
                    
                else
                    
                    i0 = I(i)-1;  %Previous index in P
                    
                end
                
                x1 = X(I(i));
                
                x0 = X(i0);
                
                if x > max(x1,x0)
                    
                    rc(i) = 0; %no right crossing
                    
                elseif x <= min(x1,x0)
                    
                    rc(i) = 1; %must be a right crossing
                    
                else  %x between x1 and x0
                    
                    a = (x - x0)/(x1 - x0); %must have 0 <= a <= 1
                    
                    y1 = Y(I(i));
                
                    y0 = Y(i0);  %corresponding y-coords
                    
                    if (y1 - y0)/(x1 - x0) >= 0 %slope of crossing is nonnegative
                    
                        if y > a*y1 + (1-a)*y0
                            
                            rc(i) = 1; %no right crossing
                            
                        else
                            
                            rc(i) = 0;
                            
                        end
                        
                    elseif y > a*y1 + (1-a)*y0 %slope of crossing is negative
                        
                        rc(i) = 0; %no right crossing
                        
                    else
                        
                        rc(i) = 1;
                        
                    end                        
                                      
                end %end if                
                
            end %end for 1 = 1:k
            
            C = sum(rc); %number of right crossings
                
            if mod(C,2) == 1
                
                hit = 1;                   
                                   
            else
                
                hit = 0;
                
            end            
            
        end %end outer if
               
        if hit == 1
            
            pts(p,:) = [x,y];
            
            flag = 1;
            
        end
        
        ctr = ctr + 1;
        
        if ctr > 100
            
            error('appears to be infinite loop');
            
        end
        
    end %end while
    
end %end construction loop

%plot results

x = poly(:,1);

y = poly(:,2);

plot(x,y,'k','LineWidth',2); %plot polygon
    
hold on;

set(gca,'DataAspectRatio',[1 1 1]);

hold on;

x = pts(:,1);

y = pts(:,2);

plot(x,y,'r.','MarkerSize',15); %plot points