in progress

film_gls.cc

 /*  film_gls.cc
 
-    Mark Woolrich, FMRIB Image Analysis Group
+    Mark Woolrich and Matthew Webster, FMRIB Image Analysis Group
 
-    Copyright (C) 1999-2000 University of Oxford  */
+    Copyright (C) 1999-2008 University of Oxford  */
 
 /*  CCOPYRIGHT  */
 
 #include <iostream>
-#include <fstream>
-#include <sstream>
 #define WANT_STREAM
 #define WANT_MATH
 
 #include "newimage/newimageall.h"
-#include "newmatap.h"
-#include "newmatio.h"
-#include "miscmaths/volumeseries.h"
-#include "miscmaths/volume.h"
-#include "miscmaths/miscmaths.h"
 #include "utils/log.h"
 #include "AutoCorrEstimator.h"
 #include "paradigm.h"
@@ -50,31 +43,25 @@ int main(int argc, char *argv[])
     int sizeTS = input_data.tsize();
 
     Matrix datam;
-    VolumeSeries x,y;
-    y.read(globalopts.inputfname);
-    x.setInfo(y.getInfo());
-    x=input_data.matrix();
-    x.thresholdSeries(globalopts.thresh, true);
+    volume<float> mask(input_data[0]);
 
-    //save_volume_dtype(input_data[int(sizeTS/2)-1],logger.getDir() + "/" + globalopts.epifname,DT_SIGNED_INT );
-    save_volume(input_data[int(sizeTS/2)-1],logger.getDir() + "/" + globalopts.epifname);
+    volume4D<float> reference(mask.xsize(),mask.ysize(),mask.zsize(),1);
+    reference[0]=input_data[int(sizeTS/2)-1];
+    copybasicproperties(input_data,reference);
+    FslSetCalMinMax(&vinfo,reference.min(),reference.max());
+    save_volume_dtype(reference[0],logger.getDir() + "/" + globalopts.epifname,DT_SIGNED_SHORT,vinfo);
 
-    volume<float> mask(input_data[0]);
-    volume4D<float> epivol2(mask.xsize(),mask.ysize(),mask.zsize(),1);
-    epivol2[0]=input_data[int(sizeTS/2)-1];
     mask=0;
     for(int t=0;t<input_data.tsize();t++) mask+=input_data[t];
     mask/=input_data.tsize();
-    for(int z=0;z<input_data.zsize();z++)
-      for(int y=0;y<input_data.ysize();y++)          
-	for(int x=0;x<input_data.xsize();x++)
-	  if ( mask(x,y,z) > globalopts.thresh )
-	    for(int t=0;t<input_data.tsize();t++) input_data(x,y,z,t)-= mask(x,y,z);
+    input_data-=mask;
     mask.binarise(globalopts.thresh,mask.max()+1,exclusive);
     input_data*=mask;
     datam=input_data.matrix(mask);
    
-    int numTS = mask.sum();
+    int numTS = datam.Ncols();
+    ColumnVector epivol = reference.matrix(mask).t();
+    //epivol =  datam.Row(int(sizeTS/2)).AsColumn();
     
     // Load paradigm:
     Paradigm parad;
@@ -106,7 +93,7 @@ int main(int argc, char *argv[])
     GlimGls glimGls(numTS, sizeTS, numParams);
 
     // Residuals container:
-    VolumeSeries residuals(sizeTS, numTS, x.getInfo(), x.getPreThresholdPositions());
+    Matrix residuals(sizeTS, numTS);
     Matrix residualsm = datam;
 
     // Setup autocorrelation estimator:
@@ -120,7 +107,7 @@ int main(int argc, char *argv[])
 	for(int i = 1; i <= numTS; i++)
 	  {						    
             glimGls.setData(datam.Column(i), parad.getDesignMatrix(), i);
-	    residuals.setSeries(glimGls.getResiduals(),i);
+	    residuals.Column(i)=glimGls.getResiduals();
 	  }
 	cout << "Completed" << endl; 
 	
@@ -144,7 +131,6 @@ int main(int argc, char *argv[])
 	    // SUSAN smooth raw estimates:
 	    if(globalopts.smoothACEst)
 	      {
-		ColumnVector epivol = epivol2.matrix(mask).t();
 		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.epith, globalopts.tukeysize);	
 	      }	    
 		
@@ -162,7 +148,6 @@ int main(int argc, char *argv[])
 	    // Smooth raw estimates:
 	    if(globalopts.smoothACEst)
 	      { 
-		ColumnVector epivol = epivol2.matrix(mask).t();
 		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.epith);
 	      }
 	    
@@ -218,7 +203,7 @@ int main(int argc, char *argv[])
 		co++;
 	      }
 	    glimGls.setData(datam.Column(i), parad.getDesignMatrix(), i);
-	    residuals.setSeries(glimGls.getResiduals(),i);
+	    residuals.Column(i)=glimGls.getResiduals();
 
 	    if(globalopts.output_pwdata || globalopts.verbose)
 	      {
@@ -236,9 +221,9 @@ int main(int argc, char *argv[])
 
     cerr << "Saving results... " << endl;
 
-    // write residuals
-    residuals.writeThresholdedSeriesAsFloat(x.getInfo(),x.getPreThresholdPositions(),logger.getDir() + "/res4d");
-    residuals.CleanUp();
+    input_data.setmatrix(residuals,mask);
+    FslSetCalMinMax(&vinfo,input_data.min(),input_data.max());
+    save_volume4D(input_data,logger.getDir() + "/res4d",vinfo);
 
     if(globalopts.output_pwdata || globalopts.verbose)
       {
@@ -251,25 +236,23 @@ int main(int argc, char *argv[])
 		
       }
 
-    x.CleanUp();
-
     // Write out threshac:
     Matrix& threshacm = acEst.getEstimates();
     int cutoff = sizeTS/2;
-    if(globalopts.tukey)
-      cutoff = globalopts.tukeysize;
-    cutoff = MISCMATHS::Max(1,cutoff);
-    threshacm = threshacm.Rows(1,cutoff); 
-    volume4D<float> threshac;
-    threshac.setmatrix(threshacm,mask);
-    copybasicproperties(input_data,threshac);
-    threshac.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
-    FslSetCalMinMax(&vinfo,threshac.min(),threshac.max());
-    save_volume4D(threshac,logger.getDir() + "/threshac1",vinfo);
+    if(globalopts.tukey) cutoff = globalopts.tukeysize;
+    //cutoff = MISCMATHS::Max(1,cutoff);
+    threshacm = threshacm.Rows(1,MISCMATHS::Max(1,cutoff)); 
+
+    input_data.setmatrix(threshacm,mask);
+    input_data.settdim(reference.tdim());
+    input_data.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
+    FslSetCalMinMax(&vinfo,input_data.min(),input_data.max());
+    save_volume4D(input_data,logger.getDir() + "/threshac1",vinfo);
+
     threshacm.CleanUp();
 
     // save gls results:
-    glimGls.Save(x.getInfo(), x.getPreThresholdPositions());
+    glimGls.Save(vinfo,mask,reference.tdim());
     glimGls.CleanUp();
 
     cerr << "Completed" << endl;