diff --git a/miscmaths.cc b/miscmaths.cc
index f38dc19fae95a60595c9b153c68076dae3ac772a..8b39669963dbeb138e466caa5816345f8c0c507f 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -1177,6 +1177,31 @@ float percentile(const ColumnVector& in, float p)
   return interp1(xx,yy,p);
 }
 
+ReturnMatrix quantile(const Matrix& in, int which)
+{
+  int num = in.Ncols();
+  Matrix res(1,num);
+  for (int ctr=1; ctr<=num; ctr++){
+    ColumnVector tmp = in.Column(ctr);
+    res(1,ctr) = quantile(tmp,which);
+  }
+  res.Release();
+  return res;
+}
+
+ReturnMatrix  percentile(const Matrix& in, float p)
+{
+  int num = in.Ncols();
+  Matrix res(1,num);
+  for (int ctr=1; ctr<=num; ctr++){
+    ColumnVector tmp = in.Column(ctr);
+    res(1,ctr) = percentile(tmp,p);
+  }
+  res.Release();
+  return res;
+}
+
+
 
 void cart2sph(const ColumnVector& dir, float& th, float& ph)
 {
@@ -1239,7 +1264,7 @@ void cart2sph(const vector<ColumnVector>& dir,ColumnVector& th,ColumnVector& ph)
     th.ReSize(dir.size());
   }
 
-  if(ph.Nrows()!=dir.size()){
+  if(ph.Nrows()!=(int)dir.size()){
     ph.ReSize(dir.size());
   }
   //double _2pi=2*M_PI;
diff --git a/miscmaths.h b/miscmaths.h
index 9a5ac2e01cfdc235b41ad1e50385b4b97e4e0f4b..1df1bcb1bc5ab91c114b6be4ad4a318fe43cc231 100644
--- a/miscmaths.h
+++ b/miscmaths.h
@@ -1,4 +1,4 @@
-/*  miscmaths.h
+*  miscmaths.h
 
     Mark Jenkinson & Mark Woolrich & Christian Beckmann & Tim Behrens, FMRIB Image Analysis Group
 
@@ -50,6 +50,7 @@ namespace MISCMATHS {
   ReturnMatrix read_ascii_matrix(const string& filename);
   ReturnMatrix read_vest(string p_fname);
   ReturnMatrix read_binary_matrix(const string& filename);
+  ReturnMatrix read_matrix(const string& filename);
 
   int write_ascii_matrix(const Matrix& mat, const string& filename, 
 			 int precision=-1);
@@ -162,13 +163,19 @@ namespace MISCMATHS {
 			     const Matrix& qform_mat, int qform_code,
 			     int& icode, int& jcode, int& kcode);
     
-  
+  // 1D lookup table with linear interpolation
+  float interp1(const ColumnVector& x, const ColumnVector& y, float xi);
+
   float quantile(const ColumnVector& in, int which);
   float percentile(const ColumnVector& in, float p);
-  float interp1(const ColumnVector& x, const ColumnVector& y, float xi);
-  inline float median(const ColumnVector& x){ return quantile(x,2);}
+  inline float median(const ColumnVector& in){ return quantile(in,2);}
   inline float iqr(const ColumnVector &in) { return quantile(in,3) - quantile(in,1); }
 
+  ReturnMatrix quantile(const Matrix& in, int which);
+  ReturnMatrix percentile(const Matrix& in, float p);
+  inline ReturnMatrix median(const Matrix& in){ return quantile(in,2);}
+  inline ReturnMatrix iqr(const Matrix& in){ Matrix res = quantile(in,3) - quantile(in,1); res.Release(); return res;}
+
   void cart2sph(const ColumnVector& dir, float& th, float& ph);// cartesian to sperical polar coordinates
   void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph);//ditto
   void cart2sph(const vector<ColumnVector>& dir,ColumnVector& th,ColumnVector& ph);// same but in a vector