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

% Goal: Calculate cortical features
% This code use BrainSuite21a. If you use other versions of BrainSuite, please change the code in the funtion "convex_hull/reconstructAndLabelConvexHull"

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

%% 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');

%% Main program starts here 
matlab_dir = fullfile(main_dir,'matlab_code');
app_dir = '/Applications';
cd(main_dir)
addpath(genpath(fullfile(main_dir,'matlab_code')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','peyre')))
addpath(genpath(fullfile(main_dir,'matlab_toolbox','alecjacobson-gptoolbox-4b473f8')))

cprintf('_blue','Step 3/3: Compute Brain Meshes Indices\n');

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

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)

    parcellationPath = fullfile(Sbj_path,'parcellation');

    if ~exist(fullfile(parcellationPath,[Sbj_ID '-0' Scan_ID '_CH.nii'])) % CH: Convex Hull
        convex_hull(matlab_dir,app_dir,parcellationPath,Sbj_ID,Scan_ID);
    end

    for ii=1:2

        %Get path name per file and filename
        Mesh_name=fullfile(parcellationPath,[Sbj_ID '-0' Scan_ID '_WM_' num2str(ii) '.off']);
        Mesh_name_CH=fullfile(parcellationPath,[Sbj_ID '-0' Scan_ID '_CH_' num2str(ii) '.off']);
                
        try
            [pathstr,name,ext] = fileparts(Mesh_name);
            [node,fc]=read_off(Mesh_name);
            [FF,I] = orient_outward(node',fc',node'); % FF is "faces"
            
            [node_CH,fc_CH]=read_off(Mesh_name_CH);
            [FF_CH,I_CH] = orient_outward(node_CH',fc_CH',node_CH');         

            %% Calcualte curvatures
            FV.faces=fc';
            FV.vertices=node';
            getderivatives=0;
            [PrincipalCurvatures,PrincipalDir1,PrincipalDir2,FaceCMatrix,VertexCMatrix,Cmagnitude]= GetCurvatures( FV ,getderivatives);
            Cmean =((PrincipalCurvatures(1,:)+PrincipalCurvatures(2,:))/2)';
            Cgauss=(PrincipalCurvatures(1,:).*PrincipalCurvatures(2,:))';
            Curvedness=(sqrt((PrincipalCurvatures(1,:).^2 + PrincipalCurvatures(2,:).^2)/2))';

            % Avoid using any of the curvature measure except mean curvature, curvedness and gauss for clinical paper as not really streamline        
            savePath_foldingMeasures = strtrim(fullfile(pathstr,[name,'_foldingMeasures']));
            foldingMeasures.Cmean=Cmean;
            foldingMeasures.Cgauss=Cgauss;
            foldingMeasures.Curvedness=Curvedness;
            save(savePath_foldingMeasures,'foldingMeasures');

            %% Compute Local gyrification index for each vertex and mean surface area for each vertex
            savePath_lGI = strtrim(fullfile(pathstr,[name,'_lGI']));  
            savePath_area_mean = strtrim(fullfile(pathstr,[name,'_area_mean'])); 
            computeGI(Mesh_name,Mesh_name_CH,savePath_lGI,savePath_area_mean);

            %% Compute Sulcal Depth
            savePath_SulcalDepth = strtrim(fullfile(pathstr,[name,'_SulcalDepth']));  
            [Sulcal_Depth] = computeSulcalDepth(Mesh_name,Mesh_name_CH,savePath_SulcalDepth);
        catch ME
            fails{ii}=Mesh_name;
            disp 'error in one of the features; check your mesh';
        end
    end
    clearvars -except main_dir matlab_dir data_path Sbj Sbj_count app_dir lGI_radius

end

