diff --git a/Makefile b/Makefile
index 97fb300085d1ec9d704e87d8a38382b29b2a4a9a..c1e4932f27809f05541a1dc60ddc62305df2989e 100644
--- a/Makefile
+++ b/Makefile
@@ -22,6 +22,7 @@ DIFF_PVM=diff_pvm
 XFIBRES=xfibres
 RV=replacevols
 MDV=make_dyadic_vectors
+FMO=fdt_matrix_ops
 TEST=testfile
 
 DTIFITOBJS=dtifit.o dtifitOptions.o
@@ -36,12 +37,13 @@ DIFF_PVMOBJS=diff_pvm.o diff_pvmoptions.o
 XFIBOBJS=xfibres.o xfibresoptions.o
 RVOBJS=replacevols.o
 MDVOBJS=make_dyadic_vectors.o
+FMOOBJS=fdt_matrix_ops.o
 TESTOBJS=testfile.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 replacevols
+FXFILES = reord_OM sausages replacevols fdt_matrix_ops
 
 
 RUNTCLS = Fdt
@@ -84,6 +86,9 @@ ${RV}:    	${RVOBJS}
 ${MDV}:    	${MDVOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${MDVOBJS} ${DLIBS}
 
+${FMO}:    	${FMOOBJS}
+		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${FMOOBJS} ${DLIBS}
+
 ${TEST}:    	${TESTOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${TESTOBJS} ${DLIBS}
 
diff --git a/fdt_matrix_ops.cc b/fdt_matrix_ops.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a71547fabffcc822813e94118dbb93f870bc6a98
--- /dev/null
+++ b/fdt_matrix_ops.cc
@@ -0,0 +1,139 @@
+/*  Copyright (C) 2004 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include "newimage/newimageall.h"
+#include <vector>
+#include <algorithm>
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace NEWMAT;
+
+
+int main ( int argc, char **argv ){
+  if(argc<5){
+    cout<<"usage: fdt_matrix_ops <matrix1> <seedmask1> <matrix2> <seedmask2> ... [inclusion mask] <output>"<<endl;
+    cout<<"creates one big uber-matrix from lots of cute wee dinky ones"<<endl;
+    cout<<"If seedmasks overlap, it just takes the values from the first of them"<<endl;
+    cout<<"If you specify an incluson mask (optional), only voxels that are"<<endl;
+    cout<<"inside this mask will be included in the output"<<endl;
+    exit(0);
+  }
+  int var=argc-1;
+  bool incmaskyn = (float(var)/2.0)==int(float(var)/2.0);
+  int Nmats;
+  if(incmaskyn) Nmats=var/2-1;
+  else Nmats=(var-1)/2;
+  
+  vector<volume<int>* > masks;
+  vector<volume<float>* > mats;
+  vector<volume<int>* > coords;
+  volume<int> incmask;
+  volume<int> totalmask;
+  for(int i=0;i<Nmats;i++){
+    volume<float> *tmp= new volume<float>;
+    volume<int> *tmpcoords= new volume<int>;
+    volume<int> *tmpmask= new volume<int>;
+    read_volume(*tmp,argv[ 2*i + 1 ]);
+    read_volume(*tmpcoords,"coords_for_" + string(argv[ 2*i + 1 ]) );
+    read_volume(*tmpmask,argv[ 2*(i + 1) ]);
+    mats.push_back(tmp);
+    masks.push_back(tmpmask);
+    if(i==0) totalmask=*tmpmask;
+    else totalmask=totalmask+ *tmpmask;
+
+}
+  if(incmaskyn){
+    read_volume(incmask,argv[var-1]);
+    incmask=incmask*totalmask;
+  }else{
+    incmask=totalmask;
+  }
+  
+  vector<volume<int>* >lookups;
+  for(int i=0;i<Nmats;i++){
+    volume<int> *lu = new volume<int>;
+    *lu=*masks[i];int conrow=0;
+    for(int z=0;z<(*lu).zsize();z++){
+      for(int y=0;y<(*lu).ysize();y++){
+	for(int x=0;x<(*lu).xsize();x++){
+	  if((*masks[i])(x,y,z)>0){
+	    (*lu)(x,y,z)=conrow;
+	    conrow++;
+	  }
+	}
+      }
+    }
+    lookups.push_back(lu);
+  }
+  
+  
+  int nvoxels=0;
+
+  for(int z=0;z<incmask.zsize();z++) {
+    for(int y=0;y<incmask.ysize();y++){
+      for(int x=0;x<incmask.xsize();x++){
+	  if(incmask(x,y,z)>0) nvoxels++;
+      }
+    }
+  }
+  
+  volume<float> output(nvoxels,(*mats[0]).ysize(),0);
+  volume<float> outcoords(nvoxels,3,0);
+  
+  int newrow=0;
+  for(int z=0;z<incmask.zsize();z++) {
+    for(int y=0;y<incmask.ysize();y++){
+      for(int x=0;x<incmask.xsize();x++){
+	if(incmask(x,y,z)>0){
+	  bool found=false;
+	  for(unsigned int i=0;i<mats.size();i++){
+	    if(!found){
+	      if((*masks[i])(x,y,z)>0){
+		int oldrow=(*lookups[i])(x,y,z);
+		for(int col=0;col<(*mats[i]).ysize();col++){
+		  output(newrow,col,0)=(*mats[i])(oldrow,col,0);
+		}
+		outcoords(newrow,0,0)=(*coords[i])(oldrow,0,0);
+		outcoords(newrow,1,0)=(*coords[i])(oldrow,1,0);
+		outcoords(newrow,2,0)=(*coords[i])(oldrow,2,0);
+		
+		found=true;
+		newrow++;
+	      }
+	    }
+	    
+	  }
+	  
+	}
+      }
+    }
+  }
+
+  save_volume(output,argv[argc-1]);
+  save_volume(outcoords,"coords_for_"+string(argv[argc-1]));
+ return 0;
+}
+ 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+