function reconstructAndLabelCorticalSurface(reduced_parcellationPath,meshFileName_1,meshFileName_2,outputLabeFileName,matlab_dir,app_dir)

if(isempty(reduced_parcellationPath))
    disp 'You need to input the parcellation file name and path';
    return
end

[meshes_path,name,~] = fileparts(reduced_parcellationPath);

if (isempty(meshFileName_1) || isempty(meshFileName_2) || isempty(outputLabeFileName))
    
    meshFileName_1=strtrim(fullfile(meshes_path,[name,'_WM_1.off']));
    meshFileName_2=strtrim(fullfile(meshes_path,[name,'_WM_2.off']));        
    outputLabeFileName=strtrim(fullfile(meshes_path,[name,'_labels.csv']));        
end

scriptClean = fullfile(matlab_dir,'quadricEdgeCollapse.mlx');
image_data_nii = load_untouch_nii(reduced_parcellationPath);
image_labels = image_data_nii.img;
image_data = (image_labels>0);

%Determine hemisphere separation
[L,num] = bwlabeln(image_data);

max_1=0;
max_2=0;
max_1_ind=0;
max_2_ind=0;
if (num>1)
    for j=1:num
        if sum(L(:)==j)>=max_1
            max_2=max_1;
            max_2_ind=max_1_ind;
            max_1=sum(L(:)==j);
            max_1_ind=j;
        else
            if sum(L(:)==j)>=max_2
                max_2=sum(L(:)==j);
                max_2_ind=j;
            end
        end
    end
end

if( max_1>1000 && max_2>1000)
    
    im1 = L==max_1_ind;
    im2 = L==max_2_ind;
    [x1 y1 z1]=find(im1>0);
    [x2 y2 z2]=find(im2>0);
    dim = mean([x1;x2]);
else
    im1 = L==max_1_ind;
    [x1 y1 z1]=find(im1>0);
    dim = mean(x1);
end

image_data(round(dim)-1:round(dim)+1,:,:)=0;

maskedFileName = strtrim(fullfile(meshes_path,[name,'_masked','.nii']));
nii = make_nii(uint8(image_data),[1 1 1],[0 0 0],2, '');
save_nii(nii,maskedFileName);

%Scrub mask
scrubbedFileName = strtrim(fullfile(meshes_path,[name,'_scrubbed','.nii']));
cmd=[fullfile(app_dir,'BrainSuite21a','bin','scrubmask -i '),maskedFileName,' -o ',scrubbedFileName,' -b 1']; % Sonia used "BrainSuite15c"
system(cmd);

%Create one mask for each hemisphere
scrubbed_image=load_nii(scrubbedFileName);
scrubbed_image=scrubbed_image.img;

[L,num] = bwlabeln(scrubbed_image);
max_1=0;
max_2=0;
max_1_ind=0;
max_2_ind=0;
if (num>1)
    for j=1:num
        if sum(L(:)==j)>=max_1
            max_2=max_1;
            max_2_ind=max_1_ind;
            max_1=sum(L(:)==j);
            max_1_ind=j;
        else
            if sum(L(:)==j)>=max_2
                max_2=sum(L(:)==j);
                max_2_ind=j;
            end
        end
    end
end


if (max_1>1000 && max_2>1000)
    scrubbed_file_hemi_1 = zeros(size(scrubbed_image));
	scrubbed_file_hemi_1(L==max_1_ind) = 1;
%     scrubbed_file_hemi_1(L==max_2_ind) = 1;  %##################################      Modified by flipped image     ###################################
    s  = regionprops(scrubbed_file_hemi_1,'centroid');
    
    scrubbed_file_hemi_2 = zeros(size(scrubbed_image));
	scrubbed_file_hemi_2(L==max_2_ind) = 1;
