function computeGI(Mesh_name,Mesh_name_CH,savePath_lGI,savePath_area_mean) %,savePath_area_init,savePath_area_smooth,savePath_gGI

%==================================================================
% Read the pial surface
%==================================================================
[vertex,face] = read_off(Mesh_name);

%==================================================================
% Read the Convex Hull surface
%==================================================================
[vertex_CH,face_CH]=read_off(Mesh_name_CH);

%==================================================================
%Compute the local gyrification index
%==================================================================

%compute area of each face for the initial mesh
face2=face.';
vertex2 = vertex.';

v0 = vertex2(face2(:,1),:).';
v1 = vertex2(face2(:,2),:).';
v2 = vertex2(face2(:,3),:).';
e1 = (v1-v0);
e2 = (v2-v0);

area_init = 0.5.*sqrt(dot(cross(e1,e2),cross(e1,e2)));
area_init =area_init.';

%compute area of each face for the smooth mesh
face2=face_CH.';
Xmu2 = vertex_CH.';

v0 = Xmu2(face2(:,1),:).';
v1 = Xmu2(face2(:,2),:).';
v2 = Xmu2(face2(:,3),:).';

e1 = (v1-v0);
e2 = (v2-v0);
area_smooth = 0.5.*sqrt(dot(cross(e1,e2),cross(e1,e2)));
area_smooth =area_smooth.';

radius=2.7.^((sqrt(sum(area_smooth)/(4*pi))/10)-1);

% Calculate areas
vertex_area_smooth=zeros(size(vertex,2),1);
vertex_area_pial=zeros(size(vertex,2),1);
vertex_area_mean_pial=zeros(size(vertex,2),1);

for seed=1:size(vertex,2)
    connected=find(vecnorm((vertex'-vertex(:,seed)')')<radius);
    face_idx=find(sum(ismember(face,connected))==3);

    [~,closest_idx]=min(vecnorm((vertex_CH'-vertex(:,seed)')'));
    connected=find(vecnorm((vertex_CH'-vertex_CH(:,closest_idx)')')<radius);
    face_idx_CH=find(sum(ismember(face_CH,connected))==3); 

    res=sum(area_init(face_idx));
    res_smooth = sum(area_smooth(face_idx_CH));
    area_mean=mean(area_init(face_idx));
   
    vertex_area_smooth(seed) = res_smooth;
    vertex_area_pial(seed) = res; 
    vertex_area_mean_pial(seed) = area_mean; 
end

lGI = vertex_area_pial./vertex_area_smooth;
lGI(isnan(lGI))=0;
save(savePath_lGI,'lGI');

area_mean=vertex_area_mean_pial;
area_mean(isnan(lGI))=0;
save(savePath_area_mean,'area_mean');

end

