octave_packages/geometry-1.5.0/geom2d/hexagonalGrid.m

   1 %% Copyright (c) 2011, INRA
   2 %% 2007-2011, David Legland <david.legland@grignon.inra.fr>
   3 %% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 %%
   5 %% All rights reserved.
   6 %% (simplified BSD License)
   7 %%
   8 %% Redistribution and use in source and binary forms, with or without
   9 %% modification, are permitted provided that the following conditions are met:
  10 %%
  11 %% 1. Redistributions of source code must retain the above copyright notice, this
  12 %%    list of conditions and the following disclaimer.
  13 %%
  14 %% 2. Redistributions in binary form must reproduce the above copyright notice,
  15 %%    this list of conditions and the following disclaimer in the documentation
  16 %%    and/or other materials provided with the distribution.
  17 %%
  18 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 %% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22 %% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 %% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 %% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 %% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 %% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 %% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 %% POSSIBILITY OF SUCH DAMAGE.
  29 %%
  30 %% The views and conclusions contained in the software and documentation are
  31 %% those of the authors and should not be interpreted as representing official
  32 %% policies, either expressed or implied, of copyright holder.
  33
  34 %% -*- texinfo -*-
  35 %% @deftypefn {Function File} {@var{pts} = } hexagonalGrid (@var{bounds}, @var{origin}, @var{size})
  36 %% Generate hexagonal grid of points in the plane.
  37 %%
  38 %%   usage
  39 %%   PTS = hexagonalGrid(BOUNDS, ORIGIN, SIZE)
  40 %%   generate points, lying in the window defined by BOUNDS (=[xmin ymin
  41 %%   xmax ymax]), starting from origin with a constant step equal to size.
  42 %%   SIZE is constant and is equals to the length of the sides of each
  43 %%   hexagon.
  44 %%
  45 %%   TODO: add possibility to use rotated grid
  46 %% @end deftypefn
  47
  48 function varargout = hexagonalGrid(bounds, origin, size, varargin)
  49
  50   size = size(1);
  51   dx = 3*size;
  52   dy = size*sqrt(3);
  53
  54
  55
  56   % consider two square grids with different centers
  57   pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
  58   pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
  59   pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
  60   pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});
  61
  62   % gather points
  63   pts = [pts1;pts2;pts3;pts4];
  64
  65
  66
  67
  68   % eventually compute also edges, clipped by bounds
  69   % TODO : manage generation of edges
  70   if nargout>1
  71       edges = zeros([0 4]);
  72       x0 = origin(1);
  73       y0 = origin(2);
  74
  75       % find all x coordinate
  76       x1 = bounds(1) + mod(x0-bounds(1), dx);
  77       x2 = bounds(3) - mod(bounds(3)-x0, dx);
  78       lx = (x1:dx:x2)';
  79
  80       % horizontal edges : first find y's
  81       y1 = bounds(2) + mod(y0-bounds(2), dy);
  82       y2 = bounds(4) - mod(bounds(4)-y0, dy);
  83       ly = (y1:dy:y2)';
  84
  85       % number of points in each coord, and total number of points
  86       ny = length(ly);
  87       nx = length(lx);
  88
  89       if bounds(1)-x1+dx<size
  90           disp('intersect bounding box');
  91       end
  92
  93       if bounds(3)-x2<size
  94           disp('intersect 2');
  95           edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
  96           x2 = x2-dx;
  97           lx = (x1:dx:x2)';
  98           nx = length(lx);
  99       end
 100
 101       for i=1:length(ly)
 102           ind = (1:nx)';
 103           tmpEdges(ind, 1) = lx;
 104           tmpEdges(ind, 2) = ly(i);
 105           tmpEdges(ind, 3) = lx+size;
 106           tmpEdges(ind, 4) = ly(i);
 107           edges = [edges; tmpEdges];
 108       end
 109
 110   end
 111
 112   % process output arguments
 113   if nargout>0
 114       varargout{1} = pts;
 115
 116       if nargout>1
 117           varargout{2} = edges;
 118       end
 119   end
 120
 121 endfunction
 122