% Developed by Dr. Sonia Dahdouh
% Modified by Dr. Yuan-Chiao Lu
% Last Modified Date: 1/26/2022
clear
clc
close all

% Goal: Visual inspection of the mesh

%% Preparation
% Create a folder called 'Data_Step2'
% Copy all results of Step 1 from 'Data_Step1' to 'Data_Step2'

%% Please change the file path
main_dir = '/data/mril/users/all/mrdata/research/processed/CNMC/chd_r01/fetal/wip/yuan/Final_Work/Brain/Cortical_features';
data_path=fullfile(main_dir,'Test_run/Data_Step2_2');

%% Main program starts here 
cd(main_dir)
addpath(genpath(fullfile(main_dir,'matlab_code')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','nifti')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','iso2mesh')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','peyre')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','toolbox_brainsuite','utilities')));
addpath(genpath(fullfile(main_dir,'matlab_toolbox','CaptureFigVid')))

cprintf('_blue','Step 2/3: Check Parcellation\n');

Sbj=dir(fullfile(data_path,'fetus*'));
csvfile = fullfile(data_path,'QA.csv');
selection_data=[];
color_map=jet(40); color_map_new=zeros(40,3);
for rand_itr=1:40
    color_map_new(1:20,:)=color_map(1:2:end,:);
    color_map_new(21:40,:)=color_map(2:2:end,:);
    color_map=color_map_new;
end

for Sbj_count=1:size(Sbj,1)

    Sbj_path=fullfile(data_path,Sbj(Sbj_count).name);
    
    nm = strfind(Sbj(Sbj_count).name,'_');
    Sbj_ID=Sbj(Sbj_count).name(nm+1:end-1);
    Scan_ID=Sbj(Sbj_count).name(end);
    cprintf('_blue',['##### Subject ID: ' num2str(Sbj_ID) ', and Scan: ' num2str(Scan_ID) ' #####\n']);
    cd(Sbj_path)

    parcellation_path = fullfile(Sbj_path,'parcellation');

    QA=[];
    list=[];
    
    cmd=['find ',parcellation_path,' -name ','*_WM_1.off'];
    [status,listing]= system(cmd);
    startIndex=regexp(listing,'\n'); 
    startIndex=[1 startIndex];

    for ii=1:(size(startIndex,2)-1)    

        meshFileName_1 = strtrim(listing(startIndex(ii):startIndex(ii+1)));
        [path_res,name,ext] = fileparts(meshFileName_1);
        nm = strfind(name,'_');
        fet_name = name(1:5);
        scan_name = name(7:8);
        list=[list;{fet_name scan_name}];
        
        meshFileName_2 = strtrim(fullfile(path_res,[name(1:(nm(1)-1)),'_WM_2.off']));

        labels_1 = strtrim(fullfile(path_res,[name(1:(nm(1)-1)),'_labels_1.csv']));
        labels_2 = strtrim(fullfile(path_res,[name(1:(nm(1)-1)),'_labels_2.csv']));

        [vertices1,faces1]=read_off(meshFileName_1);
        [vertices2,faces2]=read_off(meshFileName_2);

        labels1=csvread(labels_1);
        labels2=csvread(labels_2);

        % Rotate the hemispheres for plot
        vertices1=[vertices1(1,:);vertices1(2,:);-vertices1(3,:)];
        [vertices1,~,~] = AxelRot(vertices1, 180, [1 0 0], mean(vertices1'));
            
        vertices2=[vertices2(1,:);vertices2(2,:);-vertices2(3,:)];
        [vertices2,~,~] = AxelRot(vertices2, 180, [1 0 0], mean(vertices2')); 
        
        % Export plot
        plot_mesh_with_map(vertices1,faces1,labels1,color_map)
        plot_mesh_with_map(vertices2,faces2,labels2,color_map)

        %% Determine the quality of parcellation
        a=input('Are the parcellation satisfactory? 0/1/2 (yes/maybe/no)');
        if (a==0 || a==1 || a==2)
            QA=[QA;a];
        else
            disp 'You must input 0, 1 or 2. One last try before program stops.'
            a=input('Are the parcellation satisfactory? 0/1/2 (yes/maybe/no)');
            if(a==0 || a==1|| a==2)
                QA=[QA;a];
            else
                disp 'Input error. The program is stopping.';
                return;
            end
        end
        
        close all
    end
    selection_data=[selection_data;[list num2cell(QA)]];

    clearvars -except main_dir data_path Sbj Sbj_count csvfile selection_data color_map

end
sbj_ID=selection_data(:,1);
scan_ID=selection_data(:,2);
QA=selection_data(:,3);
T=table(sbj_ID,scan_ID,QA);
writetable(T,csvfile)
