*** empty log message ***

cluster.cc

@@ -34,28 +34,28 @@ void print_usage(char *argv[])
   //  cerr << "6 - unprocessed unique integer label"<<endl;
 }
 
-//struct tripleint { int x; int y; int z; };
-//struct tripleint { float x; float y; float z; };
+struct triplei { int x; int y; int z; };
+struct triplef { float x; float y; float z; };
+typedef struct triplei tripleint;
+typedef struct triplef triplefloat;
 
 template <class T>
 void get_stats(const volume<int>& labelim, const volume<T>& origim,
 	       std::vector<int>& size,
 	       std::vector<T>& maxvals, std::vector<T>& meanvals,
-	       std::vector<int>& maxx, std::vector<int>& maxy, 
-	       std::vector<int>& maxz, std::vector<float>& cogx,
-	       std::vector<float>& cogy, std::vector<float>& cogz,
+	       std::vector<tripleint>& max, std::vector<triplefloat>& cog,
 	       bool invert) 
 {
   int labelnum = labelim.max();
   size.resize(labelnum+1,0);
   maxvals.resize(labelnum+1, (T) 0);
   meanvals.resize(labelnum+1,(T) 0);
-  maxx.resize(labelnum+1,0);
-  maxy.resize(labelnum+1,0);
-  maxz.resize(labelnum+1,0);
-  cogx.resize(labelnum+1,0.0f);
-  cogy.resize(labelnum+1,0.0f);
-  cogz.resize(labelnum+1,0.0f);
+  tripleint zero;
+  zero.x = 0; zero.y = 0; zero.z = 0;
+  triplefloat zerof;
+  zerof.x = 0; zerof.y = 0; zerof.z = 0;
+  max.resize(labelnum+1,zero);
+  cog.resize(labelnum+1,zerof);
   std::vector<float> sum(labelnum+1,0.0);
   for (int z=labelim.minz(); z<=labelim.maxz(); z++) {
     for (int y=labelim.miny(); y<=labelim.maxy(); y++) {
@@ -63,34 +63,30 @@ void get_stats(const volume<int>& labelim, const volume<T>& origim,
 	int idx = labelim(x,y,z);
 	T oxyz = origim(x,y,z);
 	size[idx]++;
-	cogx[idx]+=((float) oxyz)*x;
-	cogy[idx]+=((float) oxyz)*y;
-	cogz[idx]+=((float) oxyz)*z;
+	cog[idx].x+=((float) oxyz)*x;
+	cog[idx].y+=((float) oxyz)*y;
+	cog[idx].z+=((float) oxyz)*z;
 	sum[idx]+=(float) oxyz;
-	if (invert) {
-	  if ((size[idx]==1) || (maxvals[idx]>oxyz)) {
+	if ( (size[idx]==1) || 
+	     ( (oxyz>maxvals[idx]) && (!invert) ) || 
+	     ( (oxyz<maxvals[idx]) && (invert) ) ) 
+	  {
 	    maxvals[idx] = oxyz;
-	    maxx[idx] = x;
-	    maxy[idx] = y;
-	    maxz[idx] = z;
+	    max[idx].x = x;
+	    max[idx].y = y;
+	    max[idx].z = z;
 	  }
-	} else {
-	  if ((size[idx]==1) || (maxvals[idx]<oxyz)) {
-	    maxvals[idx] = oxyz;
-	    maxx[idx] = x;
-	    maxy[idx] = y;
-	    maxz[idx] = z;
-	  }
-	}
       }
     }
   }
   for (int n=0; n<=labelnum; n++) {
     if (size[n]>0.0) {
       meanvals[n] = (T) (sum[n]/((float) size[n]));
-      cogx[n] /= sum[n];
-      cogy[n] /= sum[n];
-      cogz[n] /= sum[n];
+    }
+    if (sum[n]>0.0) {
+      cog[n].x /= sum[n];
+      cog[n].y /= sum[n];
+      cog[n].z /= sum[n];
     }
   }
 }
@@ -117,7 +113,6 @@ std::vector<int> get_sortindex(const std::vector<int>& vals)
 void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder) 
 {
   std::vector<int> sizecopy(size), idx;
-  sizecopy[0] = 0;    // force the background component to be in 0th posn
   idx = get_sortindex(sizecopy);
 
   // second part of pair is now the prior-index of the sorted values
@@ -152,8 +147,9 @@ int fmrib_main(int argc, char *argv[])
   int outform = atoi(argv[6]);
   float dLh = atof(argv[7]);
   int V = atoi(argv[8]);
-  std::vector<int> size, maxx, maxy, maxz;
-  std::vector<float> cogx, cogy, cogz;
+  std::vector<int> size;
+  std::vector<tripleint> max;
+  std::vector<triplefloat> cog;
   std::vector<T> maxvals, meanvals;
   float th;
   bool invert=false, fractionalthresh=false;
@@ -196,8 +192,7 @@ int fmrib_main(int argc, char *argv[])
     
     // process according to the output format
     read_volume(vin,argv[1]);  // re-read the original volume
-    get_stats(labelim,vin,size,maxvals,meanvals,
-	      maxx,maxy,maxz,cogx,cogy,cogz,invert);
+    get_stats(labelim,vin,size,maxvals,meanvals,max,cog,invert);
   }
 
   cout << "CLUSTERS " << size.size()-1 << endl;
