Sort a dataset according to the angular distance of each DW direction with a reference vector

rearrange.cc

+/*  Copyright (C) 2010 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <string>
+#include <cmath>
+#include "newimage/newimageall.h"
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace NEWMAT;
+
+//Rearranges a dataset according to the angular distance between the corresponding bvecs entry and a reference vector  
+int main ( int argc, char *argv[]){
+  if(argc<5 || argc>6){
+    cerr<<" "<<endl;
+    cerr<<"usage: rearrange data bvecs ref_dyads [brain_mask] ouput"<<endl;
+    cerr<<" "<<endl;
+    exit(1);
+  }
+  
+  ColumnVector ref_vector(3);
+  Matrix bvecs;   //Read Input 
+  volume4D<float> data, ref; volume<float> mask;
+  
+  read_volume4D(data,argv[1]);
+  bvecs=read_ascii_matrix(argv[2]);
+  read_volume4D(ref,argv[3]);
+  if (argc==6)
+    read_volume(mask,argv[4]);
+  else{ 
+    mask.reinitialize(data.xsize(),data.ysize(),data.zsize());
+    mask=1;
+  }
+  
+  volume4D<float> output(data.xsize(),data.ysize(),data.zsize(),data.tsize());
+  copybasicproperties(data,output);output=0;
+
+  if(bvecs.Nrows()>3) bvecs=bvecs.t();   //Make sure the bvecs entries are normalized unit vectors
+  for(int i=1;i<=bvecs.Ncols();i++){
+    float tmpsum=sqrt(bvecs(1,i)*bvecs(1,i)+bvecs(2,i)*bvecs(2,i)+bvecs(3,i)*bvecs(3,i));
+    if(tmpsum!=0){
+      bvecs(1,i)=bvecs(1,i)/tmpsum;
+      bvecs(2,i)=bvecs(2,i)/tmpsum;
+      bvecs(3,i)=bvecs(3,i)/tmpsum;
+    }  
+  }
+  for(int z=data.minz();z<=data.maxz();z++){   //Rearrange data for each voxel
+    for(int y=data.miny();y<=data.maxy();y++){
+      for(int x=data.minx();x<=data.maxx();x++){
+	if (mask(x,y,z)!=0){
+	  vector<pair<float,int> > dot;  //Keep in the first entry of each pair the dot product value and the in the second the index 
+	  pair<float,int> ftmp;
+	  ref_vector(1)=ref(x,y,z,0); 	  ref_vector(2)=ref(x,y,z,1); 	  ref_vector(3)=ref(x,y,z,2);
+	  for (int n=1; n<=bvecs.Ncols(); n++){     //Get the dot product of each bvecs entry with the reference orientation in this voxel
+	    if (bvecs(1,n)==0 && bvecs(2,n)==0 && bvecs(3,n)==0)
+	      ftmp.first=2.0;              //The b=0 entries are first in the sequence
+	    else
+	      ftmp.first=fabs(ref_vector(1)*bvecs(1,n)+ref_vector(2)*bvecs(2,n)+ref_vector(3)*bvecs(3,n));
+	    ftmp.second=n-1;
+	    dot.push_back(ftmp);
+	    } 
+	  sort(dot.begin(),dot.end());     //Sort the dot products in ascending order
+	  reverse(dot.begin(),dot.end());  //Reverse the ordering, so that it is in descending order
+	  for (int n=0; n<bvecs.Ncols(); n++)//Get the dot product of each bvecs entry with the reference orientation in this voxel
+	    output(x,y,z,n)=data(x,y,z,dot[n].second);
+	}
+      }
+    }
+    cout<<z+1<<" slices processed"<<endl;
+  }
+  save_volume4D(output,argv[argc-1]);
+  return 0;
+}
+
+
+
+
+
+
+
+
+