forked from fieldtrip/fieldtrip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathft_volumebiascorrect.m
194 lines (160 loc) · 6.94 KB
/
ft_volumebiascorrect.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
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
function [mri_unbias] = ft_volumebiascorrect(cfg, mri)
% FT_VOLUMEBIASCORRECT corrects the image inhomogeneity bias in an anatomical MRI
%
% Use as
% mri_unbias = ft_volumebiascorrect(cfg, mri)
% where the input mri should be a single anatomical volume organised in a structure
% as obtained from the FT_READ_MRI function
%
% The configuration structure can contain
% cfg.spmversion = string, 'spm8', 'spm12' (default = 'spm12')
% cfg.opts = struct, containing spmversion specific options.
% See the code below and the SPM-documentation for
% more information.
%
% See also FT_VOLUMEREALIGN FT_VOLUMESEGMENT FT_VOLUMENORMALISE
% Copyright (C) 2017, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
ft_revision = '$Id$';
ft_nargin = nargin;
ft_nargout = nargout;
% do the general setup of the function
% the ft_preamble function works by calling a number of scripts from
% fieldtrip/utility/private that are able to modify the local workspace
ft_defaults % this ensures that the path is correct and that the ft_defaults global variable is available
ft_preamble init
ft_preamble debug
ft_preamble loadvar mri
ft_preamble provenance mri
if ft_abort
% do not continue function execution in case the outputfile is present and the user indicated to keep it
return
end
% ensure that the input data is valid for this function
mri = ft_checkdata(mri, 'datatype', {'volume'});
% set the defaults
cfg.spmversion = ft_getopt(cfg, 'spmversion', 'spm12');
cfg.keepintermediate = ft_getopt(cfg, 'keepintermediate', 'no');
cfg.opts = ft_getopt(cfg, 'opts');
if ~isfield(cfg, 'intermediatename')
cfg.intermediatename = tempname;
end
% check that the preferred SPM version is on the path
ft_hastoolbox(cfg.spmversion, 1);
% check whether the input has an anatomy
if ~isfield(mri, 'anatomy')
ft_error('no anatomical information available, this is required for bias correction');
end
% do an approximate alignment
mri_acpc = ft_convert_coordsys(mri, 'acpc');
% flip and permute the 3D volume itself, so that the voxel and
% headcoordinates approximately correspond this improves the convergence
% of the segmentation algorithm
[mri_acpc, permutevec, flipflags] = align_ijk2xyz(mri_acpc);
% create an spm-compatible file for the anatomical volume data
V = ft_write_mri([cfg.intermediatename '_anatomy.nii'], mri_acpc.anatomy, 'transform', mri_acpc.transform, 'spmversion', cfg.spmversion, 'dataformat', 'nifti_spm');
% get the options from hte cfg
opts = cfg.opts;
switch cfg.spmversion
case 'spm8'
opts.nbins = ft_getopt(opts, 'nbins', 256);
opts.reg = ft_getopt(opts, 'reg', 0.01);
opts.cutoff = ft_getopt(opts, 'cutoff', 30);
opts.niter = ft_getopt(opts, 'niter', 128);
T = spm_bias_estimate(V, opts);
VO = spm_bias_apply(V,T);
mri_unbias = keepfields(mri, {'coordsys', 'dim', 'transform', 'unit', 'inside'});
mri_unbias.anatomy = VO.dat.*scale;
% if strcmp(cfg.keepintermediate, 'no')
% % remove the intermediate files
% for flop=1:length(files)
% [p, f, x] = fileparts(files{flop});
% delete(fullfile(p, [f, '.*']));
% [p, f, x] = fileparts(wfiles{flop});
% delete(fullfile(p, [f, '.*']));
% end
% end
case 'spm12'
if ~isfield(cfg, 'tpm') || isempty(cfg.tpm)
cfg.tpm = fullfile(spm('dir'),'tpm','TPM.nii');
end
% create the structure that is required for spm_preproc8
opts.image = V;
opts.tpm = spm_load_priors8(cfg.tpm);
opts.biasreg = ft_getopt(opts, 'biasreg', 0.0001);
opts.biasfwhm = ft_getopt(opts, 'biasfwhm', 60);
opts.lkp = ft_getopt(opts, 'lkp', [1 1 2 2 3 3 4 4 4 5 5 5 5 6 6 ]);
opts.reg = ft_getopt(opts, 'reg', [0 0.001 0.5 0.05 0.2]);
opts.samp = ft_getopt(opts, 'samp', 3);
opts.fwhm = ft_getopt(opts, 'fwhm', 1);
Affine = spm_maff8(opts.image(1),3,32,opts.tpm,eye(4),'mni');
Affine = spm_maff8(opts.image(1),3, 1,opts.tpm,Affine,'mni');
opts.Affine = Affine;
% run the segmentation
fprintf('Estimating the bias field etc..\n');
p = spm_preproc8(opts);
% this writes the 'native' segmentations
spm_preproc_write8(p, [ones(6,1) zeros(6,3)], [0 1], [0 1], 0, 0, nan(2,3), nan);
[pathname,name,ext] = fileparts([cfg.intermediatename '_anatomy.nii']);
filename = fullfile(pathname, ['m' name ext]);
VO = spm_vol(filename);
VO.dat = spm_read_vols(VO);
tmp = zeros([size(VO.dat) 6]);
for k = 1:6
cfilename = fullfile(pathname, ['c' num2str(k) name ext]);
Vtmp = spm_vol(cfilename);
tmp(:,:,:,k) = spm_read_vols(Vtmp);
end
tmp_tc = false(size(tmp));
for k = 1:6
notk = setdiff(1:6,k);
for m = notk(:)'
tmp_tc(:,:,:,k) = tmp(:,:,:,m)<tmp(:,:,:,k) & tmp(:,:,:,k);
end
end
diagn = zeros(6,5);
for k = 1:6
tmpvals = tmp(:,:,:,k);
tmpvals = tmpvals(tmp_tc(:,:,:,k));
diagn(k,:) = [sum(sum(sum(tmp_tc(:,:,:,k)))) mean(VO.dat(tmp_tc(:,:,:,k))) std(VO.dat(tmp_tc(:,:,:,k))) mean(tmpvals) std(tmpvals)];
end
mri_unbias = keepfields(mri, {'coordsys', 'dim', 'transform', 'unit', 'inside'});
mri_unbias.anatomy = VO.dat;
mri_unbias.params = p;
mri_unbias.diagn = diagn;
otherwise
ft_error('unsupported spmversion requested');
end
% flip the volumes back according to the changes introduced by align_ijk2xyz
for k = 1:3
if flipflags(k)
mri_unbias.anatomy = flipdim(mri_unbias.anatomy, k);
end
end
if ~all(permutevec == [1 2 3])
mri_unbias.anatomy = ipermute(mri_unbias.anatomy, permutevec);
mri_unbias.dim = size(mri_unbias.anatomy);
end
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous mri
ft_postamble provenance mri_unbias
ft_postamble history mri_unbias
ft_postamble savevar mri_unbias