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