diff --git a/Makefile b/Makefile
index 9016deff1a634f575cec76147b183ebffdb176b4..2251bd3f045d8379206df2fcc54e8e5e5f70bd1c 100644
--- a/Makefile
+++ b/Makefile
@@ -35,7 +35,7 @@ MDVOBJS=make_dyadic_vectors.o
 SCRIPTS = eddy_correct bedpost bedpost_proc bedpost_cleanup bedpost_kill_all bedpost_kill_pid zeropad bedpost_datacheck
 FSCRIPTS=correct_and_average ocmr_preproc
 XFILES = dtifit probtrack find_the_biggest medianfilter diff_pvm make_dyadic_vectors proj_thresh
-FXFILES = reord_OM sausages
+FXFILES = reord_OM sausages replacevols
 
 
 RUNTCLS = Fdt
diff --git a/dtifit.cc b/dtifit.cc
index 73f7532fe5d67c77775398f09095260bec0bfe45..ae84c19015d96c221bd24df0e69fea70662c2d13 100644
--- a/dtifit.cc
+++ b/dtifit.cc
@@ -135,6 +135,7 @@ void tensorfit(DiagonalMatrix& Dd,ColumnVector& evec1,ColumnVector& evec2,Column
       logS(i)=(S(i)/s0)>0.01 ? log(S(i)):log(0.01*s0);
     }
   Dvec = -pinv(Amat)*logS;
+  cout << Dvec;
   s0=exp(-Dvec(7));
   if(s0<S.Sum()/S.Nrows()){ s0=S.Sum()/S.Nrows();  }
   tens = vec2tens(Dvec);
diff --git a/pt_matrix.cc b/pt_matrix.cc
index 22a1569c549fe60a5e849dbb9f4c6f2ff04fcde2..4c78d2584f33564ebec40b7f76395810eabe69b2 100644
--- a/pt_matrix.cc
+++ b/pt_matrix.cc
@@ -626,17 +626,16 @@ void maskmatrix(){
     ColumnVector clustsums(maxclusternum),clustmaxes(maxclusternum);
     clustsums=0;clustmaxes=0;
     for(int Sz=Seeds.minz();Sz<=Seeds.maxz();Sz++){
-      cout<<Sz<<endl;
+      cout<<seedclust<<" "<<Sz<<endl;
       for(int Sy=Seeds.miny();Sy<=Seeds.maxy();Sy++){
 	for(int Sx=Seeds.minx();Sx<=Seeds.maxx();Sx++){
-	  if(Seeds(Sx,Sy,Sz)>0){
+	  if(Seeds(Sx,Sy,Sz)==seedclust){
 
 	    clustsize++;
-	    ColumnVector flags(maxclusternum);  flags=0;
+	    ColumnVector flags(maxclusternum); 
 	    ColumnVector clustcounts(maxclusternum);  clustcounts=0;
 
 
-
 	    ColumnVector xyz_seeds(3),dim_seeds(3),xyz_dti;
 	    xyz_seeds << Sx << Sy << Sz;
 	    dim_seeds <<Seeds.xdim()<<Seeds.ydim()<<Seeds.zdim();	  
@@ -645,7 +644,7 @@ void maskmatrix(){
 	    Particle part(0,0,0,0,0,0,opts.steplength.value(),mask.xdim(),mask.ydim(),mask.zdim(),false);
 	    hitcount=0;
 	    for( int p = 0; p < nparticles ; p++ ){
-	    
+	     flags=0;
 	      for(int direc=1;direc<=2;direc++){
 		x=xst;y=yst;z=zst;
 		part.change_xyz(x,y,z);