diff --git a/volume.cc b/volume.cc
index 8929716df79d622ae14cb488723c1eb9ec9a2fdd..7bb1c0319d7fb24549afd5714833c8b867baadbb 100644
--- a/volume.cc
+++ b/volume.cc
@@ -8,7 +8,7 @@
 
 #include <iostream>
 #include <cstdlib>
-#include "avwio/avwio.h"
+#include "fslio/fslio.h"
 
 #include "newmatap.h"
 #include "newmatio.h"
@@ -129,12 +129,12 @@ namespace MISCMATHS {
        
     float fmin, fmax;
     const ColumnVector& outputvol = *this;
-    AVW* OP = AvwOpen(fname.c_str(), "wc");
+    FSLIO* OP = FslOpen(fname.c_str(), "wc");
 
-    AvwSetDim(OP,volinfo.x, volinfo.y, volinfo.z, 1);
-    AvwSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, 0);
-    AvwSetDataType(OP, DT_FLOAT);
-    AvwSetOriginator(OP, volinfo.originator);
+    FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, 1);
+    FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, 0);
+    FslSetDataType(OP, DT_FLOAT);
+    FslSetOriginator(OP, volinfo.originator);
 
     int sizeVol = outputvol.Nrows();
 
@@ -152,15 +152,15 @@ namespace MISCMATHS {
 	qv[i-1] = outputvol(i);
       }
 
-    AvwSetMinMax(OP, (short)fmin, (short)fmax);
+    FslSetMinMax(OP, (short)fmin, (short)fmax);
 
     //       fwrite(qv,sizeVol*sizeof(float),1,OP->imgfp);
     //fwrite(&OP->header,sizeof(OP->header),1,OP->hdrfp);
-    AvwWriteVolumes(OP, qv, 1);
+    FslWriteVolumes(OP, qv, 1);
 
     delete [] qv;
 
-    AvwClose(OP);
+    FslClose(OP);
   }
 
   void Volume::writeAsInt(const string& fname)
