diff --git a/fsldecorr4d.cc b/fsldecorr4d.cc
index 13b49acfe64e9c32612d138159544ba9e1eaa174..70530189210615105766b97acd5f4a2caac5c356 100644
--- a/fsldecorr4d.cc
+++ b/fsldecorr4d.cc
@@ -89,7 +89,7 @@ int do_work(int argc, char *argv[])
   mask.binarise(1e-8);  // arbitrary "0" threshold
 
   Matrix tc;
-  tc = read_matrix(tcname.value());
+  tc = read_ascii_matrix(tcname.value());
 
   if (tc.Nrows() != vin.tsize()) {
     cerr << "ERROR: Different number of timepoints in timecourse file and input volume." << endl;
@@ -113,14 +113,14 @@ int do_work(int argc, char *argv[])
       for (int z=smaps[n-1].minz(); z<=smaps[n-1].maxz(); z++) {
 	for (int y=smaps[n-1].miny(); y<=smaps[n-1].maxy(); y++) {
 	  for (int x=smaps[n-1].minx(); x<=smaps[n-1].maxx(); x++) {
-	    XY(n) += smaps(x,y,z,n-1) * tc(t+1,n) * vin(x,y,z,t);
+	    XY(n) += smaps(x,y,z,n-1) * tc(t+1,n) * vin(x,y,z,t) * mask(x,y,z);
 	  }
 	}
       }
     }
     for (int m=1; m<=n; m++) {
       volume<float> tmp;
-      tmp = smaps[n-1] * smaps[m-1];
+      tmp = smaps[n-1] * smaps[m-1] * mask;
       XX(m,n) = tmp.sum();
       double tval=0.0;
       for (int t=0; t<nt; t++) {
@@ -142,7 +142,7 @@ int do_work(int argc, char *argv[])
       for (int z=smaps[n-1].minz(); z<=smaps[n-1].maxz(); z++) {
 	for (int y=smaps[n-1].miny(); y<=smaps[n-1].maxy(); y++) {
 	  for (int x=smaps[n-1].minx(); x<=smaps[n-1].maxx(); x++) {
-	    vout(x,y,z,t) -= Beta(n) * smaps(x,y,z,n-1) * tc(t+1,n);
+	    vout(x,y,z,t) -= Beta(n) * smaps(x,y,z,n-1) * tc(t+1,n) * mask(x,y,z);
 	  }
 	}
       }