@@ -211,6 +206,7 @@ int fmrib_main(int argc, char *argv[])
   } else if (outform==1) {
     volume<int> relabeledim;
     std::vector<int> sizeorder;
+    size[0] = 0;
     get_sizeorder(size,sizeorder);
     relabel_image(labelim,relabeledim,sizeorder);
     save_volume(relabeledim,argv[2]);
@@ -230,11 +226,11 @@ int fmrib_main(int argc, char *argv[])
     // process the COPE image
     volume<T> vin;
     read_volume(vin,argv[3]);
-    std::vector<int> copesize, copemaxx, copemaxy, copemaxz;
-    std::vector<float> copecogx, copecogy, copecogz;
-    std::vector<T> copemax, copemean;
-    get_stats(labelim,vin,copesize,copemax,copemean,
-	      copemaxx,copemaxy,copemaxz,copecogx,copecogy,copecogz,invert);
+    std::vector<int> copesize;
+    std::vector<tripleint> copemax;
+    std::vector<triplefloat> copecog;
+    std::vector<T> copemaxval, copemean;
+    get_stats(labelim,vin,copesize,copemaxval,copemean,copemax,copecog,invert);
 
     // re-threshold for p
     size[0] = 0;
@@ -266,12 +262,12 @@ int fmrib_main(int argc, char *argv[])
 	pthreshindex[index] = n - nozeroclust;
 	cout << pthreshindex[index] << "\t" << k << "\t" 
 	     << maxvals[index] << "\t" 
-	     << p << "\t" << -log10(p) << "\t" << maxx[index] << "\t" 
-	     << maxy[index] << "\t" << maxz[index] << "\t"
-	     << cogx[index] << "\t" << cogy[index] << "\t" 
-	     << cogz[index] << "\t" << copemax[index] << "\t"
-	     << copemaxx[index] << "\t" << copemaxy[index] << "\t" 
-	     << copemaxz[index] << "\t" << copemean[index] << endl;
+	     << p << "\t" << -log10(p) << "\t" << max[index].x << "\t" 
+	     << max[index].y << "\t" << max[index].z << "\t"
+	     << cog[index].x << "\t" << cog[index].y << "\t" 
+	     << cog[index].z << "\t" << copemaxval[index] << "\t"
+	     << copemax[index].x << "\t" << copemax[index].y << "\t" 
+	     << copemax[index].z << "\t" << copemean[index] << endl;
       }
     }
       
@@ -314,4 +310,3 @@ int main(int argc,char *argv[])
   return call_fmrib_main(dtype(argv[1]),argc,argv);
 }
 
-