1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
|
#include "correlation.h"
fmriResult fmri_eval(const Data<float,1>& timecourse, const Data<float,1>& designvec) {
Log<OdinData> odinlog("","fmri_eval");
fmriResult result;
if( !same_shape(timecourse, designvec) ) {
ODINLOG(odinlog,errorLog) << "design file size mismatch" << STD_endl;
return result;
}
int nrep=timecourse.numElements();
float maxdesign=max(designvec);
float mindesign=min(designvec);
// Calculate baseline
int nbaseline=0;
for(int i=0; i<nrep; i++) {
if(designvec(i)) break;
nbaseline++;
}
if(nbaseline) result.Sbaseline=mean(timecourse(Range(0,nbaseline-1)));
// calculate fMRI contrast
int nrest=0;
int nstim=0;
for(int i=0; i<nrep; i++) {
if(designvec(i)==mindesign) nrest++;
if(designvec(i)==maxdesign) nstim++;
}
Array<float,1> restdata(nrest);
Array<float,1> stimdata(nstim);
int irest=0;
int istim=0;
for(int i=0; i<nrep; i++) {
if(designvec(i)==mindesign) {restdata(irest)=timecourse(i); irest++;}
if(designvec(i)==maxdesign) {stimdata(istim)=timecourse(i); istim++;}
}
statisticResult reststat=statistics(restdata);
statisticResult stimstat=statistics(stimdata);
result.Srest=reststat.mean;
result.Sstim=stimstat.mean;
result.rel_diff=secureDivision(stimstat.mean-reststat.mean,reststat.mean);
result.rel_err=secureDivision(stimstat.meandev+reststat.meandev,reststat.mean);
return result;
}
//////////////////////////////////////////////////////////////////////////////////
// instantiate one special template for checking at compile time
#ifdef ODIN_DEBUG
template correlationResult correlation(const Array<float,1>&, const Array<float,1>&);
#endif
|
