function [sym_samples, Isym_vertex_samples] = segment_mesh_uni_sym(FV, Nsamples, sym_indices, options)
%SEGMENT_MESH_UNI_SYM   Symmetric Uniform Mesh Segmentation.
%   [sym_samples, Isym_vertex_samples] = segment_mesh_uni_sym(FV, Nsamples, sym_indices, options)
%   outputs a symmetric uniform segmentation of a triangular mesh by 
%   greedily selecting a specified number of vertices via an iterative 
%   farthest point strategy based on geodesic distances, and by computing 
%   the associated Voronoi tessellation.
%
%   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.
%       * Nsamples: positive integer representing the number of vertices to sample (can be set to vector of length 2 with same entries for number of symmetric vertex pairs to sample).
%       * sym_indices: array of size #symmetric vertex pairs x 2 where rows represent the symmetric vertex indices for the respective symmetric vertex pairs.
%       * options: structure of options with the following fields:
%           - seeds: array of size #symmetric seed pairs x 2 where rows represent the symmetric vertex indices used for seeding (default: empty);
%           - verbose: boolean to display the current symmetric sampling iterations (default: false).
%
%   Output arguments:
%       * sym_samples: array of size #symmetric sample pairs x 2 where rows represent the symmetric indices of sampled vertices sorted in sampling order (first rows coincide with seeds);
%       * Isym_vertex_samples: array of size #symmetric vertex pairs x #symmetric sample pairs where rows contain the closest symmetric sample indices at successive iterations for the respective symmetric vertex pairs (assumes no ties).
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Number of vertices
Nvertices       = size(FV.vertices, 1);
Nsym_vertices   = size(sym_indices, 1);

% Number of samples
switch numel(Nsamples)
    case 1
        Nsym_samples = Nsamples; % higher bound
    case 2
        if Nsamples(1) ~= Nsamples(2)
            error('The number of samples should be the same on both sides.');
        end
        Nsym_samples    = Nsamples(1);
        Nsamples        = 2 * Nsym_samples; % higher bound
    otherwise
        error('The number of samples should be of length 1 or 2.');
end

% Options
if nargin < 4
    options = struct;
end

% Symmetric seeds
sym_seeds = [];
if isfield(options, 'sym_seeds')
    sym_seeds = options.sym_seeds;
end

% Number of seeds
Nseeds      = length(unique(sym_seeds(:)));
Nsym_seeds  = size(sym_seeds, 1);
if Nseeds > Nsamples || Nsym_seeds > Nsym_samples
    error('MATLAB:samplingMoreSeedsThanSamples', 'The number of seeds should be less than the number of samples.');
end

% Distance between symmetric vertices
if isfield(options, 'Dsym_vertex') && ~isempty(options.Dsym_vertex)
    Dsym_vertex                                         = sqrt(2) * options.Dsym_vertex;
    Dsym_vertex(sym_indices(:, 1) == sym_indices(:, 2)) = inf;
else
    Dsym_vertex                                         = compute_dist_sym(FV, sym_indices, options);
    Dsym_vertex                                         = sqrt(2) * Dsym_vertex;
    Dsym_vertex(sym_indices(:, 1) == sym_indices(:, 2)) = inf;
end

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

% Initialization
sym_samples       	= zeros(Nsym_samples, 2);               % symmetric samples
Dsym_vertex_samples	= inf(Nsym_vertices, 1);                % distance of each symmetric vertex to the nearest symmetric sample
Isym_vertex_samples = zeros(Nsym_vertices, Nsym_samples);	% memory of the nearest symmetric sample for each symmetric vertex

% Build the symmetric constraint map
Lsym                     	= -inf(Nvertices, 2);
Lsym(sym_indices(:, 1), 1)  = +inf;
Lsym(sym_indices(:, 2), 2)  = +inf;

% Symmetric seeding
if Nsym_seeds
    
    % Sample the symmetric seeds
    for i = 1: Nsym_seeds
        
        if verbose
            disp(['Sampling symmetric point number ' num2str(i) '.']);
        end
        
        % Select the symmetric seed
        sym_samples(i, :) = sym_seeds(i, :);
        
        % Compute the distance of each symmetric vertex to the last symmetric sample
        Lsym(sym_indices(:, 1), 1)  = Dsym_vertex_samples;
        Dsym_vertex_lsample(:, 1)   = fast_marching_mesh(FV, sym_samples(i, 1), [], Lsym(:, 1));
        Lsym(sym_indices(:, 2), 2)	= sqrt(max(0, Dsym_vertex_samples.^2 - Dsym_vertex_lsample(sym_indices(:, 1), 1).^2));
        Dsym_vertex_lsample(:, 2)   = fast_marching_mesh(FV, sym_samples(i, 2), [], Lsym(:, 2));
        
        % Update the distance and nearest symmetric sample for each symmetric vertex
        tmp_distances                               = sqrt(Dsym_vertex_lsample(sym_indices(:, 1), 1).^2 + Dsym_vertex_lsample(sym_indices(:, 2), 2).^2);
        tmp_indices                                 = find(tmp_distances < Dsym_vertex_samples);
        Dsym_vertex_samples(tmp_indices)            = tmp_distances(tmp_indices);
        Isym_vertex_samples(tmp_indices, i: end)    = i;
        
    end 
       
