% 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
    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
      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
  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

% 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