function [grown_sym_patches, sym_patches, Isym_neigh_patches, Dsym, Dsym_max] = grow_patches_proj_sym(FV, sym_samples, Isym_vertex_samples, sym_indices, ratio, options)
%GROW_PATCHES_PROJ  Patch Growing by Projection.
%   [grown_sym_patches, sym_patches, Isym_neigh_patches, Dsym, Dsym_max] = grow_patches_proj_sym(FV, sym_samples, Isym_vertex_samples, sym_indices, ratio, options)
%   outputs symmetric 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 symmetry-aware 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.
%       * sym_samples: array of size #symmetric sample pairs x 2 where rows represent the symmetric indices of sampled vertices.
%       * Isym_vertex_samples: vector of length #symmetric vertex pairs where entries contain the closest symmetric sample indices for the respective symmetric vertex pairs.
%       * sym_indices: array of size #symmetric vertex pairs x 2 where rows represent the symmetric vertex indices for the respective symmetric vertex pairs.
%       * 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 symmetric patch growing (default: false).
%
%   Output arguments:
%       * grown_sym_patches: cell of length #symmetric sample pairs where entries are 2-column arrays containing the vertex indices that are included in the respective grown symmetric patches;
%       * sym_patches: cell of length #symmetric sample pairs where entries are 2-column arrays containing the vertex indices that are included in the respective ungrown symmetric patches;
%       * Isym_neigh_patches: cell of length #symmetric sample pairs where entries are 2-column arrays containing the indices of symmetric patches that are neighbors of the respective symmetric patches (first column) and associated types (0: both sides for self-symmetric patches, 1 or 2: same or opposite sides for non-self-symmetric patches) of neighboring relationship (second column);
%       * Dsym: array of size #vertices x 2 x #symmetric sample pairs where rows represent the distances of the respective vertices to the different symmetric patches (infinity if greater than threshold);
%       * Dsym_max: array of size #vertices x 2 x #symmetric sample pairs where rows represent the thresholds for maximum distances of the respective vertices to the different symmetric patches.
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Options
if nargin < 6
    options = struct;
end

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

% Get the number of symmetric patches
Nsym_patches = size(sym_samples, 1);

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

% Build the symmetry list
sym_list                    = zeros(size(FV.vertices, 1), 1);
sym_list(sym_indices(:, 1)) = sym_indices(:, 2);
sym_list(sym_indices(:, 2)) = sym_indices(:, 1);

% Initialization
Dsym                = inf(Nvertices, 2, Nsym_patches);
Dsym_max            = -inf(Nvertices, 2, Nsym_patches);
sym_patches         = cell(1, Nsym_patches);
grown_sym_patches   = cell(1, Nsym_patches);
Isym_neigh_patches  = cell(1, Nsym_patches);

% Build the symmetric patches
for p = 1: Nsym_patches
    
    % Find the symmetric vertices
    sym_patches{p} = sym_indices(Isym_vertex_samples == p, :);
    
    % Fix the self-symmetric patches
    if sym_samples(p, 1) == sym_samples(p, 2)
        sym_patches{p} = unique(cat(1, sym_patches{p}, sym_patches{p}(:, 2: -1: 1)), 'rows');
    end
    
    % Find the faces for the first symmetric part
    faces_onein     = FV.faces(any(ismember(FV.faces, sym_patches{p}(:, 1)), 2), :);
    faces_allin 	= FV.faces(all(ismember(FV.faces, sym_patches{p}(:, 1)), 2), :);
    faces_oneout    = setdiff(faces_onein, faces_allin, 'rows');
    
    % Find the neighboring vertices for the first symmetric part
    tmp_indices = setdiff(unique(faces_oneout(:)), sym_patches{p}(:, 1));
    
    % Find the neighboring patches in the first symmetric part
    [~, tmp_is]     	= intersect(sym_indices(:, 1), tmp_indices);
    I_neigh_patches1 	= unique(Isym_vertex_samples(tmp_is));
    
    % Find the neighboring patches in the second symmetric part
    [~, tmp_is]         = intersect(sym_indices(:, 2), tmp_indices);
    I_neigh_patches2	= unique(Isym_vertex_samples(tmp_is));
    
    % Fix the self-symmetric neighboring patches
    I_neigh_patches0 = I_neigh_patches1(sym_samples(I_neigh_patches1, 1) == sym_samples(I_neigh_patches1, 2));
    I_neigh_patches1 = setdiff(I_neigh_patches1, I_neigh_patches0);
    I_neigh_patches2 = setdiff(I_neigh_patches2, I_neigh_patches0);
    
    % Build the neighboring symmetric patches
    Isym_neigh_patches{p}   = [I_neigh_patches0; I_neigh_patches1; I_neigh_patches2];
    tmp_col                 = [ zeros(length(I_neigh_patches0), 1); ...
                                 ones(length(I_neigh_patches1), 1); ...
                                 twos(length(I_neigh_patches2), 1)];
    Isym_neigh_patches{p}   = [Isym_neigh_patches{p}, tmp_col];
    
