function [Mean_Cmean,Mean_Cgauss,Mean_Curvedness,Mean_lGI,Mean_SulcalDepth,Volume_all,Area_mean_all]=...
    feature_calc(parcellation_path,Sbj_ID,Scan_ID,WM_N)

    if WM_N==1 || WM_N==2
        labels=csvread(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_labels_' num2str(WM_N) '.csv']));

        load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(WM_N) '_foldingMeasures.mat']));
        load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(WM_N) '_lGI.mat']));
        load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(WM_N) '_area_mean.mat']));
        load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(WM_N) '_SulcalDepth.mat']));

    elseif WM_N==3
        labels1=csvread(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_labels_' num2str(1) '.csv']));
        labels2=csvread(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_labels_' num2str(2) '.csv']));
        labels=[labels1;labels2];

        [FM1]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(1) '_foldingMeasures.mat']));
        [FM2]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(2) '_foldingMeasures.mat']));
        foldingMeasures.Cmean=[FM1.foldingMeasures.Cmean;FM2.foldingMeasures.Cmean];
        foldingMeasures.Cgauss=[FM1.foldingMeasures.Cgauss;FM2.foldingMeasures.Cgauss];
        foldingMeasures.Curvedness=[FM1.foldingMeasures.Curvedness;FM2.foldingMeasures.Curvedness]; 
        
        [lGI1]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(1) '_lGI.mat']));
        [lGI2]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(2) '_lGI.mat']));
        lGI=[lGI1.lGI;lGI2.lGI];
        
        [area_mean1]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(1) '_area_mean.mat']));
        [area_mean2]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(2) '_area_mean.mat']));
        area_mean=[area_mean1.area_mean;area_mean2.area_mean];

        [SulcalDepth1]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(1) '_SulcalDepth.mat']));
        [SulcalDepth2]=load(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_' num2str(2) '_SulcalDepth.mat']));
        Sulcal_Depth=[SulcalDepth1.Sulcal_Depth;SulcalDepth2.Sulcal_Depth]; 
    end

    %########## Calculate Volume from parcellation #########
    WM_SubCort_nii = load_untouch_nii(fullfile(parcellation_path,[num2str(Sbj_ID) '-0' num2str(Scan_ID) '_WM_SubCort_labels.nii']));
    pixdimX_par=WM_SubCort_nii.hdr.dime.pixdim(2);
    pixdimY_par=WM_SubCort_nii.hdr.dime.pixdim(3);
    pixdimZ_par=WM_SubCort_nii.hdr.dime.pixdim(4);
    WM_SubCort = uint32(WM_SubCort_nii.img);

    xx=reshape(WM_SubCort,[1,size(WM_SubCort,1)*size(WM_SubCort,2)*size(WM_SubCort,3)])';
    aa=unique(xx);
    WM_SubCort_unique=[aa,histc(xx(:),aa)];
    %#######################################################

    label_ID1=[8 9 11  3  5  6  7 15 17 19];
    label_ID2=[1 2  4 10 12 13 14 16 18 20];
    if WM_N==1 % Right hemisphere
    	label_ID=label_ID1;
    elseif WM_N==2 % Left hemisphere
    	label_ID=label_ID2;
    end
    
    if WM_N==1 || WM_N==2
        for label_ID_count=1:length(label_ID)
            label_idx=labels==label_ID(label_ID_count);
            Mean_Cmean(label_ID_count)=sum(foldingMeasures.Cmean(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_Cgauss(label_ID_count)=sum(foldingMeasures.Cgauss(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_Curvedness(label_ID_count)=sum(foldingMeasures.Curvedness(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_lGI(label_ID_count)=sum(lGI(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_SulcalDepth(label_ID_count)=sum(Sulcal_Depth(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));

            %########## Calculate Volume from parcellation #########
            if sum(WM_SubCort_unique(:,1)==label_ID(label_ID_count))==0
                Volume_all(label_ID_count)=0;
            else
                Volume_all(label_ID_count)=double(WM_SubCort_unique(WM_SubCort_unique(:,1)==label_ID(label_ID_count),2))*pixdimX_par*pixdimY_par*pixdimZ_par;
            end
            %########## Calculate Area from area_mean #########
            Area_mean_all(label_ID_count)=sum(area_mean(label_idx));
            %#######################################################
        end
        label_ID_count=label_ID_count+1;

        labels_WM=labels==label_ID(1) | labels==label_ID(2) | labels==label_ID(3) | labels==label_ID(4);
        Mean_Cmean(label_ID_count)=sum(foldingMeasures.Cmean(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_Cgauss(label_ID_count)=sum(foldingMeasures.Cgauss(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_Curvedness(label_ID_count)=sum(foldingMeasures.Curvedness(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_lGI(label_ID_count)=sum(lGI(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_SulcalDepth(label_ID_count)=sum(Sulcal_Depth(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        
        %########## Calculate Volume from parcellation #########
        Volume_all(label_ID_count)=sum(Volume_all([1 2 3 4]));
        %########## Calculate Area from area_mean #########
        Area_mean_all(label_ID_count)=sum(Area_mean_all([1 2 3 4]));
        %#######################################################

    elseif WM_N==3
        for label_ID_count=1:length(label_ID1)
            label_idx=(labels==label_ID1(label_ID_count) | labels==label_ID2(label_ID_count));
            Mean_Cmean(label_ID_count)=sum(foldingMeasures.Cmean(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_Cgauss(label_ID_count)=sum(foldingMeasures.Cgauss(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_Curvedness(label_ID_count)=sum(foldingMeasures.Curvedness(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_lGI(label_ID_count)=sum(lGI(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            Mean_SulcalDepth(label_ID_count)=sum(Sulcal_Depth(label_idx).*area_mean(label_idx))/sum(area_mean(label_idx));
            
            %########## Calculate Volume from parcellation #########
            if sum(WM_SubCort_unique(:,1)==label_ID1(label_ID_count) | WM_SubCort_unique(:,1)==label_ID2(label_ID_count))==0
                Volume_all(label_ID_count)=0;
            else
                Volume_all(label_ID_count)=sum(double(WM_SubCort_unique(WM_SubCort_unique(:,1)==label_ID1(label_ID_count) | WM_SubCort_unique(:,1)==label_ID2(label_ID_count),2)))*pixdimX_par*pixdimY_par*pixdimZ_par;
            end
            %########## Calculate Area from area_mean #########
            Area_mean_all(label_ID_count)=sum(area_mean(label_idx));
            %#######################################################
        end
        label_ID_count=label_ID_count+1;

        labels_WM=labels==label_ID1(1) | labels==label_ID1(2) | labels==label_ID1(3) | labels==label_ID1(4) |...
            labels==label_ID2(1) | labels==label_ID2(2) | labels==label_ID2(3) | labels==label_ID2(4);
        Mean_Cmean(label_ID_count)=sum(foldingMeasures.Cmean(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_Cgauss(label_ID_count)=sum(foldingMeasures.Cgauss(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_Curvedness(label_ID_count)=sum(foldingMeasures.Curvedness(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_lGI(label_ID_count)=sum(lGI(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        Mean_SulcalDepth(label_ID_count)=sum(Sulcal_Depth(labels_WM).*area_mean(labels_WM))/sum(area_mean(labels_WM));
        
        %########## Calculate Volume from parcellation #########
        Volume_all(label_ID_count)=sum(Volume_all([1 2 3 4]));
        %########## Calculate Area from area_mean #########
        Area_mean_all(label_ID_count)=sum(Area_mean_all([1 2 3 4]));
        %#######################################################

    end

end