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