polylinesimplify.m

function [newnodes, len] = polylinesimplify(nodes, minangle)
%
% [newnodes, len]=polylinesimplify(nodes, minangle)
%
% Calculate a simplified polyline by removing nodes where two adjacent
% segment have an angle less than a specified limit
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    node: an N x 3 array defining each vertex of the polyline in
%          sequential order
%    minangle:(optional) minimum segment angle in radian, if not given, use
%          0.75*pi
%
% output:
%    newnodes: the updated node list; start/end will not be removed
%    len: the length of each segment between the start and the end points
%
%
% -- this function is part of brain2mesh toolbox (http://mcx.space/brain2mesh)
%    License: GPL v3 or later, see LICENSE.txt for details
%

if (nargin < 2)
    minangle = 0.75 * pi;
end

v = segvec(nodes(1:end - 1, :), nodes(2:end, :));
ang = acos(max(min(sum(-v(1:end - 1, :) .* (v(2:end, :)), 2), 1), -1));

newnodes = nodes;
newv = v;
newang = ang;

idx = find(newang < minangle);

while (~isempty(idx))
    newnodes(idx + 1, :) = [];
    newv(idx + 1, :) = [];
    newang(idx) = [];
    idx = unique(idx - (0:(length(idx) - 1))');
    idx1 = idx(idx < size(newnodes, 1));
    newv(idx1, :)  = segvec(newnodes(idx1, :), newnodes(idx1 + 1, :));
    idx1 = idx(idx < size(newv, 1));
    newang(idx1)  = acos(sum(-newv(idx1, :) .* (newv(idx1 + 1, :)), 2));
    idx0 = idx(idx > 1);
    newang(idx0 - 1) = acos(sum(-newv(idx0 - 1, :) .* (newv(idx0, :)), 2));
    idx = find(newang < minangle);
end

if (nargout > 1)
    len = newnodes(1:end - 1, :) - newnodes(2:end, :);
    len = sqrt(sum(len .* len, 2));
end

function v = segvec(n1, n2)

v = n2 - n1;
normals = sqrt(sum(v .* v, 2));
v = v ./ repmat(normals, 1, size(v, 2));