Add compatibility for complex vox2ras matrices (#33)

The previous implementation assumed a single term per column in the vox2ras matrix of the segmentationFile NIfTI. However, if a vox2ras matrix with multiple terms was used,  CATO would throw an error. This update adds compatibility for such complex vox2ras matrices. Additionally, the tests have been updated to include these complex vox2ras matrices (with mri_info as reference).

src/misc/nifti_orientation.m

@@ -0,0 +1,50 @@
+function orientationStr = nifti_orientation(vox2ras)
+% NIFTI_ORIENTATION computes the orientation string from the vox2ras
+% matrix.
+%
+% Input:
+%   vox2ras: a 4x4 matrix (typically the vox2ras parameter from the
+%   NIFTI file header).
+%
+% Output:
+%   orientation_str: a string indicating the orientation (e.g. 'RAS' or
+%   'LPI').
+
+% Notes:
+% For more information on the vox2ras transformation see the FreeSurfer
+% documentation:
+% https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems
+
+% Extract the direction cosine matrix (3x3)
+R = vox2ras(1:3, 1:3);
+
+% Compute the orientation codes
+orientationCode = zeros(1, 3);
+for i = 1:3
+    [~, maxIndx] = max(abs(R(:,i)));
+    if R(maxIndx, i) > 0
+        orientationCode(i) = maxIndx;
+    else
+        orientationCode(i) = -maxIndx;
+    end
+end
+
+% Convert orientation codes to orientation string
+orientationStr = '???';
+for i = 1:3
+    switch orientationCode(i)
+        case 1
+            orientationStr(i) = 'R';
+        case -1
+            orientationStr(i) = 'L';
+        case 2
+            orientationStr(i) = 'A';
+        case -2
+            orientationStr(i) = 'P';
+        case 3
+            orientationStr(i) = 'S';
+        case -3
+            orientationStr(i) = 'I';
+    end
+end
+end

src/structural_pipeline/Reconstruction-Fibers/reconstruction_fibers.m

@@ -39,11 +39,7 @@ function reconstruction_fibers(configParams)
 segmentationVol = data.vol;
 
 % Extracte orientation from transformation
-props = data.vox2ras(1:3, 1:3);
-[I, ~] = find(props); 
-I = (I-1)*2+ 1 + (props([0 3 6]' + I) > 0);
-orientationOptions = ['L' 'R', 'P', 'A', 'I', 'S']';
-orientation = orientationOptions(I)';
+orientation = nifti_orientation(data.vox2ras);
 
 clear data
 

tests/structural_pipeline/test_reconstruction_fibers.m

@@ -38,17 +38,21 @@ function testOrientation(testCase, subjectDir)
             cp = parseConfigParams(cp);
 
             sform = eye(3);
-            sform(:, :, 2) = [-1 0 0; 0 1 0; 0 0 1];
-            sform(:, :, 3) = [-1 0 0; 0 0 1; 0 1 0];
+            sform(:, :, 2) = [0 -1 0; 1 0 0; 0 0 1];
+            sform(:, :, 3) = ...
+                            [[   -2.0000    0.0000   -0.0000]; ...
+                            [    0.0000    1.7560   -0.9574]; ...
+                            [         0    0.9574    1.7560]];
             sform(:, :, 4) = [0 0 -1; 0 -1 0; -1 0 0];
-            ref = {'RAS', 'LAS', 'LSA', 'IPL'};
+            ref = {'RAS', 'ALS', 'LAS', 'IPL'}; 
 
             for iPerm = 1:size(sform,3)
 
                 NT = load_nifti(cp.structural_preprocessing.segmentationFile);
-                NT.srow_x = [sform(1,:, iPerm) 0];
-                NT.srow_y  = [sform(2,:, iPerm) 0];
-                NT.srow_z  = [sform(3,:, iPerm) 0];
+                NT.srow = sform(:, :,iPerm);
+                NT.srow_x = [sform(1,:, iPerm) 0]';
+                NT.srow_y  = [sform(2,:, iPerm) 0]';
+                NT.srow_z  = [sform(3,:, iPerm) 0]';
                 NT.pixdim = [1 1 1 1 1000 1 1 1]';
                 save_nifti(NT, strrep(cp.structural_preprocessing.segmentationFile, ...
                     'CATO_ref', 'DWI_processed_test'));
@@ -65,10 +69,17 @@ function testOrientation(testCase, subjectDir)
 
                 obs = strtrim(strip(header.voxel_order, char(0)));
 
-                fprintf('Check orientation: %s', ...
+                fprintf('Reference orientation: %s\n', ...
                     ref{iPerm});
+
+                % Show the orientation according to mri_info
+                fprintf('mri_info:');
+                system(sprintf(['source ~/.bashrc; ' ...
+                    'mri_info --orientation %s | tail -n1'], ...
+                    strrep(cp.structural_preprocessing.segmentationFile, ...
+                    'CATO_ref', 'DWI_processed_test')));
 
-                testCase.verifyEqual(ref{iPerm}, obs, 'Orientation not correct.')
+                testCase.verifyEqual(obs, ref{iPerm}, 'Orientation not correct.')
 
             end