diff --git a/process_results/parse_regenerating_fragments.m b/process_results/parse_regenerating_fragments.m
index 55226a2..4a2adce 100644
--- a/process_results/parse_regenerating_fragments.m
+++ b/process_results/parse_regenerating_fragments.m
@@ -1,809 +1,809 @@
 % PARSE_REGENERATING_FRAGMENT studies the properties of the vascular tissue during regeneration
 %
 % Blanchoud Group, UNIFR
 % Simon Blanchoud
 % 06/12/2021
 function parse_regenerating_fragments(frag_indx)
 
   % load the required packages
   pkg load matgeom
   pkg load image
   pkg load statistics
 
   % Absolute path to the directory containing the segmentation files
-  %fdir = '/home/guest/data/Nathalie/Quantifications/*.roi';
-  fdir = '/home/blanchou/Documents/Nathalie/Quantifications/*.roi';
+  fdir = '/home/guest/data/Nathalie/Quantifications/*.roi';
+  %fdir = '/home/blanchou/Documents/Nathalie/Quantifications/*.roi';
   files = glob(fdir);
   nfiles = length(files);
 
   % Determine which data to show
   if nargin < 1
     frag_indx = [1:nfiles];
   else
     frag_indx = frag_indx(isfinite(frag_indx) & (frag_indx <= nfiles));
   end
 
   % Structures to store the data
   all_tunic = cell(nfiles, 2, 1);
   all_ampullae = cell(nfiles, 2, 4);
   all_vessels = cell(nfiles, 2, 2);
   all_evals = cell(nfiles, 2, 4);
   all_quants = cell(nfiles, 2, 4);
 
   % Display?
   plot_segmentations = true;
 
   % Statistics
   nrepeats = 100;
   nothers = 5;
 
   disp('Parsing segmentations...')
   % Loop through all files
   for i=1:nfiles
 
     % Extract the parts of the filename
     [dir, name, ext] = fileparts (files{i});
     filename = fullfile(dir, name);
 
     % Inform the user on the progress
     disp([num2str(i) ': ' name]);
 
     % The actual files
     imfile = [filename '.tif'];
     ijfile = [filename '.zip'];
     umfile = [filename '.txt'];
     tgfile = [filename '.roi'];
 
     % Extract the segmentations
     ROIs = ReadImageJROI(ijfile);
     [props, data] = analyze_ROI(ROIs, 'Slice', 'Area', 'Length', 'Centroid');
 
     % Setup the various thresholds
     scale = csvread(umfile);
     tdist = (550/scale);
     rdist = ( 50/scale);
     pdist = ( 30/scale);
     vdist = ( 10/scale);
 
     % Measure the size of the dataset
     nframes = max(props(:,1));
     ntracks = length(data);
 
     % Prepare some temporary variables
     start = props(:,1)<nframes/2;
     first = [1 1];
     refs = [1 nframes];
 
     % Identify which ROI is what
     tunic = (~isnan(props(:,2)));
     ampullae = (~tunic & isnan(props(:,3)));
     vessels = (start & ~ampullae & ~tunic);
 
     % Loop through all the segmentatiosn
     for j=1:ntracks
 
       % Record the index of the first frame to be used for display
       if start(j) && first(1)
         refs(1) = props(j,1); 
         first(1) = false;
 
       % And the corresponding one for the regenerating fragment
       elseif ~start(j) && first(2)
         refs(2) = props(j,1); 
         first(2) = false;
       end
 
       % For the contour of the tunic, we simply copy the data
       if (tunic(j))
 
         % Check if we need to circularize the polygon
         pts = data{j};
         if (any(pts(1,:) ~= pts(end,:)))
           pts = pts([1:end 1], :);
         end
 
         % Store it
         all_tunic{i, 2-start(j)} = pts;
 
       % For the ampullae we build a neighborhood and find their periphery
       elseif (ampullae(j))
 
         % Connect the ampullae into a network
         pts = data{j};
         [links, nneigh] = span_vertices(pts, tdist);
         
         % Store them
         all_ampullae{i, 2-start(j), 1} = pts;
         all_ampullae{i, 2-start(j), 2} = links;
         all_ampullae{i, 2-start(j), 4} = nneigh;
 
       % For the vessels, we pool them and check for branchings
       elseif (vessels(j))
 
         % Pool the curves
         pts = data{j};
         all_vessels{i,2-start(j),1} = [all_vessels{i,2-start(j),1}; pts; NaN(1,2)];
 
         % Check for intersections with other curves
         for k=j+1:ntracks
           if ~vessels(k)
             continue;
           end
 
           % Intersect the two datasets
           newpts = data{k};
           cross = intersectPolylines(pts, newpts);
 
           % Also check if neighbouring points could be already close enough
           [dists, indk] = minDistancePoints(pts, newpts);
           if any(dists < vdist)
             cross = [cross; 0.5*(pts(dists<vdist,:) + newpts(indk(dists<vdist),:))];
           end
 
           % Keep everything
           all_vessels{i,2-start(j),2} = [all_vessels{i,2-start(j),2}; cross];
         end
 
       % Ignore the other types
       else
         % Nothing
       end
     end
 
     % Some post-processing
     for n=1:2
 
       % Set minimal sets
       % For the vessels
       if isempty(all_vessels{i,n,1})
         all_vessels{i,n,1} = NaN(1,2);
       end
       if isempty(all_vessels{i,n,2})
         all_vessels{i,n,2} = NaN(1,2);
       end
 
       % For the ampullae
       if isempty(all_ampullae{i,n,1})
         all_ampullae{i,n,1} = NaN(1,2);
       end
       if isempty(all_ampullae{i,n,2})
         all_ampullae{i,n,2} = [1 1];
       end
       if isempty(all_ampullae{i,n,3})
         all_ampullae{i,n,3} = NaN(1,2);
       end
       if isempty(all_ampullae{i,n,4})
         all_ampullae{i,n,4} = 0;
       end
 
       % For the tunic
       if isempty(all_tunic{i,n,1})
         all_tunic{i,n,1} = NaN(1,2);
       end
 
       % Determine the periphery of the ampullae network
       all_ampullae{i,n,3} = ampullae_periphery(all_ampullae{i,n,1}, all_ampullae{i,n,2}, all_tunic{i,n}, tdist);
 
       % Fuse the potentially duplicated crossings
       all_vessels{i,n,2} = merge_points(all_vessels{i,n,2}, pdist);
 
       % Define additional non-niches evaluation points
       [all_evals{i,n,2}, all_evals{i,n,3}, all_evals{i,n,4}] = mesh_fragment(all_tunic{i,n}, all_ampullae{i,n,3}, tdist/4);
     end
 
     % Extract the segmentation for the niches
     try
       niches = ReadImageJROI(tgfile);
 
       % Get the position of the starting points
       pt = niches.mfCoordinates(niches.vnSlices==refs(1), :);
       all_evals{i, 1, 1} = pt;
 
       % And of the regeneration niches
       pt = niches.mfCoordinates(niches.vnSlices==refs(2), :);
       all_evals{i, 2, 1} = pt;
 
     % Nothing in that file it seems
     catch
       all_evals{i, 1, 1} = NaN(1,2);
       all_evals{i, 2, 1} = NaN(1,2);
     end
 
     % All the plotting if it is required
     if plot_segmentations && any(frag_indx==i)
       figure;
       hax = NaN(1,2);
 
       % Do each part separately
       for n=1:2
         % Create the axis
         hax(n) = subplot(1,2,n); hold on;
         set(hax(n), 'position', [0.01+(n-1)*0.5 0.01 0.48 0.98]);
 
         % Plot the image
         img = imread(imfile, refs(n));
         image(img);
         axis('equal', 'ij');
 
         % Display the evaluation points
         scatter(all_evals{i,n,2}(:,1), all_evals{i,n,2}(:,2), 'k')
         scatter(all_evals{i,n,3}(:,1), all_evals{i,n,3}(:,2), 'r')
         scatter(all_evals{i,n,4}(:,1), all_evals{i,n,4}(:,2), 'w')
 
         % Display the ampullae
         pts = all_ampullae{i,n,1};
         links = all_ampullae{i,n,2};
         periph = all_ampullae{i,n,3};
         scatter(pts(:,1), pts(:,2), 'k')
         plot([pts(links(:,1),1) pts(links(:,2),1)].', [pts(links(:,1),2) pts(links(:,2),2)].', 'k')
         plot(periph(:,1), periph(:,2), 'm')
 
         % The tunic
         plot(all_tunic{i,n}(:,1), all_tunic{i,n}(:,2), 'r')
 
         % The vessels
         plot(all_vessels{i,n,1}(:,1), all_vessels{i,n,1}(:,2), 'b')
         scatter(all_vessels{i,n,2}(:,1), all_vessels{i,n,2}(:,2), 'b', 'filled')
 
         % The niches
         scatter(all_evals{i,n,1}(:,1), all_evals{i,n,1}(:,2), 'r', 'filled')
       end
 
       % Link, scale and fancy the plots
       linkaxes(hax)
       s = min(all_tunic{i,1}, [], 1);
       l = max(all_tunic{i,1}, [], 1);
       trange = [s(1) l(1) s(2) l(2)];
       axis(hax(1), trange, 'off')
       axis(hax(2), 'off')
 
       keyboard
     end
 
     % Quantify each part of the fragment separately
     for n=1:2
 
       % First the tunic
       pts = all_tunic{i,n};
 
       % Compute general statistics
       [gtunic, ctunic, tunic] = quantify_area(pts, scale, rdist);
 
       % Then the ampullae
       amps = all_ampullae{i,n,1};
       pamps = all_ampullae{i,n,3};
       namps = all_ampullae{i,n,4};
 
       % Other general statistics
       [gamps, camps] = quantify_area(pamps, scale, rdist);
       gamps = [gamps, gamps(1)/gtunic(1), gamps(1)/size(amps,1)];
 
       % And the vessels
       pts = all_vessels{i,n,1};
       crosses = all_vessels{i,n,2};
 
       % General statistics
       [gvessels, vessels] = quantify_path(pts, scale, rdist);
 
       % Compute the detailed statistics for each set of evaluation points 
       for ev=1:4
 
         % The current evaluation points
         ref = all_evals{i,n,ev};
 
         % The distance to the tunic centroid
         dcentroid = sqrt(sum((bsxfun(@minus, ref, ctunic)).^2,2))*scale;
 
         % The radial four first moments
         vtunic = central_moments(ref, tunic, scale);
 
         % Combined with the general statistics
         tvals = [repmat(gtunic, size(ref,1), 1), dcentroid, vtunic];
 
         % Distance to the centroid
         dcamps = sqrt(sum((bsxfun(@minus, ref, camps)).^2,2))*scale;
 
         % The moments of the periphery and of the neigbours
         vamps = [central_moments(ref, pamps, scale), central_moments(ref, amps, scale, namps, tdist)];
 
         % All together
         avals = [repmat(gamps, size(ref,1), 1), dcamps, vamps];
 
         % The moments of the vessels and of the crosses
         vvessl = [central_moments(ref, vessels, scale), central_moments(ref, crosses, scale)];
 
         % All together
         vvals = [repmat(gvessels, size(ref,1), 1), vvessl];
 
         % Store everything
         all_quants{i,n,ev} = [tvals, avals, vvals];
       end
     end
   end
 
   all_features = {'tunic:area', ...
           'tunic:periph', ...
           'tunic:dist2center', ...
           'tunic:min_a', 'tunic:mean_a', 'tunic:stdev_a', 'tunic:skew_a', 'tunic:kurtosis_a', ...
           'tunic:min_d', 'tunic:mean_d', 'tunic:stdev_d', 'tunic:skew_d', 'tunic:kurtosis_d', ...
           'amps:area', ...
           'amps:periph', ...
           'amps:ratio2tunic', ...
           'amps:density', ...
           'amps:dist2center', ...
           'amps_periph:min_a', 'amps_periph:mean_a', 'amps_periph:stdev_a', 'amps_periph:skew_a', 'amps_periph:kurtosis_a', ...
           'amps_periph:min_d', 'amps_periph:mean_d', 'amps_periph:stdev_d', 'amps_periph:skew_d', 'amps_periph:kurtosis_d', ...
           'amps_pos:min_a', 'amps_pos:mean_a', 'amps_pos:stdev_a', 'amps_pos:skew_a', 'amps_pos:kurtosis_a', ...
           'amps_pos:min_d', 'amps_pos:mean_d', 'amps_pos:stdev_d', 'amps_pos:skew_d', 'amps_pos:kurtosis_d', ...
           'amps_pos:nneighb', 'amps_pos:mean_nneigh', 'amps_pos:stdev_nneigh', 'amps_pos:skew_nneigh', 'amps_pos:kurtosis_nneigh', ...
           'vessels:length', ...
           'vessels:min_a', 'vessels:mean_a', 'vessels:stdev_a', 'vessels:skew_a', 'vessels:kurtosis_a', ...
           'vessels:min_d', 'vessels:mean_d', 'vessels:stdev_d', 'vessels:skew_d', 'vessels:kurtosis_d', ...
           'crosses:min_a', 'crosses:mean_a', 'crosses:stdev_a', 'crosses:skew_a', 'crosses:kurtosis_a', ...
           'crosses:min_d', 'crosses:mean_d', 'crosses:stdev_d', 'crosses:skew_d', 'crosses:kurtosis_d'};
 
   % All the plotting of the weights if it is requested
   if plot_segmentations
 
     % For all the file
     for i=1:nfiles
 
       % Display the requested ones
       if any(frag_indx==i)
         figure;
 
         % One for each picture
         for n=1:2
 
           % Create the axes
           subplot(1,2,n);
 
           % Combine all submatrices
           vals = cat(1,all_quants{i,n,:});
 
           % Display
           imagesc(abs(log(vals)));
         end
 
         keyboard
       end
     end
   end
 
   % Preparing the variables required for the statistics
   nfeatures = size(all_quants{1,1,1}, 2);
   Hs = NaN(3,nfeatures,nrepeats);
   Ps = Hs;
 
   disp('Bootstrapping statistics...')
 
   % Loop through all files to subset them
   counts = zeros(1,3);
   curr_quants = all_quants;
   for i=1:nfiles
 
     % Extract the full data
     curr_niches = all_quants{i,1,1};
     curr_ampull = all_quants{i,1,2};
     curr_tunics = all_quants{i,1,3};
     curr_outsid = all_quants{i,1,4};
 
     % Discard the non-finite data
     is_niches = all(isfinite(curr_niches),2);
     is_ampull = all(isfinite(curr_ampull),2);
     is_tunics = all(isfinite(curr_tunics),2);
     is_outsid = all(isfinite(curr_outsid),2);
 
     curr_niches = curr_niches(is_niches,:);
     curr_ampull = curr_ampull(is_ampull,:);
     curr_tunics = curr_tunics(is_tunics,:);
     curr_outsid = curr_outsid(is_outsid,:);
 
     % Identify the closest points
     maps = cat(1, all_evals{i,1,2:4});
     maps = maps([is_ampull; is_tunics; is_outsid],:);
     feat = [curr_ampull; curr_tunics; curr_outsid];
 
     dist = bsxfun(@minus, all_evals{i,1,1}(:,1), maps(:,1).').^2 + bsxfun(@minus, all_evals{i,1,1}(:,2), maps(:,2).').^2;
     [junk, dindx] = sort(min(dist, [], 1));
 
     % Extract a random subset of a size proportional to the number of niches
     nniches = size(curr_niches,1)*nothers;
 
     curr_closes = feat(dindx(1:nniches),:);
     curr_maps = maps(dindx(1:nniches),:);
 
     % Draw the sets for multiple subsets
     namp = size(curr_ampull,1);
     ntun = size(curr_tunics,1);
     nout = size(curr_outsid,1);
 
     mapa = NaN(nrepeats, nniches);
     mapt = NaN(nrepeats, nniches);
     mapo = NaN(nrepeats, nniches);
     for r=1:nrepeats
       mapa(r,:) = randperm(namp,nniches) + counts(1);
       mapt(r,:) = randperm(ntun,nniches) + counts(2);
       mapo(r,:) = randperm(nout,nniches) + counts(3);
     end
 
     % Store everything
     curr_quants{i,1,1} = curr_niches;
     curr_quants{i,1,2} = curr_ampull;
     curr_quants{i,1,3} = curr_tunics;
     curr_quants{i,1,4} = curr_closes;
     curr_quants{i,1,5} = curr_outsid;
 
     curr_quants{i,2,2} = mapa.';
     curr_quants{i,2,3} = mapt.';
     curr_quants{i,2,5} = mapo.';
 
     counts = counts + [namp, ntun, nout];
   end
 
   % Gather all the results
   all_niches = cat(1, curr_quants{:,1,1});
   all_ampull = cat(1, curr_quants{:,1,2});
   all_tunics = cat(1, curr_quants{:,1,3});
   all_closes = cat(1, curr_quants{:,1,4});
   all_outsid = cat(1, curr_quants{:,1,5});
 
   all_mapamp = cat(1, curr_quants{:,2,2});
   all_maptun = cat(1, curr_quants{:,2,3});
   all_mapout = cat(1, curr_quants{:,2,5});
 
   % Compute t-tests one-by-one for each feature
   for r=1:nrepeats
     for f=1:nfeatures
       [Hs(1,f,r),Ps(1,f,r)] = ttest2(all_niches(:,f), all_closes(:,f));
       [Hs(2,f,r),Ps(2,f,r)] = ttest2(all_niches(:,f), all_ampull(all_mapamp(:,r),f));
       [Hs(3,f,r),Ps(3,f,r)] = ttest2(all_niches(:,f), all_tunics(all_maptun(:,r),f));
       [Hs(4,f,r),Ps(4,f,r)] = ttest2(all_niches(:,f), all_outsid(all_mapout(:,r),f));
     end
   end
 
   % Average the bootstrapping
   hits = mean(Hs,3);
   pval = mean(Ps,3);
 
   % Keep and sort the best ones
   [p,best] = sort(pval(2,:));
 
   disp('Displaying the 9 best features')
   figure;
   for f=1:9
     curr = best(f);
     [mn,sn] = mymean(all_niches(:,curr));
     [mc,sc] = mymean(all_closes(:,curr));
     [ma,sa] = mymean(all_ampull(:,curr));
     [mt,st] = mymean(all_tunics(:,curr));
     [mo,so] = mymean(all_outsid(:,curr));
 
     ha=subplot(3,3,f);hold on;
     bar([mn, mc, ma, mt, mo].');
     errorbar([mn, mc, ma, mt, mo], [sn, sc, sa, st, mo], '~.');
     title({strrep(all_features{curr}, '_', '.');['p-val: ' num2str(pval(:,curr).')]})
     set(gca, 'XTickLabel', {'', 'Niches', 'Closest', 'Ampullae', 'Tunic', 'Outside', ''})
   end
 
   keyboard
 
   return
 end
 
 % Create the ampullae network
 function [outer, tissue, inner] = mesh_fragment(tunic, ampullae, dt)
 
   mins = min(tunic);
   maxs = max(tunic);
 
   [X,Y] = meshgrid([mins(1):dt:maxs(1)], [mins(2):dt:maxs(2)]);
   pts = [X(:) Y(:)];
 
   inside = isPointInPolygon(pts, ampullae);
   outside = isPointInPolygon(pts(~inside, :), tunic);
 
   outer = pts(~inside, :);
   inner = pts(inside,:);
   tissue = outer(outside,:);
   outer = outer(~outside,:);
 
   return
 end
 
 % Interpolate the polygons
 function newpts = linear_interpolation(pts, step)
 
   newpts = pts;
   if numel(pts)<4
     return;
   end
 
   dpoly = sqrt(sum((pts(1:end-1, :) - pts(2:end, :)).^2, 2));
   dpoly = [0; dpoly];
   dpoly(isnan(dpoly)) = 0;
 
   x = cumsum(dpoly);
   newx = [0:step:max(x)].';
 
   warning('off');
   newpts = interp1(x, pts, newx);
   warning('on');
 
   return
 end
 
 % Quantify a shape
 function [vals, centroid, contour] = quantify_area(pts, scale, rdist)
 
   [peri, contour] = quantify_path(pts, scale, rdist);
 
   p = pts(1:end-1,1).*pts(2:end,2) - pts(2:end,1).*pts(1:end-1,2);
   valids = isfinite(p);
   area = sum(p(valids)) / 2;
 
   if area==0
     centroid = NaN(1,2);
   else
     proj = ((pts(1:end-1,:)+pts(2:end,:)) .* [p p]) / 6 / area;
     centroid = sum(proj(valids,:));
   end
 
   vals = [abs(area*scale*scale), peri];
 
   return
 end
 
 % Quantify a path
 function [len, resampl] = quantify_path(pts, scale, rdist)
 
   len = 0;
   start = 1;
   nans = find(isnan(pts(:,1)));
   for i=1:length(nans)
     len = len + polygonLength(pts(start:nans(i)-1,:)*scale);
     start = nans(i)+1;
   end
   len = len + polygonLength(pts(start:end,:)*scale);
 
   resampl = linear_interpolation(pts, rdist);
 
   return
 end
 
 % Compute the 4 first central movements
 function vals = central_moments(ref, pts, scale, weights, norm)
 
   nrefs = size(ref, 1);
   valids = all(isfinite(pts),2);
   pts = pts(valids,:);
 
   if nargin < 4
     vals = NaN(nrefs, 10);
     if isempty(pts)
       return;
     end
 
     for i=1:nrefs
       centered = [pts(:,1) - ref(i, 1) pts(:,2) - ref(i, 2)]*scale;
       a = atan2(centered(:,2), centered(:,1));
       d = sqrt(sum(centered.^2, 2));
 
       ma = mean(a);
       sa = std(a,1);
       wa = mean((a - ma).^3) / sa.^3;
       ka = mean((a - ma).^4) / sa.^4;
 
       md = mean(d);
       sd = std(d,1);
       wd = mean((d - md).^3) / sd.^3;
       kd = mean((d - md).^4) / sd.^4;
 
       v = [min(a), ma, sa, wa, ka, min(d), md, sd, wd, kd];
       vals(i,:) = v;
     end
   else
     vals = NaN(nrefs, 15);
     if isempty(pts)
       return;
     end
 
     for i=1:nrefs
       centered = [pts(:,1) - ref(i, 1) pts(:,2) - ref(i, 2)]*scale;
       a = atan2(centered(:,2), centered(:,1));
       d = sqrt(sum(centered.^2, 2));
 
       ma = mean(a);
       sa = std(a,1);
       wa = mean((a - ma).^3) / sa.^3;
       ka = mean((a - ma).^4) / sa.^4;
 
       [id, indx] = min(d);
       md = mean(d);
       sd = std(d,1);
       wd = mean((d - md).^3) / sd.^3;
       kd = mean((d - md).^4) / sd.^4;
 
       p = exp(-d.^2 / (4*norm^2));
       dp = sum(p);
       if dp<1e-10
         norm = norm*2;
         p = exp(-d.^2 / (4*norm^2));
         dp = sum(p);
 
         if dp~=1
           p = ones(size(p));
         end
       end
       p = p / sum(p);
 
       iw = weights(indx);
       mw = sum(p .* weights);
       sw = sum(p .* (weights - mw).^2);
       ww = sum(p .* (weights - mw).^3) / sw.^3;
       kw = sum(p .* (weights - mw).^4) / sw.^4;
 
       v = [min(a), ma, sa, wa, ka, id, md, sd, wd, kd, iw, mw, sw, ww, kw];
       vals(i,:) = v;
     end
   end
 
   return
 end
 
 % Identify the periphery of the ampullae network
 function periph = ampullae_periphery(pts, edges, tunic, thresh)
 
   npts = size(pts, 1);
   sets = zeros(npts, 1);
   nsets = 0;
 
   resampl = linear_interpolation(tunic, thresh/10);
   ntunic = size(resampl, 1);
   sqthresh = 5*thresh*thresh;
 
   for i=1:npts
     first = find(sets==0, 1);
 
     if (isempty(first))
       break;
     end
 
     nsets = nsets + 1;
     sets(first) = nsets;
     curr_set = first;
     for j=1:npts
       conx = ismember(edges, curr_set);
       conx = edges(any(conx, 2), :);
       new_set = unique(conx(:));
       sets(new_set) = nsets;
 
       if numel(curr_set) == numel(new_set)
         break
       else
         curr_set = new_set;
       end
     end
   end
 
   indxs = [];
   curr_indx = NaN;
   for i=1:ntunic
     d = ((pts(:,1) - resampl(i,1)).^2 + (pts(:,2) - resampl(i,2)).^2);
     [junk, indx] = min(d);
     if indx~= curr_indx && ~(any(indxs==indx))
 
       if (numel(indxs)>0)
         l = sum((pts(indxs(end),:) - pts(indx,:)).^2);
 
         if l > sqthresh && sets(indxs(end))==sets(indx)
           p = grShortestPath(pts, edges, indxs(end), indx);
           p = p(:).';
           indxs = [indxs(1:end-1), p(1:end-1)];
         end
       end
 
       indxs(end+1) = indx;
       curr_indx = indx;
     end
   end
   indxs(end+1) = indxs(1);
   periph = pts(indxs([1:end 1]),:);
 
   %figure;plot(resampl(:,1), resampl(:,2), 'r');
   %hold on;
   %scatter(pts(:,1), pts(:,2), 'k');
   %plot([pts(edges(:,1),1) pts(edges(:,2),1)].', [pts(edges(:,1),2) pts(edges(:,2),2)].', 'k')
   %plot(periph(:,1), periph(:,2), 'm');
 
 	return;
 end
 
 % Mesh the vertices with a connectivity network
 function [edges, nneigh] = span_vertices(pts, thresh)
 
   tri = delaunay(pts);
   thresh2 = thresh*thresh;
 
   ntri = size(tri,1);
   props = NaN(ntri, 3);
   for i=1:ntri
     vertices = tri(i, [1:3 1]);
 
     pt1 = pts(vertices(1),:);
     pt2 = pts(vertices(2),:);
     pt3 = pts(vertices(3),:);
 
     area = triangleArea(pt1, pt2, pt3);
     side = sqrt(sum([pt1-pt2; pt2 - pt3; pt3 - pt1].^2, 2));
     [l,imax] = max(side);
     peri = sum(side);
     props(i, 1) = peri/area;
     props(i, 2:3) = vertices(imax:imax+1);
   end
   tiny = (abs(props(:,1)) > 1);
   tri = tri(~tiny,:);
 
   edges = [tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])];
   [edges] = unique(edges, 'rows');
 
   bads = ismember(edges, props(tiny,2:3), 'rows');
   edges = edges(~bads,:);
 
   dist = bsxfun(@minus, pts(:,1), pts(:,1).').^2 + bsxfun(@minus, pts(:,2), pts(:,2).').^2;
   ind = sub2ind(size(dist), edges(:,1), edges(:,2));
 
   edist = dist(ind);
   edges = edges(edist <= thresh2,:);
 
   nneigh = histc(edges(:),[1:size(pts,1)]);
 
   return
 end
 
 % Fuse close points
 function merged = merge_points(pts, thresh)
 
   merged = pts;
   if (size(pts, 1) < 2)
     return;
   end
 
   thresh = thresh*thresh;
   dist = bsxfun(@minus, pts(:,1), pts(:,1).').^2 + bsxfun(@minus, pts(:,2), pts(:,2).').^2;
 
   groups = {};
   npts = size(pts,1);
   candidates = [1:npts];
 
   for i=1:npts
     if isempty(candidates)
       break;
     end
 
     ncands = size(candidates, 1);
     cluster = false(ncands,1);
     cluster(1) = true;
     num = 1;
     for j=1:ncands
       d = dist(candidates(cluster), candidates);
       cluster = any(d<thresh, 1);
       new = sum(cluster);
 
       if new == num || new==ncands
         break
       end
       num = new;
     end
     groups{i} = candidates(cluster);
     candidates = candidates(~cluster);
   end
 
   ngroups = length(groups);
   merged = NaN(ngroups, 2);
   for i=1:ngroups
     merged(i,:) = mean(pts(groups{i},:),1);
   end
 
   return
 end