octave_packages/geometry-1.5.0/shape2d/curve2polyline.m

   1 %% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   2 %%
   3 %%    This program is free software: you can redistribute it and/or modify
   4 %%    it under the terms of the GNU General Public License as published by
   5 %%    the Free Software Foundation, either version 3 of the License, or
   6 %%    any later version.
   7 %%
   8 %%    This program is distributed in the hope that it will be useful,
   9 %%    but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 %%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 %%    GNU General Public License for more details.
  12 %%
  13 %%    You should have received a copy of the GNU General Public License
  14 %%    along with this program. If not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {@var{polyline} = } curve2polyline (@var{curve})
  18 %% @deftypefnx {Function File} {@var{polyline} = } curve2polyline (@dots{},@var{property},@var{value},@dots{})
  19 %% Adaptive sampling of a parametric curve.
  20 %%
  21 %% The @var{curve} is described as a 2-by-N matrix. Rows correspond to the
  22 %% polynomial (compatible with @code{polyval}) describing the respective component
  23 %% of the curve. The curve must be parametrized in the interval [0,1].
  24 %% The vertices of the polyline are accumulated in regions of the curve where
  25 %% the curvature is higher.
  26 %%
  27 %% @strong{Parameters}
  28 %% @table @samp
  29 %% @item 'Nmax'
  30 %% Maximum number of vertices. Not used.
  31 %% @item 'Tol'
  32 %% Tolerance for the error criteria. Default value @code{1e-4}.
  33 %% @item 'MaxIter'
  34 %% Maximum number of iterations. Default value @code{10}.
  35 %% @item 'Method'
  36 %% Not implemented.
  37 %% @end table
  38 %%
  39 %% @seealso{shape2polygon, curveval}
  40 %% @end deftypefn
  41
  42 %% This function is based on the algorithm described in
  43 %% L. H. de Figueiredo (1993). "Adaptive Sampling of Parametric Curves". Graphic Gems III.
  44 %% I had to remove the recursion so this version could be improved.
  45 %% Thursday, April 12 2012 -- JuanPi
  46
  47 function [polyline t bump]= curve2polyline (curve, varargin)
  48 %% TODO make tolerance relative to the "diameter" of the curve.
  49
  50   # --- Parse arguments --- #
  51   parser = inputParser ();
  52   parser.FunctionName = "curve2polyline";
  53   parser = addParamValue (parser,'Nmax', 32, @(x)x>0);
  54   parser = addParamValue (parser,'Tol', 1e-4, @(x)x>0);
  55   parser = addParamValue (parser,'MaxIter', 10, @(x)x>0);
  56   parser = parse(parser,varargin{:});
  57
  58   Nmax      = parser.Results.Nmax;
  59   tol       = parser.Results.Tol;
  60   MaxIter   = parser.Results.MaxIter;
  61
  62   clear parser toldef
  63   # ------ #
  64
  65   t = [0; 1];
  66   tf = 1;
  67   points = 1;
  68   for iter = 1:MaxIter
  69     % Add parameter values where error is still bigger than tol.
  70     t = interleave(t, tf);
  71     nt = length(t);
  72
  73     % Update error
  74     polyline = curveval (curve,t);
  75     bump = bumpyness(polyline);
  76
  77     % Check which intervals must be subdivided
  78     idx = find(bump > tol);
  79     % The position of the bumps mpas into intervals
  80     % 1 -> 1 2
  81     % 2 -> 3 4
  82     % 3 -> 5 6
  83     % and so on
  84     idx = [2*(idx-1)+1; 2*idx](:);
  85     tf = false (nt-1,1);
  86     tf(idx) = true;
  87
  88     if all (!tf)
  89       break;
  90     end
  91
  92   end
  93
  94 endfunction
  95
  96 function f = bumpyness (p)
  97 %% Check for co-linearity
  98 %% TODO implement various method for this
  99 %% -- Area of the triangle close to zero (used currently).
 100 %% -- Angle close to pi.
 101 %% -- abs(p0-pt) + abs(pt-p1) - abs(p0-p1) almost zero.
 102 %% -- Curve's tange at 0,t,1 are almost parallel.
 103 %% -- pt is in chord p0 -> p1.
 104 %% Do this in isParallel.m and remove this function
 105
 106   PL = p(1:2:end-2,:);
 107   PC = p(2:2:end-1,:);
 108   PR = p(3:2:end,:);
 109
 110   a = PL - PC;
 111   b = PR - PC;
 112
 113   f = (a(:,1).*b(:,2) - a(:,2).*b(:,1)).^2;
 114
 115 endfunction
 116
 117 function tt = interleave (t,varargin)
 118
 119  nt   = length(t);
 120  ntt = 2 * nt -1;
 121  tt = zeros(ntt,1);
 122  tt(1:2:ntt) = t;
 123  beta = 0.4 + 0.2*rand(nt-1, 1);
 124  tt(2:2:ntt) = t(1:end-1) + beta.*(t(2:end)-t(1:end-1));
 125
 126  if nargin > 1
 127   tf = true (ntt,1);
 128   tf(2:2:ntt) = varargin{1};
 129   tt(!tf) = [];
 130  end
 131
 132 endfunction
 133
 134 %!demo
 135 %! curve    = [0 0 1 0;1 -0.3-1 0.3 0];
 136 %! polyline = curve2polyline(curve,'tol',1e-8);
 137 %!
 138 %! t  = linspace(0,1,100)';
 139 %! pc = curveval(curve,t);
 140 %!
 141 %! plot(polyline(:,1),polyline(:,2),'-o',pc(:,1),pc(:,2),'-r')