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