-
Notifications
You must be signed in to change notification settings - Fork 27
/
nii_dti_clean.m
114 lines (109 loc) · 3.35 KB
/
nii_dti_clean.m
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
function fnm = nii_dti_clean(fnm, thresh);
%Detect and remove DTI volumes with substantial movement artifacts
%Warning: only for DTI: for fMRI we need to impute missing timepoints
% fnm : name of file to check and decimate if required
%Alternative
% fsl's eddy with outlier replacement (--repol)
%Examples
% nii_dti_clean; %use GUI
% nii_dti_clean('MUSC_2.7.nii')
if ~exist('fnm','var')
%fnm = 'MUSC_2.7_AVE1_MCBI_40_dir_1007_V1_22.nii';
fnm = spm_select(1,'^.*\.(gz|nii)$','Select DTI images to clean');
end
if ~exist('thresh','var')
zThresh = 4;
end
[hdr,img, bvec,bval] = loadSub(fnm);
nvol = size(img,4);
nslice = hdr(1).dim(3);
if (nvol < 6) || (nslice < 5)
error('%s requires DTI images with at least 7 volumes and 5 slices', mfilename);
end
%collapse X/Y dimension, so we have 3 dimensions (X*Y, Z (slices), volumes
img = reshape(img,size(img,1)*size(img,2),size(img,3),size(img,4));
stdev = zeros(nvol,1);
for v = 1: nvol
dx = zeros(nslice-2,1);
%mn = zeros(nslice-2,1);
%mnTB = zeros(nslice-2,1);
for s = 2: (nslice-1)
mn = mean(img(:,s,v));
mnTB = 0.5*(mean(img(:,s-1,v)+mean(img(:,s+1,v))));
dx(s-1) = mn - mnTB;
end;
stdev(v) = std(dx)/mean(mean(img(:,:,v)));
end
stdev = stdev - mean(stdev(:)); %mean = 0
z = stdev ./ std(stdev(:));
%identify outliers
reject = zeros(nvol,1);
for v = 1: nvol
if z(v) > zThresh
fprintf('Deleting volume %d, z = %g\n', v, z(v));
reject(v) = 1;
elseif z(v) > (0.5 * zThresh)
fprintf('Please visually inspect volume %d, z = %g\n', v, z(v));
end
end
if sum(reject(:)) < 1
fprintf('%s quitting: no volumes z>%g (max z=%g) in %s\n', mfilename, zThresh, max(z(:)), fnm);
return;
end
fprintf('%s found %d volumes with z>%g in %s\n', mfilename, sum(reject(:)), zThresh, fnm);
%return raw image to 4D
img = reshape(img,hdr(1).dim(1), hdr(1).dim(2), hdr(1).dim(3), nvol);
%create decimated image
imgOK = img(:,:,:, ~reject(:));
nvolOK = sum(~reject(:));
%save image
hdrOK = hdr(1);
[p,n,x] = fileparts(hdr(1).fname);
n = [n '_' num2str(nvolOK)];
hdrOK.fname = fullfile(p, [n x]);
for v=1:nvolOK
hdrOK.n(1)=v;
spm_write_vol(hdrOK,imgOK(:, :, :, v));
end;
fnm = hdrOK.fname;
%save bvec/bval
bvecOK = bvec(~reject(:), :);
bvalOK = bval(~reject(:));
dlmwrite(fullfile(p,[ n, '.bval']), bvalOK','delimiter','\t');
dlmwrite(fullfile(p,[n, '.bvec']), bvecOK','delimiter','\t');
%end nii_dti_clean()
function [hd,im, bve,bva] = loadSub(fnm)
im = []; bve = []; bva = [];
[p,n,x] = spm_fileparts(fnm);
if (length(x)==3) && min((x=='.gz')==1)
fnm = char(gunzip(fnm));
delnam = fnm;
[p,n,x] = spm_fileparts(char(fnm));
else
delnam = '';
end
hd = spm_vol(fnm); %input header
if hd(1).dt(1) == 128
fprintf('Warning: skipping RGB image %s\n', fnm);
if ~isempty(delnam), delete(delnam); end;
return;
end
im = spm_read_vols(hd);%Input image
hd = hd(1);
if ~isempty(delnam), delete(delnam); end;
bve = loadTxt(fullfile(p,[n,'.bvec']), 3);
if isempty(bve), bve = zeros( size(im,4),3); end;
bva = loadTxt(fullfile(p,[n,'.bval']));
if isempty(bva), bva = zeros( size(im,4),1); end;
%end loadSub()
function txt = loadTxt(fnm, nRow)
txt = [];
%fprintf('%s\n', fnm);
if exist(fnm,'file') == 0, return; end;
fileID = fopen(fnm);
txt = cell2mat( textscan(fileID,'%f'));
fclose(fileID);
if exist('nRow', 'var')
txt = reshape(txt, [numel(txt)/nRow nRow]);
end
%end loadTxt()