diff options
| author | Rafael Laboissiere <rafael@laboissiere.net> | 2012-06-10 08:04:00 +0000 |
|---|---|---|
| committer | Rafael Laboissiere <rafael@laboissiere.net> | 2012-06-10 08:04:00 +0000 |
| commit | 24a8d839d5a45b6b38ec6704d96af8cfad55123f (patch) | |
| tree | dffe8c6d8eab8cbc82269009cd612fb240a69405 | |
| parent | 1ba9420dac11ebcd5ef87843f1322b518d2d1071 (diff) | |
Imported Upstream version 1.5.0
45 files changed, 2259 insertions, 933 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 42f0e16..bfeb0c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,11 @@ Name: Geometry -Version: 1.4.1 -Date: 2012-03-24 -Author: David Legland <david.legland@grignon.inra.fr>, José Luis García Pallero <jgpallero@gmail.com>, Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +Version: 1.5.0 +Date: 2012-06-05 +Author: David Legland <david.legland@grignon.inra.fr>, José Luis García Pallero <jgpallero@gmail.com>, Juan Pablo Carbajal <carbajal@ifi.uzh.ch>, Simeon Simeonov <sss41@cam.ac.uk> Maintainer: Juan Pablo Carbajal <carbajal@ifi.uzh.ch> Title: Computational Geometry Description: Library for geometric computing extending MatGeom functions. Useful to create, transform, manipulate and display geometric primitives. -Depends: octave (>= 3.4.0) +Depends: octave (>= 3.4.0), general (>= 1.3.0) Autoload: yes License: GPLv3+, FreeBSD Url: http://octave.sf.net, http://matgeom.sf.net, http://davis.wpi.edu/~matt/courses/clipping/, https://bitbucket.org/jgpallero/octclip @@ -14,8 +14,10 @@ geometry >> Computational Geometry polygons2d 2D Points centroid + closed_path distancePoints drawPoint + isCounterClockwise isPointOnRay isPointInCircle isPointOnCircle @@ -96,9 +98,11 @@ geometry >> Computational Geometry intersectLineCircle radicalAxis 2D Ellipses - ellipseAsPolygon + cov2ellipse drawEllipseArc drawEllipse + ellipseAsPolygon + ellipse2cov inertiaEllipse 2D Transformations transformLine @@ -115,16 +119,34 @@ geometry >> Computational Geometry drawBezierCurve cbezier2poly 2D Polygons + curvature + distancePolygons + distancePointPolygon + distancePointPolyline drawPolygon - simplifypolygon + expandPolygon + medialAxisConvex oc_polybool + parametrize + polygonLoops + polygonSelfIntersections + polylineSelfIntersections + reversePolyline + reversePolygon + simplifypolygon + simplifypolyline + splitPolygons + supportFunction 2D Piecewise polynomial shapes polygon2shape shape2polygon shapecentroid shapeplot shapetransform + curveval + curve2polyline 2D Others + beltproblem bisector crackPattern2 crackPattern @@ -133,7 +155,6 @@ geometry >> Computational Geometry drawParabola drawShape hexagonalGrid - isCounterClockwise squareGrid triangleGrid Geometric graphs creation @@ -1,6 +1,51 @@ Summary of important user-visible changes for releases of the geometry package =============================================================================== +geometry-1.5.0 Release Date: 2012-06-05 Release Manager: Juan Pablo Carbajal +=============================================================================== + +* Added functions: + - cov2ellipse.m & ellipse2cov: transform between ellipses and covariances matrices. + + - beltproblem.m : Finds the four lines tangent to two circles with given centers and + radii. This is the solution to the belt problem in 2D. + + - curveval.m: Evaluates a polynomial curve defined as a 2-by-N matrix. + + - curve2polyline.m: Converts a polynomial curve into a polyline by the adaptive + sampling method. + + - simplifypolyline.m: Ramer-Douglas-Peucker algorithm to simplify polylines. + + - parametrize.m: Estimate a parametrization of a polygon/line based on the distance + between the points. + + - curvature.m: Estimation of the curvature of a polygon/line based on polynomial + approximation. + + - reversePolygon.m and reversePolyline.m: reverse the orders of the points in + of polygon/line. + + - supportFunction.m: Compute support function of a polygon. + + - distancePointPolygon.m, distancePointPolyline.m, distancePolygons.m , + expandPolygon.m, medialAxisConvex.m, polygonLoops.m, polygonSelfIntersections.m + polylineSelfIntersections.m, splitPolygons.m + + - close_path.m : given a set of points in the plane calculate a piecewise linear + simple path that passes through all points. + +* Changed functions: + - distancePointEdge.m: Now the function computes the distance between all points + and all edges. A third optional argument provides + backward compatibility. + +* Solved bugs: + - simplifypolygon.m returned empty polygons when points are repeated, i.e when + the polygon is not correctly formed. + - Removed installation warnings. + +=============================================================================== geometry-1.4.1 Release Date: 2012-03-24 Release Manager: Juan Pablo Carbajal =============================================================================== @@ -2,21 +2,34 @@ dirlist = {"geom2d","io","polygons2d","shape2d","octclip", "graphs"}; dirname = fileparts (canonicalize_file_name (mfilename ("fullpath"))); pp = strsplit (dirname,filesep (), true); - -%% Check if prefix was used -[pkg_folder dep_folder] = pkg ("prefix"); -pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ]; -dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ]; - -%% If we are in Architecture dependent folder add from outside arch = cstrcat (octave_config_info ("canonical_host_type"), "-", octave_config_info ("api_version")); -if strcmp (arch , pp{end}) - dirname = [strcat(filesep(),{pp{1:end-1}}){:}]; - pkg_folder = strsplit (pkg_folder,filesep (), true); - pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}]; + +%% Get the correct path +% Search installed packages +[local_packages, global_packages] = pkg("list"); +installed_pkgs_lst = {local_packages{:}, global_packages{:}}; +pkg_data = installed_pkgs_lst (cellfun(@(x) ismember (x.name, {"geometry"}), ... + installed_pkgs_lst, "unif", true)); +if isempty(pkg_data) + % The package is not installed yet + [pkg_folder dep_folder] = pkg ("prefix"); + pkg_folder = [pkg_folder strcat(filesep(),{pp{end-1:end}}){:} ]; +% dep_folder = [dep_folder strcat(filesep(),{pp{end-1:end}}){:} ]; + if strcmp (arch , pp{end}) + %% If we are in Architecture dependent folder add from outside + pkg_folder = strsplit (pkg_folder,filesep (), true); + pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}]; +% dep_folder = strsplit (dep_folder,filesep (), true); +% dep_folder = [strcat(filesep(),{dep_folder{1:end-1}}){:}]; + end +else + pkg_folder = pkg_data{1}.dir; +% dep_folder = pkg_data{1}.archprefix; end +%dep_folder = [dep_folder filesep arch]; + if (! exist (fullfile (dirname, "inst"), "dir")) %% Installing for ii=1:length (dirlist) @@ -30,6 +43,8 @@ else addpath ([ dirname "/inst/" dirlist{ii}]) endfor endif +%addpath (dep_folder) warning('off', 'Octave:fopen-file-in-path'); -clear dirlist dirname pp +clear dirlist dirname pp arch pkg_folder +%dep_folder @@ -2,22 +2,33 @@ dirlist = {"geom2d","io","polygons2d","shape2d","octclip","graphs"}; dirname = fileparts (canonicalize_file_name (mfilename ("fullpath"))); pp = strsplit (dirname,filesep (), true); - -%% Check if prefix was used -[pkg_folder dep_folder] = pkg ("prefix"); -pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ]; -dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ]; - -%% If we are not in Architecture dependent folder arch = cstrcat (octave_config_info ("canonical_host_type"), "-", octave_config_info ("api_version")); -pp = strsplit (dirname,filesep (), true); -if strcmp(arch , pp{end}) - dirname = [pkg("prefix") filesep() pp{end-1}]; - pkg_folder = strsplit (pkg_folder,filesep (), true); - pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}]; +%% Get the correct path +% Search installed packages +[local_packages, global_packages] = pkg("list"); +installed_pkgs_lst = {local_packages{:}, global_packages{:}}; +pkg_data = installed_pkgs_lst (cellfun(@(x) ismember (x.name, {"geometry"}), ... + installed_pkgs_lst, "unif", true)); +if isempty(pkg_data) + % The package is not installed yet + [pkg_folder dep_folder] = pkg ("prefix"); + pkg_folder = [pkg_folder strcat(filesep(),{pp{end-1:end}}){:} ]; +% dep_folder = [dep_folder strcat(filesep(),{pp{end-1:end}}){:} ]; + if strcmp (arch , pp{end}) + %% If we are in Architecture dependent folder add from outside + pkg_folder = strsplit (pkg_folder,filesep (), true); + pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}]; +% dep_folder = strsplit (dep_folder,filesep (), true); +% dep_folder = [strcat(filesep(),{dep_folder{1:end-1}}){:}]; + end +else + pkg_folder = pkg_data{1}.dir; +% dep_folder = pkg_data{1}.archprefix; end +%dep_folder = [dep_folder filesep arch]; + if (! exist (fullfile (dirname, "inst"), "dir")) ## Run this if the package is installed for ii=1:length (dirlist) @@ -29,5 +40,7 @@ else rmpath ([ dirname "/inst/" dirlist{ii}]) endfor endif +%rmpath (dep_folder); -clear dirlist dirname pp +clear dirlist dirname pp arch pkg_folder +%dep_folder diff --git a/inst/geom2d/beltproblem.m b/inst/geom2d/beltproblem.m new file mode 100644 index 0000000..d3207bd --- /dev/null +++ b/inst/geom2d/beltproblem.m @@ -0,0 +1,136 @@ +%% 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{tangent},@var{inner}] = } beltproblem (@var{c}, @var{r}) +%% Finds the four lines tangent to two circles with given centers and radii. +%% +%% The function solves the belt problem in 2D for circles with center @var{c} and +%% radii @var{r}. +%% +%% @strong{INPUT} +%% @table @var +%% @item c +%% 2-by-2 matrix containig coordinates of the centers of the circles; one row per circle. +%% @item r +%% 2-by-1 vector with the radii of the circles. +%%@end table +%% +%% @strong{OUPUT} +%% @table @var +%% @item tangent +%% 4-by-4 matrix with the points of tangency. Each row describes a segment(edge). +%% @item inner +%% 4-by-2 vector with the point of intersection of the inner tangents (crossed belts) +%% with the segment that joins the centers of the two circles. If +%% the i-th edge is not an inner tangent then @code{inner(i,:)=[NaN,NaN]}. +%% @end table +%% +%% Example: +%% +%% @example +%% c = [0 0;1 3]; +%% r = [1 0.5]; +%% [T inner] = beltproblem(c,r) +%% @result{} T = +%% -0.68516 0.72839 1.34258 2.63581 +%% 0.98516 0.17161 0.50742 2.91419 +%% 0.98675 -0.16225 1.49338 2.91888 +%% -0.88675 0.46225 0.55663 3.23112 +%% +%% @result{} inner = +%% 0.66667 2.00000 +%% 0.66667 2.00000 +%% NaN NaN +%% NaN NaN +%% +%% @end example +%% +%% @seealso{edges2d} +%% @end deftypefn + +function [edgeTan inner] = beltproblem(c,r) + + x0 = [c(1,1) c(1,2) c(2,1) c(2,2)]; + r0 = r([1 1 2 2]); + + middleEdge = createEdge(c(1,:),c(2,:)); + + ind0 = [1 0 1 0; 0 1 1 0; 1 1 1 0; -1 0 1 0; 0 -1 1 0; -1 -1 1 0;... + 1 0 0 1; 0 1 0 1; 1 1 0 1; -1 0 0 1; 0 -1 0 1; -1 -1 0 1;... + 1 0 1 1; 0 1 1 1; 1 1 1 1; -1 0 1 1; 0 -1 1 1; -1 -1 1 1;... + 1 0 -1 0; 0 1 -1 0; 1 1 -1 0; -1 0 -1 0; 0 -1 -1 0; -1 -1 -1 0;... + 1 0 0 -1; 0 1 0 -1; 1 1 0 -1; -1 0 0 -1; 0 -1 0 -1; -1 -1 0 -1;... + 1 0 -1 -1; 0 1 -1 -1; 1 1 -1 -1; -1 0 -1 -1; 0 -1 -1 -1; -1 -1 -1 -1]; + nInit = size(ind0,1); + ir = randperm(nInit); + edgeTan = zeros(4,4); + inner = zeros(4,2); + nSol = 0; + i=1; + + %% Solve for 2 different lines + while nSol<4 && i<nInit + ind = find(ind0(ir(i),:)~=0); + x = x0; + x(ind)=x(ind)+r0(ind); + [sol f0 nev]= fsolve(@(x)belt(x,c,r),x); + if nev~=1 + perror('fsolve',nev) + end + + for j=1:4 + notequal(j) = all(abs(edgeTan(j,:)-sol) > 1e-6); + end + if all(notequal) + nSol = nSol+1; + edgeTan(nSol,:) = createEdge(sol(1:2),sol(3:4)); + % Find innerTangent + inner(nSol,:) = intersectEdges(middleEdge,edgeTan(nSol,:)); + end + + i = i+1; + end + +endfunction + +function res = belt(x,c,r) + res = zeros(4,1); + + res(1,1) = (x(3:4) - c(2,1:2))*(x(3:4) - x(1:2))'; + res(2,1) = (x(1:2) - c(1,1:2))*(x(3:4) - x(1:2))'; + res(3,1) = sumsq(x(1:2) - c(1,1:2)) - r(1)^2; + res(4,1) = sumsq(x(3:4) - c(2,1:2)) - r(2)^2; + +end + +%!demo +%! c = [0 0;1 3]; +%! r = [1 0.5]; +%! [T inner] = beltproblem(c,r) +%! +%! figure(1) +%! clf +%! h = drawEdge(T); +%! set(h(find(~isnan(inner(:,1)))),'color','r'); +%! set(h,'linewidth',2); +%! hold on +%! drawCircle([c(1,:) r(1); c(2,:) r(2)],'linewidth',2); +%! axis tight +%! axis equal +%! +%! % ------------------------------------------------------------------- +%! % The circles with the tangents edges are plotted. The solution with +%! % crossed belts (inner tangets) is shown in red. diff --git a/inst/geom2d/closed_path.m b/inst/geom2d/closed_path.m new file mode 100644 index 0000000..4ea3839 --- /dev/null +++ b/inst/geom2d/closed_path.m @@ -0,0 +1,105 @@ +%% Copyright (C) 2012 Simeon Simeonov <simeon.simeonov.s@gmail.com> +%% +%% 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{y} =} polygon (@var{x}) +%% Returns a simple closed path that passes through all the points in @var{x}. +%% @var{x} is a vector containing 2D coordinates of the points. +%% +%% @end deftypefn + +%% Author: Simeon Simeonov <simeon.simeonov.s@gmail.com> + +function y = closed_path(x) + + if(size(x,1) > 1 && size(x,1) < size(x,2)) + x = x'; + end + + N = size(x,1); % Number of points + idx = zeros(N, 1); % ind contains the indices of the sorted coordinates + + a = find(x(:,2)==min(x(:,2))); + + if(size(a,1) > 1) + [~, i] = sort(x(a,1)); + a = a(i(1)); + end + + x_1 = x(x(:,2)==x(a,2),:); % find all x with the same y coordinate + + if(x(a,1) == min(x(:,1))) + x_2 = x(x(:,1)==x(a,1),:); % find all x with the same x coordinate + else + x_2 = x(a,:); + end + + if(size(x_1,1) > 1 || size(x_2,1) > 1) + if(size(x_1,1) > 1) + x_1 = sort(x_1); % Sort by x coordinate + y(1,:) = x(a,:); % original starting point + end + + if (size(x_2,1) > 1) + x_2 = sort(x_2, 'descend'); + end + + x_not = [x_1; x_2]; + i = ismember(x,x_not,'rows'); + x(i, :) = []; + x = [x_1(size(x_1,1),:); x]; + x_1(size(x_1, 1),:) = []; + N = size(x,1); + a = 1; + else + x_1 = []; + x_2 = x(a,:); + end + d = x - repmat(x(a,:), N, 1); + th = d(:,2)./(d(:,1) + d(:,2)); + + [~, idx0] = ismember(sort(th(th==0)), th); + [~, idx1] = ismember(sort(th(th>0)), th); + [~, idx2] = ismember(sort(th(th<0)), th); + + idx = [a; idx0; idx1; idx2]; + % I contains the indices of idx in a sorted order. [v i] = sort(idx) then + % i==I. + [~,I,J]= unique(idx); + if(size(I,1) ~= size(J,1)) + R = histc(J, 1:size(I,1)); % R(1) will always be 1? + idx_sorted = idx(I); + r = find(R>1); + for ri = r' + idx_repeated = idx_sorted(ri); + idx(idx==idx_repeated) = find(th==th(idx_sorted(ri))); + end + end + + y = [x_1; x(idx,:); x_2;]; + +endfunction + +%!demo +%! maxInt = 100; +%! N = 25; +%! +%! for i = 1:5 +%! x = randi(maxInt, N, 2); +%! y = closed_path(x); +%! plot(y(:,1), y(:,2), '*'); +%! hold on; +%! plot(y(:,1), y(:,2)); +%! end diff --git a/inst/geom2d/cov2ellipse.m b/inst/geom2d/cov2ellipse.m new file mode 100644 index 0000000..214381a --- /dev/null +++ b/inst/geom2d/cov2ellipse.m @@ -0,0 +1,71 @@ +## 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 +## (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 -*- +%% @deftypefn {Function File} {@var{ellipse} = } cov2ellipse (@var{K}) +%% @deftypefnx {Function File} {[@var{ra} @var{rb} @var{theta}] = } cov2ellipse (@var{K}) +%% @deftypefnx {Function File} {@dots{} = } cov2ellipse (@dots{}, @samp{tol},@var{tol}) +%% Calculates ellipse parameters from covariance matrix. +%% +%% @var{K} must be symmetric positive (semi)definite. The optional argument +%% @samp{tol} sets the tolerance for the verification of the +%% positive-(semi)definiteness of the matrix @var{K} (see @command{isdefinite}). +%% +%% If only one output argument is supplied a vector defining a ellipse is returned +%% as defined in @command{ellipses2d}. Otherwise the angle @var{theta} is given +%% in radians. +%% +%% Run @code{demo cov2ellipse} to see an example. +%% +%% @seealso{ellipses2d, cov2ellipse, drawEllipse} +%% @end deftypefn + +function varargout = cov2ellipse (K, varargin); + + parser = inputParser (); + parser.FunctionName = "cov2ellipse"; + parser = addParamValue (parser,'Tol', 100*eps*norm (K, "fro"), @(x)x>0); + parser = parse(parser,varargin{:}); + + if isdefinite (K,parser.Results.Tol) == -1 + print_usage + end + [R S W] = svd (K); + theta = atan (R(1,1)/R(2,2)); + v = sort (diag(S), 'ascend')'; + + if nargout == 1 + varargout{1} = [0 0 v theta*180/pi]; + elseif nargout == 3 + varargout{1} = v(1); + varargout{2} = v(2); + varargout{3} = theta; + end + +endfunction + +%!demo +%! K = [2 1; 1 2]; +%! L = chol(K,'lower'); +%! u = randn(1e3,2)*L'; +%! +%! elli = cov2ellipse (K) +%! +%! figure(1) +%! plot(u(:,1),u(:,2),'.r'); +%! hold on; +%! drawEllipse(elli,'linewidth',2); +%! hold off +%! axis tight diff --git a/inst/geom2d/distancePointEdge.m b/inst/geom2d/distancePointEdge.m index 1a689db..709a2bf 100644 --- a/inst/geom2d/distancePointEdge.m +++ b/inst/geom2d/distancePointEdge.m @@ -1,6 +1,5 @@ -%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+%% Copyright (c) 2007-2011, David Legland <david.legland@grignon.inra.fr>
+%% Copyright (c) 2012, Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
%%
%% All rights reserved.
%% (simplified BSD License)
@@ -10,8 +9,8 @@ %%
%% 1. Redistributions of source code must retain the above copyright notice, this
%% list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
+%%
+%% 2. Redistributions in binary form must reproduce the above copyright notice,
%% this list of conditions and the following disclaimer in the documentation
%% and/or other materials provided with the distribution.
%%
@@ -19,9 +18,9 @@ %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
@@ -33,53 +32,102 @@ %% -*- texinfo -*-
%% @deftypefn {Function File} {@var{dist} = } distancePointEdge (@var{point}, @var{edge})
+%% @deftypefnx {Function File} {@var{dist} = } distancePointEdge (@dots{}, @var{opt})
+%% @deftypefnx {Function File} {[@var{dist} @var{pos}]= } distancePointEdge (@dots{})
%% Minimum distance between a point and an edge
%%
-%% DIST = distancePointEdge(POINT, EDGE);
-%% Return the euclidean distance between edge EDGE and point POINT.
-%% EDGE has the form: [x1 y1 x2 y2], and POINT is [x y].
+%% Return the euclidean distance between edge @var{edge} and point @var{point}.
+%% @var{edge} has the form: [x1 y1 x2 y2], and @var{point} is [x y].
+%% If @var{edge} is Ne-by-4 and @var{point} is Np-by-2, then @var{dist} is
+%% Np-by-Ne, where each row contains the distance of each point to all the edges.
%%
-%% If EDGE is N-by-4 array, result is N-by-1 array computed for each edge.
-%% If POINT is a N-by-2 array, the result is computed for each point.
-%% If both POINT and EDGE are array, they must have the same number of
-%% rows, and the result is computed for each couple point(i,:);edge(i,:).
+%% If @var{opt} is true (or equivalent), the optput is cmpatible with the original function:
+%% @table @samp
+%% @item 1
+%% If @var{point} is 1-by-2 array, the result is Ne-by-1 array computed for each edge.
+%% @item 2
+%% If @var{edge} is a 1-by-4 array, the result is Np-by-1 computed for each point.
+%% @item 3
+%% If both @var{point} and @var{edge} are array, they must have the same number of
+%% rows, and the result is computed for each couple @code{@var{point}(i,:),@var{edge}(i,:)}.
+%% @end table
%%
-%% [DIST POS] = distancePointEdge(POINT, EDGE);
-%% Also returns the position of closest point on the edge. POS is
-%% comprised between 0 (first point) and 1 (last point).
+%% If the the second output argument @var{pos} is requested, the function also
+%% returns the position of closest point on the edge. @var{pos} is comprised
+%% between 0 (first point) and 1 (last point).
%%
-%% @seealso{edges2d, points2d, distancePoints, distancePointLine}
+%% @seealso{edges2d, points2d, distancePoints, distancePointLine}
%% @end deftypefn
-function varargout = distancePointEdge(point, edge)
+%% Rewritten to accept different numbers of points and edges.
+%% 2012 - Juan Pablo Carbajal
+
+function [dist, tp] = distancePointEdge(point, edge, opt=[])
+
+ Np = size (point,1);
+ Ne = size (edge,1);
+ edge = edge.';
+
% direction vector of each edge
- dx = edge(:, 3) - edge(:,1);
- dy = edge(:, 4) - edge(:,2);
+ dx = edge(3,:) - edge(1,:);
+ dy = edge(4,:) - edge(2,:);
% compute position of points projected on the supporting line
% (Size of tp is the max number of edges or points)
- delta = dx .* dx + dy .* dy;
- tp = ((point(:, 1) - edge(:, 1)) .* dx + (point(:, 2) - edge(:, 2)) .* dy) ./ delta;
- % ensure degenerated edges are correclty processed (consider the first
- % vertex is the closest)
- tp(delta < eps) = 0;
+ delta = dx .* dx + dy .* dy;
+ mask = delta < eps;
+ delta(mask) = 1;
+ warning ('off', 'Octave:broadcast');
+ tp = ((point(:, 1) - edge(1, :)) .* dx + (point(:, 2) - edge(2, :)) .* dy) ./ delta;
+ tp(:,mask) = 0;
% change position to ensure projected point is located on the edge
tp(tp < 0) = 0;
tp(tp > 1) = 1;
% coordinates of projected point
- p0 = [edge(:,1) + tp .* dx, edge(:,2) + tp .* dy];
+ p0x = (edge(1,:) + tp .* dx);
+ p0y = (edge(2,:) + tp .* dy);
% compute distance between point and its projection on the edge
- dist = sqrt((point(:,1) - p0(:,1)) .^ 2 + (point(:,2) - p0(:,2)) .^ 2);
+ dist = sqrt((point(:,1) - p0x) .^ 2 + (point(:,2) - p0y) .^ 2);
+
+ warning ('on', 'Octave:broadcast');
- % process output arguments
- varargout{1} = dist;
- if nargout > 1
- varargout{2} = tp;
+ %% backwards compatibility
+ if opt
+ if Np != Ne && (Ne != 1 && Np !=1)
+ error ("geometry:InvalidArgument", ...
+ "Sizes must be equal or one of the inputs must be 1-by-N array.");
+ end
+ if Np == Ne
+ dist = diag(dist)(:);
+ tp = diag(tp)(:);
+ elseif Np == 1
+ dist = dist.';
+ tp = tp.';
+ end
end
endfunction
+%% Original code
+%%tp = ((point(:, 1) - edge(:, 1)) .* dx + (point(:, 2) - edge(:, 2)) .* dy) ./ delta;
+%%% ensure degenerated edges are correclty processed (consider the first
+%%% vertex is the closest)
+%%if Ne > 1
+%% tp(delta < eps) = 0;
+%%elseif delta < eps
+%% tp(:) = 0;
+%%end
+
+%%% change position to ensure projected point is located on the edge
+%%tp(tp < 0) = 0;
+%%tp(tp > 1) = 1;
+
+%%% coordinates of projected point
+%%p0 = [edge(:,1) + tp .* dx, edge(:,2) + tp .* dy];
+
+%%% compute distance between point and its projection on the edge
+%%dist = sqrt((point(:,1) - p0(:,1)) .^ 2 + (point(:,2) - p0(:,2)) .^ 2);
diff --git a/inst/geom2d/ellipse2cov.m b/inst/geom2d/ellipse2cov.m new file mode 100644 index 0000000..276b411 --- /dev/null +++ b/inst/geom2d/ellipse2cov.m @@ -0,0 +1,92 @@ +## 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 +## (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 -*- +%% @deftypefn {Function File} {@var{K} = } ellipse2cov (@var{elli}) +%% @deftypefnx {Function File} {@var{K} = } ellipse2cov (@var{ra}, @var{rb}) +%% @deftypefnx {Function File} {@var{K} = } ellipse2cov (@dots{}, @var{theta}) +%% Calculates covariance matrix from ellipse. +%% +%% If only one input is given, @var{elli} must define an ellipse as described in +%% @command{ellipses2d}. +%% If two inputs are given, @var{ra} and @var{rb} define the half-lenght of the +%% axes. +%% If a third input is given, @var{theta} must be the angle of rotation of the +%% ellipse in radians, and in counter-clockwise direction. +%% +%% The output @var{K} contains the covariance matrix define by the ellipse. +%% +%% Run @code{demo ellipse2cov} to see an example. +%% +%% @seealso{ellipses2d, cov2ellipse, drawEllipse} +%% @end deftypefn + +function K = ellipse2cov (elli, varargin); + + ra = 1; + rb = 1; + theta = 0; + switch numel (varargin) + case 0 + %% ellipse format + if numel (elli) != 5 + print_usage (); + end + ra = elli(1,3); + rb = elli(1,4); + theta = elli(1,5)*pi/180; + + case 2 + %% ra,rb + if numel (elli) != 1 + print_usage (); + end + ra = elli; + rb = varargin{1}; + + case 3 + %% ra,rb, theta + if numel (elli) != 1 + print_usage (); + end + ra = elli; + rb = varargin{1}; + theta = varargin{2}; + + otherwise + print_usage (); + end + + T = createRotation (theta)(1:2,1:2); + K = T*diag([ra rb])*T'; + +endfunction + +%!demo +%! elli = [0 0 1 3 -45]; +%! +%! % Create 2D normal random variables with covarinace defined by elli. +%! K = ellipse2cov (elli) +%! L = chol(K,'lower'); +%! u = randn(1e3,2)*L'; +%! +%! Kn = cov (u) +%! +%! figure(1) +%! plot(u(:,1),u(:,2),'.r'); +%! hold on; +%! drawEllipse(elli,'linewidth',2); +%! hold off +%! axis tight diff --git a/inst/geom2d/geometry_Contents.m b/inst/geom2d/geometry_Contents.m deleted file mode 100644 index b87427f..0000000 --- a/inst/geom2d/geometry_Contents.m +++ /dev/null @@ -1,224 +0,0 @@ -%% Copyright (c) 2011, INRA
-%% 2007-2011, David Legland <david.legland@grignon.inra.fr>
-%% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%%
-%% All rights reserved.
-%% (simplified BSD License)
-%%
-%% Redistribution and use in source and binary forms, with or without
-%% modification, are permitted provided that the following conditions are met:
-%%
-%% 1. Redistributions of source code must retain the above copyright notice, this
-%% list of conditions and the following disclaimer.
-%%
-%% 2. Redistributions in binary form must reproduce the above copyright notice,
-%% this list of conditions and the following disclaimer in the documentation
-%% and/or other materials provided with the distribution.
-%%
-%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-%% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-%% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-%% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-%% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-%% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-%% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-%% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-%% POSSIBILITY OF SUCH DAMAGE.
-%%
-%% The views and conclusions contained in the software and documentation are
-%% those of the authors and should not be interpreted as representing official
-%% policies, either expressed or implied, of copyright holder.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} geometry_Contents ()
-%% Geometry 2D Toolbox
-%% Version 1.2.0 21-Oct-2011 .
-%%
-%% Library to handle and visualize geometric primitives such as points,
-%% lines, circles and ellipses, polygons...
-%%
-%% The goal is to provide a low-level library for manipulating geometrical
-%% primitives, making easier the development of more complex geometric
-%% algorithms.
-%%
-%% Most functions works for planar shapes, but some ones have been
-%% extended to 3D or to any dimension.
-%%
-%% Points
-%% points2d - Description of functions operating on points
-%% clipPoints - Clip a set of points by a box
-%% centroid - Compute centroid (center of mass) of a set of points
-%% midPoint - Middle point of two points or of an edge
-%% isCounterClockwise - Compute relative orientation of 3 points
-%% polarPoint - Create a point from polar coordinates (rho + theta)
-%% angle2Points - Compute horizontal angle between 2 points
-%% angle3Points - Compute oriented angle made by 3 points
-%% angleSort - Sort points in the plane according to their angle to origin
-%% distancePoints - Compute distance between two points
-%% minDistancePoints - Minimal distance between several points
-%% transformPoint - Transform a point with an affine transform
-%% drawPoint - Draw the point on the axis.
-%%
-%% Vectors
-%% vectors2d - Description of functions operating on plane vectors
-%% createVector - Create a vector from two points
-%% vectorNorm - Compute norm of a vector, or of a set of vectors
-%% vectorAngle - Angle of a vector, or between 2 vectors
-%% normalizeVector - Normalize a vector to have norm equal to 1
-%% isPerpendicular - Check orthogonality of two vectors
-%% isParallel - Check parallelism of two vectors
-%% transformVector - Transform a vector with an affine transform
-%% rotateVector - Rotate a vector by a given angle
-%%
-%% Straight lines
-%% lines2d - Description of functions operating on planar lines
-%% createLine - Create a straight line from 2 points, or from other inputs
-%% medianLine - Create a median line between two points
-%% cartesianLine - Create a straight line from cartesian equation coefficients
-%% orthogonalLine - Create a line orthogonal to another one.
-%% parallelLine - Create a line parallel to another one.
-%% intersectLines - Return all intersection points of N lines in 2D
-%% lineAngle - Computes angle between two straight lines
-%% linePosition - Position of a point on a line
-%% lineFit - Fit a straight line to a set of points
-%% clipLine - Clip a line with a box
-%% reverseLine - Return same line but with opposite orientation
-%% transformLine - Transform a line with an affine transform
-%% drawLine - Draw the line on the current axis
-%%
-%% Edges (line segments between 2 points)
-%% edges2d - Description of functions operating on planar edges
-%% createEdge - Create an edge between two points, or from a line
-%% edgeToLine - Convert an edge to a straight line
-%% edgeAngle - Return angle of edge
-%% edgeLength - Return length of an edge
-%% midPoint - Middle point of two points or of an edge
-%% edgePosition - Return position of a point on an edge
-%% clipEdge - Clip an edge with a rectangular box
-%% reverseEdge - Intervert the source and target vertices of edge
-%% intersectEdges - Return all intersections between two set of edges
-%% intersectLineEdge - Return intersection between a line and an edge
-%% transformEdge - Transform an edge with an affine transform
-%% drawEdge - Draw an edge given by 2 points
-%% drawCenteredEdge - Draw an edge centered on a point
-%%
-%% Rays
-%% rays2d - Description of functions operating on planar rays
-%% createRay - Create a ray (half-line), from various inputs
-%% bisector - Return the bisector of two lines, or 3 points
-%% clipRay - Clip a ray with a box
-%% drawRay - Draw a ray on the current axis
-%%
-%% Relations between points and lines
-%% distancePointEdge - Minimum distance between a point and an edge
-%% distancePointLine - Minimum distance between a point and a line
-%% projPointOnLine - Project of a point orthogonally onto a line
-%% pointOnLine - Create a point on a line at a given position on the line
-%% isPointOnLine - Test if a point belongs to a line
-%% isPointOnEdge - Test if a point belongs to an edge
-%% isPointOnRay - Test if a point belongs to a ray
-%% isLeftOriented - Test if a point is on the left side of a line
-%%
-%% Circles
-%% circles2d - Description of functions operating on circles
-%% createCircle - Create a circle from 2 or 3 points
-%% createDirectedCircle - Create a directed circle
-%% intersectCircles - Intersection points of two circles
-%% intersectLineCircle - Intersection point(s) of a line and a circle
-%% circleAsPolygon - Convert a circle into a series of points
-%% circleArcAsCurve - Convert a circle arc into a series of points
-%% isPointInCircle - Test if a point is located inside a given circle
-%% isPointOnCircle - Test if a point is located on a given circle.
-%% enclosingCircle - Find the minimum circle enclosing a set of points.
-%% radicalAxis - Compute the radical axis (or radical line) of 2 circles
-%% drawCircle - Draw a circle on the current axis
-%% drawCircleArc - Draw a circle arc on the current axis
-%%
-%% Ellipses
-%% ellipses2d - Description of functions operating on ellipses
-%% inertiaEllipse - Inertia ellipse of a set of points
-%% isPointInEllipse - Check if a point is located inside a given ellipse
-%% ellipseAsPolygon - Convert an ellipse into a series of points
-%% drawEllipse - Draw an ellipse on the current axis
-%% drawEllipseArc - Draw an ellipse arc on the current axis
-%%
-%% Geometric transforms
-%% transforms2d - Description of functions operating on transforms
-%% createTranslation - Create the 3*3 matrix of a translation
-%% createRotation - Create the 3*3 matrix of a rotation
-%% createScaling - Create the 3*3 matrix of a scaling in 2 dimensions
-%% createHomothecy - Create the the 3x3 matrix of an homothetic transform
-%% createBasisTransform - Compute matrix for transforming a basis into another basis
-%% createLineReflection - Create the the 3x3 matrix of a line reflection
-%% fitAffineTransform2d - Fit an affine transform using two point sets
-%%
-%% Angles
-%% angles2d - Description of functions for manipulating angles
-%% normalizeAngle - Normalize an angle value within a 2*PI interval
-%% angleAbsDiff - Absolute difference between two angles
-%% angleDiff - Difference between two angles
-%% deg2rad - Convert angle from degrees to radians
-%% rad2deg - Convert angle from radians to degrees
-%%
-%% Boxes
-%% boxes2d - Description of functions operating on bounding boxes
-%% intersectBoxes - Intersection of two bounding boxes
-%% mergeBoxes - Merge two boxes, by computing their greatest extent
-%% randomPointInBox - Generate random point within a box
-%% drawBox - Draw a box defined by coordinate extents
-%%
-%% Various drawing functions
-%% drawBezierCurve - Draw a cubic bezier curve defined by 4 control points
-%% drawParabola - Draw a parabola on the current axis
-%% drawOrientedBox - Draw centered oriented rectangle
-%% drawRect - Draw rectangle on the current axis
-%% drawArrow - Draw an arrow on the current axis
-%% drawLabels - Draw labels at specified positions
-%% drawShape - Draw various types of shapes (circles, polygons...)
-%%
-%% Other shapes
-%% squareGrid - Generate equally spaces points in plane.
-%% hexagonalGrid - Generate hexagonal grid of points in the plane.
-%% triangleGrid - Generate triangular grid of points in the plane.
-%% crackPattern - Create a (bounded) crack pattern tessellation
-%% crackPattern2 - Create a (bounded) crack pattern tessellation
-%%
-%%
-%% Credits:
-%% * function 'enclosingCircle' rewritten from a file from Yazan Ahed
-%% , available on Matlab File Exchange
-%%
-%% @end deftypefn
-
-function geometry_Contents ()
-
- help('Contents');
-
- %% Deprecated functions
-
- % createMedian - create a median line
- % minDistance - compute minimum distance between a point and a set of points
- % homothecy - create a homothecy as an affine transform
- % rotation - return 3*3 matrix of a rotation
- % translation - return 3*3 matrix of a translation
- % scaling - return 3*3 matrix of a scaling in 2 dimensions
- % lineSymmetry - create line symmetry as 2D affine transform
- % vecnorm - compute norm of vector or of set of vectors
- % normalize - normalize a vector
- % onCircle - test if a point is located on a given circle.
- % inCircle - test if a point is located inside a given circle.
- % onEdge - test if a point belongs to an edge
- % onLine - test if a point belongs to a line
- % onRay - test if a point belongs to a ray
- % invertLine - return same line but with opposite orientation
- % clipLineRect - clip a line with a polygon
- % formatAngle - Ensure an angle value is comprised between 0 and 2*PI
-
-
- %% Others...
- % drawRect2 - Draw centered rectangle on the current axis
-
-endfunction
-
diff --git a/inst/io/data2geo.m b/inst/io/data2geo.m index 94ff611..d349a5e 100644 --- a/inst/io/data2geo.m +++ b/inst/io/data2geo.m @@ -16,7 +16,7 @@ %% -*- texinfo -*- %% @deftypefn {Function File} {@var{fileStr} =} data2geo (@var{data}, @var{lc}) %% @deftypefnx {Function File} {@var{fileStr} =} data2geo (@dots{}, @var{param}, @var{value}) -%% Builds a file compatible with gmsh form data. +%% Uses data to build a file compatible with Gmsh. %% %% @var{data} is assumed to describe a polygon in @code{polygon2d} format. %% The argument @var{lc} specifies the edge size. @@ -97,6 +97,7 @@ endfunction %! ids = dc.pathid(); %! P = dc.path2polygon(ids{1},12)(1:end-1,:); %! P = bsxfun(@minus, P, centroid(P)); +%! P = simplifypolygon(P,'tol',5e-1); %! filename = tmpnam (); %! meshsize = sqrt(mean(sumsq(diff(P,1,1),2)))/2; %! data2geo (P, meshsize, 'output', [filename '.geo']); diff --git a/inst/io/private/SVGstrPath2SVGpath.m b/inst/io/private/SVGstrPath2SVGpath.m deleted file mode 100644 index 71c0dea..0000000 --- a/inst/io/private/SVGstrPath2SVGpath.m +++ /dev/null @@ -1,103 +0,0 @@ -%% Copyright (c) 2011 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/>. - -function SVGpath = SVGstrPath2SVGpath (SVGstrPath) - - nPaths = numel (SVGstrPath); - SVGpath = repmat (struct('coord', [], 'closed', [], 'id', []), 1, nPaths); - - for ip = 1:nPaths - path = SVGstrPath{ip}; - - % Match data - [s e te m] = regexpi (path, 'd="(?:(?!").)*'); - data=strtrim (m{1}); - % parse data - d = parsePathData (data); - SVGpath(ip).coord = d.coord; - SVGpath(ip).closed = d.closed; - - % Match id - [s e te m] = regexp (path, 'id="(?:(?!").)*'); - if ~isempty (m) - SVGpath(ip).id = strtrim (m{1}(5:end)); - end - end - -end - -function d = parsePathData (data) - - d = struct ('coord', [], 'closed', []); - - [s e te comm] = regexp (data, '[MmLlHhVvCcSsQqTtAaZz]{1}'); - % TODO - % This info could be used to preallocate d - [zcomm zpos] = ismember ({'Z','z'}, comm); - - if any (zcomm) - d.closed = true; - else - d.closed = false; - s(end+1) = length (data); - end - comm(zpos(zcomm)) = []; - ncomm = size (comm, 2); - for ic = 1:ncomm - - switch comm{ic} - case {'M','L'} - [x y] = strread (data(s(ic) + 1 : s(ic + 1) - 1), ... - '%f%f', 'delimiter', ', '); - coord = [x y]; - - case 'm' - [x y] = strread( data(s(ic) + 1 : s(ic + 1) - 1), ... - '%f%f', 'delimiter', ', '); - nc = size (x, 1); - coord = [x y]; - if ic == 1 - % relative moveto at begining of data. - % First coordinates are absolute, the rest are relative. - coord = cumsum (coord); - else - % Relative moveto. - % The coordinates are relative to the last one loaded - coord(1,:) =coord(1,:) + d.coord(end,:); - coord = cumsum(coord); - warning('svg2octWarning',['Relative moveto in path data.'... - ' May mean that the orginal path was not a simple polygon.' ... - ' If that is the case, the path will not display equally.']) - end - - case 'l' - % Relative lineto, coordinates are relative to last point loaded. - [x y] = strread( data(s(ic) + 1 : s(ic + 1) - 1), ... - '%f%f', 'delimiter', ', '); - nc = size (x, 1); - coord = [x y] - coord(1,:) =coord(1,:) + d.coord(end,:); - coord = cumsum(coord); - - otherwise - warning('svg2oct:Warning',... - 'Path data command "%s" not implemented yet.',comm{ic}); - end - - nc = size(coord,1); - d.coord(end+1:end+nc,:) = coord; - - end -end diff --git a/inst/io/private/_parsePath.py b/inst/io/private/_parsePath.py deleted file mode 100644 index a6738ae..0000000 --- a/inst/io/private/_parsePath.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python - -import inkex, simplepath -import sys -#import getopt - -def parsePaths (filen=None): - - svg = inkex.Effect () - svg.parse (filen) - - paths = svg.document.xpath ('//svg:path', namespaces=inkex.NSS) - for path in paths: - D = simplepath.parsePath (path.attrib['d']) - cmdlst = []; - parlst = []; - for cmd,params in D: - cmdlst.append(cmd) - parlst.append(params) - - print 'svgpath = struct("cmd","{0}","data",{{{1}}});' \ - .format(''.join(cmdlst),str(parlst).replace('[[','[').replace(']]',']')) - - print 'svgpathid = "{0}"; $'.format(path.attrib['id']) - - - -# ---------------------------- - -if __name__=="__main__": - ''' - try: - optlist,args = getopt.getopt(sys.argv[1:],"thdp") - except getopt.GetoptError: - usage() - sys.exit(2) - - doHelp = 0 - c = Context() - c.doPrint = 1 - for opt in optlist: - if opt[0] == "-d": c.debug = 1 - if opt[0] == "-p": c.plot = 1 - if opt[0] == "-t": c.triangulate = 1 - if opt[0] == "-h": doHelp = 1 - - if not doHelp: - pts = [] - fp = sys.stdin - if len(args) > 0: - fp = open(args[0],'r') - for line in fp: - fld = line.split() - x = float(fld[0]) - y = float(fld[1]) - pts.append(Site(x,y)) - if len(args) > 0: fp.close() - - if doHelp or len(pts) == 0: - usage() - sys.exit(2) - ''' - svg = sys.argv[1] - parsePaths(svg) diff --git a/inst/io/private/formatSVGstr.m b/inst/io/private/formatSVGstr.m deleted file mode 100644 index a1310e6..0000000 --- a/inst/io/private/formatSVGstr.m +++ /dev/null @@ -1,24 +0,0 @@ -%% Copyright (c) 2011 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/>. - -function svgF = formatSVGstr(svg) - - % Remove all newlines and tabs - svgF = strrep(svg,"\n",' '); - svgF = strrep(svgF,"\t",' '); - % Remove consecutive blanks - svgF = regexprep(svgF,' +',' '); - -end diff --git a/inst/io/private/getSVGPaths_py.m b/inst/io/private/getSVGPaths_py.m deleted file mode 100644 index cca08ee..0000000 --- a/inst/io/private/getSVGPaths_py.m +++ /dev/null @@ -1,113 +0,0 @@ -%% Copyright (c) 2011 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/>. - -function Paths = getSVGPaths_py (svg, varargin) - - %% Call python script - if exist (svg,'file') - % read from file - [st str]=system (sprintf ('python parsePath.py %s', svg)); - - else - % inline SVG - [st str]=system (sprintf ('python parsePath.py < %s', svg)); - end - - %% Parse ouput - strpath = strsplit (str(1:end-1), '$', true); - - npaths = numel (strpath); - - %% Convert path data to polygons - for ip = 1:npaths - - eval (strpath{ip}); - %% FIXME: intialize struct with cell field - svgpath2.cmd = svgpath(1).cmd; - svgpath2.data = {svgpath.data}; - - nD = length(svgpath2.cmd); - pathdata = cell (nD-1,1); - - point_end=[]; - %% If the path is closed, last command is Z and we set initial point == final - if svgpath2.cmd(end) == 'Z' - nD -= 1; - point_end = svgpath2.data{1}; - end - - %% Initial point - points(1,:) = svgpath2.data{1}; - - for jp = 2:nD - switch svgpath2.cmd(jp) - case 'L' - %% Straigth segment to polygon - points(2,:) = svgpath2.data{jp}; - pp = [(points(2,:)-points(1,:))' points(1,:)']; - clear points - points(1,:) = [polyval(pp(1,:),1) polyval(pp(2,:),1)]; - - case 'C' - %% Cubic bezier to polygon - points(2:4,:) = reshape (svgpath2.data{jp}, 2, 3).'; - pp = cbezier2poly (points); - clear points - points(1,:) = [polyval(pp(1,:),1) polyval(pp(2,:),1)]; - end - - pathdata{jp-1} = pp; - end - - if ~isempty(point_end) - %% Straight segmet to close the path - points(2,:) = point_end; - pp = [(points(2,:)-points(1,:))' points(1,:)']; - - pathdata{end} = pp; - end - - Paths.(svgpathid).data = pathdata; - end -endfunction - -%!test -%! figure(1) -%! hold on -%! paths = getSVGPaths_py ('../drawing.svg'); -%! -%! % Get path ids -%! ids = fieldnames(paths); -%! npath = numel(ids); -%! -%! t = linspace (0, 1, 64); -%! -%! for i = 1:npath -%! x = []; y = []; -%! data = paths.(ids(i)).data; -%! -%! for j = 1:numel(data) -%! x = cat (2, x, polyval (data{j}(1,:),t)); -%! y = cat (2, y, polyval (data{j}(2,:),t)); -%! end -%! -%! plot(x,y,'-'); -%! end -%! axis ij -%! if strcmpi(input('You should see drawing.svg [y/n] ','s'),'n') -%! error ("didn't get what was expected."); -%! end -%! close - diff --git a/inst/io/private/getSVGdata.m b/inst/io/private/getSVGdata.m deleted file mode 100644 index 21a1dd8..0000000 --- a/inst/io/private/getSVGdata.m +++ /dev/null @@ -1,33 +0,0 @@ -%% Copyright (c) 2011 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/>. - -function svgData = getSVGdata(svg) - - svgData = struct('height',[],'width',[]); - attr = fieldnames(svgData); - nattr = numel(attr); - - [s e te data] = regexpi(svg,'<[ ]*svg(?:(?!>).)*>'); - data=strtrim(data{1}); - - for a = 1:nattr - pattr =sprintf('%s="(?:(?!").)*',attr{a}); - [s e te m] = regexpi(data,pattr); - m=strtrim(m{1}); - [dummy value] = strread(m,'%s%f','delimiter','"'); - svgData.(attr{a}) = value; - end - -end diff --git a/inst/io/private/parsePath.py b/inst/io/private/parsePath.py deleted file mode 100644 index a6738ae..0000000 --- a/inst/io/private/parsePath.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python - -import inkex, simplepath -import sys -#import getopt - -def parsePaths (filen=None): - - svg = inkex.Effect () - svg.parse (filen) - - paths = svg.document.xpath ('//svg:path', namespaces=inkex.NSS) - for path in paths: - D = simplepath.parsePath (path.attrib['d']) - cmdlst = []; - parlst = []; - for cmd,params in D: - cmdlst.append(cmd) - parlst.append(params) - - print 'svgpath = struct("cmd","{0}","data",{{{1}}});' \ - .format(''.join(cmdlst),str(parlst).replace('[[','[').replace(']]',']')) - - print 'svgpathid = "{0}"; $'.format(path.attrib['id']) - - - -# ---------------------------- - -if __name__=="__main__": - ''' - try: - optlist,args = getopt.getopt(sys.argv[1:],"thdp") - except getopt.GetoptError: - usage() - sys.exit(2) - - doHelp = 0 - c = Context() - c.doPrint = 1 - for opt in optlist: - if opt[0] == "-d": c.debug = 1 - if opt[0] == "-p": c.plot = 1 - if opt[0] == "-t": c.triangulate = 1 - if opt[0] == "-h": doHelp = 1 - - if not doHelp: - pts = [] - fp = sys.stdin - if len(args) > 0: - fp = open(args[0],'r') - for line in fp: - fld = line.split() - x = float(fld[0]) - y = float(fld[1]) - pts.append(Site(x,y)) - if len(args) > 0: fp.close() - - if doHelp or len(pts) == 0: - usage() - sys.exit(2) - ''' - svg = sys.argv[1] - parsePaths(svg) diff --git a/inst/io/svgload.m b/inst/io/svgload.m deleted file mode 100644 index a9e9963..0000000 --- a/inst/io/svgload.m +++ /dev/null @@ -1,62 +0,0 @@ -%% Copyright (c) 2011 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{SVG} = loadSVG (@var{filename}) -%% Reads the plain SVG file @var{filename} and returns an SVG structure. -%% -%% In the current version only SVG path elements are parsed and stored in the field -%% path of the @var{SVG} structure. -%% The coordinates of the path are not modified in anyway. This produces the path -%% to be ploted vertically reflected. -%% -%% @seealso{svgnormalize, svgpath2polygon} -%% @end deftypefn - - -function SVG = svgload (filename) - - SVG = struct ('height', [], 'width', [], 'path', [], 'normalized', []); - - SVG.normalized = false; - - fid = fopen (filename); - svg = char (fread (fid, "uchar")'); - fclose (fid); - svgF = formatSVGstr (svg); - - % Get SVG Data - data = getSVGdata (svgF); - SVG.height = data.height; - SVG.width = data.width; - - % Get SVG Paths - SVGstrPath = getSVGstrPath (svgF); - SVGpath = SVGstrPath2SVGpath (SVGstrPath); - SVG.path = SVGpath; - -end - -%!demo -%! file = 'tmp__.svg'; -%! fid = fopen (file,'w'); -%! svgfile = '<html><body><svg xmlns="http://www.w3.org/2000/svg" version="1.1" height="250" width="250"><path d="M150,0 75,200 225,200 Z" /></svg></body></html>'; -%! fprintf (fid,"%s\n",svgfile); -%! fclose (fid); -%! SVG = svgload (file); -%! SVG -%! plot([SVG.path.coord(:,1); SVG.path.coord(1,1)], ... -%! [SVG.path.coord(:,2); SVG.path.coord(1,2)]); -%! delete (file); diff --git a/inst/io/svgnormalize.m b/inst/io/svgnormalize.m deleted file mode 100644 index 381d29a..0000000 --- a/inst/io/svgnormalize.m +++ /dev/null @@ -1,68 +0,0 @@ -%% Copyright (c) 2011 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{SVGn} = loadSVG (@var{SVG}) -%% Scales and reflects the @var{SVG} structure and returns a modified @var{SVGn} -%% structure. -%% -%% The height and width of the SVG are scaled such that the diagonal of the -%% bounding box has length 1. Coordinates are trasnfomed such that a plot of the -%% paths coincides with the visualization of the original SVG. -%% -%% @seealso{svgload, svgpath2polygon} -%% @end deftypefn - -function SVGn = svgnormalize (SVG) - - SVGn = SVG; - if ~SVG.normalized - % Translate - TransV = [0, SVG.height/2]; - % Scale - Dnorm = sqrt (SVG.width ^ 2 + SVG.height ^ 2); - S = (1 / Dnorm) * eye (2); - for i = 1:numel (SVG.path) - nc = size (SVG.path(i).coord, 1); - % Translate to middle - T = repmat (TransV,nc, 1); - SVGn.path(i).coord = SVG.path(i).coord - T; - % Reflect - SVGn.path(i).coord(:, 2) = -SVGn.path(i).coord(:,2); - T(:,2) = -T(:,2); - % Translate to bottom - SVGn.path(i).coord = SVGn.path(i).coord - T; - %Scale - SVGn.path(i).coord = SVGn.path(i).coord * S; - end - SVGn.height = SVG.height / Dnorm; - SVGn.width = SVG.width / Dnorm; - SVGn.normalized = true; - end - -end - -%!demo -%! file = 'tmp__.svg'; -%! fid = fopen (file,'w'); -%! svgfile = '<html><body><svg xmlns="http://www.w3.org/2000/svg" version="1.1" height="250" width="250"><path d="M150,0 75,200 225,200 Z" /></svg></body></html>'; -%! fprintf (fid,"%s\n",svgfile); -%! fclose (fid); -%! SVG = svgload (file); -%! SVGn = svgnormalize (SVG); -%! SVGn -%! plot([SVGn.path.coord(:,1); SVGn.path.coord(1,1)], ... -%! [SVGn.path.coord(:,2); SVGn.path.coord(1,2)]); -%! delete (file); diff --git a/inst/io/svgpath2polygon.m b/inst/io/svgpath2polygon.m deleted file mode 100644 index e03a78f..0000000 --- a/inst/io/svgpath2polygon.m +++ /dev/null @@ -1,60 +0,0 @@ -%% Copyright (c) 2011 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{P} = svgpath2polygon (@var{SVGpath}) -%% Converts the SVG path structure @var{SVGpath} to an array of polygons -%% compatible with the geometry package and matGeom (@url{http://matgeom.sf.net}). -%% -%% @var{SVGpath} is a substructure of the SVG structure output by loadSVG. This -%% function extracts the field named "coord" if there is only one path. If there -%% are more than oe path it puts the "coord" field of each path in the same -%% array, separated by nans. -%% -%% @seealso{svgnormalize, svgload} -%% @end deftypefn - -function P = svgpath2polygon (SVGpath) - - P = SVGpath(1).coord; - for ip = 2:numel (SVGpath) - P = [P; nan(1,2); SVGpath(ip).coord]; - end - -end - -%!demo -%! file = 'tmp__.svg'; -%! fid = fopen (file,'w'); -%! svgfile = '<html><body><svg xmlns="http://www.w3.org/2000/svg" version="1.1" height="250" width="250"><path d="M150,0 75,200 225,200 Z" /></svg></body></html>'; -%! fprintf (fid,"%s\n",svgfile); -%! fclose (fid); -%! SVG = svgload (file); -%! SVG = svgnormalize (SVG); -%! P = svgpath2polygon (SVG.path); -%! plot (P(:,1),P(:,2)); -%! delete (file); - -%!demo -%! file = 'tmp__.svg'; -%! fid = fopen (file,'w'); -%! svgfile = '<html><body><svg xmlns="http://www.w3.org/2000/svg" version="1.1" height="250" width="250"><path d="M150,0 75,200 225,200 Z" /><path d="M0,0 100,0 100,100 0,100 Z" /></svg></body></html>'; -%! fprintf (fid,"%s\n",svgfile); -%! fclose (fid); -%! SVG = svgload (file); -%! SVG = svgnormalize (SVG); -%! P = svgpath2polygon (SVG.path); -%! plot (P(:,1),P(:,2)); -%! delete (file); diff --git a/inst/octclip/oc_polybool.m b/inst/octclip/oc_polybool.m index 05c9f96..259f057 100644 --- a/inst/octclip/oc_polybool.m +++ b/inst/octclip/oc_polybool.m @@ -45,7 +45,7 @@ ## ## For the matrices @var{sub} and @var{clip}, the first point is not needed to ## be repeated at the end (but is permitted). Pairs of (NaN,NaN) coordinates in -## @var{sub} and/or @var{clip} are ommitted. +## @var{sub} and/or @var{clip} are omitted. ## ## @var{X} is a column vector containing the X coordinates of the vertices for. ## resultant polygon(s). diff --git a/inst/octclip/src/_oc_polybool.cc b/inst/octclip/src/_oc_polybool.cc index b40e7e5..db696a5 100644 --- a/inst/octclip/src/_oc_polybool.cc +++ b/inst/octclip/src/_oc_polybool.cc @@ -49,7 +49,7 @@ Operation of @var{clip} - @var{sub}.\n\ \n\ For the matrices @var{sub} and @var{clip}, the first point is not needed to\n\ be repeated at the end (but is permitted). Pairs of (NaN,NaN) coordinates in\n\ -@var{sub} and/or @var{clip} are ommitted.\n\ +@var{sub} and/or @var{clip} are omitted.\n\ \n\ @var{X} is a column vector containing the X coordinates of the vertices for.\n\ resultant polygon(s).\n\n\ diff --git a/inst/polygons2d/curvature.m b/inst/polygons2d/curvature.m new file mode 100644 index 0000000..e14b14c --- /dev/null +++ b/inst/polygons2d/curvature.m @@ -0,0 +1,177 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{kappa} = } curvature (@var{t}, @var{px}, @var{py},@var{method},@var{degree}) +%% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{t}, @var{poly},@var{method},@var{degree}) +%% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{px}, @var{py},@var{method},@var{degree}) +%% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{points},@var{method},@var{degree}) +%% @deftypefnx {Function File} {[@var{kappa} @var{poly} @var{t}] = } curvature (@dots{}) +%% Estimate curvature of a polyline defined by points. +%% +%% First compute an approximation of the curve given by PX and PY, with +%% the parametrization @var{t}. Then compute the curvature of approximated curve +%% for each point. +%% @var{method} used for approximation can be only: 'polynom', with specified degree. +%% Further methods will be provided in a future version. +%% @var{t}, @var{px}, and @var{py} are N-by-1 array of the same length. The points +%% can be specified as a single N-by-2 array. +%% +%% If the argument @var{t} is not given, the parametrization is estimated using +%% function @code{parametrize}. +%% +%% If requested, @var{poly} contains the approximating polygon evlauted at the +%% parametrization @var{t}. +%% +%% @seealso{parametrize, polygons2d} +%% @end deftypefn + +function [kappa, varargout] = curvature(varargin) + + % default values + degree = 5; + t=0; % parametrization of curve + tc=0; % indices of points wished for curvature + + + % ================================================================= + + % Extract method and degree ------------------------------ + + nargin = length(varargin); + varN = varargin{nargin}; + varN2 = varargin{nargin-1}; + + if ischar(varN2) + % method and degree are specified + method = varN2; + degree = varN; + varargin = varargin(1:nargin-2); + elseif ischar(varN) + % only method is specified, use degree 6 as default + method = varN; + varargin = varargin{1:nargin-1}; + else + % method and degree are implicit : use 'polynom' and 6 + method = 'polynom'; + end + + % extract input parametrization and curve. ----------------------- + nargin = length(varargin); + if nargin==1 + % parameters are just the points -> compute caracterization. + var = varargin{1}; + px = var(:,1); + py = var(:,2); + elseif nargin==2 + var = varargin{2}; + if size(var, 2)==2 + % parameters are t and POINTS + px = var(:,1); + py = var(:,2); + t = varargin{1}; + else + % parameters are px and py + px = varargin{1}; + py = var; + end + elseif nargin==3 + var = varargin{2}; + if size(var, 2)==2 + % parameters are t, POINTS, and tc + px = var(:,1); + py = var(:,2); + t = varargin{1}; + else + % parameters are t, px and py + t = varargin{1}; + px = var; + py = varargin{3}; + end + elseif nargin==4 + % parameters are t, px, py and tc + t = varargin{1}; + px = varargin{2}; + py = varargin{3}; + tc = varargin{4}; + end + + % compute implicit parameters -------------------------- + + % if t and/or tc are not computed, use implicit definition + if t==0 + t = parametrize(px, py, 'norm'); + end + + % if tc not defined, compute curvature for all points + if tc==0 + tc = t; + else + % else convert from indices to parametrization values + tc = t(tc); + end + + + % ================================================================= + % compute curvature for each point of the curve + + if strcmp(method, 'polynom') + % compute coefficients of interpolation functions + x0 = polyfit(t, px, degree); + y0 = polyfit(t, py, degree); + + % compute coefficients of first and second derivatives. In the case of a + % polynom, it is possible to compute coefficient of derivative by + % multiplying with a matrix. + derive = diag(degree:-1:0); + xp = circshift(x0*derive, [0 1]); + yp = circshift(y0*derive, [0 1]); + xs = circshift(xp*derive, [0 1]); + ys = circshift(yp*derive, [0 1]); + + % compute values of first and second derivatives for needed points + xprime = polyval(xp, tc); + yprime = polyval(yp, tc); + xsec = polyval(xs, tc); + ysec = polyval(ys, tc); + + % compute value of curvature + kappa = (xprime.*ysec - xsec.*yprime)./ ... + power(xprime.*xprime + yprime.*yprime, 3/2); + + if nargout > 1 + varargout{1} = [polyval(x0,tc(:)) polyval(y0,tc(:))]; + if nargout > 2 + varargout{2} = tc; + end + end + else + error('unknown method'); + end + +endfunction diff --git a/inst/polygons2d/distancePointPolygon.m b/inst/polygons2d/distancePointPolygon.m new file mode 100755 index 0000000..6290b91 --- /dev/null +++ b/inst/polygons2d/distancePointPolygon.m @@ -0,0 +1,51 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{dist} = } distancePointPolygon (@var{point},@var{poly}) +## Compute shortest distance between a point and a polygon +## +## @seealso{polygons2d, points2d, distancePointPolyline, distancePointEdge, projPointOnPolyline} +## @end deftypefn + +function varargout = distancePointPolygon(point, poly) + + % eventually copy first point at the end to ensure closed polygon + if sum(poly(end, :)==poly(1,:))~=2 + poly = [poly; poly(1,:)]; + end + + % call to distancePointPolyline + minDist = distancePointPolyline(point, poly); + + % process output arguments + if nargout<=1 + varargout{1} = minDist; + end + +endfunction diff --git a/inst/polygons2d/distancePointPolyline.m b/inst/polygons2d/distancePointPolyline.m new file mode 100755 index 0000000..c4fcd4a --- /dev/null +++ b/inst/polygons2d/distancePointPolyline.m @@ -0,0 +1,71 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{dist} = } distancePointPolyline (@var{point},@var{poly}) +## Compute shortest distance between a point and a polyline +## Example: +## @example +## pt1 = [30 20]; +## pt2 = [30 5]; +## poly = [10 10;50 10;50 50;10 50]; +## distancePointPolyline([pt1;pt2], poly) +## ans = +## 10 +## 5 +## @end example +## +## @seealso{polygons2d, points2d,distancePointEdge, projPointOnPolyline} +## @end deftypefn + +function varargout = distancePointPolyline(point, poly, varargin) + + # number of points + Np = size(point, 1); + + # allocate memory for result + minDist = inf * ones(Np, 1); + + # process each point + for p = 1:Np + # construct the set of edges + edges = [poly(1:end-1, :) poly(2:end, :)]; + + # compute distance between current each point and all edges + dist = distancePointEdge(point(p, :), edges); + + # update distance if necessary + minDist(p) = min(dist); + end + + # process output arguments + if nargout<=1 + varargout{1} = minDist; + end + +endfunction diff --git a/inst/polygons2d/distancePolygons.m b/inst/polygons2d/distancePolygons.m new file mode 100755 index 0000000..53dfc9b --- /dev/null +++ b/inst/polygons2d/distancePolygons.m @@ -0,0 +1,43 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{dist} = } distancePolygons (@var{poly1},@var{poly2}) +## Compute the shortest distance between 2 polygons +## +## @end deftypefn +function dist = distancePolygons(poly1, poly2) + + # compute distance of each vertex of a polygon to the other polygon + dist1 = min(distancePointPolygon(poly1, poly2)); + dist2 = min(distancePointPolygon(poly2, poly1)); + + # keep the minimum of the two distances + dist = min(dist1, dist2); + +endfunction diff --git a/inst/polygons2d/expandPolygon.m b/inst/polygons2d/expandPolygon.m new file mode 100755 index 0000000..d009454 --- /dev/null +++ b/inst/polygons2d/expandPolygon.m @@ -0,0 +1,85 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{loops} = } expandPolygon (@var{poly}, @var{dist}) +## Expand a polygon by a given (signed) distance +## +## Associates to each edge of the polygon @var{poly} the parallel line located +## at distance @var{dist} from the current edge, and compute intersections with +## neighbor parallel lines. The resulting polygon is simplified to remove +## inner "loops", and can be disconnected. +## The result is a cell array, each cell containing a simple linear ring. +## +## This is a kind of dilation, but behaviour on corners is different. +## This function keeps angles of polygons, but there is no direct relation +## between length of 2 polygons. +## +## It is also possible to specify negative distance, and get all points +## inside the polygon. If the polygon is convex, the result equals +## morphological erosion of polygon by a ball with radius equal to the +## given distance. +## +## @seealso{polygons2d} +## @end deftypefn + +function loops = expandPolygon(poly, dist) + + # eventually copy first point at the end to ensure closed polygon + if sum(poly(end, :)==poly(1,:))~=2 + poly = [poly; poly(1,:)]; + end + + # number of vertices of the polygon + N = size(poly, 1)-1; + + # find lines parallel to polygon edges located at distance DIST + lines = zeros(N, 4); + for i=1:N + side = createLine(poly(i,:), poly(i+1,:)); + lines(i, 1:4) = parallelLine(side, dist); + end + + # compute intersection points of consecutive lines + lines = [lines;lines(1,:)]; + poly2 = zeros(N, 2); + for i=1:N + poly2(i,1:2) = intersectLines(lines(i,:), lines(i+1,:)); + end + + # split result polygon into set of loops (simple polygons) + loops = polygonLoops(poly2); + + # keep only loops whose distance to original polygon is correct + distLoop = zeros(length(loops), 1); + for i=1:length(loops) + distLoop(i) = distancePolygons(loops{i}, poly); + end + loops = loops(abs(distLoop-abs(dist))<abs(dist)/1000); + +endfunction diff --git a/inst/polygons2d/medialAxisConvex.m b/inst/polygons2d/medialAxisConvex.m new file mode 100755 index 0000000..ac7a1fa --- /dev/null +++ b/inst/polygons2d/medialAxisConvex.m @@ -0,0 +1,153 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{n} @var{e}] = } medialAxisConvex (@var{polygon}) +## Compute medial axis of a convex polygon +## +## @var{polygon} is given as a set of points [x1 y1;x2 y2 ...], returns +## the medial axis of the polygon as a graph. +## @var{n} is a set of nodes. The first elements of @var{n} are the vertices of the +## original polygon. +## @var{e} is a set of edges, containing indices of source and target nodes. +## Edges are sorted according to order of creation. Index of first vertex +## is lower than index of last vertex, i.e. edges always point to newly +## created nodes. +## +## Notes: +## - Is not fully implemented, need more development (usually crashes for +## polygons with more than 6-7 points...) +## - Works only for convex polygons. +## - Complexity is not optimal: this algorithm is O(n*log n), but linear +## algorithms exist. +## +## @seealso{polygons2d, bisector} +## @end deftypefn + +function [nodes, edges] = medialAxisConvex(points) + + # eventually remove the last point if it is the same as the first one + if points(1,:) == points(end, :) + nodes = points(1:end-1, :); + else + nodes = points; + end + + # special case of triangles: + # compute directly the gravity center, and simplify computation. + if size(nodes, 1)==3 + nodes = [nodes; mean(nodes, 1)]; + edges = [1 4;2 4;3 4]; + return + end + + # number of nodes, and also of initial rays + N = size(nodes, 1); + + # create ray of each vertex + rays = zeros(N, 4); + rays(1, 1:4) = bisector(nodes([2 1 N], :)); + rays(N, 1:4) = bisector(nodes([1 N N-1], :)); + for i=2:N-1 + rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :)); + end + + # add indices of edges producing rays (indices of first vertex, second + # vertex is obtained by adding one modulo N). + rayEdges = [[N (1:N-1)]' (1:N)']; + + pint = intersectLines(rays, rays([2:N 1], :)); + #ti = linePosition(pint, rays); + #ti = min(linePosition(pint, rays), linePosition(pint, rays([2:N 1], :))); + ti = distancePointLine(pint, ... + createLine(points([N (1:N-1)]', :), points((1:N)', :))); + + # create list of events. + # terms are : R1 R2 X Y t0 + # R1 and R2 are indices of involved rays + # X and Y is coordinate of intersection point + # t0 is position of point on rays + events = sortrows([ (1:N)' [2:N 1]' pint ti], 5); + + # initialize edges + edges = zeros(0, 2); + + + # ------------------- + # process each event until there is no more + + # start after index of last vertex, and process N-3 intermediate rays + for i=N+1:2*N-3 + # add new node at the rays intersection + nodes(i,:) = events(1, 3:4); + + # add new couple of edges + edges = [edges; events(1,1) i; events(1,2) i]; + + # find the two edges creating the new emanating ray + n1 = rayEdges(events(1, 1), 1); + n2 = rayEdges(events(1, 2), 2); + + # create the new ray + line1 = createLine(nodes(n1, :), nodes(mod(n1,N)+1, :)); + line2 = createLine(nodes(mod(n2,N)+1, :), nodes(n2, :)); + ray0 = bisector(line1, line2); + + # set its origin to emanating point + ray0(1:2) = nodes(i, :); + + # add the new ray to the list + rays = [rays; ray0]; + rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2]; + + # find the two neighbour rays + ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0; + ir = unique(events(ind, 1:2)); + ir = ir(~ismember(ir, events(1,1:2))); + + # create new intersections + pint = intersectLines(ray0, rays(ir, :)); + #ti = min(linePosition(pint, ray0), linePosition(pint, rays(ir, :))) + events(1,5); + ti = distancePointLine(pint, line1); + + # remove all events involving old intersected rays + ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0; + events = events(ind, :); + + # add the newly formed events + events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)]; + + # and sort them according to 'position' parameter + events = sortrows(events, 5); + end + + # centroid computation for last 3 rays + nodes = [nodes; mean(events(:, 3:4))]; + edges = [edges; [unique(events(:,1:2)) ones(3, 1)*(2*N-2)]]; + +endfunction diff --git a/inst/polygons2d/parametrize.m b/inst/polygons2d/parametrize.m new file mode 100644 index 0000000..0e36196 --- /dev/null +++ b/inst/polygons2d/parametrize.m @@ -0,0 +1,96 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{par} = } parametrize (@var{poly}) +%% @deftypefnx {Function File} {@var{par} = } parametrize (@var{px},@var{py}) +%% @deftypefnx {Function File} {@var{par} = } parametrize (@dots{},@var{normalize}) +%% +%% Parametrization of a curve, based on edges length +%% +%% Returns a parametrization of the curve defined by the serie of points, +%% based on euclidean distance between two consecutive points. +%% POLY is a N-by-2 array, representing coordinates of vertices. The +%% result PAR is N-by-1, and contains the cumulative length of edges until +%% corresponding vertex. The function also accepts the polygon as two inputs +%% @var{px} and @var{py} representinx coordinates x and y respectively. +%% When optional argument @var{normalize} is non-empty @var{par} is rescaled +%% such that the last element equals 1, i.e. @code{@var{par}(end)==1}. +%% +%% Example +%% @example +%% % Parametrize a circle approximation +%% poly = circleToPolygon([0 0 1], 200); +%% p = parametrize(poly); +%% p(end) +%% ans = +%% 6.2829 +%% p = parametrize(poly,'norm'); +%% p(end) +%% ans = +%% 1 +%% p = parametrize(poly,true); +%% p(end) +%% ans = +%% 1 +%% @end example +%% +%% @seealso{polygons2d, polylineLength} +%% @end deftypefn +function par = parametrize(varargin) + %% Process inputs + + % extract vertex coordinates + if size(varargin{1}, 2) > 1 + % vertices in a single array + pts = varargin{1}; + normalize = numel(varargin) > 1; + elseif size(varargin{1}, 2) == 1 && numel(varargin) >= 2 + % points as separate arrays + pts = [varargin{1} varargin{2}]; + normalize = numel(varargin) > 2; + end + + %% Parametrize polyline + + % compute cumulative sum of euclidean distances between consecutive + % vertices, setting distance of first vertex to 0. + if size(pts, 2) == 2 + % process points in 2D + par = [0 ; cumsum(hypot(diff(pts(:,1)), diff(pts(:,2))))]; + else + % process points in arbitrary dimension + par = [0 ; cumsum(sqrt(sum(diff(pts).^2, 2)))]; + end + + % eventually rescale between 0 and 1 + if normalize + par = par / par(end); + end + +endfunction diff --git a/inst/polygons2d/polygonLoops.m b/inst/polygons2d/polygonLoops.m new file mode 100755 index 0000000..e024f6b --- /dev/null +++ b/inst/polygons2d/polygonLoops.m @@ -0,0 +1,143 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{loops} = } polygonLoops (@var{poly}) +## Divide a possibly self-intersecting polygon into a set of simple loops +## +## @var{poly} is a polygone defined by a series of vertices, +## @var{loops} is a cell array of polygons, containing the same vertices of the +## original polygon, but no loop self-intersect, and no couple of loops +## intersect each other. +## +## Example: +## @example +## poly = [0 0;0 10;20 10;20 20;10 20;10 0]; +## loops = polygonLoops(poly); +## figure(1); hold on; +## drawPolygon(loops); +## polygonArea(loops) +## @end example +## +## @seealso{polygons2d, polygonSelfIntersections} +## @end deftypefn + +function loops = polygonLoops(poly) + + ## Initialisations + + # compute intersections + [inters pos1 pos2] = polygonSelfIntersections(poly); + + # case of a polygon without self-intersection + if isempty(inters) + loops = {poly}; + return; + end + + # array for storing loops + loops = cell(0, 1); + + positions = sortrows([pos1 pos2;pos2 pos1]); + + + ## First loop + + pos0 = 0; + loop = polygonSubcurve(poly, pos0, positions(1, 1)); + pos = positions(1, 2); + positions(1, :) = []; + + while true + # index of next intersection point + ind = find(positions(:,1)>pos, 1, 'first'); + + # if not found, break + if isempty(ind) + break; + end + + # add portion of curve + loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW> + + # look for next intersection point + pos = positions(ind, 2); + positions(ind, :) = []; + end + + # add the last portion of curve + loop = [loop;polygonSubcurve(poly, pos, pos0)]; + + # remove redundant vertices + loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; + if sum(diff(loop([1 end], :))==0)==2 + loop(end, :) = []; + end + + # add current loop to the list of loops + loops{1} = loop; + + + ## Other loops + + Nl = 1; + while ~isempty(positions) + + loop = []; + pos0 = positions(1, 2); + pos = positions(1, 2); + + while true + # index of next intersection point + ind = find(positions(:,1)>pos, 1, 'first'); + + # add portion of curve + loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW> + + # look for next intersection point + pos = positions(ind, 2); + positions(ind, :) = []; + + # if not found, break + if pos==pos0 + break; + end + end + + # remove redundant vertices + loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; ##ok<AGROW> + if sum(diff(loop([1 end], :))==0)==2 + loop(end, :) = []; ##ok<AGROW> + end + + # add current loop to the list of loops + Nl = Nl + 1; + loops{Nl} = loop; + end + +endfunction diff --git a/inst/polygons2d/polygonSelfIntersections.m b/inst/polygons2d/polygonSelfIntersections.m new file mode 100755 index 0000000..ec5f43f --- /dev/null +++ b/inst/polygons2d/polygonSelfIntersections.m @@ -0,0 +1,87 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pts} = } polygonSelfIntersections (@var{poly}) +## @deftypefnx {Function File} {[@var{pts} @var{pos1} @var{pos2}] = } polygonSelfIntersections (@var{poly}) +## Find-self intersection points of a polygon +## +## Return the position of self intersection points +## +## Also return the 2 positions of each intersection point (the position +## when meeting point for first time, then position when meeting point +## for the second time). +## +## Example +## @example +## # use a '8'-shaped polygon +## poly = [10 0;0 0;0 10;20 10;20 20;10 20]; +## polygonSelfIntersections(poly) +## ans = +## 10 10 +## @end example +## +## @seealso{polygons2d, polylineSelfIntersections} +## @end deftypefn + +function varargout = polygonSelfIntersections(poly, varargin) + + tol = 1e-14; + + # ensure the last point equals the first one + if sum(abs(poly(end, :)-poly(1,:)) < tol) ~= 2 + poly = [poly; poly(1,:)]; + end + + # compute intersections by calling algo for polylines + [points pos1 pos2] = polylineSelfIntersections(poly); + + # It may append that first vertex of polygon is detected as intersection, + # the following tries to detect this + n = size(poly, 1) - 1; + inds = (pos1 == 0 & pos2 == n) | (pos1 == n & pos2 == 0); + points(inds, :) = []; + pos1(inds) = []; + pos2(inds) = []; + + # remove multiple intersections + [points I J] = unique(points, 'rows', 'first'); ##ok<NASGU> + pos1 = pos1(I); + pos2 = pos2(I); + + + ## Post-processing + + # process output arguments + if nargout <= 1 + varargout = {points}; + elseif nargout == 3 + varargout = {points, pos1, pos2}; + end + +endfunction diff --git a/inst/polygons2d/polylineSelfIntersections.m b/inst/polygons2d/polylineSelfIntersections.m new file mode 100755 index 0000000..d759be9 --- /dev/null +++ b/inst/polygons2d/polylineSelfIntersections.m @@ -0,0 +1,153 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pts} = } polylineSelfIntersections (@var{poly}) +## Find self-intersections points of a polyline +## +## Return the position of self intersection points +## +## Also return the 2 positions of each intersection point (the position +## when meeting point for first time, then position when meeting point +## for the second time). +## +## Example +## @example +## # use a gamma-shaped polyline +## poly = [0 0;0 10;20 10;20 20;10 20;10 0]; +## polylineSelfIntersections(poly) +## ans = +## 10 10 +## +## # use a 'S'-shaped polyline +## poly = [10 0;0 0;0 10;20 10;20 20;10 20]; +## polylineSelfIntersections(poly) +## ans = +## 10 10 +## @end example +## +## @seealso{polygons2d, intersectPolylines, polygonSelfIntersections} +## @end deftypefn + +function varargout = polylineSelfIntersections(poly, varargin) + + ## Initialisations + + # determine whether the polyline is closed + closed = false; + if ~isempty(varargin) + closed = varargin{1}; + if ischar(closed) + if strcmp(closed, 'closed') + closed = true; + elseif strcmp(closed, 'open') + closed = false; + end + end + end + + # if polyline is closed, ensure the last point equals the first one + if closed + if sum(abs(poly(end, :)-poly(1,:))<1e-14)~=2 + poly = [poly; poly(1,:)]; + end + end + + # arrays for storing results + points = zeros(0, 2); + pos1 = zeros(0, 1); + pos2 = zeros(0, 1); + + # number of vertices + Nv = size(poly, 1); + + + ## Main processing + + # index of current intersection + ip = 0; + + # iterate over each couple of edge ( (N-1)*(N-2)/2 iterations) + for i=1:Nv-2 + # create first edge + edge1 = [poly(i, :) poly(i+1, :)]; + for j=i+2:Nv-1 + # create second edge + edge2 = [poly(j, :) poly(j+1, :)]; + + # check conditions on bounding boxes + if min(edge1([1 3]))>max(edge2([1 3])) + continue; + end + if max(edge1([1 3]))<min(edge2([1 3])) + continue; + end + if min(edge1([2 4]))>max(edge2([2 4])) + continue; + end + if max(edge1([2 4]))<min(edge2([2 4])) + continue; + end + + # compute intersection point + inter = intersectEdges(edge1, edge2); + + if sum(isfinite(inter))==2 + # add point to the list + ip = ip + 1; + points(ip, :) = inter; + + # also compute positions on the polyline + pos1(ip, 1) = i+edgePosition(inter, edge1)-1; + pos2(ip, 1) = j+edgePosition(inter, edge2)-1; + end + end + end + + # if polyline is closed, the first vertex was found as an intersection, so + # we need to remove it + if closed + dist = distancePoints(points, poly(1,:)); + [minDist ind] = min(dist); ##ok<ASGLU> + points(ind,:) = []; + pos1(ind) = []; + pos2(ind) = []; + end + + ## Post-processing + + # process output arguments + if nargout<=1 + varargout{1} = points; + elseif nargout==3 + varargout{1} = points; + varargout{2} = pos1; + varargout{3} = pos2; + end + +endfunction diff --git a/inst/polygons2d/reversePolygon.m b/inst/polygons2d/reversePolygon.m new file mode 100644 index 0000000..b0febf2 --- /dev/null +++ b/inst/polygons2d/reversePolygon.m @@ -0,0 +1,44 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{poly2} = } reversePolygon (@var{poly}) +%% Reverse a polygon, by iterating vertices from the end +%% +%% POLY2 = reversePolygon(POLY) +%% POLY2 has same vertices as POLY, but in different order. The first +%% vertex of the polygon is still the same. +%% +%% @seealso{reversePolyline} +%% @end deftypefn + +function rev = reversePolygon(poly) + + rev = poly([1 end:-1:2], :); + +endfunction diff --git a/inst/polygons2d/reversePolyline.m b/inst/polygons2d/reversePolyline.m new file mode 100644 index 0000000..bc16338 --- /dev/null +++ b/inst/polygons2d/reversePolyline.m @@ -0,0 +1,43 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{poly2} = } reversePolyline (@var{poly}) +%% Reverse a polyline, by iterating vertices from the end +%% +%% POLY2 = reversePolyline(POLY) +%% POLY2 has same vertices as POLY, but POLY2(i,:) is the same as +%% POLY(END-i+1,:). +%% +%% @seealso{reversePolygon} +%% @end deftypefn +function rev = reversePolyline(poly) + + rev = poly(end:-1:1, :); + +endfunction diff --git a/inst/polygons2d/simplifypolygon.m b/inst/polygons2d/simplifypolygon.m index 4a197a6..5993cfc 100644 --- a/inst/polygons2d/simplifypolygon.m +++ b/inst/polygons2d/simplifypolygon.m @@ -15,37 +15,32 @@ %% -*- texinfo -*- %% @deftypefn {Function File} {@var{spoly} = } simplifypolygon (@var{poly}) -%% Filter colinear vertex from a 2D polygon. +%% Simplify a polygon using the Ramer-Douglas-Peucker algorithm. %% %% @var{poly} is a N-by-2 matrix, each row representing a vertex. %% -%% @seealso{shape2polygon} +%% @seealso{simplifypolyline, shape2polygon} %% @end deftypefn -function polygonsimp = simplifypolygon (polygon) +function polygonsimp = simplifypolygon (polygon, varargin) - # Filter colinear points - edges = diff(polygon([1:end 1],:)); - ned = size(edges,1); - nxt = [2:ned 1]; + polygonsimp = simplifypolyline (polygon,varargin{:}); - # check if consecutive edges are parallel - para = edges(:,1).*edges(nxt,2) - edges(:,2).*edges(nxt,1); - ind = abs(para) > sqrt(eps); - - polygonsimp = polygon(circshift (ind,1),:); - - if isempty(polygonsimp) - warning('simplifypolygon:devel',"The simplification gives an empty polygon. Returning original\n"); - polygonsimp = polygon; - end + # Remove parrallel consecutive edges + PL = polygonsimp(1:end-1,:); + PC = polygonsimp(2:end,:); + PR = polygonsimp([3:end 1],:); + a = PL - PC; + b = PR - PC; + tf = find(isParallel(a,b))+1; + polygonsimp (tf,:) = []; endfunction %!test %! P = [0 0; 1 0; 0 1]; %! P2 = [0 0; 0.1 0; 0.2 0; 0.25 0; 1 0; 0 1; 0 0.7; 0 0.6; 0 0.3; 0 0.1]; -%! assert(P,simplifypolygon (P2)) +%! assert(simplifypolygon (P2),P,min(P2(:))*eps) %!demo %! diff --git a/inst/polygons2d/simplifypolyline.m b/inst/polygons2d/simplifypolyline.m new file mode 100644 index 0000000..0be0efd --- /dev/null +++ b/inst/polygons2d/simplifypolyline.m @@ -0,0 +1,156 @@ +%% 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{pline2} @var{idx}] = } simplifypolyline (@var{pline}) +%% @deftypefnx {Function File} {@dots{} = } simplifypolyline (@dots{},@var{property},@var{value},@dots{}) +%% Simplify or subsample a polyline using the Ramer-Douglas-Peucker algorithm, +%% a.k.a. the iterative end-point fit algorithm or the split-and-merge algorithm. +%% +%% The @var{pline} as a N-by-2 matrix. Rows correspond to the +%% verices (compatible with @code{polygons2d}). The vector @var{idx} constains +%% the indexes on vetices in @var{pline} that generates @var{pline2}, i.e. +%% @code{pline2 = pline(idx,:)}. +%% +%% @strong{Parameters} +%% @table @samp +%% @item 'Nmax' +%% Maximum number of vertices. Default value @code{1e3}. +%% @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 +%% +%% Run @code{demo simplifypolyline} to see an example. +%% +%% @seealso{curve2polyline, curveval} +%% @end deftypefn + +function [pline idx] = simplifypolyline (pline_o, varargin) +%% TODO do not print warnings if user provided Nmax or MaxIter. + + # --- Parse arguments --- # + parser = inputParser (); + parser.FunctionName = "simplifypolyline"; + parser = addParamValue (parser,'Nmax', 1e3, @(x)x>0); + toldef = 1e-4;%max(sumsq(diff(pline_o),2))*2; + parser = addParamValue (parser,'Tol', toldef, @(x)x>0); + parser = addParamValue (parser,'MaxIter', 100, @(x)x>0); + parser = parse(parser,varargin{:}); + + Nmax = parser.Results.Nmax; + tol = parser.Results.Tol; + MaxIter = parser.Results.MaxIter; + + clear parser toldef + msg = ["simplifypolyline: Maximum number of points reached with maximum error %g." ... + " Increase %s if the result is not satisfactory."]; + # ------ # + + [N dim] = size(pline_o); + idx = [1 N]; + + for iter = 1:MaxIter + % Find the point with the maximum distance. + [dist ii] = maxdistance (pline_o, idx); + + tf = dist > tol; + n = sum(tf); + if all(!tf); + break; + end + + idx(end+1:end+n) = ii(tf); + idx = sort(idx); + + if length(idx) >= Nmax + %% TODO remove extra points + warning('geometry:MayBeWrongOutput', sprintf(msg,max(dist),'Nmax')); + break; + end + + end + if iter == MaxIter + warning('geometry:MayBeWrongOutput', sprintf(msg,max(dist),'MaxIter')); + end + + pline = pline_o(idx,:); +endfunction + +function [dist ii] = maxdistance (p, idx) + + %% Separate the groups of points according to the edge they can divide. + func = @(x,y) x:y; + idxc = arrayfun (func, idx(1:end-1), idx(2:end), "UniformOutput",false); + points = cellfun (@(x)p(x,:), idxc, "UniformOutput",false); + + %% Build the edges + edges = [p(idx(1:end-1),:) p(idx(2:end),:)]; + edges = mat2cell (edges, ones(1,size(edges,1)), 4)'; + + %% Calculate distance between the points and the corresponding edge + [dist ii] = cellfun(@dd, points,edges,idxc); + +endfunction + +function [dist ii] = dd (p,e,idx) + [d pos] = distancePointEdge(p,e); + [dist ii] = max(d); + ii = idx(ii); +endfunction + +%!demo +%! t = linspace(0,1,100).'; +%! y = polyval([1 -1.5 0.5 0],t); +%! pline = [t y]; +%! +%! figure(1) +%! clf +%! plot (t,y,'-r;Original;','linewidth',2); +%! hold on +%! +%! tol = [8 2 1 0.5]*1e-2; +%! colors = jet(4); +%! +%! for i=1:4 +%! pline_ = simplifypolyline(pline,'tol',tol(i)); +%! msg = sprintf('-;%g;',tol(i)); +%! h = plot (pline_(:,1),pline_(:,2),msg); +%! set(h,'color',colors(i,:),'linewidth',2,'markersize',4); +%! end +%! hold off +%! +%! % --------------------------------------------------------- +%! % Four approximations of the initial polyline with decreasing tolerances. + +%!demo +%! P = [0 0; 3 1; 3 4; 1 3; 2 2; 1 1]; +%! func = @(x,y) linspace(x,y,5); +%! P2(:,1) = cell2mat( ... +%! arrayfun (func, P(1:end-1,1),P(2:end,1), ... +%! 'uniformoutput',false))'(:); +%! P2(:,2) = cell2mat( ... +%! arrayfun (func, P(1:end-1,2),P(2:end,2), ... +%! 'uniformoutput',false))'(:); +%! +%! P2s = simplifypolyline (P2); +%! +%! plot(P(:,1),P(:,2),'s',P2(:,1),P2(:,2),'o',P2s(:,1),P2s(:,2),'-ok'); +%! +%! % --------------------------------------------------------- +%! % Simplification of a polyline in the plane. diff --git a/inst/polygons2d/splitPolygons.m b/inst/polygons2d/splitPolygons.m new file mode 100644 index 0000000..e373668 --- /dev/null +++ b/inst/polygons2d/splitPolygons.m @@ -0,0 +1,64 @@ +## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## All rights reserved. +## +## Redistribution and use in source and binary forms, with or without +## modification, are permitted provided that the following conditions are met: +## +## 1 Redistributions of source code must retain the above copyright notice, +## this list of conditions and the following disclaimer. +## 2 Redistributions in binary form must reproduce the above copyright +## notice, this list of conditions and the following disclaimer in the +## documentation and/or other materials provided with the distribution. +## +## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +## +## The views and conclusions contained in the software and documentation are +## those of the authors and should not be interpreted as representing official +## policies, either expressed or implied, of the copyright holders. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{polygons} = } splitPolygons (@var{polygon}) +## Convert a NaN separated polygon list to a cell array of polygons. +## +## @var{polygon} is a N-by-2 array of points, with possibly couples of NaN values. +## The functions separates each component separated by NaN values, and +## returns a cell array of polygons. +## +## @seealso{polygons2d} +## @end deftypefn +function polygons = splitPolygons(polygon) + + if iscell(polygon) + # case of a cell array + polygons = polygon; + + elseif sum(isnan(polygon(:)))==0 + # single polygon -> no break + polygons = {polygon}; + + else + # find indices of NaN couples + inds = find(sum(isnan(polygon), 2)>0); + + # number of polygons + N = length(inds)+1; + polygons = cell(N, 1); + + # iterate over NaN-separated regions to create new polygon + inds = [0;inds;size(polygon, 1)+1]; + for i=1:N + polygons{i} = polygon((inds(i)+1):(inds(i+1)-1), :); + end + end + +endfunction diff --git a/inst/polygons2d/supportFunction.m b/inst/polygons2d/supportFunction.m new file mode 100644 index 0000000..7936aae --- /dev/null +++ b/inst/polygons2d/supportFunction.m @@ -0,0 +1,72 @@ +%% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr> +%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +%% All rights reserved. +%% +%% Redistribution and use in source and binary forms, with or without +%% modification, are permitted provided that the following conditions are met: +%% +%% 1 Redistributions of source code must retain the above copyright notice, +%% this list of conditions and the following disclaimer. +%% 2 Redistributions in binary form must reproduce the above copyright +%% notice, this list of conditions and the following disclaimer in the +%% documentation and/or other materials provided with the distribution. +%% +%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' +%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR +%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +%% +%% The views and conclusions contained in the software and documentation are +%% those of the authors and should not be interpreted as representing official +%% policies, either expressed or implied, of the copyright holders. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{h} = } suppportFunction (@var{polygon}) +%% @deftypefnx {Function File} {@var{h} = } suppportFunction (@var{polygon}, @var{n}) +%% @deftypefnx {Function File} {@var{h} = } suppportFunction (@var{polygon}, @var{v}) +%% Compute support function of a polygon +%% +%% H = supportFunction(POLYGON, N) +%% uses N points for suport function approximation +%% +%% H = supportFunction(POLYGON) +%% assume 24 points for approximation +%% +%% H = supportFunction(POLYGON, V) +%% where V is a vector, uses vector V of angles to compute support +%% function. +%% +%% @seealso{convexification} +%% @end deftypefn + +function h = supportFunction(polygon, varargin) + N = 24; + u = (0:2*pi/N:2*pi*(1-1/N)).'; + + if length(varargin)==1 + var = varargin{1}; + if length(var)==1 + N = var; + u = (0:2*pi/N:2*pi*(1-1/N)).'; + else + u = var(:); + end + end + + + h = zeros(size(u)); + + for i=1:length(u) + + v = repmat([cos(u(i)) sin(u(i))], [size(polygon, 1), 1]); + + h(i) = max(dot(polygon, v, 2)); + end + +endfunction diff --git a/inst/shape2d/curve2polyline.m b/inst/shape2d/curve2polyline.m new file mode 100644 index 0000000..389aec0 --- /dev/null +++ b/inst/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') diff --git a/inst/io/private/getSVGstrPath.m b/inst/shape2d/curveval.m index ae58654..b47ec77 100644 --- a/inst/io/private/getSVGstrPath.m +++ b/inst/shape2d/curveval.m @@ -1,5 +1,5 @@ -%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch> -%% +%% 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 @@ -13,8 +13,18 @@ %% You should have received a copy of the GNU General Public License %% along with this program. If not, see <http://www.gnu.org/licenses/>. -function strPath = getSVGstrPath(svg) +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{p} = } curveval (@var{curve}, @var{t}) +%% Evaluates parametric @var{curve} at @var{t}. +%% +%% @end deftypefn +function p = curveval (curve, t) + + dim = size (curve,1); + p = zeros (length(t), dim); - [s e te strPath] = regexpi(svg,'<[ ]*path(?:(?!/>).)*/>'); + for i = 1:dim + p(:,i) = polyval (curve(i,:), t); + end -end +endfunction diff --git a/inst/shape2d/shape2polygon.m b/inst/shape2d/shape2polygon.m index 1d7a987..86fb143 100644 --- a/inst/shape2d/shape2polygon.m +++ b/inst/shape2d/shape2polygon.m @@ -15,24 +15,22 @@ %% -*- texinfo -*- %% @deftypefn {Function File} {@var{polygon} = } shape2polygon (@var{shape}) -%% @deftypefnx {Function File} {@var{polygon} = } shape2polygon (@var{shape},@var{N}) +%% @deftypefnx {Function File} {@var{polygon} = } shape2polygon (@dots{},@var{property},@var{value},@dots{}) %% Transforms a 2D shape described by piecewise smooth polynomials into a polygon. %% %% @var{shape} is a n-by-1 cell where each element is a pair of polynomials %% compatible with polyval. %% @var{polygon} is a k-by-2 matrix, where each row represents a vertex. -%% @var{N} defines the number of points to be used in non-straigth edges. +%% The property-value pairs are passed to @code{curve2polyline}. %% -%% @seealso{polygon2shape, drawPolygon} +%% @seealso{polygon2shape, curve2poyline} %% @end deftypefn -function polygon = shape2polygon (shape, N=16) +function polygon = shape2polygon (shape, varargin) polygon = cell2mat ( ... - cellfun (@(x) func (x,N), shape,'UniformOutput',false) ); - - %% TODO simply the polygon based on curvature -% polygon = unique (polygon, 'rows'); + cellfun (@(x) curve2polyline(x,varargin{:}), shape,'UniformOutput',false) ); + polygon = simplifypolygon(polygon); if size (polygon, 1) == 1 polygon(2,1) = polyval (shape{1}(1,:), 1); @@ -41,14 +39,27 @@ function polygon = shape2polygon (shape, N=16) endfunction -function y = func (x,N) - - if size (x,2) > 2 - t = linspace (0,1-1/N,N).'; - y(:,1) = polyval (x(1,:), t); - y(:,2) = polyval (x(2,:), t); - else - y = x(:,2).'; - end - -endfunction +%!demo +%! shape = {[-93.172 606.368 -476.054 291.429; ... +%! -431.196 637.253 11.085 163.791]; ... +%! [-75.3626 -253.2337 457.1678 328.5714; ... +%! 438.7659 -653.6278 -7.9953 380.9336]; ... +%! [-89.5841 344.9716 -275.3876 457.1429; ... +%! -170.3613 237.8858 1.0469 158.0765];... +%! [32.900 -298.704 145.804 437.143; ... +%! -243.903 369.597 -34.265 226.648]; ... +%! [-99.081 409.127 -352.903 317.143; ... +%! 55.289 -114.223 -26.781 318.076]; ... +%! [-342.231 191.266 168.108 274.286; ... +%! 58.870 -38.083 -89.358 232.362]}; +%! +%! % Estimate a good tolerance +%! n = cell2mat(cellfun(@(x)curveval(x,rand(1,10)), shape, 'uniformoutput',false)); +%! dr = (max(n(:,1))-min(n(:,1)))*(max(n(:,2))-min(n(:,2)))*40; +%! p = shape2polygon (shape,'tol',dr); +%! +%! figure(1) +%! shapeplot(shape,'-b'); +%! hold on; +%! drawPolygon (p,'-or'); +%! hold off diff --git a/inst/shape2d/shapecentroid.m b/inst/shape2d/shapecentroid.m index a7ad8d5..8094b0b 100644 --- a/inst/shape2d/shapecentroid.m +++ b/inst/shape2d/shapecentroid.m @@ -66,9 +66,8 @@ endfunction %! 0 0]}; %! CoM = shapecentroid (boomerang) %! Gcentroid = centroid(shape2polygon(boomerang)) -%! %! figure(1); clf; -%! shapeplot(boomerang,10,'-o'); +%! shapeplot(boomerang,'-o'); %! hold on %! drawPoint(CoM,'xk;shape centroid;'); %! drawPoint(Gcentroid,'xr;point centroid;'); @@ -85,7 +84,7 @@ endfunction %! CoM = shapecentroid (Lshape) %! Gcentroid = centroid (shape2polygon (Lshape)) %! -%! shapeplot(Lshape,10,'-o'); +%! shapeplot(Lshape,'-o'); %! hold on %! drawPoint(CoM,'xk;shape centroid;'); %! drawPoint(Gcentroid,'xr;point centroid;'); diff --git a/inst/shape2d/shapeplot.m b/inst/shape2d/shapeplot.m index adeae7d..b820464 100644 --- a/inst/shape2d/shapeplot.m +++ b/inst/shape2d/shapeplot.m @@ -27,9 +27,11 @@ %% @seealso{drawPolygon, shape2polygon} %% @end deftypefn -function h = shapeplot(shape, N = 16, varargin) +function h = shapeplot(shape, varargin) - p = shape2polygon(shape, N); + n = cell2mat(cellfun(@(x)curveval(x,rand(1,5)), shape, 'uniformoutput',false)); + dr = (max(n(:,1))-min(n(:,1)))*(max(n(:,2))-min(n(:,2)))/100; + p = shape2polygon(shape,'tol', dr); h = drawPolygon(p,varargin{:}); endfunction diff --git a/inst/shape2d/shapetransform.m b/inst/shape2d/shapetransform.m index 6cf387e..ca04f9d 100644 --- a/inst/shape2d/shapetransform.m +++ b/inst/shape2d/shapetransform.m @@ -98,14 +98,12 @@ endfunction %! shape = shapetransform (shape,-T + [2; 0]); %! %! close -%! shapeplot (shape,10,'-r','linewidth',2) +%! shapeplot (shape,'-r','linewidth',2); %! hold on %! for i = 1:9 %! T = createRotation (i*pi/5)(1:2,1:2)/exp(0.3*i); -%! shapeplot (shapetransform(shape, T), 10, 'color',rand(1,3),'linewidth',2); +%! shapeplot (shapetransform(shape, T), 'color',rand(1,3),'linewidth',2); %! end %! hold off %! axis tight %! axis square - - |