-
Notifications
You must be signed in to change notification settings - Fork 25
/
hrf.pas
executable file
·166 lines (152 loc) · 6.36 KB
/
hrf.pas
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
unit hrf;
interface
uses
define_types, utypes,metagraph;
const
kHRFdur = 24000; //ms for 'full' HRF - window size for HRF
function CreateHRF (lTRsec: double; var lKernelBins: integer; lDefaultsStatsFmriT: integer; var lHRFra, lTDra: doublep): boolean;
function ConvolveTimeCourse(var lTimeCourse: PMatrix; var lKernel: doublep; var l4DTrace: T4DTrace;
lCond,lCondOut,lnVol,lKernelBins,lDefaultsStatsFmriT,lDefaultsStatsFmriT0: integer;
lTRSec: single; lSliceTime: boolean): boolean;
implementation
uses math {power}, ugamma {gamma},sysutils,dialogs;
const
kHRFkernelSec = 32;
//kDefaultsStatsFmriT = 16; //each TR is supersampled at 16x resolution
//from SPM's hrf.m
function spm_Gpdf(x,h,l: double): double;
//emulates spm_Gpdf
begin
result := power(l,h)*power(x,(h-1))* exp(-l*x);
result := result / gamma(h);
end;
function fHRF (u,dt: double): double;
//emulates spm_hrf.m
const
//TR = 1;
p1= 6; //delay of response
p2=16;//delay of undershoot (relative to onset)
p3=1; //dispersion of response
p4=1; //dispersion of undershoot
p5=6; //ratio of response to undershoot
p7=kHRFkernelSec;//length of kernel (seconds)
begin
if u <= 0 then
result := 0
else
result := spm_Gpdf(u,p1/p3,dt/p3) - spm_Gpdf(u,p2/p4,dt/p4)/p5;
end;
function CreateHRF (lTRsec: double; var lKernelBins: integer; lDefaultsStatsFmriT: integer; var lHRFra, lTDra: doublep): boolean;
//NOTE: if this returns TRUE, you MUST freemem lHRFra, lTDra
//returns lHRFra and lTDra with lBins of data - equal to 32sec convolution kernel for
//hemodynamic response (HRF) and the HRF's temporal derivative
var
lDT,lSum,l1sec: double;
lI: integer;
begin
result := false;
if lDefaultsStatsFmriT < 1 then exit;
lDT := (lTRsec / lDefaultsStatsFmriT); //DeltaTime - width of each sample in sec
lKernelBins := round ( kHRFkernelSec / lDT);
if lKernelBins < 1 then
exit;
getmem(lHRFra,lKernelBins*sizeof(double));
//generate whole HRF kernel
for lI := 1 to lKernelBins do
lHRFra^[lI] := fHRF (lI-1,lDT);
//find sum
lSum := 0;
for lI := 1 to lKernelBins do
lSum := lSum + lHRFra^[lI];
//normalize - so sum = 1
for lI := 1 to lKernelBins do
lHRFra^[lI] := lHRFra^[lI]/lsum;
//next temporal derivative
getmem(ltdra,lKernelBins*sizeof(double));
l1sec := 1/lDT;
for lI := 1 to lKernelBins do
ltdra^[lI] := fHRF((lI-1)-l1sec,lDT); //tdHRF (lI-1,lDT);
//find sum
lSum := 0;
for lI := 1 to lKernelBins do
lSum := lSum + ltdra^[lI];
//normalize - so sum = 1
for lI := 1 to lKernelBins do
ltdra^[lI] := ltdra^[lI]/lsum;
//temporal derivative is difference between normalized TD and normalized HRF
for lI := 1 to lKernelBins do
ltdra^[lI] := lHRFra^[lI]- ltdra^[lI];
result := true;
end;
function Convolve(var lTimeCoursePrecise,lKernel: doublep; lEventBin,lnVolPrecise,lKernelBins: integer): boolean;
var
lVol,lStart,lEnd: integer;
begin
result := false;
if (lEventBin > lnVolPrecise) then exit; //event too late to influence timecourse
if ((lEventBin+lKernelBins)< 1) then exit;//event too early to influence timecourse
lStart := lEventBin;
if lStart < 1 then
lStart := 1;
lEnd := (lEventBin+lKernelBins-1);
if lEnd > lnVolPrecise then
lEnd := lnVolPrecise;
//lOffset := lEventBin;
for lVol := lStart to lEnd do begin
lTimeCoursePrecise^[lVol] := lTimeCoursePrecise^[lVol] + lKernel^[lVol -lEventBin+1];
end;
result := true;
end;
function ConvolveTimeCourse(var lTimeCourse: PMatrix; var lKernel: doublep; var l4DTrace: T4DTrace; lCond,lCondOut, lnVol,lKernelBins,lDefaultsStatsFmriT,lDefaultsStatsFmriT0: integer;
lTRSec: single; lSliceTime: boolean): boolean;
var
lnVolPrecise,lEvent,lVol,lVolx,lEventBin,lEventEnd: integer;
lDT: double;
lTimeCoursePrecise: doublep;//supersampled by kDefaultsStatsFmriT
lAllEvents: boolean;
begin
result := false;
if (l4DTrace.Conditions[lCond].Events < 1) or (lnVol < 1) or (lTRSec <= 0) then exit;
lnVolPrecise := lnVol * lDefaultsStatsFmriT;
getmem(lTimeCoursePrecise,lnVolPrecise * sizeof(double));
for lVol := 1 to lnVolPrecise do
lTimeCoursePrecise^[lVol] := 0;
lDT := (lTRsec / lDefaultsStatsFmriT); //DeltaTime - width of each sample in sec
//spm_fmri_design
//X is supersampled at 16 times (fMRI_T) the number of volumes - with (32 bin offset)
//k = SPM.nscan(s);
//X = X([0:(k - 1)]*fMRI_T + fMRI_T0 + 32,:);
for lEvent := 1 to l4DTrace.Conditions[lCond].Events do begin
lEventBin := round((l4DTrace.Conditions[lCond].EventRA^[lEvent])/lDT);
//incorrect: same dur will have different number of bins due to rounding:
//lEventEnd := round((l4DTrace.Conditions[lCond].EventRA^[lEvent]+l4DTrace.Conditions[lCond].DurRA^[lEvent])/lDT);
//correct: all stimuli of same duration will have identical number of bins
lEventEnd := lEventBin+round(l4DTrace.Conditions[lCond].DurRA^[lEvent]/lDT);
//if lEvent = 1 then fx(lEventBin,lEventEnd,l4DTrace.Conditions[lCond].DurRA^[lEvent]);
repeat
if (lEventBin > 0) and (lEventBin <= lnVolPrecise) then
Convolve(lTimeCoursePrecise,lKernel,lEventBin,lnVolPrecise,lKernelBins);
inc(lEventBin);
until lEventBin > lEventEnd;
end; //for each event
//output - scaled by reciprocal of DT: e.g. if TR=2, DT=0.125, Scale = 8
//if TR=2.2, DT=0.1375 Scale = 7.2727
//this linear scaling does not change any effects - it simply clones SPM2
lAllEvents := true;
for lEvent := 1 to l4DTrace.Conditions[lCond].Events do
if l4DTrace.Conditions[lCond].DurRA^[lEvent] > lDT then
lAllEvents := false;
if lAllEvents then
lDT := 1/lDT
else
lDT := 1;
lVolx := lDefaultsStatsFmriT0;
for lVol := 1 to lnVol do begin
if (lVolx > 0) and (lVolx < lnVolPrecise) then
lTimeCourse^[lCondOut]^[lVol] := lDT * lTimeCoursePrecise^[lVolx];
inc(lVolx,lDefaultsStatsFmriT);
end;
freemem(lTimeCoursePrecise);
result := true;
end;//func ConvolveTimeCourse
end.