Added -rank and -ranknorm options

fslmaths.cc

@@ -5,7 +5,7 @@
 //     
 //     CCOPYRIGHT  
 //     
-
+#define EXPOSE_TREACHEROUS
 #include "newimage/newimageall.h"
 #include "miscmaths/miscmaths.h"
 #include "utils/fsl_isfinite.h"
@@ -109,6 +109,8 @@ int printUsage(const string& programName)
   cout << " -cpval   : Same as -pval, but gives FWE corrected P-values" << endl;
   cout << " -ztop    : Convert Z-stat to (uncorrected) P" << endl;
   cout << " -ptoz    : Convert (uncorrected) P to Z" << endl;
+  cout << " -rank    : Convert data to ranks (over T dim)" << endl;
+  cout << " -ranknorm: Transform to Normal dist via ranks" << endl;
 
   cout << "\nMulti-argument operations:" << endl;
   cout << " -roi <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> : zero outside roi (using voxel coordinates)" << endl;
@@ -729,6 +731,70 @@ if (!separatenoise)
               input_volume.value(x,y,z,t)=(T)((v<=0 || v>=1)?0:ndtri(1-v));
             }
     }
+    /***** Convert to ranks ********/
+    else if (string(argv[i])=="-rank")
+    {
+      ColumnVector val(input_volume.tsize()),rank(input_volume.tsize()),sortval(input_volume.tsize());
+      for(int z=0;z<input_volume.zsize();z++)
+	for(int y=0;y<input_volume.ysize();y++)           
+	  for(int x=0;x<input_volume.xsize();x++) {
+	    for(int t=0;t<input_volume.tsize();t++)
+	      val(t+1)=input_volume.value(x,y,z,t);
+	    
+	    
+	    /* Take sortval and 'unsort' it, finding ranks */
+	    sortval=val;
+	    SortAscending(sortval);
+	    for(int k=1;k<=val.Nrows();k++)	       
+	      rank(k)=k;
+	    for(int k=1;k<=val.Nrows();k++)	       
+	      if(val(k)!=sortval(k))
+		for(int l=k+1;l<=val.Nrows();l++) 
+		  if(sortval(l)==val(k))
+		    {
+		      swap(rank(l),rank(k));
+		      swap(sortval(l),sortval(k));
+		    }
+	    for(int t=0;t<input_volume.tsize();t++)
+	      input_volume.value(x,y,z,t)=(T)rank(t+1);
+	  }
+
+    }
+    /***** Convert to normal distribution via ranks ********/
+    else if (string(argv[i])=="-ranknorm")
+    {
+      ColumnVector val(input_volume.tsize()),rank(input_volume.tsize()),sortval(input_volume.tsize());
+      volume<double> valv(input_volume.tsize(),1,1);
+      for(int z=0;z<input_volume.zsize();z++)
+	for(int y=0;y<input_volume.ysize();y++)           
+	  for(int x=0;x<input_volume.xsize();x++) {
+	    for(int t=0;t<input_volume.tsize();t++) {
+	      val(t+1)=input_volume.value(x,y,z,t); 
+	      valv(t,0,0)=val(t+1);
+	    }
+	    
+	    /* Take sortval and 'unsort' it, finding ranks */
+	    sortval=val;
+	    SortAscending(sortval);
+	    for(int k=1;k<=val.Nrows();k++)	       
+	      rank(k)=k;
+	    for(int k=1;k<=val.Nrows();k++)	       
+	      if(val(k)!=sortval(k))
+		for(int l=k+1;l<=val.Nrows();l++) 
+		  if(sortval(l)==val(k))
+		    {
+		      swap(rank(l),rank(k));
+		      swap(sortval(l),sortval(k));
+		    }
+	    /* Transform to expected order statisics of a Uniform */
+	    rank = (rank-0.5)/rank.Nrows();
+	    /* Transform to expected order statisics of a Normal with same Mean & Sd */
+	    for(int t=0;t<input_volume.tsize();t++)
+	      input_volume.value(x,y,z,t)=(T)(valv.mean()+ndtri(rank(t+1))*valv.stddev()); 
+	    
+	  }
+
+    }
     /***************************************************************/
     else if (string(argv[i])=="-abs") input_volume=abs(input_volume);
     /***************************************************************/