diff --git a/install_SiBoRG.m b/install_SiBoRG.m
index 1ab8c92..97c6071 100755
--- a/install_SiBoRG.m
+++ b/install_SiBoRG.m
@@ -1,180 +1,183 @@
 function install_SiBoRG(recompile)
   % INSTALL_SIBORG adds the package to the work path, creates the required directories,
   % handles the dependencies and compiles the necessary MEX libraries.
   %
   % INSTALL_SIBORG('UNINSTALL') removes this package from the work path.
   %
   % INSTALL_SIBORG(RECOMPILE) if RECOMPILE is TRUE, recompiles all MEX files.
   %
   % Blanchoud group, University of Fribourg
   % Simon Blanchoud
   % 04.11.2019
 
   % By default, do not recompile everything
   if (nargin == 0)
     recompile = false;
 
     % In that particular case, delete the directories
   elseif (strcmp(recompile, 'uninstall'))
     cell_folder = which('install_SiBoRG');
     if (~isempty(cell_folder))
       [current_dir, junk, junk] = fileparts(cell_folder);
       folders = strsplit(path(), ':');
 
       % Remove the used directories from the work path
       for i=1:length(folders)
         if (strncmp(folders{i}, current_dir, length(current_dir)))
           rmpath(folders{i});
         end
       end
       savepath;
 
       % Remove the added java libraries if any
       fname = tilde_expand('~/.octaverc');
       fnew = tilde_expand('~/octave.rc');
 
       % We open the file and a new new
       fid = fopen(fname, 'rt');
       if fid >= 0
         fnw = fopen(fnew, 'wt');
         if fnw >= 0
 
           javastring = ['javaaddpath(''' current_dir];
 
           % And we copy everything that isn't in our directory
           line = fgetl(fid);
           while ischar(line)
             if (~strncmp(line, javastring, length(javastring)))
               fwrite(fnw, line);
             end
             line = fgetl(fid);
           end
 
           fclose(fnw);
           movefile(fnew, fname);
         end
 
         fclose(fid);
       end
 
       % Move to the root of SiBoRG
       [root_dir, junk, junk] = fileparts(current_dir);
       cd(root_dir);
 
       % Actually delete this package if required
       btn = questdlg ("Do you want to delete SiBoRG completely (cannot be undone)?", "DELETE FILES??", "Yes", "No", "No");
       if (strcmp (btn, "Yes"))
         disp('(not really) DELETED')
         % rmdir(current_dir, 's');
       end
     end
 
     return;
   end
 
   % Start by moving inside the containing folder
   cell_folder = which('install_SiBoRG');
   [current_dir, junk, junk] = fileparts(cell_folder);
   [root_dir, junk, junk] = fileparts(current_dir);
   cd(current_dir);
 
   % Add the proper directories to the work path
   subdirs = ls();
   folders = strsplit(path(), ':');
   for i=1:size(subdirs, 1)
     dir_name = strtrim(subdirs(i,:));
     if (isdir(dir_name))
       if (dir_name(1) ~= '.' && dir_name(1) ~= '_')
         full_path = fullfile(current_dir, dir_name);
         if (~any(strncmp(folders, full_path, length(full_path))))
           addpath(full_path);
         end
       end
     end
   end
   addpath(current_dir);
   savepath;
 
   % Added the required java libraries if any to octaverc
   fname = tilde_expand('~/.octaverc');
 
   % We open the file to add our files
   fid = fopen(fname, 'at');
   if fid >= 0
 
     % And we copy everything that is a jar in our libraries
     libraries = dir(fullfile('libraries', '*.jar'));
     fdisp(fid, ' ');
 
     % We both run the command and store it for the next boot
     for i=1:length(libraries)
       library = fullfile(pwd, 'libraries', libraries(i).name);
       fdisp(fid, ['javaaddpath(''' library ''');']);
       javaaddpath(library);
     end
 
     fclose(fid);
   end
 
   % Install the required packages
   if isempty(pkg('list', 'image'))
     pkg install -forge image
   end
   if isempty(pkg('list', 'io'))
     pkg install -forge io
   end
   if isempty(pkg('list', 'statistics'))
     pkg install -forge statistics
   end
   if isempty(pkg('list', 'control'))
     pkg install -forge control
   end
   if isempty(pkg('list', 'signal'))
     pkg install -forge signal
   end
   if isempty(pkg('list', 'tisean'))
     pkg install -forge tisean
   end
+  if isempty(pkg('list', 'matgeom'))
+    pkg install -forge matgeom
+  end
 
   % Try to compile the necessary MEX files
   cd('MEX');
   files = ls('*_mex.c*');
   troubles = false;
   for i=1:size(files, 1)
     fname = strtrim(files(i,:));
     [junk, no_ext, ext] = fileparts(fname);
     if (recompile || exist(no_ext) ~= 3)
       failure = mex(fname);
 
       if (failure == 1)
         troubles = true;
         wtharning('SiBoRG:MEX', {'Could not compile the required MEX function!' ME.message});
       end
     end
   end
 
   cd(root_dir);
 
   % These folders are required as well
   if (~exist('TmpData', 'dir'))
     mkdir('TmpData');
   end
   if (~exist('export', 'dir'))
     mkdir('export');
   end
 
   % Confirm to the user that everything went fine
   if (troubles)
     disp('Installation (almost) successful...');
   else
     disp('Installation successful !');
   end
 
   % Gnu GPL notice
   fprintf(1, ['\nSimon''s Botrylloides Regeneration Group plateform,  Copyright (C) 2019  Simon Blanchoud\n', ...
     'This program comes with ABSOLUTELY NO WARRANTY;\n', ...
     'This is free software, and you are welcome to redistribute it\n', ...
     'under certain conditions; read licence.txt for details.\n\n']);
 
   return;
 end
diff --git a/process_results/parse_regenerating_fragments.m b/process_results/parse_regenerating_fragments.m
index 8cb94f4..ce6c809 100644
--- a/process_results/parse_regenerating_fragments.m
+++ b/process_results/parse_regenerating_fragments.m
@@ -1,638 +1,648 @@
 % PARSE_REGENERATING_FRAGMENT studies the properties of the vascular tissue during regeneration
 %
 % Blanchoud Group, UNIFR
 % Simon Blanchoud
 % 06/12/2021
 function parse_regenerating_fragments
 
   % load the required packages
   pkg load matgeom
   pkg load image
 
   % Absolute path to the directory containing the segmentation files
   fdir = '/home/guest/data/Nathalie/Quantifications/*.roi';
   files = glob(fdir);
   nfiles = length(files);
 
   % 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;
 
   % 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(['Processing: ' 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 = (450/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, periph, 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), 3} = periph;
         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
 
       % 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/2);
     end
 
     % Extract the segmentation for the niches
-    niches = ReadImageJROI(tgfile);
+    try
+      niches = ReadImageJROI(tgfile);
+
+      % Get the position of the starting points
+      pt = niches.mfCoordinates(niches.vnSlices==refs(1), :);
+      all_evals{i, 1, 1} = pt;
 
-    % 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;
 
-    % 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
       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')
     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 the plotting of the weights if it is requested
   if plot_segmentations
     figure;
 
     % For all the file
     for i=1:nfiles
 
       % One for each picture
       for n=1:2
 
         % Create the axes
         subplot(nfiles,2,2*(i-1)+n);
 
         % Combine all submatrices
         vals = cat(1,all_quants{i,n,:});
 
         % Display
         imagesc(abs(log(vals)));
       end
     end
   end
 
   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
 
 % Mesh the vertices with a connectivity network
 function [edges, periph, 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,:);
 
   periph = outer_edges(pts, edges, thresh);
   nneigh = histc(edges(:),[1:size(pts,1)]);
 
   return
 end
 
 % Identify the periphery of a network
 function periph = outer_edges(pts, edges)
 
   npts = size(pts, 1);
   sets = zeros(npts, 1);
   nsets = 0;
   dist = sqrt((pts(edges(:,1),1) - pts(edges(:,2),1)).^2 + (pts(edges(:,1),2) - pts(edges(:,2),2)).^2);
   periph = [];
 
   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
 
   pindx = [1:npts].';
   for i=1:nsets
     curr = (sets==i);
     tmp = pts(curr,:);
 
     if (size(tmp,1)<=3)
       full = pindx(curr);
     else
       pcurr = pindx(curr);
       convs = pcurr(convhull(tmp(:,1), tmp(:,2)));
 
       full = [];
       for j=1:length(convs)-1
         if ~any(all(edges==convs(j) | edges==convs(j+1), 2))
 
           [d, p] = distancePointEdge(tmp, [pts(convs(j),:) pts(convs(j+1),:)]);
           d(d==0) = NaN;
           [v,dindx] = min(d);
           midpt = pcurr(dindx);
 
           p1 = grShortestPath(pts, edges, convs(j), midpt, dist);
           p2 = grShortestPath(pts, edges, midpt, convs(j+1), dist);
 
           full = [full; p1(1:end-1); p2(1:end-1)];
         else
           full = [full; convs(j)];
         end
       end
     end
 
     periph = [periph; pts(full([1:end 1]),:); NaN(1,2);];
   end
 
   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