function[F_stats,detvtx,sided_pval]=runlme3(hemi,sel,contrast,out_dir,qdecname)
sel=6;



dir1=sprintf('~/ANALYSIS/lme/%s',out_dir);
%
%http://www.mathworks.com/help/techdoc/ref/exist.html
if exist(dir1,'dir')
    display('folder exist')
else
    mkdir(dir1)
end    

        
        



[Y,mri]=fs_read_Y(sprintf('%s.long_thickness.30.mgh',hemi));

% The fsaverage's spherical surface (lh.sphere) and cortex label (lh.cortex.label) were read with:
sphere=fs_read_surf(sprintf('$FREESURFER_HOME/subjects/fsaverage/surf/%s.sphere',hemi));
cortex=fs_read_label(sprintf('$FREESURFER_HOME/subjects/fsaverage/label/%s.cortex.label',hemi));


% read Qdec dat file
Qdec=fReadQdec(qdecname);
Qdec = rmQdecCol(Qdec,1);
sID = Qdec(2:end,1);

Qdec = rmQdecCol(Qdec,1);



label=Qdec(1,:);





X = [ones(length(M),1), M(:,1),M(:,1).^2,M(:,2),M(:,2).*M(:,1),M(:,3) ];
% interception, time, time^2, group, group*time, gender
%contrast
field = 'C';
value=[0 0 0 0 1 0];
% interception, time, time^2, group, group*time, gender
CM=struct(field,value);











[M,Y,ni]=sortData(M,1,Y,sID);



%reads in mgh file
%;
%Thus, our design matrix X was built from matrix M to represent the previous model.Then initial vertex-wise temporal covariance estimates were computed with:

[Th0,Re] = lme_mass_fit_EMinit(X,[1 2],Y,ni,cortex,3);

%These covariance estimates were segmented into homogeneous regions using:

[Rgs,RgMeans] = lme_mass_RgGrow(sphere,Re,Th0,cortex,2,95);


%Here both lhTh0 and lhRgMeans maps were overlaid onto lhsphere and visually compared each other to ensure that they were similar enough (the essential spatial organization of the initial covariance estimates was not lost after the segmentation procedure, otherwise the above input value 2 must be reduced to 1.8 and so on):

surf.faces = sphere.tri;

surf.vertices = sphere.coord';

figure; p1 = patch(surf);

set(p1,'facecolor','interp','edgecolor','none','facevertexcdata',Th0(1,:)');

figure; p2 = patch(surf); set(p2,'facecolor','interp','edgecolor','none','facevertexcdata',RgMeans(1,:)');

%The spatiotemporal LME model was fitted using:

stats = lme_mass_fit_Rgw(X,[1 2],Y,ni,Th0,Rgs,sphere);










%i) Fit the model with one random effect using the segmentation obtained from the previous model:

Th0_1RF = lme_mass_fit_EMinit(X,[1],Y,ni,cortex,3);

stats_1RF = lme_mass_fit_Rgw(X,[1],Y,ni,Th0_1RF,Rgs,sphere);

%ii) Apply the likelihood ratio test:

LR_pval = lme_mass_LR(stats,stats_1RF,1);

%iii) Correct for multiple comparisons:

dvtx = lme_mass_FDR2(LR_pval,ones(1,length(LR_pval)),cortex,0.05,0);

compare=100*length(dvtx)/length(cortex);
disp(compare);

fprintf('Lenght of fdr compare to total length %f',compare);

txt=contrasttxt(CM.C);

       

F_stats = lme_mass_F(stats,CM); 

[detvtx,sided_pval,pth]  = lme_mass_FDR2(F_stats.pval,F_stats.sgn,cortex,0.05,0); 
fs_write_fstats(F_stats,mri,sprintf('%s/%s_%s_%s_sig.mgh',dir1,'all',hemi,txt),'sig'); 
fs_write_fstats(F_stats,mri,sprintf('%s/%s_%s_%s_F.mgh',dir1,'all',hemi,txt),'fval'); 
fs_write_fstats(F_stats,mri,sprintf('%s/%s_%s_%s_p.mgh',dir1,'all',hemi,txt),'pval');