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.

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 
0012 % ------
0013 % Author: David Legland
0014 % E-mail: david.legland@inrae.fr
0015 % Created: 2005-08-06
0016 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0017 
0018 size = size(1);
0019 dx = 3*size;
0020 dy = size*sqrt(3);
0021 
0022 % consider two rectangular grids with shifted centers
0023 pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
0024 pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
0025 pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
0026 pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});
0027 
0028 % gather points
0029 pts = [pts1;pts2;pts3;pts4];
0030 
0031 
0032 % eventually compute also edges, clipped by bounds
0033 if nargout > 1
0034     edges = zeros([0 4]);
0035     x0 = origin(1);
0036     y0 = origin(2);
0037 
0038     % find all x coordinate
0039     x1 = bounds(1) + mod(x0-bounds(1), dx);
0040     x2 = bounds(3) - mod(bounds(3)-x0, dx);
0041     lx = (x1:dx:x2)';
0042 
0043     % horizontal edges: first find y's
0044     y1 = bounds(2) + mod(y0-bounds(2), dy);
0045     y2 = bounds(4) - mod(bounds(4)-y0, dy);
0046     ly = (y1:dy:y2)';
0047     
0048     % number of points in each coord, and total number of points
0049     ny = length(ly);
0050     nx = length(lx);
0051  
0052     if bounds(1)-x1+dx < size
0053         disp('intersect bounding box');
0054     end
0055     
0056     if bounds(3)-x2 < size
0057         disp('intersect 2');
0058         edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
0059         x2 = x2-dx;
0060         lx = (x1:dx:x2)';
0061         nx = length(lx);
0062     end
0063   
0064     for i = 1:length(ly)
0065         ind = (1:nx)';
0066         tmpEdges = zeros(length(ind), 4);
0067         tmpEdges(ind, 1) = lx;
0068         tmpEdges(ind, 2) = ly(i);
0069         tmpEdges(ind, 3) = lx+size;
0070         tmpEdges(ind, 4) = ly(i);
0071         edges = [edges; tmpEdges]; %#ok<AGROW>
0072     end
0073     
0074 end
0075 
0076 % process output arguments
0077 if nargout > 0
0078     varargout{1} = pts;
0079     
0080     if nargout > 1
0081         varargout{2} = edges;
0082     end
0083 end