function pts = rand_points_plot(N,bnd);

% RAND_POINTS_PLOT generates N random points in the polygon defined
% by bnd, and plots the results.
%
%INPUTS:      N = number of points to generate
%           bnd = boundary file in 'polyform.m' format 
%                 [1,n
%                  x1,y1
%                  x2,y2
%                  .
%                  .
%                  x1,y1]; last point same as first.

% GOOD EXAMPLE:    P = [1 7
%                       -1 -1
%                        1  0
%                       -1  1
%                        0  2
%                        2  1
%                        2 -1
%                        0 -2
%                       -1 -1];
 
Z = polyform(bnd);

P = bnd(Z(1):Z(2),:);

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

poly_plot(bnd,pts);