function [lH_image, rH_image, lH_mask, rH_mask, lH_scouts, rH_scouts] = flatten_source(source_path, surface_path, size_x, size_y, freq_name)

% From 3D spherical representation to 2D projection source images
%
% Code to project registered source images from a spherical representation 
% to a mollweide 2D projection representation.
%
% INPUT:
%       - source_path (char): path to the source image .mat file
%       - surface_path (char): path the to surface file .mat file
%       - size_x (integer): x dimension of the resulting image
%       - size_y (integer): y dimension of the resulting image
%       - freq_name (string): the desired frequency band. 
%       (delta, theta, alpha, beta1, beta2, gamma1, gamma2)
%
% OUTPUT:
%       - lH_image (double matrix): the resulting 2D projection of the left
%       hemisphere
%       - rH_image (double matrix): the resulting 2D projection of the
%       right hemisphere
%       - lH_mask (double matrix): matrix the same size as the image with
%       1's if the image contains a number and 0 if the value is NaN
%       - rH_mask (double matrix): matrix the same size as the image with
%       1's if the image contains a number and 0 if the value is NaN
%       - lH_scouts (structure): contains information about scouts called
%       Painareas
%
% FUNCTION DEPENDENCIES:
%       - tess_hemisplit: This function is part of the Brainstorm software:
%         https://neuroimage.usc.edu/brainstorm
%       - mollweide_proj: Code written by Saskia van Heumen
%
%  Written by: Saskia van Heumen, September 2019
%  Written for: Internship at Centre of Pain Medicine, Erasmus MC,
%  Rotterdam

    freq_bands = {'delta', 'theta', 'alpha', 'beta1', 'beta2', 'gamma1', 'gamma2'};
    freq_no = find(ismember(freq_bands, freq_name));
    
    if isempty(freq_no)
        error('Input frequency band is not correct');
    end

    % Load the values of the sourceimages
    sourcemap = load(source_path);

    % Load the surface file containing the spherical coordinates
    surface = load(surface_path);

    % Find the indexes of the right and left hemispheres
    [rH_idx, lH_idx, ~, ~, ~, ~] = tess_hemisplit(surface);

    % Get the values of the source image

    if isfield(sourcemap, 'ImageGridAmp')
        sourcevalues = sourcemap.ImageGridAmp(:, 1);
    elseif isfield(sourcemap, 'TF')
        sourcevalues = squeeze(sourcemap.TF);
    else
        error('The sourcemap structure must contain an ImageGridAmp or TF field');
    end

    % Get spherical vertices of the right and left hemisphere
    lHvertices = double(surface.Reg.Sphere.Vertices(lH_idx, :));
    rHvertices = double(surface.Reg.Sphere.Vertices(rH_idx, :));

    % Get the latidude and longitude coordinates of the vertices
    [lH_lon, lH_lat, lH_radius] = cart2sph(lHvertices(:, 1), ...
        lHvertices(:, 2), lHvertices(:, 3));
    [rH_lon, rH_lat, rH_radius] = cart2sph(rHvertices(:, 1), ...
        rHvertices(:, 2), rHvertices(:, 3));

    % Get the 2D projection coordinates using mollweide projection
    [lH_x, lH_y] = mollweide_proj(lH_lat, lH_lon, mean(lH_radius(:)));
    [rH_x, rH_y] = mollweide_proj(rH_lat, rH_lon, mean(rH_radius(:)));

    % Get a equally spaced coordinate system for the image (uniform grid)
    lH_xgrid = linspace(min(lH_x), max(lH_x), size_x);
    lH_ygrid = linspace(min(lH_y), max(lH_y), size_y);
    [lH_xnew, lH_ynew] = meshgrid(lH_xgrid, lH_ygrid);

    rH_xgrid = linspace(min(rH_x), max(rH_x), size_x);
    rH_ygrid = linspace(min(rH_y), max(rH_y), size_y);
    [rH_xnew, rH_ynew] = meshgrid(rH_xgrid, rH_ygrid);

    % Interpolate the scattered data on the uniform grid
    lH_image = griddata(lH_x, lH_y, sourcevalues(lH_idx, freq_no), ...
        lH_xnew, lH_ynew, 'linear');
    rH_image = griddata(rH_x, rH_y, sourcevalues(rH_idx, freq_no), ...
        rH_xnew, rH_ynew, 'linear');
    
    % Create the masks
    lH_mask = double(~isnan(lH_image));
    rH_mask = double(~isnan(rH_image));
    
    % Scouts
    for i = 1: length(surface.Atlas)
        name_atlas = surface.Atlas(i);
        if strcmp(name_atlas.Name, 'PainAreas')
            index = i;
            break
        end
    end

    % Initialise structures for scouts
    lH_scouts = struct;
    rH_scouts = struct;
    num_lH = 1;
    num_rH = 1;

    % Get the scouts names and vertices number and vertices coordinates and
    % projection coordinates for both hemispheres
    for i = 1: length(surface.Atlas(index).Scouts)
        scouts_vertices_no = surface.Atlas(index).Scouts(i).Vertices;
        
        if mean(ismember(scouts_vertices_no, lH_idx)) == 1
            lH_scouts(num_lH).name = surface.Atlas(index).Scouts(i).Label;
            lH_scouts(num_lH).vertices_no = surface.Atlas(index).Scouts(i).Vertices;
            lH_scouts(num_lH).vertices = double(surface.Reg.Sphere.Vertices ...
                (scouts_vertices_no, :));
            [~, idx] = ismember(surface.Atlas(index).Scouts(i).Vertices, lH_idx);
            lH_scouts(num_lH).projection = [lH_x(idx, :), lH_y(idx, :)];
            num_lH = num_lH + 1;
        elseif mean(ismember(scouts_vertices_no, rH_idx)) == 1
            rH_scouts(num_rH).name = surface.Atlas(index).Scouts(i).Label;
            rH_scouts(num_rH).vertices_no = surface.Atlas(index).Scouts(i).Vertices;
            rH_scouts(num_rH).vertices = double(surface.Reg.Sphere.Vertices ...
                (scouts_vertices_no, :));
            [~, idx] = ismember(surface.Atlas(index).Scouts(i).Vertices, rH_idx);
            rH_scouts(num_rH).projection = [rH_x(idx, :), rH_y(idx, :)];
            num_rH = num_rH + 1;
        else
            disp(strcat('Scout: ', {' '},  surface.Atlas(index).Scouts(i).Label , {' '},'covered both the left and the right hemisphere and was skipped'))
        end
        
    end

end


