From 9ced373966b75bd99a2537ff9d36ad944ae7c619 Mon Sep 17 00:00:00 2001
From: Stephen Smith <steve@fmrib.ox.ac.uk>
Date: Wed, 8 Feb 2006 20:10:01 +0000
Subject: [PATCH]  CVS:
 ----------------------------------------------------------------------

---
 volumeseries.cc | 135 ++++++++++++++++++++++++++----------------------
 volumeseries.h  |   1 +
 2 files changed, 75 insertions(+), 61 deletions(-)

diff --git a/volumeseries.cc b/volumeseries.cc
index 1638665..9bb49a3 100644
--- a/volumeseries.cc
+++ b/volumeseries.cc
@@ -154,76 +154,54 @@ namespace MISCMATHS {
       volinfo.miscinfo = FslInit();
       FslCloneHeader(volinfo.miscinfo,IP);
 
-      size_t imagesize=x*y*z*v;
+      size_t imagesize=x*y*z;
       FslGetDataType(IP,&type);
   
       output.ReSize(v,x*y*z);
   
-      switch(type)
+      for(size_t i=1; i<=(size_t)v; i++)
 	{
-	case DT_SIGNED_SHORT:
-	  {
-	    short* sbuffer=new short[imagesize];
-	
-	    FslReadVolumes(IP, sbuffer,v);
-	
-	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
-	    size_t volstart = 1;
-
-	    for(size_t i=1; i<=(size_t)v; i++)
+	  switch(type)
+	    {
+	    case DT_SIGNED_SHORT:
 	      {
-		volstart = (i-1)*volsize;
-		for(size_t j = 1; j<=(size_t)x*y*z; j++)
+		short* sbuffer=new short[imagesize];
+		FslReadVolumes(IP, sbuffer,1);
+		for(size_t j = 1; j<=imagesize; j++)
 		  {
-		    if (doscaling==0) { output(i,j)=sbuffer[volstart+j-1]; }
-		    else { output(i,j)=(slope * sbuffer[volstart+j-1]) + intercept;}
+		    if (doscaling==0) { output(i,j)=sbuffer[j-1]; }
+		    else { output(i,j)=(slope * sbuffer[j-1]) + intercept;}
 		  }
+		delete[] sbuffer;
 	      }
-	  
-	    delete[] sbuffer;
-	  }
-	  break;
-	case DT_FLOAT:
-	  {
-	    float* fbuffer=new float[imagesize];
-	    FslReadVolumes(IP,fbuffer,v);
-
-	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
-	    size_t volstart = 1;
-	    for(size_t i=1; i<=(size_t)v; i++)
+	      break;
+	    case DT_FLOAT:
 	      {
-		volstart = (i-1)*volsize;
-		for(size_t j = 1; j<=(size_t)x*y*z; j++)
+		float* fbuffer=new float[imagesize];
+		FslReadVolumes(IP,fbuffer,1);
+		for(size_t j = 1; j<=imagesize; j++)
 		  {
-		    if (doscaling==0) { output(i,j)=fbuffer[volstart+j-1]; }
-		    else { output(i,j)=(slope * fbuffer[volstart+j-1]) + intercept;}
-
+		    if (doscaling==0) { output(i,j)=fbuffer[j-1]; }
+		    else { output(i,j)=(slope * fbuffer[j-1]) + intercept;}
 		  }
+		delete[] fbuffer;
 	      }
-	    delete[] fbuffer;
-	  }
-	  break;
-	case DT_UNSIGNED_CHAR:
-	  {
-	    unsigned char* cbuffer=new unsigned char[imagesize];
-	    FslReadVolumes(IP,cbuffer,v);
-
-	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
-	    size_t volstart = 1;
-	    for(size_t i=1; i<=(size_t)v; i++)
+	      break;
+	    case DT_UNSIGNED_CHAR:
 	      {
-		volstart = (i-1)*volsize;
-		for(size_t j = 1; j<=(size_t)x*y*z; j++)
+		unsigned char* cbuffer=new unsigned char[imagesize];
+		FslReadVolumes(IP,cbuffer,1);
+		for(size_t j = 1; j<=imagesize; j++)
 		  {
-		    if (doscaling==0) { output(i,j)=cbuffer[volstart+j-1]; }
-		    else { output(i,j)=(slope * cbuffer[volstart+j-1]) + intercept;}
+		    if (doscaling==0) { output(i,j)=cbuffer[j-1]; }
+		    else { output(i,j)=(slope * cbuffer[j-1]) + intercept;}
 		  }
+		delete[] cbuffer;
 	      }
-	    delete[] cbuffer;
-	  }
-	  break;
-	default:
-	  perror("FslRead: DT not supported");
+	      break;
+	    default:
+	      perror("FslRead: DT not supported");
+	    }
 	}
 
       FslClose(IP);
@@ -252,23 +230,58 @@ namespace MISCMATHS {
        FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2, 
 		    volinfo.intent_p3);
 
-       int volStart = 1;
        int volSize = getNumSeries();
        int volNum = getNumVolumes();
 
-       float *qv = new float[volSize*volNum];
+       FslWriteHeader(OP);
+
+       float *qv = new float[volSize];
      
        for(int i = 1; i<= volNum; i++)
 	 {
-	   volStart = (i-1)*volSize;
 	    for(int j = 1; j <= volSize; j++)
-	     {
-	        qv[volStart+j-1] = (*this)(i,j);
-	     }
+	      qv[j-1] = (*this)(i,j);
+	    FslWriteVolumes(OP, qv, 1);
 	 }
-            
+
+       delete [] qv;
+
+       FslClose(OP);
+       
+    }
+
+  void VolumeSeries::writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const string& fname)
+    {
+       volinfo = pvolinfo;
+       preThresholdPositions = in;
+
+       Time_Tracer ts(string("VolumeSeries::writeThresholdedSeriesAsFloat" + fname).c_str());
+
+       FSLIO* OP = FslOpen(fname.c_str(), "wb");
+
+       FslCloneHeader(OP,volinfo.miscinfo);
+       
+       FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
+       FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
+       FslSetDataType(OP, DT_FLOAT);
+       FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2, volinfo.intent_p3);
+
+       int volSize = getUnthresholdNumSeries();
+       int numThresholdedSeries = getNumSeries();
+       int volNum = getNumVolumes();
+
        FslWriteHeader(OP);
-       FslWriteVolumes(OP, qv, volNum);
+
+       float *qv = new float[volSize];
+      
+       for(int i = 1; i<= volNum; i++)
+	 {
+	   for(int j = 1; j <= volSize; j++)
+	     qv[j-1]=0;
+	   for(int j = 1; j <= numThresholdedSeries; j++)
+	     qv[int(preThresholdPositions(j))-1] = (*this)(i,j);
+	   FslWriteVolumes(OP, qv, 1);
+	 }
 
        delete [] qv;
 
diff --git a/volumeseries.h b/volumeseries.h
index 855fd83..bde1d7d 100644
--- a/volumeseries.h
+++ b/volumeseries.h
@@ -98,6 +98,7 @@ namespace MISCMATHS {
       void read(const string& fname);
       void writeAsInt(const string& fname);
       void writeAsFloat(const string& fname);
+      void writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const string& fname);
       
       void replaceMeans();
 
-- 
GitLab