*** empty log message ***

fsl_dualreg.cc

+/*  test.cc
+  
+    Christian F. Beckmann, FMRIB Analysis Group
+  
+    Copyright (C) 1999-20013 University of Oxford */
+
+/*  CCOPYRIGHT  */
+
+#include "libvis/miscplot.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "utils/options.h"
+#include <vector>
+#include <ctime>
+#include "newimage/newimageall.h"
+#include "melhlprfns.h"
+#include <iostream>
+
+#ifdef __APPLE__
+#include <mach/mach.h>
+#define memmsg(msg) { \
+  struct task_basic_info t_info; \
+  mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+  if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+	{ \
+		cout << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+		cout.flush(); \
+	} \
+}
+#else
+#define memmsg(msg) { \
+   cout << msg; \
+}
+#endif
+
+// a simple message macro that takes care of cout and log
+#define message(msg) { \
+  cout << msg; \
+  cout.flush(); \
+}
+
+#define outMsize(msg,Mat) { \
+    cerr << "     " << msg << "  " <<Mat.Nrows() << " x " << Mat.Ncols() << endl;	\
+}
+
+
+using namespace MISCPLOT;
+using namespace MISCMATHS;
+using namespace Utilities;
+using namespace std;
+using namespace Melodic;
+
+// GLOBALS
+
+clock_t tictime;
+
+
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+
+  string title=string("fsl_BLAH")+
+		string("\nAuthor: Christian F. Beckmann \nCopyright(c) 2008-2013 University of Oxford\n")+
+		string(" \n \n")+
+		string(" \n");
+  string examples="fsl_BLAH [options]";
+
+//Command line Options {
+    Option<string> fnin(string("-i,--in"), string(""),
+		string("input file name (matrix 3D or 4D image)"),
+		false, requires_argument);
+    Option<string> fnmask(string("-m"), string(""),
+			string("mask file name "),
+			false, requires_argument);
+	Option<int> help(string("-h,--help"), 0,
+		string("display this help text"),
+		false,no_argument);
+	Option<int> xdim(string("-x,--xdim"), 0,
+		string("xdim"),
+		false,requires_argument);
+	Option<int> ydim(string("-y,--ydim"), 0,
+		string("ydim"),
+		false,requires_argument);
+	Option<int> econ(string("-e,--econ"), 0,
+		string("econ: how to liump stuff"),
+		false,requires_argument);
+		/*
+}
+*/
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+
+void tic(){
+	tictime = clock();
+}
+
+void toc(){
+	cerr << endl << "TOC: " << float(clock()-tictime)/CLOCKS_PER_SEC << " seconds" << endl<<endl;
+}
+
+Matrix calccorr(const Matrix& in, int econ)
+  { 
+    Matrix Res;
+	int nrows=in.Nrows();
+	int ncols=in.Ncols();    
+    Res = zeros(nrows,nrows);
+
+    if(econ>0){
+      RowVector colmeans(ncols);
+	  for (int n=1; n<=ncols; n++) {
+        colmeans(n)=0;
+        for (int m=1; m<=nrows; m++) {
+          colmeans(n)+=in(m,n);
+        }
+        colmeans(n)/=nrows;
+      }
+      int dcol = econ;
+	  Matrix suba; 
+
+      for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
+	    suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
+		int scolmax = suba.Ncols();
+
+		for (int n=1; n<=scolmax; n++) {
+	        double cmean=colmeans(ctr + n - 1);
+	        for (int m=1; m<=nrows; m++) {
+	          suba(m,n)-=cmean;
+	        }
+	    }
+		
+	    Res += suba*suba.t() / ncols;
+      }
+    }
+    else
+      Res = cov(in.t());
+    return Res;
+  }  //Matrix calccorr
+
+void convert_to_mat(volume4D<float>& in, volume<float>& mask, Matrix& out)
+{
+	out = in[0].vec(mask).t();
+	in.deletevolume(0);
+	while(in.tsize()>0){
+		out &= in[0].vec(mask).t();
+		in.deletevolume(0);	
+	}
+}
+
+int do_work(int argc, char* argv[]) {
+	
+
+	tic();
+	Matrix MatrixData;
+	volume<float> Mean;
+  	
+	if(xdim.value()==0 && ydim.value()==0)
+    {
+		volume4D<float> RawData;
+		volume<float> theMask;
+		toc();
+    	//read data
+    	message("Reading data file " << (string)fnin.value() << "  ... ");
+    	read_volume4D(RawData,fnin.value());
+    	message(" done" << endl);
+
+		Mean = meanvol(RawData);
+		toc();
+    	message("Reading mask file " << (string)fnmask.value() << "  ... ");
+      	read_volume(theMask,fnmask.value());
+
+		MatrixData = RawData.matrix();
+
+		memmsg("start");
+		Matrix res;
+		convert_to_mat(RawData,theMask, res);
+
+		memmsg(" Before reshape ");
+    
+		write_ascii_matrix("mat1", res);
+		write_ascii_matrix("mat2", MatrixData);
+		
+
+    }
+	else{
+		Matrix data = unifrnd(xdim.value(),ydim.value());
+	    outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+		
+		data = unifrnd(10,100);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+		
+		data = unifrnd(100,1000);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+		
+		data = unifrnd(100,10000);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+
+		data = unifrnd(300,200000);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();		
+		
+		data = unifrnd(500,20000);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+		
+		data = unifrnd(500,200000);
+		outMsize("data", data);
+		tic(); calccorr(data,econ.value()); toc();
+		
+		
+	}
+	
+	
+	return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+	int main(int argc,char *argv[]){
+	  Tracer tr("main");
+	  OptionParser options(title, examples);
+	  try{
+	    // must include all wanted options here (the order determines how
+	    //  the help message is printed)
+	
+			options.add(fnin);		
+			options.add(fnmask);		
+			options.add(help);
+			options.add(xdim);
+			options.add(ydim);
+			options.add(econ);
+	    options.parse_command_line(argc, argv);
+
+	    // line below stops the program if the help was requested or 
+	    //  a compulsory option was not set
+	    if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+				options.usage();
+				exit(EXIT_FAILURE);
+	    }else{
+	  		// Call the local functions
+	  		return do_work(argc,argv);
+			}
+		}catch(X_OptionError& e) {
+			options.usage();
+	  	cerr << endl << e.what() << endl;
+	    exit(EXIT_FAILURE);
+	  }catch(std::exception &e) {
+	    cerr << e.what() << endl;
+	  } 
+	}
+