end

clear faces_allin faces_onein faces_oneout tmp_indices tmp_is tmp_col I_neigh_patches0 I_neigh_patches1 I_neigh_patches2

% Grow the symmetric patches
if ratio <= 0
    grown_sym_patches = sym_patches;
    return
end

for p = 1: Nsym_patches
    
    if verbose
        disp(['Growing symmetric patch number ' num2str(p) '.']);
    end
    
    % Get the symmetric patch
    sym_patch = sym_patches{p};
    
    % Update the grown symmetric patch
    grown_sym_patches{p} = sym_patch;
    
    % If the patch is self-symmetric
    if sym_samples(p, 1) == sym_samples(p, 2)
        
        % For each neighboring symmetric patch
        for np = 1: size(Isym_neigh_patches{p}, 1)
            
            % Prune the symmetric neighbors that have already been treated
            if Isym_neigh_patches{p}(np, 2) == 2
                break
            end
            
            % Get the neighboring symmetric patch
            neigh_patch = sym_patches{Isym_neigh_patches{p}(np, 1)};
                        
            % Build the symmetric constraint map
            L                       = -inf(size(FV.vertices, 1), 2);
            L(sym_patch(:, 1), 1)   = +inf;
            L(neigh_patch(:, 1), 1) = +inf;
            L(sym_patch(:, 2), 2)   = +inf;
            L(neigh_patch(:, 2), 2) = +inf;
            
            % Compute the maximum symmetric geodesic distance (reuse if already computed)
            dmax = Dsym_max(sym_samples(p, 1), 1, Isym_neigh_patches{p}(np, 1));
            if dmax == -inf
                if sym_samples(Isym_neigh_patches{p}(np, 1), 1) == sym_samples(Isym_neigh_patches{p}(np, 1), 2)
                    d       = fast_marching_mesh(FV, sym_samples(p, 1), sym_samples(Isym_neigh_patches{p}(np, 1), 1), L(:, 1));
                    dmax    = ratio * sqrt(2) * d(sym_samples(Isym_neigh_patches{p}(np, 1), 1));
                else
                    d1      = fast_marching_mesh(FV, sym_samples(p, 1), sym_samples(Isym_neigh_patches{p}(np, 1), 1), L(:, 1));
                    d2      = fast_marching_mesh(FV, sym_samples(p, 2), sym_samples(Isym_neigh_patches{p}(np, 1), 2), L(:, 2));
                    dmax    = ratio * sqrt(                                         ...
                              d1(sym_samples(Isym_neigh_patches{p}(np, 1), 1))^2 +	...
                              d2(sym_samples(Isym_neigh_patches{p}(np, 1), 2))^2  	...
                              );
                end
            end
            
            % Build the constraint map
            L                       = -inf(size(FV.vertices, 1), 2);
            L(neigh_patch(:, 1), 1) = dmax;
            
            % Compute the symmetric geodesic distances to the symmetric patch
            ds1                     = fast_marching_mesh(FV, sym_patch(:, 1), [], L(:, 1));
            L(neigh_patch(:, 2), 2) = sqrt(max(0, dmax^2 - ds1(neigh_patch(:, 1)).^2));
            ds2                     = fast_marching_mesh(FV, sym_patch(:, 2), [], L(:, 2));
            
            % Update the grown patch
            tmp_indices             = find(isfinite(ds2));
            grown_sym_patches{p}    = cat(1, grown_sym_patches{p}, [sym_list(tmp_indices), tmp_indices]);
            
            % Update the distances
            dsym                                = sqrt(ds1(neigh_patch(:, 1)).^2 + ds2(neigh_patch(:, 2)).^2);
            Dsym(neigh_patch(:, 1), 1, p)       = dsym;
            Dsym(neigh_patch(:, 2), 2, p)       = dsym;
            Dsym_max(neigh_patch(:, 1), 1, p)   = dmax;
            Dsym_max(neigh_patch(:, 2), 2, p)   = dmax;
            
        end
        
        % Fix the self-symmetry
        grown_sym_patches{p} = unique(cat(1, grown_sym_patches{p}, grown_sym_patches{p}(:, 2: -1: 1)), 'rows');
        
    % Otherwise
    else
        
        % For each neighboring symmetric patch
        for np = 1: size(Isym_neigh_patches{p}, 1)
            
            % Get the neighboring symmetric patch
            if Isym_neigh_patches{p}(np, 2) == 2
                neigh_patch = sym_patches{Isym_neigh_patches{p}(np, 1)}(:, 2: -1: 1);
            else
                neigh_patch = sym_patches{Isym_neigh_patches{p}(np, 1)};
            end
            
            % Build the symmetric constraint map
            L                       = -inf(size(FV.vertices, 1), 2);
            L(sym_patch(:, 1), 1)   = +inf;
            L(neigh_patch(:, 1), 1) = +inf;
            L(sym_patch(:, 2), 2)   = +inf;
            L(neigh_patch(:, 2), 2) = +inf;
            
            % Compute the maximum symmetric geodesic distance (reuse if already computed)
            dmax = Dsym_max(sym_samples(p, 1), max(1, Isym_neigh_patches{p}(np, 2)), Isym_neigh_patches{p}(np, 1));
            if dmax == -inf
                if sym_samples(p, 2) == sym_samples(Isym_neigh_patches{p}(np, 1), max(1, Isym_neigh_patches{p}(np, 2)))
                    d       = fast_marching_mesh(FV, sym_samples(p, 1), sym_samples(p, 2), L(:, 1));
                    dmax    = ratio * sqrt(2) * d(sym_samples(p, 2));
                else
                    d1      = fast_marching_mesh(FV, sym_samples(p, 1), sym_samples(Isym_neigh_patches{p}(np, 1), max(1,     Isym_neigh_patches{p}(np, 2))), L(:, 1));
                    d2      = fast_marching_mesh(FV, sym_samples(p, 2), sym_samples(Isym_neigh_patches{p}(np, 1), min(2, 3 - Isym_neigh_patches{p}(np, 2))), L(:, 2));
                    dmax    = ratio * sqrt(                                                                                	...
                              d1(sym_samples(Isym_neigh_patches{p}(np, 1), max(1,     Isym_neigh_patches{p}(np, 2))))^2 +	...
                              d2(sym_samples(Isym_neigh_patches{p}(np, 1), min(2, 3 - Isym_neigh_patches{p}(np, 2))))^2  	...
                              );
                end
            end
            
            % Build the constraint map
            L                       = -inf(size(FV.vertices, 1), 2);
            L(neigh_patch(:, 1), 1) = dmax;
            
            % Compute the symmetric geodesic distances to the symmetric patch
            ds1                     = fast_marching_mesh(FV, sym_patch(:, 1), [], L(:, 1));
            L(neigh_patch(:, 2), 2) = sqrt(max(0, dmax^2 - ds1(neigh_patch(:, 1)).^2));
            ds2                     = fast_marching_mesh(FV, sym_patch(:, 2), [], L(:, 2));
            
            % Update the grown patch
            tmp_indices             = find(isfinite(ds2));
            grown_sym_patches{p}    = cat(1, grown_sym_patches{p}, [sym_list(tmp_indices), tmp_indices]);
            
            % Update the distances
            dsym                                = sqrt(ds1(neigh_patch(:, 1)).^2 + ds2(neigh_patch(:, 2)).^2);
            Dsym(neigh_patch(:, 1), 1, p)       = dsym;
            Dsym(neigh_patch(:, 2), 2, p)       = dsym;
            Dsym_max(neigh_patch(:, 1), 1, p)   = dmax;
            Dsym_max(neigh_patch(:, 2), 2, p)   = dmax;
            
        end
        
        % Fix the common vertices
        grown_sym_patches{p} = unique(grown_sym_patches{p}, 'rows');
        
    end
    
    % Update the distances
    Dsym(sym_patch(:, 1), 1, p)     = 0;
    Dsym(sym_patch(:, 2), 2, p)     = 0;
    Dsym_max(sym_patch(:, 1), 1, p) = 0;
    Dsym_max(sym_patch(:, 2), 2, p) = 0;
    
end

end

function x = twos(varargin)

x = 2 * ones(varargin{:});

end