diff --git a/Makefile b/Makefile
index 91d688321197f50c24832e91c24fb7b81fada624..744fb651452b52d246487b1c5311cb894aaac5d9 100644
--- a/Makefile
+++ b/Makefile
@@ -22,6 +22,8 @@ FSL_SBCA_OBJS = melhlprfns.o fsl_sbca.o
 
 FSL_MVLM_OBJS = melhlprfns.o fsl_mvlm.o
 
+FSL_BWTF_OBJS = melhlprfns.o fsl_bwtf.o
+
 FSL_REGFILT_OBJS = melhlprfns.o fsl_regfilt.o
 
 FSL_SCHURPROD_OBJS = melhlprfns.o fsl_schurprod.o
@@ -30,7 +32,7 @@ MELODIC_OBJS =  meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o
 
 TESTXFILES = test
 
-XFILES = fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod melodic
+XFILES = fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod fsl_bwtf melodic
 
 OTHERFILES =
 
@@ -44,6 +46,11 @@ ifeq ($(BUILD_MVLM), 0)
         OTHERFILES+=fsl_mvlm
 endif
 
+BUILD_BWTF = $(shell [ -f fsl_bwtf.cc ]; echo $$?)
+ifeq ($(BUILD_BPTF), 0)
+        OTHERFILES+=fsl_bwtf
+endif
+
 BUILD_SCHUR = $(shell [ -f fsl_schurprod.cc ]; echo $$?)
 ifeq ($(BUILD_SCHUR), 0)
 	OTHERFILES+=fsl_schurprod
@@ -79,6 +86,9 @@ fsl_schurprod: ${FSL_SCHURPROD_OBJS}
 fsl_mvlm: ${FSL_MVLM_OBJS}
 		$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_MVLM_OBJS} ${LIBS}
 
+fsl_bwtf: ${FSL_BWTF_OBJS}
+		$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_BWTF_OBJS} ${LIBS}
+
 fsl_regfilt: ${FSL_REGFILT_OBJS}
 	$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_REGFILT_OBJS} ${LIBS}
 
diff --git a/meldata.cc b/meldata.cc
index ce24816b0eccacceb8e6079f79bd02dd94cc0fc9..86b6f067e11c7f6c4e0cc61e8822793241a36e59 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -518,7 +518,7 @@ namespace Melodic{
 
   	//update mask
     if(opts.update_mask.value()){
-      message(endl<< "Excluding voxels with constant value ...");
+      message(endl<<  "Excluding voxels with constant value ...");
       update_mask(Mask, Data); 
       message(" done" << endl);
     }
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 87a1c46012575f187be441cc0eea22bf521aa131..2b8e769a3074b9786b548f316cae0898d57024c0 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -219,13 +219,13 @@ namespace Melodic{
     Res = zeros(nrows,nrows);
 
     if(econ>0){
-      RowVector colmeans(ncols);
-	  for (int n=1; n<=ncols; n++) {
-        colmeans(n)=0;
-        for (int m=1; m<=nrows; m++) {
-          colmeans(n)+=in(m,n);
+      ColumnVector rowmeans(nrows);
+	  for (int n=1; n<=nrows; n++) {
+        rowmeans(n)=0;
+        for (int m=1; m<=ncols; m++) {
+          rowmeans(n)+=in(n,m);
         }
-        colmeans(n)/=nrows;
+        rowmeans(n)/=ncols;
       }
       int dcol = econ;
 	  Matrix suba; 
@@ -234,13 +234,12 @@ namespace Melodic{
 	    suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
 		int scolmax = suba.Ncols();
 
-		for (int n=1; n<=scolmax; n++) {
-	        double cmean=colmeans(ctr + n - 1);
-	        for (int m=1; m<=nrows; m++) {
-	          suba(m,n)-=cmean;
+		for (int n=1; n<=nrows; n++) {
+			for (int m=1; m<=scolmax; m++) {
+	          suba(n,m)-=rowmeans(n);
 	        }
 	    }
-		
+
 	    Res += suba*suba.t() / ncols;
       }
     }
@@ -260,32 +259,32 @@ namespace Melodic{
       localweights = ones(1,ncols);
     else
       localweights = weights;
-    if(econ>0){	
-		RowVector colmeans(ncols);
-		  for (int n=1; n<=ncols; n++) {
-	        colmeans(n)=0;
-	        for (int m=1; m<=nrows; m++) {
-	          colmeans(n)+=in(m,n);
+    if(econ>0){
+      ColumnVector rowmeans(nrows);
+	  for (int n=1; n<=nrows; n++) {
+        rowmeans(n)=0;
+        for (int m=1; m<=ncols; m++) {
+          rowmeans(n)+=in(n,m);
+        }
+        rowmeans(n)/=ncols;
+      }
+      int dcol = econ;
+	  Matrix suba; 
+
+      for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
+	    suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
+		int scolmax = suba.Ncols();
+
+		for (int n=1; n<=nrows; n++) {
+			double locw =localweights(ctr + n - 1);
+			for (int m=1; m<=scolmax; m++) {
+	          suba(n,m)-=rowmeans(n);
+			suba(n,m)*=locw;
 	        }
-	        colmeans(n)/=nrows;
-	      }
-	      int dcol = econ;
-		  Matrix suba; 
-
-	      for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
-		    suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
-			int scolmax = suba.Ncols();
-
-			for (int n=1; n<=scolmax; n++) {
-		        double cmean=colmeans(ctr + n - 1);
-				double locw =localweights(ctr + n - 1);
-		        for (int m=1; m<=nrows; m++) {
-		          suba(m,n)-=cmean;
-			      suba(m,n)*=	locw;
-		        }
-		    }
-		    Res += suba*suba.t() / ncols;
-	      }
+	    }
+
+	    Res += suba*suba.t() / ncols;
+      }
     }
     else{
       Res = SP2(in,localweights);