% ----------------------------------------------------------------------- %
% Author: Mireia Leon Dasi
% European Space Agency - ESAC
% Last modified: 19/10/2020
% Description: matlab script to produce individual QuACK projections for       
% each region including the surface features
% ----------------------------------------------------------------------- %
clc; clear all; close all;

% Define save path 
%save_path='';

% Color definitions
Colors_beige = {[249,248,235]./255,[237,233,218]./255,[224,218,201]./255,...
                [212,203,184]./255,[200,188,168]./255};

%-------------------- LOAD REGION DEFINITIONS ----------------------------%
% Load triangulation matrix
load('triangulation_2M.mat');      % Triangulation matrix defining regions
T_reg = Triangulation;
load('triangulation_sub.mat');     % Triangulation matrix defining regions
T_sub = Triangulation;

% Load point mesh
fileID_xy = fopen('Mesh2M_out.txt','r');
points_xy = fscanf(fileID_xy,'%f %f %f',[3,Inf])';
x = points_xy(:,1);
y = points_xy(:,2);
z = zeros(length(x),1);

%-------------------- SELECT REGION(S) TO PLOT ---------------------------%
load('region_order.mat');
load('Colors.mat');
load('idx_sub.mat');

% Select region in array. Example: regions = {'Anuket','Apis'}
regions = {'Anuket'};

for k = 1:length(regions)
    % Determine if map contains features (Northern hemisphere, logic = 1) or 
    % the region has not been mapped yet (Southern hemisphere, logic = 0) 
    idx = ismember(region_order,regions{k});
    index = find(idx);
    logic = ismember(index,[1 3 4 5 7 9 10 13 14 15 18 19 22 23 25 26]);

    % Add directories
    if k == 26
        index = 27;
    end
    index_sub = idx_sub{index};

    % For mapped regions, read data from SBMT files

    if logic == 1

        % Load inout SBMT files
        Paths_input = strcat(regions{k},'_Paths');
        Ellipses_input = strcat(regions{k},'_Ellipses');

        % Open and read Input file
        fid1 = fopen(Paths_input);
        data_paths=textscan(fid1,'%s','delimiter','\n', 'HeaderLines', 2 );
        fclose(fid1);
        data_paths = data_paths{1,1};
        n_paths = length(data_paths)-1;

        fid2 = fopen(Ellipses_input);
        data_ellipses = textscan(fid2,'%d%s%f%f%f%f%f%f%s%s%s%s%f%f%f%s%s%s');
        fclose(fid2);

        % Extract path color, length and ellipse size data
        latlong = [];
        for i = 1:n_paths
            text = data_paths{i};

            % Read color
            [SC,EC] = regexp(text, 'color="\w*\W*\w*\W*\w*"');
            c = textscan( text(SC+7:EC-1), '%f', 'Delimiter',',' );
            c = c{1,1};
            colors(i,1:3) = c./255;

            %Read latitude and longitude
            [S,E] = regexp(text, 'vertices="');
            E = size(text,2);
            vertices = textscan( text(S+10:E-3), '%f', 'Delimiter',' ' );
            vertices = vertices{1,1};
            latlong = [latlong; reshape(vertices,3,size(vertices,1)/3)'];
            n_points(i,1) = size(vertices,1)/3;
        end
        path_coords = llr2cart(latlong')';
        ellipse_coords = [data_ellipses{1,3} data_ellipses{1,4} data_ellipses{1,5}];  

        % Read path and ellipse coordinates in QuACK layout

        paths_xy = readtable(strcat(regions{k},'_paths_q.txt'));
        paths_xy = table2array(paths_xy);

        opts = detectImportOptions(strcat(regions{k},'_ellipses_q.txt'));
        opts = setvartype(opts,3,'double');
        opts = setvaropts(opts,3,'DecimalSeparator','.');
        fid3 = fopen(strcat(regions{k},'_ellipses_q.txt'));
        ellipses_xy= textscan(fid3,'%f%f%f');
        fclose(fid3);

        ellipses_xy = [ellipses_xy{1} ellipses_xy{2}];
        n_ellipses = size(ellipses_xy,1);
        radius_ellipses = data_ellipses{1,13}/2;
        color_ellipses = data_ellipses{1,16};
    end
    %---------------------- PLOT GEOLOGICAL MAP(S) ---------------------------%
    f1 = figure;  

    hold on; axis off;axis equal

    % Plot region background
    if k ~=26
        trisurf(T_reg{1,index},x,y,z,'EdgeColor','k','FaceColor','k');
        trisurf(T_reg{1,index},x,y,z,'EdgeColor','none','FaceColor',Colors_beige{rem(k,5)+1},'FaceAlpha',1);

        trisurf(T_sub{1,index_sub(i)},x,y,z,'EdgeColor','k','FaceColor','k');
        trisurf(T_sub{1,index_sub(i)},x,y,z,'EdgeColor','none','FaceColor',Colors_beige{rem(k,5)+1},'FaceAlpha',1);
    end
    if k == 26
        for i = 1:length(index_sub)
            trisurf(T_sub{1,index_sub(i)},x,y,z,'EdgeColor','k','FaceColor','k');
            trisurf(T_sub{1,index_sub(i)},x,y,z,'EdgeColor','none','FaceColor',Colors_beige{rem(k,5)+1},'FaceAlpha',1);
        end
    end

    % Plot paths and ellipses if regions are mapped

    if logic == 1
        init = 1;
        for i = 1:n_paths
            if ((colors(i,1)~=0 )||(colors(i,2)~=0 )||(colors(i,3)~=0 ))
                plot(paths_xy(init:init+n_points(i)-1,1),paths_xy(init:init+n_points(i)-1,2),'Color',colors(i,:),'LineWidth',1)
            end
            init = init + n_points(i);
        end
        for i = 1:n_ellipses
           c =  textscan(color_ellipses{i},'%f', 'Delimiter',',');
           c = c{1,1};
           c = c';
           c = c./255;
           if ((c(1)~=0 )||(c(2)~=0 )||(c(3)~=0 ))
            circle(ellipses_xy(i,1),ellipses_xy(i,2),radius_ellipses(i)*0.11355,c);
           end
        end
    end
end

