randomPointInPolygon

RANDOMPOINTINPOLYGON Generate random point(s) in a polygon.

   POINT = randomPointInPolygon(POLY)
   POINTS = randomPointInPolygon(POLY, NPTS)
   Generate random point(s) within the specified polygon. If the number of
   points is not specified, only one point is generated.
   
   POINT = randomPointInPolygon(POLY, NPTS, QRS)
   Specifies a Quasi-random number generator QRS used to generate point.
   coordinates (requires the statistics toolbox).

   Example
     % generate an ellipse-like polygon, and fill it with points
     elli = [50 50 50 30 30];
     poly = ellipseToPolygon(elli, 200);
     pts = randomPointInPolygon(poly, 500);
     figure; hold on;
     drawPolygon(poly, 'k');
     drawPoint(pts, 'b.');
     axis equal; axis([0 100 0 100]);

     % use halton sequence to distribute points within the polygon
     qrs = haltonset(2, 'Skip', 1e3, 'Leap', 1e2);
     pts = randomPointInPolygon(poly, 500, qrs);
     figure; hold on;
     drawPolygon(poly, 'k');
     drawPoint(pts, 'b.');
     axis equal; axis([0 100 0 100]);

   See also
     polygons2d, randomPointInBox, drawPolygon

0001 function points = randomPointInPolygon(poly, varargin)
0002 %RANDOMPOINTINPOLYGON Generate random point(s) in a polygon.
0003 %
0004 %   POINT = randomPointInPolygon(POLY)
0005 %   POINTS = randomPointInPolygon(POLY, NPTS)
0006 %   Generate random point(s) within the specified polygon. If the number of
0007 %   points is not specified, only one point is generated.
0008 %
0009 %   POINT = randomPointInPolygon(POLY, NPTS, QRS)
0010 %   Specifies a Quasi-random number generator QRS used to generate point.
0011 %   coordinates (requires the statistics toolbox).
0012 %
0013 %   Example
0014 %     % generate an ellipse-like polygon, and fill it with points
0015 %     elli = [50 50 50 30 30];
0016 %     poly = ellipseToPolygon(elli, 200);
0017 %     pts = randomPointInPolygon(poly, 500);
0018 %     figure; hold on;
0019 %     drawPolygon(poly, 'k');
0020 %     drawPoint(pts, 'b.');
0021 %     axis equal; axis([0 100 0 100]);
0022 %
0023 %     % use halton sequence to distribute points within the polygon
0024 %     qrs = haltonset(2, 'Skip', 1e3, 'Leap', 1e2);
0025 %     pts = randomPointInPolygon(poly, 500, qrs);
0026 %     figure; hold on;
0027 %     drawPolygon(poly, 'k');
0028 %     drawPoint(pts, 'b.');
0029 %     axis equal; axis([0 100 0 100]);
0030 %
0031 %   See also
0032 %     polygons2d, randomPointInBox, drawPolygon
0033 %
0034  
0035 % ------
0036 % Author: David Legland
0037 % e-mail: david.legland@inra.fr
0038 % Created: 2017-08-28,    using Matlab 9.1.0.441655 (R2016b)
0039 % Copyright 2017 INRA - Cepia Software Platform.
0040 
0041 % determine number of points to generate
0042 nPts = 1;
0043 if nargin > 1
0044     nPts = varargin{1};
0045     varargin(1) = [];
0046 end
0047 
0048 % if additional input arg is provided, use it as quasi-random generator
0049 stream = [];
0050 if ~isempty(varargin)
0051     stream = qrandstream(varargin{1});
0052 end
0053 
0054 % polygon extreme coordinates
0055 box = polygonBounds(poly);
0056 xmin = box(1);  xmax = box(2);
0057 ymin = box(3);  ymax = box(4);
0058 
0059 % compute size of box
0060 boxSizeX = xmax - xmin;
0061 boxSizeY = ymax - ymin;
0062 
0063 % allocate memory for result
0064 points = zeros(nPts, 2);
0065 
0066 % contains indices of remaining points to process
0067 ind = (1:nPts)';
0068 
0069 % iterate until all points have been sampled within the polygon
0070 if isempty(stream)
0071     % use default random number generator
0072     while ~isempty(ind)
0073         NI = length(ind);
0074         x = rand(NI, 1) * boxSizeX + xmin;
0075         y = rand(NI, 1) * boxSizeY + ymin;
0076         points(ind, :) = [x y];
0077 
0078         ind = ind(~polygonContains(poly, points(ind, :)));
0079     end
0080 
0081 else
0082     % use specified quasi-random generator
0083     while ~isempty(ind)
0084         NI = length(ind);
0085         pts = qrand(stream, NI);
0086         x = pts(:, 1) * boxSizeX + xmin;
0087         y = pts(:, 2) * boxSizeY + ymin;
0088         points(ind, :) = [x y];
0089 
0090         ind = ind(~polygonContains(poly, points(ind, :)));
0091     end
0092     
0093 end