% load dataset
t1 = xff('*.hdr', 'Please select the dataset to be corrected...');

% apply two rounds of correction with default parameters and PO := 3
% - no mask image is available, auto-detection required!
% - first round will correct basic inhomogeneity (crude mask)
% - second round will correct global inhomogeneity (full mask)
ihc = pmbfilter(t1.VoxelData(:, :, :), 3);
ihc = pmbfilter(ihc, 3);

% apply one last round with
% - slightly lower cutoff (retaining a larger portion of the data)
% - forcing a contiguous mask (only voxels of white matter)
% - using robust regression (takes only slightly longer)
ihc = pmbfilter(ihc, 3, [], struct( ...
    'bcutoff', 0.1, ...
    'cmask',   true, ...
    'robust',  true));

% set back into data
t1.VoxelData = int16(round(ihc));

% save as
t1.SaveAs(strrep(t1.FilenameOnDisk, '.hdr', '_IHC.hdr'));
t1.SaveVoxelData;