function [grown_patches, patches, I_neigh_patches, D, Dmax] = grow_patches_proj(FV, samples, I_vertex_samples, ratio, options)
%GROW_PATCHES_PROJ  Patch Growing by Projection.
%   [grown_patches, patches, I_neigh_patches, D, Dmax] = grow_patches_proj(FV, samples, I_vertex_samples, ratio, options)
%   outputs grown patches by projection of vertices from neighbor patches  
%   with respect to geodesic distances on a triangular mesh. The threshold  
%   to include a vertex from a neighbor patch to grow a given patch is 
%   given by the distance between patch centers scaled by an overlap ratio.
%
%   Input arguments:
%       * FV: structure representing a mesh with the following fields:
%           - faces: array of size #faces x 3 where rows contain the 3 vertex indices of the respective faces;
%           - vertices: array of size #vertices x 3 where rows contain the 3 spatial coordinates of the respective vertices.
%       * samples: vector of length #samples where entries represent the indices of sampled vertices.
%       * I_vertex_samples: vector of length #vertices where entries contain the closest sample indices for the respective vertices.
%       * ratio: non-negative number representing the overlap ratio (can be set to infinity).
%       * options: structure of options with the following fields:
%           - verbose: boolean to display the current patch growing (default: false).
%
%   Output arguments:
%       * grown_patches: cell of length #samples where entries are vectors containing the vertex indices that are included in the respective grown patches;
%       * patches: cell of length #samples where entries are vectors containing the vertex indices that are included in the respective ungrown patches;
%       * I_neigh_patches: cell of length #samples where entries are vectors containing the indices of patches that are neighbors of the respective patches;
%       * D: array of size #vertices x #samples where rows represent the distances of the respective vertices to the different patches (infinity if greater than threshold);
%       * Dmax: array of size #vertices x #samples where rows represent the thresholds for maximum distances of the respective vertices to the different patches.
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Options
if nargin < 5
    options = struct;
end

% Verbose
verbose = false;
if isfield(options, 'verbose')
    verbose = options.verbose;
end

% Get the number of patches
Npatches = length(samples);

% Get the number of vertices
Nvertices = size(FV.vertices, 1);

% Initialization
D               = inf(Nvertices, Npatches);
Dmax            = -inf(Nvertices, Npatches);
patches         = cell(1, Npatches);
grown_patches   = cell(1, Npatches);
I_neigh_patches = cell(1, Npatches);

% Build the patches
for p = 1: Npatches
    
    % Find the vertices
    patches{p} = find(I_vertex_samples == p);
    
    % Find the faces
    faces_onein     = FV.faces(any(ismember(FV.faces, patches{p}), 2), :);
    faces_allin     = FV.faces(all(ismember(FV.faces, patches{p}), 2), :);
    faces_oneout    = setdiff(faces_onein, faces_allin, 'rows');
    
    % Find the neighboring patches
    I_neigh_patches{p} = unique(I_vertex_samples(setdiff(unique(faces_oneout(:)), patches{p})));
    
end

clear faces_allin faces_onein faces_oneout

% Grow the patches
if ratio <= 0
    grown_patches = patches;
    return
end

for p = 1: Npatches
    
    if verbose
        disp(['Growing patch number ' num2str(p) '.']);
    end
    
    % Get the patch
    patch = patches{p};
    
    % Update the grown patch
    grown_patches{p} = patch;
    
    % For each neighboring patch
    for np = 1: length(I_neigh_patches{p})
        
        % Get the neighboring patch
        neigh_patch = patches{I_neigh_patches{p}(np)};
        
        % Build the constraint map
        L           	= -inf(size(FV.vertices, 1), 1);
        L(patch)      	= +inf;
        L(neigh_patch)  = +inf;
        
        % Compute the maximum geodesic distance (reuse if already computed)
        dmax = Dmax(samples(p), I_neigh_patches{p}(np));
        if dmax == -inf
            d   	= fast_marching_mesh(FV, samples(p), samples(I_neigh_patches{p}(np)), L);
            dmax    = ratio * d(samples(I_neigh_patches{p}(np)));
        end
        
        % Build the constraint map
        L               = -inf(size(FV.vertices, 1), 1);
        L(neigh_patch)  = dmax;
        
        % Compute the geodesic distances to the patch
        ds = fast_marching_mesh(FV, patch, [], L);
        
        % Update the grown patch
        grown_patches{p} = cat(1, grown_patches{p}, find(isfinite(ds)));
        
        % Update the distances
        D(neigh_patch, p)       = ds(neigh_patch);
        Dmax(neigh_patch, p)    = dmax;
        
    end
    
    % Fix the common vertices
    grown_patches{p} = unique(grown_patches{p});
    
    % Update the distances
    D(patch, p)     = 0;
    Dmax(patch, p)  = 0;
    
end

end