[FIX] Fix the way the voxel size is computed (#1053)

* change the way the voxel size is conputed

* check voxel size before running preproc

* use tmp dir

Author: Remi-Gau

bidspm.m

@@ -278,7 +278,7 @@ function preprocess(args)
                       'dummyScans', opt.dummy_scans, ...
                       'force', false);
   end
-
+  bidsCheckVoxelSize(opt);
   if opt.useFieldmaps && ~opt.anatOnly
     bidsCreateVDM(opt);
   end

src/utils/getFuncVoxelDims.m

@@ -28,13 +28,9 @@
     return
   else
     hdr = spm_vol(fullfile(subFuncDataDir, fileName));
-    voxDim = diag(hdr(1).mat);
-    % Voxel dimensions are not pure integers before reslicing, therefore
-    % round the dimensions of the functional files to the 1st decimal point
-    voxDim = abs(voxDim(1:3)');
-    voxDim = round(voxDim * 10) / 10;
-    % Add it to opt.funcVoxelDims to have the same value for
-    % all subjects and sessions
+    M = hdr(1).mat;
+    voxDim  = sqrt(sum(M(1:3, 1:3).^2));
+    % Add it to opt.funcVoxelDims to have the same value for all sessions
     opt.funcVoxelDims = voxDim;
   end
 

src/workflows/bidsCheckVoxelSize.m

@@ -0,0 +1,75 @@
+function bidsCheckVoxelSize(opt)
+  %
+  % Check that all file to preprocess have the same voxel size
+  %
+  % USAGE::
+  %
+  %  status = bidsCheckVoxelSize(opt)
+  %
+  % :param opt: Options chosen for the analysis.
+  %             See also: ``checkOptions()`` and ``loadAndCheckOptions()``.
+  % :type opt: structure
+  %
+
+  % (C) Copyright 2023 bidspm developers
+
+  % If no normalisation is to be done,
+  % checking that all files have the same res is not as important
+  if ~ismember('IXI549Space', opt.space)
+    return
+  end
+
+  if ~isempty(opt.funcVoxelDims)
+    return
+  end
+
+  opt.pipeline.type = 'preproc';
+
+  [BIDS, opt] = getData(opt, opt.dir.preproc);
+
+  allFiles = {};
+  allVoxDim = [];
+
+  for iSub = 1:numel(opt.subjects)
+
+    subLabel = opt.subjects{iSub};
+
+    for iTask = 1:numel(opt.taskName)
+
+      opt.query.task = opt.taskName{iTask};
+
+      filter = fileFilterForBold(opt, subLabel, 'stc');
+
+      files = bids.query(BIDS, 'data', filter);
+
+      for iFile = 1:numel(files)
+        allFiles{end + 1, 1} = files{iFile}; %#ok<*AGROW>
+        [subFuncDataDir, boldFilename, ext] = spm_fileparts(files{iFile});
+        voxDim = getFuncVoxelDims(opt, subFuncDataDir, [boldFilename, ext]);
+        allVoxDim(end + 1, :) = voxDim;
+      end
+
+    end
+
+  end
+
+  differentFromFirstImageRes = (allVoxDim - allVoxDim(1, :)) > 10^-3;
+  if any(differentFromFirstImageRes(:))
+    msg = sprintf('%s\n', 'Not all input files have the same resolution.');
+    for i = 1:numel(allFiles)
+      msg = [msg, sprintf('\n%s: %0.3f, %0.3f, %0.3f', ...
+                          strrep(allFiles{i}, fullfile(BIDS.pth, filesep), ''), ...
+                          allVoxDim(i, 1), ...
+                          allVoxDim(i, 2), ...
+                          allVoxDim(i, 3))];
+    end
+    msg = [msg, sprintf('\n\n%s\n', ...
+                        ['Use the "opt.funcVoxelDims" to force the voxel dimensions', ...
+                         'to have the same resolution after normalisation.'])];
+    logger('ERROR', msg, ...
+           'options', opt, ...
+           'filename', mfilename(), ...
+           'id', 'differentVoxelSize');
+  end
+
+end

tests/tests_batches/test_saveMatlabBatch.m

@@ -13,43 +13,43 @@ function test_saveMatlabBatch_basic()
   subLabel = '01';
   opt = setOptions('dummy', subLabel);
   opt.dryRun = false;
-  opt.dir.jobs = pwd;
+  tmpDir = tempdir;
+  spm_mkdir(tmpDir);
+  opt.dir.jobs = tmpDir;
 
   matlabbatch = struct('test', 1);
   saveMatlabBatch(matlabbatch, 'test', opt, subLabel);
 
-  expectedOutput = fullfile(pwd, ['sub-' subLabel], ...
+  expectedOutput = fullfile(tmpDir, ['sub-' subLabel], ...
                             ['batch_test_' datestr(now, 'yyyy-mm-ddTHH-MM') '.mat']);
 
   assertEqual(exist(expectedOutput, 'file'), 0);
   assertEqual(exist(spm_file(expectedOutput, 'ext', '.json'), 'file'), 2);
 
-  expectedOutput = fullfile(pwd, ['sub-' subLabel], ...
+  expectedOutput = fullfile(tmpDir, ['sub-' subLabel], ...
                             ['batch_test_' datestr(now, 'yyyy_mm_ddTHH_MM') '.m']);
 
   assertEqual(exist(expectedOutput, 'file'), 2);
 
-  cleanUp(fullfile(pwd, ['sub-' subLabel]));
-
 end
 
 function test_saveMatlabBatch_group()
 
   subLabel = '01';
   opt = setOptions('dummy', subLabel);
   opt.dryRun = false;
-  opt.dir.jobs = pwd;
+  tmpDir = tempdir;
+  spm_mkdir(tmpDir);
+  opt.dir.jobs = tmpDir;
 
   matlabbatch = struct('test', 1);
 
-  expectedOutput = fullfile(pwd, 'group', ...
+  expectedOutput = fullfile(tmpDir, 'group', ...
                             ['batch_groupTest_' datestr(now, 'yyyy-mm-ddTHH-MM') '.mat']);
 
   saveMatlabBatch(matlabbatch, 'groupTest', opt);
 
   assertEqual(exist(expectedOutput, 'file'), 0);
   assertEqual(exist(spm_file(strrep(expectedOutput, '-', '_'), 'ext', '.m'), 'file'), 2);
 
-  cleanUp(fullfile(pwd, 'group'));
-
 end