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

% Goal: Export cortical features to Excel file

%% Preparation
% Create a folder called 'Data_Step4'

%% 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,'Attempt3/Data_Step1');
result_path=fullfile(main_dir,'Attempt3/Data_Step4');

%% Main program starts here
% mat_path=fullfile(result_path,'Mat');
% if ~exist(mat_path)
%     mkdir(mat_path)
% end
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','GEEQBOX')));

%% Data cleaning
Sbj_info=[];
Mean_Cmean=[];
Mean_Cgauss=[];
Mean_Curvedness=[];
Mean_lGI=[];
Mean_SulcalDepth=[];
Msur_vol=[];
Msur_area=[];

Sbj=dir(fullfile(data_path,'fetus*'));

for Sbj_count=1:size(Sbj,1)

    Sbj_path=fullfile(data_path,Sbj(Sbj_count).name);
    cd(Sbj_path)    
    parcellation_path=fullfile(Sbj_path,'parcellation');
    segmentation_path=fullfile(Sbj_path,'segmentation');

    nm = strfind(Sbj(Sbj_count).name,'_');
    Sbj_ID=Sbj(Sbj_count).name(nm+1:end-1);
    Sbj_ID_str=Sbj_ID;

    Scan_ID=Sbj(Sbj_count).name(end);
    cprintf('_blue',['##### Subject ID: ' num2str(Sbj_ID_str) ', and Scan: ' num2str(Scan_ID) ' #####\n']);
      
    sbj_info={Sbj_ID_str str2num(Scan_ID)};

    % Plot cortical features
%     feature_plot(parcellation_path,Sbj_ID,Scan_ID)

    try
        % Calculate mean cortical features
        [Mean_Cmean_1,Mean_Cgauss_1,Mean_Curvedness_1,Mean_lGI_1,Mean_SulcalDepth_1,vol_1,area_1]=...
            feature_calc(parcellation_path,Sbj_ID_str,Scan_ID,1); % 1: Right hemisphere
        [Mean_Cmean_2,Mean_Cgauss_2,Mean_Curvedness_2,Mean_lGI_2,Mean_SulcalDepth_2,vol_2,area_2]=...
            feature_calc(parcellation_path,Sbj_ID_str,Scan_ID,2); % 2: Left hemisphere
        [Mean_Cmean_3,Mean_Cgauss_3,Mean_Curvedness_3,Mean_lGI_3,Mean_SulcalDepth_3,vol_3,area_3]=...
            feature_calc(parcellation_path,Sbj_ID_str,Scan_ID,3); % 3: Both hemisphere

        Mean_Cmean=[Mean_Cmean;[sbj_info num2cell(Mean_Cmean_1) num2cell(Mean_Cmean_2) num2cell(Mean_Cmean_3)]];
        Mean_Cgauss=[Mean_Cgauss;[sbj_info num2cell(Mean_Cgauss_1) num2cell(Mean_Cgauss_2) num2cell(Mean_Cgauss_3)]];
        Mean_Curvedness=[Mean_Curvedness;[sbj_info num2cell(Mean_Curvedness_1) num2cell(Mean_Curvedness_2) num2cell(Mean_Curvedness_3)]];
        Mean_lGI=[Mean_lGI;[sbj_info num2cell(Mean_lGI_1) num2cell(Mean_lGI_2) num2cell(Mean_lGI_3)]];
        Mean_SulcalDepth=[Mean_SulcalDepth;[sbj_info num2cell(Mean_SulcalDepth_1) num2cell(Mean_SulcalDepth_2) num2cell(Mean_SulcalDepth_3)]];
        Msur_vol=[Msur_vol;[sbj_info num2cell(vol_1) num2cell(vol_2) num2cell(vol_3)]];
        Msur_area=[Msur_area;[sbj_info num2cell(area_1) num2cell(area_2) num2cell(area_3)]];    
    catch ME
        disp 'Segmentation or parcellation labels are not correct!'
    end

end

%% export_to_excel_noGA and save the data
filename=fullfile(result_path,'Cortical_features.xls');
export_to_excel(filename,Mean_Cmean,'Cmean')
export_to_excel(filename,Mean_Cgauss,'Cgauss')
export_to_excel(filename,Mean_Curvedness,'Curvedness')
export_to_excel(filename,Mean_lGI,'lGI')
export_to_excel(filename,Mean_SulcalDepth,'SulcalDepth')
export_to_excel(filename,Msur_vol,'Volume')
export_to_excel(filename,Msur_area,'Area')
