*** empty log message ***

test.cc

@@ -72,18 +72,43 @@ void usage(void)
   exit(1);
 }
 
+ReturnMatrix normpdf3(const RowVector& vals, const RowVector& mu, const RowVector& var)
+{
+  Matrix res(mu.Ncols(),vals.Ncols());
+  for (int mc=1; mc<=res.Ncols(); mc++){
+    for (int mr=1; mr<=res.Nrows(); mr++){
+      res(mr,mc) = std::exp(-0.5*(std::pow(vals(mc)-mu(mr),2)/var(mr)))*std::pow(2*M_PI*var(mr),-0.5);
+    }
+  }
+
+  res.Release();
+  return res;
+}
+
+
+
+
 int main(int argc, char *argv[])
 {
 
- double tmptime = time(NULL);
-    cerr << tmptime << endl << endl;;
-    srand((unsigned int) tmptime);
-   
-    Matrix test(5,5);
-    test = unifrnd(5,5);
-    cerr<< test;
-    exit(1);
+  RowVector mutmp(1);
+  RowVector vartmp(1);
+
+  RowVector valtmp(1);
+
+  mutmp(1) =0.0;
+  vartmp(1) =1.0;
+
+  float mintmp=-1;
+  float inttmp=0.1;
 
+  for(int ctr=1; ctr<=5; ctr++){
+    valtmp(1) = mintmp + ctr*inttmp;
+    cerr << valtmp << "    " << normpdf(valtmp,mutmp,vartmp) << endl;
+    // cerr << valtmp << "    " << normpdf(valtmp,mutmp,vartmp) << endl;
+  }
+
+  //exit(1);
 
   if (argc<3)
     usage();
@@ -104,9 +129,15 @@ int main(int argc, char *argv[])
   if (argc>3)
     MIXfname = string(argv[3]);
 
-  cerr << argc << "  " << RAWfname << " " <<ICfname << " " << MIXfname << endl;
+  string MASKfname;
+  if (argc>4)
+    MASKfname = string(argv[4]);
+
+
+
+  cerr << argc << "  " << RAWfname << " " <<ICfname << " " << MIXfname << MASKfname << endl;
  
-  {
+  /* {
     volume4D<float> RawData;
     cout << " Reading orig. data " << RAWfname << " ... ";
     read_volume4D(RawData,RAWfname,ICvolInfo);
@@ -137,7 +168,7 @@ int main(int argc, char *argv[])
 
 
     cout << " done! " << endl;
-  }
+    }*/
 
    {
     volume4D<float> RawIC;
@@ -146,10 +177,13 @@ int main(int argc, char *argv[])
     cout << " done" << endl;
 
     cout << "Creating mask   ... ";
-    Mask = binarise(RawIC[0],float(RawIC[0].min()),float(RawIC[0].max()));
+
+    read_volume(Mask,MASKfname);
+
+    //  Mask = binarise(RawIC[0],float(RawIC[0].min()),float(RawIC[0].max()));
 
     ICs = RawIC.matrix(Mask);
-    if(ICs.Nrows()>1){
+    /*if(ICs.Nrows()>1){
       Matrix DStDev=stdev(ICs);
       
       volume4D<float> tmpMask;
@@ -170,9 +204,9 @@ int main(int argc, char *argv[])
       Mask = binarise(RawIC[0],float(0.001),float(RawIC[0].max())) 
 	+ binarise(RawIC[0],float(RawIC[0].min()),float(-0.001));
       ICs = RawIC.matrix(Mask);
-    }
+      }*/
 
-    //cerr << "ICs : " << ICs.Ncols() << ICs.Nrows() << endl;
+    cerr << "ICs : " << ICs.Ncols() << ICs.Nrows() << endl;
     cout << " done" << endl;
   }
 
@@ -187,8 +221,9 @@ int main(int argc, char *argv[])
   }
 
   cout << " ICs: " << ICs.Nrows() << " x " << ICs.Ncols() << endl;
-  if(ICs.Nrows()>0){
-    cout << " Plotting histogram for map 1" << endl;
+for(int ctr=1; ctr <= ICs.Nrows(); ctr++){ 
+
+  cout << " Plotting histogram for map " << ctr << endl;
   
     miscplot newplot;
     // newplot.add_label("legend1");
@@ -201,153 +236,19 @@ int main(int argc, char *argv[])
     Matrix mu(1,3);
     Matrix pi(1,3);
     Matrix std(1,3);
-    mu(1,1)=0;mu(1,2)=2.10;mu(1,3)=-1.99388;
-    pi(1,1)=0.911806;pi(1,2)=0.066722;pi(1,3)=0.021472;
-    std(1,1)=1;std(1,2)=0.831590;std(1,3)=0.639299;
-    newplot.gmmfit(ICs.Row(1),mu,std,pi,string("./res/g1.png"),string("Title"));
-
-    mu(1,1)=0;mu(1,2)=2.209024;mu(1,3)=-2.016637;
-    pi(1,1)=0.873335;pi(1,2)=0.086185;pi(1,3)=0.040480;
-    std(1,1)=1;std(1,2)=1.012970;std(1,3)=0.616374;
-    newplot.gmmfit(ICs.Row(2),mu,std,pi,string("./res/g2.png"),string("Title"));
-
-    mu(1,1)=0;mu(1,2)=2.694458;mu(1,3)=-3.418763;
-    pi(1,1)=0.9779397;pi(1,2)=0.017768;pi(1,3)=0.004293;
-    std(1,1)=1;std(1,2)=2.848577;std(1,3)=0.385187;
-    newplot.gmmfit(ICs.Row(3),mu,std,pi,string("./res/g3.png"),string("Title"));
-
-    // newplot.histogram(ICs.Row(1),string("./res/h1.png"),string("Title"));
-  }
+    mu(1,1)=0;mu(1,2)=1.10;mu(1,3)=-1.1;
+    pi(1,1)=0.96;pi(1,2)=0.02;pi(1,3)=0.02;
+    std(1,1)=1;std(1,2)=0.1;std(1,3)=0.1;
 
-  cout << endl << endl;
-  /*
-  Matrix data;
- 
-  cerr<< "Reading mixing matrix ";
-  data = read_ascii_matrix(string("mixmat"));
-  if (data.Storage()<=0) {
-    cerr <<" Please specify the mixing matrix correctly" << endl;
-  }
-  cerr << "done " << endl;
-  cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;
-
-
-  //set up data
-
-  if(data.Ncols()==0){
-
-    cerr << " create new data " << M_PI <<endl;
-
-    ColumnVector X(1000);
-   
-    for(int i=1; i<=X.Nrows(); i++){
-      X(i) = (float((i-1))/(X.Nrows()-1))*M_PI;
-    }
-
-    RowVector Y(10);
-    for(int i=1; i<=Y.Ncols(); i++)
-      Y(i)=3*i;
-
-    data = X*Y;
-
-    for(int i=1; i<=data.Ncols(); i++)
-      for(int j=1; j<=data.Nrows(); j++)
-	data(j,i) = sin(data(j,i));
+    //   Matrix returnM;
+    // returnM = normpdf(ICs.Row(ctr)*10,10,11);
 
-   cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;
-
-   Log mmat;
-mmat.setDir(string("./res/"),string("mix.txt"));
- mmat << data <<endl;
+    newplot.gmmfit(ICs.Row(ctr),mu,std,pi,string("./res/g"+num2str(ctr)+".png"),string("Title"));
+    //newplot.histogram(ICs.Row(ctr),string("./res/g"+num2str(ctr)+".png"),string("Title"));
 
   }
 
-  Log IChtml;
-  IChtml.setDir(string("./res/"),string("index.html"));
-
-  IChtml << "<HTML> " << endl
-	     << "<TITLE>MELODIC Component " 
-	     << "</TITLE>" << endl
-	     << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR") 
-	     << "/doc/images/fsl-bg.jpg\">" << endl 
-	     << "<hr><CENTER><H1>MELODIC Component "
-	     << "</H1>"<< endl;
-
-
-  float tr = 5.0;
-
-  Matrix data2;
-  data2 = calcFFT(data);
-
-  for (int i=1; i<=data.Ncols(); i++)
-    {
-      
-      {//plot time course
-	miscplot newplot;
-	//newplot.add_label("legend1");
-	//newplot.add_label("legend2");
-	//        newplot.add_ylabel(string("yll1"));
-        //newplot.add_ylabel(string("ylabel2guiygyiasyugfuyigfuigasgfuiguiyasgfuigasuigfagifgiuygasisgfuigasgfgasgfigiagfgaggfigaif"));
-	//newplot.add_xlabel(string("xl1"));
-	//newplot.add_xlabel(string("xlabel2"));
-	//newplot.set_xysize(600,200);
-	if(tr>0.0)
-	  newplot.timeseries(data.Column(i).t(),
-			     string("./res/t")+num2str(i)+".png",
-			     string("Timecourse (in seconds); TR = ")+
-			     float2str(tr,0,2,0)+" s", 
-			     tr);
-	else
-	  newplot.timeseries(data.Column(i).t(),
-			     string("./res/t")+num2str(i)+".png",
-			     string("Timecourse (in TRs)"));
-	//	write_ascii_matrix(string("./res/t")+num2str(i)+".txt"),
-	//	   data.Column(i));
-	Log tfile;
-	tfile.setDir(string("./res/"),string("t")+num2str(i)+".txt");
-	tfile << data.Column(i) << endl;
-
-	IChtml << "<A HREF=\"" << string("t")
-	  +num2str(i)+".txt" << "\"> ";
-	IChtml << "<img BORDER=0 SRC=\"" 
-	  +string("t")+num2str(i)+".png\"></A><p>" << endl;
-      }//time series plot
-      
-      {//plot frequency 
-	miscplot newplot;
-	int fact = int(std::pow(10.0,int(std::log10(float(data.Nrows())))));
-
-	if(tr>0.0)
-	  newplot.timeseries(data2.Column(i).t(),
-			     string("./res/f")+num2str(i)+".png",
-			     string("FFT of timecourse (in Hz / ") +num2str(fact)+")",
-			     fact/(tr*data2.Nrows()));
-	else
-	  newplot.timeseries(data2.Column(i).t(),
-			     string("./res/f")+num2str(i)+".png",
-			     string(string("FFT of timecourse (in cycles); ")
-				  +"frequency(Hz)=cycles/("
-				  +num2str(data2.Nrows())
-				  +"* TR); period(s)=("
-				  +num2str(data2.Nrows())
-				  +"* TR)/cycles"));
-
-	Log ffile;
-	ffile.setDir(string("./res/"),string("f")+num2str(i)+".txt");
-	ffile << data2.Column(i) << endl;
-
-	IChtml << "<A HREF=\"" << string("f")
-	  +num2str(i)+".txt" << "\"> ";
-	IChtml << "<img BORDER=0 SRC=\"" 
-	  +string("f")+num2str(i)+".png\"></A><p>" << endl;
-      }//frequency plot
-
-    }
-
+  cout << endl << endl;
 
-  IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "	
-	<< "FMRIB Software Library</A>.</FONT>" << endl
-	<< "</BODY></HTML>" << endl;
-	cerr << endl;  */
   return 0;
 }