From 67ee073272c771211eecb0565fae8f9fe9316659 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Tue, 4 May 2004 14:20:41 +0000
Subject: [PATCH] *** empty log message ***
---
Makefile | 13 +++-
reord_OM.cc | 191 +++++++++++++++++++++++++++++++++++++++++++---------
2 files changed, 170 insertions(+), 34 deletions(-)
diff --git a/Makefile b/Makefile
index 0f67ade..4bd0b77 100644
--- a/Makefile
+++ b/Makefile
@@ -14,18 +14,22 @@ PT=probtrack
FTB=find_the_biggest
PJ=proj_thresh
MED=medianfilter
+ROM=reord_OM
+SAUS=sausages
DTIFITOBJS=dtifit.o dtifitOptions.o
PTOBJS=probtrack.o probtrackOptions.o
FTBOBJS=find_the_biggest.o
PJOBJS=proj_thresh.o
MEDOBJS=medianfilter.o
+ROMOBJS=reord_OM.o
+SAUSOBJS=sausages.o
SCRIPTS = eddy_correct
XFILES = dtifit probtrack find_the_biggest medianfilter
RUNTCLS = Fdt
-all: ${XFILES}
+all: ${XFILES} reord_OM sausages
${PT}: ${PTOBJS}
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${PTOBJS} ${DLIBS}
@@ -42,6 +46,13 @@ ${MED}: ${MEDOBJS}
${DTIFIT}: ${DTIFITOBJS}
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${DTIFITOBJS} ${DLIBS}
+${ROM}: ${ROMOBJS}
+ ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${ROMOBJS} ${DLIBS}
+
+${SAUS}: ${SAUSOBJS}
+ ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${SAUSOBJS} ${DLIBS}
+
+
diff --git a/reord_OM.cc b/reord_OM.cc
index ef3418e..35d5c00 100644
--- a/reord_OM.cc
+++ b/reord_OM.cc
@@ -182,18 +182,104 @@ void rem_zrowcol(const Matrix& myOMmat,const Matrix& coordmat,const Matrix& trac
}
}
-
}
+
+
+
+
+
+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){
+ cerr<<"Updating Tract Coordinates"<<endl;
+ 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 ){
if(argc<4){
- cerr<<"usage: reord_OM <matrix> <output base> <bin>"<<endl;
+ cerr<<"usage: reord_OM <matrix> <output base> <bin> [excl_mask]"<<endl;
cerr<<" matrix in 2D analyze format"<<endl;
cerr<<" bin=0 gives no binarising at all"<<endl;
cerr<<" any other <bin> binarises above that number "<<endl;
+ cerr<<""<<endl;
+ cerr<<"Optional"<<endl;
+ cerr<<" excl_mask - exclude contribution from "<<endl;
+ cerr<<" voxels in this mask - this needs to be in the same low res"<<endl;
+ cerr<<" space as the initial tract space analysis"<<endl;
+
exit(0);
}
string ip=argv[1];
@@ -268,18 +354,53 @@ int main ( int argc, char **argv ){
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(argc==5){
+ volume<int> lookup_tract;
+ volume<int> excl;
+ read_volume(lookup_tract,"lookup_tractspace_"+ip);
+ string exname=argv[4];
+ 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){
+ excl_cols.push_back(lookup_tract(i,j,k)+1);
+ }
+ }
+ }
+ }
+
+ 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);
-
+//write_ascii_matrix("ncm",newcoordmat);
+// write_ascii_matrix("nctm",newtractcoordmat);
+ string base=argv[2];
+ make_basename(base);
+ volume<float> outCCvol(newOMmat.Nrows(),newOMmat.Nrows(),1);
+ volume<int> outcoords(newcoordmat.Nrows(),3,1);
cerr<<"Starting first Reordering"<<endl;
SymmetricMatrix CtCt;
@@ -287,24 +408,46 @@ cerr<<"Starting first Reordering"<<endl;
CtCt << CtCt+1;
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);
+
+
+
+
cerr<<"Starting Second Reordering"<<endl;
SymmetricMatrix CC;
CC << corrcoef(newOMmat);
CC<<CC+1;
spect_reord(CC,r2,y2);
-
- string base=argv[2];
- make_basename(base);
- write_ascii_matrix(r1,base+"r1");
- write_ascii_matrix(y1,base+"y1");
write_ascii_matrix(r2,base+"r2");
write_ascii_matrix(y2,base+"y2");
+
volume<int> outvol(newOMmat.Nrows(),newOMmat.Ncols(),1);
- volume<float> outCCvol(newOMmat.Nrows(),newOMmat.Nrows(),1);
- volume<int> outcoords(newcoordmat.Nrows(),3,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++){
@@ -312,33 +455,15 @@ cerr<<"Starting first Reordering"<<endl;
}
}
- cerr<<"Seed CC"<<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(r1(i+1),1);
- outcoords(i,1,0)=(int)newcoordmat(r1(i+1),2);
- outcoords(i,2,0)=(int)newcoordmat(r1(i+1),3);
- }
- }
if(tractcoordbool){
cerr<<"Permuting Tract Coordinates"<<endl;
for(int i=0;i<outtractcoords.xsize();i++){
- outtractcoords(i,0,0)=(int)newtractcoordmat(r2(i+1),1);
- outtractcoords(i,1,0)=(int)newtractcoordmat(r2(i+1),2);
- outtractcoords(i,2,0)=(int)newtractcoordmat(r2(i+1),3);
+ 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(outCCvol,"reord_CC_"+base);
- save_volume(outcoords,"coords_for_reord_"+base);
save_volume(outtractcoords,"tract_space_coords_for_reord_"+base);
return 0;
--
GitLab