From c0c3bcefe5edb1db84deabe92f06bc58a45e2327 Mon Sep 17 00:00:00 2001 From: Tim Behrens <behrens@fmrib.ox.ac.uk> Date: Tue, 8 Nov 2005 16:17:04 +0000 Subject: [PATCH] *** empty log message *** --- Makefile | 7 +- ccops.cc | 513 ++++++++++++++++++++++++++++++++++++++++++++++++ ccopsOptions.cc | 47 +++++ ccopsOptions.h | 111 +++++++++++ 4 files changed, 677 insertions(+), 1 deletion(-) create mode 100644 ccops.cc create mode 100644 ccopsOptions.cc create mode 100644 ccopsOptions.h diff --git a/Makefile b/Makefile index c1e4932..8489c03 100644 --- a/Makefile +++ b/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} diff --git a/ccops.cc b/ccops.cc new file mode 100644 index 0000000..1985a3a --- /dev/null +++ b/ccops.cc @@ -0,0 +1,513 @@ +/* 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; +} + + + + + + + + + + + + + + + + + + + diff --git a/ccopsOptions.cc b/ccopsOptions.cc new file mode 100644 index 0000000..5a8272d --- /dev/null +++ b/ccopsOptions.cc @@ -0,0 +1,47 @@ +/* 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; +} + +} + + + + + + + diff --git a/ccopsOptions.h b/ccopsOptions.h new file mode 100644 index 0000000..a86b77c --- /dev/null +++ b/ccopsOptions.h @@ -0,0 +1,111 @@ +/* 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 + + + + + -- GitLab