From 1c46c0a80142c896f0ae94ee3870d3cba9cb5fb0 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Wed, 2 Oct 2013 09:32:08 +0000
Subject: [PATCH] *** empty log message ***
---
fsl_dualreg.cc | 254 +++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 254 insertions(+)
create mode 100644 fsl_dualreg.cc
diff --git a/fsl_dualreg.cc b/fsl_dualreg.cc
new file mode 100644
index 0000000..bc0fb91
--- /dev/null
+++ b/fsl_dualreg.cc
@@ -0,0 +1,254 @@
+/* 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;
+ }
+ }
+
--
GitLab