removing marks volume: In progress

film_gls.cc

@@ -12,6 +12,7 @@
 #define WANT_STREAM
 #define WANT_MATH
 
+#include "newimage/newimageall.h"
 #include "newmatap.h"
 #include "newmatio.h"
 #include "miscmaths/volumeseries.h"
@@ -29,6 +30,7 @@ using namespace NEWMAT;
 using namespace FILM;
 using namespace Utilities;
 using namespace MISCMATHS;
+using namespace NEWIMAGE;
 
 int main(int argc, char *argv[])
 {
@@ -42,10 +44,32 @@ int main(int argc, char *argv[])
     globalopts.parse_command_line(argc, argv, logger);
 
     // load data
-    VolumeSeries x;
-    x.read(globalopts.inputfname);
-
-    int sizeTS = x.getNumVolumes();
+    volume4D<double> input_data;
+    read_volume4D(input_data,globalopts.inputfname);
+
+    Matrix datam;
+    VolumeSeries x,y;
+    y.read(globalopts.inputfname);
+    x.setInfo(y.getInfo());
+    x=input_data.matrix();
+
+    volume<double> mask(input_data[0]);
+    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);
+    mask.binarise(globalopts.thresh,mask.max()+1,exclusive);
+    input_data*=mask;
+    save_volume4D(input_data,"foo2");
+    datam=input_data.matrix(mask);
+
+    int sizeTS = input_data.tsize();
+    int numTS = mask.sum();
+    cerr << datam.Nrows() << " " << datam.Ncols() << endl;
 
     // if needed for SUSAN spatial smoothing, output the halfway volume for use later
     Volume epivol;
@@ -58,6 +82,7 @@ int main(int argc, char *argv[])
 
     // This also removes the mean from each of the time series:
     x.thresholdSeries(globalopts.thresh, true);
+    x.writeThresholdedSeriesAsFloat(x.getInfo(),x.getPreThresholdPositions(),"foo3");
 
     // if needed for SUSAN spatial smoothing, also threshold the epi volume
     if(globalopts.smoothACEst)
@@ -66,8 +91,6 @@ int main(int argc, char *argv[])
 	epivol.threshold();
       }
     
-    int numTS = x.getNumSeries();
-
     // Load paradigm:
     Paradigm parad;
     if(!globalopts.ac_only)
@@ -108,7 +131,8 @@ int main(int argc, char *argv[])
 	cout << "Calculating residuals..." << endl; 
 	for(int i = 1; i <= numTS; i++)
 	  {						    
-	    glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
+	    //glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
+            glimGls.setData(datam.Column(i), parad.getDesignMatrix(), i);
 	    residuals.setSeries(glimGls.getResiduals(),i);
 	  }
 	cout << "Completed" << endl;