Implemented efficient sort

cluster.cc

@@ -9,17 +9,18 @@
 #include "fmribmain.h"
 #include "newimageall.h"
 #include <vector>
+#include <algorithm>
 
 using namespace NEWIMAGE;
 
 void print_usage(char *argv[])
 {
   cerr << "Usage: " << argv[0] << " <input file> <output file> "  
-       << "<threshold> <output_format> [-fe]" << endl;
+       << "<threshold> <output_format> [options]" << endl;
   cerr << "[-f] : threshold is fractional (of robust range) instead of "
        << "absolute"<<endl;
-  cerr << "[-e] : used only edge connections (6-way) instead of "
-       << "corner connections (26-way)"<<endl;
+  cerr << "[-6,18,26] : number of connections at each point (default is 26)"
+       <<endl;
   cerr << "If threshold is -ve then binarization is inverted"<<endl;
   cerr << "Output formats: the output intensity signifies:"<<endl;
   cerr << "0 - size"<<endl;
@@ -67,21 +68,19 @@ void get_stats(const volume<int>& labelim, const volume<T>& origim,
 
 void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder) 
 {
-  // A BAD O(N^2) SORTING ALGORITHM...
-  int labelnum=size.size()-1, maxidx=0, maxsize=0;
-  sizeorder.resize(labelnum,0);
-  std::vector<int> sizecopy(size);
-  for (int n=labelnum; n>=1; n--) {
-    maxidx=0;
-    maxsize=0;
-    for (int m=1; m<=labelnum; m++) {
-      if (sizecopy[m]>maxsize) {
-	maxidx = m;
-	maxsize = sizecopy[m];
-      }
-    }
-    sizecopy[maxidx]=0;
-    sizeorder[maxidx] = n;
+  // find order of labels so that sizes are in *ascending* order
+  int length=size.size()-1;
+  std::vector<pair<int, int> > sortlist(length+1);
+  for (int n=0; n<=length; n++) {
+    sortlist[n] = pair<int, int>(size[n],n);
+  }
+  sortlist[0].first = 0;  // force the background component to be in 0th posn
+  sort(sortlist.begin(),sortlist.end());  // O(N.log(N))
+
+  // second part of pair is now the prior-index of the sorted values
+  sizeorder.resize(length+1,0);
+  for (int n=0; n<=length; n++) {
+    sizeorder[sortlist[n].second] = n;    // maps old index to new
   }
 }
 
@@ -113,29 +112,32 @@ int fmrib_main(int argc, char *argv[])
     volume<T> vin;
     read_volume(vin,argv[1]);
     
-    // threshold the volume
-    bool invert=false;
-    float th;
-    string opt1="", opt2="", opt="";
-    if (argc>=6) { opt = argv[5]; }
-    if (opt=="-f") { opt1="-f"; }
-    else if (opt=="-e") { opt2="-e"; }
-    else if ( (opt=="-fe") || (opt=="-ef") ) { opt1="-f"; opt2="-e"; }
-    else if (opt!="") { 
-      cerr << "Unrecognised option "<<opt<<endl<<endl;
-      print_usage(argv); 
-      return -1;
+    // parse the options
+    bool invert=false, fractionalthresh=false;
+    int numconnected=26;
+
+    for (int n=5; n<argc; n++) {
+      string opt=argv[n];
+      if (opt=="-f") { fractionalthresh=true; }
+      else if (opt=="-6") { numconnected=6; }
+      else if (opt=="-18") { numconnected=18; }
+      else if (opt=="-26") { numconnected=26; }
+      else { 
+	cerr << "Unrecognised option "<<opt<<endl<<endl;
+	print_usage(argv); 
+	return -1;
+      }
     }
 
-    bool use_corners=true;
-    if (opt2=="-e") { use_corners=false; }
-
+    // threshold the volume
+    float th;
     th = atof(argv[3]);
     if (th<0.0) { 
       invert=true;
       th=-th;
     }
-    if ( opt1=="-f" ) {
+
+    if ( fractionalthresh ) {
       float frac = th;
       th = frac*(vin.robustmax() - vin.robustmin()) + vin.robustmin();
     }
@@ -143,7 +145,7 @@ int fmrib_main(int argc, char *argv[])
     if (invert) { vin = ((T) 1) - vin; }
     
     // Find the connected components
-    labelim = connected_components(vin,use_corners);
+    labelim = connected_components(vin,numconnected);
     
     // process according to the output format
     read_volume(vin,argv[1]);  // re-read the original volume