reordering dyads - first sketch

reorder_dyadic_vectors.cc

+/*  reorder_dyadic_vectors.cc
+
+    Saad Jbabdi, FMRIB Image Analysis Group
+
+    Copyright (C) 2006 University of Oxford */
+
+/*  CCOPYRIGHT */
+
+#include "utils/options.h"
+#include "heap.h"
+#include "newimage/newimageall.h"
+#include "miscmaths/miscprob.h"
+#include "stdlib.h"
+
+#define AP 2
+#define NB 1
+#define FA 0
+
+
+using namespace Utilities;
+using namespace std;
+using namespace NEWIMAGE;
+using namespace MISCMATHS;
+
+
+string title="reorder_dyadic_vectors (Version 1.0)\nReordering of the dyadic vectors with direction preservation";
+string examples="reorder_dyadic_vectors -b <dirname> [-m mask]\n<dirname> is supposed to contain mean_fisamples and dyadsi for i=1:n";
+
+Option<bool> verbose(string("-v,--verbose"),false,
+		       string("switch on diagnostic messages"),
+		       false,no_argument);
+Option<bool> help(string("-h,--help"),false,
+		       string("display this message"),
+		       false,no_argument);
+Option<string> basename(string("-b,--basename"),string(""),
+		       string("output basename"),
+		       true,requires_argument);
+Option<string> mask(string("-m,--mask"),string(""),
+		       string("filename of brain mask"),
+		       false,requires_argument);
+////////////////////////////////////////////////////////
+
+class FM{
+ private:
+  volume<float> mask,tmap,state;
+  vector< volume4D<float> > dyads;
+  vector< volume<float> >  meanf;
+  int nx,ny,nz;
+  Matrix P;
+  Matrix neighbours;
+
+  Option<string>& obasename;
+  Option<string>& omask;
+
+ public:
+  FM(Option<string>& optbasename,Option<string>& optmask):obasename(optbasename),omask(optmask){
+
+    int fib=1;
+    bool fib_existed=true;
+    while(fib_existed){
+      if(fsl_imageexists(obasename.value()+"mean_f"+num2str(fib)+"samples")){
+	volume4D<float> *tmpdptr=new volume4D<float>;
+	volume<float> *tmpfptr=new volume<float>;
+	read_volume4D(*tmpdptr,obasename.value()+"dyads"+num2str(fib));
+	dyads.push_back(*tmpdptr);
+	read_volume(*tmpfptr,obasename.value()+"mean_f"+num2str(fib)+"samples");
+	meanf.push_back(*tmpfptr);
+	fib++;
+      }
+      else
+	fib_existed=false;
+    }
+
+    // read input volumes
+    nx=dyads[0].xsize();
+    ny=dyads[0].ysize();
+    nz=dyads[0].zsize();
+    // create other volumes
+    if(omask.set()){
+      volume<float> mask;
+      read_volume(mask,omask.value());
+    }
+    else{
+      volume<float> maskvol(nx,ny,nz);
+      copybasicproperties(dyads[0],mask);
+      mask=1;
+    }
+
+    state.reinitialize(nx,ny,nz);
+    tmap.reinitialize(nx,ny,nz);
+
+    P=MISCMATHS::perms(dyads.size());
+
+    neighbours.ReSize(3,26);
+    neighbours << 1 << 0 << 0 << -1 << 0 << 0 << 1 << 1 <<-1 <<-1 << 1 <<-1 << 1 <<-1 << 0 << 0 << 0 << 0 << 1 <<-1 << 1 << 1 <<-1 <<-1 << 1 <<-1 
+	       << 0 << 1 << 0 << 0  <<-1 << 0 << 1 <<-1 << 1 <<-1 << 0 << 0 << 0 << 0 << 1 <<-1 << 1 <<-1 << 1 << 1 <<-1 << 1 <<-1 << 1 <<-1 <<-1
+	       << 0 << 0 << 1 << 0  << 0 <<-1 << 0 << 0 << 0 << 0 << 1 << 1 <<-1 <<-1 << 1 << 1 <<-1 <<-1 << 1 << 1 << 1 <<-1 << 1 <<-1 <<-1 <<-1;
+  }
+  
+  ////////////////////// functions
+  void do_fast_marching(){
+    int i0,j0,k0;                      /* current point coords */
+    float tval; 
+    
+    /********************************/
+    /*** performing fmt algorithm ***/
+    /********************************/
+    
+    /*** initialization ***/
+    
+    /* look for bigger f1+f2 point as a seed */
+    float maxf=0,curf;
+    for(int z=0;z<nz;z++)
+      for(int y=0;y<ny;y++)
+	for(int x=0;x<nx;x++){
+	  curf=0;
+	  for(unsigned int f=0;f<meanf.size();f++)
+	    curf += meanf[f](x,y,z);
+	  if(curf>maxf){
+	    i0=x;j0=y;k0=z;
+	    maxf=curf;
+	  }
+	  state(x,y,z)=FA;
+	  tmap(x,y,z)=1;
+	}
+    state(i0,j0,k0)=AP;
+    tmap(i0,j0,k0)=0;
+    
+     /* create heap sort structure */
+    Heap h(nx,ny,nz);
+    
+    /* and all points of the ROIs as Alive Points */
+    updateNBvalue(i0,j0,k0,h);
+    /*--------------------------------------------------------*/
+    /*** big loop ***/
+    int STOP = 0;
+    int counter = 0;
+    while(STOP==0){
+      /*break;*/
+      counter++;
+      
+      /*** uses the heap sort structure to find the NB point with the smallest T-value ***/
+      h.heapRemove(tval,i0,j0,k0);
+      
+      /*** add this point to the set of alive points ***/
+      state(i0,j0,k0)=AP;
+      
+      if(h.get_N()==0)
+	break;
+      
+      /*** update narrow band's T-value ***/
+      updateNBvalue(i0,j0,k0,h);
+      
+    }
+    
+  }
+  int isInside(int i,int j,int k){
+    return((i>=0) && (i<nx)
+	   && (j>=0) && (j<ny) 
+	   && (k>=0) && (k<nz) 
+	   && (mask(i,j,k)!=0));
+  }
+  void computeT(int x,int y,int z){
+    Matrix V(P.Ncols(),3);
+    Matrix nV(P.Ncols(),3);
+    ColumnVector mF(P.Ncols());
+
+    V=get_vector(x,y,z);
+    mF=get_f(x,y,z);
+
+    int nbx,nby,nbz;
+    int opt_perm=1;
+    float opt_sperm=0;
+    for(int per=1;per<P.Nrows();per++){
+      
+      float opt_nb=0;int nnb=0;
+      for(int nb=1;nb<=neighbours.Ncols();nb++){
+	nbx=x+neighbours(1,nb);nby=y+neighbours(2,nb);nbz=z+neighbours(3,nb);
+	if(!isInside(nbx,nby,nbz))
+	  continue;
+	if(state(nbx,nby,nbz)==AP){
+	  nV=get_vector(nbx,nby,nbz);
+	  
+	  float opt_s=0,s;
+	  for(int f=1;f<=V.Nrows();f++){
+	    s=dot(nV.Row(f).t(),V.Row(P(per,f)).t());
+	    if(s>opt_s){
+	      opt_s=s;
+	    }
+	  }//f
+	  opt_nb+=opt_s;
+	  nnb++;
+	}//nb
+	opt_nb/=nnb;
+	
+	if(opt_nb>opt_sperm){
+	  opt_sperm=opt_nb;
+	  opt_perm=per;
+	}
+	
+      }//endif
+    }//perm
+    tmap(x,y,z)=1-opt_sperm; // store optimal mean scalar product
+
+    // reorder dyads
+    for(unsigned int f=0;f<dyads.size();f++){
+      dyads[f](x,y,z,0)=V(P(opt_perm,f+1),1);
+      dyads[f](x,y,z,1)=V(P(opt_perm,f+1),2);
+      dyads[f](x,y,z,2)=V(P(opt_perm,f+1),3);
+
+      meanf[f](x,y,z)=mF(P(opt_perm,f+1));
+    }
+  }
+
+  ReturnMatrix get_vector(int x,int y,int z){
+    Matrix V(P.Ncols(),3);
+    for(unsigned int f=0;f<dyads.size();f++)
+      for(int i=1;i<=3;i++)
+	V(f+1,i)=dyads[f](x,y,z,i-1);
+
+    V.Release();
+    return V;
+  }
+
+  ReturnMatrix get_f(int x,int y,int z){
+    ColumnVector V(P.Ncols());
+    for(unsigned int f=0;f<dyads.size();f++)
+      V(f+1)=meanf[f](x,y,z);
+
+    V.Release();
+    return V;
+  }
+
+  void updateNBvalue(int i,int j,int k,Heap& h){    
+    int ni,nj,nk;
+    int pos;
+    double val;
+   
+    for(int nb=1;nb<neighbours.Ncols();nb++){
+      ni=i+neighbours(1,nb);nj=j+neighbours(2,nb);nk=k+neighbours(3,nb);
+      if(isInside(ni,nj,nk)){
+	if(state(ni,nj,nk)==AP)continue;
+	computeT(ni,nj,nk);
+	val=tmap(ni,nj,nk);
+	/* update value if in NB */
+	if(state(ni,nj,nk)==NB){
+	  pos=h.get_bpr(ni,nj,nk);
+	  h.set_val(pos,MIN(h.get_val(pos),val));
+	  h.heapUp(pos);
+	}
+	/* insert value if in FA */
+	else{
+	  state(ni,nj,nk)=NB;
+	  h.heapInsert(val,ni,nj,nk);
+	}
+      }
+      
+    }
+  }
+  
+  void save_results(){
+    int fib=1;
+    for(unsigned int f=0;f<dyads.size();f++){
+      save_volume4D(dyads[f],obasename.value()+"reordered_dyads"+num2str(fib));
+      save_volume(meanf[f],obasename.value()+"reordered_meanf"+num2str(fib)+"samples");
+      fib++;
+    }
+  }
+
+};
+
+
+
+
+int main(int argc,char *argv[]){
+
+  Tracer tr("main");
+  OptionParser options(title,examples);
+
+  try{
+    options.add(verbose);
+    options.add(help);
+    options.add(basename);
+    options.add(mask);
+
+    options.parse_command_line(argc,argv);
+
+    
+    if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+      options.usage();
+      exit(EXIT_FAILURE);
+    }
+
+
+    FM fm(basename,mask);
+  
+    fm.do_fast_marching();
+
+    fm.save_results();
+
+
+
+  }
+  catch(X_OptionError& e) {
+    options.usage();
+    cerr << endl << e.what() << endl;
+    exit(EXIT_FAILURE);
+  } 
+  catch(std::exception &e) {
+    cerr << e.what() << endl;
+  } 
+
+  return 0;
+  
+  
+}