*** empty log message ***

Makefile

@@ -11,6 +11,7 @@ DLIBS = -lmeshclass -lbint -lnewimage -lutils -lmiscmaths  -lnewmat -lfslio -lni
 
 
 DTIFIT=dtifit
+CCOPS=ccops
 PT=probtrack
 PTX=probtrackx
 FTB=find_the_biggest
@@ -26,6 +27,7 @@ FMO=fdt_matrix_ops
 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
 PTXOBJS=probtrackx.o probtrackxOptions.o streamlines.o ptx_simple.o ptx_seedmask.o ptx_twomasks.o
 FTBOBJS=find_the_biggest.o
@@ -42,7 +44,7 @@ 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
+XFILES = dtifit ccops probtrack find_the_biggest medianfilter diff_pvm make_dyadic_vectors proj_thresh
 FXFILES = reord_OM sausages replacevols fdt_matrix_ops
 
 
@@ -68,6 +70,9 @@ ${MED}:    	${MEDOBJS}
 ${DTIFIT}:    	${DTIFITOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${DTIFITOBJS} ${DLIBS}
 
+${CCOPS}:    	${CCOPSOBJS}	
+		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${CCOPSOBJS} ${DLIBS}
+
 ${ROM}:    	${ROMOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${ROMOBJS} ${DLIBS}
 

ccops.cc

+/*  Copyright (C) 2004 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include "newimage/newimageall.h"
+#include "ccopsOptions.h"
+#include <vector>
+#include <algorithm>
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace NEWMAT;
+using namespace CCOPS;
+  void spect_reord(SymmetricMatrix& A,ColumnVector& r,ColumnVector& y){
+    SymmetricMatrix Q=-A;
+      DiagonalMatrix t(Q.Nrows());
+      t=0;
+      DiagonalMatrix D;
+      Matrix V;
+      for(int i=1;i<=Q.Nrows();i++){
+        float rowsum1=0, rowsum2=0;
+        for(int j=1;j<=Q.Ncols();j++){
+          if(i!=j) rowsum1+=Q(i,j);
+          rowsum2+=A(i,j);
+        }
+          Q(i,i)=-rowsum1;
+          t(i)=1/sqrt(rowsum2);
+      }
+      Q << t*Q*t;
+      EigenValues(Q,D,V);
+      vector<pair<float,int> > myvec;
+      vector<pair<float,int> > myvec2;
+      
+      for(int i=1;i<=D.Nrows();i++){
+	pair<float,int> mypair;
+	mypair.first=D(i);
+	mypair.second=i;
+	myvec.push_back(mypair);
+      }
+      sort(myvec.begin(),myvec.end());
+      int ind=myvec[1].second; // index for second eigenval
+      
+      ColumnVector v2scale(V.Nrows());
+      for(int i=1;i<=V.Nrows();i++){
+	v2scale(i)=V(i,ind); //second eigvec
+      }
+      v2scale=t*v2scale; //scale it
+      
+      
+      for(int i=1;i<=D.Nrows();i++){
+	pair<float,int> mypair;
+	mypair.first=v2scale(i);
+	mypair.second=i;
+	myvec2.push_back(mypair);
+      }
+      //myvec2 contains scaled second eigenvector and index for sorting.
+      
+      
+      sort(myvec2.begin(),myvec2.end());
+      r.ReSize(D.Nrows());
+      y.ReSize(D.Nrows());
+      
+      for(int i=1;i<=D.Nrows();i++){
+	y(i)=myvec2[i-1].first;
+	r(i)=myvec2[i-1].second;
+      }
+      
+  } 
+
+
+void rem_zrowcol(const Matrix& myOMmat,const Matrix& coordmat,const Matrix& tractcoordmat,const bool coordbool,const bool tractcoordbool,Matrix& newOMmat,Matrix& newcoordmat, Matrix& newtractcoordmat)
+{
+ 
+ vector<int> zerorows;
+ vector<int> zerocols;
+ int dimsum =0;
+
+
+
+ cerr<< "Checking for all zero rows"<<endl;
+ for(int i=1;i<=myOMmat.Nrows();i++){
+   dimsum=0;
+   for(int j=1;j<=myOMmat.Ncols();j++){
+     dimsum+=(int)myOMmat(i,j);
+   }
+   if(dimsum==0){zerorows.push_back(i);}
+ }
+
+
+
+ cerr<< "Checking for all zero cols"<<endl;
+ for(int j=1;j<=myOMmat.Ncols();j++){
+   dimsum=0;
+   for(int i=1;i<=myOMmat.Nrows();i++){
+     dimsum+=(int)myOMmat(i,j);
+   }
+   if(dimsum==0){zerocols.push_back(j);}
+ }
+
+
+ newOMmat.ReSize(myOMmat.Nrows()-zerorows.size(),myOMmat.Ncols()-zerocols.size());
+ if(coordbool){
+   newcoordmat.ReSize(coordmat.Nrows()-zerorows.size(),3);
+ }
+ if(tractcoordbool){
+   newtractcoordmat.ReSize(tractcoordmat.Nrows()-zerocols.size(),3);
+ }
+
+ int zrowcounter=0,zcolcounter=0,nzrowcounter=1,nzcolcounter=1;
+ cerr<<"Forming New Matrix"<<endl;
+ for(int j=1;j<=myOMmat.Ncols();j++){
+   zrowcounter=0;
+   nzrowcounter=1;
+
+   if(zerocols.size()>0){ //Are there any Zero Columns
+     if(zerocols[zcolcounter]!=j){  // Only add a col if it's not the next zcol
+       for(int i=1;i<=myOMmat.Nrows();i++){
+	 if(zerorows.size()>0){ //Are there any zero rows?
+	   if(zerorows[zrowcounter]!=i){ //Only add a row if it's not the next zrow
+	     newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);
+	     nzrowcounter++;
+	   }
+	   else{zrowcounter++;}
+	 }
+	 else{newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);nzrowcounter++;}//No Zero Rows
+       }
+       nzcolcounter++;
+     }
+     else{zcolcounter++;} //move onto next z col
+   }
+   
+   else{  //No zero Columns
+     for(int i=1;i<=myOMmat.Nrows();i++){
+       if(zerorows.size()>0){
+	 if(zerorows[zrowcounter]!=i){ //Only add a row if it's not the next zrow
+	   newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);
+	   nzrowcounter++;
+	 }
+	 else{zrowcounter++;}
+       }
+       else{newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);nzrowcounter++;}
+     }
+     nzcolcounter++;
+   }
+
+}
+
+
+ if(coordbool){
+   cerr<<"Updating Seed Coordinates"<<endl;
+   zrowcounter=0;nzrowcounter=1;
+   if(zerorows.size()>0){//Are there any zero rows?
+     for(int i=1;i<=coordmat.Nrows();i++){
+       if(zerorows[zrowcounter]!=i){
+	 newcoordmat(nzrowcounter,1)=coordmat(i,1);
+	 newcoordmat(nzrowcounter,2)=coordmat(i,2);
+	 newcoordmat(nzrowcounter,3)=coordmat(i,3);
+	 nzrowcounter++;
+       }
+       else{zrowcounter++;}
+     }
+   } 
+   else{//No zero Rows
+     newcoordmat=coordmat;
+   }
+}
+
+ if(tractcoordbool){
+   cerr<<"Updating Tract Coordinates"<<endl;
+   zcolcounter=0;nzcolcounter=1;
+   if(zerocols.size()>0){//Are there any zero cols?
+     for(int i=1;i<=tractcoordmat.Nrows();i++){
+       if(zerocols[zcolcounter]!=i){
+	 newtractcoordmat(nzcolcounter,1)=tractcoordmat(i,1);
+	 newtractcoordmat(nzcolcounter,2)=tractcoordmat(i,2);
+	 newtractcoordmat(nzcolcounter,3)=tractcoordmat(i,3);
+	 nzcolcounter++;
+       }
+       else{zcolcounter++;}
+     }
+   }
+   else{//No zero Cols 
+     newtractcoordmat=tractcoordmat;
+   }
+ }
+
+}
+
+
+
+
+
+
+
+
+
+void rem_cols(Matrix& myOMmat,Matrix& tractcoordmat,const bool tractcoordbool,const vector<int>& excl_cols)
+{
+ 
+
+ Matrix newOMmat,newtractcoordmat;
+ newOMmat.ReSize(myOMmat.Nrows(),myOMmat.Ncols()-excl_cols.size());
+
+ if(tractcoordbool){
+   newtractcoordmat.ReSize(tractcoordmat.Nrows()-excl_cols.size(),3);
+ }
+
+ int zrowcounter=0,zcolcounter=0,nzcolcounter=1,nzrowcounter=1;
+ vector<int> zerorows;
+
+ for(int j=1;j<=myOMmat.Ncols();j++){
+   zrowcounter=0;
+   nzrowcounter=1;   
+
+   if(excl_cols.size()>0){ //Are there any excl Columns
+     if(excl_cols[zcolcounter]!=j){  // Only add a col if it's not the next zcol
+       for(int i=1;i<=myOMmat.Nrows();i++){
+	 if(zerorows.size()>0){ //Are there any excl rows?
+	   if(zerorows[zrowcounter]!=i){ //Only add a row if it's not the next zrow
+	     newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);
+	     nzrowcounter++;
+	   }
+	   else{zrowcounter++;}
+	 }
+	 else{newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);nzrowcounter++;}//No Zero Rows
+       }
+       nzcolcounter++;
+     }
+     else{zcolcounter++;} //move onto next z col
+   }
+   else{  //No zero Columns
+     for(int i=1;i<=myOMmat.Nrows();i++){
+       if(zerorows.size()>0){
+	 if(zerorows[zrowcounter]!=i){ //Only add a row if it's not the next zrow
+	   newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);
+	   nzrowcounter++;
+	 }
+	 else{zrowcounter++;}
+       }
+       else{newOMmat(nzrowcounter,nzcolcounter)=myOMmat(i,j);nzrowcounter++;}
+     }
+     nzcolcounter++;
+   }
+}
+
+ 
+ if(tractcoordbool){
+   zcolcounter=0;nzcolcounter=1;
+   if(excl_cols.size()>0){//Are there any zero cols?
+     for(int i=1;i<=tractcoordmat.Nrows();i++){
+       if(excl_cols[zcolcounter]!=i){
+	 newtractcoordmat(nzcolcounter,1)=tractcoordmat(i,1);
+	 newtractcoordmat(nzcolcounter,2)=tractcoordmat(i,2);
+	 newtractcoordmat(nzcolcounter,3)=tractcoordmat(i,3);
+	 nzcolcounter++;
+       }
+       else{zcolcounter++;}
+     }
+   }
+   else{//No zero Cols 
+     newtractcoordmat=tractcoordmat;
+   }
+ }
+
+ myOMmat = newOMmat;
+ tractcoordmat = newtractcoordmat;
+  
+}
+
+
+
+int main ( int argc, char **argv ){
+  ccopsOptions& opts = ccopsOptions::getInstance();
+  int success=opts.parse_command_line(argc,argv);
+  string ip=opts.inmatrix.value();
+  make_basename(ip);
+  
+  ColumnVector y1,r1,y2,r2;
+  volume<int> myOM;
+  volume<int> coordvol;
+  volume<int> tractcoordvol;
+  bool coordbool=false,tractcoordbool=false;
+  read_volume(myOM,ip);
+  
+  Matrix myOMmat(myOM.xsize(),myOM.ysize());
+  Matrix mycoordmat,mytractcoordmat;
+  Matrix newOMmat,newcoordmat,newtractcoordmat;
+  
+  for(int j=0;j<myOM.ysize();j++){
+    for(int i=0;i<myOM.xsize();i++){
+      if(opts.bin.value()==0)
+	myOMmat(i+1,j+1)=float(myOM(i,j,0)); 
+      else{
+	if(myOM(i,j,0)>opts.bin.value()){
+	  myOMmat(i+1,j+1)=1.0f;
+	}
+	else{
+	  myOMmat(i+1,j+1)=0.0f;
+	}
+      }
+    }
+  }
+  // write_ascii_matrix(newOMmat.t(),"preprecock");
+  //Checking for and loading up Seed Coordinates
+
+  string coordname="coords_for_"+ip;
+  if(fsl_imageexists(coordname)){
+    read_volume(coordvol,coordname);
+    coordbool=true;
+    mycoordmat.ReSize(coordvol.xsize(),coordvol.ysize());
+    for(int j=0;j<coordvol.ysize();j++){
+      for(int i=0;i<coordvol.xsize();i++){
+	mycoordmat(i+1,j+1)=coordvol(i,j,0);
+      }
+    }
+  }
+  else{
+    cerr<<"Seed Space Coordinate File Not present - Ignoring"<<endl;
+  }
+
+  //Checking For and Loading Up Tract coordinates
+  string trcoordname="tract_space_coords_for_"+ip;
+  if(fsl_imageexists(trcoordname)){
+    read_volume(tractcoordvol,trcoordname);
+    tractcoordbool=true;
+    mytractcoordmat.ReSize(tractcoordvol.xsize(),tractcoordvol.ysize());
+    for(int j=0;j<tractcoordvol.ysize();j++){
+      for(int i=0;i<tractcoordvol.xsize();i++){
+	mytractcoordmat(i+1,j+1)=tractcoordvol(i,j,0);
+      }
+    }
+  }
+  else{
+    cerr<<"Tract Space Coordinate File Not present - Ignoring"<<endl;
+  }
+ 
+
+
+  // If user specifies an exclusion mask in tract space. 
+  // work out which columns in the matrix to remove 
+  // This only works if there is a lookup matrix available
+
+
+  if(opts.excl_mask.value()!=""){
+    volume<int> lookup_tract;
+    volume<int> excl;
+    read_volume(lookup_tract,"lookup_tractspace_"+ip);
+    string exname=opts.excl_mask.value();
+    make_basename(exname);
+    read_volume(excl,exname);
+    if(!samesize(excl,lookup_tract)){
+      cerr<<"Whoops - your exlusion mask does not appear to be in the same space as your original low resolution mask - sorry"<<endl;
+      return(-1);
+    }
+    vector<int> excl_cols;
+    for(int k=0;k<=excl.zsize();k++){
+      for(int j=0;j<=excl.ysize();j++){
+	for(int i=0;i<=excl.xsize();i++){
+	  if(excl(i,j,k)==1){
+	    if(lookup_tract(i,j,k)!=0){
+	     
+	      if(lookup_tract(i,j,k)<=mytractcoordmat.Nrows()){
+		excl_cols.push_back(lookup_tract(i,j,k)+1);
+	      }
+	      else{
+		cerr<<"Something a bit dodgy has happened here"<<endl;
+		cerr<<"Have you already run a reord_OM on this matrix"<<endl;
+		cerr<<"If so you can't use an exclusion mask as the"<<endl;
+		cerr<<"tractspace_lookup volume is not valid for this matrix"<<endl;
+	      }
+	     
+	     
+	    }
+	  }
+	}
+      }
+    }
+
+    rem_cols(myOMmat,mytractcoordmat,tractcoordbool,excl_cols);
+  }
+
+ 
+  rem_zrowcol(myOMmat,mycoordmat,mytractcoordmat,coordbool,tractcoordbool,newOMmat,newcoordmat,newtractcoordmat);
+  //   cerr<<"NOW"<<endl;
+  //   cerr<<myOMmat.MaximumAbsoluteValue()<<endl;
+  //   cerr<<newOMmat.MaximumAbsoluteValue()<<endl;
+ 
+  //write_ascii_matrix("ncm",newcoordmat);
+  // write_ascii_matrix("nctm",newtractcoordmat);
+
+
+
+  string base=opts.basename.value();
+  make_basename(base);
+  volume<float> outCCvol(newOMmat.Nrows(),newOMmat.Nrows(),1);
+  volume<int> outcoords(newcoordmat.Nrows(),3,1);
+
+
+
+  cerr<<"Computing correlation"<<endl;
+  SymmetricMatrix CtCt;
+  CtCt << corrcoef(newOMmat.t());
+  CtCt << CtCt+1;
+  if(opts.power.value()!=1){
+    CtCt << pow(CtCt,opts.power.value());
+  }
+ 
+  if(!opts.reord1.value()){
+     
+    for(int j=0;j<outCCvol.ysize();j++){
+      for(int i=0;i<outCCvol.xsize();i++){
+	outCCvol(i,j,0)=CtCt(i+1,j+1);
+      }
+    }
+   
+   
+    if(coordbool){
+      for(int i=0;i<outcoords.xsize();i++){
+	outcoords(i,0,0)=(int)newcoordmat(i+1,1);
+	outcoords(i,1,0)=(int)newcoordmat(i+1,2);
+	outcoords(i,2,0)=(int)newcoordmat(i+1,3);
+      } 
+    }
+   
+    save_volume(outCCvol,"CC_"+base);
+    save_volume(outcoords,"coords_for_"+base);
+
+  }
+  else{
+    cerr<<"Starting First Reordering"<<endl;
+    spect_reord(CtCt,r1,y1);
+   
+   
+    cerr<<"Permuting seed CC matrix"<<endl;
+    for(int j=0;j<outCCvol.ysize();j++){
+      for(int i=0;i<outCCvol.xsize();i++){
+	outCCvol(i,j,0)=CtCt((int)r1(i+1),(int)r1(j+1));
+      }
+    }
+   
+   
+    if(coordbool){
+      cerr<<"Permuting Seed Coordinates"<<endl;
+      for(int i=0;i<outcoords.xsize();i++){
+	outcoords(i,0,0)=(int)newcoordmat(int(r1(i+1)),1);
+	outcoords(i,1,0)=(int)newcoordmat(int(r1(i+1)),2);
+	outcoords(i,2,0)=(int)newcoordmat(int(r1(i+1)),3);
+      } 
+    }
+   
+    write_ascii_matrix(r1,base+"r1");
+    write_ascii_matrix(y1,base+"y1");
+    save_volume(outCCvol,"reord_CC_"+base);
+    save_volume(outcoords,"coords_for_reord_"+base);
+ 
+  }
+
+  if(opts.reord2.value()){
+    cerr<<"Starting Second Reordering"<<endl;
+    SymmetricMatrix CC;
+    CC << corrcoef(newOMmat);
+    CC<<CC+1;
+    spect_reord(CC,r2,y2);
+ 
+    write_ascii_matrix(r2,base+"r2");
+    write_ascii_matrix(y2,base+"y2");
+
+    volume<int> outvol(newOMmat.Nrows(),newOMmat.Ncols(),1);
+    volume<int> outtractcoords(newtractcoordmat.Nrows(),3,1);
+
+
+    cerr<<"Permuting Matrix"<<endl;
+    for(int j=0;j<outvol.ysize();j++){
+      for(int i=0;i<outvol.xsize();i++){
+	outvol(i,j,0)=(int)newOMmat((int)r1(i+1),(int)r2(j+1));
+      }
+    }
+
+    if(tractcoordbool){
+      cerr<<"Permuting Tract Coordinates"<<endl;
+      for(int i=0;i<outtractcoords.xsize();i++){
+	outtractcoords(i,0,0)=(int)newtractcoordmat(int(r2(i+1)),1);
+	outtractcoords(i,1,0)=(int)newtractcoordmat(int(r2(i+1)),2);
+	outtractcoords(i,2,0)=(int)newtractcoordmat(int(r2(i+1)),3);
+      } 
+    }
+    save_volume(outvol,"reord_"+base);
+    save_volume(outtractcoords,"tract_space_coords_for_reord_"+base);
+  }
+  return 0;
+}
+ 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

ccopsOptions.cc

+/*  ccopsOptions.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#define WANT_STREAM
+#define WANT_MATH
+
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "ccopsOptions.h"
+#include "utils/options.h"
+//#include "newmat.h"
+using namespace Utilities;
+
+namespace CCOPS {
+
+ccopsOptions* ccopsOptions::gopt = NULL;
+
+bool ccopsOptions::parse_command_line(int argc, char** argv)
+{
+
+  
+  for(int a = options.parse_command_line(argc, argv); a < argc; a++)
+    ;
+  if(help.value() || ! options.check_compulsory_arguments())
+    {
+      options.usage();
+      //throw NEWMAT::Exception("Not all of the compulsory arguments have been provided");
+      return false;
+    }      
+  return true;
+}
+
+}
+
+
+
+
+
+
+

ccopsOptions.h

+/*  BpmOptions.h
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(ccopsOptions_h)
+#define ccopsOptions_h
+
+#include <string>
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "utils/options.h"
+//#include "newmatall.h"
+using namespace Utilities;
+
+namespace CCOPS {
+
+class ccopsOptions {
+ public:
+  static ccopsOptions& getInstance();
+  ~ccopsOptions() { delete gopt; }
+  
+  Option<bool> help;
+  Option<string> inmatrix;
+  Option<string> basename;
+  Option<string> excl_mask;
+  Option<bool> reord1;
+  Option<bool> reord2;
+  Option<int> bin;
+  Option<float> power;
+  bool parse_command_line(int argc, char** argv);
+  
+ private:
+  ccopsOptions();  
+  const ccopsOptions& operator=(ccopsOptions&);
+  ccopsOptions(ccopsOptions&);
+
+  OptionParser options; 
+      
+  static ccopsOptions* gopt;
+  
+};
+
+ inline ccopsOptions& ccopsOptions::getInstance(){
+   if(gopt == NULL)
+     gopt = new ccopsOptions();
+   
+   return *gopt;
+ }
+
+ inline ccopsOptions::ccopsOptions() :
+   help(string("-h,--help"), false,
+	string("display this message"),
+	false, no_argument),
+   inmatrix(string("-i,--in"), string(""),
+	       string("input matrix"),
+	       true, requires_argument),  
+   basename(string("-b,--basename"), string(""),
+	       string("Output basename"),
+	       true, requires_argument),
+   excl_mask(string("-x"), string(""),
+	     string("exclusion mask"),
+	     false, requires_argument),  
+   reord1(string("--r1"), bool(false),
+	     string("do seedspace reordering (default no)"),
+	     false, no_argument), 
+   reord2(string("--r2"), bool(false),
+	     string("do tractspace reordering (default no)"),
+	     false, no_argument), 
+   bin(string("--bin"), 0, 
+	 string("binarise at (default 0 - no binarisation)"), 
+	 false, requires_argument),
+   power(string("-p,--power"), 1, 
+	 string("power to raise the correlation matrix to (default 1)"), 
+	 false, requires_argument),
+   options("ccops","")
+   {
+     
+    
+     try {
+       options.add(help);
+       options.add(inmatrix);
+       options.add(excl_mask);
+       options.add(reord1);
+       options.add(reord2);
+       options.add(bin);
+       options.add(power);
+       
+     }
+     catch(X_OptionError& e) {
+       options.usage();
+       cerr << endl << e.what() << endl;
+     } 
+     catch(std::exception &e) {
+       cerr << e.what() << endl;
+     }    
+     
+   }
+}
+
+#endif
+
+
+
+
+