source: misclib/fmri_smooth.m @ 119

Last change on this file since 119 was 119, checked in by tjniskan, 9 years ago
  • Added support for triple-reference EPI data.
  • Added two new properties for aedes_readfid:
    • FileChunkSize? for reducing memory requirements when reading large fid-files (default FileChunkSize?=500 MB)
    • EPIBlockSize for reducing memory requirements when reading large multireceiver EPI (default=100 volumes)
  • Tried to speed up smoothing in fMRI analysis by performing the

operation in frequency space. Still needs some work...

M misclib/fmri_smooth.m
M misclib/fmri_analysis.m
M aedes_readfid.m
M aedes_revision.m

File size: 2.6 KB
Line 
1function smooth_data_out = fmri_smooth(data,fwhm_sz,voxsize,out_file)
2
3if nargin==4
4  writeSmoothedData = true;
5else
6  writeSmoothedData = false;
7end
8
9if nargin<3 || isempty(voxsize)
10  fwhmInPixels = true;
11else
12  fwhmInPixels = false;
13end
14
15% Check if data is an Aedes structure
16if isstruct(data)
17  data = data.FTDATA;
18elseif ischar(data)
19  data = aedes_data_read(data);
20  data = data.FTDATA;
21end
22
23UseImFilter = true;
24
25smooth_data = zeros(size(data));
26
27% Check data dimensions
28if any(ndims(data)==[4,3]) && size(data,3)~=1
29  use3Dkernel = true;
30else
31  use3Dkernel = false;
32end
33
34% Calculate standard deviations using FWHM
35if fwhmInPixels
36   stds = fwhm_sz/sqrt(8*log(2));
37else
38  stds = (fwhm_sz/sqrt(8*log(2)))./voxsize;
39end
40
41% Calculate kernel size using STDs
42kernel_sz = round(6*stds);
43if kernel_sz(3)>floor((size(data,3)/2))
44  kernel_sz(3)=floor((size(data,3)-1)/2);
45end
46kernel_sz
47
48% Construct the smoothing kernel
49if 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(:));
57else
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(:));
63end
64
65% Smooth the image data -----------------------
66nVols = size(data,4);
67if ~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;
76end
77for 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
97end
98fprintf(1,'\n');
99
100if nargout>0
101  smooth_data_out = smooth_data;
102end
103
104% Write output files
105if writeSmoothedData
106  aedes_write_nifti(smooth_data,out_file)
107end
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008