%     scrubbed_file_hemi_2(L==max_1_ind) = 1;  %##################################      Modified by flipped image     ###################################
    s2  = regionprops(scrubbed_file_hemi_2,'centroid');
    
    if (norm(s(:,1).Centroid - [0 0 0],2)<norm(s2(:,1).Centroid - [0 0 0],2))
        
        Hemi1FileName = strtrim(fullfile(meshes_path,[name,'_hemi1','.nii']));
        nii = make_nii(uint8(scrubbed_file_hemi_1),[1 1 1],[0 0 0],2, '');
        save_nii(nii,Hemi1FileName);
        
        Hemi2FileName = strtrim(fullfile(meshes_path,[name,'_hemi2','.nii']));
        nii = make_nii(uint8(scrubbed_file_hemi_2),[1 1 1],[0 0 0],2, '');
        save_nii(nii,Hemi2FileName);
        
    else
        Hemi1FileName = strtrim(fullfile(meshes_path,[name,'_hemi1','.nii']));
        nii = make_nii(uint8(scrubbed_file_hemi_2),[1 1 1],[0 0 0],2, '');
        save_nii(nii,Hemi1FileName);
        
        Hemi2FileName = strtrim(fullfile(meshes_path,[name,'_hemi2','.nii']));
        nii = make_nii(uint8(scrubbed_file_hemi_1),[1 1 1],[0 0 0],2, '');
        save_nii(nii,Hemi2FileName);
    end
    
    
    %For hemisphere 1
    %topology correction mask
    tcaFileName = strtrim(fullfile(meshes_path,[name,'_tca','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','tca -i '),Hemi1FileName,' -o ',tcaFileName,' -m 200'];
    system(cmd);
    
    %dewisp
    dewispFileName = strtrim(fullfile(meshes_path,[name,'_dewisp','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dewisp -i '),tcaFileName,' -o ',dewispFileName];
    system(cmd);
    
    % create surface
    surfaceFileName1 = strtrim(fullfile(meshes_path,[name,'_surface_WM_1','.dfs']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dfs -i '),dewispFileName,' -o ',surfaceFileName1];
    system(cmd);
    
    %save to off format
    [NFV,hdr]=readdfs(surfaceFileName1);
    write_off(meshFileName_1,NFV.vertices', NFV.faces');
    
    cmd=[fullfile(app_dir,'meshlab.app','Contents','MacOS','meshlabserver -i '),meshFileName_1,' -o ',meshFileName_1,' -s ',scriptClean];
    system(cmd);
    
    %For hemisphere 2
    %topology correction mask
    tcaFileName = strtrim(fullfile(meshes_path,[name,'_tca','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','tca -i '),Hemi2FileName,' -o ',tcaFileName,' -m 200'];
    system(cmd);
    
    %dewisp
    dewispFileName = strtrim(fullfile(meshes_path,[name,'_dewisp','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dewisp -i '),tcaFileName,' -o ',dewispFileName];
    system(cmd);
    
    % create surface
    surfaceFileName2 = strtrim(fullfile(meshes_path,[name,'_surface_WM_2','.dfs']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dfs -i '),dewispFileName,' -o ',surfaceFileName2];
    system(cmd);
    
    %save to off format
    [NFV,hdr]=readdfs(surfaceFileName2);
    write_off(meshFileName_2,NFV.vertices', NFV.faces');
    
    cmd=[fullfile(app_dir,'meshlab.app','Contents','MacOS','meshlabserver -i '),meshFileName_2,' -o ',meshFileName_2,' -s ',scriptClean];
    system(cmd);
    
    %cmd=['rm ',maskedFileName,' ',scrubbedFileName,' ',tcaFileName,' ',dewispFileName,' ',surfaceFileName1,' ',surfaceFileName2,' ',Hemi1FileName,' ',Hemi2FileName,' ',fractionFileName,' ',brainMaskFileName];
    cmd=['rm ',maskedFileName,' ',scrubbedFileName,' ',tcaFileName,' ',dewispFileName,' ',surfaceFileName1,' ',surfaceFileName2,' ',Hemi1FileName,' ',Hemi2FileName];
    system(cmd);
    
else
    
    %topology correction mask
    tcaFileName = strtrim(fullfile(meshes_path,[name,'_tca','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','tca -i '),scrubbedFileName,' -o ',tcaFileName,' -m 200'];
    system(cmd);
    
    %dewisp
    dewispFileName = strtrim(fullfile(meshes_path,[name,'_dewisp','.nii']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dewisp -i '),tcaFileName,' -o ',dewispFileName];
    system(cmd);
    
    % create surface
    surfaceFileName = strtrim(fullfile(meshes_path,[name,'_surface_WM','.dfs']));
    cmd=[fullfile(app_dir,'BrainSuite21a','bin','dfs -i '),dewispFileName,' -o ',surfaceFileName];
    system(cmd);
    
    %save to off format
    [NFV,hdr]=readdfs(surfaceFileName);
    meshFileName_2 = strtrim(fullfile(meshes_path,[name,'_surface_WM','.off']));
    write_off(meshFileName_2,NFV.vertices', NFV.faces');
    
    %cmd=['rm ',scrubbedFileName,' ',maskedFileName,' ',tcaFileName,' ',dewispFileName,' ',surfaceFileName,' ',brainMaskFileName,' ',fractionFileName];
    cmd=['rm ',scrubbedFileName,' ',maskedFileName,' ',tcaFileName,' ',dewispFileName,' ',surfaceFileName];
    system(cmd);
end

clear im1 image_data image_data_nii L j kk NFV nii s s2 scrubbed_file_hemi_1 scrubbed_file_hemi_2 scrubbed_image x1 y1 z1;
%% Label vertices
% Create volumetric mesh
image_labels2=image_labels;
image_labels2(image_labels==12)=0;
image_labels2(image_labels==5)=0;

image_insula_1=(image_labels==5);
se = strel('cube',2);
image_insula_1 = imerode(image_insula_1, se);
CC=bwconncomp(image_insula_1);
numPixels=cellfun(@numel,CC.PixelIdxList);
[biggest,idx]=max(numPixels);
image_labels2(CC.PixelIdxList{idx})=5;


image_insula_2=(image_labels==12);
se = strel('cube',2);
image_insula_2 = imerode(image_insula_2, se);
CC=bwconncomp(image_insula_2);
numPixels=cellfun(@numel,CC.PixelIdxList);
[biggest,idx]=max(numPixels);
image_labels2(CC.PixelIdxList{idx})=12;


ISO2MESH_TEMP=meshes_path;

opt.radbound=1.5;
[node_raw,face,~]=v2m(uint8(image_labels2),unique(image_labels2),opt,100,'cgalmesh');  % (Yuan) This step will show  "surface and volume meshes complete" % Elapsed time is 19.517804 seconds.


%For each hemisphere determine nodes labels
%Hem 1
[vertices,faces]=read_off(meshFileName_1);

node=node_raw(node_raw(:,1)<ceil(max(vertices(1,:))) & node_raw(:,1)>floor(min(vertices(1,:))),:); % Added by Yuan on 9/28/2018

%  dists=pdist2(node(:,1:3),vertices');
dists=pdist2(vertices',node(:,1:3)); % Elapsed time is ~4 seconds.
[~,ids]=min(dists'); % Elapsed time is 8.363290 seconds.
labels=node(ids,4);

% clf
% hold on
% scatter3(vertices(1,:),vertices(2,:),vertices(3,:),'y')
% scatter3(node(:,1),node(:,2),node(:,3),'r')
% axis equal
% hold off

ids=find(((labels==0)|(labels==-1)));
vertices_red = vertices';
vertices_red(ids,:)=[];
labels_red=labels;
labels_red(ids)=[];

dists_ref=pdist2(vertices_red,vertices(:,ids)');
[~,ids2]=min(dists_ref);
labels(ids)=labels_red(ids2);

%Correct errors due to pb in registration for insula by changing the
%islands label to the temporal pole label
labelInsula = 5; % 12; %##################################      Modified by flipped image     ###################################
AdjMatrix = triangulation2adjacency_Sparse(faces);

[labeledComponents,sizeOfLabeledComponents,~,~]=getConnectedComponentLabeledMesh(vertices',AdjMatrix, labels,labelInsula);

[~,sortedLab]=sort(cell2mat(sizeOfLabeledComponents),'descend');

for kk=2:max(size(sortedLab))
    labels(labeledComponents{sortedLab(kk)})=11;
end

%Hem 2
[vertices2,faces2]=read_off(meshFileName_2);

node=node_raw(node_raw(:,1)<ceil(max(vertices2(1,:))) & node_raw(:,1)>floor(min(vertices2(1,:))),:); % Added by Yuan on 9/28/2018

dists2=pdist2(vertices2',node(:,1:3));
[~,ids2]=min(dists2'); % (Yuan) This step takes time.
labels2=node(ids2,4);

ids3=find(((labels2==0)|(labels2==-1)));
vertices_red = vertices2';
vertices_red(ids3,:)=[];
labels_red=labels2;
labels_red(ids3)=[];

dists_ref=pdist2(vertices_red,vertices2(:,ids3)');
[~,ids2]=min(dists_ref);
labels2(ids3)=labels_red(ids2);

%Correct errors due to pb in registration for insula by changing the
%islands label to the temporal pole label
labelInsula = 12; %5;   %##################################      Modified by flipped image     ###################################
AdjMatrix = triangulation2adjacency_Sparse(faces2);
[labeledComponents,sizeOfLabeledComponents,~,~]=getConnectedComponentLabeledMesh(vertices2',AdjMatrix, labels2,labelInsula);

[~,sortedLab]=sort(cell2mat(sizeOfLabeledComponents),'descend');

for kk=2:max(size(sortedLab))
    labels2(labeledComponents{sortedLab(kk)})=4;
end

%Remove the elements in common in both hemisphere due to mismatch in
%registration
sums_1=[];
sums_2=[];
labs_uniques=unique(image_labels);
for k=1:max(size(labs_uniques))
    sums_1=[sums_1; sum(labels==labs_uniques(k))];
    sums_2=[sums_2; sum(labels2==labs_uniques(k))];
end

for k=1:max(size(labs_uniques))
    if(sums_1(k)>0 && sums_2(k)>0 && labs_uniques(k)~=15 && labs_uniques(k)~=16)%%Don't take into account subcortical as common to both hemispheres
        if(sums_1(k)>sums_2(k))
            labels2(labels2==labs_uniques(k))=0;
            
        else
            labels(labels==labs_uniques(k))=0;
        end
        
    end
end

ids=find((labels==0));
vertices_red = vertices';
vertices_red(ids,:)=[];
labels_red=labels;
labels_red(ids)=[];

dists_ref=pdist2(vertices_red,vertices(:,ids)');
[~,ids2]=min(dists_ref);
labels(ids)=labels_red(ids2);

ids=find((labels2==0));
vertices_red = vertices2';
vertices_red(ids,:)=[];
labels_red=labels2;
labels_red(ids)=[];
dists_ref=pdist2(vertices_red,vertices2(:,ids)');
[~,ids2]=min(dists_ref);
labels2(ids)=labels_red(ids2);

[labels_path,name,~] = fileparts(outputLabeFileName);
outputLabeFileName_1=strtrim(fullfile(labels_path,[name,'_1.csv'])); % Right hemisphere
outputLabeFileName_2=strtrim(fullfile(labels_path,[name,'_2.csv'])); % Left hemisphere

csvwrite(outputLabeFileName_1,labels);
csvwrite(outputLabeFileName_2,labels2);
clear node vertices2 faces2 AdjMatrix labels labels2 dists dists2 labels_red labeledComponents sizeOfLabeledComponents dists_ref opt ids3;
clear AdjMatrix vertices faces labeledComponents sizeOfLabeledComponents dists vertices_red labels_red ids ids2 sortedLab image_labels hdr;
end