Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / geometry-1.5.0 / shape2d / curve2polyline.m

diff --git a/octave_packages/geometry-1.5.0/shape2d/curve2polyline.m b/octave_packages/geometry-1.5.0/shape2d/curve2polyline.m
new file mode 100644 (file)
index 0000000..389aec0
--- /dev/null
+++ b/octave_packages/geometry-1.5.0/shape2d/curve2polyline.m
@@ -0,0 +1,141 @@
+%% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%%
+%%    This program is free software: you can redistribute it and/or modify
+%%    it under the terms of the GNU General Public License as published by
+%%    the Free Software Foundation, either version 3 of the License, or
+%%    any later version.
+%%
+%%    This program is distributed in the hope that it will be useful,
+%%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%%    GNU General Public License for more details.
+%%
+%%    You should have received a copy of the GNU General Public License
+%%    along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{polyline} = } curve2polyline (@var{curve})
+%% @deftypefnx {Function File} {@var{polyline} = } curve2polyline (@dots{},@var{property},@var{value},@dots{})
+%% Adaptive sampling of a parametric curve.
+%%
+%% The @var{curve} is described as a 2-by-N matrix. Rows correspond to the
+%% polynomial (compatible with @code{polyval}) describing the respective component
+%% of the curve. The curve must be parametrized in the interval [0,1].
+%% The vertices of the polyline are accumulated in regions of the curve where
+%% the curvature is higher.
+%%
+%% @strong{Parameters}
+%% @table @samp
+%% @item 'Nmax'
+%% Maximum number of vertices. Not used.
+%% @item 'Tol'
+%% Tolerance for the error criteria. Default value @code{1e-4}.
+%% @item 'MaxIter'
+%% Maximum number of iterations. Default value @code{10}.
+%% @item 'Method'
+%% Not implemented.
+%% @end table
+%%
+%% @seealso{shape2polygon, curveval}
+%% @end deftypefn
+
+%% This function is based on the algorithm described in
+%% L. H. de Figueiredo (1993). "Adaptive Sampling of Parametric Curves". Graphic Gems III.
+%% I had to remove the recursion so this version could be improved.
+%% Thursday, April 12 2012 -- JuanPi
+
+function [polyline t bump]= curve2polyline (curve, varargin)
+%% TODO make tolerance relative to the "diameter" of the curve.
+
+  # --- Parse arguments --- #
+  parser = inputParser ();
+  parser.FunctionName = "curve2polyline";
+  parser = addParamValue (parser,'Nmax', 32, @(x)x>0);
+  parser = addParamValue (parser,'Tol', 1e-4, @(x)x>0);
+  parser = addParamValue (parser,'MaxIter', 10, @(x)x>0);
+  parser = parse(parser,varargin{:});
+
+  Nmax      = parser.Results.Nmax;
+  tol       = parser.Results.Tol;
+  MaxIter   = parser.Results.MaxIter;
+
+  clear parser toldef
+  # ------ #
+
+  t = [0; 1];
+  tf = 1;
+  points = 1;
+  for iter = 1:MaxIter
+    % Add parameter values where error is still bigger than tol.
+    t = interleave(t, tf);
+    nt = length(t);
+
+    % Update error
+    polyline = curveval (curve,t);
+    bump = bumpyness(polyline);
+
+    % Check which intervals must be subdivided
+    idx = find(bump > tol);
+    % The position of the bumps mpas into intervals
+    % 1 -> 1 2
+    % 2 -> 3 4
+    % 3 -> 5 6
+    % and so on
+    idx = [2*(idx-1)+1; 2*idx](:);
+    tf = false (nt-1,1);
+    tf(idx) = true;
+
+    if all (!tf)
+      break;
+    end
+
+  end
+
+endfunction
+
+function f = bumpyness (p)
+%% Check for co-linearity
+%% TODO implement various method for this
+%% -- Area of the triangle close to zero (used currently).
+%% -- Angle close to pi.
+%% -- abs(p0-pt) + abs(pt-p1) - abs(p0-p1) almost zero.
+%% -- Curve's tange at 0,t,1 are almost parallel.
+%% -- pt is in chord p0 -> p1.
+%% Do this in isParallel.m and remove this function
+
+  PL = p(1:2:end-2,:);
+  PC = p(2:2:end-1,:);
+  PR = p(3:2:end,:);
+
+  a = PL - PC;
+  b = PR - PC;
+
+  f = (a(:,1).*b(:,2) - a(:,2).*b(:,1)).^2;
+
+endfunction
+
+function tt = interleave (t,varargin)
+
+ nt   = length(t);
+ ntt = 2 * nt -1;
+ tt = zeros(ntt,1);
+ tt(1:2:ntt) = t;
+ beta = 0.4 + 0.2*rand(nt-1, 1);
+ tt(2:2:ntt) = t(1:end-1) + beta.*(t(2:end)-t(1:end-1));
+
+ if nargin > 1
+  tf = true (ntt,1);
+  tf(2:2:ntt) = varargin{1};
+  tt(!tf) = [];
+ end
+
+endfunction
+
+%!demo
+%! curve    = [0 0 1 0;1 -0.3-1 0.3 0];
+%! polyline = curve2polyline(curve,'tol',1e-8);
+%!
+%! t  = linspace(0,1,100)';
+%! pc = curveval(curve,t);
+%!
+%! plot(polyline(:,1),polyline(:,2),'-o',pc(:,1),pc(:,2),'-r')