diff --git a/Makefile b/Makefile
index ee2a88733d90ca61d32987b4dd21c850a1bf7755..e19ee8c9ed57b227ae0922d4f0340dee4faa860d 100644
--- a/Makefile
+++ b/Makefile
@@ -24,8 +24,10 @@ XFIBRES=xfibres
 RV=replacevols
 MDV=make_dyadic_vectors
 FMO=fdt_matrix_ops
+INDEXER=indexer
 TEST=testfile
 
+
 DTIFITOBJS=dtifit.o dtifitOptions.o
 CCOPSOBJS=ccops.o ccopsOptions.o
 PTOBJS=probtrack.o probtrackOptions.o pt_alltracts.o pt_matrix.o pt_seeds_to_targets.o pt_simple.o pt_twomasks.o pt_matrix_mesh.o
@@ -40,12 +42,14 @@ XFIBOBJS=xfibres.o xfibresoptions.o
 RVOBJS=replacevols.o
 MDVOBJS=make_dyadic_vectors.o
 FMOOBJS=fdt_matrix_ops.o
+INDEXEROBJS=indexer.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 ccops probtrack find_the_biggest medianfilter diff_pvm make_dyadic_vectors proj_thresh
-FXFILES = reord_OM sausages replacevols fdt_matrix_ops probtrackx xfibres
+FXFILES = reord_OM sausages replacevols fdt_matrix_ops probtrackx xfibres indexer
 
 
 RUNTCLS = Fdt
@@ -94,6 +98,9 @@ ${MDV}:    	${MDVOBJS}
 ${FMO}:    	${FMOOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${FMOOBJS} ${DLIBS}
 
+${INDEXER}:    	${INDEXEROBJS}
+		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${INDEXEROBJS} ${DLIBS}
+
 ${TEST}:    	${TESTOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${TESTOBJS} ${DLIBS}
 
diff --git a/indexer.cc b/indexer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..60823cf40c2ae4aa430339132f5a82a8d311529f
--- /dev/null
+++ b/indexer.cc
@@ -0,0 +1,112 @@
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+
+
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace Utilities;
+//using namespace NEWMAT;
+//////////////////////////
+/////////////////////////
+
+void read_masks(vector<string>& masks, const string& filename){
+  ifstream fs(filename.c_str());
+  string tmp;
+  if(fs){
+    fs>>tmp;
+    while(!fs.eof()){
+      masks.push_back(tmp);
+      fs>>tmp;
+    }
+  }
+  else{
+    cerr<<filename<<" does not exist"<<endl;
+    exit(0);
+  }
+}
+
+
+ColumnVector vox_to_vox(const ColumnVector& xyz1,const ColumnVector& dims1,const ColumnVector& dims2,const Matrix& xfm){
+  ColumnVector xyz1_mm(4),xyz2_mm,xyz2(3);
+  xyz1_mm<<xyz1(1)*dims1(1)<<xyz1(2)*dims1(2)<<xyz1(3)*dims1(3)<<1;
+  xyz2_mm=xfm*xyz1_mm;
+  xyz2_mm=xyz2_mm/xyz2_mm(4);
+  xyz2<<xyz2_mm(1)/dims2(1)<<xyz2_mm(2)/dims2(2)<<xyz2_mm(3)/dims2(3);
+  return xyz2;
+}
+
+
+int main ( int argc, char **argv ){
+  if(argc<6){
+    cerr<<"Usage: indexer <mask_for_going_ahead> <file_of_volume_names> <x> <y> <z>"<<endl;
+    cerr<<endl;
+    cerr<< "<x> <y> <z> in mni mm coords"<<endl;
+    return 0;
+  }
+ 
+  vector<string> masknames;
+  volume<int> inside_mask;
+  read_volume(inside_mask,argv[1]);
+  read_masks(masknames,argv[2]);
+  vector<volume<float> > cortex_masks;
+  volume<float> tmpcort;
+  for( unsigned int m = 0; m < masknames.size(); m++ ){
+    read_volume(tmpcort,masknames[m]);
+    cortex_masks.push_back(tmpcort);
+  }
+  
+  float x_mm=atof(argv[3]);
+  float y_mm=atof(argv[4]);
+  float z_mm=atof(argv[5]);
+  float x_orig_mm=90;
+  float y_orig_mm=124;
+  float z_orig_mm=72;
+  int x_roi_start_vox=45;
+  int y_roi_start_vox=70;
+  int z_roi_start_vox=50;
+  int numsubjects=11;
+  float xvox=(x_mm+x_orig_mm)/cortex_masks[0].xdim()-x_roi_start_vox;
+  float yvox=(y_mm+y_orig_mm)/cortex_masks[0].ydim()-y_roi_start_vox;
+  float zvox=(z_mm+z_orig_mm)/cortex_masks[0].zdim()-z_roi_start_vox;
+  
+  if(inside_mask((int)round(xvox),(int)round(yvox),(int)round(zvox))==0){
+    cout<<"Sorry - Your input is not in our defined thalamus"<<endl;
+    for(unsigned int i=0;i<masknames.size();i++){
+      cout<<0<<endl;
+    }
+  }  
+  else{
+    cout<<"Input inside thalamus"<<endl;
+    for(unsigned int i=0;i<masknames.size();i++){
+      if(cortex_masks[i].interpolate(xvox,yvox,zvox)==0){
+	cout<<std::setprecision(0)<<std::fixed;
+      }
+      else{
+	cout<<std::setprecision(2)<<std::fixed;
+      }
+      cout<<cortex_masks[i].interpolate(xvox,yvox,zvox)/numsubjects<<endl;
+    }
+  }
+  return 0;
+}
+
+
+ 
+
+
+
+
+
+
+
+
+
+
+
+
+
+