added dof2: adjusted dofs to GLM cleanup

meldata.cc

@@ -620,7 +620,7 @@ namespace Melodic{
     message("Creating new PEs, DOF and sigmasquareds ... " << endl);
     // First, save all the old files
     volume4D<float> tmpvol, old_dof;
-    Matrix pemaps;
+    Matrix pemaps, old_ssq;
     volumeinfo tmpvolinfo;
     int pe_ctr = 1;
     while(fsl_imageexists(opts.glmcleanup.value() + "/stats/pe" + num2str(pe_ctr)))
@@ -631,8 +631,11 @@ namespace Melodic{
 	else pemaps = pemaps & tmpvol.matrix(Mask);
 	pe_ctr++;
       }
+
     read_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/sigmasquareds",tmpvolinfo);
     save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/old_sigmasquareds",tmpvolinfo);
+    old_ssq = tmpvol.matrix(Mask);
+
     if(fsl_imageexists(opts.glmcleanup.value() + "/stats/dof"))
       {
 	read_volume4D(old_dof,opts.glmcleanup.value() + "/stats/dof",tmpvolinfo);
@@ -696,10 +699,9 @@ namespace Melodic{
     
     allmaps = pemaps & threshmaps;
 
-    outMsize("    allmaps ",allmaps);
+    outMsize(" size of pes and thresholded maps: ",allmaps);
     //outMsize("    orig_data ",orig_data);
 
-    allmaps = pemaps & threshmaps;
     alltcs = ( orig_data * allmaps.t() ) * pinv (allmaps * allmaps.t() ); 
     
     //outMsize("    alltcs",alltcs);
@@ -713,8 +715,11 @@ namespace Melodic{
     sigsq = sum(SP(resi,resi)) / R.Trace(); 
     //outMsize("  ... new sigmasquareds :",sigsq);
     
-    Matrix new_dof;
-    new_dof = old_dof.matrix(Mask);
+    Matrix new_dof, new_dof2;
+    new_dof  = old_dof.matrix(Mask);
+    new_dof2 = old_dof.matrix(Mask);
+
+    // Calculate new_dof
     for(int j_ctr=1;j_ctr<=threshmaps.Ncols();j_ctr++)
       for(int i_ctr=1;i_ctr<=threshmaps.Nrows();i_ctr++)
 	if(abs(threshmaps(i_ctr,j_ctr)) > 0)
@@ -722,6 +727,12 @@ namespace Melodic{
     //outMsize("  ... new DOFs :", new_dof);
     message(" ... done " << endl); 
     
+
+    // Calculate new_dof2
+    //    for(int j_ctr=1;j_ctr<=threshmaps.Ncols();j_ctr++)
+    //  for(int i_ctr=1;i_ctr<=threshmaps.Nrows();i_ctr++)
+    new_dof2 = SP(new_dof2,SP(sigsq,pow(old_ssq,-1)));
+
     //save out new GLM results
     message(" Saving results ... ");
     tmpvol.setmatrix(sigsq,Mask);
@@ -729,13 +740,16 @@ namespace Melodic{
     
     for(int tmp_ctr=1;tmp_ctr<=pemaps.Nrows(); tmp_ctr++){
       tmpvol.setmatrix(pemaps.Row(tmp_ctr),Mask);
-      save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/pe" + num2str(tmp_ctr));
+      save_volume4D(tmpvol,opts.glmcleanup.value()+"/stats/pe"+num2str(tmp_ctr));
     }
     tmpvol.setmatrix(new_dof,Mask);
     save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/dof");
+   
+    tmpvol.setmatrix(new_dof2,Mask);
+    save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/dof2");
+   
     message(" ... done " << endl);
     message(endl << " finished GLM cleanup ! " << endl << endl);
-   
   }
 
   Matrix MelodicData::calc_FFT(const Matrix& Mat)