Code covered by the BSD License

### Highlights fromToolbox Fast Marching

from Toolbox Fast Marching by Gabriel Peyre
A toolbox for the computation of the Fast Marching algorithm in 2D and 3D.

compute_shape_boundary(M)
```function bound = compute_shape_boundary(M)

% compute_shape_boundary - extract boundary points
%
% bound = compute_shape_boundary(M);
%
%   bound is oriented counter clockwise.
%   This is the 8 connectivity boundary.
%
%   Copyright (c) 2007 Gabriel Peyre

n = size(M,1);

% enforce only 1 connected component
epsilon = 1e-10;
options.nb_iter_max = Inf;
options.Tmax = 2*n;
W = epsilon+M;
I = find(M(:)==1); x = linspace(-1,1,n);
[Y,X] = meshgrid(x,x);
[tmp,c] = min( X(I).^2 + Y(I).^2 );
[a b] = ind2sub([n n],I(c));
options.constraint_map = [];
[d,tmp] = perform_fast_marching(W, [a;b], options);
I = find(d>1e5 | d<0); M(I) = 0;

% points on the boundary
if 0
% outside the shape
h = [0 1 0; 1 1 1; 0 1 0];
Mh = perform_convolution(M,h);
A = Mh>0 & M==0;
else
% inside the shape
h = ones(3);
Mh = perform_convolution(M,h);
A = Mh>0 & Mh<9 & M==1;
end

% 8 connectivity
wx = [1 -1 0  0 -1 1 -1  1];
wy = [0  0 1 -1 -1 1  1 -1];

% order the points
a = find(A(:)==1); a = a(1);
[x,y] = ind2sub(size(A),a);

boundx = x;
boundy = y;
A(x,y) = 0;

while sum(A(:))>0
% find closest point
for k=1:8
if x+wx(k)>0 && x+wx(k)<=n && y+wy(k)>0 && y+wy(k)<=n && A(x+wx(k),y+wy(k))>0
x = x+wx(k);
y = y+wy(k);
break;
end
end
if A(x,y)==0
% try to find the existing point the closest to the last one
I = find(A>0);
[xlist,ylist] = ind2sub(size(A),I);
d = (xlist-x).^2 + (ylist-y).^2;
[tmp,I] = min(d(:));
x = xlist(I(1)); y = ylist(I(1));
end
A(x,y) = A(x,y)-1;
boundx(end+1) = x;
boundy(end+1) = y;
end
bound = cat(2,boundx(:), boundy(:));

% reorient the curve so that it is clockwise oriented
p = size(bound,1);
c = mean(bound,1);
v = bound - repmat(c, [p 1]);
% cross product
z = v(1:end-1,1).*v(2:end,2) - v(1:end-1,2).*v(2:end,1);
if sum(z<0)>sum(z<0)
% should flip the curve
bound = bound(end:-1:1,:);
end```