@@ -169,12 +169,12 @@ namespace MISCMATHS {
     int fmin, fmax;
 
     const ColumnVector& outputvol = *this;
-    AVW* OP = AvwOpen(fname.c_str(), "wc");
+    FSLIO* OP = FslOpen(fname.c_str(), "wc");
  
-    AvwSetDim(OP, volinfo.x, volinfo.y, volinfo.z, 1);  
-    AvwSetVoxDim(OP, volinfo.vx, volinfo.vy, volinfo.vz, 0);  
-    AvwSetDataType(OP, DT_SIGNED_SHORT);
-    AvwSetOriginator(OP, volinfo.originator);
+    FslSetDim(OP, volinfo.x, volinfo.y, volinfo.z, 1);  
+    FslSetVoxDim(OP, volinfo.vx, volinfo.vy, volinfo.vz, 0);  
+    FslSetDataType(OP, DT_SIGNED_SHORT);
+    FslSetOriginator(OP, volinfo.originator);
 
     int sizeVol = outputvol.Nrows();
 
@@ -192,34 +192,34 @@ namespace MISCMATHS {
 	qv[i-1] = (short)outputvol(i);
       }
 
-    AvwSetMinMax(OP, (short)fmin, (short)(fmax+0.9999));
+    FslSetMinMax(OP, (short)fmin, (short)(fmax+0.9999));
 
-    AvwWriteVolumes(OP, qv, 1);
+    FslWriteVolumes(OP, qv, 1);
 
     delete [] qv;
 
-    AvwClose(OP);
+    FslClose(OP);
   }
  
   void Volume::read(const string& fname)
   {
     Time_Tracer ts(string("Volume::read" + fname).c_str());
       
-    AVW* IP = AvwOpen(fname.c_str(), "r");
+    FSLIO* IP = FslOpen(fname.c_str(), "r");
     ColumnVector& output = *this;
        
     short x,y,z,v,type;
     float vx,vy,vz,tr;
     
-    AvwGetDim(IP,&x,&y,&z,&v);
-    AvwGetVoxDim(IP,&vx,&vy,&vz,&tr);
-    AvwGetOriginator(IP,volinfo.originator);
+    FslGetDim(IP,&x,&y,&z,&v);
+    FslGetVoxDim(IP,&vx,&vy,&vz,&tr);
+    FslGetOriginator(IP,volinfo.originator);
       
     volinfo.x = x; volinfo.y = y; volinfo.z = z; volinfo.v = v;
     volinfo.vx = vx; volinfo.vy = vy; volinfo.vz = vz; volinfo.tr = tr;
 
     size_t imagesize=x*y*z;
-    AvwGetDataType(IP,&type);
+    FslGetDataType(IP,&type);
   
     output.ReSize(x*y*z);
   
@@ -228,7 +228,7 @@ namespace MISCMATHS {
       case DT_SIGNED_SHORT:
 	{
 	  short* sbuffer=new short[imagesize];
-	  AvwReadVolumes(IP,sbuffer,v);
+	  FslReadVolumes(IP,sbuffer,v);
 	     
 	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
 	    {
@@ -241,7 +241,7 @@ namespace MISCMATHS {
       case DT_FLOAT:
 	{
 	  float* fbuffer=new float[imagesize];
-	  AvwReadVolumes(IP,fbuffer,v);
+	  FslReadVolumes(IP,fbuffer,v);
 
 	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
 	    {
@@ -254,7 +254,7 @@ namespace MISCMATHS {
       case DT_UNSIGNED_CHAR:
 	{
 	  unsigned char* cbuffer=new unsigned char[imagesize];
-	  AvwReadVolumes(IP,cbuffer,v);
+	  FslReadVolumes(IP,cbuffer,v);
 
 	  for(size_t j = 1; j<=(size_t)x*y*z; j++)
 	    {
@@ -265,10 +265,10 @@ namespace MISCMATHS {
 	}
 	break;
       default:
-	perror("AvwRead: DT not supported");
+	perror("FslRead: DT not supported");
       }
 
-    AvwClose(IP);
+    FslClose(IP);
     
     return;
   }
diff --git a/volumeseries.cc b/volumeseries.cc
index 4587407fd86bde3b23602be39b534ee968b78ff5..8e2a7552ecf2bde3b25c1b058b8076b5af8b4a06 100644
--- a/volumeseries.cc
+++ b/volumeseries.cc
@@ -8,7 +8,7 @@
 
 #include <iostream>
 #include <cstdlib>
-#include "avwio/avwio.h"
+#include "fslio/fslio.h"
 
 #include "newmatap.h"
 #include "newmatio.h"
@@ -97,7 +97,8 @@ namespace MISCMATHS {
 	       cerr << "j = " << j+1 << endl;
 	     }
 	   */
-	   if(m > thresh && s > 0.0)
+
+	   if(m > thresh && s > 1e-10)
 	     {
 	       j++;	
 	       preThresholdPositions(j) = i;
@@ -133,21 +134,21 @@ namespace MISCMATHS {
     {
       Time_Tracer ts(string("VolumeSeries::read-" + fname).c_str());
 
-      AVW* IP = AvwOpen(fname.c_str(), "r");
+      FSLIO* IP = FslOpen(fname.c_str(), "r");
       Matrix& output = *this;
 
       short x,y,z,v,type;
       float vx,vy,vz,tr;
 
-      AvwGetDim(IP,&x,&y,&z,&v);
-      AvwGetVoxDim(IP,&vx,&vy,&vz,&tr);
-      AvwGetOriginator(IP,volinfo.originator);
+      FslGetDim(IP,&x,&y,&z,&v);
+      FslGetVoxDim(IP,&vx,&vy,&vz,&tr);
+      FslGetOriginator(IP,volinfo.originator);
 
       volinfo.x = x; volinfo.y = y; volinfo.z = z; volinfo.v = v;
       volinfo.vx = vx; volinfo.vy = vy; volinfo.vz = vz; volinfo.tr = tr;
 
       size_t imagesize=x*y*z*v;
-      AvwGetDataType(IP,&type);
+      FslGetDataType(IP,&type);
   
       output.ReSize(v,x*y*z);
   
@@ -157,7 +158,7 @@ namespace MISCMATHS {
 	  {
 	    short* sbuffer=new short[imagesize];
 	
-	    AvwReadVolumes(IP, sbuffer,v);
+	    FslReadVolumes(IP, sbuffer,v);
 	
 	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
 	    size_t volstart = 1;
@@ -177,7 +178,7 @@ namespace MISCMATHS {
 	case DT_FLOAT:
 	  {
 	    float* fbuffer=new float[imagesize];
-	    AvwReadVolumes(IP,fbuffer,v);
+	    FslReadVolumes(IP,fbuffer,v);
 
 	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
 	    size_t volstart = 1;
@@ -195,7 +196,7 @@ namespace MISCMATHS {
 	case DT_UNSIGNED_CHAR:
 	  {
 	    unsigned char* cbuffer=new unsigned char[imagesize];
-	    AvwReadVolumes(IP,cbuffer,v);
+	    FslReadVolumes(IP,cbuffer,v);
 
 	    size_t volsize = volinfo.x*volinfo.y*volinfo.z;
 	    size_t volstart = 1;
@@ -211,10 +212,10 @@ namespace MISCMATHS {
 	  }
 	  break;
 	default:
-	  perror("AvwRead: DT not supported");
+	  perror("FslRead: DT not supported");
 	}
 
-      AvwClose(IP);
+      FslClose(IP);
 
       return;
     }
@@ -230,12 +231,12 @@ namespace MISCMATHS {
     {
        Time_Tracer ts(string("VolumeSeries::writeAsFloat" + fname).c_str());
 
-       AVW* OP = AvwOpen(fname.c_str(), "wc");
+       FSLIO* OP = FslOpen(fname.c_str(), "wc");
 
-       AvwSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
-       AvwSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
-       AvwSetDataType(OP, DT_FLOAT);
-       AvwSetOriginator(OP, volinfo.originator);
+       FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
+       FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
+       FslSetDataType(OP, DT_FLOAT);
+       FslSetOriginator(OP, volinfo.originator);
 
        int volStart = 1;
        int volSize = getNumSeries();
@@ -252,11 +253,11 @@ namespace MISCMATHS {
 	     }
 	 }
             
-       AvwWriteVolumes(OP, qv, volNum);
+       FslWriteVolumes(OP, qv, volNum);
 
        delete [] qv;
 
-       AvwClose(OP);
+       FslClose(OP);
        
     }