Mercurial > forge

--- a/main/geometry/doc/NEWS	Fri Oct 21 18:58:15 2011 +0000
+++ b/main/geometry/doc/NEWS	Fri Oct 21 19:10:17 2011 +0000
@@ -131,6 +131,12 @@
  inertiaEllipse.m
  changelog.txt
  readme.txt
-
+ hexagonalGrid.m
+ squareGrid.m
+ triangleGrid.m
+ intersectCircles.m
+ intersectEdges.m
+ intersectLineCircle.m
+
 ===============================================================================
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/hexagonalGrid.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,122 @@
+%% Copyright (c) 2011, INRA
+%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{pts} = } hexagonalGrid (@var{bounds}, @var{origin}, @var{size})
+%% 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
+%% @end deftypefn
+
+function varargout = hexagonalGrid(bounds, origin, size, varargin)
+
+  size = size(1);
+  dx = 3*size;
+  dy = size*sqrt(3);
+
+
+
+  % consider two square grids with different centers
+  pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
+  pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
+  pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
+  pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});
+
+  % gather points
+  pts = [pts1;pts2;pts3;pts4];
+
+
+
+
+  % eventually compute also edges, clipped by bounds
+  % TODO : manage generation of edges
+  if nargout>1
+      edges = zeros([0 4]);
+      x0 = origin(1);
+      y0 = origin(2);
+
+      % find all x coordinate
+      x1 = bounds(1) + mod(x0-bounds(1), dx);
+      x2 = bounds(3) - mod(bounds(3)-x0, dx);
+      lx = (x1:dx:x2)';
+
+      % horizontal edges : first find y's
+      y1 = bounds(2) + mod(y0-bounds(2), dy);
+      y2 = bounds(4) - mod(bounds(4)-y0, dy);
+      ly = (y1:dy:y2)';
+
+      % number of points in each coord, and total number of points
+      ny = length(ly);
+      nx = length(lx);
+
+      if bounds(1)-x1+dx<size
+          disp('intersect bounding box');
+      end
+
+      if bounds(3)-x2<size
+          disp('intersect 2');
+          edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
+          x2 = x2-dx;
+          lx = (x1:dx:x2)';
+          nx = length(lx);
+      end
+
+      for i=1:length(ly)
+          ind = (1:nx)';
+          tmpEdges(ind, 1) = lx;
+          tmpEdges(ind, 2) = ly(i);
+          tmpEdges(ind, 3) = lx+size;
+          tmpEdges(ind, 4) = ly(i);
+          edges = [edges; tmpEdges];
+      end
+
+  end
+
+  % process output arguments
+  if nargout>0
+      varargout{1} = pts;
+
+      if nargout>1
+          varargout{2} = edges;
+      end
+  end
+
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/intersectCircles.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,129 @@
+%% Copyright (c) 2011, INRA
+%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{points} = } intersectCircles (@var{circle1}, @var{circle2})
+%% Intersection points of two circles.
+%%
+%%   POINTS = intersectCircles(CIRCLE1, CIRCLE2)
+%%   Computes the intersetion point of the two circles CIRCLE1 and CIRCLE1.
+%%   Both circles are given with format: [XC YC R], with (XC,YC) being the
+%%   coordinates of the center and R being the radius.
+%%   POINTS is a 2-by-2 array, containing coordinate of an intersection
+%%   point on each row.
+%%   In the case of tangent circles, the intersection is returned twice. It
+%%   can be simplified by using the 'unique' function.
+%%
+%%   Example
+%%     % intersection points of two distant circles
+%%     c1 = [0  0 10];
+%%     c2 = [10 0 10];
+%%     pts = intersectCircles(c1, c2)
+%%     pts =
+%%         5   -8.6603
+%%         5    8.6603
+%%
+%%     % intersection points of two tangent circles
+%%     c1 = [0  0 10];
+%%     c2 = [20 0 10];
+%%     pts = intersectCircles(c1, c2)
+%%     pts =
+%%         10    0
+%%         10    0
+%%     pts2 = unique(pts, 'rows')
+%%     pts2 =
+%%         10    0
+%%
+%%   References
+%%   http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
+%%   http://mathworld.wolfram.com/Circle-CircleIntersection.html
+%%
+%%   @seealso{circles2d, intersectLineCircle, radicalAxis}
+%% @end deftypefn
+
+function points = intersectCircles(circle1, circle2)
+
+  % adapt sizes of inputs
+  n1 = size(circle1, 1);
+  n2 = size(circle2, 1);
+  if n1 ~= n2
+      if n1 > 1 && n2 == 1
+          circle2 = repmat(circle2, n1, 1);
+      elseif n2 > 1 && n1 == 1
+          circle1 = repmat(circle1, n2, 1);
+      else
+          error('Both input should have same number of rows');
+      end
+  end
+
+  % extract center and radius of each circle
+  center1 = circle1(:, 1:2);
+  center2 = circle2(:, 1:2);
+  r1 = circle1(:,3);
+  r2 = circle2(:,3);
+
+  % allocate memory for result
+  nPoints = length(r1);
+  points = NaN * ones(2*nPoints, 2);
+
+  % distance between circle centers
+  d12 = distancePoints(center1, center2, 'diag');
+
+  % get indices of circle couples with intersections
+  inds = d12 >= abs(r1 - r2) & d12 <= (r1 + r2);
+
+  if sum(inds) == 0
+      return;
+  end
+
+  % angle of line from center1 to center2
+  angle = angle2Points(center1(inds,:), center2(inds,:));
+
+  % position of intermediate point, located at the intersection of the
+  % radical axis with the line joining circle centers
+  d1m  = d12(inds) / 2 + (r1(inds).^2 - r2(inds).^2) ./ (2 * d12(inds));
+  tmp = polarPoint(center1(inds, :), d1m, angle);
+
+  % distance between intermediate point and each intersection point
+  h   = sqrt(r1(inds).^2 - d1m.^2);
+
+  % indices of valid intersections
+  inds2 = find(inds)*2;
+  inds1 = inds2 - 1;
+
+  % create intersection points
+  points(inds1, :) = polarPoint(tmp, h, angle - pi/2);
+  points(inds2, :) = polarPoint(tmp, h, angle + pi/2);
+
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/intersectEdges.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,175 @@
+%% Copyright (c) 2011, INRA
+%% 2003-2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{point} = } intersectEdges (@var{edge1}, @var{edge2})
+%% Return all intersections between two set of edges
+%%
+%%   P = intersectEdges(E1, E2);
+%%   returns the intersection point of lines L1 and L2. E1 and E2 are 1-by-4
+%%   arrays, containing parametric representation of each edge (in the form
+%%   [x1 y1 x2 y2], see 'createEdge' for details).
+%%
+%%   In case of colinear edges, returns [Inf Inf].
+%%   In case of parallel but not colinear edges, returns [NaN NaN].
+%%
+%%   If each input is [N*4] array, the result is a [N*2] array containing
+%%   intersections of each couple of edges.
+%%   If one of the input has N rows and the other 1 row, the result is a
+%%   [N*2] array.
+%%
+%%   @seealso{edges2d, intersectLines}
+%% @end deftypefn
+
+function point = intersectEdges(edge1, edge2)
+  %% Initialisations
+
+  % ensure input arrays are same size
+  N1  = size(edge1, 1);
+  N2  = size(edge2, 1);
+
+  % ensure input have same size
+  if N1~=N2
+      if N1==1
+          edge1 = repmat(edge1, [N2 1]);
+          N1 = N2;
+      elseif N2==1
+          edge2 = repmat(edge2, [N1 1]);
+      end
+  end
+
+  % tolerance for precision
+  tol = 1e-14;
+
+  % initialize result array
+  x0  = zeros(N1, 1);
+  y0  = zeros(N1, 1);
+
+
+  %% Detect parallel and colinear cases
+
+  % indices of parallel edges
+  %par = abs(dx1.*dy2 - dx1.*dy2)<tol;
+  par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2));
+
+  % indices of colinear edges
+  % equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps
+  col = abs(  (edge2(:,1)-edge1(:,1)).*(edge1(:,4)-edge1(:,2)) - ...
+              (edge2(:,2)-edge1(:,2)).*(edge1(:,3)-edge1(:,1)) )<tol & par;
+
+  % Parallel edges have no intersection -> return [NaN NaN]
+  x0(par & ~col) = NaN;
+  y0(par & ~col) = NaN;
+
+
+  %% Process colinear edges
+
+  % colinear edges may have 0, 1 or infinite intersection
+  % Discrimnation based on position of edge2 vertices on edge1
+  if sum(col)>0
+      % array for storing results of colinear edges
+      resCol = Inf*ones(size(col));
+
+      % compute position of edge2 vertices wrt edge1
+      t1 = edgePosition(edge2(col, 1:2), edge1(col, :));
+      t2 = edgePosition(edge2(col, 3:4), edge1(col, :));
+
+      % control location of vertices: we want t1<t2
+      if t1>t2
+          tmp = t1;
+          t1  = t2;
+          t2  = tmp;
+      end
+
+      % edge totally before first vertex or totally after last vertex
+      resCol(col(t2<-eps))  = NaN;
+      resCol(col(t1>1+eps)) = NaN;
+
+      % set up result into point coordinate
+      x0(col) = resCol;
+      y0(col) = resCol;
+
+      % touches on first point of first edge
+      touch = col(abs(t2)<1e-14);
+      x0(touch) = edge1(touch, 1);
+      y0(touch) = edge1(touch, 2);
+
+      % touches on second point of first edge
+      touch = col(abs(t1-1)<1e-14);
+      x0(touch) = edge1(touch, 3);
+      y0(touch) = edge1(touch, 4);
+  end
+
+
+  %% Process non parallel cases
+
+  % process intersecting edges whose interecting lines intersect
+  i = find(~par);
+
+  % use a test to avoid process empty arrays
+  if sum(i)>0
+      % extract base parameters of supporting lines for non-parallel edges
+      x1  = edge1(i,1);
+      y1  = edge1(i,2);
+      dx1 = edge1(i,3)-x1;
+      dy1 = edge1(i,4)-y1;
+      x2  = edge2(i,1);
+      y2  = edge2(i,2);
+      dx2 = edge2(i,3)-x2;
+      dy2 = edge2(i,4)-y2;
+
+      % compute intersection points of supporting lines
+      delta = (dx2.*dy1-dx1.*dy2);
+      x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta;
+      y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta;
+
+      % compute position of intersection points on each edge
+      % t1 is position on edge1, t2 is position on edge2
+      t1  = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1);
+      t2  = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2);
+
+      % check position of points on edges.
+      % it should be comprised between 0 and 1 for both t1 and t2.
+      % if not, the edges do not intersect
+      out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol;
+      x0(i(out)) = NaN;
+      y0(i(out)) = NaN;
+  end
+
+
+  %% format output arguments
+
+  point = [x0 y0];
+
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/intersectLineCircle.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,106 @@
+%% Copyright (c) 2011, INRA
+%% 2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{points} = } intersectLineCircle (@var{line}, @var{circle})
+%% Intersection point(s) of a line and a circle
+%%
+%%   INTERS = intersectLineCircle(LINE, CIRCLE);
+%%   Returns a 2-by-2 array, containing on each row the coordinates of an
+%%   intersection point. If the line and circle do not intersect, the result
+%%   is filled with NaN.
+%%
+%%   Example
+%%   % base point
+%%   center = [10 0];
+%%   % create vertical line
+%%   l1 = [center 0 1];
+%%   % circle
+%%   c1 = [center 5];
+%%   pts = intersectLineCircle(l1, c1)
+%%   pts =
+%%       10   -5
+%%       10    5
+%%   % draw the result
+%%   figure; clf; hold on;
+%%   axis([0 20 -10 10]);
+%%   drawLine(l1);
+%%   drawCircle(c1);
+%%   drawPoint(pts, 'rx');
+%%   axis equal;
+%%
+%%   @seealso{lines2d, circles2d, intersectLines, intersectCircles}
+%% @end deftypefn
+
+function points = intersectLineCircle(line, circle)
+
+  % local precision
+  eps = 1e-14;
+
+  % center parameters
+  center = circle(:, 1:2);
+  radius = circle(:, 3);
+
+  % line parameters
+  dp = line(:, 1:2) - center;
+  vl = line(:, 3:4);
+
+  % coefficient of second order equation
+  a = sum(line(:, 3:4).^2, 2);
+  b = 2*sum(dp .* vl, 2);
+  c =  sum(dp.^2, 2) - radius.^2;
+
+  % discriminant
+  delta = b .^ 2 - 4 * a .* c;
+
+  if delta > eps
+      % find two roots of second order equation
+      u1 = (-b - sqrt(delta)) / 2 ./ a;
+      u2 = (-b + sqrt(delta)) / 2 ./ a;
+
+      % convert into 2D coordinate
+      points = [line(1:2) + u1 * line(3:4) ; line(1:2) + u2 * line(3:4)];
+
+  elseif abs(delta) < eps
+      % find unique root, and convert to 2D coord.
+      u = -b / 2 ./ a;
+      points = line(1:2) + u*line(3:4);
+
+  else
+      % fill with NaN
+      points = NaN * ones(2, 2);
+      return;
+  end
+
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/squareGrid.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,82 @@
+%% Copyright (c) 2011, INRA
+%% 2005-2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{pts} = } squaregrid (@var{bounds}, @var{origin}, @var{size})
+%% Generate equally spaces points in plane.
+%%
+%%   usage
+%%   PTS = squareGrid(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.
+%%
+%%   Example
+%%   PTS = squareGrid([0 0 10 10], [3 3], [4 2])
+%%   will return points :
+%%   [3 1;7 1;3 3;7 3;3 5;7 5;3 7;7 7;3 9;7 9];
+%%
+%%   TODO: add possibility to use rotated grid
+%%
+%% @end deftypefn
+
+function varargout = squareGrid(bounds, origin, size)
+
+  % find all x coordinate
+  x1 = bounds(1) + mod(origin(1)-bounds(1), size(1));
+  x2 = bounds(3) - mod(bounds(3)-origin(1), size(1));
+  lx = (x1:size(1):x2)';
+
+  % find all y coordinate
+  y1 = bounds(2) + mod(origin(2)-bounds(2), size(2));
+  y2 = bounds(4) - mod(bounds(4)-origin(2), size(2));
+  ly = (y1:size(2):y2)';
+
+  % number of points in each coord, and total number of points
+  ny = length(ly);
+  nx = length(lx);
+  np = nx*ny;
+
+  % create points
+  pts = zeros(np, 2);
+  for i=1:ny
+      pts( (1:nx)'+(i-1)*nx, 1) = lx;
+      pts( (1:nx)'+(i-1)*nx, 2) = ly(i);
+  end
+
+  % process output
+  if nargout>0
+      varargout{1} = pts;
+  end
+
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/geometry/geom2d/inst/triangleGrid.m	Fri Oct 21 19:10:17 2011 +0000
@@ -0,0 +1,69 @@
+%% Copyright (c) 2011, INRA
+%% 2005-2011, David Legland <david.legland@grignon.inra.fr>
+%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%% All rights reserved.
+%% (simplified BSD License)
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1. Redistributions of source code must retain the above copyright notice, this
+%%    list of conditions and the following disclaimer.
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%    this list of conditions and the following disclaimer in the documentation
+%%    and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+%% POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of copyright holder.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{pts} = } triangleGrid (@var{bounds}, @var{origin}, @var{size})
+%% Generate triangular grid of points in the plane.
+%%
+%%   usage
+%%   PTS = triangleGrid(BOUNDS, ORIGIN, SIZE)
+%%   generate points, lying in the window defined by BOUNDS, given in form
+%%   [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
+%%   triangles.
+%%
+%%   TODO: add possibility to use rotated grid
+%%
+%% @end deftypefn
+
+function varargout = triangleGrid(bounds, origin, size, varargin)
+
+  dx = size(1);
+  dy = size(1)*sqrt(3);
+
+  % consider two square grids with different centers
+  pts1 = squareGrid(bounds, origin, [dx dy], varargin{:});
+  pts2 = squareGrid(bounds, origin + [dx dy]/2, [dx dy], varargin{:});
+
+  % gather points
+  pts = [pts1;pts2];
+
+
+  % process output
+  if nargout>0
+      varargout{1} = pts;
+  end
+
+endfunction
+
--- a/main/geometry/matGeom_raw/geom2d/hexagonalGrid.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,122 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function varargout = hexagonalGrid(bounds, origin, size, varargin)
-%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
-%
-%   ---------
-%
-%   author : David Legland
-%   INRA - TPV URPOI - BIA IMASTE
-%   created the 06/08/2005.
-%
-size = size(1);
-dx = 3*size;
-dy = size*sqrt(3);
-
-
-
-% consider two square grids with different centers
-pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
-pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
-pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
-pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});
-
-% gather points
-pts = [pts1;pts2;pts3;pts4];
-
-
-
-
-% eventually compute also edges, clipped by bounds
-% TODO : manage generation of edges
-if nargout>1
-    edges = zeros([0 4]);
-    x0 = origin(1);
-    y0 = origin(2);
-
-    % find all x coordinate
-    x1 = bounds(1) + mod(x0-bounds(1), dx);
-    x2 = bounds(3) - mod(bounds(3)-x0, dx);
-    lx = (x1:dx:x2)';
-
-    % horizontal edges : first find y's
-    y1 = bounds(2) + mod(y0-bounds(2), dy);
-    y2 = bounds(4) - mod(bounds(4)-y0, dy);
-    ly = (y1:dy:y2)';
-
-    % number of points in each coord, and total number of points
-    ny = length(ly);
-    nx = length(lx);
-
-    if bounds(1)-x1+dx<size
-        disp('intersect bounding box');
-    end
-
-    if bounds(3)-x2<size
-        disp('intersect 2');
-        edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
-        x2 = x2-dx;
-        lx = (x1:dx:x2)';
-        nx = length(lx);
-    end
-
-    for i=1:length(ly)
-        ind = (1:nx)';
-        tmpEdges(ind, 1) = lx;
-        tmpEdges(ind, 2) = ly(i);
-        tmpEdges(ind, 3) = lx+size;
-        tmpEdges(ind, 4) = ly(i);
-        edges = [edges; tmpEdges];
-    end
-
-end
-
-% process output arguments
-if nargout>0
-    varargout{1} = pts;
-
-    if nargout>1
-        varargout{2} = edges;
-    end
-end
--- a/main/geometry/matGeom_raw/geom2d/intersectCircles.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,135 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function points = intersectCircles(circle1, circle2)
-%INTERSECTCIRCLES Intersection points of two circles
-%
-%   POINTS = intersectCircles(CIRCLE1, CIRCLE2)
-%   Computes the intersetion point of the two circles CIRCLE1 and CIRCLE1.
-%   Both circles are given with format: [XC YC R], with (XC,YC) being the
-%   coordinates of the center and R being the radius.
-%   POINTS is a 2-by-2 array, containing coordinate of an intersection
-%   point on each row.
-%   In the case of tangent circles, the intersection is returned twice. It
-%   can be simplified by using the 'unique' function.
-%
-%   Example
-%     % intersection points of two distant circles
-%     c1 = [0  0 10];
-%     c2 = [10 0 10];
-%     pts = intersectCircles(c1, c2)
-%     pts =
-%         5   -8.6603
-%         5    8.6603
-%
-%     % intersection points of two tangent circles
-%     c1 = [0  0 10];
-%     c2 = [20 0 10];
-%     pts = intersectCircles(c1, c2)
-%     pts =
-%         10    0
-%         10    0
-%     pts2 = unique(pts, 'rows')
-%     pts2 =
-%         10    0
-%
-%   References
-%   http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
-%   http://mathworld.wolfram.com/Circle-CircleIntersection.html
-%
-%   See also
-%   circles2d, intersectLineCircle, radicalAxis
-%
-%
-% ------
-% Author: David Legland
-% e-mail: david.legland@grignon.inra.fr
-% Created: 2011-01-20,    using Matlab 7.9.0.529 (R2009b)
-% Copyright 2011 INRA - Cepia Software Platform.
-
-%   HISTORY
-%   2011-06-06 add support for multiple inputs
-
-
-% adapt sizes of inputs
-n1 = size(circle1, 1);
-n2 = size(circle2, 1);
-if n1 ~= n2
-    if n1 > 1 && n2 == 1
-        circle2 = repmat(circle2, n1, 1);
-    elseif n2 > 1 && n1 == 1
-        circle1 = repmat(circle1, n2, 1);
-    else
-        error('Both input should have same number of rows');
-    end
-end
-
-% extract center and radius of each circle
-center1 = circle1(:, 1:2);
-center2 = circle2(:, 1:2);
-r1 = circle1(:,3);
-r2 = circle2(:,3);
-
-% allocate memory for result
-nPoints = length(r1);
-points = NaN * ones(2*nPoints, 2);
-
-% distance between circle centers
-d12 = distancePoints(center1, center2, 'diag');
-
-% get indices of circle couples with intersections
-inds = d12 >= abs(r1 - r2) & d12 <= (r1 + r2);
-
-if sum(inds) == 0
-    return;
-end
-
-% angle of line from center1 to center2
-angle = angle2Points(center1(inds,:), center2(inds,:));
-
-% position of intermediate point, located at the intersection of the
-% radical axis with the line joining circle centers
-d1m  = d12(inds) / 2 + (r1(inds).^2 - r2(inds).^2) ./ (2 * d12(inds));
-tmp = polarPoint(center1(inds, :), d1m, angle);
-
-% distance between intermediate point and each intersection point
-h   = sqrt(r1(inds).^2 - d1m.^2);
-
-% indices of valid intersections
-inds2 = find(inds)*2;
-inds1 = inds2 - 1;
-
-% create intersection points
-points(inds1, :) = polarPoint(tmp, h, angle - pi/2);
-points(inds2, :) = polarPoint(tmp, h, angle + pi/2);
--- a/main/geometry/matGeom_raw/geom2d/intersectEdges.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,185 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function point = intersectEdges(edge1, edge2)
-%INTERSECTEDGES Return all intersections between two set of edges
-%
-%   P = intersectEdges(E1, E2);
-%   returns the intersection point of lines L1 and L2. E1 and E2 are 1-by-4
-%   arrays, containing parametric representation of each edge (in the form
-%   [x1 y1 x2 y2], see 'createEdge' for details).
-%
-%   In case of colinear edges, returns [Inf Inf].
-%   In case of parallel but not colinear edges, returns [NaN NaN].
-%
-%   If each input is [N*4] array, the result is a [N*2] array containing
-%   intersections of each couple of edges.
-%   If one of the input has N rows and the other 1 row, the result is a
-%   [N*2] array.
-%
-%   See also:
-%   edges2d, intersectLines
-%
-%   ---------
-%   author : David Legland
-%   INRA - TPV URPOI - BIA IMASTE
-%   created the 31/10/2003.
-%
-
-%   HISTORY
-%   19/02/2004 add support for multiple edges.
-%   15/08/2005 rewrite a lot, and create unit tests
-%   08/03/2007 update doc
-%   21/09/2009 fix bug for edge arrays (thanks to Miquel Cubells)
-%   03/08/2010 fix another bug for edge arrays (thanks to Reto Zingg)
-
-%% Initialisations
-
-% ensure input arrays are same size
-N1  = size(edge1, 1);
-N2  = size(edge2, 1);
-
-% ensure input have same size
-if N1~=N2
-    if N1==1
-        edge1 = repmat(edge1, [N2 1]);
-        N1 = N2;
-    elseif N2==1
-        edge2 = repmat(edge2, [N1 1]);
-    end
-end
-
-% tolerance for precision
-tol = 1e-14;
-
-% initialize result array
-x0  = zeros(N1, 1);
-y0  = zeros(N1, 1);
-
-
-%% Detect parallel and colinear cases
-
-% indices of parallel edges
-%par = abs(dx1.*dy2 - dx1.*dy2)<tol;
-par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2));
-
-% indices of colinear edges
-% equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps
-col = abs(  (edge2(:,1)-edge1(:,1)).*(edge1(:,4)-edge1(:,2)) - ...
-            (edge2(:,2)-edge1(:,2)).*(edge1(:,3)-edge1(:,1)) )<tol & par;
-
-% Parallel edges have no intersection -> return [NaN NaN]
-x0(par & ~col) = NaN;
-y0(par & ~col) = NaN;
-
-
-%% Process colinear edges
-
-% colinear edges may have 0, 1 or infinite intersection
-% Discrimnation based on position of edge2 vertices on edge1
-if sum(col)>0
-    % array for storing results of colinear edges
-    resCol = Inf*ones(size(col));
-
-    % compute position of edge2 vertices wrt edge1
-    t1 = edgePosition(edge2(col, 1:2), edge1(col, :));
-    t2 = edgePosition(edge2(col, 3:4), edge1(col, :));
-
-    % control location of vertices: we want t1<t2
-    if t1>t2
-        tmp = t1;
-        t1  = t2;
-        t2  = tmp;
-    end
-
-    % edge totally before first vertex or totally after last vertex
-    resCol(col(t2<-eps))  = NaN;
-    resCol(col(t1>1+eps)) = NaN;
-
-    % set up result into point coordinate
-    x0(col) = resCol;
-    y0(col) = resCol;
-
-    % touches on first point of first edge
-    touch = col(abs(t2)<1e-14);
-    x0(touch) = edge1(touch, 1);
-    y0(touch) = edge1(touch, 2);
-
-    % touches on second point of first edge
-    touch = col(abs(t1-1)<1e-14);
-    x0(touch) = edge1(touch, 3);
-    y0(touch) = edge1(touch, 4);
-end
-
-
-%% Process non parallel cases
-
-% process intersecting edges whose interecting lines intersect
-i = find(~par);
-
-% use a test to avoid process empty arrays
-if sum(i)>0
-    % extract base parameters of supporting lines for non-parallel edges
-    x1  = edge1(i,1);
-    y1  = edge1(i,2);
-    dx1 = edge1(i,3)-x1;
-    dy1 = edge1(i,4)-y1;
-    x2  = edge2(i,1);
-    y2  = edge2(i,2);
-    dx2 = edge2(i,3)-x2;
-    dy2 = edge2(i,4)-y2;
-
-    % compute intersection points of supporting lines
-    delta = (dx2.*dy1-dx1.*dy2);
-    x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta;
-    y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta;
-
-    % compute position of intersection points on each edge
-    % t1 is position on edge1, t2 is position on edge2
-    t1  = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1);
-    t2  = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2);
-
-    % check position of points on edges.
-    % it should be comprised between 0 and 1 for both t1 and t2.
-    % if not, the edges do not intersect
-    out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol;
-    x0(i(out)) = NaN;
-    y0(i(out)) = NaN;
-end
-
-
-%% format output arguments
-
-point = [x0 y0];
-
--- a/main/geometry/matGeom_raw/geom2d/intersectLineCircle.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,116 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function points = intersectLineCircle(line, circle)
-%INTERSECTLINECIRCLE Intersection point(s) of a line and a circle
-%
-%   INTERS = intersectLineCircle(LINE, CIRCLE);
-%   Returns a 2-by-2 array, containing on each row the coordinates of an
-%   intersection point. If the line and circle do not intersect, the result
-%   is filled with NaN.
-%
-%   Example
-%   % base point
-%   center = [10 0];
-%   % create vertical line
-%   l1 = [center 0 1];
-%   % circle
-%   c1 = [center 5];
-%   pts = intersectLineCircle(l1, c1)
-%   pts =
-%       10   -5
-%       10    5
-%   % draw the result
-%   figure; clf; hold on;
-%   axis([0 20 -10 10]);
-%   drawLine(l1);
-%   drawCircle(c1);
-%   drawPoint(pts, 'rx');
-%   axis equal;
-%
-%   See also
-%   lines2d, circles2d, intersectLines, intersectCircles
-%
-%   References
-%   http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
-%   http://mathworld.wolfram.com/Circle-LineIntersection.html
-%
-% ------
-% Author: David Legland
-% e-mail: david.legland@grignon.inra.fr
-% Created: 2011-01-14,    using Matlab 7.9.0.529 (R2009b)
-% Copyright 2011 INRA - Cepia Software Platform.
-
-%   HISTORY
-%   2011-06-06 fix bug in delta test
-
-
-% local precision
-eps = 1e-14;
-
-% center parameters
-center = circle(:, 1:2);
-radius = circle(:, 3);
-
-% line parameters
-dp = line(:, 1:2) - center;
-vl = line(:, 3:4);
-
-% coefficient of second order equation
-a = sum(line(:, 3:4).^2, 2);
-b = 2*sum(dp .* vl, 2);
-c =  sum(dp.^2, 2) - radius.^2;
-
-% discriminant
-delta = b .^ 2 - 4 * a .* c;
-
-if delta > eps
-    % find two roots of second order equation
-    u1 = (-b - sqrt(delta)) / 2 ./ a;
-    u2 = (-b + sqrt(delta)) / 2 ./ a;
-
-    % convert into 2D coordinate
-    points = [line(1:2) + u1 * line(3:4) ; line(1:2) + u2 * line(3:4)];
-
-elseif abs(delta) < eps
-    % find unique root, and convert to 2D coord.
-    u = -b / 2 ./ a;
-    points = line(1:2) + u*line(3:4);
-
-else
-    % fill with NaN
-    points = NaN * ones(2, 2);
-    return;
-end
-
--- a/main/geometry/matGeom_raw/geom2d/squareGrid.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function varargout = squareGrid(bounds, origin, size)
-%SQUAREGRID Generate equally spaces points in plane.
-%
-%   usage
-%   PTS = squareGrid(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.
-%
-%   Example
-%   PTS = squareGrid([0 0 10 10], [3 3], [4 2])
-%   will return points :
-%   [3 1;7 1;3 3;7 3;3 5;7 5;3 7;7 7;3 9;7 9];
-%
-%
-%
-%   TODO: add possibility to use rotated grid
-%
-%   ---------
-%
-%   author : David Legland
-%   INRA - TPV URPOI - BIA IMASTE
-%   created the 06/08/2005.
-%
-
-% find all x coordinate
-x1 = bounds(1) + mod(origin(1)-bounds(1), size(1));
-x2 = bounds(3) - mod(bounds(3)-origin(1), size(1));
-lx = (x1:size(1):x2)';
-
-% find all y coordinate
-y1 = bounds(2) + mod(origin(2)-bounds(2), size(2));
-y2 = bounds(4) - mod(bounds(4)-origin(2), size(2));
-ly = (y1:size(2):y2)';
-
-% number of points in each coord, and total number of points
-ny = length(ly);
-nx = length(lx);
-np = nx*ny;
-
-% create points
-pts = zeros(np, 2);
-for i=1:ny
-    pts( (1:nx)'+(i-1)*nx, 1) = lx;
-    pts( (1:nx)'+(i-1)*nx, 2) = ly(i);
-end
-
-% process output
-if nargout>0
-    varargout{1} = pts;
-end
--- a/main/geometry/matGeom_raw/geom2d/triangleGrid.m	Fri Oct 21 18:58:15 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%%    list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%%    this list of conditions and the following disclaimer in the documentation
-%%    and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-
-function varargout = triangleGrid(bounds, origin, size, varargin)
-%TRIANGLEGRID Generate triangular grid of points in the plane.
-%
-%   usage
-%   PTS = triangleGrid(BOUNDS, ORIGIN, SIZE)
-%   generate points, lying in the window defined by BOUNDS, given in form
-%   [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
-%   triangles.
-%
-%   TODO: add possibility to use rotated grid
-%
-%   ---------
-%
-%   author : David Legland
-%   INRA - TPV URPOI - BIA IMASTE
-%   created the 06/08/2005.
-%
-
-dx = size(1);
-dy = size(1)*sqrt(3);
-
-% consider two square grids with different centers
-pts1 = squareGrid(bounds, origin, [dx dy], varargin{:});
-pts2 = squareGrid(bounds, origin + [dx dy]/2, [dx dy], varargin{:});
-
-% gather points
-pts = [pts1;pts2];
-
-
-% process output
-if nargout>0
-    varargout{1} = pts;
-end