hexagonalGrid

HEXAGONALGRID Generate hexagonal grid of points in the plane.

   usage:
   PTS = hexagonalGrid(BOUNDS, ORIGIN, SIZE)
   generate points, lying in the window defined by BOUNDS (=[xmin ymin
   xmax ymax]), starting from origin with a constant step equal to size.
   SIZE is constant and is equals to the length of the sides of each
   hexagon. 

   TODO: add possibility to use rotated grid

0001 function varargout = hexagonalGrid(bounds, origin, size, varargin)
0002 %HEXAGONALGRID Generate hexagonal grid of points in the plane.
0003 %
0004 %   usage:
0005 %   PTS = hexagonalGrid(BOUNDS, ORIGIN, SIZE)
0006 %   generate points, lying in the window defined by BOUNDS (=[xmin ymin
0007 %   xmax ymax]), starting from origin with a constant step equal to size.
0008 %   SIZE is constant and is equals to the length of the sides of each
0009 %   hexagon.
0010 %
0011 %   TODO: add possibility to use rotated grid
0012 %
0013 
0014 %   ---------
0015 %
0016 %   author : David Legland
0017 %   INRA - TPV URPOI - BIA IMASTE
0018 %   created the 06/08/2005.
0019 %
0020 
0021 size = size(1);
0022 dx = 3*size;
0023 dy = size*sqrt(3);
0024 
0025 % consider two square grids with different centers
0026 pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
0027 pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
0028 pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
0029 pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});
0030 
0031 % gather points
0032 pts = [pts1;pts2;pts3;pts4];
0033 
0034 
0035 % eventually compute also edges, clipped by bounds
0036 % TODO : manage generation of edges
0037 if nargout > 1
0038     edges = zeros([0 4]);
0039     x0 = origin(1);
0040     y0 = origin(2);
0041 
0042     % find all x coordinate
0043     x1 = bounds(1) + mod(x0-bounds(1), dx);
0044     x2 = bounds(3) - mod(bounds(3)-x0, dx);
0045     lx = (x1:dx:x2)';
0046 
0047     % horizontal edges : first find y's
0048     y1 = bounds(2) + mod(y0-bounds(2), dy);
0049     y2 = bounds(4) - mod(bounds(4)-y0, dy);
0050     ly = (y1:dy:y2)';
0051     
0052     % number of points in each coord, and total number of points
0053     ny = length(ly);
0054     nx = length(lx);
0055  
0056     if bounds(1)-x1+dx < size
0057         disp('intersect bounding box');
0058     end
0059     
0060     if bounds(3)-x2 < size
0061         disp('intersect 2');
0062         edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
0063         x2 = x2-dx;
0064         lx = (x1:dx:x2)';
0065         nx = length(lx);
0066     end
0067   
0068     for i = 1:length(ly)
0069         ind = (1:nx)';
0070         tmpEdges = zeros(length(ind), 4);
0071         tmpEdges(ind, 1) = lx;
0072         tmpEdges(ind, 2) = ly(i);
0073         tmpEdges(ind, 3) = lx+size;
0074         tmpEdges(ind, 4) = ly(i);
0075         edges = [edges; tmpEdges]; %#ok<AGROW>
0076     end
0077     
0078 end
0079 
0080 % process output arguments
0081 if nargout > 0
0082     varargout{1} = pts;
0083     
0084     if nargout > 1
0085         varargout{2} = edges;
0086     end
0087 end