diff --git a/test.cc b/test.cc
index 8cc9a075b7ee344ffa917500f953d172fb4aa656..88f43211c05736a8be7d845f0e9cc73c1e2c622a 100644
--- a/test.cc
+++ b/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;
}