1 | function maps_out = fmri_analysis(DATA,varargin) |
---|
2 | % |
---|
3 | % This function does a quick fMRI analysis in SPM style. No motion |
---|
4 | % correction or slice timing at this time... |
---|
5 | % |
---|
6 | |
---|
7 | % Defaults |
---|
8 | TR = 2.039; |
---|
9 | contr = [1 0]; |
---|
10 | onset = [30 75 120].'; |
---|
11 | durat = [15 15 15].'; |
---|
12 | smooth_kernel = [2 2 1]; |
---|
13 | |
---|
14 | |
---|
15 | if rem(length(varargin),2)~=0 |
---|
16 | error('parameters/values missing!') |
---|
17 | end |
---|
18 | |
---|
19 | % Parse varargin |
---|
20 | for ii=1:2:length(varargin) |
---|
21 | param = lower(varargin{ii}); |
---|
22 | value = varargin{ii+1}; |
---|
23 | switch param |
---|
24 | case 'TR' |
---|
25 | TR = value; |
---|
26 | case {'onset','onsets'} |
---|
27 | onset = value; |
---|
28 | onset = onset(:); |
---|
29 | case {'durat','durations'} |
---|
30 | durat = value; |
---|
31 | durat = durat(:); |
---|
32 | case 'smooth' |
---|
33 | smooth_kernel = value; |
---|
34 | case 'contrast' |
---|
35 | contr = value; |
---|
36 | end |
---|
37 | end |
---|
38 | |
---|
39 | if length(durat)==1 |
---|
40 | durat = repmat(durat,1,length(onset)); |
---|
41 | elseif length(durat)~=length(onset) |
---|
42 | error('Mismatch in the number of elements in onset and duration!') |
---|
43 | end |
---|
44 | |
---|
45 | |
---|
46 | % Check data |
---|
47 | if isstruct(DATA) && isfield(DATA,'FTDATA') |
---|
48 | % if isfield(DATA,'PROCPAR') && isfield(DATA.PROCPAR,'readres') |
---|
49 | % data = DATA.FTDATA(:,:,:,2:end); |
---|
50 | % else |
---|
51 | data = DATA.FTDATA; |
---|
52 | % end |
---|
53 | elseif isnumeric(DATA) && ndims(DATA)>2 |
---|
54 | data=DATA; |
---|
55 | else |
---|
56 | error('Input data has to be 3D or 4D numeric matrix!') |
---|
57 | end |
---|
58 | |
---|
59 | if ndims(data)==3 |
---|
60 | data = permute(data,[1 2 4 3]); |
---|
61 | end |
---|
62 | |
---|
63 | % Load Rat HRF (8 seconds long) |
---|
64 | bf = [ |
---|
65 | 0 |
---|
66 | 0.000008876945146 |
---|
67 | 0.000387621003097 |
---|
68 | 0.003012441669394 |
---|
69 | 0.011548214895067 |
---|
70 | 0.030056522080168 |
---|
71 | 0.061233628208021 |
---|
72 | 0.105349995559045 |
---|
73 | 0.160158500633323 |
---|
74 | 0.221528149631163 |
---|
75 | 0.284405279593501 |
---|
76 | 0.343761753554314 |
---|
77 | 0.395323573307322 |
---|
78 | 0.436002661868442 |
---|
79 | 0.464045998242252 |
---|
80 | 0.478965423658889 |
---|
81 | 0.481327079129854 |
---|
82 | 0.472473786978796 |
---|
83 | 0.454237528221769 |
---|
84 | 0.428680153860941 |
---|
85 | 0.397883140401152 |
---|
86 | 0.363793656845973 |
---|
87 | 0.328124931280142 |
---|
88 | 0.292303457117353 |
---|
89 | 0.257453130446147 |
---|
90 | 0.224406054242129 |
---|
91 | 0.193730677047637 |
---|
92 | 0.165769521647024 |
---|
93 | 0.140680556701329 |
---|
94 | 0.118477987483032 |
---|
95 | 0.099069734765398 |
---|
96 | 0.082290068881472 |
---|
97 | 0.067926763712502 |
---|
98 | 0.055742762013807 |
---|
99 | 0.045492744488986 |
---|
100 | 0.036935220196777 |
---|
101 | 0.029840852217657 |
---|
102 | 0.023997740489123 |
---|
103 | 0.019214335904674 |
---|
104 | 0.015320580887648 |
---|
105 | 0.012167779438633 |
---|
106 | 0.009627606069476 |
---|
107 | 0.007590575511228 |
---|
108 | 0.005964217696074 |
---|
109 | 0.004671136982604 |
---|
110 | 0.003647081075751 |
---|
111 | 0.002839102788446 |
---|
112 | 0.002203865389816 |
---|
113 | 0.001706118264261 |
---|
114 | 0.001317352436400 |
---|
115 | 0.001014633776781 |
---|
116 | 0.000779604140824 |
---|
117 | 0.000597636251764 |
---|
118 | 0.000457125954290 |
---|
119 | 0.000348904854607 |
---|
120 | 0.000265756797153 |
---|
121 | 0.000202022712087 |
---|
122 | 0.000153279812313 |
---|
123 | 0.000116082720651 |
---|
124 | 0.000087755728843 |
---|
125 | 0.000066226941806 |
---|
126 | 0.000049896490408 |
---|
127 | 0.000037532277385 |
---|
128 | ]; |
---|
129 | |
---|
130 | % Create stimulus function (32 bin offset) |
---|
131 | k = size(data,4); % Number of scans |
---|
132 | T = 16; |
---|
133 | dt = TR/T; |
---|
134 | u = onset.^0; |
---|
135 | if ~any(durat) |
---|
136 | u = u/dt; |
---|
137 | end |
---|
138 | ton = round(onset*TR/dt) + 32; % onsets |
---|
139 | tof = round(durat*TR/dt) + ton + 1; % offset |
---|
140 | sf = zeros((k*T + 128),size(u,2)); |
---|
141 | ton = max(ton,1); |
---|
142 | tof = max(tof,1); |
---|
143 | for j = 1:length(ton) |
---|
144 | if numel(sf)>ton(j), |
---|
145 | sf(ton(j),:) = sf(ton(j),:) + u(j,:); |
---|
146 | end; |
---|
147 | if numel(sf)>tof(j), |
---|
148 | sf(tof(j),:) = sf(tof(j),:) - u(j,:); |
---|
149 | end; |
---|
150 | end |
---|
151 | sf = cumsum(sf); % integrate |
---|
152 | |
---|
153 | % Convolve stimulus with the HRF |
---|
154 | conv_sf = conv(sf,bf); |
---|
155 | |
---|
156 | % Resample the convolved stimulus |
---|
157 | R = conv_sf([0:(k-1)]*16+1+32); |
---|
158 | R=R(:); |
---|
159 | H = [R ones(size(R))]; |
---|
160 | |
---|
161 | % Smooth data |
---|
162 | if all(smooth_kernel) |
---|
163 | %fprintf(1,'Smoothing data...'); |
---|
164 | smooth_data = fmri_smooth(data,smooth_kernel); |
---|
165 | %fprintf(1,'...done\n'); |
---|
166 | else |
---|
167 | fprintf(1,'fMRI analysis warning: Could not smooth data!\n'); |
---|
168 | smooth_data = data; |
---|
169 | end |
---|
170 | |
---|
171 | % Calculate parametric map(s) |
---|
172 | nPlanes = size(smooth_data,3); |
---|
173 | nCols = size(smooth_data,2); |
---|
174 | pmap1 = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3)); |
---|
175 | pmap2 = zeros(size(pmap1)); |
---|
176 | tmap = zeros(size(pmap1)); |
---|
177 | c = repmat(contr,nCols,1); |
---|
178 | HTH = pinv(H'*H); |
---|
179 | R = eye(size(H,1))-H*HTH*H'; |
---|
180 | for ii=1:nPlanes |
---|
181 | fprintf(1,'Processing plane %d/%d\n',ii,nPlanes); |
---|
182 | for kk=1:nCols |
---|
183 | if kk==1 |
---|
184 | fprintf(1,'.....Processing column %03d/%03d',kk,nCols); |
---|
185 | else |
---|
186 | fprintf(1,'\b\b\b\b\b\b\b'); |
---|
187 | fprintf(1,'%03d/%03d',kk,nCols); |
---|
188 | end |
---|
189 | drawnow |
---|
190 | col_data = squeeze(smooth_data(:,kk,ii,:)).'; |
---|
191 | col_means = mean(col_data,1); |
---|
192 | col_data = detrend(col_data,'linear'); |
---|
193 | col_data = col_data+repmat(col_means,k,1); |
---|
194 | th=H\col_data; |
---|
195 | pmap1(:,kk,ii)=th(1,:).'; |
---|
196 | pmap2(:,kk,ii)=th(2,:).'; |
---|
197 | |
---|
198 | r = col_data-H*th; |
---|
199 | rr=diag(r'*r); |
---|
200 | rr=rr(:); |
---|
201 | sig2 = rr./trace(R); |
---|
202 | sig2 = repmat(sig2,1,2); |
---|
203 | T = diag(c*th)./sqrt(c.*sig2*HTH*H'*H*HTH*contr'); |
---|
204 | tmap(:,kk,ii)=T.'; |
---|
205 | end |
---|
206 | fprintf(1,'\n'); |
---|
207 | end |
---|
208 | |
---|
209 | maps_out.pmap1 = pmap1; |
---|
210 | maps_out.pmap2 = pmap2; |
---|
211 | maps_out.tmap = tmap; |
---|
212 | |
---|
213 | % Calculate effective degrees of freedom |
---|
214 | %v = (trace(R).^2)/trace(R*R); |
---|
215 | |
---|
216 | |
---|
217 | |
---|
218 | % % Calculate T-maps and contrasts |
---|
219 | % c = contr; |
---|
220 | % |
---|
221 | % % Projection and residual |
---|
222 | % HTH = pinv(H'*H); |
---|
223 | % R = eye(size(H,1))-H*HTH*H'; |
---|
224 | % r = Z-H*th; |
---|
225 | % sig2 = (r'*r)./trace(R); |
---|
226 | % |
---|
227 | % % Calculate T-value |
---|
228 | % T = (c*th)/sqrt(c*sig2*HTH*H'*H*HTH*c'); |
---|
229 | % |
---|
230 | |
---|
231 | |
---|
232 | |
---|
233 | |
---|