else
    
    if verbose
        disp('Sampling symmetric point number 1.');
    end
    
    % Sample the first symmetric point
    sym_samples(1, :) = sym_indices(randi(Nsym_vertices), :);
        
    % Replace it with one of the farthest symmetric vertices
    Dsym_vertex_lsample(:, 1)   = fast_marching_mesh(FV, sym_samples(1, 1), [], Lsym(:, 1));
    Dsym_vertex_lsample(:, 2)   = fast_marching_mesh(FV, sym_samples(1, 2), [], Lsym(:, 2));
    tmp_distances               = min(Dsym_vertex, sqrt(Dsym_vertex_lsample(sym_indices(:, 1), 1).^2 + Dsym_vertex_lsample(sym_indices(:, 2), 2).^2));
    tmp_indices                 = find(tmp_distances == max(inf2zero(tmp_distances))); % due to numerical issues in fast marching, we need to remove the vertices that appear to be at infinite distance even if the mesh is connected
    tmp_length                	= length(tmp_indices);
    if tmp_length == 1
        sym_samples(1, :) = sym_indices(tmp_indices, :);
    else
        sym_samples(1, :) = sym_indices(tmp_indices(randi(tmp_length)), :);
    end
    
    % Compute the distance of each symmetric vertex to the last symmetric sample
    Dsym_vertex_lsample(:, 1) = fast_marching_mesh(FV, sym_samples(1, 1), [], Lsym(:, 1));
    Dsym_vertex_lsample(:, 2) = fast_marching_mesh(FV, sym_samples(1, 2), [], Lsym(:, 2));
    
    % Update the distance and nearest symmetric sample for each symmetric vertex
    Dsym_vertex_samples     = sqrt(Dsym_vertex_lsample(sym_indices(:, 1), 1).^2 + Dsym_vertex_lsample(sym_indices(:, 2), 2).^2);
    Isym_vertex_samples(:)	= 1;
   
    % Update the number of symmetric seeds
    Nsym_seeds = 1;

end

% Check end of sampling
samples = sym_samples(1: Nsym_seeds, :);
if length(unique(samples(:))) >= Nsamples
    sym_samples(Nsym_seeds + 1: Nsym_samples, :)            = [];
    Isym_vertex_samples(:, Nsym_seeds + 1: Nsym_samples)    = [];
    return
end

% Next symmetric points
for i = Nsym_seeds + 1: Nsym_samples
    
    if verbose
        disp(['Sampling symmetric point number ' num2str(i) '.']);
    end

    % Select one of the farthest symmetric vertices from the symmetric samples
    tmp_distances   = min(Dsym_vertex, Dsym_vertex_samples);
    tmp_indices     = find(tmp_distances == max(inf2zero(tmp_distances))); % due to numerical issues in fast marching, we need to remove the vertices that appear to be at infinite distance even if the mesh is connected
    tmp_length      = length(tmp_indices);
    if tmp_length == 1
        sym_samples(i, :) = sym_indices(tmp_indices, :);
    else
        sym_samples(i, :) = sym_indices(tmp_indices(randi(tmp_length)), :);
    end
       
    % Compute the distance of each symmetric vertex to the last symmetric sample
    Lsym(sym_indices(:, 1), 1)  = Dsym_vertex_samples;
    Dsym_vertex_lsample(:, 1)	= fast_marching_mesh(FV, sym_samples(i, 1), [], Lsym(:, 1));
    Lsym(sym_indices(:, 2), 2)  = sqrt(max(0, Dsym_vertex_samples.^2 - Dsym_vertex_lsample(sym_indices(:, 1), 1).^2));
    Dsym_vertex_lsample(:, 2)   = fast_marching_mesh(FV, sym_samples(i, 2), [], Lsym(:, 2));
    
    % Update the distance and nearest symmetric sample for each symmetric vertex
    tmp_distances                               = sqrt(Dsym_vertex_lsample(sym_indices(:, 1), 1).^2 + Dsym_vertex_lsample(sym_indices(:, 2), 2).^2);
    tmp_indices                                 = find(tmp_distances < Dsym_vertex_samples);
    Dsym_vertex_samples(tmp_indices)            = tmp_distances(tmp_indices);
    Isym_vertex_samples(tmp_indices, i: end)    = i;
    
    % Check end of sampling
    samples = sym_samples(1: i, :);
    if length(unique(samples(:))) >= Nsamples
        sym_samples(i + 1: Nsym_samples, :)         = [];
        Isym_vertex_samples(:, i + 1: Nsym_samples) = [];
        return
    end
    
end

end

function y = inf2zero(x)

y           = x;
y(isinf(y)) = 0;

end