function [samples, I_vertex_samples] = segment_mesh_uni(FV, Nsamples, options)
%SEGMENT_MESH_UNI	Uniform Mesh Segmentation.
%   [samples, I_vertex_samples] = segment_mesh_uni(FV, Nsamples, options)
%   outputs a 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.
%       * options: structure of options with the following fields:
%           - seeds: vector of length #seeds where entries represent the vertex indices used for seeding (default: empty);
%           - verbose: boolean to display the current sampling iterations (default: false).
%
%   Output arguments:
%       * samples: vector of length Nsamples where entries represent the indices of sampled vertices sorted in sampling order (first entries coincide with seeds);
%       * I_vertex_samples: array of size #vertices x Nsamples where rows contain the closest sample indices at successive iterations for the respective vertices (assumes no ties).
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

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

% Options
if nargin < 3
    options = struct;
end

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

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

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

% Initialization
samples             = zeros(Nsamples, 1);           % samples
D_vertex_samples	= inf(Nvertices, 1);            % distance of each vertex to the nearest sample
I_vertex_samples    = zeros(Nvertices, Nsamples);	% memory of the nearest sample for each vertex

% Seeding
if Nseeds
    
    % Sample the seeds
    for i = 1: Nseeds
        
        if verbose
            disp(['Sampling point number ' num2str(i) '.']);
        end
        
        % Select the seed
        samples(i) = seeds(i);
        
        % Compute the distance of each vertex to the last sample
        D_vertex_lsample = fast_marching_mesh(FV, samples(i), [], D_vertex_samples);
        
        % Update the distance and nearest sample for each vertex
        tmp_indices                         	= find(D_vertex_lsample < D_vertex_samples);
        D_vertex_samples(tmp_indices)           = D_vertex_lsample(tmp_indices);
        I_vertex_samples(tmp_indices, i: end)   = i;
        
    end
    
else
    
    if verbose
        disp('Sampling point number 1.');
    end
    
    % Sample the first point
    samples(1) = randi(Nvertices);
    
    % Replace it with one of the farthest vertices
    D_vertex_lsample	= fast_marching_mesh(FV, samples(1), [], []);
    tmp_indices         = find(D_vertex_lsample == max(inf2zero(D_vertex_lsample))); % 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
        samples(1) = tmp_indices;
    else
        samples(1) = tmp_indices(randi(tmp_length));
    end
    
    % Compute the distance of each vertex to the last sample
    D_vertex_lsample = fast_marching_mesh(FV, samples(1), [], []);
    
    % Update the distance and nearest sample for each vertex
    D_vertex_samples    = D_vertex_lsample;
    I_vertex_samples(:) = 1;
    
    % Update the number of seeds
    Nseeds = 1;

end

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

    % Select one of the farthest vertices from the samples
    tmp_indices	= find(D_vertex_samples == max(inf2zero(D_vertex_samples))); % 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
        samples(i) = tmp_indices;
    else
        samples(i) = tmp_indices(randi(tmp_length));
    end
    
    % Compute the distance of each vertex to the last sample
    D_vertex_lsample = fast_marching_mesh(FV, samples(i), [], D_vertex_samples);
    
    % Update the distance and nearest sample for each vertex
    tmp_indices                             = find(D_vertex_lsample < D_vertex_samples);
    D_vertex_samples(tmp_indices)           = D_vertex_lsample(tmp_indices);
    I_vertex_samples(tmp_indices, i: end)   = i;
    
end

end

function y = inf2zero(x)

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

end