From 41a3ae563441a365d78db85eca180e1d1a4f10d8 Mon Sep 17 00:00:00 2001 From: Siemon de Lange Date: Wed, 19 Apr 2023 13:36:38 +0200 Subject: [PATCH] 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 | 50 +++++++++++++++++++ .../reconstruction_fibers.m | 6 +-- .../test_reconstruction_fibers.m | 27 +++++++--- 3 files changed, 70 insertions(+), 13 deletions(-) create mode 100644 src/misc/nifti_orientation.m diff --git a/src/misc/nifti_orientation.m b/src/misc/nifti_orientation.m new file mode 100644 index 0000000..9ed18ed --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/structural_pipeline/Reconstruction-Fibers/reconstruction_fibers.m b/src/structural_pipeline/Reconstruction-Fibers/reconstruction_fibers.m index 206b3d8..4968031 100644 --- a/src/structural_pipeline/Reconstruction-Fibers/reconstruction_fibers.m +++ b/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 diff --git a/tests/structural_pipeline/test_reconstruction_fibers.m b/tests/structural_pipeline/test_reconstruction_fibers.m index 2e7d3d1..a843e23 100644 --- a/tests/structural_pipeline/test_reconstruction_fibers.m +++ b/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