summaryrefslogtreecommitdiff
path: root/inst/curve2polyline.m
blob: a44c3515d7eb708b281cd71b7f3742bb1fe3153e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
## Copyright (C) 2012-2019 Juan Pablo Carbajal
##
## 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 (at your option) 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 -*-
## @defun {@var{polyline} = } curve2polyline (@var{curve})
## @defunx {@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}.
##
## Currently the area of the smaller triangle formed by three consecutive points
## on the curve is considered. When this area is smaller than tolerance the
## points are colinear, and hence no more sampling between these points is
## needed.
##
## @item 'MaxIter'
## Maximum number of iterations. Default value @code{10}.
## @item 'Method'
## Not implemented.
## @end table
##
## This function is based on the algorithm described in
## L. H. de Figueiredo (1993). "Adaptive Sampling of Parametric Curves". Graphic Gems III.
## @seealso{shape2polygon, curveval}
## @end defun

## 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 ('Nmax', 32, @(x)x>0);
  parser.addParamValue ('Tol', 1e-4, @(x)x>0);
  parser.addParamValue ('MaxIter', 10, @(x)x>0);
  parser.parse(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 maps 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')