*** empty log message ***

make_dyadic_vectors.cc

@@ -7,16 +7,28 @@ using namespace NEWIMAGE;
 
 int main ( int argc, char **argv ){
   if(argc<4){
-    cerr<<"usage: make_dyadic_vectors <theta_vol> <phi_vol> <output>"<<endl;
+    cerr<<"usage: make_dyadic_vectors <theta_vol> <phi_vol> [mask] <output>"<<endl;
+    cerr<<"[mask] is optional"<<endl;
     exit(0);
   }
+
   volume4D<float> ths,phs;
   volumeinfo tempinfo;
   read_volume4D(ths,argv[1],tempinfo);
   read_volume4D(phs,argv[2]);
+  volume<float> mask;
+  string oname;
+  if(argc==5){
+    oname=argv[4];
+    read_volume(mask,argv[3]);
+  }
+  else{
+    mask=ths[0]*0+1;
+    oname=argv[3];
+  }
   volume4D<float> dyadic_vecs(ths.xsize(),ths.ysize(),ths.zsize(),3);
   dyadic_vecs=0;
-  copybasicproperties(ths,dyadic_vecs);
+  copybasicproperties(ths[0],dyadic_vecs[0]);
   SymmetricMatrix dyad(3);dyad=0;
   ColumnVector dir(3);
   
@@ -26,35 +38,42 @@ int main ( int argc, char **argv ){
   for(int k=ths.minz();k<=ths.maxz();k++){
     for(int j=ths.miny();j<=ths.maxy();j++){
       for(int i=ths.minx();i<=ths.maxx();i++){
-	dyad=0;
-	for(int s=ths.mint();s<=ths.maxt();s++){		   
-	  float th=ths(i,j,k,s);
-	  float ph=phs(i,j,k,s);
-	  dir(1)=sin(th)*cos(ph);
-	  dir(2)=sin(th)*sin(ph);
-	  dir(3)=cos(th);
-	  dyad << dyad+dir*dir.t();
-	}
-	
-	EigenValues(dyad,dyad_D,dyad_V);
-	int maxeig;
-	if(dyad_D(1)>dyad_D(2)){
-	  if(dyad_D(1)>dyad_D(3)) maxeig=1;
-	  else maxeig=3;
+	if(mask(i,j,k)>0){
+	  dyad=0;
+	  for(int s=ths.mint();s<=ths.maxt();s++){		   
+	  
+	    float th=ths(i,j,k,s);
+	    float ph=phs(i,j,k,s);
+	    dir(1)=sin(th)*cos(ph);
+	    dir(2)=sin(th)*sin(ph);
+	    dir(3)=cos(th);
+	    dyad << dyad+dir*dir.t();
+	  }
+	  
+	  EigenValues(dyad,dyad_D,dyad_V);
+	  int maxeig;
+	  if(dyad_D(1)>dyad_D(2)){
+	    if(dyad_D(1)>dyad_D(3)) maxeig=1;
+	    else maxeig=3;
+	  }
+	  else{
+	    if(dyad_D(2)>dyad_D(3)) maxeig=2;
+	    else maxeig=3;
+	  }
+	  dyadic_vecs(i,j,k,0)=dyad_V(1,maxeig);
+	  dyadic_vecs(i,j,k,1)=dyad_V(2,maxeig);
+	  dyadic_vecs(i,j,k,2)=dyad_V(3,maxeig);
 	}
 	else{
-	  if(dyad_D(2)>dyad_D(3)) maxeig=2;
-	  else maxeig=3;
+	  dyadic_vecs(i,j,k,0)=0;
+	  dyadic_vecs(i,j,k,1)=0;
+	  dyadic_vecs(i,j,k,2)=0;
 	}
-	dyadic_vecs(i,j,k,0)=dyad_V(1,maxeig);
-	dyadic_vecs(i,j,k,1)=dyad_V(2,maxeig);
-	dyadic_vecs(i,j,k,2)=dyad_V(3,maxeig);
       }
-      
     }
   }
   
-  save_volume4D(dyadic_vecs,argv[3],tempinfo);
+  save_volume4D(dyadic_vecs,oname,tempinfo);
 }
 
 

old_bedpost

@@ -197,6 +197,8 @@ ${FSLDIR}/bin/avwmaths ${subjdir}.bedpost/merged_thsamples -Tmean ${subjdir}.bed
 ${FSLDIR}/bin/avwmaths ${subjdir}.bedpost/merged_phsamples -Tmean ${subjdir}.bedpost/mean_phsamples
 ${FSLDIR}/bin/avwmaths ${subjdir}.bedpost/merged_fsamples -Tmean ${subjdir}.bedpost/mean_fsamples
 
+${FSLDIR}/bin/make_dyadic_vectors ${subjdir}.bedpost/merged_thsamples ${subjdir}.bedpost/merged_phsamples ${subjdir}.bedpost/nodif_brain_mask ${subjdir}.bedpost/dyadic_vectors
+
 if [ ${finished} -eq 1 ];then
     if [ `imtest ${subjdir}.bedpost/merged_thsamples` -eq 1 ];then
 	if [ `imtest ${subjdir}.bedpost/merged_phsamples` -eq 1 ];then