1 | function smooth_data_out = fmri_smooth(data,fwhm_sz,voxsize,out_file) |
---|
2 | |
---|
3 | if nargin==4 |
---|
4 | writeSmoothedData = true; |
---|
5 | else |
---|
6 | writeSmoothedData = false; |
---|
7 | end |
---|
8 | |
---|
9 | if nargin<3 || isempty(voxsize) |
---|
10 | fwhmInPixels = true; |
---|
11 | else |
---|
12 | fwhmInPixels = false; |
---|
13 | end |
---|
14 | |
---|
15 | % Check if data is an Aedes structure |
---|
16 | if isstruct(data) |
---|
17 | data = data.FTDATA; |
---|
18 | elseif ischar(data) |
---|
19 | data = aedes_data_read(data); |
---|
20 | data = data.FTDATA; |
---|
21 | end |
---|
22 | |
---|
23 | UseImFilter = true; |
---|
24 | |
---|
25 | smooth_data = zeros(size(data)); |
---|
26 | |
---|
27 | % Check data dimensions |
---|
28 | if any(ndims(data)==[4,3]) && size(data,3)~=1 |
---|
29 | use3Dkernel = true; |
---|
30 | else |
---|
31 | use3Dkernel = false; |
---|
32 | end |
---|
33 | |
---|
34 | % Calculate standard deviations using FWHM |
---|
35 | if fwhmInPixels |
---|
36 | stds = fwhm_sz/sqrt(8*log(2)); |
---|
37 | else |
---|
38 | stds = (fwhm_sz/sqrt(8*log(2)))./voxsize; |
---|
39 | end |
---|
40 | |
---|
41 | % Calculate kernel size using STDs |
---|
42 | kernel_sz = round(6*stds); |
---|
43 | if kernel_sz(3)>floor((size(data,3)/2)) |
---|
44 | kernel_sz(3)=floor((size(data,3)-1)/2); |
---|
45 | end |
---|
46 | kernel_sz |
---|
47 | |
---|
48 | % Construct the smoothing kernel |
---|
49 | if use3Dkernel |
---|
50 | [x,y,z] = meshgrid(-kernel_sz(2):kernel_sz(2),... |
---|
51 | -kernel_sz(1):kernel_sz(1),... |
---|
52 | -kernel_sz(3):kernel_sz(3)); |
---|
53 | s_kernel = exp(-(x).^2/(2*(stds(1)).^2)... |
---|
54 | -(y).^2/(2*(stds(2)).^2)... |
---|
55 | -(z).^2/(2*(stds(3)).^2)); |
---|
56 | s_kernel = s_kernel/sum(s_kernel(:)); |
---|
57 | else |
---|
58 | [x,y] = meshgrid(-kernel_sz(2):kernel_sz(2),... |
---|
59 | -kernel_sz(1):kernel_sz(1)); |
---|
60 | s_kernel = exp(-(x).^2/(2*(stds(1)).^2)... |
---|
61 | -(y).^2/(2*(stds(2)).^2)); |
---|
62 | s_kernel = s_kernel/sum(s_kernel(:)); |
---|
63 | end |
---|
64 | |
---|
65 | % Smooth the image data ----------------------- |
---|
66 | nVols = size(data,4); |
---|
67 | if ~UseImFilter |
---|
68 | tmp_sz = size(data(:,:,:,1)); |
---|
69 | %tmp_sz(3)=tmp_sz(3)+fwhm_sz(3); |
---|
70 | ind_hi=ceil(tmp_sz/2)+floor(size(s_kernel)/2); |
---|
71 | ind_lo=ind_hi-(size(s_kernel)-1); |
---|
72 | ind_lo(ind_lo<1)=1; |
---|
73 | tmp_kernel = zeros(tmp_sz); |
---|
74 | tmp_kernel((ind_lo(1):ind_hi(1))+1,(ind_lo(2):ind_hi(2))+1,... |
---|
75 | (ind_lo(3):ind_hi(3)))=s_kernel; |
---|
76 | end |
---|
77 | for ii=1:nVols |
---|
78 | if ii==1 |
---|
79 | fprintf(1,'Smoothing volume %d/%d',ii,nVols); |
---|
80 | bsz = length(sprintf('%d/%d',ii,nVols)); |
---|
81 | else |
---|
82 | fprintf(1,repmat('\b',1,bsz)); |
---|
83 | bsz = length(sprintf('%d/%d',ii,nVols)); |
---|
84 | fprintf(1,'%d/%d',ii,nVols); |
---|
85 | end |
---|
86 | if UseImFilter |
---|
87 | tmp_data = data(:,:,:,ii); |
---|
88 | smooth_data(:,:,:,ii) = imfilter(tmp_data,s_kernel); |
---|
89 | else |
---|
90 | tmp_data = double(data(:,:,:,ii)); |
---|
91 | %pad_sz=size(tmp_data)+fwhm_sz*2; |
---|
92 | %tmp_data(pad_sz(1),pad_sz(2),pad_sz(3))=0; |
---|
93 | %tmp_data = cat(3,tmp_data(:,:,1),tmp_data,tmp_data(:,:,end)); |
---|
94 | tmp_smooth = fftshift(ifftn(fftn(tmp_data).*fftn(tmp_kernel))); |
---|
95 | smooth_data(:,:,:,ii) = tmp_smooth(:,:,1:size(data,3)); |
---|
96 | end |
---|
97 | end |
---|
98 | fprintf(1,'\n'); |
---|
99 | |
---|
100 | if nargout>0 |
---|
101 | smooth_data_out = smooth_data; |
---|
102 | end |
---|
103 | |
---|
104 | % Write output files |
---|
105 | if writeSmoothedData |
---|
106 | aedes_write_nifti(smooth_data,out_file) |
---|
107 | end |
---|