function reconstructAndLabelConvexHull(reduced_parcellationPath,meshFileName_1,meshFileName_2,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)
    
    meshFileName_1=strtrim(fullfile(meshes_path,[name,'_WM_1.off']));
    meshFileName_2=strtrim(fullfile(meshes_path,[name,'_WM_2.off']));               
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;
    
    s  = regionprops(scrubbed_file_hemi_1,'centroid');
    scrubbed_file_hemi_2 = zeros(size(scrubbed_image));
    scrubbed_file_hemi_2(L==max_2_ind) = 1;
    
    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;
